# Statistical analysis of the local strut thickness of open cell foams.

INTRODUCTIONFoams 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 ([micro]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 [micro]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.sup.+]. A three-dimensional binary image is defined as a subset VI of the cubic grid [Z.sup.3] together with a function f : [V.sub.I][right arrow]{0,1}. We refer to X [subset][V.sub.I] mapped to {1} as foreground and to its complement [bar.X], that is mapped to {0}, as background. A point x = ([x.sub.1], [x.sub.2], [x.sub.3])[member of][V.sub.I] with [x.sub.i][member of]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) = [square root of ([[summation].sub.i][([x.sub.i] - [y.sub.i]).sup.2])]. We denote the ball of radius r > 0 centered in x [member of] [V.sub.I] by

[B.sub.r](x) = {y [member of] [V.sub.I] | d(x, y)[less than or equal to]r}. (1)

The Euclidean distance of X [subset] VI and x [member of] VI is defined as d(x, X) = [min.sub.y[member of]X] d(x, y). The mapping [D.sub.X] (x) = d(x, [bar.X]) that assigns to each point x [member of] X its distance to the complement [bar.X] of X is called Euclidean distance map. Note that for all x [member of] X, the value [D.sub.X] (x) is the radius of the maximum ball centered in x that is completely contained in X. [D.sub.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 [member of] [V.sub.I] is defined as the set of all points in [V.sub.I] that are adjacent to x with an adjacency notion depending on some distance. In the remainder of this work, we consider the neighborhoods

[N.sub.6](x) = {y [member of] [V.sub.I] |[summation over (i)][absolute value of [x.sub.i] - [y.sub.i]] [less than or equal to] 1} (2)

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3)

that are induced by the city-block and the chessboard distance. The neighborhood obtained from [N.sub.26](x) by removing its outer corners is given by

[N.sub.18](x) = {y [member of] [V.sub.I]|[summation over (i)][absolute value of [x.sub.i] - [y.sub.i]] [less than or equal to]2} [intersection] [N.sub.26](x). (4)

We define [N.sup.*.sub.n](x) = [N.sub.n](x)\{x}. Note that the number of points in [N.sub.n](x) is n + 1.

Let x and y be two points in VI. We say that y is nadjacent to x iff y [member of] [N.sup.*.sub.n] (x). A set Y [subset] [V.sub.I] is n-adjacent to x [member of] [V.sub.I] iff there exists a point y [member of] Y which is n-adjacent to x. We call a sequence of points [x.sub.0], ..., [x.sub.k] in [V.sub.I] an n-path, if [x.sub.i-1] is n-adjacent to [x.sub.i] for all 1[less than or equal to]i[less than or equal to]k. Given a nonempty subset X [subset] [V.sub.I], two points x, y [member of] 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 [subset] [V.sub.I] 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 [C.sub.n][x, X].

