How to calculate geometric tortuosity of a 3D binary image in all directions (x,y,z)

Hi,

I'd like to calculate the geometric tortuosity of a 3d binary image (2d .tiff stack) in in the x, y, and z directions using MATLAB. What I have is a 3d binary image (obtained from xray u-CT) of a porous material and I would like to calculate the tortuosity in all directions of the sample. Please help?

Answer 1

[Sean de Wolski]

Here's the general approach I would take:

• For every pore, identify the end-points. We'll get to how to do this later.
• Once you have the end-points, select one, take the geodesic distance transform with this point selected and then calculate the maximum value in that pore (which I believe should always be the other end point).

doc bwdistgeodesic  %will give you this transform in three dimensions.

• Divide the geodesic distance by the hypotenuse between the two end points to calculate tortuosity.

The second module of this webinar should help you get an understanding of the general idea as it is exactly the same thing in two dimensions:

https://www.mathworks.com/company/events/webinars/wbnr66293.html

Now, how to calculate endpoints in three dimensions? This is non-trivial. There are no prebuilt 3d skeletonization or thinning routines ( bwmorph for 2d) to my knowledge. If this was my task, I would probably try something like the following as a starter.

• Use bwperim to get just the perimeter voxels of the pores (they will contain endpoints).
• Then loop over the voxels in the perimeter taking the geodesic distance transform with respect to each voxel.
• Calculate the maximum value and store it only if the max is greater than the previous max. This is a brute force approach to finding the endpoints (i.e. the two points that maximize distance along the pore. I can't think of a cheap way around this though because a weird shape could throw off anything general.
• Once you have the endpoints, the geodesic distance between them, calculate the distance between them using the distance formula. Divide these two and then relax.

Maybe Image Analyst will have some additional thoughts as well.

[Dami]

Thanks Sean!

[Image Analyst]

Well the common definition for tortuosity is the path distance divided by the straight line distance like Sean said. And like Sean says, the skeleton in 3D is not trivial. And it depends on what your binary object looks like and what you want. The skeleton of an n-dimensional image is a (n-1) object. That makes sense for 2D where a 2D blob will give a bunch of 1D lines, though possibly connected to each other. Let's say your 3D object is like staghorn coral. You might reasonably expect or want your skeleton to be the main spine/backbone of the coral, perhaps with little stray struts going off of it to every little tiny bump on the surface of the coral. (Aside: I've asked the IP team for a routine that would give just the main, longest line and not the little stray spurs sticking off of it, but it's really not so easy and it's ill defined, especially if there are loops. Anyway, that's a side topic.) Back to the main point, you might want a 1D line in 3 dimensions for that coral, but remember what I said about the skeleton being (n-1)D? So the skeleton in 3D is actually a 2D surface, not a 1D backbone. Once I wanted the 3D skeleton of a 3D microCT object with void spaces (like bubbles) in it (think of an object like bread or some other closed cell foam). I wanted the sheets surrounding each void space so I could multiply that by the Euclidean Distance Transform to get the average thickness of the walls between the voids. (Sorry if I've lost you, I know it's complicated). First I tried the skeletonization routine in ITK and it didn't do what I want (I've forgotten why). I didn't know how to do it so I asked Steve Eddins. Steve ended up writing my a 3D skeletonization routine using a clever algorithm that works like it theoretically should. Even an object like a staghorn coral will give 2D surfaces for the skeleton? Why? Well let's imagine a very simple object like a rectangular block. What would the skeleton of that look like? A line? No! How about 9 lines (a centerline and 4 lines going to each corner)? Again, no. Well then what the heck does it look like? Imagine this. The 3D skeleton is the surface the center of a ball would trace out as it ran around the inside of the object, changing its radius to be the biggest it can be. So what about when the ball is in the center of a flat rectangular slab? With it's diameter the thickness of the narrowest dimension of the slab, the center of the ball would trace out a 2D sheet. And each end of the block would have 4 sheets attached to the main sheet along the main body of the slab/block. Sorry I don't have an image of that. The closest I could find is http://www.inf.u-szeged.hu/~palagyi/skel/voroskel.gif

So, the skeleton of a 3D object is a 2D surface. To get the 2D tortuosity you take the path length and divide it by the straight line distance. What's the equivalent in 3D? Is it the surface area divided by the volume or something? Does a 3D equivalent even exist or make sense? Does the answer depend on what Dami's object looks like? The desired 3D skeleton of a skull my look different than the desired 3D skeleton of a finger bone. Can you get a 1D skeleton from a 3D object? Yes but you have to make some assumptions.

Anyway, I guess that's just a long winded answer of saying "Ask Steve!" Here's the start of his code he wrote after he and I discussed it:

%------------------------
% Compute the 3D skiz using Steve Eddins proposed method.
% There is a way of using a 3-D watershed transform that I think matches closely your inflating balls metaphor.
% Specifically: 
% 1.       Compute the distance transform of your 3-D blobs image using bwdist. 
%          (Fortunately, bwdist runs much faster and takes up less memory for the multidimensional case in R2009b than in earlier releases.)
% 2.       Compute the watershed transform of the distance transform using watershed.  
%          The zero-valued pixels in the output of watershed will be at the interfaces between the inflated balls.
% I think this is what morphology people call the “SKIZ,” or “skeleton by influence zones.”  The interior of each inflated ball is an “influence zone.”
% Okay, so here we go......

And if anyone followed everything I said, I'd be really impressed! Because it is confusing, I know. Luckily I know some really smart people on the Image Processing team at the Mathworks who can help me out with the most challenging problems.