STATISTICAL ANALYSIS OF THE LOCAL STRUT THICKNESS OF OPEN CELL FOAMS

Open cell foams are formed by an interconnected network of struts whose thickness varies locally. These variations are known to have an impact on the elastic and thermal properties of the foam. In this paper we quantify the local strut thickness by means of micro computed tomography (μCT) imaging. We develop a fully automatic algorithm to extract the foam’s skeleton from a binary image and its topological decomposition into vertices and struts. This allows to estimate the thickness of individual strut segments by the Euclidean distance transform, where an appropriate correction for struts with nonspherical cross-sectional shape is applied. Conflating these estimates based on the strut lengths results in a strut thickness profile for the entire foam. Based on this profile we give a statistical justification that a strut thickness model should depend on the strut length. Furthermore, the investigation of polynomial models for the strut thickness profile by means of regression analysis leads to a new three-parameter strut thickness model.


INTRODUCTION
Foams are nowadays used in a wide range of application areas including heat exchangers, filters, insulators or sound absorbers (Banhart, 2001).In this work we are interested in so-called open cell foams which are formed by a continuous network of struts.The macroscopic properties of a foam such as thermal conductivity, permeability, elasticity or acoustic absorption are highly influenced by the microstructure.Therefore, an understanding of the reaction of these properties to changes of the microstructure is crucial for the optimization of foams for given applications.
Three-dimensional images obtained by micro computed tomography (µCT) are a valuable source of information on the microstructure geometry of foams.Geometric characteristics which can be estimated from image data include the volume fraction, the specific surface area, distributions of cell size or shape as well as mean value characteristics of the strut system, e.g., the mean strut thickness (Ohser and Schladitz, 2009).
Using these characteristics models based on random tessellations can be fit to the observed foam structure.By change of the model parameters, virtual foams with altered microstructure can be generated.Application of finite element techniques then allows to study the influence of certain geometric microstructure characteristics on macroscopic properties of the foam (Roberts and Garboczi, 2002;Li et al., 2006;Kanaun and Tkachenko, 2006;Redenbach, 2009;Tekoglu et al., 2011;Liebscher et al., 2012).
A typical feature of open foams is that the strut thickness varies locally.Usually, the struts are thicker at the vertices than at the centers (Fig. 1).This local thickness variation was shown to have an impact on the elastic and thermal properties of the foam (Jang et al., 2008;Kanaun and Tkachenko, 2008).Therefore, it should be included in the model.In Lautensack et al. (2008), a foam model based on locally adaptable dilation of the edges of a random tessellation was introduced.The local strut thickness was modeled by a quadratic function depending on the distance from the vertices.The function was parametrized using mean strut characteristics which were estimated from µCT image data.Owing to the parametrization with mean values, the model did not reproduce the variability of the real foam.
In the same year, Jang et al. (2008) presented an analysis of the strut length distribution and the cross-section area in polymeric and aluminum foams.Cross-section areas were measured from sections of the foam's struts leading to a polynomial strut thickness model.They further revealed a nonlinear dependence of the mid-strut thickness on the strut length.This was incorporated into their microstructure model by forming two length classes and scaling the normalized strut thickness by a function depending on the strut length.Hence, modelling local strut thickness disregarding the strut length is insufficient.
In this paper, we extend the approach of Jang et al. (2008) by a systematic analysis of the complete strut system of the foam.For our study we develop a fully automatic extraction of the foam skeleton adopting the skeletonization technique of Chaussard et al. (2010).Further, we introduce a novel algorithm for decomposing the skeleton into individual curve segments.
Based on the decomposition we establish the notion of the local spherical contact profile as the radius of the maximal inscribed ball at any point of a curve segment.A regression analysis of the extracted profile reveals that the strut thickness depends on the strut length in a more complex way than could be expressed by scaling.As result we introduce a new three-parameter model for the strut thickness that is superior compared to existing models (van Merkerk, 2002;Jang et al., 2008;Kanaun and Tkachenko, 2008).
The paper is organized as follows: We start with a survey on digital images and digital topology.These notions are essential for the following two sections.In these we present the skeletonization technique that forms the core of our algorithm and explain how the skeleton can be decomposed into its topological components.Based on this decomposition we introduce the local spherical contact profile to describe the varying strut thickness.Finally, the algorithm is applied to a tomographic image of an open cell aluminum foam, for which we conduct a regression analysis of its spherical contact profile.We conclude with a discussion of our results and an outlook to future work.

