Tortuosimetric operator for complex porous media characterization

Geometric tortuosity is one of the foremost topological characteristics of porous media. Despite the various deﬁnitions in the literature, to our knowledge, they are all linked to an arbitrary propagation direction. This paper proposes a novel topological descriptor, named M-tortuosity , by giving a more straightforward deﬁnition, describing the data regardless of physicochemical processes. M-tortuosity , based on the concept of geometric tortuosity, is a scalable descriptor, meaning that information of several dimensions (scalar, histograms, 3D maps) is available. It is applicable on complex disconnected structures without any arbitrary deﬁnition of entry and exit. Topological information can be represented by aggregation into a unique scalar descriptor for classiﬁcation purposes. It is extended by iterative erosions to take into account porous structure narrowness, especially bottleneck effects. This new descriptor, called M-tortuosity-by-iterative-erosions , describes tortuosity of the porous part as seen by a spherical particle of given size walking along the network. Boolean models are used to simulate different porous media structures in order to test the proposed characterization.


INTRODUCTION
Analysis and characterization of the topology of biphasic materials, in particular porous media (formed by solid and voids), is of paramount importance.As their usage properties depend on their microstructure, intuitive features help to design the porous phase topology targeting a specific application.Their complex interconnected capillary network renders their accurate description, a difficult task.Statistical modeling of materials by random morphological models (Matheron , 1975) is an essential tool for behavior prediction.In particular, Boolean models (Serra , 1982) and multi-scale Cox Boolean models (Jeulin , 1996;2010) are considered in this paper.Image processing operators quantify porous media volumes by geometrical and topological features, as specific descriptors.This paper addresses a topological description based on the concept of tortuosity, first introduced by Carman (1937) in the study of permeability.The term "tortuosity" is polysemic (Clennell , 1997), even in the recent literature (Ghanbarian et al., 2013a).Therefore, it is important to accurately define which tortuosity is considered to avoid misunderstanding.
The focus here is on the geometric approach of tortuosity characterization, quantifying porous volume sinuosity (Clennell , 1997), without taking into account any physicochemical phenomenon.More precisely a specific definition is used; geometric tortuosity is defined between two points as the ratio of their geodesic and Euclidean distances (Lantuéjoul and Beucher , 1981;Decker et al., 1998).Despite the broad range of definitions, up to our knowledge, geometric tortuosity has always been defined with respect to a certain propagation direction.This can be explained by the effort aiming to connect it with percolation theory (Ghanbarian et al., 2013b).Theoretical definitions from basic modelizations of porous structures (Yu and Li , 2004;Yun et al., 2006) exist.In Lindquist et al. (1996), the probability density function of the geometric tortuosity is fit by a gamma distribution for specific porous media.Some original definitions of morphological tortuosity (Decker et al., 1998), which stands for geometric tortuosity, have also been proposed.Morphological dilations, with different structuring elements, are used for tortuosity computation in Peyrega and Jeulin (2013).Gommes et al. (2009) use the morphological recontruction operator to assess the tortuosity of a porous volume.

INTRODUCTION
Analysis and characterization of the topology of biphasic materials, in particular porous media (formed by solid and voids), is of paramount importance.As their usage properties depend on their microstructure, intuitive features help to design the porous phase topology targeting a specific application.Their complex interconnected capillary network renders their accurate description, a difficult task.Statistical modeling of materials by random morphological models (Matheron , 1975) is an essential tool for behavior prediction.In particular, Boolean models (Serra , 1982) and multi-scale Cox Boolean models (Jeulin , 1996;2010) are considered in this paper.Image processing operators quantify porous media volumes by geometrical and topological features, as specific descriptors.This paper addresses a topological description based on the concept of tortuosity, first introduced by Carman (1937) in the study of permeability.The term "tortuosity" is polysemic (Clennell , 1997), even in the recent literature (Ghanbarian et al., 2013a).Therefore, it is important to accurately define which tortuosity is considered to avoid misunderstanding.
The focus here is on the geometric approach of tortuosity characterization, quantifying porous volume sinuosity (Clennell , 1997), without taking into account any physicochemical phenomenon.More precisely a specific definition is used; geometric tortuosity is defined between two points as the ratio of their geodesic and Euclidean distances (Lantuéjoul and Beucher , 1981;Decker et al., 1998).Despite the broad range of definitions, up to our knowledge, geometric tortuosity has always been defined with respect to a certain propagation direction.This can be explained by the effort aiming to connect it with percolation theory (Ghanbarian et al., 2013b).Theoretical definitions from basic modelizations of porous structures (Yu and Li , 2004;Yun et al., 2006) exist.In Lindquist et al. (1996), the probability density function of the geometric tortuosity is fit by a gamma distribution for specific porous media.Some original definitions of morphological tortuosity (Decker et al., 1998), which stands for geometric tortuosity, have also been proposed.Morphological dilations, with different structuring elements, are used for tortuosity computation in Peyrega and Jeulin (2013).Gommes et al. (2009) use the morphological recontruction operator to assess the tortuosity of a porous volume.

