Estimation of tortuosity and reconstruction of geodesic paths in 3D

The morphological tortuosity of a geodesic path in a medium can be deﬁned as the ratio between its geodesic length and the Euclidean distance between its two extremities. Thus, the minimum tortuosity of all the geodesic paths into a medium in 2D or in 3D can be estimated by image processing methods using mathematical morphology. Considering a medium, the morphological tortuosities of its internal paths are estimated according to one direction, which is perpendicular to both starting and ending opposite extremities of the geodesic paths. The used algorithm estimates the morphological tortuosities from geodesic distance maps, which are obtained from geodesic propagations. The shape of the propagated structuring element used to estimate the geodesic distance maps on a discrete grid has a direct inﬂuence on the morphological tortuosity and has to be chosen very carefully. The results of our algorithm is an image with pixels p having a value equal to the length of the shortest path containing p and connected to two considered opposite boundaries A and B of the image. The analysis of the histogram of the morphological tortuosities gives access to their statistical distribution. Moreover, for each tortuosity the paths can be extracted from the original image, which highlights the location of them into the sample. However, these geodesic paths have to be reconstructed for further processing. The extraction, because applying a threshold on the tortuosities, results in disconnected components, especially for highly tortuous paths. This reconstruction consists in reconnecting these components to the geodesic path linking the two opposite faces, by means of a backtracking algorithm


INTRODUCTION
The morphological tortuosity is one of the measurements which have been processed on 3D images to characterize a fibrous material called Thermisorel™, and which is made of wooden fibres.This medium was used as a reference material in the Silent Wall project for its good properties in thermal and acoustical insulating (Peyrega et al., 2009;2010;Peyrega, 2010).The main objective of the Silent Wall ANR1 project consisted in optimizing the acoustic properties of fibrous media by modifying their microstructures.Since the morphological tortuosity of the microscopic pores has a direct influence on the acoustic absorption on such a material at the macroscopic scale, this morphological parameter was studied for this project.This paper focuses on methods and algorithms to estimate the morphological tortuosity of a medium and to reconstruct on discrete grids in 2D and 3D geodesic paths it contains.After presenting an algorithm (Decker et al., 1998;Peyrega, 2010) to estimate the morphological tortuosity of any phase into 2D and 3D images from geodesic dilations (Lantuejoul and Beucher, 1981;Serra, 1988;Soille, 2003), we will focus on the influence of the shape of the propagating structuring elements on this estimation.Then, the studied phases can be separately characterized by histograms of tortuosities in the corresponding directions of propagation.Moreover, since geodesic propagations are handled by the algorithm, it is possible to extract the shortest paths in any direction, i.e., the geodesic paths, linking any point to two separate components which are in the same phase.Finally, these methods were applied to 3D X-Ray CT images of fibrous materials used for acoustic absorption.

Geodesic distance
In order to define the morphological tortuosity, we must start with the definition of the geodesic distance (Lantuejoul and Beucher, 1981;Serra, 1988;Soille, 2003).Considering two points x and y belonging to a set X, the shortest paths connecting them, while constrained to remain in X, are called the geodesic paths in the geodesic mask X.Thus, using the definition (Eq. 1) proposed by (Soille, 2003), with X a connected component of an image, the geodesic distance Geo Dist X (xy) between two points x and y in X is the minimum length L of all the paths P = (p 1 , p 2 , p 3 , ..., p n ) linking x and y included in X.We note p i the i th pixel of the path P between x and y composed of n connected pixels belonging to X.
Geo Dist X (xy) = min {L (P) |p 1 = x, p n = y, P ⊆ X} . (1) The geodesic distance between two points x and y belonging to a medium X is illustrated by Fig. 1.According to this definition, the distance between x and w is infinite because both points belong to two different connected components of the set X.One of the challenges of this work consists in estimating a geodesic distance map on a discrete grid.Previous publications proposed methods to estimate it in different contexts (Borgefors, 1986;Soille, 1992;1994;2003;Coeurjolly et al., 2004).
From a morphological point of view, this can be processed with geodesic dilations from x to y into X.In Fig. 2, the propagating fronts are circular, which is intuitive to estimate a Euclidean distance map in a continuum framework.However, in this paper, working on digitized images on a grid of points, we show to what extent the shape of the structuring element influences the estimation of the geodesic distance map.