BASIC NOTATIONS
In this section we summarize the basic definitions this work is founded on.We adhere to the notations given in Couprie et al. (2007).For a more detailed survey we refer to Rosenfeld (1981), Nakamura and Aizawa (1985) and Kong and Rosenfeld (1989).

NOTIONS FOR DIGITAL IMAGES
We denote by Z the set of integers, by N the set of strictly positive integers and by R the set of real numbers.The set of positive real numbers is denoted by R + .A three-dimensional binary image is defined as a subset V I of the cubic grid Z 3 together with a function f : V I → {0, 1}.We refer to X ⊂ V I mapped to {1} as foreground and to its complement X, that is mapped to {0}, as background.A point x = (x 1 , x 2 , x 3 ) ∈ V I with x i ∈ Z is called a voxel if we want to emphasize its volumetric extent.The Euclidean distance of two points x and y in V I is given by d The Euclidean distance of X ⊂ V I and x ∈ V I is defined as d(x, X) = min y∈X d(x, y).The mapping D X (x) = d(x, X) that assigns to each point x ∈ X its distance to the complement X of X is called Euclidean distance map.Note that for all x ∈ X, the value D X (x) is the radius of the maximum ball centered in x that is completely contained in X. D X can be computed in linear time using any of the algorithms presented in Meijster et al. (2002), Maurer et al. (2003) or Schouten et al. (2006).

SURVEY ON DIGITAL TOPOLOGY
The neighborhood of a point x ∈ V I is defined as the set of all points in V I that are adjacent to x with an adjacency notion depending on some distance.In the remainder of this work, we consider the neighborhoods and that are induced by the city-block and the chessboard distance.The neighborhood obtained from N 26 (x) by removing its outer corners is given by We define N * n (x) = N n (x) \ {x}.Note that the number of points in N n (x) is n + 1.
Let x and y be two points in V I .We say that y is nadjacent to x iff y ∈ N * n (x).A set Y ⊂ V I is n-adjacent to x ∈ V I iff there exists a point y ∈ Y which is n-adjacent to x.We call a sequence of points x 0 , . . ., x k in V I an n-path, if x i−1 is n-adjacent to x i for all 1 ≤ i ≤ k.Given a nonempty subset X ⊂ V I , two points x, y ∈ X are called n-connected iff x and y are joined by an npath in X.This definition fulfills the three axioms of an equivalence relation (reflexive, symmetric, transitive) and hence its equivalence classes induce a partition of X into its n-connected components, in the following abbreviated as n-components.We call X ⊂ V I nconnected if it consists of exactly one n-component.
The set of all n-components of X that are n-adjacent to a point x is denoted by C n [x, X].
To have a correspondence between the topology of the foreground and the background of V I -and thus to avoid connectivity paradoxes (Kong and Rosenfeld, 1989) -we must assign them with complementary neighborhoods.For instance, if we use N 26 for the foreground of V I then we have to use N 6 for the background and vice versa.In the sequel, we assign N 26 to the foreground and N 6 to the background.

EUCLIDEAN SKELETON
The first step of our procedure to estimate the local strut thickness is to compute the skeleton of the foam structure.Cenens et al. (1994) were the first who introduced a suitable algorithm for three-dimensional foams based on the watershed transform (Soille, 2002).However, this algorithm must be carefully parametrized to avoid over-or undersegmentation, respectively.Even with an optimal parametrization we could not avoid undersegmentation at the boundaries completely.Hence the watershed transform seems not suitable for our fully automated approach.
A parameter free concept for skeletonization is the sequential removal of points from the original shape that are not necessary to preserve the topology of the structure.As it is important for our application that the skeleton is centered within the foam structure, the main challenge is to determine the order in which these points shall be removed.But before we devote our attention to this problem, we give a formal definition of skeletonization as used throughout this work.The following sections are based on the work of Couprie et al. (2007).

SIMPLE POINTS
The notion of simple points is essential for the definition of topology preserving skeletons.Informally speaking, a point of X ⊂ V I is called simple if its removal from X does not change the number of connected components and tunnels of X and X.
In the literature various methods for characterizing simple points have been proposed (Kong and Rosenfeld, 1989;Couprie and Bertrand, 2009;Evako, 2011).We adopt the approach based on counting the 26-and 6-components in the neighborhood of a point (Bertrand, 1994;Bertrand and Malandain, 1994).More precisely, the two topological numbers for a point x ∈ X ⊂ V I are defined as where #X denotes the cardinality of X.The point x is then simple (for X) iff T 26 (x, X) = 1 and T 6 (x, X) = 1.Note that we use N * 18 in Eq. 6 to avoid explicit checking for holes.For a detailed explanation we refer to Bertrand (1994).