1
Notwithstanding that tortuosity is defined with respect to a path, our goal is to look for a single scalar value, gathering topological information, able to describe the whole structure of a porous medium, for discrimination purposes.This is motivated by the rationale that the flow direction inside a porous solid is non predictable in practice.Descriptor scalability, information of different dimensions (scalar, histograms, 3D maps), is necessary for a complete characterization, since a unique scalar is not sufficient.The resulting operator named M-tortuosity checks all criteria.An extension is proposed, named Mtortuosity-by-iterative-erosions, probing the porous structure with a spherical percolating particle with given radius.The first part of this paper deals with theoretical definitions, determinist and estimator definitions, and M-tortuosity estimator properties.Numerical methods are then explained in detail.Finally, the results obtained on toy cases and Boolean models are exposed and discussed.

METHODS
Notations used for M-tortuosity descriptor are first presented.Then theoretical foundations for the different operators are presented.

Notations
We base the notations on the notation system proposed by Criminisi et al. (2010) and Lohou and Bertrand (2005).Let I be a binary function, I : R 3 → {0, 1}.Feature points are defined by the set Let ∂ I be the convex hull of X, with I a bounded subset of R 3 , defined as the smallest convex subset such that X ⊂ I .Let I be a 3D binary image defined by I and I , I : I → {0, 1}.In other words, I is the limitation of is the 3D geodesic distance map and ∀x ∈ X, D G (x, S; X) is the geodesic distance transform value at the point x from S restricted to X.

Geodesic Distance Transform
The geodesic distance transform is the equivalent of the distance transform (Rosenfeld and Pfaltz , 1968;Borgefors , 1986) in a non-Euclidean space, typically the distance transform on gray-level images.The grayvalues are seen as "heights", a functional point of view introduced by Rutovitz (1968); the distance is propagated on the resulting surface.Geodesy could also be seen as the restriction of the distance transform to a subset of the image.This is the viewpoint taken in this work.The geodesic distance transform D G (., S; X) (Criminisi et al., 2010) of each voxel x ∈ X, from the set S, restricted to the support X is defined as: where γ x,y;X is the set of all possible paths in R 3 constrained by X, between y ∈ X (starting point) and x ∈ X (ending point).Γ is one of these paths, and s ∈ [0, L(Γ)] its curvilinear abscissa, with L(Γ) the length of the path Γ. D G (x, y; X) (S = {y}) "is the greatest lower bound of the lengths of the arcs in X ending at points x and y, if such arcs exist, and +∞ if not" (Lantuéjoul and Beucher , 1981).

Geometric tortuosity
For each (x, y) ∈ X 2 , such that x = y and x = c, the geometric tortuosity between x and y, named τ x,y , is defined as the ratio of the geodesic distance (length of the shortest path) and the Euclidean distance (length of the straight path) between these two points (Decker et al., 1998).
with D(y, x) the Euclidean distance between x, the starting point, and y, the ending point.

M-coefficient
Given x ∈ X, x = c, the M-coefficient associated to x, C x , is defined as the arithmetic mean of the geometric tortuosities, τ x,y , where the respective weights are the geodesic distances D G (y, x; X).
The idea behind the weighting using geodesic distances (Berrocal et al., 2016) is that, the longer is the path the more representative it is of the overall porous structure tortuosity.

M-scalar
The M-scalar, τ M , is defined as the arithmetic mean of {C x } x∈X\{c} , where the respective weights are the Euclidean distances from c, D(x, c).
The idea behind this second weighting is close to the previous one: the more a point is away from the center of mass of the network, the higher the probability to get long paths in its corresponding M-coefficient value.
At this step, such a formulation does not take into account the case of multiple connected components.It would lead to consider each of them separately.In fact, it is possible to process simultaneously all the connected components using the following formulation.

Redefinition of M-coefficient and M-scalar
With this perspective, we replace Eq. 4 with Eq. 6 and Eq. 5 with Eq. 7. The M-coefficient associated to x is defined, for a given x ∈ X, x = c, using the harmonic mean of the geometric tortuosities τ x,y , weighted by the inverse of the respective geodesic distance D G (y, x; X).
The harmonic mean is defined as the reciprocal of the arithmetic mean of the reciprocals of the values, here the geometric tortuosities.The denominator is equal to zero only if x is isolated, i.e. connected to no point of X.In this case, we impose C x = 0.For the same purpose, τ M , the M-scalar, is defined as the harmonic mean of {C −1 x } x∈X\{c} , weighted by the inverse of the respective Euclidean distance from c, D(x, c).
This formulation handles the disconnections corresponding to an infinite geodesic distance.The integral of the denominator is equal to zero only if each point of X \ {c} is disconnected from the others.
Unfortunately, this deterministic way to define the M-tortuosity descriptor, T M , is difficult to apply on large volumes in practice.Therefore, an estimator of the M-tortuosity, TM , is necessary, and this will be done using a skeletonization step and a point sampling approach.The next section introduces the various steps which define TM , the M-tortuosity estimator, and its extension, the M-tortuosity-by-iterative-erosions estimator.

