Image Anal Stereol 2013;32:27-43 Original Research Paper
ESTIMATION OF TORTUOSITY AND RECONSTRUCTION OF GEODESIC PATHS IN 3D
Charles Peyregab and Dominique Jeulin
Centre de Morphologie Mathématique, Mathématiques et Systèmes, Mines ParisTech, 35 rue Saint Honoré, 77300 Fontainebleau, France
e-mail: charles.peyrega@mines-paristech.fr, dominique.jeulin@ensmp.fr (Received June 6, 2012; revised January 6, 2013; accepted January 22, 2013)
ABSTRACT
The morphological tortuosity of a geodesic path in a medium can be defined 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 influence 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.
Keywords: 3D images, fibrous media, geodesic paths, mathematical morphology, tortuosity.
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
1 http://us2b.pierroton.inra.fr/Projets/Silent_Wall/description.htm
(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.
ESTIMATION OF THE MORPHOLOGICAL TORTUOSITY
DEFINITIONS 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_DistX (xy) between two points x and y in X is the minimum length L of all the paths P = (p1,p2,p3,...,pn) linking x and y included in X. We note pi the ith pixel of the path P between x and y composed of n connected pixels belonging to X.
Geo_DistX (xy) =
min {L (P) |p1 = x,pn = y,P Q 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 G X to a subset A is defined by Lantuejoul and Beucher (1981), and by Soille (2003) in Eq. 2.
Geo_DistX (x, A) =
min {Geo_DistX (xy),x G X,y G A} (2)
The geodesic paths connecting x to A are the paths corresponding to Geo_DistX (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 JDist (AB). For every point x G X, the morphological tortuosity corresponding to the sources A and B is defined as the ratio Eq. 3.
Fig. 1. Geodesic distance in the set X between x and y.
Fig. 2. Geodesic distance in the set X between x and y estimated by geodesic dilations.
TortuosityXAB (x) =
Geo_DistX (x, A) + Geo_DistX (x, B) Euclidean_Dist (AB)
(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 a (to) has a physical interpretation linked to both the frequency to of the acoustic wave, and to the geometry of the microstructure, through which the wave front is propagating. Moreover, we note a^ 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 a (to) and the acoustic velocity u of the propagating wave through the microstructure into Eq. 4.
ax = lim a (to ) = lim

to
<u2) <u>2
(4)
Thus, the physical interpretation of a (a) 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: (u2) = (u)2.
THE ALGORITHM
The algorithm makes use of the geodesic dilation (Lantuejoul and Beucher, 1981; Serra, 1988; Soille, 2003), defined as follows.
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 C4 and the C8 connectivities starting from a square grid and the 4 or the 8 nearest neighbours, the C6 and the C26 connectivities on a cubic grid by means of the 6 or the 26 nearest neighbours. Noting 8 (A) the dilation of a set A by the elementary structuring element, and 8X(-1)(A) the geodesic dilation of size 1 defined from the set X , we have:
8XX1)(A) = 8(A) nX .
5)	Extraction of the percolating paths: IF ImInfFwdBwd > 0: ImTortuosity = ImAddFwdBwd, ELSE:
ImTortuosity = 0
6)	Normalization 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.
n
J\
I Direction of Forward I propagation
The geodesic dilation of A of size n in X is obtained after n iterations of 8X(-1)(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) 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) ^ ImInfFwdBwd
Fig. 3. Estimation of the tortuosity into the white set between the top (blue) and the bottom (red) faces.
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.
w
v\ r .........
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 etal., 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 3-dimensional 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.
(a)	(b)
Fig. 4. Pattern of points on the faces and at the center of a cube. (a) 3D view, (b) 2D projection of the point-pattern (red: visible points ; blue: invisible points ; green: central point invisible on the 3D view.)
(a)	(b)
Fig. 5. 3D structuring elements of size 15 propagating in the Oz+ direction on the point-pattern of Fig. 4. (a) Half cubes, (b) Half spheres.
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.
(a)	(b)
Tortuosity toward Oy into a path (by SQUARE)
,2 45 r
CN
a 40
g 35
I*
ji 30
in
£ 25
re a
a) 20 sz
o 15
Sä 10
0
1	5 .n
ï 0
Tort Y PATH ; AVG: 2.721 ; std dev: 0.036
2.7	2.8	2.9	3	3.1	3.2
Tortuosity
(c)
Fig. 6. Estimation, from propagations of squares, of the tortuosities of the 2D white part of Fig. 3. (a) Forward propagation, (b) Image of the tortuosities (minimum tortuosity: 2.70 in yellow), (c) Histogram of the tortuosities.
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.
(a)
(b)
(c)
Tortuosity along a COMPOSITE TUBE (by CUBE)
90 80 70 60 50 40 30 20 10 0
Tort_TUBE ; AVG: 1.003 ; std_dev: 0.008
" 1	1.01 1.02 1.03 1.04 1.05 1.06
Tortuosity
(d)
Fig. 7. Tortuosity along a 3D discontinuous composite tube processed with cubes. (a) Forward propagation, (b) 2D projection of the forward propagation, (c) 2D projection of the tortuosities (minimum tortuosity: 1 in yellow), (d) Histogram of the tortuosities.
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 Po (T(po) = 0), and f the cost function related to the speed of propagation:
|VT|| = f .
(5)
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 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: v(p) = -¡(pp) = jjpj. 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 fi,j and fi,jk, the cost functions respectively at the pixel (i, j) and at the voxel (i, j, k). We note Tij and Tij,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 6-connectivity neighborhood in 3D of the voxels into the narrow band.
{max{ [Tj - T-u] , [Tij - T+i,j] ,0})2+ 0max {[Tj - T',j-i] , [Tij - Tij+i] , °})2 = fi2
(6)
i,j
°max{ [Tijk - Ti-i,j,k\ , [Ti,j,k - Ti+i,j,k] ,°})2+ (max{ [Ti,jk - Tij-i,k] , [Ti,j,k - TiJ+i,k] ,°})2+
°max{ [Ti,j,k - Ti,j,k-i] , [Ti,j,k - Ti,j,k+i] ,°})2 = j
(7)
t; D fc l i
r f y
JS
4 if
(a)	(b)
Tortuosity toward Oy into a path (by DISC)
« 30
CM
<D
^ 25 in
20
s th
15
v th
£ 10
w
§ 5
D th
ii 0
Tort	_Y_PATH	AVG: 3.04	2 ; std_de ♦	v: 0.040	♦
					
					
					
			♦		
					
3.2
^ 2.7	2.8	2.9	3	3.1
Tortuosity
(c)
Fig. 8. Estimation, from propagations of discs, of the tortuosities of the 2D white part of Fig. 3. (a) Forward propagation, (b) Image of the tortuosities (minimum tortuosity: 3.01 in yellow), (c) Histogram of the tortuosities.
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.
(a)
(b)
(c)
Tortuosity along a COMPOSITE TUBE (by SPHERE)
" 1	1.01 1.02 1.03 1.04 1.05 1.06
Tortuosity
(d)
Fig. 9. Tortuosity along a 3D discontinuous composite tube processed with spheres. (a) Forward propagation along the tube in 3D, (b) 2D projection of the forward propagation, (c) 2D projection of the tortuosities (minimum tortuosity: 1 in yellow), (d) Histogram of the tortuosities.
0
A
as
Direction of
Forward propagation
A
te
I
Direction of Forward propagation
B
B
(a)	(b)
Fig. 10. Influence of isotropic and half discs on the paths in X1. (a) Isotropic discs, (b) Half discs.
(a)	(b)
Fig. 11. Influence of isotropic and half discs on the paths in X2. (a) Isotropic discs, (b) Half discs.
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.
(max{ [Ti,j - Ti+i,j] ,0})2+ {max {[T,j - Tj-i] , [Ti,j - T',j+i] , °})2 = f2
(8)
i,j
(max{ [Ti,j,k - Ti+ij,k] ,0})2+
(max{ [Tijk - Ti,j-i,k] , [T,j,k - T,j+i,k] ,0})2+
(max{ [Ti,j,k - Ti,j,k-i] , [Ti,j,k - Ti,j,k+1] ,^)2 = ff,j,k
(9)
Let XI and X2 be two sets respectively illustrated in Fig. 10 and 11. The set XI 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 XI (Fig. 10b and 11b) leads to non-percolating 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.
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.
INPUT images: ImMark, ImSpeed OUTPUT images: ImSrc4C_6C, ImDistFM
1.1	Initialization()
1.2	{
1.3	Declaration of the Priority Queue FM-PQ;
1.4	Allocate and initialize all the voxels v of ImLabels at the value FAR;
1.5	Initialize all the voxels v of ImDistFM at
1.6	Initialize all the voxels v of ImSrc4C_6C at 0;
1.7	foreach voxel v of ImMark do
1.8	if ImMark[v] = Value_BG then
1.9	ImLabels[v] ^ ALIVE;
1.10	ImDistFM[v] ^ 0;
1.11	foreach voxel vn of the 6C-Neighborhood of v do
1.12	ImSrc4C_6C[vn] ^ RP[v, vn]+1;
1.13	if [(ImLabels[vn] = FAR )& (ImMark[vn] = Value_BG)] then
o ir 7 /ImSpeed[vn\s
1.14	Speed[vn] ^ - vaLMAX
1.15	f[vn] ^ Speed[vn\ ;
1.16	ImLabels[vn] ^ TRIAL;
1.17	Add vn into FM_PQ with the priority f[vn];
1.18	ImDistFM[vn] ^ f[vn];
1.19	end
1.20	end
1.21	end
1.22	end
1.23
2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8 2.9
2.10
2.11 2.12 2.13
2.14
2.15
2.16
2.17
2.18
2.19
2.20
2.21 2.22
2.23
2.24
2.25
2.26
2.27
2.28
FMLoop()
{
while FMPQ is NOT_EMPTY do
foreach voxel v of the 1 st plateau TP of FMPQ do
if ImLabels[v] = ALIVE then
Top .Priority ^ priority of TP[v] into FMPQ;
ImLabels[v] ^ ALIVE; ImDistFM[v] ^ Top_Priority; foreach voxel vn in the 6C-Neighborhood of v do
ImSrc4C_6C[vn] ^ RP[v, vn]+1;
if (ImLabels[vn] = FAR) then ImLabels[vn] ^ TRIAL; Calculate a 1 st estimation of the arrival time T[vn\0 by solving the Eq. 7 from f[vn] and the ALIVE 6C-Neighbors of vn in ImDistFM;
Add vn into FM-PQ with the priority T[vn\0; ImDistFM[vn] ^ T [vn\ 0;
else
if (ImLabels[vn] = TRIAL) then
Update the arrival time Ti,j,k[vn\UPDATE by solving the Eq. 7 from f[vn] and the ALIVE 6C-Neighbors of vn in ImDistFM; Add vn again into FM-PQ with its new priority T[vu\update ; ImDistFM[vn] ^ T [vn\UPDATE;
end
end
end
end
end
Remove the 1st plateau TP from FMPQ;
end
Fig. 13. Iterative loop of the Fast Marching algorithm.
Table 1. Coding square for the 8-connectivity neighbors vn of v in 2D.
5	6	7
3	v	4
0	1	2
Fig. 12. Initialization of the Fast Marching algorithm.
}
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 4-connectivity (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.
			17 {-1,-1,1}	18 {0,-1,1}	19 {1,-1,1}			
			9 {-1,-1,0}	10 {0,-1,0}	11 {1,-1,0}			
			0 {-1,-1,-1}	1 {0,-1,-1}	2 {1,-1,-1}			
Z 17 {-1,-1,1}	9 {-1,-1,0}	0 {-1,-1,-1}	0 {-1,-1,-1}	MM 1 {0,-1,-1}	X 2 {1,-1,-1}	2 {1,-1,-1}	11 {1,-1,0}	19 {1,-1,1}
20 {-1,0,1i	12 {-1,0,0}	3 {-1,0,-1}	3 {-1,0,-1}	4 {0,0,-1}	5 {1,0,-1}	5 {1,0,-1}	13 {1,0,0}	22 {1,0,1}
23 {-1,1,1}	14 {-1,1,0}	6 {-1,1,-1}	Y|6 {-1,1,-1}	7 {0,1,-1}	8 {1,1,-1}	8 {1,1,-1}	16 {1,1,0}	25 {1,1,1}
			6 {-1,1,-1}	7 {0,1,-1}	8 {1,1,-1}			
			14 {-1,1,0}	15 {0,1,0}	16 {1,1,0 }			
			23 {-1,1,1}	24 {0,1,1}	25 {1,1,1}			
			23 {-1,1,1}	24 {0,1,1}	25 {1,1,1}			
			20 {-1,0,1}	21 {0,0,1}	22 {1,0,1}			
			17 {-1,-1,1}	18 {0,-1,1}	19 {1,-1,1}			
Fig. 14. Pattern of the coding cube for the 26-connectivity neighbors vn of v in 3D.
6-connectivity (3D) backtracking
Based on the 26-connectivity cubic neighbor pattern illustrated in Fig. 14, the identification numbers of the relative positions (25 — vn) of the voxels v are determined from v and vn in 3D. Note that the voxel v is supposed to be at the center of the square in 2D (Table 1) and of the cube in 3D (Fig. 14). The "+1" instructions (lines 1.12 and line 2.10) are handled because all the voxels of ImSrc4CJ6C which are not reached by the propagating front are set to zero. Thus, to sum up, the value of each voxel of the image ImSrc4C-6C corresponds to the relative position of its nearest 4-connectivity (2D) or 6-connectivity (3D) neighbor from the marker seeds. More details about this method of source backtracking can be found in (Peyrega, 2010).
8-connectivity (2D) and 26-connectivity (3D) backtracking
Another method can be handled to process a source backtracking in 8-connectivity for 2D images, or in 26-connectivity for 3D images. However, this induces additional processing time. In order to extract the nearest voxel from the marker seed among the 8 neighbors vn of any voxel v into the geodesic mask, it is necessary to scan the full 8-connectivity (2D) or 26-connectivity (3D) neighborhood around v to find the voxel vn with the minimum value into the image ImDistFM of the geodesic distance map (algorithms 12 and 13). Both Table 1 and Fig. 14 are used to fill the corresponding values of the voxels of the image ImSrc8CJ26C of the sources in 8-connectivity (2D) or 26-connectivity (3D).
PATH RECONSTRUCTION
After having defined the source of each voxel of the geodesic mask for both forward and backward propagations of the algorithm estimating the morphological tortuosity, it is then possible to reconnect any components into this mask to both marker faces, even if they have been disconnected by thresholding the tortuosity. The method to process the path reconstruction consists in reconnecting the isolated components to the marker face by using the source images (forward and backward) which indicate the path to follow until the marker faces.
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. 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.
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.
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. 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.
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 etal. 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.
Tortuosity of the FIBRES of Thermisorel
(a)
(b)
Fig. 19. 3D X-Ray CT image of Thermisorel™™. (a)-Gray level image. (b)-Binarized image; (Dimensions: 5.6 x 5.6 x 3.4 mm3. Resolution: 9.36 jjm/voxel. Source: US2B).
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.
Tort_X_FIBRES ; AVG: 1.039 ; std_dev: 0.015 ♦ Tort_Y_FIBRES ; AVG: 1.027 ; std_dev: 0.014 ■ Tort Z FIBRES ; AVG: 1.197 ; std dev: 0.063 •
= 4
t 0
1.4 1.5 1.6 1.7 1.8 1.9 Tortuosity
(a)
Tortuosity of the PORES of Thermisorel
* 3
Ö 2
x o
* 0
1.1	1.15
Tortuosity
(b)
Fig. 20. Histograms of the morphological tortuosities (a) into the fibres and (b) into the pores of ThermisorelM in Fig. 19 in the Ox, Oy and Oz directions.
6
5
3
2
1
6
5
4
1
PATH RECONSTRUCTIONS IN 3D IMAGES OF THERMISOREL
Table 2. Minimal, maximal, and average values, and standard deviation of the morphological tortuosities into the fibres of Thermisorel™ in Fig. 19 in the Ox, Oy and Oz directions.
Direction	Ox	Oy	Oz
TMinimal	1.005	1.002	1.050
TMaximal	1.092	1.084	1.418
^ Average	1.039	1.027	1.197
OT	0.015	0.014	0.063
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 26-connectivity 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.
Table 3. Minimal, maximal, and average values, and standard deviation of the morphological tortuosities into the pores of Thermisorel™ in Fig. 19 in the Ox, Oy and Oz directions.
Direction	Ox	Oy	Oz
TMinimal	1.002	1.003	1.004
TMaximal	1.095	1.060	1.247
TAverage	1.028	1.022	1.082
OT	0.012	0.010	0.039
In the fibres
(c) TMax = 1.092
Fig. 21. Reconstruction of the paths by backtracking in 26-connectivity after thresholding the morphological tortuosity in the Ox direction for the fibres of Thermisorel™.
(c) TMax = 1.084
Fig. 22. Reconstruction of the paths by backtracking in 26-connectivity after thresholding the morphological tortuosity in the Oy direction for the fibres of Thermisorel™.
(c) TMax = 1.418
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™.
In the pores
(c) TMax = 1.095
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™.
(a) TMin = 1.003
(a) TMin = 1.004
(b) TAverage = 1 .°22
(b) TAverage = 1.082
(c) TMax = 1.060
(c) TMax = 1.247
Fig. 25. Reconstruction of the paths by backtracking in	Fig. 26. Reconstruction of the paths by backtracking in
26-connectivity after thresholding the morphological	26-connectivity after thresholding the morphological
tortuosity in the Oy direction for the pores of	tortuosity in the Oz direction for the pores of
Thermisorel™M.	Thermisorel™.
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.
ACKNOWLEDGEMENTS
We gratefully acknowledge the Agence Nationale de la Recherche (ANR) for financially supporting the Silent Wall project (Ref.: ANR-06-MAPR-00-18) in which this study takes place.
REFERENCES
Allard J (1993). Propagation of sound in porous media. Modelling sound absorbing materials. London: Chapman & Hall.
Allard J, Castagnede B, Henry M, Lauriks W (1994). Evaluation of tortuosity in acoustic porous materials saturated by air. Rev Sci Inst 65:754-5.
Borgefors G (1986). Distance transformations in digital images. Comput Vision Graph 34(3):344-71.
Coeurjolly D, Miguet S, Tougne L (2004). 2D and 3D visibility in discrete geometry: an application to discrete geodesic paths. Pat Rec Let 25:561-70.
Cohen L, Kimmel R (1997). Global minimum for active contour models: a minimal path approach. Int J Comput Vis 24(1):57-78.
Decker L, Jeulin D, Tovena I (1998). 3D morphological analysis of the connectivity of a porous medium. Acta Stereol 17(1):107-12.
Lantuejoul C, Beucher S (1981). On the use of the geodesic metric in image analysis. J Microsc 121(1):39-49.
Petres C, Pailhas Y, Evans J, Petillot Y, Lane D (2005). Underwater path planning using fast marching algorithms. In: Proc IEEE Oceans 2005 Eur Conf. June 20-23, 2005. Brest, France. Vol 2.
Peyrega C (2010). Prediction of the acoustic properties of heterogeneous fibrous materials from their 3D microstructures. PhD thesis, Mines ParisTech.
Peyrega C, Jeulin D, Delisee C, Malvestio J (2009). 3D morphological modelling of a random fibrous network. Image Anal Stereol 28:129-41.
Peyrega C, Jeulin D, Delisee C, Malvestio J (2010). 3D morphological characterization of phonic insulation fibrous media. Adv Eng Mat 13(3):156-64.
Serra J (1988). Image analysis and mathematical morphology vol 2: Theoretical advances. London: Academic Press.
Sethian J (1996). Level set methods. Evolving interfaces in geometry, fluid mechanics, computer vision, and materials science. Cambridge: Cambridge University Press. pp 81-95.
Soille P (1992). Morphologie mathematique: du relief a la dimensionalite-Algorithmes et methodes. PhD thesis, Universite Catholique de Louvain; in collaboration with the Ecole des Mines of Paris.
Soille P (1994). Generalized geodesy via geodesic time. Pattern Recogn Lett 15(12):1235-40.
Soille P (2003). Morphological image analysis. Principles and applications. 2nd Ed. Berlin: Springer-Verlag. pp 219-40.