ULTIMATE HOMOTOPIC SKELETON
Using this notion we can now define the skeleton of a finite subset X ⊂ V I .We call Y ⊂ V I a homotopic thinning of X if Y = X or Y was obtained from X by iterative deletion of simple points.An ultimate homotopic skeleton of X is a homotopic thinning Y of X with no simple points left in Y .
The definition of the homotopic thinning itself already provides an algorithm to compute the ultimate homotopic skeleton.However, the result of the skeletonization will depend on the order in which points are deleted.This order is guided by a priority function.Points with a low priority shall be removed before those with a higher priority.The ultimate homotopic skeleton can be computed with time complexity O(n log(n)) independently of the priority function (Couprie et al., 2007).

ROBUST EUCLIDEAN SKELETON
The remaining question is how to choose the priority function to ensure that the skeleton is centered with respect to the Euclidean distance.Using the Euclidean distance map may lead to unnatural branching of the skeleton (Talbot and Vincent, 1992) or even spurious branches (Chaussard et al., 2010).
On the other hand, being centered implies the inclusion of the ridge points of the Euclidean distance map in the skeleton.Formally, for X ⊂ V I , ridge points are either (a) the centers of all balls in X that are not strictly included in any other ball contained in X (Blum, 1967) or (b) all points x ∈ X that have at least two closest points in X.Both definitions differ only by a negligible set of points (Matheron, 1988).The set comprised of all ridge points is known as the medial axis.
By construction, the medial axis is centered in the object with respect to the distance used in its definition.If we weight each point of an object X with the radius corresponding to its medial axis point we obtain a function that guides the thinning process towards the innermost medial axis points (Chaussard et al., 2010).However, the medial axis is rather sensitive to noise on the object boundary and thus in applications often some sort of filtering has to be applied (Attali et al., 2009).
To alleviate the stability problem, Chazal and Lieutier (2005) extended notion (b) to the λ -medial axis, which was recently transferred to the digital framework by Chaussard et al. (2010).It was shown to be more stable under small perturbations of the boundary and also to be less sensitive to rotations than the classical medial axis.These properties are inherited by the skeleton.
Let S be a subset of R n .We denote by the radius of the smallest ball enclosing S. Furthermore, let X be a subset of V I .The extended projection of x on X is given by where denotes the plain projection of x ∈ X on the complement X of X. Couprie et al. (2007) then define the projection radius map as the function F X (x) = R(Π e X (x)) that assigns to each point x ∈ X the radius of the minimum ball enclosing all extended projections of x on X.Using the extended projection in the discrete framework avoids the reduction of the projection to a singleton in ill-conditioned cases (Chaussard et al., 2010, Sec. 3).
Thresholding F X with λ > 0 yields the λ -medial axis at level λ .By definition all λ -medial axes at level λ ′ > λ are included in the axis at level λ .The computational cost for computing the extended projection is O(n) (Couprie et al., 2007).A stochastic algorithm for computing enclosing spheres in expected linear time was proposed by Welzl (1991).A slightly tuned version was given by Gärtner (1999).
In our application we do not use the λ -medial axis itself.We rather apply the radius projection map F X as priority function during the thinning process that automatically guides the process to the topological kernel of the several nested λ -medial axes.That is, the smallest subset of all λ -medial axes having the same topology as the original shape.As a consequence the skeleton is curvilinear (Fig. 7a).
During homotopic thinning topological artifacts may evolve if a semi-continuous function such as F X is used (Passat et al., 2007;2008).To avoid this behavior we enforce the thinning process to start with the outermost points of the shape by preordering all points with respect to D X .Note that a partial order is sufficient as points with the same priority may be removed in arbitrary order.

TOPOLOGICAL DECOMPOSITION OF SKELETONS
The second step of our procedure consists of the decomposition of the skeleton into its topological components.As we are dealing with open cell foams we work under the model assumption that the skeleton consists of curves and curve junctions, only.Hence, each skeleton point has to be labeled as either a curve or a curve junction point.
In the context of characterizing simple points two methods were proposed for such a classification.Note that both methods are not restricted to curves and curve junctions.In the following two sections we give a short survey on these methods and their limitations.Subsequently, we will propose a new algorithm that overcomes these limitations.