To have a correspondence between the topology of the foreground and the background of [V.sub.I]--and thus to avoid connectivity paradoxes (Kong and Rosenfeld, 1989)--we must assign them with complementary neighborhoods. For instance, if we use [N.sub.26] for the foreground of [V.sub.I] then we have to use [N.sub.6] for the background and vice versa. In the sequel, we assign [N.sub.26] to the foreground and [N.sub.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 [subset] [V.sub.I] is called simple if its removal from X does not change the number of connected components and tunnels of X and [bar.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 [member of] X [subset] [V.sub.I] are defined as

[T.sub.26](x, X) = #[C.sub.26] [x, [N.sup.*.sub.26](x)[intersection]X] (5)

and

[T.sub.6](x,[bar.X]) = #[C.sub.6][x, [N.sup.*.sub.18](x)[intersection][bar.X]], (6)

where #X denotes the cardinality of X. The point x is then simple (for X) iff [T.sub.26](x,X) = 1 and [T.sub.6](x, [bar.X]) = 1. Note that we use [N.sup.*.sub.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 [subset] [V.sub.I]. We call Y [subset] [V.sub.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(nlog(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 [subset] [V.sub.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 [member of] X that have at least two closest points in [bar.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 [lambda]-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.sup.n]. We denote by

R(S) = min{r [member of] [R.sup.+]|[there exists]x [member of] [R.sup.n]: [B.sub.r](x)[contains or equal to]S} (7)

the radius of the smallest ball enclosing S. Furthermore, let X be a subset of VI. The extended projection of x on [bar.X] is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (8)

where

[[PI].sub.[bar.X]](x) = {y [member of] [bar.X]|[for all]z [member of] [bar.X] : d(y,x)[less than or equal to]d(z,x)} (9)

denotes the plain projection of x [member of] X on the complement [bar.X] of X. Couprie et al. (2007) then define the projection radius map as the function [F.sub.X] (x) = R([[PI].sup.e.sub.[bar.X]](x)) that assigns to each point x [member of] X the radius of the minimum ball enclosing all extended projections of x on [bar.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.sub.X] with [lambda] > 0 yields the [lambda]-medial axis at level [lambda]. By definition all [lambda]-medial axes at level [lambda]' > [lambda] are included in the axis at level [lambda]. 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 [lambda]-medial axis itself. We rather apply the radius projection map [F.sub.X] as priority function during the thinning process that automatically guides the process to the topological kernel of the several nested [lambda]-medial axes. That is, the smallest subset of all [lambda]-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.sub.X] is used (Passatetal., 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.sub.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 Malandainetal. (1993), x can be classified topologically by evaluating [T.sub.26](x, X) and [T.sub.6](x, [bar.X]). As the skeleton contains no surfaces, [T.sub.6](x, [bar.X]) is always one and therefore only [T.sub.26](x, X) is meaningful. So x is part of a curve in X if [T.sub.26](x, X) = 2 and part of a curve junction if [T.sub.26](x, X) > 2, respectively. Note that Malandain etal. (1993) used [N.sub.18] instead of [N.sup.*.sub.18] in the definition of [T.sub.6](x, [bar.X]). However, this does not make a difference since [N.sub.18(x)][intersection]X = [N.sup.*.sub.18](x)[intersection][bar.X] for any x [member of] X [subset] [V.sub.I].

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.sub.1], the number of tunnels [B.sub.2], and the number of cavities [B.sub.3]. Let x be a point of a skeleton X. Then the first Betti number of [N.sup.*.sub.26](x)[intersection]X is given by [B.sub.1] (x,X) = [T.sub.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.sup.*.sub.18](x)[intersection][bar.X] that contain at least one 6neighbor of x. More precisely, the number of tunnels of x [member of] X is defined as

[B.sub.2](x, X)= 0 (10)

if #{[N.sup.*.sub.6](x)[intersection]X} = 6 and

[B.sub.2](x,X)= #{Y [member of] [C.sub.6] [x, [N.sup.*.sub.18] (x)[intersection][bar.X]] | Y[intersection] [N.sup.*.sub.6] (x)[not equal to] 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 [B.sub.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 X [subset] [V.sub.I] by [N.sub.n](X) = [[union].sub.x[member of]]X [N.sub.n](x). The set composed of all n-components of X [subset] [V.sub.I] that are n-adjacent to Y [subset] [V.sub.I] is denoted by [C.sub.n][Y,X]. Let X be a skeleton and C [subset] X the set of all 26-components of X that form a curve junction. We now iterate over all curve junctions J [subset] X and remove all simple points j from J (in decreasing order w. r. t. [D.sub.X]) for which the number of 26-components in [N.sub.26](J')[intersection]X\J' does not change, i.e.,

#[C.sub.26] [J, [N.sub.26](J)[intersection]X\J'] = #[C.sub.26] [J, [N.sub.26](J)[intersection]X\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.

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 [member of] Y which is assumed to be of length l. The local strut thickness at point x [member of] 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.sub.0]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (13)

yielding a function [p.sub.Z]([xi]), with [xi] [member of] [-l/2, l/2] (see Fig. 3 for an illustration). The expected spherical contact profile of the entire foam is defined as p([xi]| I) = E([p.sub.Z]([xi]) | l), that is the mean thickness at distance [xi] 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, [x.sub.0]) + d (y, [x.sub.0]).

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.sub.Z] for a given curve segment Z by an individual volume correction factor [c.sub.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.sub.r] (x) with radius r = [D.sub.X] (x) centered in the curve junction x.