Morphological tortuosity
We consider two separate subsets A and B of a medium X, that will be used as propagation sources.The geodesic distance of a point x ∈ X to a subset A is defined by Lantuejoul and Beucher (1981), and by Soille (2003) in Eq. 2.

Geo Dist
The geodesic paths connecting x to A are the paths corresponding to Geo Dist X (x, A) (Soille, 2003).In the remaining part of the paper, we consider two parallel planar sources A and B in 3D, and their separation noted Euclidean Dist (AB).For every point x ∈ X, the morphological tortuosity corresponding to the sources A and B is defined as the ratio Eq. 3.
In what follows, the morphological tortuosity of every point of a set X (namely fibres or pores) will be estimated by means of geodesic dilations (Lantuejoul and Beucher, 1981;Serra, 1988;Soille, 2003) implemented on a grid of points, since we are using 3D digitized images of a real material.

Acoustic tortuosity in physics
In acoustics, the tortuosity α(ω) has a physical interpretation linked to both the frequency ω of the acoustic wave, and to the geometry of the microstructure, through which the wave front is propagating.Moreover, we note α ∞ the acoustic tortuosity for sound waves having an infinite frequency.According to Allard (1993) and Allard et al. (1994), we can write the relationship between the tortuosity α(ω) and the acoustic velocity u of the propagating wave through the microstructure into Eq.4.
(4) Thus, the physical interpretation of α(ω) is directly linked to both the average acoustic velocity of the wave and to the morphological tortuosity of the microstructure.For straight paths, both morphological and acoustic tortuosities are equal to 1 because in this case, we can write:
Consider an elementary structuring element defined on the grid of points from a point and its nearest neighbours.It defines a connectivity on the graph generated by the grid, for instance the C 4 and the C 8 connectivities starting from a square grid and the 4 or the 8 nearest neighbours, the C 6 and the C 26 connectivities on a cubic grid by means of the 6 or the 26 nearest neighbours.Noting δ (A) the dilation of a set A by the elementary structuring element, and δ (1) X (A) the geodesic dilation of size 1 defined from the set X, we have: The geodesic dilation of A of size n in X is obtained after n iterations of δ (1) X (A).From its definition, the result of the geodesic dilation depends on the choice of the elementary structuring element and of its corresponding connectivity.In the rest of the paper, we will compare geodesic paths and geodesic distances obtained on images with different elementary structuring elements.
In order to estimate the morphological tortuosities of every path linking a voxel v of a set X to two opposite faces A and B of X (Eq.3), the algorithm proposed in 3D by Decker et al. (1998) has been used.It consists of the following steps using morphological operations such as geodesic dilations (Serra, 1988)   The steps 1 and 2 calculate the geodesic distance of each pixel of the set X respectively to the top and to the bottom faces.After these two steps, the components in pink are eliminated because they are not connected to a marker-face and consequently not reached by the propagating fronts in the vertical Oy direction.
The results of step 3, ImAddFwdBwd, is an image where the value of each pixel of the white component of X is the length of the geodesic path linking the two marker-faces and going through the corresponding pixel.However, the pixels in the green components of X are just linked to one marker-face and are then eliminated by steps 4 and 5, which isolate all the pixels belonging to the paths linking the two marker-faces.The pixels of the background are set to zero elsewhere.Finally the resulting image ImTortuosity, giving for each pixel the tortuosity of the shortest paths in which it is contained, is then normalized by the Euclidean distance between both marker-faces.

INFLUENCE OF THE PROPAGATING STRUCTURING ELEMENT TWO DIFFERENT STRUCTURING ELEMENTS
The geodesic dilations on a digitized image are morphological operations and then directly depend on the shape of the structuring element used (Serra, 1988).Previous publications (Soille, 1992;1994;2003;Coeurjolly et al., 2004) highlighted the influence of the propagating structuring element on the geodesic distance map.The images processed in the framework of the Silent Wall project are in 3D.That is why two 3dimensional structuring elements have been compared, on the one hand the cube with 26 neighbors, and on the other hand the sphere.
Considering the problem of acoustics handled in the Silent Wall project, the assumption is made that the wave fronts propagate only outwards in the direction perpendicular to the marker-faces and do not return to the source face.Thus the structuring elements are not isotropic.On the contrary, they are half cubes or half spheres.The Fig. 5 illustrates structuring elements with a radius equal to 15 voxels obtained from 3D dilations of the point-pattern shown in Fig. 4.