CLASSIFICATION BY TOPOLOGICAL NUMBERS
Denote by x a point of a skeleton X.According to Malandain et al. (1993), x can be classified topologically by evaluating T 26 (x, X) and T 6 (x, X).As the skeleton contains no surfaces, T 6 (x, X) is always one and therefore only T 26 (x, X) is meaningful.So x is part of a curve in X if T 26 (x, X) = 2 and part of a curve junction if T 26 (x, X) > 2, respectively.Note that Malandain et al. (1993) used N 18 instead of N * As a consequence of discretization effects "thick junctions" may be misclassified as depicted in Fig. 2 (top row), where all points are classified as curve points.Such misclassifications are corrected by postprocessing the preliminary classification: For each curve point the neighbors that are curve or curve junction points are counted.If the number of such neighbors exceeds two, the point is classified as curve junction.
However, the algorithm also allows for the classification of points of other topological types, e.g., surface points or surface junctions.In contrast to our model assumption there may be point configurations of the skeleton that are recognized accordingly.An example is shown in Fig. 2 (bottom row), where the central point is falsely classified as a single surface point.Such misclassifications are ignored by the correction scheme.

CLASSIFICATION BY BETTI NUMBERS
Another classification scheme was introduced by Saha and Chaudhuri (1996).It is based on the binary versions of the Betti numbers in three-dimensional space: the number of connected components B 1 , the number of tunnels B 2 , and the number of cavities B 3 .Let x be a point of a skeleton X.Then the first Betti number of N * 26 (x) ∩ X is given by B 1 (x, X) = T 26 (x, X).Following Saha and Chaudhuri (1996) the number of tunnels is zero if all 6-neighbors of x are also contained in X. Otherwise it is one less than the number of 6components in N * 18 (x) ∩ X that contain at least one 6neighbor of x.More precisely, the number of tunnels of x ∈ X is defined as otherwise.Note that only 6-connected tunnels are considered by this definition.As the skeleton contains no surfaces the number of cavities is always zero.If B 2 (x, X) is zero, we get the same result as by topological number classification without correction.Otherwise there are several possible topological classes for a point x and postprocessing is necessary.For each such point, the topological type of its neighbors in the skeleton is checked: If all neighbors of x are curve points or curve junctions, then x belongs to a curve junction (Fig. 2, bottom row) for an example.Though, for the "thick junction" as depicted in Fig. 2 (top row) the algorithm falsely classifies the junction points as curve points.

CLASSIFICATION BY NEIGHBOR COUNTING
For our application we have two requirements on a point classification scheme.Firstly, it is crucial to reliably detect curve junctions.Secondly, if the curve junction points are removed from the skeleton, its curves should become segmented.The approaches introduced above fail on both requirements even on rather simple point configurations.Inspired by the postprocessing method used in the topological numbers approach, we propose a new classification method that fulfills both requirements.Note that a similar approach has also been used by Montminy (2001).
The basic idea of our approach is to count the number of neighbors of each skeleton point.If this number is greater than two, the point belongs to a curve junction.Otherwise it is a curve point.If possible, the size of the junction is reduced by relabeling junction points which are not necessary to ensure the segmentation of the skeleton's curves as curve points.Fig. 2 illustrates one out of four ("thick junction") or six ("thin junction") possible configurations found by our algorithm.
To determine the order in which points are removed a priority function is necessary.We suggest the Euclidean distance map as priority function as it guarantees that points close to the boundary of the original shape are removed first and the most centered points remain as curve junction.In addition we must assure that no topological changes occur during the thinning step.Therefore, we allow only points to be removed from a junction that are simple for it.
To make this more rigorous let us introduce some definitions.We denote the neighborhood of Let X be a skeleton and C ⊂ X the set of all 26-components of X that form a curve junction.We now iterate over all curve junctions J ⊂ X and remove all simple points j from J (in decreasing order w.r. t.D X ) for which the number of 26-components in N 26 (J ′ )∩X \J ′ does not change, i.e., (12) where J ′ = J \ { j}.Eq. 12 guarantees that the curves of the skeleton become separated from each other if all curve junction points are removed from the image.