The strut segment between any two junctions x, y [member of] Y is computed as follows: Let [L.sub.xy] be the vector from x to y whose normalized direction is denoted by r. Further, define two planes [p.sub.1] and [p.sub.2] that pass orthogonally to [L.sub.xy] through x + [D.sub.X](x)r and y - [D.sub.X](y)r, respectively. The connecting strut segment is composed of all points of the strut between [p.sub.1] and [p.sub.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.sub.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 [member of] Y. Then [c.sub.Z] is computed as

[c.sub.z] = v(C)/[[summation].sub.z[member of]Z[LAMBDA]x + [D.sub.X](x)r<z<y-[D.sub.x](y)r][pi][D.sub.X][(z).sup.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 [p.sub.1] and [p.sub.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 [mu]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).

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([xi]|l) 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 [bar.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([xi]|l). It shows the typical flat form in the strut center ([xi] = 0) which steeply rises towards the nodes (increasing absolute value of [xi]).

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 | l) depicted in Fig. 9. We subdivided the struts into the percentile ranges [q.sub.i,j] = ([x.sub.i], [x.sub.j]], where [x.sub.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.sub.0,5] stands almost on top of the other ones. This indicates a nonlinear dependence of p(0 | l) on l as reported by Jang et al. (2008).

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 | l) and strut length l. Let us denote by [??] = [xi]/ l the distance from the strut center [xi] normalized by l. We then assume a polynomial of degree of at most eight as initial model:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15)

[a.sub.1], ..., [a.sub.4][member of]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.sub.0,10]), medial ([q.sub.45,55]) and longest ([q.sub.90,100]) 10% of the struts normalized by p(0 | l) and l along with the initial model [[??].sub.2468] fitted to the individual profiles. The resulting curves differ heavily from each other, which indicates a dependence of p([xi]| l) on l. Therefore, we subdivided the data into 10% length percentile ranges [q.sub.i,j] for our analysis.

Table 1 summarizes the results of the all-subset regression of [[??].sub.2468] guided by the [MC.sub.p]-statistic (Fujikoshi and Satoh, 1997). [MC.sub.p] is a modification of Mallows' (1973) [C.sub.p]-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 [MC.sub.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 p-values of a t-test for the significance of the individual parameters in [[??].sub.2468] are given (null hypothesis [H.sub.0] : [a.sub.i] = 0 against the alternative [H.sub.1]: [a.sub.i][not equal to]0). The p-values for [a.sub.1] were left out as they are all zero and thus [[??].sup.2] is mandatory.

No general model for strut lengths smaller than [x.sub.50] could be found, as in this region the central plateaus start to grow. We also observed, that only struts smaller than [x.sub.30] could be reliably modeled by a two-parameter model. However, for all but [q.sub.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.sub.p]-statistic. For the three-parameter models, the smallest [MC.sub.p]-values were mainly obtained by [[??].sub.248]. The only exception is [q.sub.40,50], but with a negligible difference.

For strut lengths greater than [x.sub.50], with the exception of [q.sub.90,100], [[??].sub.248] 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 [[??].sup.6] instead of [[??].sup.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 etal. (2008, Eq. 1)--Eq. 18) the following subsets of [[??].sub.2468] have already been proposed to describe the varying strut thickness:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (17)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (18)

All these models assume that the strut thickness scales with length. Here, however, [[??].sub.2] and [[??].sub.6] are disqualified: The first one by its high [MC.sub.p]-values and the second one by the mandatory but missing [[??].sub.2]. Among these models only [[??].sub.24] showed a reasonable performance. For strut lengths greater than [x.sub.50] its [MC.sub.p]-values were at least two times smaller than the ones reached by the other two-parameter models (with exception of [q.sub.50,60] and [q.sub.60,70] for [[??].sub.26] 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 [[??].sub.248] 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.sub.i,j] and compared its [R.sup.2]-statistics with the ones of the initial parametric model [[??].sub.2468]. 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 [micro]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 [xi] and l avoiding the consideration of different length classes is topic of our ongoing research.