SQUARES AND CUBES
The 8-connectivity square in 2D We first used the 8-connectivity square and 26-connectivity cube to estimate the morphological tortuosity of fibrous media respectively for 2D and 3D images.However their shapes involve problems of inaccuracy because of their square corners, which do not correctly suit to estimate the geodesic distance between two faces of a sample.From a physical point of view, the cubic propagating fronts do not estimate a geodesic distance faithful to the definition illustrated in Fig. 2.
The Fig. 6 shows how rough the estimation of the tortuosities in the white part represented in Fig. 3 can be, when it is processed with squares as structuring elements.The histogram of the tortuosities of this white part (Fig. 6c) shows that 45% of the pixels of the path have a tortuosity equal to 2.70, with a maximum tortuosity equal to 2.89.Moreover according to the Fig. 6b, the paths with a minimum tortuosity equal to 2.70 (in yellow) are physically incorrect because they do not follow the curves of the white part in Fig. 3.
To sum up, the 2D square propagating front underestimates the geodesic distances and then underestimate the tortuosities.This is observable on the histogram which has consequently a high number of pixels with a tortuosity close to 1.This problem is observable in 3D with cubes as well.The 26-connectivity cube in 3D In 3D the cubic fronts involve a too high number of pixels having a tortuosity close to 1. Consider a straight composite tube composed of concatenated big and small cylinders (Fig. 7).The Fig. 7a represents a cropped 3D view of the real object in order to see the fronts inside it.
Thus discontinuities of the fronts are introduced at the interfaces between big and small cylinders.Since the structuring element is cubic, these discontinuities do not deform the propagating fronts enough, that is why most of them look like plane waves especially in the first input cylinder on the left (Fig. 7b).
As shown in Fig. 7c and Fig. 7d, 84% of the voxels have a tortuosity equal to 1 which makes the cube blind to the discontinuities of the composite tube, when it is used as propagating structuring element.

FAST MARCHING DISCRETIZED DISCS AND SPHERES
The square in 2D and the cube in 3D are irrelevant to estimate the morphological tortuosity from geodesic dilations.Intuitively, the disc in 2D and the sphere in 3D are the solutions to this problem.However both must be discretized on 2D or 3D lattices to be handled in the image processing algorithm.

The fast marching method
The fast marching was originally introduced by Sethian (1996).It consists of iteratively solving the Eikonal equation (Eq.5) into a narrow band around the marker seed, with T (p) the arrival time at a point p of the propagating front starting from the marker seed p 0 (T (p 0 ) = 0), and f the cost function related to the speed of propagation: (5) The cost function can be compared to the refractive index n in geometrical optics (Cohen and Kimmel, 1997;Petres et al., 2005).Thus, the speed v(p) of the front at the point p can be written as a function of the celerity of light: In the case of this study, the cost function is homogeneous and equal to 1, which leads to fronts propagating at the speed of one pixel (in 2D) and one voxel (in 3D) per iteration.The fronts travel outwards and do not return to the source of the propagation.
As explained by Sethian (1996), solving the Eikonal equation (Eq.5) with the fast marching algorithm is made by solving the second order Eq.6 in 2D and Eq. 7 in 3D for T(p), with f i, j and f i, j,k , the cost functions respectively at the pixel (i, j) and at the voxel (i, j, k).We note T i, j and T i, j,k the arrival time T(p) at the point p respectively (i, j) in 2D and (i, j, k) in 3D.Thus, the algorithm only handles the 4-connectivity neighborhood in 2D and the 6connectivity neighborhood in 3D of the voxels into the narrow band.
The algorithm itself is precisely described and illustrated by Sethian (1996), Cohen and Kimmel (1997) and Petres et al. (2005).An implementation using hierarchical priority queues (algorithms 12 and 13) is proposed by (Peyrega, 2010).Moreover, if the cost function f = 1 inside the geodesic mask and if f = ∞ outside, then the speed is equal to 1 voxel per iteration, and the time of arrival T (p) at a point p inside the geodesic mask is equivalent to the geodesic distance U(p) between p and the marker seed inside the mask.

