Image Anal Stereol 2013;32:1-12 Original Research Paper
STATISTICAL ANALYSIS OF THE LOCAL STRUT THICKNESS OF
*
OPEN CELL FOAMS
Andre Liebscherb and Claudia Redenbach
University of Kaiserslautern, Department of Mathematics, Erwin-Schrödinger Straße, 67663 Kaiserslautern, Germany
e-mail: liebscher@mathematik.uni-kl.de, redenbach@mathematik.uni-kl.de (Received June 13, 2012; revised November 16, 2012; accepted November 16, 2012)
ABSTRACT
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 (uCT) 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.
Keywords: digital topology, Euclidean skeleton, image analysis, regression.
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 etal., 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.
1 2 3 4 5 6 7
1	2 3 4 5 6 7
Fig. 1: Visualization of a strut of an open cell aluminum foam showing the locally varying strut thickness.
In Lautensack etal. (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 topic of this paper was presented at the S4G Conference, June 25-28, 2012 in Prague, Czech Republic.
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 VI of the cubic grid Z3 together with a function f : Vi ^ {0,1}. We refer to X c Vi_mapped to {1} as foreground and to its complement X, that is mapped to {0}, as background. A point x = (xi,x2,x3) e VI with xi e Z is called a voxel if we want to emphasize its volumetric extent. The Euclidean distance of two points x and y in VI is given by d(x, y) = a/£i(xi — yi)2. We denote the ball of radius r > 0 centered in x e VI by
Br(x) = {y e Vi | d(x, y) < r} .	(1)
The Euclidean distance of X c VI and x e VI is defined as d(x, X) = minyeXd(x, y). The mapping DX (x) = d(x, X) that assigns to each point x e X its distance to the complement X of X is called Euclidean distance map. Note that for all x e X, the value DX (x) is the radius of the maximum ball centered in x that is completely contained in X. DX 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 e Vi is defined as the set of all points in Vi that are adjacent to x with an adjacency notion depending on some distance. In the remainder of this work, we consider the neighborhoods
N6(x) = {y e Vi | £ |x,- — yi| < 1} (2) i
and
N26(x) = {y e Vi | max |xi — yi| < 1} , (3) i
that are induced by the city-block and the chessboard distance. The neighborhood obtained from N26(x) by removing its outer corners is given by
N18(x) = {y e Vi | £ |x,- — yi| < 2} nN26(x) . (4) i
We define N*(x) = Nn (x) \ {x}. Note that the number of points in Nn(x) is n + 1.
Let x and y be two points in VI. We say that y is nadjacent to x iff y € N* (x). A set Y C VI is n-adjacent to x € VI iff there exists a point y € Y which is n-adjacent to x. We call a sequence of points xo,...,xk in VI an n-path, if xi_1 is n-adjacent to xi for all 1 < i < k. Given a nonempty subset X C VI, two points x,y € X are called n-connected iff x and y are joined by an n-path 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 C VI n-connected 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 Cn [x,X].
To have a correspondence between the topology of the foreground and the background of VI - and thus to avoid connectivity paradoxes (Kong and Rosenfeld, 1989) - we must assign them with complementary neighborhoods. For instance, if we use N26 for the foreground of VI then we have to use N6 for the background and vice versa. In the sequel, we assign N26 to the foreground and N6 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 C VI 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 C VI are defined as
T26(x, X) = #C26 [x, N26(x) n X]	(5)
and
T(x,X) = #C6[x,Nj8(x) nX] ,	(6)
where #X denotes the cardinality of X. The point x is then simple (for X) iff T26(x,X) = 1 and T6(x,X) = 1. Note that we use Nj8 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 C VI. We call Y C VI 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 C VI, 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 X-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 Rn. We denote by
R(S) = min{r € R+ | 3x € Rn : Br(x) D S} (7)
the radius of the smallest ball enclosing S. Furthermore, let X be a subset of VI. The extended projection of x on X is given by
nX (x) =
u
% (y)
(8)
y€N6(x)nX :d(y,X )<d(x,X)
where
%(x) = {y € X | Vz € X : d(y,x) < d(z,x)} (9)
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 FX (x) = R(nX(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 FX with X > 0 yields the X-medial axis at level X. By definition all X-medial axes at level X' > X are included in the axis at level X. 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 Gartner (1999).
In our application we do not use the X-medial axis itself. We rather apply the radius projection map FX as priority function during the thinning process that automatically guides the process to the topological kernel of the several nested X-medial axes. That is, the smallest subset of all X-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 Fx is used (Passat etal., 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 DX. 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 Malandainetal. (1993), x can be classified topologically by evaluating T26(x,X) and T6(x,X). As the skeleton contains no surfaces, T6(x,X) is always one and therefore only T26(x,X) is meaningful. So x is part of a curve in X if T26(x, X) = 2 and part of a curve junction if T26(x, X) > 2, respectively. Note that Malandain etal. (1993) used Ni8 instead of N{8 in the definition of T6(x,X). However, this does not make a difference since Ni8(x) nX = Nf8(x) nX for any x € X c VI.
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 B1, the number of tunnels B2, and the number of cavities B3. Let x be a point of a skeleton X. Then the first Betti number of N|6(x) nX is given by B1 (x,X) = T26(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 6-components in Nj8(x) n X that contain at least one 6-neighbor of x. More precisely, the number of tunnels of x € X is defined as
B2(x, X )= 0	(10)
if #{N6 (x) n X } = 6 and
B2(x,X)= #{Y € C6 [x, Nig (x) nX] | Y nK(x)= 0} -1,
(11)
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 B2(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 X c VI by Nn(X) = UxeXNn(x). The set composed of all n-components of X c VI that are n-adjacent to Y c VI is denoted by Cn[Y,X]. Let X be a skeleton and C c X the set of all 26-components of X that form a curve junction. We now iterate over all curve junctions J c X and remove all simple points j from J (in decreasing order w. r. t. DX) for which the number of 26-components in N26(J) nX \ J' does not change, i.e.,
#C26 [J', N26(J) nx\ j'] = #C26 [J, N26(J) nX\ J] ,
(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.
rt
o
rt —j
O
,13 H
A
o rt
—J
s
¿3 H
o □
Curve point Surface point Curve junction point
Change to particular type after correction
II—				—II
				
				
				
•-				—II
II—				—II
				
				
				
II—				—II
	4	»	
	ft r	n i	
	' L	J '	
	i	»	
		»	
	■ ■	■ I	■ ii
• ■	I 1	I 1	■ •
	4	»	
Topological numbers
Betti numbers
Neighbor counting
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.
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 I. 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 x0
—d(x, x0) x < x0, %(x) = < 0	x = x0,
d(x, x0) x > x0
(13)
yielding a function pZ(%), with % € [—f, f] (see Fig. 3 for an illustration). The expected spherical contact profile of the entire foam is defined as p(% | I) = E(pZ(%) 11), that is the mean thickness at distance % of all struts Z with length I.
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 x0 of the struts' curve in the skeleton to the adjacent curve junctions x and y, that is l(Z) = d (x, xo)+ d (y, xo).
Dx (x)

from x to y whose normalized direction is denoted by r. Further, define two planes p1 and p2 that pass orthogonally to Lxy through x + DX(x)r and y — DX (y)r, respectively. The connecting strut segment is composed of all points of the strut between p1 and p2. See Fig. 5 for an illustration along with a visualization of a foam cell whose nodes were removed.
We define the volume correction cZ 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 e Y. Then cZ is computed as
Fig. 3: Illustration of the spherical contact profile for an individual strut.
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 pZ for a given curve segment Z by an individual volume correction factor cZ 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 Br (x) with radius r = DX (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 e Y is computed as follows: Let Lxy be the vector
cz =
v(C)
LzeZhx+Dx (x)r<z<y-Dx (y)r nDx(z)2 '
(14)
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 p1 and p2. We assume here, that the volume contribution of the nodes adjacent to C are of the same magnitude as for C.
r i
Xf-*
r 2
f—? y
i
(a)
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).
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 |j,m. This corresponds to 8.9 mm x 19.4 mm x 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).
H
Yt
!
Fig. 6: Cell of the aluminum foam and its skeleton.
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(E, 11) 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 15494 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(E, 11). 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 | I) depicted in Fig. 9. We subdivided the struts into the percentile ranges q^j = (xi; xj], where x, 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 q0,5 stands almost on top of the other ones. This indicates a nonlinear dependence of p(0 | I) on I as reported by Jang et al. (2008).
►HL,
(a)
Fig. 7: Structural pecularities found in the aluminum foam and their effect on the skeleton: Closed face (a) and eight-fold vertex (b)
( / ßm
-400	-200	0	200	400
Too 150 200 250 300 350^
p(f | I) / ßm
Fig. 8: Plot of the measured spherical contact profile p(% |1) for the entire aluminum foam.
Q0,5 Q5,10 Q10,20 Q20,30 Q30,40 q40,50 Q50,60 Q60,70 q70,80 Q80,90 Q90,100 Length percentile ranges
Fig. 9: Boxplot of the estimated mid-strut thickness p(0 | I) subdivided into its length percentile ranges qij. The median of p(0 | I) is depicted as gray line in the background.
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 | I) and strut length I. Let us denote by f = % / I the distance from the strut center % normalized by I. We then assume a polynomial of degree of at most eight as initial model:
pWf ) = a4|8 + as!6 + a2^4 + a^2 + 1 , (15)
a1,..., a4 € R, where the indices denote the exponents of f 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 f = 0 and thus an intercept of one was used.
Fig. 10 shows a plot of the measurements for the shortest (40,10), medial (445,55) and longest (490,100) 10% of the struts normalized by p(0 | I) and I 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( % | I) on I. Therefore, we subdivided the data into 10% length percentile ranges qij for our analysis.
Table 1 summarizes the results of the all-subset regression of p2468 guided by the MCp-statistic (Fujikoshi and Satoh, 1997). MCp is a modification of Mallows' (1973) Cp-statistic that was shown to be the minimum variance unbiased estimator of the expected overall Gauss discrepancy (Davies etal., 2006). Roughly speaking: The smaller the MCp-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 p-values of a t-test for the significance of the individual parameters in p2468 are given (null hypothesis H0 : ai = 0 against the alternative H1 : ai = 0). The p-values for a1 were left out as they are all zero and thus |2 is mandatory.
- P2468 fitted to q90,100
— P2468 fitted to 545,55 - P2468 fitted to qo,10
^-04 -02 0.0 02 04
f = i/t
Fig. 10: Plot of the measurements for the shortest (foreground points, black), medial (intermediate points, red) and longest (backgroundpoints, blue) 10% of the struts normalized by p(0 11) and I. Model p2468 was fitted individually to shortest, medial and longest struts using weighted least squares.
No general model for strut lengths smaller than x50 could be found, as in this region the central plateaus start to grow. We also observed, that only struts smaller than x30 could be reliably modeled by a two-parameter model. However, for all but q0,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 MCp-statistic. For the three-parameter models, the smallest
Table 1: Results of the all-subset regression analysis of model p2468. The best model for each qij is highlighted.
		p-value							MCp					
#0,10	0.01	0.01	0	51.49		9.96		6.97	5.39	9.95	8.02	8.39		
#10,20	0.23	0.37	0.19	105.53		5.25		-0.29	0.91	2.71	1.81	2.44		
q20,30	0	0.01	0.01	182.96		5.26		10.45	20.47	7.86	7.16	11.93		
#30,40	0.41	0.76	0.69	34.88		17.24		8.86	3.91	1.16	1.10	1.69		
#40,50	0.53	0.51	0.12	30.72		31.12		33.58	33.38	3.45	1.44		1.40	
#50,60	0	0.93	0.10	268.99		94.58		153.27	194.32	3.70	1.01	9.50		
#60,70	0	0.49	0.13	655.29	177.45			316.58	419.45	3.25	1.47	27.63		
#70,80	0	0.98	0.01	1500.08	251.10			530.68	759.55	7.15	1.00	35.94		
#80,90	0	0.65	0.03	2741.50	270.72			713.66	1106.94	5.57	1.21	58.42		
#90,100	0.42	0	0	5786.65	247.43			876.88	1616.74	94.69	59.40	1.65		
MCp-values were mainly obtained by p248. The only exception is q40,50, but with a negligible difference.
For strut lengths greater than x50, with the exception of q90,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 f6 instead of f4. 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:
p^f )= a^~2 + 1,	(16)
p6(f )= asf6 + 1,	(17)
p24(f )= a2f4 + a1 f2 +1.	(18)
All these models assume that the strut thickness scales with length. Here, however, p2 and p6 are disqualified: The first one by its high MCp-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 x50 its MCp-values were at least two times smaller than the ones reached by the other two-parameter models (with exception of q50,60 and q60,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 qi, j and compared its R2-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 I avoiding the consideration of different length classes is topic of our ongoing research.
Table 2: Parameters of the polynomial model p248 individually fitted to the length percentile ranges q^ j.
	a1	a2	a4
#0,10	1.4500	-0.2886	-19.4497
#10,20	2.5887	-0.7984	-24.3497
#20,30	3.2839	-2.3636	-10.1886
#30,40	3.0614	1.5318	-44.6902
#40,50	2.8883	3.7503	-57.7422
#50,60	2.3003	8.8237	-97.2233
#60,70	1.7494	12.8970	-126.1603
#70,80	1.2383	15.6366	-132.4259
#80,90	0.6412	18.8818	-136.8497
#90,100	0.3144	19.2832	-95.6320
ACKNOWLEDGMENTS
This work was supported by the Deutsche Forschungsgemeinschaft grant RE3002/1-1. We thank M. Couprie for providing the source code of his skeletonization algorithms on the Internet (Couprie, 2012) which was very useful to verify our own implementation.
REFERENCES
Attali D, Boissonnat JD, Edelsbrunner H (2009). Stability and computation of medial axes - a state-of-the-art report. In: Moller T, Hamann B, Russell RD, eds. Mathematical foundations of scientific visualization, computer graphics, and massive data exploration, mathematics and visualization. Math Visual 109-25. Banhart J (2001). Manufacture, characterisation and application of cellular metals and metal foams. Prog Mater Sci 46:559-632. Bertrand G (1994). Simple points, topological numbers and geodesic neighborhoods in cubic grids. Pattern Recogn Lett 15:1003-11. Bertrand G, Malandain G (1994). A new characterization of three-dimensional simple points. Pattern Recogn Lett 15:169-75.
Blum H (1967). A transformation for extracting new descriptors of shape. In: Models for the perception of speech and visual form. MIT Press.
Cenens V, Huis R, Chauvaux B, Dereppe JM, Gratin C, Meyer F (1994). 3d cellular structure characterization of flexible polyurethane foam. In: Kumar V, Seeler KA, eds., Cellular and microcellular materials, vol. 53. ASME, 29-44.
Chaussard J, Couprie M, Talbot H (2010). Robust skeletonization using the discrete A-medial axis. Pattern Recogn Lett 32:1384-94.
Chazal F, Lieutier A (2005). The A-medial axis. Graph Models 67:304-31.
Couprie M (2012). Euclidean skeletons. http://www.esiee. fr/^coupriem/es/. Accessed 10 June 2012.
Couprie M, Bertrand G (2009). New characterizations of simple points in 2D, 3D, and 4D discrete spaces. IEEE T Pattern Anal 31:637-48.
Couprie M, Coeurjolly D, Zrour R (2007). Discrete bisector function and Euclidean skeleton in 2D and 3D. Image Vision Comput 25:1543-56.
Davies SL, Neath AA, Cavanaugh JE (2006). Estimation optimality of corrected AIC and modified Cp in linear regression. Int Stat Rev 74:161-8.
Evako AV (2011). Characterizations of simple points, simple edges and simple cliques of digital spaces: One method of topology-preserving transformations of digital spaces by deleting simple points and edges. Graph Models 73:1-9.
Fraunhofer ITWM (2012). MAVI - Modular algorithms for volume images. http://www.mavi-3d.de/. Accessed 10 June 2012.
Fujikoshi Y, Satoh K (1997). Modified AIC and Cp in multivariate linear regression. Biometrika 84:707-16.
Gartner B (1999). Fast and robust smallest enclosing balls. In: Nesetril J, ed. Algorithms-ESA '99. Lect Not Comput Sci 1643:325-38.
Gibson LJ, Ashby MF (1999). Cellular Solids: Structure and properties. Cambridge University Press.
Green PJ, Silverman BW (1995). Nonparametric regression and generalized linear models. Chapman & Hall.
Jang WY, Kraynik A, Kyriakides S (2008). On the microstructure of open-cell foams and its effect on elastic properties. Int J Solids Struct 45:1845-75.
Jang WY, Kyriakides S (2009). On the crushing of aluminum open-cell foams: Part II Analysis. Int J Solids Struct 46:635-50.
Kanaun S, Tkachenko O (2006). Mechanical properties of open cell foams: Simulations by laguerre tesselation procedure. Int J Fracture 140:305-12.
Kanaun S, Tkachenko O (2008). Effective conductive properties of open-cell foams. Int J Eng Sci 46:551-71.
Kong TY, Rosenfeld A (1989). Digital topology: Introduction and survey. Comput Vision Graph 48:357-93.
Lautensack C, Giertzsch M, Godehardt M, Schladitz K (2008). Modelling a ceramic foam using locally adaptable morphology. J Microsc 230:396-404.
Li K, Gao XL, Subhash G (2006). Effects of cell shape and strut cross-sectional area variations on the elastic properties of three-dimensional open-cell foams. J Mech Phys Solids 54:783-806.
Liebscher A, Proppe C, Redenbach C, Schwarzer D (2012). Uncertainty quantification for metal foam structures by means of image analysis. Probabilist Eng Mech 28:143-51.
Liebscher A, Redenbach C (2012). 3D image analysis and stochastic modelling of open foams. Int J Mater Res 103:155-61.
Malandain G, Betrand G, Ayache N (1993). Topological segmentation of discrete surfaces. Int J Comput Vision 10:183-97.
Mallows CL (1973). Some comments on Cp. Technometrics 15:661-75.
Matheron G (1988). Example of topological properties of skeletons. In: Serra J, ed., Image analysis and mathematical morphology, vol. 2. Academic Press, 239-56.
Maurer CR, Qi R, Raghavan V (2003). A linear time algorithm for computing exact Euclidean distance transforms of binary images in arbitrary dimensions. IEEE T Pattern Anal 25:265-70.
Meijster A, Roerdink JBTM, Hesselink WH (2002). A general algorithm for computing distance transforms in linear time. In: Viergever M, Goutsias J, Vincent L, Bloomberg DS, eds. Mathematical morphology and its applications to image and signal processing. Comput Imaging Vision 18:331-40.
Montminy MD (2001). Complete structural characterization of foam using three-dimensional images. Ph.D. thesis, Department of Control Science and Dynamical Systems, University of Minnesota.
Nakamura A, Aizawa K (1985). On the recognition of properties of three-dimensional pictures. IEEE T Pattern Anal 7:708-13.
Ohser J, Schladitz K (2009). 3D images of materials structures: Processing and analysis. Wiley-VCH.
Otsu N (1979). A threshold selection method from gray-level histograms. IEEE T Syst Man Cyb 9:62-6.
Passat N, Couprie M, Bertrand G (2007). Topological
monsters in Z3: A non-exhaustive bestiary. In: Bannon GJF, Barrera J, Braga-Neto UM, Hirata NST, eds., Proc 8th Int Symp Math Morphol, vol. 2.
Passat N, Couprie M, Bertrand G (2008). Minimal simple pairs in the 3-D cubic grid. J Math Imaging Vis 32:239-49.
Redenbach C (2009). Modelling foam structures using random tessellations. In: Capasso V, Aletti G, Micheletti A, eds. Proc 10th Eur Congr Stereol Image Anal. Bologna: Esculapio.
Roberts AP, Garboczi EJ (2002). Elastic properties of model random three-dimensional open-cell solids. J Mech Phys Solids 50:33-55.
Rosenfeld A (1981). Three-dimensional digital topology. Inform Control 50:119-27.
Saha PK, Chaudhuri BB (1996). 3d digital topology under binary transformation with applications. Comput Vis Image Und 63:418-29.
Schouten TE, Kuppens HC, van den Broek EL (2006). Three-dimensional fast exact Euclidean distance (3D-FEED) maps. In: Proceedings of the Vision Geometry XIV, vol. 6066 of Proceedings ofSPIE.
Soille P (2002). Morphological image analysis. Heidelberg: Springer.
Talbot H, Vincent L (1992). Euclidean skeletons and conditional bisectors. In: Proc Visual Commun Image Proc. 1818:862-76.
Tekoglu C, Gibson LJ, Pardoen T, Onck PR (2011). Size effects in foams: Experiments and modeling. Prog Mater Sci 56:109-38.
van Merkerk R (2002). Fracture of metal foams. Master's thesis, Department of Applied Physics, Rijksuniversiteit Groningen.
Welzl E (1991). Smallest enclosing disks (balls and ellipsoids). In: Maurer H, ed. New results and new trends in computer science. Lect Not Comput Sci 555:359-70.