Image Anal Stereol 2013;32:45-55 Original Research Paper
IMPROVED ESTIMATION OF FIBER LENGTH FROM 3-DIMENSIONAL
*
IMAGES*
Joachim Ohser^'1, Konrad Sandau2, Jürgen Kampf3, Irene Vecchio4 and Ali Moghiseh5
1Univ. Appl. Sci. Darmstadt, Dept. Math & Nat. Sci., Schofferstr. 3, D-64295 Darmstadt, Germany; 2Univ. Appl. Sci. Darmstadt, Dept. Math & Nat. Sci., Schofferstr. 3, D-64295 Darmstadt, Germany; 3Techn. University Kaiserslautern, Dept. Math., PB 3049, D-67653 Kaiserslautern; 4Fraunhofer ITWM Kaiserslautern, Fraunhofer-Platz 1, D-67663 Kaiserslautern, Germany; 5Univ. Appl. Sci. Darmstadt, Dept. Math & Nat. Sci., Schofferstr. 3, D-64295 Darmstadt, Germany
e-mail: joachim.ohser@h-da.de, konrad.sandau@h-da.de, kampf@mathematik.uni-kl.de, irene.vecchio@itwm.fhg.de, ali.moghiseh@h-da.de
(Received August 27, 2012; revised December 22, 2012; accepted January 7, 2013) ABSTRACT
A new method is presented for estimating the specific fiber length from 3D images of macroscopically homogeneous fiber systems. The method is based on a discrete version of the Crofton formula, where local knowledge from 3 x 3 x 3-pixel configurations of the image data is exploited. It is shown that the relative error resulting from the discretization of the outer integral of the Crofton formula amounts at most 1.2 %. An algorithmic implementation of the method is simple and the runtime as well as the amount of memory space are low. The estimation is significantly improved by considering 3 x 3 x 3-pixel configurations instead of 2 x 2 x 2, as already studied in literature.
Keywords: Crofton formula, image analysis, integral of mean curvature, systematic error.
INTRODUCTION
Fibers in fiber-reinforced materials, fleeces, felts and paper form systems of non-overlapping and more or less curved fibers. The characterization of these fiber systems is a topic of materials science and technology, where quantities estimated from image data are used in industrial quality control, as parameters for geometric models of fiber systems, and as input for simulating materials properties. Currently, the development in fiber system characterization is mainly driven by computed microtomography (uCT) which has been established as an appropriate source of three-dimensional (3D) image data of sufficiently high lateral resolution and contrast between fibers and the matrix constituents or air in pore spaces (Weitkamp etal., 2011). Often, simple thresholding (e.g., using a global threshold level or a hysteresis) can be used to segment the image data, where the fiber system is assigned to the foreground and the matrix is background. Throughout the present article we assume that the foreground pixels form a sampling of the fiber system on a cubic primitive point lattice.
Besides the fiber direction distribution, the length density of the fibers is probably the most important characteristic of a macroscopically homogeneous (i.e., stationary) fiber system. An intuitive method of
estimating fiber lengths is based on skeletonization of the fibers and counting the number of skeleton pixels, where the counts are weighted with the lattice distance a if two neighboring pixels are connected along an edge of the cubic lattice cells, with \/2a if the pixels are connected via a face diagonal, and with V3a if the pixels are connected via a space diagonal, respectively. This intuitive estimator can be seen as a 3D analogue of the boundary length estimation based on the Freeman encoding of a 2D object. It is not surprising that, depending on the adjacency chosen for the pixels of the 3D image, the systematic error of the estimated fiber length can be considerably large, in particular in the isotropic case.
One should keep in mind that the length of a fiber (with small cross-section) is approximately proportional to the integral of the mean curvature of the fiber. Furthermore, the integral of the mean curvature of a convex 3D body is (up to a constant factor) one of the four intrinsic volumes. A further one is (up to a constant factor) the surface area. Therefore, the principle behind the estimation of fiber length is conceptionally the same as for the surface area.
First, we refer to Lindblad and Nystrom (2002), who suggested to estimate the surface area of a body as the sum of the areas of the surface patches obtained
* The topic of this paper was presented at the S4G Conference, June 25-28, 2012 in Prague, Czech Republic.
by ordinary surface rendering (Ohser and Schladitz, 2009, sec. 3.6.2). In this approach, the areas of the surface patches serve as weights for computing the surface area from local knowledge. Analogously, the integral of the mean curvature (and thus the fiber length) can be estimated from a tessellation of the foreground set into polytopes induced by rendering data. The corresponding weights are computed via the inclusion-exclusion principle and Hadwiger's formula for the mean width of polytopes (e.g., Hadwiger, 1957, p. 215; Santalo, 1976, p. 226). However, this estimator is not multigrid convergent, even not in the isotropic case, i.e., for decreasing lattice distance the estimated integral of the mean curvature does not converge to the true value.
There are alternative approaches of estimating curve length (or circumference of a 2D set) based on an explicit approximation of the curve (or the boundary of the set) but ensuring multigrid convergence, see Klette and Rosenfeld (2004) for an overview. In Klette etal. (1999) digital curves are approximated by piecewise straight lines. However, such approximations are computationally demanding since they do not claim locality.
The surface area can be measured directly from the image data without need to approximate the surface (e.g., by a surface rendering). For each 2 x 2 x 2-pixel configuration we count its occurrences in the image and approximate the surface area by a linear combination of these counts. There are several proposals how to choose the coefficients, called weights, of the linear combination. The weights suggested by Lindblad (2005) minimize the estimation variance of the surface area of a plane with random normal direction uniformly distributed on the unit sphere. This idea goes back to Mullikin and Verbeek (1993), see also the discussion in Windreichetal. (2003). A further approach is presented in Chapter 4 of Ohser and Mucklich (2000) and in Lang et al. (2001) where the weights are computed using a discrete version of one of Crofton's intersection formulas (Schladitz et al., 2006a; Ohser and Schladitz, 2009, Chapter 5). We remark that the method of Ohser and Mucklich (2000) is designed for the general case of cuboidal lattices, while the restriction on the particular case of cubic primitive lattices allows to exploy symmetry properties (Ohser and Schladitz, 2009), and to present the weights in a very condensed form (depending on representatives of the 22 equivalence classes, as in Table 1).
A comprehensive treatment of the subject of surface area estimation is given in Ziegel and Kiderlen (2010). The weights suggested in this article minimize the worst case asymptotic error for surface area
estimation where asymptotics is understood with respect to decreasing lattice distances. This approach is based on a general asymptotic result shown in Kiderlen and Rataj (2006). It also allows a comparison of the various methods of surface area estimation. Unfortunately, Ziegel's methods can not be extended to the estimation of the integral of the mean curvature, since the estimator based on optimal weights completely ignores the image data (Kampf, 2012). In Svane (2012) the weights for estimating the intrinsic volumes are designed in such a way that their estimation is unbiased for Boolean models or isotropic random closed sets fulfilling certain smoothness assumptions.
Further approaches of estimating the intrinsic volumes are based on solving systems of linear equations deduced from the extended Steiner formula or the principal kinematic formula (Schmidt and Spodarev, 2005; Klenk etal., 2006; Mrkvicka and Rataj, 2008; 2009; Meschenmoser and Spodarev, 2012). These methods are studied in detail for the 2D case but work in principle in arbitrary dimensions. However, so far none of these algorithms has been proven to work in practice for dimensions > 3. The authors are claiming multigrid convergence for the resulting estimators of the intrinsic volumes. Comparisons for 2D Boolean models in Guderlei et al. (2007) and Mrkvicka and Rataj (2008) show, that the accuracy is sometimes higher than that of the estimators given in Section 4.2 of Ohser and Mucklich (2000).
We follow the approach of Ohser and Schladitz (2009) based on a Crofton formula which boils down the estimation of the integral of the mean curvature to computing Euler numbers in virtual planar sections. Discretisation of this formula combined with an efficient calculation of the Euler numbers in the planar sections yields a fast algorithm for determining the integral of the mean curvature from tomographic images, see, e.g., Vogel et al. (2010) for an application in the field of characterizing soil structure. The backbone of the Euler number calculation are thorough investigations of digital connectivity and consistency from Nagel et al. (2000); Ohser et al. (2002; 2003). In the present article it will be shown that the method of Ohser and Schladitz (2009) can be improved considerably by finer discretization of the outer integral in the Crofton formula.
The above approach has certain advantages over others. First, the method can simply be extended to arbitrary homogeneous lattices so that we are no longer restricted to cubic primitive lattices (Ohseretal., 2009; Ohser and Schladitz, 2009). This is an important fact since many 3D image acquisition techniques
produce images on non-cubic lattices. Furthermore, our approach can be extended to cases of hollow fibers (or partially hollow fibers, e.g., kinds of cellulose fibers) where a high fraction of the cross-section is not simply connected. In fact, the computation of the Euler number in virtual planar sections can be replaced by counting connected components or tangent points, as suggested in Schladitz et al. (2006b). Finally, we remark that the approach presented in this article has the capability to involve the estimation of the fiber direction distribution (Ohser and Schladitz, 2009, sec. 5.4).
The present article is organized as follows: First we recall the principle of estimating of the integral of the mean curvature from numbers of 2 x 2 x 2-pixel configurations. Then we compute the systematic error related to the discretization of the outer integral of Crofton's intersection formula. It will be shown that the systematic error can be reduced considerably by a finer discretization of the outer integral. Finally, we compute the optimal weights for the quadrature rule of the applied numerical integration.
H1 exists for each point on the surface dy of y and, using the 2-dimensional Hausdorff measure H2, the integral of the mean curvature is defined as
M(y) = J Hi(s) dH 2(s)
dy
From Theorem 5.9 in Federer (1959) it follows that M(y) ^ nl(0) as r 10. In other words, M(y)/n is an appropriate estimate of the fiber length l as r ^ l.
We remark that in the special case K = Br, the fibers are similar to the tubular neighborhood of 0 considered in Baddeley and Averback (1983). The only difference is that the fibers end with half-balls, while the tubular neighborhoods are cut off by the planes through the endpoints of the curve.
In order to apply Crofton's formula (Schneider, 1993, Theorem 4.5.5) we first show that every fiber can be approximated by a polyconvex set. Let us consider a piecewise straight approximation 0m of the curve 0 by the segments [f ((k — 1)£/m), f (kl/m)] between the curve points f ((k — 1)£/m) and f (kl/m),
FIBERS AND FIBER LENGTH
In our setting a fiber y in r3 with a convex cross-section of nonempty interior can be defined as follows: Let f : r ^ r3 be an arclength parametrization of a finite immersed curve 0 C r3 with f (0) = f (l), where f is twice continuously differentiable and l is the curve length,
0 = {f (s) : 0 < s < l} .
The first derivative of f is the unit tangent vector of the curve 0 at the point f (s). More precisely, let S+ denote the positive unit half-sphere of r3, then the non-orientated direction of the curve at the point f (s) is defined as d(s) = f(s) if f '(s) € S+, and 6(s) = —f'(s) otherwise.
Let Br denote the 3D ball of radius r centered at the origin, and let © be the Minkowski addition. The curve 0 may be shaped in such a way that the parallel set 0 ©Br of 0 is morphologically regular, i.e., there is an e > 0 such that 0 © Br is morphologically closed as well as morphologically open with respect to Be. A necessary condition for the morphological regularity of 0 © Br is that 0 is smooth enough, more precisely, the curvature k(s) = || f"(s) || of 0 must be smaller than the inverse radius, k(s) < 1/r for all s € [0, l].
Let K c Br be a compact, convex and morphologically regular set. Then a fiber y is defined as the set y = 0 © K. As a consequence of the morphological regularity of K, the mean curvature
=U [f {(k - 1)l/m), f(k£/m)} , m > 1 , k=i
fulfilling £(0m) ^ £($) as m ^ The set ym = 0m © K is polyconvex and it follows that M(ym) ^ n£(0) as m ^ ^ and r 10.
Let now L2 be the set of 2-dimensional linear subspaces in r3 (i.e., the set of planes hitting the origin). By ±L, v±L and ß we denote the orthogonal space of a plane L G L2, the Lebesgue measure on ±L, and the rotation invariant probability measure on the unit sphere S2 of r3, respectively. Then, using Crofton's intersection formula, Schneider (1993), the integral of the mean curvature M(X) of a polyconvex set X can be written in the form
M(X) = 2nJ Jx(Xn (L + x)) v±L(dx)p(dL), (1)
L2 XL
p(X, L)
where x(X n (L+x)) is the Euler number of the planar section X n (L + x) of X and the inner integral p(X, L) can be seen as the length of the orthogonal projection of X onto the straight line ±L.
ESTIMATION OF FIBER LENGTH
The sampling of the set X on a homogeneous point lattice implies a discretization of the Crofton
m
formula (Eq. 1), where the plane L is replaced with a 2-dimensional section lattice, the translations L + x over the orthogonal space lL are replaced with a translation of the section lattice over the corresponding translation lattice, and the integrals are numerically computed by quadrature rules.
Let l3 = az3 be a homogeneous lattice in r3 with the lattice distance a > 0, the unit cell C = [0,a]3 and the set {0, a}3 of vertices of C. We consider a plane L G L2. Assume now that there is a translation y G {0, a}3 such that the intersection of the translated plane hits at least three vertices of C,
#((L + y) n{0,a}3) > 3 .
Then l2 = L n l3 forms a homogeneous 2-dimensional section lattice of l3. This means that there are linearly independent vectors u, v G r3 with u + y, v + y G {0, a}3 generating the lattice l2,
l2 = {iu + jv : i, j G z} .
Without loss of generality we assume that the normal direction 0 = u x v/||u x vll is on S+ of
r3. Furthermore, there exists a vector w G r3 with | det(u, v, w)| = a3 generating a 1-dimensional section lattice ^l2 = {kw : k G z} with the property that l3 is the union of l2 + z over all translations z G ^l2,
l3 = J (l2 + z) ,
zG±L2
(Ohser and Schladitz, 2009). The lattice ^l2 is called the translation lattice of l2. The above results are summarized in the following lemma:
Lemma 1 Let be given a lattice l3 = az3, a plane L G L2 and a vector y G {0, a}3 such that #((L + y) n {0, a}3) > 3. Then there are three linearly independent vectors u,v,w G {-a,0,a}3 with l3 = {iu + jv + kw : i, j, k G z} and l3 n L = {iu + jv : i, j G z}.
As has been pointed out in Ohser and Schladitz (2009), p. 157, there exist m = 13 section lattices l|, k = 1,..., m, with pairwise different normal directions 0k of the corresponding section planes Lk = span l|. Let wk be the basis vector of translation lattice of l2. Then dk = |0kwk | is the orthogonal distance of the neighboring section planes Lk and Lk + wk. Notice that dk = |C|/|Ck|, where |C| is the volume of C and |Ck| denotes the area of the unit cell Ck = [0, uk] © [0, vk] of lf. Now, using a simple rectangular quadrature rule, the inner integral p(X, L) of the Crofton formula (Eq. 1) can be approximated with
p(X,Lk) » dk £ x{Xn (Lk + z)) ,
where the dk serve as the weights of the quadrature rule (Ohser and Schladitz, 2009, p. 154).
When using a 2-dimensional analog of the rectangular quadrature rule, Eq. 1 yields
M(X) » £jkp(X, Lk)
(3)
k=1
(Ohser and Schladitz, 2009, p. 156), and summarizing Eq. 2 and Eq. 3, one obtains an approximation of the integral of the mean curvature,
M(X) » £ Ykdk £ X{X n (Lk + z))
(4)
k=1
ze±L2
(2)
The weights of the quadrature rule applied in Eq. 3 are defined as follows: Consider the Voronoi tessellation of the unit sphere S2 with respect to the point field {0i,...,dm, — 0i,..., —0m}. By yk we denote the area of the Voronoi cell with respect to the direction 0k, k = 1,..., m. Numerical calculation yields yk = 0.575263 for the three directions parallel to the edges of the unit cell C, Yk = 0.464711 for the directions of the face diagonals of C, and yk = 0.442280 for the directions of the space diagonals of C.
In order to give an estimator of the Euler number %{Xn (Lk + z)) we consider a 6-adjacency system fk on the 2-dimensional section lattice l|. A 6-adjacency system is generated by a tessellation of the unit cell Ck into two closed triangles Tk1 and Tk,2 with Tk1 U Tk,2 = Ck. Formally, this can be done as follows: Let conv Z denote the convex hull of a subset Z ^ F°(Ck), where F0(Ck) denotes the set of vertices of Ck. Then the local adjacency system f0 C {conv Z : Z Q F°(Ck)} with respect to the tessellation mentioned above consists of the unit cell Ck, the triangles Tk1 and Tk,2, their edges, their vertices, and the empty set. Finally, the adjacency fk is the union of f0 + x over all lattice translations,
fk = J (f0 + x) ,
xehf
(Ohser and Schladitz, 2009, p. 50).
Now we introduce a discretization X n fk of the planar section X n span lk with respect to the adjacency system fk,
X n fk = U {F G fk : F0 (F) C X } .
It is important to see that the discretization X n fk can be obtained from the sampling X n lk2 of X on the section lattice l|, that is X n fk and X n l| carry the same information on X. The additivity of the Euler
number ensures that it can be computed locally. Using the Euler-Poincare formula we get
2
x(x nfk) = £ (-1)j £ 1(F0(F) c x), (5)
j=0	F eFj(Fk)
(Ohser and Schladitz, 2009, Eq. 3.5). Here Fj(fk) denotes the set of the j-dimensional polytopes belonging to fk, and 1 is the indicator function.
The Euler number x(X n fk) can be seen as an approximation of x(X n Lk), i.e., x(X n fk) ^ X(X n Lk) as a 10 (Ohser et al., 2002). Finally, when replacing the expression x (X n Lk) in Eq. 4 with x (X n fk) one gets
m
M(X) » £ Ykdk £ XX n (fk + z)) . (6)
k=1
ze±L2
The right-hand side is a discretization of the Crofton formula induced by the sampling of X on l3.
ALGORITHMIC IMPLEMENTATION
Our aim is now to present a factorization of the estimator Eq. 6 in the form
M(X ) « aqh ,
where the vector q depends on the kind of discretization of the Crofton formula and h is the vector of 2 x 2 x 2-pixel configurations in a binary image with the foreground pixels X n L3.
As in the previous section, we consider a 6-adjacency system fk on l2 generated from a
k2
tessellation of the unit cell Ck of l2 into two triangles Tk,i and Tk,2. Then the system T = {T1,T2,...} of all triangles belonging to fk forms a face-to-face tessellation of Lk, i.e., for i = j the intersection Tin Tj is either an edge of T, a vertex of Ti or empty. Notice that there are two possible tessellations of Ck. For each index k we choose one of them arbitrarily.
For each triangle T belonging to the tessellation T we define fo(T) = U2=0 Fj(T) and X n fo(T) = U{F G fo(T) : Fo(F) C X}. Then the (local) contribution Xo(X, T) of X n fo(T) to the Euler number can be expressed as
*0(X, T ) = £
(-iy = (3 - j)!
#Fj (X n f0(T ))
(7)
1
6
1,
-1, 0,
if #(F 0(T ) n X ) = 1,
if #(f 0(t ) n x ) = 2,
otherwise,
where the weights 1/(3 — j)! in Eq. 7 are due to the fact that each vertex of X n fo(T) hits 6 triangles of T and each edge hits 2 triangles, see Fig. 1 for a sketch. Notice that Eq. 7 is a special case of Eq. 3.io in Ohser and Schladitz (2oo9). Now, using the inclusion-exclusion principle for additive functionals one can prove the following lemma:
Lemma2 Let T = {Ti,T2,...} be the tessellation of Lk generated from a 6-adjacency system fk. Then E/=1 Xo(K, TÎ) = 1 for all compact and convex sets K' c Lk with the property that K' n fk is connected and nonempty.
It follows that the Euler-Poincare formula (Eq. 5) can be rewritten as
X(Xnfk) = £ (Xo(X — y, Tk,i)+ Xo(X — y, Tk,2)) .
yGL2
Furthermore, due to the local structure of this formula, it is sufficient to consider pixel configurations. By £o,...,^n—1 C {o,a}3 we denote the n = 255 local configurations on l3, where = {o, a}3 \ fy is the complementary configuration of . The numbers hl of occurrences of the in the sampling X n l3 are
hi = £ 1 (fy + x C X) ■ l(fylc + x C Xc) ,
xeL3
for I = o,..., n — 1. Let yk;i G r3 be a translation of the triangle Tki such that F°(Tk)i + yk>i) c {o,a}3. With the coefficients
1 ■ 6
2 ■ 1 -1 ■ 2
3 ■ 6 - 3 ■ 2 + 1
Fig. 1. A sketch of the computation for the local contribution Xo(X, T) of X n fo(T) to the Euler number depending on the number of vertices of T hitting X (full discs).
Pi (Tk,i)= X0(^i, Tk,i + yk,i) , i = 1,2 , and using the convention 0 • ~ = 0, one obtains £ x(xn(fk+z))
ze±L2
£ (^0(X - x, Tk,1)+ ^0(X - x, Tk,2))
xeL3
^ 1
£ CPi(Tk,1)+ Pl(Tk,i))hi . l=0
0
Finally, we get
M(X) =
m	n—1
E Ykdk E (pi(Tk,i)+ Pi(Tw))hi k=1 l=0 n—1
= a E qihi ,	(8)
l=0
as an estimator of M(X) with the weights
1 m
qi = a E Ykdk{pi(Tk,i)+ Pi(Tk,2)) .
k=1
It is important to see that dk/a does not depend on a. Eq. 8 shows that the integral of the mean curvature M(X) of X can be estimated from the sampling X n l3 as a scalar product qh of the vectors q = (qi ) and h = (hi ) where q is independent of both the data and the lattice distance a. Numerical values for the qi are given in Table 1, where the pictograms (full disc for foreground, empty disc for background) stand for the representatives of the 22 equivalence classes of the 2 x 2 x 2-pixel configurations fy. In this context, two pixel configurations are said to be equivalent if there is a rotation or a reflection such that one of them can be transformed to the other.
Table 1. The weights qi for the estimation of the integrai of the mean curvature M of a set X sampied on a homogeneous iattices l3 = az3 (Ohser and Schiaditz, 2009, p. 160).
conf.
qi
conf.
qi
0
0.589 849 0.727 768 0.616 231 0.686 796 0.446 035 0.425 548 0.334497 0 0 0
0 0 0
-0.334497 -0.425 548 -0.446035 -0.686796 -0.616231 -0.727768 -0.589 849 0
The algorithmic core of estimating M(X) is the computation of the vector h from the image data by a marching-cube algorithm and, therefore, the computation time is linear in the pixel number. The computation time on an Intel Xeon 5 148 at 2.33 GHz is about 5.2 ■ 10-9 s per pixel.
ERRORS OF ESTIMATION
We consider asymptotic errors of estimating the fiber length using the estimator Eq. 8 for vanishing
lattice distance, a ^ 0, and for a width of the fiber cross sections much smaller than the fiber length. An estimator 1(0) of the length 1(0) of the curve 0 can be deduced from Eq. 4,
1 m
1(0) = n E Ykdk E n (Lk + z)).
^ k=1 ze^Lj
First, we consider the particular case where 0 is a straight line of length l(0) = 1, direction d e S+ and having a random offset x0 uniformly distributed on the unit cell C,
0 = {sd : 0 < s < 1} + X0 . Then it follows that
\edk |
e E x{0 n (Lk + z))
ze ±L?
dk
where the expectation is taken with respect to the distribution of the random offset x0. From Eq. 4 we get the difference 8 of the estimated and the true lengths of 0,
1 m
8(o) = - 1, 6 G4 •
k=1
This error is due to the approximation in Eq. 6, which is asymptotically identical to the approximation in Eq. 4, since x(X n fk) = x(X n Lk) for sufficiently small a. Notice that ¿2 I66k|H2(d6) = n for all
6k € S+ which implies that fS2 8(6) H2(d6) = 0. This means that in the isotropic case where 6 is a random direction uniformly distributed on S22, the estimator 1(0) is unbiased for the fiber length £($). Numerical integration yields
var 1(0) = I 82(6) H2(d6) = 0.261875 ■ 10-3, 2n ,/s+
(as a I 0) and thus the standard deviation of the estimated length of 0 is about 1.6%. Furthermore, we have computed the minimum and maximum differences
min 8(d) = —7.34%
de S2
max8(d) = 2.28% .
deS2
The minimum is taken for directions parallel to the edges of the unit cell C, which might be disadvantageous for applications. For example, in many fiber-reinforced materials such as carbon-fiber-reinforced polymers (CRP) or glass-fiber-reinforced polymers (GRP), the fibers are nearly parallel to the
direction of the expected principal load, and for image acquisition this direction is usually chosen as one of the coordinate axes.
As a consequence of the above results, for each curve 0 of length 1(0), the error of the estimator l(0) is between -7.34% and 2.28 %,
0.9266l(0) < 1(0) < 1.02281(0) .
The large error is a consequence of the rough quadrature of the outer integral in Crofton's intersection formula (Eq. 1) induced by considering 2 x 2 x 2-pixel configurations fy in binary images.
Let now be a random system of non-overlapping curves, where the realizations are a. s. locally finite unions of curves and 5 is a parameter of its direction distribution. We assume that is macroscopically homogeneous, i.e., the distribution of is invariant with respect to translations. Then the specific curve length lV of is defined as the mean total length of curves in the unit cube, lV(&p) = el(n [0,1]3). We assume that Iv(0p) < The specific curve length can be estimated using lV = l(05 n W)/|W|, where W c r3 is a closed and convex window of volume |W | > 0.
The direction distribution of &p is the distribution of the curve direction d(s) at the typical point s of <Pp. Let
fß(ü, 9)
1
ß sin û
2n
;i+(ß2 -1) cos2 û3/2 '
ß > 0
be a class of probability density functions of the direction distribution, where the direction d = (sin $ cos y, sin $ sin y, cos $) is given in spherical polar coordinates, i.e., $ and y are the altitude and the longitude, respectively (Schladitz etal., 2006b; Ohseretal., 2009, sec. 7.5.2). It is easy to see that for ¡5 = 1 the fiber directions are uniformly distributed on S+ (isotropy), for 5 ^ 0 the curves tend to be straight and parallel to the z-axis, and for 5 ^ ~ the fiber directions are uniformly distributed on the positive unit half-circle in the xy-plane. Notice that 5 ^ 0 for glass fiber systems in GRP, and 5 » 1 for cellulose fiber systems in paper. The systematic error of estimating the specific curve length is E£V(&p) — lV(0p) and the relative systematic error of lV is 8 := ElV(&p)/£V— 1. Fig. 2 (red curve) shows 8 as a function of the anisotropy parameter ¡.
IMPROVED ESTIMATION
8 [%] 0- -0.31J	
	7l2^10 20 50 100 200 ß