Discs in 2D
To be compared to Fig. 6a, the Fig. 8a represents the forward propagation from the top to the bottom of the white part represented in Fig. 3, obtained by successive geodesic dilations by discs.The fronts have circular shapes approximating the geodesic distance on the 2D discretized grid.The resulting minimum tortuosity path in yellow in Fig. 8b is more acceptable that the one in Fig. 6b, because it correctly follows the interior of the curves of the white part in Fig. 3.The minimum tortuosity reached is about 3.01.
The histogram of the tortuosities estimated with discs (Fig. 8c) is completely different that the one obtained with a square.It is totally shifted to higher tortuosities between 3.01 and 3.23 and the reached tortuosities are more regularly distributed since no prominent maximum is observable like in Fig. 6c.

Spheres in 3D
Three dimensional propagations with spheres in the straight composite tube are illustrated in Fig. 9a.Contrarily to the Fig. 7b, the spherical propagating fronts in Fig. 9b are influenced by the discontinuities between the large cylinders, even if the fronts in the first input cylinder on the left look like plane waves.Moreover, as soon as the first discontinuity after the first input cylinder has been passed, the propagating fronts do not recover their plane wave shapes.
On the contrary, according to Fig. 9c and Fig. 9d the sphere as structuring element is sensitive to the discontinuities of the paths through which the fronts propagate.On the one hand the yellow path of minimum tortuosity equal to 1 in Fig. 9c is logically exclusively located around the axis of the composite tube but does not touch the sides of the cylinder as it was observable in Fig. 7c.On the other hand the histogram of Fig. 9d still highlights a high number of voxels with a tortuosity equal to 1 because of the global straightness of the tube, but it is not blind to higher tortuosities, as the cube is.
Finally estimating the tortuosities of the paths through a medium is highly dependent on the shape of the structuring element used to process the geodesic distances.Thus the best choice are the disc in 2D and the sphere in 3D which are more consistent from a physical point of view to simulate propagating fronts.Moreover they are more sensitive to the discontinuities and changes of directions of the paths through which they propagate.This makes them consequently more accurate in estimating the tortuosities and especially high tortuosities.

HALF DISCS AND HALF SPHERES
Considering acoustics, we make the assumption that the outward propagating fronts globally travel in a specific direction perpendicular to the origin faces.Thus half spheres are finally used to estimate the morphological tortuosity of fibrous media in both phases (fibres and pores) for the Silent Wall project.Using half spheres highlights the direct paths linking the opposite faces.For example, in order to propagate half discs or half spheres in the increasing Ox direction (noted Ox + ), the Eikonal equation solved into the Fast Marching (algorithm 12 and 13) is modified as Eq. 8 in 2D, and as 9 in 3D.
Let X1 and X2 be two sets respectively illustrated in Fig. 10 and 11.The set X1 is the same as in Fig. 3, and X2 is slightly different.Like in Fig. 3, the pink zones are not reached by the propagating fronts in the vertical Oy direction between both faces A and B. The green areas are just linked to one marker-face (A or B).Finally, the white zones are the percolating paths between both marker-faces in the direction of propagation.Thus, using half spheres inside the set X1 (Fig. 10b and 11b) leads to nonpercolating paths in the vertical Oy direction, and then prevents us from estimating the morphological tortuosity of it with such a structuring element.On the contrary, isotropic spheres should be used to process the geodesic distance maps (Fig. 10a and 11a ).

RECONSTRUCTION OF GEODESIC PATHS THRESHOLDING THE TORTUOSITY
In Fig. 8b, the value of each pixel p is equal to the morphological tortuosity of the geodesic paths containing p, connecting it to the opposite marker faces A and B in the vertical Oy direction.Thresholding this image at a value higher than the minimal tortuosity inside the yellow path, equal to 3.01 (Fig. 8), will result in components disconnected from the marker faces A and B (for example the red parts).Thus these disconnected components must be reconnected to the marker faces in order to extract the full geodesic paths from their tortuosities, for further morphological characterization or for further processing like flow simulation.These components are reconnected to the paths with the minimum tortuosity.These reconstructions are processed from source backtracking.