M-TORTUOSITY ESTIMATOR Skeletonization
Most numerical methods have to balance accuracy and computation time efficiency.Skeletonization, used as a pre-processing of X, allows first to decrease the computation time.A skeleton can be defined in several ways (Tagliasacchi et al., 2016).In this paper, the thinning method of Lohou and Bertrand (2005), using the notion of P-simple point, is chosen.The obtained skeleton, named Sk, is the smallest homotopic subset of X. Obviously, skeleton computation biases the results, but this bias is intentional.Let us recall that we attempt to define a global geometric tortuosity, for discrimination of complex networks.Fig. 1 displays a basic situation.Intuitively, we would like a significant difference of geometric tortuosity between the straight pore and the sinuous pore in Fig. 1.Computing the geometric tortuosity on the skeleton enhances the difference between the two pores of Fig. 1.In the section "Explanation on the use of Monte Carlo method on the skeleton", a specific case illustrates this assumption by comparing the distribution of distances of both cases, with and without skeletonization.Skeletonization provides a less scattered distribution of distances too (Fig. 2).Moreover, homotopic skeleton allows to have a stability with respect to morphological erosion, defined below.This characteristic matters, as we would like an equal geometric tortuosity coefficient for two sinuous pores of identical shape, but with distinct narrowness.This is of particular interest for the M-tortuosity-byiterative-erosions operator.For all this reasons, M-tortuosity operator TM is defined on Sk, and not on X.