doi: 10.5566/ias.v32.p1-12

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 graylevel 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 (3DFEED) 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.

ANDRE LIEBSCHER [mail] and CLAUDIA REDENBACH

University of Kaiserslautern, Department of Mathematics, Erwin-Schrodinger StraBe, 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)

* The topic of this paper was presented at the S4G Conference, June 25-28, 2012 in Prague, Czech Republic.

Table 1: Results of the all-subset regression analysis of model [[??].sub.2468]. The best model for each [q.sub.i,j] is highlighted. p-value [MC.sub.p] [a.sub.2] [a.sub.3] [a.sub.4] [[??].sub.2] [q.sub.0,10] 0.01 0.01 0 51.49 [q.sub.10,20] 0.23 0.37 0.19 105.53 [q.sub.20,30] 0 0.01 0.01 182.96 [q.sub.30,40] 0.41 0.76 0.69 34.88 [q.sub.40,50] 0.53 0.51 0.12 30.72 [q.sub.50,60] 0 0.93 0.10 268.99 [q.sub.60,70] 0 0.49 0.13 655.29 [q.sub.70,80] 0 0.98 0.01 1500.08 [q.sub.80,90] 0 0.65 0.03 2741.50 [q.sub.90,100] 0.42 0 0 5786.65 [MC.sub.p] [[??].sub.24] [[??].sub.26] [[??].sub.28] [q.sub.0,10] 9.96 6.97 5.39# [q.sub.10,20] 5.25 -0.29# 0.91 [q.sub.20,30] 5.26# 10.45 20.47 [q.sub.30,40] 17.24 8.86 3.91 [q.sub.40,50] 31.12 33.58 33.38 [q.sub.50,60] 94.58 153.27 194.32 [q.sub.60,70] 177.45 316.58 419.45 [q.sub.70,80] 251.10 530.68 759.55 [q.sub.80,90] 270.72 713.66 1106.94 [q.sub.90,100] 247.43 876.88 1616.74 [MC.sub.p] [[??].sub.246] [[??].sub.248] [[??].sub.268] [q.sub.0,10] 9.95 8.02 8.39 [q.sub.10,20] 2.71 1.81 2.44 [q.sub.20,30] 7.86 7.16 11.93 [q.sub.30,40] 1.16 1.10# 1.69 [q.sub.40,50] 3.45 1.44 1.40# [q.sub.50,60] 3.70 1.01# 9.50 [q.sub.60,70] 3.25 1.47# 27.63 [q.sub.70,80] 7.15 1.00# 35.94 [q.sub.80,90] 5.57 1.21# 58.42 [q.sub.90,100] 94.69 59.40 1.65# Note: The best model for each [q.sub.i,j] is highlighted are indicated with #. Table 2: Parameters of the polynomial model [[??].sub.248] individually fitted to the length percentile ranges [q.sub.i,j]. [a.sub.1] [a.sub.2] [a.sub.4] [q.sub.0,10] 1.4500 -0.2886 -19.4497 [q.sub.10,20] 2.5887 -0.7984 -24.3497 [q.sub.20,30] 3.2839 -2.3636 -10.1886 [q.sub.30,40] 3.0614 1.5318 -44.6902 [q.sub.40,50] 2.8883 3.7503 -57.7422 [q.sub.50,60] 2.3003 8.8237 -97.2233 [q.sub.60,70] 1.7494 12.8970 -126.1603 [q.sub.70,80] 1.2383 15.6366 -132.4259 [q.sub.80,90] 0.6412 18.8818 -136.8497 [q.sub.90,100] 0.3144 19.2832 -95.632

Printer friendly Cite/link Email Feedback | |

Title Annotation: | Original Research Paper |
---|---|

Author: | Liebscher, Andre; Redenbach, Claudia |

Publication: | Image Analysis and Stereology |

Article Type: | Report |

Geographic Code: | 4EXCZ |

Date: | Mar 1, 2013 |

Words: | 7922 |

Previous Article: | Extraction of curved fibers from 3d data. |

Next Article: | Spatial dependence of random sets and its application to dispersion of bark beetle infestation in a natural forest. |

Topics: |