IMPROVED ESTIMATION OF FIBER LENGTH FROM 3-DIMENSIONAL IMAGES

A new method is presented for estimating the specific fiber len gth from 3D images of macroscopically homogeneous fiber systems. The method is based on a discrete v e sion of the Crofton formula, where local knowledge from 3× 3× 3-pixel configurations of the image data is exploited. It is s hown that the relative error resulting from the discretization of the outer integr al of the Crofton formula amounts at most 1.2 %. An algorithmic implementation of the method is simple and the r untime as well as the amount of memory space are low. The estimation is significantly improved by conside ring 3× 3× 3-pixel configurations instead of 2×2×2, as already studied in literature.


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 of materials properties.Currently, the development in fiber system characterization is mainly driven by computed microtomography (µCT) which has been established as an appropriate source of threedimensional (3D) image data of sufficiently high lateral resolution and contrast between fibers and the matrix constituents or air in pore spaces, Weitkamp et al. (2011).In most cases, 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 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.
Except 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 √ 3a 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, there are certain similarities of estimating length and surface area.
First, we refer on Lindblad and Nyström (2002), who suggested to estimate the surface area of a body as the sum of the areas of the surface patches obtained by ordinary surface rendering, see Section 3.6.2 in Ohser and Schladitz (2009).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, see e. g.Hadwiger (1957), p. 215 andSantaló (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 et al. (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.For each 2 × 2 × 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 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 Windreich et al. (2003).A further approach is presented in Chapter 4 of Ohser and Mücklich (2000) and in Lang et al. (2001) where the weights are computed using a discrete version of one of Crofton's intersection formulas, see also Schladitz et al. (2006a) and Chapter 5 in Ohser and Schladitz (2009).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 estimator based on optimal weights completely ignore the image data, see 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, see Schmidt and Spodarev (2005), Klenk et al. (2006), Mrkvička and Rataj (2008), Mrkvička and Rataj (2009), and 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 Boolean models in 2D in Guderlei et al. (2007) and Mrkvička and Rataj (2008) show, that the accuracy is sometimes higher than of that of the estimators given in Chapter 4 of Ohser and Mücklich (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 yield 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 on cubic primitive lattices, see Ohser et al. (2009); Ohser and Schladitz (2009).This is an important fact since many 3D image acquisition techniques produce images on non-cubic lattices.Furthermore, methods based on discrete versions of the Crofton formula yield reasonable results even for low lateral resolution, where the pixel size is in the range of the fiber radius.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, see Section 5.4 in Ohser and Schladitz (2009).

FIBERS AND FIBER LENGTH
In our setting a fiber Φ in R 3 with a convex crosssection of nonempty interior can be defined as follows: Let f : R → R 3 be an arclength parametrization of a finite immersed curve φ ⊂ R 3 with f (0) = f (ℓ), where f is twice continuously differentiable and ℓ is the curve length, The first derivative of f is the normal tangent vector of the curve φ at the point f (s).More precisely, let S 2 Let B r denote the ball of radius r and ⊕ be the Minkowski addition.The curve φ may be shaped in such a way that the parallel set φ ⊕ B r of φ is morphologically regular, i. e. there is an ε > 0 such that φ ⊕ B r is morphologically closed as well as morphologically open with respect to B ε .A necessary condition for the morphological regularity of φ ⊕ B r is that φ is smooth enough, more precisely, the curvature κ(s) = f ′′ (s) of φ must be smaller than the inverse radius, κ(s) < 1/r for all s ∈ [0, ℓ].
Let K ⊆ B r be a compact, convex and morphologically regular set.Then a fiber Φ is defined as the set Φ = φ ⊕ K.As a consequence of the morphological regularity of K, the mean curvature H 1 exists for each point on the surface ∂ Φ of Φ and, using the 2-dimensional Hausdorff measure H 2 , the integral of the mean curvature is defined as From Theorem 5.9 in Federer (1959) it follows that M(Φ) → πℓ(φ ) as r ↓ 0. In other words, M(Φ)/π is an appropriate estimate of the fiber length ℓ as r ≪ ℓ.
We consider a piecewise straight approximation φ m of the curve φ , Let now L 2 be the set of 2-dimensional linear subspaces in R 3 (i.e. the set of planes hitting the origin).By ⊥ L, ν⊥ L and µ we denote the orthogonal space of a plane L ∈ L 2 , the Lebesgue measure on ⊥ L, and the rotation invariant probability measure on the unit sphere S 2 of R 3 , respectively.Then, using Crofton's intersection formula, the integral of the mean curvature M(X ) of a polyconvex set X can be written in the form where χ X ∩ (L +x) is the Euler number of the planar section X ∩ (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 formula (1), where the plane L is replaced with a 2dimensional section lattice, the translations L + x over the orthogonal space ⊥ L are replaced with a translation of the section lattice over the corresponding translation lattice, and the integrals are numerically computed by quadrature rules.
Let L 3 = aZ 3 be a homogeneous lattice in R 3 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 ∈ L 2 .Assume now that there is a translation y ∈ {0, a} 3 such that the intersection of the translated plane hits at least three vertices of C, Then L 2 = L ∩ L 3 forms a homogeneous 2dimensional section lattice of L 3 .This means that there are linearly independent vectors u, v ∈ R 3 with u + y, v + y ∈ {0, a} 3 generating the lattice L 2 , Without loss of generality we assume that the normal direction θ see Ohser and Schladitz (2009).The lattice ⊥ L 2 is called the translation lattice of L 2 .The above results are summarized in the following lemma: Lemma 1 Let be given a lattice L 3 = aZ 3 , a plane L ∈ L 2 and a vector y ∈ {0, a} 3 such that # (L + y) ∩ {0, a} 3 ≥ 3. Then there are three linearly independent vectors u, v, w ∈ {−a, 0, a} 3 with L 3 = {iu + jv + kw : i, j, k ∈ Z} and L 3 ∩ L = {iu + jv : i, j ∈ Z}.
As has been pointed out Ohser and Schladitz (2009), p. 157, there exist m = 13 section lattices L 2 k , k = 1, . . ., m, with pairwise different normal directions θ k of the corresponding section planes k .Now, using a simple rectangular quadrature rule, the inner integral p(X , L) of the Crofton formula (1) can be approximated with where the d k serve as the weights of the quadrature rule.
Consider the Voronoi tessellation of the unit sphere S 2 with respect to the point field {θ 1 , . . . ,θ m , −θ 1 , . . . ,−θ m }.By γ k we denote the area of the Voronoi cell with respect to the direction θ k , k = 1, . . ., m. Numerical calculation yields γ k = 0.575 263 for the three directions parallel to the edges of the unit cell C, γ k = 0.464 711 for the directions of the face diagonals of C, and γ k = 0.442 280 for the directions of the space diagonals of C. When using a 2-dimensional analog of the rectangular quadrature rule, Eq. ( 1) yields and summarizing ( 2) and (3), one obtains an approximation of the integral of the mean curvature, In order to give an estimator of the Euler number χ X ∩ (L k + z) we consider a 6-adjacency system F k on the 2-dimensional section lattice L 2 k .The 6adjacency system is generated by a tessellation of the unit cell C k into two closed triangles T k,1 and T k,2 with T k,1 ∪ T k,2 = C k .Formally, this can be done as follows: Let conv ζ denote the convex hull of a subset ζ ⊆ F 0 (C k ), where F 0 (C k ) denotes the set of vertices of C k .Then the local adjacency system with respect to the tessellation mentioned above consists of the unit cell C k , the triangles T k,1 and T k,2 , their edges, their vertices, and the empty set.Finally, the adjacency F k is the union of F 0 k + x over all lattice translations, see Ohser and Schladitz (2009), p. 50.Now we introduce a discretization X ⊓ F k of the planar section X ∩ spanL k with respect to the adjacency system F k , k carry the same information on X .The additivity of the Euler number ensures that it can be computed locally.Using the Euler-Poincaré formula we get where F j (F k ) denotes the set of the j-dimensional polytopes belonging to F k , and 1 is the indicator function.
The Euler number χ(X ⊓ F k ) can be seen as an approximation of χ(X ∩ L k ) which yields the estimator (5) The right-hand side is a discretization of the Crofton formula induced by the sampling of X on L 3 .

ALGORITHMIC IMPLEMENTATION
Our aim is now to present a factorization of the estimator (5) in the form where the vector q depends on the kind of discretization of the Crofton formula and h is the vector of 2 × 2 × 2-pixel configurations in a binary image with the foreground pixels X ∩ L 3 .
Let T be a triangle in R 3 , we define the local Euler number χ 0 (X , T ) as As in the previous section, we consider an 6-adjacency system F k on L 2 k generated from a tessellation of the unit cell C k of L 2 k into two triangles T k,1 and T k,2 .Then the system T = {T 1 , T 2 , . ..} of all triangles belonging to F k forms a face-to-face tesselation of L k , i. e. for i = j the intersection T i ∩ T j is either an edge of T i , a vertex of T i or empty.Notice that there are two possible tessellations of C k .Depending on the index k, we choose one of them arbitrarily.
Using the Euler-Poincaré formula and the inclusion-exclusion principle for additive functionals one can prove the following lemma: From the inclusion-exclusion principle it follows that the Euler-Poincaré formula can be rewritten as Furthermore, by ξ 0 , . . . ,ξ n−1 ⊆ {0, a} 3 we denote the n = 255 local configurations on L 3 , where ξ c l = {0, a} 3 \ ξ l is the complementary configuration of ξ l .The numbers h l of occurrences of the ξ l in the sampling X ∩ L 3 are and using the convention 0 • ∞ = 0, one obtains Finally, we get as an estimator of M(X ) with the weights It is important to see that d k /a does not depent on a.
Eq. ( 6) shows that the integral of the mean curvature M(X ) of X can be estimated from the sampling X ∩ L 3 as a scalar product qh of the vectors q = (q l ) and h = (h l ) where q is independent of the data as well as on the lattice distance a. Numerical values for the q l given in Table 1, where the pictograms stand for the representatives of the 22 equivalence classes of the 2 × 2 × 2-pixel configurations ξ l (full disc for foreground, empty disc for background).Notice that this weights are equivalent to those given in Ohser and Schladitz (2009) 1.The weights q l for the estimation of the integral of the mean curvature M of a set X sampled on a homogeneous lattices L 3 = aZ 3 .
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 (6) for vanishing lattice distance, a → 0, and for a width of the fiber cross-sections much smaller than the fiber length, r → 0. An estimator l(φ ) of the length ℓ(φ ) of the underlying curve φ can be deduced from (4), First, we consider the particular case where φ is a straight line of length ℓ(φ ) = 1, direction θ ∈ S 2 + and having a random offset x 0 uniformly distributed on the unit cell C, Then it follows that and from (4) we get the difference δ of the estimated and the true lengths of φ , This means that in the isotropic case where θ is a random direction uniformly distributed on S 2 + , the estimator l(φ ) is unbiased for the fiber length ℓ(φ ).Numerical integration yields var l(φ ) = 1 4π and thus the standard deviation of the estimated length of φ is about 1.6 %.Furthermore, we have computed the minimum and maximum differences min 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-fiberreinforced 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.
The large error is a consequence of the rough quadrature of the outer integral in Crofton's intersection formula (1) induced by considering 2 × 2 × 2-pixel configurations ξ l in binary images.
Let now φ β be a random system of nonoverlapping curves where the realizations are a.s. locally finite unions of curves, where β 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 ℓ V of φ β is defined as the mean total length of curves in the unit cube, The direction distribution of φ β is the distribution of the curve direction θ (s) at the typical point s of φ .Let be a class of probability density functions of the direction distribution, where the direction θ = (sin ϑ cos ϕ, sin ϑ sin ϕ, cos ϑ ) with 0 ≤ ϑ ≤ π and 0 ≤ ϕ < 2π is given in spherical polar coordinates, i. e. ϑ and ϕ are the altitude and the longitude, respectively, see Schladitz et al. (2006b) and Section 7.5.2 in Ohser et al. (2009).It is easy to see that for β = 1 the fiber directions are uniformly distributed on S 2 + (isotropy of φ ), for β → 0 the curves tend to be straight and parallel to the z-axis, and for β → ∞ the fiber directions are uniformly distributed on the positive unit half-circle in the xy-plane.Notice that β 0 for glass fiber systems in GRP, and β ≫ 1 for cellulose fiber systems in paper.Fig. 1 (red) shows the relative systematic error δ = E lV (φ β )/ℓ V (φ β ) − 1 of the estimator depending on the anisotropy parameter β .(13 directions) −0.31 −2.09 Fig. 1.The relative systematic error δ of estimating the specific curve length ℓ V (φ β ) of a macroscopically homogeneous system φ β of curves with the density f β of the direction distribution.

IMPROVED ESTIMATION
The estimation of the fiber length can considerably be improved by a finer discretization of the outer integral of (1).For this purpose we consider the 3 × 3 × 3-instead of the 2 × 2 × 2-pixel configurations of 3D imges, i. e. the set {0, a} 3 is now replaced with {0, a, 2a} 3 .Fig. 2. A sketch for the shift of the triangles T k,i .One of the vertices of the unit cell C k does not belong to the set {0, a, 2a} 3 (left).Hence, C k 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 L k ∈ L 2 for which a translation y k ∈ {0, a, 2a} 3 exists such that L k + y k hits at least three affinely independent points of {0, a, 2a} 3 .It turns out that the number of subspaces L k with pairwise different normal directions θ k increases from m = 13 to m = 145.For all section lattices L 2 k , the basis vectors u k and v k of the section lattice L 2 k = L 3 ∩ L k can be chosen such that there is a tessellation of the unit cell C k of L 2 k into triangles T k,1 and T k,2 and there are translations y k,i ∈ R 3 with F 0 (T k,i + y k,i ) ∈ {0, a, 2a} 3 , i = 1, 2. Let F k be the 6-adjacency system of L 2 k generated by the tessellation of C k into the triangles T k,1 and T k,2 .Then the integral of the mean curvature M(X ) of the set X can be estimated using the right-hand side of (5) but with m = 145.
As a consequence of the discretization of the outer integral with respect to the section planes L 1 , . . ., L 145 , the asymptotic error for the estimator l(φ ) is reduced considerably.For a segment φ of length 1, random direction θ uniformly distributed on S 2 + and random offset uniformly distributed on [0, a] 3 we get var l(Φ) = 0, 021 381 •10 −3 , i. e. the standard deviation is about 0.46 %.Furthermore, for a curve φ with random offset it follows, 0.987 2 ℓ(φ ) ≤ l(φ ) ≤ 1.011 9 ℓ(φ ), 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 δ of estimating the specific curve length ℓ V (φ β ) of a macroscopically homogeneous system φ β of curves is shown in Fig. 1 (blue).Notice that the error vanishes for β → ∞.
It should be noted that most translations y k,i are not uniquely determined.This means that the coefficients q l depend on specific choices for the y k,i , but the estimator M(X ) is independent of the y k,i .
The step from the 2 × 2 × 2-to the 3 × 3 × 3-pixel configurations involves a huge increase of the configurations' number from 2 8 to 2 27 .This has consequences for the algorithmic implementation.It does not seem very reasonable to follow the algorithmic approach based on (6) which needs to allocate memory space for q h.Alternatively, it is suggested to compute the vector g = (g k ) with 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 ( 8) 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 ∈ L 3 with 0 < # (X − y) ∩ {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 θ k on S 2 + is not very uniform, even when computing the outer integral of (1) based on the discrete set {L 1 , . . ., L 145 } of section planes (i.e. induced by the 3 × 3 × 3-pixel configurations).Therefore, the question is as follows: Is it possible to select a subset of planes L ′ 1 , . . ., L ′ m without a substantial loss of accuracy of the estimator l(φ )?For example, we chose all planes L ′ k from {L 1 , . . ., L 145 } with the normal directions θ ′ k = ϕ k − ψ k ∈ S 2 + and ϕ k , ψ k ∈ {0, a, 2a} 3 .There exist only m = 49 section planes L ′ k with pairwise different θ ′ k .The θ ′ k are scattered 'more uniformly' on S 2 + than the θ k .Therefore, there is only a slight loss of accuracy of the estimator l(Φ) obtained from a discretization of (1) with respect to L ′ 1 , . . ., L ′ 49 , 0.979 1 ℓ(φ ) ≤ l(φ ) ≤ 1.013 5 ℓ(φ ), 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 δ of estimating the specific curve length ℓ V (φ β ) of a macroscopically homogeneous system φ β of curves is shown in Fig. 1 (green).

OPTIMAL CHOICE OF THE γ k
Until now, the weights γ k for the quadrature rule for computing the outer integral of the Crofton formula (1) are chosen more or less empirically as the Voronoi areas.But are these weights really the best ones?As above, let φ be a segment of length 1, a direction θ uniformly distributed on S 2 + and with an offset x 0 uniformly distributed in C. Then optimal weights γ i can be defined as those values that minimize the variance var l subject to the normalization constraint.In other words, we are solving the system var l(φ This can be done by applying Lagrange multiplier method.Define where the coefficients are computed numerically.Using this, Eq. ( 11) is forming the linear equation system for the unknown vector (γ 1 , . . ., γ m , λ ).
Numerical values of the optimal γ k for the particular case m = 13 are γ k = 0.584 340 for the directions of the edges of C, γ k = 0.474 428 for the directions of the face diagonals, and γ k = 0.420 899 for the directions of the space diagonals.Clearly, the corresponding minimum variance var l(φ ) = 0.260 182 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 × 2 × 2-pixel configurations the Voronoi areas are close to optimum.pixel m var l(φ ) min δ (θ ) max δ (θ ) configs.
•10 For 3 × 3 × 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 ∩ L k ) with χ(X ⊓ F k ) used in the right-hand side of (5) depends on the lateral resolution of the image.But in what 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 F k , then χ(X ∩ L k ) = χ(X ⊓ F k ), i = 1, . . ., m.Of course, this condition is very strong and never fulfilled in practical applications.For systems Φ 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 ρ > 0 such that K ⊖ B ρ = / 0 and the minimum distance between neighboring fibers is larger than ρ, with ρ = √ 3a if the discretization of the Crofton formula (1) is induced by the 2 × 2 × 2pixel configurations, and ρ = 2 √ 3a for 3 × 3 × 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 Φ, the systematic error of l(Φ) is independent of the number m of the section planes L k and, hence, it is sufficient to estimate the fiber length based on the 2 × 2 × 2pixel 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 × 3 × 3-pixel configurations.
We recall that the basis vectors u k and v k of the section lattices L k are not uniquely defined by the rule formulated in Lemma 1. Their lengths u k and v k determine the lateral resolution of estimating the Euler number χ(X ∩ L k ).Thus, for further reducing the error resulting from the approximation of χ(X ∩ L k ) with χ(X ⊓ F k ), the basis vectors u k and v k of the section lattice L k should be chosen such that the spectral norm of the matrix (u k , v k ) is minimal, where the minimum is taken over all possible choices of the lattice basis.
The number χ(φ ∩ L) of intersections carries information about the anisotropy of a system φ of curves.More precisely, the density of the direction distribution (i.e. the rose of directions) of φ is the inverse cosine transform of For a set of subspaces L 1 , . . ., L k the quantity p(φ , L) can be estimated from the corresponding fiber system Φ by p(φ , L k ) ≈ d k g k , where in (7) the set X is replaced with Φ. Discrete versions of the inverse cosine transform can be based on a series expansion of p(φ , L) using spherical harmonics, see Section 5.4.1 in Ohser and Schladitz (2009).For a numerically stable version of the inverse cosine transform see e. g. Louis et al. (2011).Using this, the p(φ , L k ) can be applied to estimate the rose of directions, in particular for m = 145.

Table 2 .
Numerical values of the variance var l(φ ) and the relative deviations min δ (θ ) and max δ (θ ) for the optimal weights (bold face) as well as for the Voronoi areas (normal font).