LOCAL SPHERICAL CONTACT PROFILE
In the third step of our procedure we estimate the length and the local thickness of the individual strut segments based on the decomposition of the skeleton introduced in the previous section.These estimates are then combined to a thickness profile of the entire foam.
In Jang et al. (2008) the strut thickness is defined as the cross-sectional area.To compute it we would have to slice the strut, which is computationally expensive and owing to discretization errors not robust.Therefore, we employ the radius of the struts' inscribed ball as measure for the strut thickness.This corresponds to the incircle of a cross-section and can be computed in linear time using the Euclidean distance transform.
Note that for modelling there is no practical difference between both approaches.Knowing the cross-sectional shape in conjunction with the radius of the incircle provides the same information as combining the cross-sectional shape and its area.As the cross-sectional shape is generally known both definitions can be converted into each other.

COMPUTING THE SPHERICAL CONTACT PROFILE
Let us denote by Y the set of curve segments of the decomposed skeleton of an open cell foam X.Each strut of the foam contributes a curve segment Z ∈ Y which is assumed to be of length ℓ.The local strut thickness at point x ∈ Z is defined as the radius of the largest ball with center x inscribed in X.It may be parametrized using the Euclidean distance of x to the strut center x 0 x > x 0 (13) yielding a function p Z (ξ ), with ξ ∈ − ℓ 2 , ℓ 2 (see Fig. 3 for an illustration).The expected spherical contact profile of the entire foam is defined as p(ξ | ℓ) = E(p Z (ξ ) | ℓ), that is the mean thickness at distance ξ of all struts Z with length ℓ.
The struts of real foam samples are often slightly curved.We therefore approximate the length of a curve segment Z by the sum of the Euclidean distances from the center x 0 of the struts' curve in the skeleton to the adjacent curve junctions x and y, that is ℓ

VOLUME CORRECTION
According to Gibson and Ashby (1999), the volume density (also known as relative density in physics) is the single most important parameter in mechanical and thermal properties of foams.Hence it should be captured accurately in microstructure models.As the spherical contact profile measures the radius of the inscribed ball of a strut, the volume of struts with nonspherical cross-sectional shape will be underestimated.Therefore, when generating microstructure models, we scale all measurements p Z for a given curve segment Z by an individual volume correction factor c Z to yield the same volume as the corresponding strut.
To compute the volume correction for an individual strut we must first separate it from its adjacent nodes.In a real foam the nodes are formed by complex minimal surfaces which join smoothly with their connecting strut segments.Two nodes from an aluminum foam are shown in Fig. 4. As the centers of the nodes coincide with the junctions of the skeleton we define a node as the region occupied by the ball B r (x) with radius r = D X (x) centered in the curve junction x.

Fig. 4: Two nodes from a 26-ppi aluminum foam
The strut segment between any two junctions x, y ∈ Y is computed as follows: Let L xy be the vector from x to y whose normalized direction is denoted by r.Further, define two planes p 1 and p 2 that pass orthogonally to L xy through x + D X (x)r and y − D X (y)r, respectively.The connecting strut segment is composed of all points of the strut between p 1 and p 2 .See Fig. 5 for an illustration along with a visualization of a foam cell whose nodes were removed.
We define the volume correction c Z for a strut as the ratio of the cross sectional area to the incircle averaged over all slices of the corresponding strut segment.Let C be the strut segment whose curve is given by Z ∈ Y .Then c Z is computed as where v(C) denotes the volume of C as estimated by voxel counting from the image.The denominator corresponds to the sum of the circle areas (computed in voxels) between the planes p 1 and p 2 .We assume here, that the volume contribution of the nodes adjacent to C are of the same magnitude as for C.

APPLICATION TO OPEN METAL FOAM
To evaluate our algorithm we applied the technique to an 26-ppi open cell aluminum foam.The analysis is based on a tomographic image of size 630 x 1370 x 1370 voxels with an isotropic voxel edge length of 14.16 µm.This corresponds to 8.9 mm × 19.4 mm × 19.4 mm of material.Owing to the image's high contrast we could binarize it using Otsu's (1979) threshold without any preprocessing.Furthermore, some small holes in the strut system were morphologically closed to prepare the binarized image for the analysis.