-1	
-2_ -2.09	
-3-	
-4-	— 2 x 2 x 2-pixel configs. (13 directions)
-5-	— 3 x 3 x 3-pixel configs. (145 directions)
-6	3 x 3 x 3-pixel configs. (49 directions)
-7-7.33-1	
Fig. 2. The relative systematic error 8 of estimating the specific curve length lV (&ß) of a macroscopically homogeneous system &ß of curves with the density fß of the direction distribution.
In this section we will show that the estimation of the fiber length can considerably be improved by a finer discretization of the outer integral of Eq. 1. For this purpose we consider the 3 x 3 x 3- instead of the 2 x 2 x 2-pixel configurations of 3D images, i.e., the set {0, a}3 is now replaced with {0, a, 2a}3.
Fig. 3. A sketch for the shift of the triangles Tkii. One of the vertices of the unit cell Ck does not belong to the set {0, a, 2a}3 (left). Hence, Ck is tessellated into two triangles (bold). The upper triangle is shifted (right) in order to compute the Euler number locally.
SELECTION OF 145 SUBSPACES
We consider linear subspaces Lk e L2 for which a translation yk e {0, a, 2a}3 exists such that Lk + yk hits at least three affinely independent points of {0, a, 2a}3. It turns out that the number of subspaces Lk with
1
pairwise different normal directions dk increases from m = 13 to m = 145.
For all section lattices l|, the basis vectors uk and vk of the section lattice l| = l3 n Lk can be chosen such that there is a tessellation of the unit cell Ck of l| into triangles 7k, 1 and Tk,2 and there are translations yki G r3 with F0(Tk,i + yki) C {0, a, 2a}3, i = 1,2. Let fk be a 6-adjacency system of l2 generated by the tessellation of Ck into the triangles Tk1 and Tk 2 (Fig. 3). Then the integral of the mean curvature M(X) of the set X can be estimated using the right-hand side of Eq. 6 but with m = 145.
As a consequence of the discretization of the outer integral with respect to the section planes L1,..., L145, the asymptotic error for the estimator l(0) is reduced considerably. For a segment 0 of length 1, random direction d uniformly distributed on S+ and random offset uniformly distributed on [0,a]3 we get var l(0) = 0.021381 ■ 10-3, i.e., the standard deviation is about 0.46%. Furthermore, for a curve 0 with random offset it follows,
0.9872l(0) < 1(0) < 1.0119l(0) ,
and thus the error is now between -1.28% and 1.19%. Finally, we remark that in the case of parallel straight fibers, neither the minimal nor the maximum error is taken for a fiber direction equal to that of one of the coordinate axes. The relative systematic error 8 of estimating the specific curve length lV(&p) of a macroscopically homogeneous system of curves is shown in Fig. 2 (blue). Notice that the error vanishes for ¡5 ^ œ.
It should be noted that most translations yk;i are not uniquely determined. This means that the coefficients qi depend on specific choices for the yk i, but the estimator M(X) is independent of the yk i.
The step from the 2 x 2 x 2-pixel to the 3 x 3 x 3-pixel configurations involves a huge increase of the configurations' number from 28 to 227. This has consequences for the algorithmic implementation. It does not seem very reasonable to follow the algorithmic approach based on Eq. 8 which needs to allocate memory space for q and h. Alternatively, it is suggested to compute the vector g = (gk) with
2
gk = Z Z X0(X - X, Tk,i + yk,i)
xe L3 i=1
(9)
for k = 1,..., m. Using this, one obtains an estimate of the integral of the mean curvature from
Obviously, the computation of the vector g from 3D image data is more time consuming than that of h. This is the price we have to pay. Nevertheless, the algorithm based on Eq. 10 is linear in the pixel number, too. The runtime of the algorithm on an Intel Xeon 5 148 at 2.33 GHz is about 0.125 ■ 10—6 s per pixel.
Notice that the measure M(X) is concentrated on the surface of the set X. This means, that only those pixels y € L3 with 0 < #( (X - y) n {0,a,2a}3) < 27 contribute to M(X). When taking this into consideration, the runtime of the algorithm is linear in the number of surface pixels of the image. This can lead to further speeding-up of the algorithmic computing of M(X). For example, the typical volume fraction of the fibers in CRP and GRP is 20 %. If the fiber diameter is larger than 10 pixels, then at most 20 % of the pixels are surface pixels, i.e., the maximum computation time is 0.025 ■ 10-6 s per pixel.
SELECTION OF 49 SUBSPACES
Finally, we remark that the scattering of the directions dk on S+ is not very uniform, even when computing the outer integral of Eq. 1 based on the discrete set {Li,...,Li45} of section planes (i.e., induced by the 3 x 3 x 3-pixel configurations). Therefore, the question is as follows: Is it possible to select a subset of planes L1,..., L'm without a substantial loss of accuracy of the estimator 1(0)? For example, we chose all planes L'k from {L1;...,L145} with the normal directions dk = (pk — $k € S+ and (pk, $k € {0, a, 2a}3. There exist only m = 49 section planes Lk with pairwise different dk. The dk are scattered 'more uniformly' on S+ than the dk. Therefore, there is only a slight loss of accuracy of the estimator 1(0) obtained from a discretization of Eq. 1 with respect to L1 , .
L
, l49,
M(X) = Z Ykdkgk
(10)
k=1
0.97911(0) < 1(0) < 1.01351(0) ,
i.e., the error ranges from —2.09% to 1.35%. On the other hand, the runtime of an algorithm computing M(X) is reduced considerably (by about 1/3). The relative systematic error 8 of estimating the specific curve length lV(&p) of a macroscopically homogeneous system of curves is shown in Fig. 2 (green).
OPTIMAL CHOICE OF THE Yk
Until now, the weights yk for the quadrature rule for computing the outer integral of the Crofton formula (Eq. 1) are chosen more or less empirically as the
Voronoi areas. But are these weights really the best ones? As above, let 0 be a segment of length 1, a direction 6 uniformly distributed on S+ and with an offset x0 uniformly distributed in C. Then optimal weights yk can be defined as those values that minimize the variance
var£($) =
1 mm
2n3
2n3 LLW 2 66l|66k|H2(d6) - 1, (11) j=ik=1 JS+
subject to the normalization constraint. In other words, we are solving the system
for the unknown vector (y1,..., ym, A).
Numerical values of the optimal yk for the particular case m = 13 are yk = 0.584340 for the directions of the edges of C, Y = 0.474428 for the directions of the face diagonals, and yk = 0.420899 for the directions of the space diagonals. Clearly, the corresponding minimum variance var 1(0) = 0.260182 is smaller than that for the Voronoi areas. Nonetheless, the results presented in Table 2 are ambiguous: the estimation variances are reduced only slightly, while the lower bounds of the estimates decrease. Thus, we think that for 2 x 2 x 2-pixel configurations the Voronoi areas are close to optimum.
var £($) ^ min
m
LY = 2%
i=1
(12)
This can be done by applying the Lagrange multiplier method. Define
. m
A (y, A) = var £(0) + A[£y - 2n
i=i
then solving Eq. 12 is equivalent to solving
\iA(y, A) = 0,
that is
d var £($) + X = 0 dy1
dyn
var £($)+ X = 0
m
Ly = 2n
i=1
(13)
It holds
d	d mm	m
d- var £(0) = dY^L^L YjYkajk = 2 £ Ykm , Y Y j=1 k=1 k=1
where the coefficients
ajk = ¿/166j Neek|H2(de), j, k = 1,..., m
S+
are computed numerically. Using this, Eq. 13 is forming the linear equation system
( 2ö11
2am1 1
2a1m ^ /H\ ( 0 \

Ym \ X J
0
V 2n J
Table 2. Numerical values of the variance var I and the relative deviations min 5(6) and max5(6) for the optimal weights (bold face) as well as for the Voronoi areas (normal font).
pixel	m	var I($)	min 8(6)	max 8(6)
configs.		•10-3	[%]	[%]
2 x 2 x 2	13	0.260182	-7.74	2.24
		0.261 871	-7.34	2.28
3 x 3 x 3	145	0.006 189	-1.76	0.39
		0.021381	-1.28	1.19
3 x 3 x 3	49	0.015572	-3.26	0.73
		0.037 577	-2.09	1.35
For 3 x 3 x 3-pixel configurations one sees again a decrease in all three numbers (in both cases m = 145 and m = 49), which is now much more pronounced. This indicates that in these cases the Voronoi areas are far from the optimum. The reason for the decrease in the minimum in all three cases is not clear.
DISCUSSION
The approximation of x(X n Lk) with x(X n fk) used in the right-hand side of Eq. 6 depends on the lateral resolution of the image. But under which conditions for X is the approximation good enough? The following sufficient condition is given in Ohser et al. (2002): If the set X is morphologically open as well as morphologically closed with respect to all elements of the adjacency system fk, then X(X nLk) = x(X nfk), i = 1,...,m. Of course, this condition is very strong and never fulfilled in practical applications. For systems W of slightly curved fibers, however, it is sufficient to suppose that the fibers' thickness and spacing are large enough. More precisely, there must be a lower bound p > 0 such that K is morphologically open with respect to Bp, and the minimum distance between neighboring fibers is larger than p, with p = \f3a if the discretization of the
d
1
Crofton formula (Eq. 1) is induced by the 2 x 2 x 2-pixel configurations, and p = 2\/3a for 3 x 3 x 3. This means that the lateral resolution decreases for increasing angular resolution (increasing m), and vice versa. Obviously, it is necessary to find a balance between high lateral and high angular resolution. In cases of (approximately) isotropic fiber systems W, the systematic error of M(W) is independent of the number m of the section planes Lk and, hence, it is sufficient to estimate the fiber length based on the 2 x 2 x 2-pixel configurations. In the contrary case, for strongly anisotropic fiber systems it is necessary to scan the 3D image at a suitably high lateral resolution (small a), which allows to explore local knowledge from the 3 x 3 x 3-pixel configurations.
We recall that the basis vectors uk and vk of the section lattices lk are not uniquely defined by the rule formulated in Lemma 1. Their lengths ||uk || and ||vk|| determine the lateral resolution of estimating the Euler number x(X nLk). Thus, for further reducing the error resulting from the approximation of x(X n Lk) with X(X n fk), the basis vectors uk and vk of the section lattice lk should be chosen such that the spectral norm of the matrix (uk, vk) is minimal, where the minimum is taken over all possible choices of the lattice basis.
Let R be the direction distribution function of a macroscopically homogeneous system 0 of curves, and let W C r3 be a compact window with nonempty interior. The quantity
pV (L) = lim
e p(0 n pW, L) |pW |
L G L2 ,
(known as the rose of intersections of 0) is the cosine transform of R, that is pV(L) = lV y \uLd\R(d6), where uL is the normal vector of the section plane L. Vice versa, if R has a continuous differentiable density, this density is the inverse cosine transform of pV(L) (Mecke and Nagel, 1980). For a set of subspaces Li,...,Lm, the rose of intersections pV(L) can be estimated by pV (L) ^ dkgk/\W \, where in Eq. 9 the set X is replaced with a realization of W n W of the fiber system W observed through a window W. Discrete versions of the inverse cosine transform can be based on a series expansion of pV using spherical harmonics (Ohser and Schladitz, 2009, sec. 5.4.1). A numerically stable version of the inverse cosine transform is presented in Louis etal. (2011), which works well if m is large enough.
ACKNOWLEDGMENT
This work was supported by the German Federal Ministry of Education and Research (BMBF) under grant MNT/05 M10RCA/AMiNa.
REFERENCES
Baddeley A, Averback P (1983). Stereology of tubular structures. J Microsc 131:323-40.
Federer H (1959). Curvature measures. Trans Amer Math Soc 93:418-91.
Guderlei R, Klenk S, Mayer J, Schmidt V, Spodarev E (2007). Algorithms for the computation of Minkowski functionals of deterministic and random polyconvex sets. Image Vision Comput 25:464-74.
Hadwiger H (1957). Vorlesungen uber Inhalt, Oberflache und Isoperimetrie. Berlin: Springer-Verlag.
Kampf J (2012). A limitation of the estimation of intrinsic volumes via pixel configuration counts. WiMa Report, 144.
Kiderlen M, Rataj J (2006). On infinitesimal increase of volumes of morphological transforms. Mathematika 53:103-27.
Klenk S, Schmidt V, Spodarev E (2006). A new algorithmic approach to the computation of Minkowski functionals of germ-grain models. Comp Geom Th Appl 34:127-48.
Klette R, Kovalevsky VV, Yip B (1999). Length estimation of digital curves. In: Latecki LJ, Melter RA, Mount DM, Wu AY, eds. Vision Geometry VIII. Proc SPIE 3811:117-29.
Klette R, Rosenfeld A (2004). Digital geometry. Amsterdam: Morgan & Kaufman.
Lang C, Ohser J, Hilfer R (2001). On the analysis of spatial binary images. J Microsc 203:303-13.
Lindblad J (2005). Surface area estimation of digitized 3D objects using weighted local configurations. Image Vision Comput 23:111-22.
Lindblad J, Nystrom I (2002). Surface area estimation of digitized 3D objects using local computations. In: Proc 10th Int Conf Discrete Geom Computer Imagery. Lect Not Comput Sci 2301:267-78.
Louis AK, Riplinger M, Spiess M, Spodarev E (2011). Inversion algorithms for the spherical Radon and cosine transform. Inverse Problems 27:035015.
Mecke J, Nagel W (1980). Stationäre räumliche Faserprozesse und ihre Schnittzahlrosen. Elektron Informationsverarb Kyb 16:475-83..
Meschenmoser D, Spodarev E (2012). On the computation of intrinsic volumes. Preprint.
Mrkvicka T, Rataj J (2008). On the estimation of intrinsic volume densities of stationary random closed sets. Stochastic Proc Appl 118:213-31.
Mrkvicka T, Rataj J (2009). On the estimation of intrinsic volume densities of stationary random closed sets via parallel sets in the plane. Kybernetika 45:931-45.
Mullikin J, Verbeek P (1993). Surface estimation of digital planes. Bioimaging 1:6-16.
Nagel W, Ohser J, Pischang K (2000). An integral-geometric approach for the Euler-Poincare characteristic of spatial images. J Microsc 198:54-62.
Ohser J, Mucklich F (2000). Statistical analysis of microstructures in materials science. Chichester: J Wiley & Sons.
Ohser J, Nagel W, Schladitz K (2002). The Euler number of discretized sets - on the choice of adjacency in homogeneous lattices. In: Mecke KR, Stoyan D, eds. Morphology of condensed matter. Lect Not Phys 600:275-98.
Ohser J, Nagel W, Schladitz K (2003). The Euler number of discretised sets - surprising results in three dimensions. Image Anal Stereol 22:11-9.
Ohser J, Nagel W, Schladitz K (2009). Miles formulae for Boolean models observed on lattices. Image Anal Stereol 28:77-92.
Ohser J, Schladitz K (2009). 3D images of materials structures - Processing and analysis. Weinheim, Berlin: Wiley VCH.
Santalo LA (1976). Integral geometry and geometric probability. Reading Mass.: Addison-Wesley.
Schladitz K, Ohser J, Nagel W (2006a). Measurement of intrinsic volumes of sets observed on lattices. In: Kuba A, Nyul LG, Palagyi K, eds. Proc 13th Int Conf Discrete Geom Comput Imagery. Lect Not Comput Sci
4245:247-58.
Schladitz K, Peters S, Reinel-BitzerD, Wiegmann A, Ohser J (2006b). Design of aoustic trim based on geometric modeling and flow simulation for non-woven. Comp Materials Sci 38:56-66.
Schmidt V, Spodarev E (2005). Joint estimators for the specific intrinsic volumes of stationary random sets. Stoch Proc Appl 115:959-81.
Schneider R (1993). Convex bodies: The Brunn-Minkowski theory. Encyclopedia of mathematics and its application Vol. 44. Cambridge: Cambridge University Press.
Svane AM (2012). Local digital estimators of intrinsic volumes for Boolean models and in the design based setting :preprint.
Vogel HJ, Weller U, Schluter S (2010). Quantification of soil structure based on Minkowski functions. Comput Geosci 36:1236-46.
Weitkamp T, Haas D, Wegrzynek D, Rack A (2011). Ankaphase: software for single-distance phase-retrieval from inline X-ray phase contrast radiographs. J Synchrotron Radiat 18:617-29.
Windreich G, Kiryati N, Lohmann G (2003). Surface area estimation in practice. In: Nystrom I, di Baja GS, Svensson S, eds. Proc 11th Int Conf Discrete Geom Comput Imagery. Lect Not Comput Sci 2886:358-67.
Ziegel J, Kiderlen M (2010). Estimation of the surface area and surface area measure of three-dimensional sets from digitizations. Image Vision Comput 28:64-77.