The algorithm
Since the fast marching algorithm processes the voxels in the narrow band around the propagating front, it is possible to extract the source voxel in the 4-connectivity (2D) or 6-connectivity (3D) neighborhood of each voxel of the narrow band.This source voxel v is the last one from which the time of arrival at a voxel vn in the narrow band is updated, and consequently is the nearest voxel from the starting markers.The source backtracking is directly processed during the propagation, which does not extend the processing time.

1.3
Allocate and initialize all the voxels v of ImLabels Initialize all the voxels v of ImSrc4C 6C at 0; Add vn into FM PQ with the 1.17   Table 1.Coding square for the 8-connectivity neighbors vn of v in 2D.
The instructions in line 1.12 and line 2.10 (respectively in algorithms in Figs. 12 and 13) mean that the value of a voxel vn of the image ImSrc4C 6C is equal to the identification number of the relative position (RP) of the current voxel v compared to its 4connectivity (2D) or 6-connectivity (3D) neighboring voxel vn.

4-connectivity (2D) backtracking
In 2D, the relative positions of the neighbors vn of v in 4-and 8-connectivities are coded from the Table 1.Thus, the global principle of our source coding, consists in setting the value of each pixel v of the geodesic path, to the relative position (7 − vn) of the neighboring pixel vn which has the smallest geodesic distance to the starting image boundary.

4-Connectivity reconstruction
On both Fig. 15 and 16 the disconnected components in green are reconnected to both top and bottom faces using the previous method in 4-connectivity.Since the source backtracking is processed only in a neighborhood of size 1, and since this processing is made on discrete cubic grids, the reconstructed paths are approximations of the minimum geodesic paths linking the disconnected components to the marker faces.

8-Connectivity reconstruction
As shown in Figs. 17 and 18, the path reconstructions in 8-connectivity are fairly better approximations of the minimum geodesic paths than in 4-connectivity.In 8-connectivity, the reconnected paths better follow the curves of the geodesic mask.

APPLICATION TO 3D IMAGES OF FIBROUS MEDIA TORTUOSITIES OF THERMISOREL
For illustration of this work, we show an application to 3D X-Ray CT images of a fibrous network.In the framework of the Silent Wall project, the reference fibrous material is the Thermisorel™, which has very good thermal and acoustic insulation properties.The Fig. 19 illustrates a 3D X-Ray CT-image of this material.Among the different morphological measurements made on such 3D images of Thermisorel™ (Peyrega et al. 2010;Peyrega 2010), the tortuosity was estimated with the previously described algorithm using Fast Marching propagations of half spheres in the Ox, Oy and Oz directions.
The histograms of tortuosities in Fig. 20 gives access to the distribution of the connectivity into the considered media (fibrous and porous phases).The Tables 2 and 3 give more details about the extremal and average morphological tortuosities, and about their standard deviations respectively inside the fibres and the pores of Thermisorel™.A first look at these histograms shows that the tortuosity is higher for both media in the Oz direction than in the longitudinal ones (Ox and Oy).This is a direct consequence of the production method of this material, which is made from a paper like process.The board is compressed in the Oz direction, thus the fibres are globally isotropically oriented in the xOy planes perpendicularly to this axis (Peyrega et al. 2010;Peyrega 2010).This explains why the morphological tortuosities in both media are close to 1 in the xOy planes.Moreover the Fig. 20b shows that the paths inside the pores have fairly small tortuosities close to 1, which proves that the opposite faces are mostly linked by direct paths in the three directions.

PATH RECONSTRUCTIONS IN 3D IMAGES OF THERMISOREL
As shown in Fig. 8, thresholding the 3D images of tortuosities of Thermisorel™ leads to disconnected components.The previously described method of source backtracking used to reconnect these isolated components to the minimum geodesic paths between the opposite faces, in the corresponding directions, has been handled.The results of such reconnections in 26connectivity are illustrated in Fig. 21, 22 and 23 for the fibrous phase, and in Fig. 24, 25 and 26 for the porous phase, where the images of morphological tortuosity have been thresholded at the minimal, average and maximal values.In the pores   Finally, we can see in Fig. 21 to 26 that the reconstructed paths in both phases (fibres and pores) are homogeneously distributed in the 3D sample in the three directions for thresholds equal to the average tortuosities.On the contrary, extremal tortuosities are located in specific zones of the material.Thus these results show that the average tortuosities of the geodesic paths in both phases of a fibrous material such as the Thermisorel™ are geometrically representative enough to be handled as macroscopic characteristic parameters of the material.