SKELETONIZATION OF THE FOAM
The skeleton was computed as stated previously.A subvolume of the foam and its corresponding skeleton can be seen in Fig. 6.We assessed the topological correctness of the skeleton by comparing the Euler number of the skeleton and the binarized image of the foam using the MAVI software package (Fraunhofer ITWM, 2012).Both yielded the value −8898 (computed in the 26-neighborhood, Ohser and Schladitz, 2009).The drawback of topology preservation is that structural peculiarities of the foam are also captured by the skeleton.Two common ones are shown in Fig. 7.In the case of closed faces, one curve attached to the appropriate face is missing in the skeleton.Its influence on the estimated spherical contact profile is negligible.On the other hand, the short curves that evolve within eight-fold vertices may bias the estimate of p(ξ | ℓ) and thus must be ignored.Hence, we only considered curves that contain at least 10 voxels in our study.
The resulting skeleton was decomposed into 15 494 curves.For approximately 84% of those the volume correction factor could be determined.The remaining struts were either part of a closed or a partially closed face, for which we used the mean volume correction factor c = 1.58.

ANALYSIS OF THE SPHERICAL CONTACT PROFILE
After minus-sampling edge correction 8170 struts remained.From those we removed nonsymmetric outliers that were caused by broken or heavily bent struts, leaving 7058 struts for the analysis.Fig. 8 depicts the volume corrected estimate of the spherical contact profile p(ξ | ℓ).It shows the typical flat form in the strut center (ξ = 0) which steeply rises towards the nodes (increasing absolute value of ξ ).
An extremal behavior is shown below the fifth percentile of the strut length distribution where there is almost no increase of strut thickness while approaching the node.This finding is supported by the boxplot of the estimated mid-strut thickness p(0 | ℓ) depicted in Fig. 9.We subdivided the struts into the percentile ranges q i, j = (x i , x j ], where x i denotes the ith percentile of the strut length distribution.The median and the variation of the mid-strut thickness decrease with increasing strut length, while the box of q 0,5 stands almost on top of the other ones.This indicates a nonlinear dependence of p(0 | ℓ) on ℓ as reported by Jang et al. (2008).q 0,5 q 5,10 q 10,20 q 20,30 q 30,40 q 40,50 q 50,60 q 60,70 q 70,80 q 80,90 q 90,100 To study the spherical contact profile in more detail we conducted a regression analysis using weighted least squares.Following Jang et al. (2008) we first normalized the data with respect to mid-strut thickness p(0 | ℓ) and strut length ℓ.Let us denote by ξ = ξ / ℓ the distance from the strut center ξ normalized by ℓ.We then assume a polynomial of degree of at most eight as initial model: a 1 , . . ., a 4 ∈ R, where the indices denote the exponents of ξ that are included.The odd exponents were neglected due to the symmetric nature of the strut profile.Caused by the normalization the polynomial model shall equal one at ξ = 0 and thus an intercept of one was used.
Fig. 10 shows a plot of the measurements for the shortest (q 0,10 ), medial (q 45,55 ) and longest (q 90,100 ) 10% of the struts normalized by p(0 | ℓ) and ℓ along with the initial model p2468 fitted to the individual profiles.The resulting curves differ heavily from each other, which indicates a dependence of p(ξ | ℓ) on ℓ.Therefore, we subdivided the data into 10% length percentile ranges q i, j for our analysis.
Table 1 summarizes the results of the all-subset regression of p2468 guided by the MC p -statistic (Fujikoshi and Satoh, 1997).MC p is a modification of Mallows' (1973) C p -statistic that was shown to be the minimum variance unbiased estimator of the expected overall Gauss discrepancy (Davies et al., 2006).Roughly speaking: The smaller the MC p -value the better is the model as the disparity compared to the initial model is small while, simultaneously, less parameters are needed.In the first three columns the pvalues of a t-test for the significance of the individual parameters in p2468 are given (null hypothesis H 0 : a i = 0 against the alternative H 1 : a i = 0).The p-values for a 1 were left out as they are all zero and thus ξ 2 is mandatory.No general model for strut lengths smaller than x 50 could be found, as in this region the central plateaus start to grow.We also observed, that only struts smaller than x 30 could be reliably modeled by a two-parameter model.However, for all but q 0,10 the difference between the best two-parameter and the best three-parameter model is close to two.Thus it is only caused by the parameter penalization of the MC pstatistic.For the three-parameter models, the smallest MC p -values were mainly obtained by p248 .The only exception is q 40,50 , but with a negligible difference.
For strut lengths greater than x 50 , with the exception of q 90,100 , p248 yielded the smallest deviations from the initial model.This behavior is backed by the p-values and the observation that long struts tend to have the same shape: a wide flat plateau in the center which steeply rises towards the nodes.Only for the 10% longest struts the plateau was captured better by including ξ 6 instead of ξ 4 .This behavior is expected as the coefficients of higher order variables are negative and thus serve as penalization to fit wide plateaus better (cf.Table 2).
Note that in the literature (Kanaun and Tkachenko (2008, Eq. 51) -Eq. 16;van Merkerk (2002, Eq. 4.1) -Eq.17; Jang et al. (2008, Eq. 1) -Eq.18) the following subsets of p2468 have already been proposed to describe the varying strut thickness: p6 ( ξ ) = a 3 ξ 6 + 1, ( 17) All these models assume that the strut thickness scales with length.Here, however, p2 and p6 are disqualified: The first one by its high MC p -values and the second one by the mandatory but missing ξ 2 .Among these models only p24 showed a reasonable performance.
For strut lengths greater than x 50 its MC p -values were at least two times smaller than the ones reached by the other two-parameter models (with exception of q 50,60 and q 60,70 for p26 were it is about 1.6 and 1.8, respectively).Thus, it can be considered as the best two-parameter model.
In summary, the three-parameter model p248 yielded the best representation of the data.Its coefficients are given in Table 2.For comparison we also fitted a nonparametric smoothed spline model (Green and Silverman, 1995) to each q i, j and compared its R 2 -statistics with the ones of the initial parametric model p2468 .In all cases the values lay within 0.94 and 0.99 with no observable differences among the two models.As no model selection procedure is necessary for the spline model it is a feasible alternative to the polynomial model.