Sampling 1D sampling:
This sampling strategy corresponds to a uniform sampling in the pore phase skeleton.Sk can be rewritten as Sk = {x i } i∈ [[0,] , which can be seen as an indexing of Sk, with |Sk| its number of elements.Uniform distribution U([0, |Sk| − 1]) allows to randomly draw N distinct points in Sk \ {c} defining S.

Stratified sampling:
This sampling strategy corresponds to a uniform density of sampled points in Sk (Neyman , 1934) and 3 corresponds to the dimension of the space.If Sk ∩ W i = / 0, as before, one point is uniformly drawn in Sk ∩W i \{c}.N points are sampled in Sk \ {c}, with N ≤ K, defining S.

Geometric tortuosity
Given S, the geometric tortuosity (Decker et al., 1998) between p n ∈ S and p m ∈ S, m = n, named τ n,m , is defined as with D(p m , p n ) the Euclidean distance between p n , starting point, and p m , ending point.

M-coefficient
Given More details about the weighting are given in the "Numerical Methods" section ("Geodesic Distance Transform").

M-scalar
The estimated M-scalar, τM , is defined as The M-tortuosity descriptor, TM , applied to a porous network skeleton Sk, is defined by the previous steps; skeletonization, sampling, geodesic distance transforms, Euclidean distance transforms, geometric tortuosities, M-coefficients and M-scalar.A numerical definition of this descriptor is given below.

Morphological erosion
Let ε r (X) be the eroded set (Serra , 1982) of X by a sphere of radius r ∈ N defined as with D(x, X c ) the value of the distance transform (Rosenfeld and Pfaltz , 1968;Borgefors , 1986), D(., X c ), to the complementary set of X, at point x. with, where γ x,y is the set of all possible paths in R 3 , between y (starting point) and x (ending point).Γ is one of these paths, and s ∈ [0, L(Γ)] its curvilinear abscissa, with L(Γ) the length of the path Γ.
Morphological erosion is used for defining the M-tortuosity-by-iterative-erosions estimator, TM,r , assessing tortuosity for a percolating particle of a given size.
In the next section, a realization of a Boolean model is used to illustrate the use of homotopic skeletonization.We show that the conditions of the use of Monte Carlo method are fulfilled, justifying the definition of TM .

EXPLANATION ON THE USE OF MONTE CARLO METHOD ON THE SKELETON
A realization of a Boolean model of size 300 3 , with volume fraction V v = 0.7 of spheres with constant radius r = 10, is generated to illustrate our purpose (see section "Results" for details about Boolean models).
The following calculations are done on the spheres network and on the skeleton of this network, for comparison purposes.Let us consider two locations belonging to the pore phase skeleton of the given microstructure.The two values d and g, their Euclidean distance and their geodesic distance, respectively, can be represented as a point (d, g) in a coordinate system.Doing that for several point pairs (see example in Fig. 2), we obtain a kind of directional "jet" above neutral line d = g, corresponding to nearby locations.The results with and without skeletonization step are presented for exactly the same point pairs (Fig. 2).Differences are highlighted for 3 point pairs (black spheres, squares and triangles in Fig. 2).Skeletonization allows to increase the geometric tortuosity value, and especially allows to increase the difference between two microstructures, by being more sensitive to small variations.Clearly, the more tortuous is the network, the steepest will be this cloud of points; a first candidate for being a "coefficient of tortuosity" could be the average of the different slopes g / d.In fact, as we have seen above, it is more interesting to work with the inverse ratios d / g.The problem is that these calculations rely on samplings.Will coefficients obtained with two different samplings be reasonably close?Instead of considering an underlying "abstract" bivariate density, that could serve as a reference on which various samplings are operated, we prefered to operate on simulations with a very good degree of reproductibility, as we are going to see it.Fig. 2 (b) displays a cloud of points which can be seen as a probability density function.A digital operator, based on Monte Carlo method, can be a good estimator of such a density, if and only if some conditions are fulfilled (Caflisch , 1998).First, we are in a bounded set, Sk.Secondly, our estimator can be integrated on this subset if and only if at least two points are connected.Practically, we generate 100 realizations of a Boolean schemes of size 200 3 and of volume fraction V v = 0.7 of spheres of constant radius r = 3.We compute the homotopic skeleton (Lohou and Bertrand , 2005) of the complementary set of the set of spheres, for prooving the decreasing of the variance as 1/N, with N the number of sampled points, using 1Dsampling.Fig. 3 displays the variation of Nσ 2 , with σ 2 the variance, as a function of N. The stabilisation of this curve validates the fact that we are in an asymptotic domain and the convergence is reached.

M-TORTUOSITY DESCRIPTOR PROPERTIES
M-tortuosity estimator, as a topological descriptor, should verify some invariance properties.Due to its definition, the geodesic distance transform is already invariant by translation and rotation.In this section, M-tortuosity estimator invariance properties and stability on a periodic structure, are established.

Translational invariance
Let T be the translation defined by vector t ∈ R 3 .
Eq. 9 and 14 allow us to say that the M-coefficients are translation invariant.As Euclidean distance is translation invariant itself, τ(t) M is equal to τM .Finally the M-tortuosity descriptor, TM , is translation invariant.

Rotational invariance
Let R be the rotation transform with Euler angles (θ , φ ) defined by matrix M θ ,φ ∈ M 3,3 (R), Let I θ ,φ be the rotated version of angles (θ , φ ) of I.It follows, S θ ,φ and Sk θ ,φ , representing the rotated sets of S and Sk, respectively.Due to its definition, the geodesic distance transform, D G (., S; Sk), is invariant with respect to rotation.As before, there is no difficulty proving that M-tortuosity descriptor, TM , is rotation invariant.

Homothety invariance
Let H be the homothety transform defined by its center O ∈ R 3 and its ratio λ ∈ R * .
Let I O,λ be the homothetic set of I by an homothety of center O and ratio λ .Let S O,λ and Sk O,λ be the homothetic sets of S and Sk, respectively.Whatever Euclidean distance respects this relation too.Then the geometric tortuosity between x O,λ and y O,λ , (16) Eq. 9 and 16 allow us to say that the M-coefficients are homothety invariant.Then the M-scalar, Finally TM is homothety invariant.

Periodic pattern repetition stability
Periodic pattern is here defined by the function I 0 .The image I 0 , defined by I 0 and I 0 , with X 0 the feature points set, is seen in a more practical viewpoint; the convex hull of X 0 , ∂ I 0 , is a cube, as shown Fig. 4 and  5.This stability is of paramount interest for crystalline structures characterization.
Proposition 1.Let I 0 be a periodic function in all three directions.I 0 is named here the unit cell.(T x , T y , T z ) ∈ N 3 are the respective periods which are the size of I 0 too.Let I be a repetition of q times I 0 in the three space directions of R 3 .Then, if q tends to infinity and N is high enough, the estimated M-scalar, τM , will converge.In other words, there exists a τ ∈ R such that The proof of this proposition of stability of TM uses the fact that the integrality of the topological information is contained in the unit cell, I 0 .For didactic purposes, I is displayed here as a concatenation of 3 times I 0 in the x direction only (Fig. 4).For these same purposes, illustrations are in 2D (Fig. 4 and 5) and calculations are done in a 2D plan of R 3 .Moreover, we will focus on very specific points, as shown below, for illustration purposes.The results can be extended to any I, defined as a repetition of q times I 0 in every space directions of R 3 and for any points pair.Let (p, p 0 ) ∈ Sk 2 0 and p 1 and p 2 , two points of Sk, be defined by Let p k 0 ∈ Sk 0 ∩ ∂ I 0 be defined by Redefining geodesic distance transform in the specific case of a periodic image gives with p l 0 ∈ Sk 0 ∩ ∂ I 0 , as shown in Fig. 4.
Let us now consider the geodesic distance between p and p 2 .Using the same scheme we obtain: Fig. 4: Sketch of I, concatenation of 3 times I 0 (first "square"), the black broken lines represent the porous network skeleton.Possible positions of p and p 0 , and the corresponding positions of p 1 , p 2 , p k 0 and p l 0 .Geodesic paths between p and p 1 and between p and p 2 , are represented (red broken lines).
Euclidean distance is needed too, for geometric tortuosity definition.Considering (x, y) ∈ R 2 the orthogonal projections of pp 0 , as represented in Fig. 5, the Euclidean distances between p and p 1 and between p and p 2 can be written respectively, (23) Fig. 5: Sketch of I, concatenation of 3 times I 0 (first "square"), and representation of some possible positions of p and p 0 , and the corresponding positions of p 1 , p 2 , p k 0 and p l 0 .Euclidean distances between p and p 1 and between p and p 2 , are represented with their orthogonal projections on x and y directions using (x, y) ∈ R 2 and T x ∈ R.
Let us generalize these concepts.I is now a concatenation of q times I 0 in the x direction, and p q ∈ X the point defined by The geodesic distance and Euclidean distance between p and p q are respectively, D G (p q , p; Sk) = D G (p k 0 , p; Sk) + (q − 1).D G (p k 0 , p l 0 ; Sk) + D G (p 0 , p l 0 ; Sk) D(p q , p) = y 2 + (x + q.T x ) 2 . (25) We have the following limits, Therefore τ q , the geometric tortuosity between p and p q , is such that, This means that for q high enough, using the previous equivalences, Ĉn,q can be expressed using the geometric tortuosities of each percolating path of I 0 , the unit cell, here in the x direction.
where A i, j,q ∈ N is a coefficient depending on the ending point p q and on the percolating path (p l i , p k j ), and p l i and p k j belonging to two different faces of I 0 , in this case two opposite ones.Let I and I be two images generated by concatenation of q and q−1 times I 0 , respectively.Let C n,q and C n,q−1 be two estimated M-coefficients connected to the same starting points p and the same set of end points.Using the previous approximations, lim q→∞ Ĉn,q−1 Ĉn,q = 1.( 29) There exists a C ∈ R * such that As said before, all these calculations are extended to any image I defined as a concatenation of q times the unit cell I 0 in each space direction of R 3 and for any points pair.Then using Eq.30 and 25, for N high enough, there exists a τ ∈ R such that meaning that TM is stable by periodic pattern repetition.In practice, N is chosen with respect to the variance of the estimator (cf., Fig. 3).Summary: the M-tortuosity descriptor is invariant for all desirable transformations and stable for periodic structures.The part describes the numerical methods of computation of M-tortuosity and Mtortuosity-by-iterative-erosions estimators, to assess the defined characteristics from digital images.

NUMERICAL METHODS
M-tortuosity and M-tortuosity-by-iterativeerosions can be described in a discrete space.The descriptors are here the numerical version of the estimators TM and TM,r (r ∈ N) defined below (i.e., the algorithms or computational methods).Two main aspects of our methods are presented too: connectivity, being a central issue of any topological analysis, and geodesic distance transform, on which our descriptors are based.

M-tortuosity
M-tortuosity has a certain sampling as its N random points are sampled in Sk, using one of the two previous sampling methods, defining the set S = {p i } i∈ [[0, N−1]] .The geodesic distance transform, D G (., p i ; Sk), is computed for each point of S; it means that each p i will be, only once, the unique starting point for the geodesic distance map computation.In the following algorithms (algo. 1 and 3), D G (., p i ; Sk) and D(., p i ) will be written D G (p i ) and D(p i ) for brevity.Then given p n , the starting point, n ∈ [[0, N − 1]], the geometric tortuosity is computed for each p m , m ∈ [[0, N − 1]], m = n, by using Eq. 8. Then the N − 1 values τ n,m allow the computation of Ĉn , the estimated M-coefficient associated to p n (Eq.9).These operations are iterated N times on S, in order to get N M-coefficients, each one corresponding to a point of S. Finally, the M-scalar, τM , is computed using Eq.10.This computation process defines the estimator TM , called M-tortuosity descriptor.The following pseudocode (algo.1) gives a detailed description of the whole computation process.Details about M-tortuosity main step, geodesic distance transform, are given below.Sketches are done in 2D for didactic purposes.

M-tortuosity-by-iterative-erosions
Extension of TM , by using iterative erosions of X, allows to take into account local narrowness of X, and especially the bottleneck effect.Mtortuosity-by-iterative-erosions is then connected to the constrictivity concept (Petersen (1958); Holzer et al. (2013) cited in Neumann et al. (2018)).The idea behind this new descriptor, TM,r , is to assess the tortuosity of the visible porous part for a given probe size, the radius r of the spherical percolating particle used for erosions.The principle is the following; the M-tortuosity descriptor is applied on ε r (X), the eroded set of X by a sphere of radius r (Eq.11).TM,r is defined for any integer r by, TM,r (X) = TM (ε r (X)). (32) In the following pseudo-code (Algo.2), D(X) stands for the distance map to the complementary set D(., X c ) (Eq. 12) for all x ∈ X, for brevity.The stopping condition could well be defined using the number of connected points in S, or the number of connected components, or the volume of ε r (X) etc.

Algorithm 2: M-tortuosity-by-iterativeerosions
Result: M-scalar vector τM,r Computation of distance map D(X); while Condition is false do M-tortuosity computation: TM (X); Erosion by sphere of radius r: X = ε r (X); Condition computation; r = r + 1; end

Connectivity
The connectivity degree, as it has been considered, has a non-negligible impact on digital image processing computation methods.It is well known that 3D connectivity considers either 6, 18 or 26 neighbor voxels (sharing either a face, at least an edge or at least a vertex only).Using the definitions from (Lohou and Bertrand , 2005), let x ∈ Z 3 defined by (x 1 , x 2 , x 3 ) be the current point, three neighborhoods can be defined: The connected components are then directly defined by the choice of the foreground and background neighborhoods.In this paper, we take the most usual choice: N 26neighborhood for the foreground and N 6 -neighborhood for the background.

Geodesic Distance Transform
The geodesic distance transform is a very powerful tool for connectivity issues (Lantuéjoul and Beucher , 1981) and is at the very basis of our descriptors.Fig. 6 shows the difference between Euclidean distance transform on X and geodesic distance transform for the same set X and starting point s.X has two connected components and the second one, the disk, is not reached by the geodesic propagation from s unlike the Euclidean propagation.Moreover, the propagation itself is completely different too.Fig. 6: Differences between Euclidean distance transform (left) on X and geodesic distance transform (right) restricted to X. s is the starting point used for both maps.
The purpose of the weighting of Eq. 9, as said above, is to promote long paths which are more representative of the porous structure for a global tortuosity assessment.Moreover, the inverse of the geometric tortuosity and the respective geodesic distance allows to deal with disconnected components.This comes from a geodesic distance transform characteristic; the geodesic distance of two nonconnected points of X (i.e., belonging to two distinct connected components, as in Fig. 6) is infinite.Indeed, in such a case, there is no path, totally included in X, connecting the two points.Using the inverse of the geodesic distance annihilates the contribution of such paths.Then, given n ∈ N, Ĉn , the M-coefficient associated to the starting point p n , will take into account only paths belonging to the same connected component of p n ; nonconnected points will not interfere in the computation.Raster-scanning algorithm, used for geodesic distance computation, is presented next.

Raster-scanning algorithm:
This algorithm iterates image scans until stabilization (i.e., idempotence).A similar algorithm can be found in Toivanen (1996) for grayscale images.An iteration is composed of two raster scans; the forward scan -from top to bottom, and from left to right-and the backward scan -from bottom to top, and from right to left-.For each scan a specific mask is used, as shown in Fig. 7. First, voxels belonging to S, the starting voxels, are initialized to 0 and all the others to infinity.The forward scan starts from the top left voxel.Computation and update, if necessary, of the geodesic distance value of each voxel of X is done if with x the current voxel (gray pixel in Fig. 7), n i stands for a neighbor (colored pixels in Fig. 7) and W (n i ) is the weight corresponding to the neighbor.The backward scan is done in a very same way with an appropriate neighborhood.These scans are iterated until idempotence is reached, which means that no update has been done in D G at the final iteration.We have gathered these explanations in the pseudo-code (algo.3).
Fig. 7: The two scans used for the raster method, the starting voxel is represented by a black point.

RESULTS AND DISCUSSION
M-tortuosity describes the sinuosity of capillary networks.Information at various dimensions can be extracted from this scalable topological descriptor.For instance, distributions of geometric tortuosities τ n,m and M-coefficients Ĉn , characterize the porous structure.Heterogeneity can be described by scatterplots of geometric tortuosities or geodesic distances, at the points of S, according to Euclidean distances.3D geometric tortuosity maps can also be extracted from M-tortuosity estimator and its extension.In this paper, the focus is mainly on τM,r , for classification purposes.Before proceeding to our results, we are going to present some numerical details and choices about the descriptors.

PRACTICAL DETAILS ABOUT NUMERICAL ESTIMATORS Sampling
As said above, two sampling methods are considered, 1D-sampling and stratified sampling.Fig. 8 exhibits the difference between the two sampling methods.Due to its faster convergence (Baddeley and Jensen , 2005), stratified sampling has been chosen for the results exposed below.The targeted number of points in Sk is 50 then k = 4 in our case.This choice is motivated by the convergence study (Fig. 3).Although 10 Fig. 11: τM values according to the index of (φ 2 lim,down , φ 2 lim,up ) for the P 1 case.Confidence intervals, with confidence level at 95 % (see definition below), are represented by vertical bars.
As expected, τM increases when limited angles increase.The generation of the second toy case, P 2 , is based on 3D-crosses made of 3 pores.The same limiting parameters are used for pores' generation.Each pore of a cross (see Fig. 12) starts and ends at the center of the two opposite faces perpendicular to its propagation direction, and passes through the center of the cube I of size 300 3 voxels.P 2 is composed of 8 cubes or images I; 4 of them defined by φ 1 , θ 1 and L 1 , the other 4 cubes by φ 2 , θ 2 and L 2 .The 8 cubes are organized following a simple rule; a cube cannot have a common face with an identical cube.An example is given in Fig. 12.The number of points N is between 42 and 57; N increases with tortuosity.The behavior of τM corresponds to the expectations for both cases (P 1 and P 2 ).The only difference between the two cases is that the values taken by τM are higher in the second case, P 2 .This comes from the fundamental structure of the "porous networks", in other words, their lowest scale.For P 1 , this is simply a straight line whereas it is 8 3D-crosses for P 2 .

Boolean models
Boolean models (Matheron , 1975;Serra , 1982) are based on a Poisson Point process of intensity θ (Kingman , 1993).The number of points to be placed, U, is a random variable, following a Poisson distribution of parameter θV (W ), U ∼ Poi(θV (W )), W being the bounded domain and V (W ) its volume.Random primary grains A (overlapping allowed) are located at Poisson points x k .A and intensity θ define Boolean model A (Eq. 35).
Eq. 36 shows dependency of θ with respect to the average volume V (A ) of the primary grains of the material A and on the volume fraction V v . 1 For more details about Boolean models, please refer to Matheron (1975); Serra (1982); Chiu et al. (2013).Image Anal Stereol 2019;38:25-41 Image Anal Stereol ?? (Please use \volume):1-17

Multi-scale Cox Boolean models
A multi-scale microstructure can be modelized by using multi-scale Cox Boolean model (Jeulin , 1996).This process is defined by intersections and unions between objects and points, generated by several Poisson point processes, named Cox point processes (Jeulin , 1997).This way to model multi-scale complex microstructures gives more realistic results than intersections and unions of several Boolean models (Savary et al., 1999), in particular no grains are cut.This difference is illustrated in Fig. 14.  (Moreaud et al., 2018).
In this paper, the focus is mostly on threescale models defined by three volume fractions; V v,inc volume fraction of inclusion areas (defining aggregates), V v volume fraction of grains inside inclusion areas and V v,out volume fraction of grains outside of inclusion areas.Exclusion zones, free of any grain, can be used for the characterization of more complex microstructures (Moreaud , 2006), but these models are not considered here.Results on Boolean models and on multi-scale Cox Boolean models are presented in the next part.The stopping condition is "when each point of S is disconnected to the others".This is the least restrictive condition.In practice, we stop the display of curves when uncertainty exceeds 2. The confidence level is at 95 %; therefore the uncertainty is defined as twice the standard deviation divided by the square root of N.

BOOLEAN MODELS SPHERES AND SPHEROCYLINDERS
The degree of impact of grains' morphology has been studied by comparing Boolean models of spheres and spherocylinders with random orientations (cf., Fig. 15).A spherocylinder is a cylinder with two hemispherical caps at each end, thus defined by two parameters; L the length of the cylinder and R the radius of the hemispheres.The first Boolean model of spheres is defined by parameters R = 10 and V v = 0.7.The Boolean model of spherocylinders is defined by R = 5, L = 47 and V v = 0.7.Grains parameters are chosen to have similar grain average volume V (A ).n = 20 realizations for each scheme are generated and M-tortuosity-by-iterative-erosions is applied on them.The focus is on the M-scalar, τM,r , as a function of r, radius of the percolating sphere.Results highlight the shape effects (Fig. 16).Indeed, τM,r increases faster in the case of spherocylinders than in the case of spheres.The percolation threshold of Boolean models of 13 spherocylinders is lower than for Boolean models of spheres (Jeulin and Moreaud , 2006), the porous volume is disconnected faster, which explains that τM,r increases faster, with respect to r.

MULTI-SCALE COX BOOLEAN MODELS OF SPHERES
Let us now consider multi-scale Cox Boolean models of spheres with constant radius.Fig. 17 displays a realization of each model.Two scales are used; aggregate scale, defining non-empty areas, and grain scale inside aggregates, inclusion spheres.Space outside inclusion spheres is kept empty.Radius of spheres, i.e., grains, is given, R = 3, and is part of constant parameters, as well as grain volume fraction V v = 0.7 and V v,out = 0.Only aggregate sphere radii R inc and volume fractions V v,inc vary; V v,inc = 0.4 and R inc = 10 for the first model, and V v,inc = 0.6 and R inc = 20 for the second (Fig. 17  M-tortuosity-by-iterative-erosions results are presented in Fig. 18.As expected the model with a higher volume fraction V v,inc is more tortuous than the other one.Indeed, the larger the radius r of the spherical percolating particule, the more the model's low scales prevails.Hence, the gap between two models increases with r.

MULTI-SCALE COX BOOLEAN MODELS OF PLATELETS
The final application focuses on a specific modelization using multi-scale Cox Boolean models of platelets, simulating alumina catalysts (Wang et al., 2015) (cf., Fig. 19).Platelets (Chiche et al., 2008) are defined by the parameters H, L, l and l (Fig. 20).Two models are compared with constant parameters; platelets' parameters (H = 3, L = 105, l = 5, l = 7), V v = 0.4, V v,out = 0.2, and different parameters; V v,inc = 0. M-tortuosity-by-iterative-erosions, especially τM,r , allows to separate very similar multi-scale Cox Boolean models (cf.,Fig. 18 and 21).This is partially done by consideration of narrowness, thanks to the use of morphological erosion.Moreover, shape effect is also highlighted by our descriptor (see Fig. 16).

CONCLUSION AND PERSPECTIVES
We have proposed a new descriptor, TM , named M-tortuosity, based on geometric tortuosity.It characterizes porous media sinuosity.TM is a scalable descriptor which can have a proper dimensionality and is applicable on very complex and disconnected porous microstructures.Unlike definitions found in the litterature, our method, using stochastic points sampling (Monte Carlo method), does not use any arbitrary entry and exit planes or points.Translational, rotational and homothety invariances have been proven together with periodic pattern repetition stability.Extension of M-tortuosity, named M-tortuosity-byiterative-erosions, improves the discriminative power by taking into account narrowness, especially the bottleneck effect, connecting it conceptually to constrictivity.Application to basic toy cases validates M-tortuosity behavior.Boolean models of spheres and spherocylinders highlight consideration of shape effect characteristic by M-tortuosity-by-iterative-erosions.Finally, discriminations of multi-scale Cox Boolean models of spheres first, then platelets, account for the good discriminative power of our descriptor.M-tortuosity is characterized, among other things, by its scalability.This notion means here that information of different dimensions can be extracted from our descriptor.In this paper, the focus is on the final M-scalar, τM , for discrimation purposes.Such a description is not exhaustive.Information of higher dimensions, for instance M-coefficients distribution or the 3D map of mean geometric tortuosity, can improve the accuracy of the characterization.The 3D map of mean geometric tortuosity is defined thanks to the N geometric tortuosity 3D maps (one per starting point), on which an arithmetic mean is computed for each point of Sk over the N maps.This complete description using all available information of the M-tortuosity descriptor will be discussed in a further paper.Moreover, characterization of a large set of zeolites (catalysts with crystalline structure) will be performed and 3D nano volumes of alumina catalyst supports, obtained by nanotomography (Tran et al., 2014), too.Correlations with experimental physicochemical characteristics will be looked for.Finally based on this method, two novel topological descriptors will be proposed by first extending M-tortuosity to graylevel images, in order to avoid a segmentation step, which can be very delicate in some applications, and second, adapting the M-tortuosity formalism to quantify heterogeneity at small scale.

Fig. 1 :
Fig. 1: Two pores (black lines); a straight one (dashed lines) and a sinuous one (curved lines).The skeleton of both pores is represented: the dashed straight blue line for the straight pore, and the curved red line for the sinuous one.

Fig. 2 :
Fig. 2: Cases without skeletonization (a) and with skeletonization (b).Scatterplot of coordinates (d, g): d (resp.g) is the Euclidean distance (resp.geodesic distance) of a point pair.520 locations are sampled, each pair of locations is processed.The linear regression of the scatterplot is displayed (orange line) with its equation, and the reference line g = d corresponds to geometric tortuosity equal to one.Three point pairs are represented by a black square (resp.circle and triangle).

Fig. 3 :
Fig. 3: Evolution of N times the variance of τM (Eq.10), σ 2 , over the 100 realizations of Boolean models of size 200 3 and of volume fraction V v = 0.7 of spheres of constant radius r = 3, as a function of N.

Fig. 8 :
Fig. 8: Examples of skeleton Sk, represented by the broken lines, with the (a) 1D-sampling and (b) stratified sampling methods.S is the set of red or blue dots.

Fig. 13 :
Fig. 13: τM values according to the index of (φ 2 lim,down , φ 2 lim,up ) for the P 2 case.Confidence intervals, with confidence level at 95 % (see definition below), are represented by vertical bars.

Fig. 14 :
Fig. 14: Top figures: Boolean model of platelets and Boolean model of spheres denoting aggregate scale.Bottom figures: intersection of the two Boolean models, by standard operation or using Cox point process yielding a more realistic microstructure(Moreaud et al., 2018).
3 and R inc = 15 for the first model, and V v,inc = 0.4 and R inc = 25 for the second.The M-scalar τM,r increases faster in the higher volume fraction model again, as shown in Fig. 21.The final radius, r = 2, is lower than in the previous application because of the non-empty space outside inclusion spheres.
Fig.20: The shape of a platelet grain and its parameters(Chiche et al., 2008).
the translated set of I by vector t.Let S t and Sk t be the translated sets of S and Sk, respectively.By definition, the geodesic distance transform, D G (., S; Sk), is invariant with respect to translation.This result is obviously true for the Euclidean distance transform too.The geometric tortuosity is then translation invariant.Whatever x ∈ Sk and t ∈ R 3 ,