CONCLUSION
Applied to 3D X-Ray CT images, our method to reconnect geodesic paths to isolated connected components lets us locate the geodesic paths within the fibres and the pores of fibrous material according to their tortuosities.Our original algorithm is based on source backtracking, which uses the Fast Marching algorithm itself to draw a source map of the propagating front.This source map was then used to reconnect any connected component to the nearest geodesic path to the starting boundaries.

Fig. 1 .
Fig. 1.Geodesic distance in the set X between x and y.

Fig. 2 .
Fig. 2. Geodesic distance in the set X between x and y estimated by geodesic dilations.
to estimate the geodesic distance map: 1) Estimation of the geodesic distance of each voxel of the medium X to the face A: → ImGeoDistForward 2) Estimation of the geodesic distance of each voxel of the medium X to the face B: → ImGeoDistBackward 3) Addition of both images: ImGeoDistForward + ImGeoDistBackward → ImAddFwdBwd 4) Infinum of both images: INF(ImGeoDistForward, ImGeoDistBackward) of ImTortuosity by the Euclidean distance between A and B to get the tortuosity This algorithm is illustrated in 2D by the Fig. 3 but its principle is the same in 3D.Let X be the set composed of every connected components (white, green and pink) different from the background in black.Let the top A and bottom B faces be the so-called marker-faces used as markers to be dilated into the set.

Fig. 3 .
Fig. 3. Estimation of the tortuosity into the white set between the top (blue) and the bottom (red) faces.

Fig. 4 .Fig. 5 .
Fig. 4. Pattern of points on the faces and at the center of a cube.(a) 3D view, (b) 2D projection of the pointpattern (red: visible points ; blue: invisible points ; green: central point invisible on the 3D view.)

Fig. 10 .
Fig. 10.Influence of isotropic and half discs on the paths in X1.(a) Isotropic discs, (b) Half discs.

Fig. 11 .
Fig. 11.Influence of isotropic and half discs on the paths in X2.(a) Isotropic discs, (b) Half discs.
FAR;Initialize all the voxels v of ImDistFM at ∞; 1.5

Fig. 15 .
Fig. 15.Reconnections of the disconnected components in green into the loops to the top and bottom faces.2D vertical 4-connectivity reconstructions along Oy.

Fig. 16 .
Fig. 16.Reconnections of the disconnected components in green into the white part (Fig. 3) to the top and bottom faces.2D vertical 4-connectivity reconstructions along Oy.

Fig. 17 .
Fig. 17.Reconnections of the disconnected components in green into the loops to the top and bottom faces.2D vertical 8-connectivity reconstructions along Oy.

Fig. 18 .
Fig. 18.Reconnections of the disconnected components in green into the white part (Fig. 3) to the top and bottom faces.2D vertical 8-connectivity reconstructions along Oy.

Fig. 20 .
Fig. 20.Histograms of the morphological tortuosities (a) into the fibres and (b) into the pores of Thermisorel™ in Fig. 19 in the Ox, Oy and Oz directions.
Fig. 23.Reconstruction of the paths by backtracking in 26-connectivity after thresholding the morphological tortuosity in the Oz direction for the fibres of Thermisorel™.
Fig. 24.Reconstruction of the paths by backtracking in 26-connectivity after thresholding the morphological tortuosity in the Ox direction for the pores of Thermisorel™.
Fig. 25.Reconstruction of the paths by backtracking in 26-connectivity after thresholding the morphological tortuosity in the Oy direction for the pores of Thermisorel™.
Fig. 26.Reconstruction of the paths by backtracking in 26-connectivity after thresholding the morphological tortuosity in the Oz direction for the pores of Thermisorel™.
Remove the 1 st plateau TP from FM PQ;

Table 2 .
Minimal, maximal, and average values, and standard deviation of the morphological tortuosities into the fibres of Thermisorel™ in Fig.19in the Ox, Oy and Oz directions.

Table 3 .
Minimal, maximal, and average values, and standard deviation of the morphological tortuosities into the pores of Thermisorel™ in Fig.19in the Ox, Oy and Oz directions.