CONCLUSIONS
In this work we presented a complete method to quantify the local strut thickness of open cell foams.Our approach is based on the skeletonization of a binarized µCT image of the foam structure that is decomposed into its struts and vertices.The spherical contact profile is obtained by evaluating the Euclidean distance transform on every point of the skeleton that is identified as strut.
To accurately capture the volume fraction for struts with a nonspherical cross-sectional shape we incorporated a correction factor for the thickness measurements.The computation of this factor is based on a separation of the strut segments from the nodes in the foam.This could be avoided by employing a distance transformation whose structuring element is adapted to the cross-sectional shape of the struts.
Using the fitted polynomials a foam can be modeled with the observed strut thickness (Liebscher and Redenbach, 2012).Also the microstructure models presented in Jang and Kyriakides (2009) or Liebscher et al. (2012) could benefit from the new thickness model.As they use circular cross-section shape the new model could be readily applied.The development of a suitable two-dimensional model in ξ and ℓ avoiding the consideration of different length classes is topic of our ongoing research.

Fig. 1 :
Fig. 1: Visualization of a strut of an open cell aluminum foam showing the locally varying strut thickness.

Fig. 2 :
Fig. 2: Two types of curve junctions that are hard to recognize in three dimensions (illustrated in two dimensions -the upper and lower slice are empty).The columns show the classification results of different algorithms.

Fig. 3 :
Fig. 3: Illustration of the spherical contact profile for an individual strut.

Fig. 5 :
Fig. 5: (a) Illustration of the model used to separate strut segments from the adjacent node.(b) Foam cell whose nodes (dark parts) were removed using the model shown in the illustration (a).

Fig. 6 :
Fig. 6: Cell of the aluminum foam and its skeleton.

Fig. 7 :Fig. 8 :
Fig. 7: Structural pecularities found in the aluminum foam and their effect on the skeleton: Closed face (a) and eight-fold vertex (b)

Fig. 9 :
Fig.9: Boxplot of the estimated mid-strut thickness p(0 | ℓ) subdivided into its length percentile ranges q i, j .The median of p(0 | ℓ) is depicted as gray line in the background.

Fig. 10 :
Fig. 10: Plot of the measurements for the shortest (foreground points, black), medial (intermediate points, red) and longest (background points, blue) 10% of the struts normalized by p(0 | ℓ) and ℓ.Model p2468 was fitted individually to shortest, medial and longest struts using weighted least squares.

Table 1 :
Results of the all-subset regression analysis of model p2468 .The best model for each q i, j is highlighted.

Table 2 :
Parameters of the polynomial model p248 individually fitted to the length percentile ranges q i, j .