Voronoi-based estimation of Minkowski tensors from finite point samples
aa r X i v : . [ m a t h . M G ] J a n Voronoi-based estimation of Minkowski tensorsfrom finite point samples
Daniel Hug, Markus Kiderlen, and Anne Marie SvaneMarch 28, 2018
Abstract
Intrinsic volumes and Minkowski tensors have been used to describethe geometry of real world objects. This paper presents an estimatorthat allows to approximate these quantities from digital images. It isbased on a generalized Steiner formula for Minkowski tensors of setsof positive reach. When the resolution goes to infinity, the estimatorconverges to the true value if the underlying object is a set of positivereach. The underlying algorithm is based on a simple expression interms of the cells of a Voronoi decomposition associated with the image.
Intrinsic volumes, such as volume, surface area, and Euler characteristic,are widely-used tools to capture geometric features of an object; see, forinstance, [1, 24, 22]. Minkowski tensors are tensor valued generalizationsof the intrinsic volumes, associating with every sufficiently regular compactset in R d a symmetric tensor, rather than a scalar. They carry informationabout geometric features of the set such as position, orientation, and eccen-tricity. For instance, the volume tensor – defined formally in Section 2 – ofrank 0 is just the volume of the set, while the volume tensors of rank 1 and2 are closely related to the center of gravity and the tensor of inertia, re-spectively. For this reason, Minkowski tensors are used as shape descriptorsin materials science [29, 31], physics [14], and biology [3, 36].The main purpose of this paper is to present estimators that approximateall the Minkowski tensors of a set K when only weak information on K isavailable. More precisely, we assume that a finite set K which is close to K in the Hausdorff metric is known. The estimators are based on the Voronoidecomposition of R d associated with the finite set K , following an idea ofM´erigot et al. [21]. What makes these estimators so interesting is that theyare consistent; that is, they converge to the respective Minkowski tensorsof K when applied to a sequence of finite approximations converging to K in the Hausdorff metric. We emphasize that the notion of ‘estimator’ is1sed here in the sense of digital geometry [17] meaning ‘approximation ofthe true value based on discrete input’ and should not be confused with thestatistical concept related to the inference from data with random noise.The main application we have in mind is the case where K is a digitizationof K . This is detailed in the following.As data is often only available in digital form, there is a need for es-timators that allow us to approximate the Minkowski tensors from digitalimages. In a black-and-white image of a compact geometric object K ⊆ R d ,each pixel (or voxel) is colored black if its midpoint belongs to K and whiteotherwise. Thus, the information about K contained in the image is the setof black pixel (voxel) midpoints K = K ∩ a L , where L is the lattice formedby all pixel (voxel) midpoints and a − is the resolution. A natural criterionfor the reliability of a digital estimator is that it yields the correct tensorwhen a → + . If this property holds for all objects in a given family of sets,for instance, for all sets with smooth boundary, then the estimator is called multigrid convergent for this class.Digital estimators for the scalar Minkowski tensors, that is, for theintrinsic volumes, are widespread in the digital geometry literature; see,e.g., [17, 24, 25] and the references therein. For Minkowski tensors up torank two, estimators based on binary images are given in [28] for the two-dimensional and in [30] for the three-dimensional case. Even for the class ofconvex sets, multigrid convergence has not been proven for any of the abovementioned estimators. The only exception are volume related quantities.Most of the above mentioned estimators are n -local for some given fixed n ∈ N . We call an estimator n -local if it depends on the image only throughthe histogram of all n × · · · × n configurations of black and white points. Forinstance, a natural surface area estimator [19] in three-dimensional spacescans the image with a voxel cube of size 2 × × n -local estimators is that they are intuitive, easy to implement, and thecomputation time is linear in the number of pixels or voxels.However, many n -local estimators are not multigrid convergent for con-vex sets; see [32] and the detailed discussion in Section 6. This implies thatmany established estimators, like the mentioned one in [19] cannot be multi-grid convergent for convex sets. All the estimators of 2D-Minkowski tensorsin [28] are 2-local. By the results in [32], the estimators for the perimeterand the Euler characteristic can thus not be multigrid convergent for convexsets. The multigrid convergence of the other estimators has not been inves-tigated. The algorithms for 3D-Minkowski tensors in [30] have as input atriangulation of the object’s boundary, and the way this triangulation is ob-tained determines whether the resulting estimators are n -local or not. Thereare no known results on multigrid convergence for these estimators either.Summarizing, to the best of our knowledge, this paper presents for the first2ime estimators of all Minkowski tensors of arbitrary rank that come witha multigrid convergence proof for a class of sets that is considerably largerthan the class of convex sets.The present work is inspired by [21], and we therefore start by recallingsome basic notions from this paper. For a nonempty compact set K , theauthors of [21] define a tensor valued measure, which they call the Voronoicovariance measure , defined on a Borel set A ⊆ R d by V R ( K ; A ) = Z K R A ( p K ( y ))( y − p K ( y ))( y − p K ( y )) ⊤ dy. Here, K R is the set of points at distance at most R > K and p K is the metric projection on K : the point p K ( x ) is the point in K closest to x , provided that this closest point is unique. The metric projection of K iswell-defined on R d with the possible exception of a set of Lebesgue-measurezero; see, e.g., [7].The paper [21] uses the Voronoi covariance measure to determine localfeatures of surfaces. It is proved there that if K ⊆ R is a smooth surface,then V R ( K ; B ( x, r )) ≈ π R r (cid:18) u ( x ) u ( x ) ⊤ + r X i =1 , k i ( x ) P i ( x ) P i ( x ) ⊤ (cid:19) , (1)where B ( x, r ) is the Euclidean ball with midpoint x ∈ K and radius r , u ( x ) is one of the two surface unit normals at x ∈ K , P ( x ) , P ( x ) are theprincipal directions and k ( x ) , k ( x ) the corresponding principal curvatures.Hence, the eigenvalues and -directions of the Voronoi covariance measurecarry information about local curvatures and normal directions.Assuming that a compact set K approximates K , [21] suggests to es-timate V R ( K ; · ) by V R ( K ; · ). It is shown in that paper that V R ( K ; · )converges to V R ( K ; · ) in the bounded Lipschitz metric when K → K in theHausdorff metric. Moreover, if K is a finite set, then the Voronoi covariancemeasure can be expressed in the form V R ( K ; A ) = X x ∈ K ∩ A Z B ( x,R ) ∩ V x ( K ) ( y − x )( y − x ) ⊤ dy. Here, V x ( K ) is the Voronoi cell of x in the Voronoi decomposition of R d associated with K . Thus, the estimator which is used to approximate V R ( K ; A ) is easily computed. Given the Voronoi cells of K , each Voronoicell contributes with a simple integral. Figure 1 (a) shows the Voronoi cellsof a finite set of points on an ellipse. The Voronoi cells are elongated inthe normal direction. This is the intuitive reason why they can be used toapproximate (1).The Voronoi covariance measure V R ( K ; A ) can be identified with a sym-metric 2-tensor. In the present work, we explore how natural extensions of3 a ) ( b )Figure 1: (a). The Voronoi cells of a finite set of points on a surface. (b).A digital image and the associated Voronoi cells.the Voronoi covariance measure can be used to estimate general Minkowskitensors. The generalizations of the Voronoi covariance measure, which wewill introduce, will be called Voronoi tensor measures . We will then showhow the Minkowski tensors can be recovered from these. When we applythe results to digital images, we will work with full-dimensional sets K , andthe finite point sample K is obtained from the representation K = K ∩ a L of a digital image of K . The Voronoi cells associated with K = K ∩ a L are sketched in Figure 1 (b). Taking point samples from K with increasingresolution, convergence results will follow from an easy generalization of theconvergence proof in [21].The paper is structured as follows: In Section 2, we recall the definitionof Minkowski tensors and the classical as well as a local Steiner formula forsets of positive reach. In Section 3, we define the Voronoi tensor measures,discuss how they can be estimated from finite point samples, and explainhow the Steiner formula can be used to connect the Voronoi tensor measureswith the Minkowski tensors. Section 4 is concerned with the convergenceof the estimator. The results are specialized to digital images in Section 5.Finally, the estimator is compared with existing approaches in Section 6. We work in Euclidean space R d with scalar product h· , ·i and norm | · | . TheEuclidean ball with center x ∈ R d and radius r ≥ B ( x, r ),and we write S d − for the unit sphere in R d . Let ∂A and int A be theboundary and the interior of a set A ⊆ R d , respectively. The k -dimensionalHausdorff-measure in R d is denoted by H k , 0 ≤ k ≤ d . Let C d be the familyof nonempty compact subsets of R d and K d ⊆ C d the subset of nonemtpycompact convex sets. For two compact sets K, M ∈ C d , we define their Hausdorff distance by d H ( K, M ) = inf { ε > | K ⊆ M ε , M ⊆ K ε } . T p denote the space of symmetric p -tensors (tensors of rank p ) over R d . Identifying R d with its dual (via the scalar product), a symmetric p -tensor defines a symmetric multilinear map ( R d ) p → R . Letting e , . . . , e d be the standard basis in R d , a tensor T ∈ T p is determined by its coordinates T i ...i p = T ( e i , . . . , e i p )with respect to the standard basis, for all choices of i , . . . , i p ∈ { , . . . , d } .We use the norm on T p given by | T | = sup (cid:8) | T ( v , . . . , v p ) | | | v | = · · · = | v p | = 1 (cid:9) for T ∈ T p . The same definition is used for arbitrary tensors of rank p .The symmetric tensor product of y , . . . , y m ∈ R d is given by the sym-metrization y ⊙ · · · ⊙ y m = ( m !) − P ⊗ mi =1 y σ ( i ) , where the sum extends overall permutations σ of { , . . . , m } and ⊗ is the usual tensor product. We write x r for the r -fold tensor product of x ∈ R d . For two symmetric tensors of theform T = y ⊙ · · · ⊙ y r and T = y r +1 ⊙ · · · ⊙ y r + s , where y , . . . , y r + s ∈ R d ,the symmetric tensor product T ⊙ T of T and T , which we often ab-breviate by T T , is the symmetric tensor product of y , . . . , y r + s . This isextended to general symmetric tensors T and T by linearity. Moreover, itfollows from the preceding definitions that | y ⊙ · · · ⊙ y m | ≤ | y | · · · | y m | ,y , . . . , y m ∈ R d .For any compact set K ⊆ R d , we can define an element of T r called the r th volume tensor Φ r, d ( K ) = 1 r ! Z K x r dx. For s ≥ r,sd ( K ) = 0. Some of the volume tensors have well-known physical interpretations. For instance, Φ , d ( K ) is the usual volumeof K , Φ , d ( K ) is up to normalization the center of gravity, and Φ , d ( K ) isclosely related to the tensor of inertia. All three tensors together can be usedto find the best approximating ellipsoid of a particle [36]. The sequence ofall volume tensors (Φ r, d ( K )) ∞ r =0 determines the compact set K uniquely. Forconvex sets in the plane even the following stability result [10, Remark 4.4.]holds: If K, L ∈ K are contained in the unit square and have coincidingvolume tensors up to rank r , then their distance, measured in the symmetricdifference metric H (cid:0) ( K \ L ) ∪ ( L \ K ) (cid:1) , is of order O ( r − / ) as r → ∞ .We will now define Minkowski surface tensors . These can also be used tocharacterize the shape of an object or the structure of a material as in [3, 14].They require stronger regularity assumptions on K . Usually, like in [26,Section 5.4.2], the set K is assumed to be convex. However, as Minkowskitensors are tensor-valued integrals with respect to the generalized curvature5easures (also called support measures) of K , they can be defined wheneverthe latter are available. We will use this to define Minkowski tensors for setsof positive reach.First, we recall the definition of a set of positive reach and explain howcurvature measures of such sets are determined (see [8, 35]). For a compactset K ∈ C d , we let d K ( x ) denote the distance from x ∈ R d to K . Then,for R ≥ K R = { x ∈ R d | d K ( x ) ≤ R } is the R -parallel set of K . The reach Reach( K ) of K is defined as the supremum over all R ≥ x ∈ R d with d K ( x ) < R there is a unique closest point p K ( x ) in K . We say that K has positive reach if Reach( K ) >
0. Smooth surfaces(of class C , ) are examples of sets of positive reach, and compact convexsets are characterized by having infinite reach. By definition, the map p K is defined everywhere on K R if R <
Reach( K ). Let K ⊆ R d be a (compact)set of positive reach. The (global) Steiner formula for sets with positivereach states that for all R <
Reach( K ) the R -parallel volume of K is apolynomial, that is, H d ( K R ) = d X k =0 κ d − k R d − k Φ , k ( K ) . (2)Here κ j is the volume of the unit ball in R j and the numbers Φ , ( K ) , . . . , Φ , d ( K ) are the so-called intrinsic volumes of K . They are special cases ofthe Minkowski tensors to be defined below. Some of them have well-knowninterpretations. As mentioned, Φ , d ( K ) is the volume of K . Moreover,2Φ , d − ( K ) is the surface area, Φ , d − ( K ) is proportional to the total meancurvature, and Φ , ( K ) is the Euler characteristic of K . For convex sets, (2)is the classical Steiner formula which holds for all R ≥ generalized curvature measures Λ k ( K ; · ) of K , for k = 0 , . . . , d − R d × S d − . They are determined by the following local Steiner formulawhich holds for all
R <
Reach( K ) and all Borel set B ⊆ Σ: H d (cid:16)n x ∈ K R \ K | (cid:16) p K ( x ) , x − p K ( x ) | x − p K ( x ) | (cid:17) ∈ B o(cid:17) = d − X k =0 R d − k κ d − k Λ k ( K ; B ) . (3)The coefficients Λ k ( K ; B ) on the right side of (3) are signed Borel measuresΛ k ( K ; · ) evaluated on B ⊆ Σ. These measures are called the generalizedcurvature measures of K . Since the pairs of points in B on the left sideof (3) always consist of a boundary point of K and an outer unit normalof K at that point, each of the measures Λ k ( K, · ) is concentrated on theset of all such pairs. For this reason, the generalized curvature measures6 k ( K ; · ), k ∈ { , . . . , d − } , are also called support measures . They describethe local boundary behavior of the part of ∂K that consists of points x with an outer unit normal u such that ( x, u ) ∈ B . A description of thegeneralized curvature measures Λ k ( K, · ) by means of generalized curvaturesliving on the normal bundle of K was first given in [35] (see also [26, § k ( K, Σ)are the intrinsic volumes.Based on the generalized curvature measures, for every k ∈ { , . . . , d − } , r, s ≥ K ⊆ R d with positive reach, we define the Minkowskitensor Φ r,sk ( K ) = 1 r ! s ! ω d − k ω d − k + s Z Σ x r u s Λ k ( K ; d ( x, u ))in T r + s . Here ω k is the surface area of the unit sphere S k − in R k . Moreinformation on Minkowski tensors can for instance be found in [13, 20,27, 15]. As in the case of volume tensors, the Minkowski tensors carrystrong information on the underlying set. For instance, already the sequence(Φ ,s ( K )) ∞ s =0 determines any K ∈ K d up to a translation. A stability re-sult also holds: if K and L are both contained in a fixed ball and havethe same tensors Φ ,s of rank s ≤ s , then a translation of K is close to L in the Hausdorff metric and the distance is O ( s − β ) as s → ∞ for any0 < β < / ( n + 1); see [18, Theorem 4.9].One can define local Minkowski tensors in a similar way (see [11]). Fora Borel set B ⊆ Σ, for k ∈ { , . . . , d − } , r, s ≥ K ⊆ R d withpositive reach, we putΦ r,sk ( K ; B ) = 1 r ! s ! ω d − k ω d − k + s Z B x r u s Λ k ( K ; d ( x, u ))and, for a Borel set A ⊆ R d ,Φ r, d ( K ; A ) = 1 r ! Z K ∩ A x r dx. In order to avoid a distinction of cases, we also write Φ r, d ( K ; A × S d − )instead of Φ r, d ( K ; A ). Moreover, we define Φ r,sd ( K ; · ) = 0 if s ≥
1. Thelocal Minkowski tensors can be used to describe local boundary properties.For instance, local 1- and 2-tensors are used for the detection of sharp edgesand corners on surfaces in [6]. They also carry information about normaldirections and principal curvatures as explained in the introduction.We conclude this section with a general remark on continuity proper-ties of the Minkowski tensors. Although the functions K Φ r,sk ( K ) arecontinuous when considered in the metric space ( K d , d H ), they are not con-tinuous on C d . (For instance, the volume tensors of a finite set are alwaysvanishing, but finite sets can be used to approximate any compact set inthe Hausdorff metric.) This is the reason why our approach requires an7pproximation argument with parallel sets as outlined below. The consis-tency of our estimator is mainly based on a continuity result for the metricprojection map. We quote this result [4, Theorem 3.2] in a slightly differentformulation which is symmetric in the two bodies involved. Let k f k L ( E ) bethe usual L -norm of the restriction of f to a Borel set E ⊆ R d . Proposition 2.1.
Let ρ > and let E ⊆ R d be a bounded measurable set.Then there is a constant C = C ( d, diam( E ∪ { } ) , ρ ) > such that k p K − p K k L ( E ) ≤ C d H ( K, K ) for all K, K ∈ C d with K, K ⊆ B (0 , ρ ) .Proof. Let E ′ be the convex hull of E and observe that k p K − p K k L ( E ) ≤ k p K − p K k L ( E ′ ) . It is shown in [4, Lemma 3.3] (see also [8, Theorem 4.8]) that the map v K : R d → R given by v K ( x ) = | x | − d K ( x ) is convex and that its gradientcoincides almost everywhere with 2 p K . Since E ′ has rectifiable boundary,[4, Theorem 3.5] implies that k p K − p K k L ( E ′ ) ≤ c ( d )( H d ( E ′ ) + ( c + k d K − d K k ∞ ,E ′ ) H d − ( ∂E ′ )) × k d K − d K k ∞ ,E ′ . Here c = diam(2 p K ( E ′ ) ∪ p K ( E ′ )) ≤ K ∪ K ) ≤ ρ and thesupremum-norm k · k ∞ ,E ′ on E ′ can be estimated by k d K − d K k ∞ ,E ′ ≤ E ′ ∪ K ∪ K ) k d K − d K k ∞ ,E ′ ≤ (cid:2) diam( E ′ ∪ { } ) + 2 ρ (cid:3) d H ( K, K ) . Moreover, intrinsic volumes are increasing on the class of convex sets, so H d ( E ′ ) ≤ H d ( B (0 , diam( E ′ ∪ { } ))) H d − ( ∂E ′ ) ≤ H d − ( ∂B (0 , diam( E ′ ∪ { } ))) . Together with the trivial estimate d H ( K, K ) ≤ ρ and with the equalitydiam( E ∪ { } ) = diam( E ′ ∪ { } ), this yields the claim.The authors of [4] argue that the exponent 1 / In Section 3.1 below, we define the Voronoi tensor measures and show howthe Minkowski tensors can be obtained from these. We then explain in Sec-tion 3.2 how the Voronoi tensor measures can be estimated from finite pointsamples. As a special case, we obtain estimators for all intrinsic volumes.This is detailed in Section 3.3. 8 .1 The Voronoi tensor measures
Let K be a compact set. Here and in the following subsections, we let r, s ∈ N and R ≥
0. Define the T r + s -valued measures V r,sR ( K ; · ) given on aBorel set A ⊆ R d by V r,sR ( K ; A ) = Z K R A ( p K ( x )) p K ( x ) r ( x − p K ( x )) s dx. (4)When K is a smooth surface, V , R ( K ; · ) corresponds to the Voronoi covari-ance measure in [21]. We will refer to the measures defined in (4) as the Voronoi tensor measures . Note that if f : R d → R is a bounded Borelfunction, then Z R d f ( x ) V r,sR ( K ; dx ) = Z K R f ( p K ( x )) p K ( x ) r ( x − p K ( x )) s dx ∈ T r + s . (5)Suppose now that K has positive reach with Reach( K ) > R . Then aspecial case of the generalized Steiner formula derived in [12] (or an extensionof (3)) implies the following version of the local Steiner formula for theVoronoi tensor measures: V r,sR ( K ; A ) = d X k =1 ω k Z Σ Z R A ( x ) t s + k − x r u s dt Λ d − k ( K ; d ( x, u ))+ { s =0 } Z K ∩ A x r dx = r ! s ! d X k =0 κ k + s R s + k Φ r,sd − k ( K ; A × S d − ) , (6)where A ⊆ R d is a Borel set. In particular, the total measure is V r,sR ( K ) = V r,sR ( K ; R d ) = r ! s ! d X k =0 κ k + s R s + k Φ r,sd − k ( K ) . Note that the special case r = s = 0 is the Steiner formula (2) for sets withpositive reach.Equation (6), used for different parallel distances R , can be solved for theMinkowski tensors. More precisely, choosing d + 1 different values 0 < R <. . . < R d < Reach( K ) for R , we obtain a system of d + 1 linear equations: V r,sR ( K ; A )... V r,sR d ( K ; A ) = r ! s ! κ s R s . . . κ s + d R s + d ... ... κ s R sd . . . κ s + d R s + dd Φ r,sd ( K ; A × S d − )...Φ r,s ( K ; A × S d − ) . (7)9ince the Vandermonde-type matrix A r,sR ,...,R d = r ! s ! κ s R s . . . κ s + d R s + d ... ... κ s R sd . . . κ s + d R s + dd ∈ R ( d +1) × ( d +1) (8)in (7) is invertible, the system can be solved for the tensors, and thus weget Φ r,sd ( K ; A × S d − )...Φ r,s ( K ; A × S d − ) = (cid:16) A r,sR ,...,R d (cid:17) − V r,sR ( K ; A )... V r,sR d ( K ; A ) . (9)If s >
0, then Φ r,sd ( K ; A × S d − ) = 0 by definition, so we may omit one ofthe equations in the system (7). Let K be a compact set of positive reach. Suppose that we are given a com-pact set K that is close to K in the Hausdorff metric. In the applicationswe have in mind, K is a finite subset of K , but this is not necessary for thealgorithm to work. Based on K , we want to estimate the local Minkowskitensors of K . We do this by approximating V r,sR k ( K ; A ) in Formula (9) by V r,sR k ( K ; A ), for k = 0 , . . . , d and A ⊆ R d a Borel set. This leads to thefollowing set of estimators for Φ r,sk ( K ; A × S d − ), k ∈ { , . . . , d } : ˆΦ r,sd ( K ; A × S d − )...ˆΦ r,s ( K ; A × S d − ) = (cid:16) A r,sR ,...,R d (cid:17) − V r,sR ( K ; A )... V r,sR d ( K ; A ) (10)with A r,sR ,...,R d given by (8). Setting A = R d in (10), we obtain estimatorsˆΦ r,sk ( K ) = ˆΦ r,sk ( K ; R d × S d − )of the intrinsic volumes. Note that this approach requires an estimate forthe reach of K because we need to choose 0 < R < · · · < R d < Reach( K ).The idea to invert the Steiner formula is not new. It was used in [4] toapproximate curvature measures of sets of positive reach. In [16] and [23] itwas used to estimate intrinsic volumes but without proving convergence forthe resulting estimator.We now consider the case where K is finite. Let V x ( K ) = { y ∈ R d | p K ( y ) = x } denote the Voronoi cell of x ∈ K with respect to the set K . Since R d isthe union of the finitely many Voronoi cells of K , it follows that K R is the10igure 2: The Voronoi decomposition (blue lines) and R -parallel set (redcurve) associated with a digital image.union of the R -bounded parts B ( x, R ) ∩ V x ( K ), x ∈ K , of the Voronoi cells V x ( K ), x ∈ K , which have pairwise disjoint interiors. Thus (4) simplifiesto V r,sR ( K ; A ) = X x ∈ K ∩ A x r Z B ( x,R ) ∩ V x ( K ) ( y − x ) s dy. (11)Like the Voronoi covariance measure, the Voronoi tensor measure V r,sR ( K ; A )is a sum of simple contributions from the individual Voronoi cells.An example of a Voronoi decomposition associated with a digital imageis sketched in Figure 2. The original set K is the disk bounded by the innerblack circle, and the disk bounded by the outer black circle is its R -parallelset K R . The finite point sample is K = K ∩ Z , which is shown as the setof red dots in the picture, and the red curve is the boundary of its R -parallelset. The Voronoi cells of K are indicated by blue lines. The R -boundedpart of one of the Voronoi cells is the part that is cut off by the red arc. Recall that Φ , k ( K ) = Λ k ( K ; R d ) is the k th intrinsic volume. Thus, Section3.2 provides estimators for all intrinsic volumes as a special case. This caseis particularly simple. The measure V , R ( K ; A ) is simply the volume of alocal parallel set V , R ( K ; A ) = H d (cid:0) { x ∈ K R | p K ( x ) ∈ A } (cid:1) , V , R ( K ) = H d ( K R ) .
11n particular, if K ⊆ R d is a compact set with Reach( K ) > R , then Equa-tion (6) reduces to the usual local Steiner formula H d ( { x ∈ K R | p K ( x ) ∈ A } ) = d X k =0 κ k R k Λ d − k ( K ; A × S d − ) , and to the (global) Steiner formula (2) if A = R d .In this case, our algorithm approximates the parallel volume H d ( K R )by H d ( K R ). In the example in Figure 2, this corresponds to approximatingthe volume of the larger black disk by the volume of the region bounded bythe red curve. This volume is again the sum of the volumes of the regionsbounded by the red and blue curves. In other words, it is the sum of volumesof the R -bounded Voronoi cells on the right-hand side of the equation V , R ( K ; A ) = X x ∈ K ∩ A H d ( B ( x, R ) ∩ V x ( K )) . In Section 3.2 we have only considered estimators for local tensors of theform Φ r,sk ( K ; A × S d − ), where K ⊆ R d is a set with positive reach. Thenatural way to estimate Φ r,sk ( K ; B ), for a measurable set B ⊆ Σ, would beto copy the idea in Section 3.2 with V r,sR ( K ; A ) replaced by the followinggeneralization of the Voronoi tensor measures, W r,sR ( K ; B ) = Z K R \ K B ( p K ( x ) , u K ( x )) p K ( x ) r ( x − p K ( x )) s dx, (12)where u K ( x ) = ( x − p K ( x )) / | x − p K ( x ) | estimates the normal direction. Ofcourse, this definition works for any K ∈ C d . Moreover, we could defineestimators related to (12) whenever we have a set K which approximates K . However, even if K has positive reach, the map x u K ( x ) is notLipschitz on K R \ K , and therefore the convergence results in Section 4 willnot work with this definition. Since the map x u K ( x ) is Lipschitz on K R \ K R/ , it is natural to proceed as follows. For any K ∈ C d , we define V r,sR ( K ; B ) = Z K R \ K R/ B ( p K ( x ) , u K ( x )) p K ( x ) r ( x − p K ( x )) s dx. (13)Note that V r,sR ( K ; · ) = W r,sR ( K ; · ) − W r,sR/ ( K ; · ) , (14)where W r,sR ( K ; · ) is defined as in (12). We will not use the notation W r,sR ( K ; · )in the following. If K has positive reach and 0 < R < reach( K ), then thegeneralized Steiner formula yields V r,sR ( K ; B ) = r ! s ! d X k =1 κ s + k R s + k (1 − − ( s + k ) )Φ r,sd − k ( K ; B ) . < R < . . . < R d < reach( K ), we can recover theMinkowski tensors from Φ r,sd − ( K ; B )...Φ r,s ( K ; B ) = (cid:0) A r,sR ,...,R d (cid:1) − V r,sR ( K ; B )... V r,sR d ( K ; B ) where A r,sR ,...,R d = 1 r ! s ! κ s +1 (1 − − ( s +1) ) R s +11 . . . κ s + d (1 − − ( s + d ) ) R s + d ... ... κ s +1 (1 − − ( s +1) ) R s +1 d . . . κ s + d (1 − − ( s + d ) ) R s + dd is a regular matrix. Using this, we can define estimators for Φ r,sk ( K ; B ), for0 ≤ k ≤ d −
1, by Φ r,sd − ( K ; B )...Φ r,s ( K ; B ) = (cid:0) A r,sR ,...,R d (cid:1) − V r,sR ( K ; B )... V r,sR d ( K ; B ) , where K is a compact set which approximates K . Convergence of thesemodified estimators will be discussed in Section 4.The estimators Φ r,sk can be used to approximate local tensors of the formΦ r,sk ( K ; B ) where the set B ⊆ Σ involves normal directions. Thus, they aremore general than ˆΦ r,sk . However, (14) shows that estimating V r,sR ( K ; B )requires an approximation of two parallel sets, rather than one. We thereforeexpect more severe numerical errors for Φ r,sk . In this section we prove the main convergence results. This is an immediategeneralization of [21, Theorem 5.1].
For a bounded Lipschitz function f : R d → R , we let | f | ∞ denote the usualsupremum norm, | f | L = sup (cid:26) | f ( x ) − f ( y ) || x − y | | x = y (cid:27) the Lipschitz semi-norm, and | f | bL = | f | L + | f | ∞ d bL be the bounded Lipschitz metric onthe space of bounded T p -valued Borel measures on R d . For any two suchmeasures µ and ν on R d , the distance with respect to d bL is defined by d bL ( µ, ν ) = sup (cid:26)(cid:12)(cid:12)(cid:12)(cid:12) Z f dµ − Z f dν (cid:12)(cid:12)(cid:12)(cid:12) | | f | bL ≤ (cid:27) , where the supremum extends over all bounded Lipschitz functions f : R d → R with | f | bL ≤
1. The following theorem shows that the map K
7→ V r,sR ( K ; · )is H¨older continuous with exponent with respect to the Hausdorff metric on C d (restricted to compact subsets of a fixed ball) and the bounded Lipschitzmetric. In the proof, we use the symmetric difference A ∆ B = ( A \ B ) ∪ ( B \ A )of sets A, B ⊆ R d . Theorem 4.1.
Let
R, ρ > and r, s ∈ N be given. Then there is a positiveconstant C = C ( d, R, ρ, r, s ) such that d bL ( V r,sR ( K ; · ) , V r,sR ( K ; · )) ≤ C d H ( K, K ) for all compact sets K, K ⊆ B (0 , ρ ) .Proof. Let f with | f | bL ≤ (cid:12)(cid:12)(cid:12)(cid:12) Z R d f ( x ) V r,sR ( K ; dx ) − Z R d f ( x ) V r,sR ( K ; dx ) (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) Z K R f ( p K ( x )) p K ( x ) r ( x − p K ( x )) s dx − Z K R f ( p K ( x )) p K ( x ) r ( x − p K ( x )) s dx (cid:12)(cid:12)(cid:12)(cid:12) ≤ I + II, (15)where I is the integral Z K R ∩ K R | f ( p K ( x )) p K ( x ) r ( x − p K ( x )) s − f ( p K ( x )) p K ( x ) r ( x − p K ( x )) s | dx and II = ρ r R s H d ( K R ∆ K R ) . By [4, Corollary 4.4], there is a constant c = c ( d, R, ρ ) > H d ( K R ∆ K R ) ≤ c d H ( K, K ) (16)14hen d H ( K, K ) ≤ R/
2. Replacing c by a possibly even bigger constant,we can ensure that (16) also holds when R/ ≤ d H ( K, K ) ≤ ρ . Hence, II ≤ c d H ( K, K ) (17)with some constant c = c ( d, R, ρ, r, s ) > (cid:12)(cid:12)(cid:12)(cid:12) m K i =1 y i − m K i =1 z i (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12)(cid:12) m O i =1 y i − m O i =1 z i (cid:12)(cid:12)(cid:12)(cid:12) ≤ m X j =1 | y j − z j | j − Y i =1 | y i | m Y i = j +1 | z i | , (18)with m = r + s and the rank-one tensors y = . . . = y r = p K ( x ) , y r +1 = . . . = y r + s = x − p K ( x ) ,z = . . . = z r = p K ( x ) , z r +1 = . . . = z r + s = x − p K ( x ) , we get | f ( p K ( x )) p K ( x ) r ( x − p K ( x )) s − f ( p K ( x )) p K ( x ) r ( x − p K ( x )) s |≤ | f ( p K ( x )) − f ( p K ( x )) || p K ( x ) | r | x − p K ( x ) | s + | f ( p K ( x )) | r X j =1 | p K ( x ) − p K ( x ) || p K ( x ) | j − | p K ( x ) | r − j | x − p K ( x ) | s ++ | f ( p K ( x )) | s X j =1 | p K ( x ) − p K ( x ) || p K ( x ) | r | x − p K ( x ) | j − | x − p K ( x ) | s − j . Since we assumed that | f | bL ≤
1, we get I ≤ ( r + s + 1) max { ρ, } r max { R, } s Z K R ∩ K R | p K ( x ) − p K ( x ) | dx ≤ c d H ( K, K ) . (19)The existence of the constant c = c ( d, R, ρ, r, s ) in the last inequality isguaranteed by Proposition 2.1 with K R ∩ K R as the set E , because thischoice of E satisfies diam( E ∪ { } ) ≤ ρ + R ).When r = s = 0 and f = 1, the above proof simplifies to Inequality (16)as I vanishes. Hence we obtain the following strengthening of the theorem,which is relevant for the estimation of intrinsic volumes. Theorem 4.2.
Let
R, ρ > . Then there is a constant C = C ( d, R, ρ ) > such that (cid:12)(cid:12)(cid:12) V , R ( K ) − V , R ( K ) (cid:12)(cid:12)(cid:12) ≤ C d H ( K, K ) for all compact sets K, K ⊆ B (0 , ρ ) . Theorem 4.3.
Let r, s ∈ N and R > . If K i → K with respect to theHausdorff metric on C d , as i → ∞ , then V r,sR ( K i ; A ) → V r,sR ( K ; A ) in thetensor norm, for every Borel set A ⊆ R d which satisfies H d ( p − K ( ∂A ) ∩ K R ) = 0 . (20) Proof.
Convergence of tensors is equivalent to coordinate-wise convergence.Hence, it is enough to show that the coordinates satisfy V r,sR ( K i ; A ) i ...i r + s → V r,sR ( K ; A ) i ...i r + s as i → ∞ , for all choices of indices i . . . i r + s ; see the notation at the beginning ofSection 2.We write T K ( x ) = p K ( x ) r ( x − p K ( x )) s . Then V r,sR ( K ; A ) i ...i r + s = Z K R A ( p K ( x )) T K ( x ) i ...i r + s dx is a signed measure. Let T K ( x ) + i ...i r + s and T K ( x ) − i ...i r + s denote the positiveand negative part of T K ( x ) i ...i r + s , respectively. Then V r,sR ( K ; A ) ± i ...i r + s = Z K R A ( p K ( x )) T K ( x ) ± i ...i r + s dx are non-negative measures such that V r,sR ( K ; · ) i ...i r + s = V r,sR ( K ; · ) + i ...i r + s − V r,sR ( K ; · ) − i ...i r + s . The proof of Theorem 4.1 can immediately be generalized to show that V r,sR ( K i ; · ) ± i ...i r + s converges to V r,sR ( K ; · ) ± i ...i r + s in the bounded Lipschitz norm(as i → ∞ ), and hence the measures converge weakly. In particular, theyconverge on every continuity set of V r,sR ( K ; · ) ± i ...i r + s . If H d ( p − K ( ∂A ) ∩ K R ) =0, then A is such a continuity set. Remark 4.4.
Though relatively mild, the condition H d ( p − K ( ∂A ) ∩ K R ) = 0can be hard to control if K is unknown. It is satisfied if, for instance, K and A are smooth and their boundaries intersect transversely. A specialcase of this is when K is a smooth surface and A is a small ball centered onthe boundary of K . This is the case in the application from [21] that wasdescribed in the introduction. Examples where it is not satisfied are when A = K or when K is a polytope intersecting ∂A at a vertex. Remark 4.5.
Let f : R d → R be a bounded measurable function. Wedefine V r,sR ( K ; f ) := Z R d f ( x ) V r,sR ( K ; dx ) . V r,sR ( K ; A ) = V r,sR ( K ; A ) for every Borel set A ⊆ R d . Then, Theorem4.3 is equivalent to saying that, for all continuous test functions f : R d → R , V r,sR ( K i ; f ) → V r,sR ( K ; f ) , as i → ∞ , in the tensor norm, whenever K i → K with respect to the Hausdorff metricon C d , as i → ∞ . Thus, if one is interested in the local behaviour of Φ r,sk ( K ; · )at a neighborhood A , like in [21], then one can studyΦ r,sk ( K ; f ) := Z Σ f ( x ) x r u s Λ k ( K ; d ( x, u )) , where f is a continuous function with support in A . This avoids the extracondition (20).As the matrix A r,sR ,...,R d in the definition (10) of ˆΦ r,sk ( K ; A × S d − ) doesnot depend on the set K , the above results immediately yield a consistencyresult for the estimation of the Minkowski tensors. We formulate this onlyfor A = R d . Corollary 4.6.
Let ρ > and K be a compact subset of B (0 , ρ ) of positivereach such that Reach( K ) > R d > . . . > R > . Let K ⊆ B (0 , ρ ) be acompact set. Then there is a constant C = C ( d, R , . . . , R d , ρ ) such that (cid:12)(cid:12)(cid:12) ˆΦ , k ( K ) − Φ , k ( K ) (cid:12)(cid:12)(cid:12) ≤ C d H ( K , K ) , for all k ∈ { , . . . , d } .For r, s ∈ N there is a constant C = C ( d, R , . . . , R d , ρ, r, s ) such that (cid:12)(cid:12)(cid:12) ˆΦ r,sk ( K ) − Φ r,sk ( K ) (cid:12)(cid:12)(cid:12) ≤ C d H ( K , K ) , for all k ∈ { , . . . , d − } . Finally, we state the convergence results for the modified estimators forΦ r,sk ( K ; B ), where B ⊆ Σ is a Borel set, that were defined in Section 3.4. Themap x x/ | x | is Lipschitz on R d \ int( B (0 , R/ /R , and therefore the mapping u K , which was defined after (12), satisfies | u K ( x ) − u K ( x ) | ≤ R | p K ( x ) − p K ( x ) | , for x ∈ ( K R \ K R/ ) ∩ ( K R \ K R/ ). Moreover, (cid:16) K R \ K R/ (cid:17) ∆ (cid:16) K R \ K R/ (cid:17) ⊆ (cid:0) K R ∆ K R (cid:1) ∪ (cid:16) K R/ ∆ K R/ (cid:17) . Using this, it is straightforward to generalize the proofs of Theorems 4.1 and4.3 to obtain the following result. 17 heorem 4.7.
Let
R, ρ > and r, s ∈ N be given. Then there is a positiveconstant C = C ( d, R, ρ, r, s ) such that d bL ( V r,sR ( K ; · ) , V r,sR ( K ; · )) ≤ C d H ( K, K ) for all compact sets K, K ⊆ B (0 , ρ ) . This in turn leads to the next convergence result.
Theorem 4.8.
Let r, s ∈ N and R > . If K, K i ∈ C d are compact setssuch that K i → K in the Hausdorff metric, as i → ∞ , then V r,sR ( K i ; B ) converges to V r,sR ( K ; B ) in the tensor norm, for any measurable set B ⊆ Σ satisfying H d ( { x ∈ K R | ( p K ( x ) , u K ( x )) ∈ ∂B } ) = 0 . Here ∂B is the boundary of B as a subset of Σ .If B satisfies this condition and Reach ( K ) > R d , then lim i → Φ r,sk ( K i ; B ) = Φ r,sk ( K ; B ) . Remark 4.9.
We can argue as in Remark 4.5 to see that if
K, K i ∈ C d arecompact sets such that K i → K in the Hausdorff metric, as i → ∞ , then V r,sR ( K i ; g ) → V r,sR ( K ; g ) , as i → ∞ , whenever g : Σ → R is a continuous test function and V r,sR ( K ; g ) is definedsimilarly as before.If K satisfies Reach( K ) > R d , we get Φ r,sk ( K i ; g ) → Φ r,sk ( K ; g ), as i → ∞ . Our main motivation for this paper is the estimation of Minkowski tensorsfrom digital images. Recall that we model a black-and-white digital imageof K ⊆ R d as the set K ∩ a L , where L ⊆ R d is a fixed lattice and a >
0. Werefer to [2] for basic information about lattices.The lower dimensional parts of K are generally invisible in the digitalimage. When dealing with digital images, we will therefore always assumethat the underlying set is topologically regular, which means that it is theclosure of its own interior.In digital stereology, the underlying object K is often assumed to belongto one of the following two set classes: • K is called δ -regular if it is topologically regular and the reach of itsclosed complement cl( R d \ K ) and the reach of K itself are both atleast δ >
0. This is a kind of smoothness condition on the boundary,ensuring in particular that ∂K is a C manifold (see the discussionafter Definition 1 in [34]). 18 K is called polyconvex if it is a finite union of compact convex sets.While convex sets have infinite reach, note that polyconvex sets dogenerally not have positive reach. Also note that for a compact convexset K ⊆ R d , the set cl( R d \ K ) need not have positive reach.It should be observed that for a compact set K ⊆ R d both assumptionsimply that the boundary of K is a ( d − ∂K is the image of a bounded subset of R d − under a Lipschitz map),which is a much weaker property that will be sufficient for the analysis inSection 5.1. Simple and efficient estimators for the volume tensors Φ r, d ( K ) of a (topo-logically regular) compact set K are already known and are usually basedon the approximation of K by the union of all pixels (voxels) with midpointin K . This leads to the estimator φ r, d ( K ∩ a L ) = 1 r ! X z ∈ K ∩ a L Z z + aV ( L ) x r dx, where V ( L ) is the Voronoi cell of 0 in the Voronoi decomposition generatedby L . This, in turn, can be approximated byˆ φ r, d ( K ∩ a L ) = a d r ! H d ( V ( L )) X z ∈ K ∩ a L z r . When r ∈ { , } , we even have φ r, d ( K ∩ a L ) = ˆ φ r, d ( K ∩ a L ).Choose C > V ( L ) ⊆ B (0 , C ). Then K ∆ [ z ∈ K ∩ a L ( z + aV ( L )) ⊆ ( ∂K ) aC . In fact, if x ∈ (cid:2)S z ∈ K ∩ a L ( z + aV ( L )) (cid:3) \ K , then there is some z ∈ K ∩ a L such that x ∈ z + aV ( L ) and x / ∈ K . Since z ∈ K and x / ∈ K , we have[ x, z ] ∩ ∂K = ∅ . Moreover, x − z ∈ aV ( L ) ⊆ B (0 , aC ), and hence | x − z | ≤ aC . This shows that x ∈ ( ∂K ) aC . Now assume that x ∈ K and x / ∈ ( ∂K ) aC .Then B ( x, ρ ) ⊆ K for some ρ > aC . Since S z ∈ a L ( z + aV ( L )) = R d , thereis some z ∈ a L such that x ∈ z + aV ( L ). Hence x − z ∈ aV ( L ) ⊆ B (0 , aC ).We conclude that z ∈ B ( x, aC ) ⊆ K , therefore z ∈ K ∩ a L and thus x ∈ S z ∈ K ∩ a L ( z + aV ( L )).Hence | φ r, d ( K ∩ a L ) − Φ r, d ( K ) | ≤ r ! Z ( ∂K ) aC | x | r dx. (21)19f H d ( ∂K ) = 0, then the integral on the right-hand side goes to zero bymonotone convergence, solim a → + φ r, d ( K ∩ a L ) = Φ r, d ( K ) . (22)If ∂K is ( d − ∂K is the image of a bounded subset of R d − under a Lipschitz map, then H d ( ∂K ) = 0. Since ∂K is compact, [9, Theorem 3.2.39] implies thatlim a → + H d (( ∂K ) aC ) /a exists and equals a fixed multiple of H d − ( ∂K ) whichis finite. Hence, (21) shows that the speed of convergence in (22) is O ( a ) as a → + .Inequality (18) yields that | x r − z r | ≤ aCr ( | x | + aC ) r − whenever x ∈ z + aV ( L ) and r ≥
1. Therefore, | ˆ φ r, d ( K ∩ a L ) − φ r, d ( K ∩ a L ) | ≤ aC ( r − X z ∈ K ∩ a L Z z + aV ( L ) ( | x | + aC ) r − dx ≤ aC ( r − Z K aC ( | x | + aC ) r − dx, which shows that lim a → + ˆ φ r, d ( K ∩ a L ) = Φ r, d ( K ) , provided that H d ( ∂K ) = 0. If ∂K is ( d − O ( a ).Hence, we suggest to simply use the estimators ˆ φ r, d ( K ∩ a L ) for thevolume tensors. This estimator can be computed much faster and moredirectly than ˆΦ r, d ( K ∩ a L ). Moreover, it does not require an estimate forthe reach of K , and it converges for a much larger class of sets than thoseof positive reach. For the estimation of the remaining tensors we suggest to use the Voronoitensor measures. Choosing K = K ∩ a L in (11), we obtain V r,sR ( K ∩ a L ; A ) = X x ∈ K ∩ a L ∩ A x r Z B ( x,R ) ∩ V x ( K ∩ a L ) ( y − x ) s dy, (23)where A ⊆ R d is a Borel set.To show some convergence results in Corollary 5.2 below, we first notethat the digital image converges to the original set in the Hausdorff metric. Lemma 5.1. If K is compact and topologically regular, then lim a → + d H ( K, K ∩ a L ) = 0 . If K is δ -regular, then d H ( K, K ∩ a L ) is of order O ( a ) . The same holds if K is topologically regular and polyconvex. roof. Recall from [2, p. 311] that µ ( L ) = max x ∈ R d dist( x, L ) is well definedand denotes the covering radius of L .Let ε > K is compact, there are points x , . . . , x m ∈ K such that K ⊆ m [ i =1 B ( x i , ε ) . Using the fact that K is topologically regular, we conclude that there arepoints y i ∈ int( K ) ∩ int( B ( x i , ε )) for i = 1 , . . . , m . Hence, there are ε i ∈ (0 , ε ) such that B ( y i , ε i ) ⊆ K ∩ B ( x i , ε ) for i = 1 , . . . , m . Let0 < a < min { ε i /µ ( L ) | i = 1 , . . . , m } . Since ε i /a > µ ( L ) it followsthat a L ∩ B ( y i , ε i ) = ∅ , for i = 1 , . . . , m . Thus we can choose z i ∈ a L ∩ B ( y i , ε i ) ⊆ a L ∩ K for i = 1 , . . . , m . By the triangle inequality, wehave | z i − x i | ≤ ε i + 2 ε ≤ ε , and hence x i ∈ ( K ∩ a L ) + B (0 , ε ), for i = 1 , . . . , m . Therefore, K ⊆ ( K ∩ a L ) + B (0 , ε ) if a > K is δ -regular, for some δ >
0. We choose 0 < a <δ/ (2 µ ( L )). Since aµ ( L ) < δ/
2, for any x ∈ K there is a ball B ( y, aµ ( L )) ofradius aµ ( L ) such that x ∈ B ( y, aµ ( L )) ⊆ K . From a L ∩ B ( y, aµ ( L )) = ∅ we conclude that there is a point z ∈ K ∩ a L with | x − z | ≤ aµ ( L ). Hence x ∈ ( K ∩ a L ) + B (0 , aµ ( L )), and therefore d H ( K, K ∩ a L ) ≤ aµ ( L ).Finally, we assume that K is topologically regular and polyconvex. Then K is the union of finitely many compact convex sets with interior points.Hence, for the proof we may assume that K is convex with B (0 , ρ ) ⊆ K for a fixed ρ >
0. Choose 0 < a < ρ/ (2 µ ( L )) and put r = 2 aµ ( L ) < ρ . If x ∈ K , then B ((1 − r/ρ ) x, r ) ⊆ K and B ((1 − r/ρ ) x, r ) contains a point z ∈ a L . Since | x − z | ≤ r + ( r/ρ ) | x | ≤ aµ ( L ) (1 + diam( K ) /ρ ) , we get K ⊆ ( K ∩ a L ) + B (cid:0) , aµ ( L ) (1 + diam( K ) /ρ ) (cid:1) , which completes the argument.Thus Theorems 4.1 and 4.2 and Corollary 4.6 together with Lemma 5.1yield the following result. Corollary 5.2. If K is compact and topologically regular, then lim a → + d bL ( V r,sR ( K ; · ) , V r,sR ( K ∩ a L ; · )) = 0 , lim a → + V r,sR ( K ∩ a L ) = V r,sR ( K ) . If, in addition, K has positive reach, then lim a → + ˆΦ r,sk ( K ∩ a L ) = Φ r,sk ( K ) . (24)21 f K is δ -regular or a topologically regular convex set, then the speed ofconvergence is O ( a ) when r = s = 0 and O ( √ a ) otherwise. The property (24) means that ˆΦ r,sk ( K ∩ a L ) is multigrid convergent forthe class of sets of positive reach as defined in the introduction. A similarstatement about local tensors, but without the speed of convergence, can bemade. We omit this here. We first describe how the number of necessary radii R < R < . . . < R d in(10) can be reduced by one if s = 0 and A = R d . Setting s = 0 and A = R d and subtracting ( r !)Φ r, d ( K ) on both sides of Equation (6) yields Z K R \ K p K ( x ) r dx = V r, R ( K ) − ( r !)Φ r, d ( K ) = ( r !) d X k =1 κ k R k Φ r, d − k ( K ) . (25)As mentioned in Section 5.1, the volume tensor Φ r, d ( K ) can be estimated byˆ φ r, d ( K ∩ a L ). We may take V r, R ( K ∩ a L ) − ( r !) ˆ φ r, d ( K ∩ a L ) as an improvedestimator for (25). This corresponds to replacing the integration domains B ( x, R ) ∩ V x ( K ∩ a L ) in (23) by( B ( x, R ) ∩ V x ( K ∩ a L )) \ V x ( a L ) . This makes sense since V x ( a L ) is likely to be contained in K while the left-hand side of (25) is an integral over K R \ K . The Minkowski tensors can nowbe isolated from only d equations of the form (25) with d different values of R . We now suggest a slightly modified estimator for the Minkowski tensorssatisfying the same convergence results as ˆΦ r,sk ( K ∩ a L ) but where the numberof summands in (23) is considerably reduced. As the volume tensors caneasily be estimated with the estimators in Section 5.1, we focus on thetensors with k < d .Let K be a compact set. We define the Voronoi neighborhood N L (0) of0 to be the set of points y ∈ L such that the Voronoi cells V ( L ) and V y ( L )of 0 and y , respectively, have exactly one common ( d − z ∈ L the Voronoi neighborhood N L ( z ) of z is defined, andthus clearly N L ( z ) = z + N L (0). When L ⊂ R is the standard lattice, N L ( z ) consists of the four points in L that are neighbors of z in the usual4-neighborhood [24]. Define I ( K ∩ a L ) to be the set of points z ∈ K ∩ a L such that N a L ( z ) ⊆ K ∩ a L . The relative complement B ( K ∩ a L ) = ( K ∩ a L ) \ I ( K ∩ a L ) of I ( K ∩ a L ) can be considered as the set of lattice pointsin K ∩ a L that are close to the boundary of the given set K .22e modify (23) by removing contributions from I ( K ∩ a L ) and define˜ V r,sR ( K ∩ a L ; A ) = X x ∈ B ( K ∩ a L ) ∩ A x r Z B ( x,R ) ∩ V x ( K ∩ a L ) ( y − x ) s dy. (26)Assuming that K has positive reach, let 0 < R < R < . . . < R d < Reach( K ). We write again K for K ∩ a L . Then we obtain the estimators ˜Φ r,sd ( K ; A × S d − )...˜Φ r,s ( K ; A × S d − ) = (cid:16) A r,sR ,...,R d (cid:17) − ˜ V r,sR ( K ; A )...˜ V r,sR d ( K ; A ) (27)with A r,sR ,...,R d given by (8).Working with ˜ V r,sR ( K ∩ a L ; A ) reduces the workload considerably. Forinstance, when K is δ -regular or polyconvex and topologically regular, thenumber of elements in I ( K ∩ a L ) increases with a − d , whereas the numberof elements in B ( K ∩ a L ) only increases with a − ( d − as a → + . The set I ( K ∩ a L ) can be obtained from the digital image of K in linear time usinga linear filter. Moreover, we have the following convergence result. Proposition 5.3.
Let K be a topologically regular compact set with positivereach and let C be such that V ( L ) ⊆ B (0 , C ) . If A is a Borel set in R d and aC < R < R < . . . < R d < Reach( K ) and K = K ∩ a L , then ˜Φ r,sk ( K ; A × S d − ) = ˆΦ r,sk ( K ; A × S d − ) for all k ∈ { , . . . , d − } , whenever s = 0 or s is odd. If s is even and k ∈ { , . . . , d − } , then lim a → + ˜Φ r,sk ( K ; A × S d − ) = lim a → + ˆΦ r,sk ( K ; A × S d − ) . Proof.
Let aC < R <
Reach( K ). For x ∈ I ( K ∩ a L ), we have B ( x, R ) ∩ V x ( K ∩ a L ) = V x ( a L ) , so the contribution of x to the sum in (23) is ( s !) x r Φ s, d ( V ( a L )). It followsthat V r,sR ( K ∩ a L ; A ) − ˜ V r,sR ( K ∩ a L ; A ) = ( s !)Φ s, d ( V ( a L )) X x ∈ I ( K ∩ a L ) ∩ A x r . (28)For odd s we have Φ s, d ( V ( a L )) = 0, so the claim follows. For s = 0 theright-hand side of (28) does not vanish, but it is independent of R . Acombination of (cid:16) A r, R ,...,R d (cid:17) − = ( r !) − , s >
0, we have that Φ s, d ( V ( a L )) = a d + s Φ s, d ( V ( L )), while (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) X x ∈ I ( K ∩ a L ) ∩ A x r (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ X x ∈ I ( K ∩ a L ) | x | r ≤ sup x ∈ K | x | r X x ∈ I ( K ∩ a L ) (cid:16) a d H d ( V ( L )) (cid:17) − H d ( V x ( a L )) ≤ sup x ∈ K | x | r · a − d · H d ( V ( L )) − · H d ( K aC ) . Therefore, the expression on the right-hand side of (28) converges to 0.It should be noted that a similar modification for Φ r,sk is not necessary.In fact the modified Voronoi tensor measure (13) with K = K has theadvantage that small Voronoi cells that are completely contained in the R / K ∩ a L do not contribute. In particular, contributionsfrom I ( K ∩ a L ) are automatically ignored when a is sufficiently small. Most existing estimators of intrinsic volumes [17, 19, 24] and Minkowskitensors [28, 30] are n -local for some n ∈ N . The idea is to look at all n ×· · ·× n pixel blocks in the image and count how many times each of the 2 n d possible configurations of black and white points occur. Each configurationis weighted by an element of T r + s and Φ r,sk ( K ) is estimated as a weightedsum of the configuration counts. It is known that estimators of this type forintrinsic volumes other than ordinary volume are not multigrid convergent,even when K is known to be a convex polytope; see [32]. It is not difficultto see that there cannot be a multigrid convergent n -local estimator for the(even rank) tensors Φ , sk ( K ) with k = 0 , . . . , d − s ∈ N , for polytopes K ,either. In fact, repeatedly taking the trace of such an estimator would leadto a multigrid convergent n -local estimator of the k th intrinsic volume, incontradiction to [32].The algorithm presented in this paper is not n -local for any n ∈ N . Itis required in the convergence proof that the parallel radius R is fixed whilethe resolution a − goes to infinity. The non-local operation in the definitionof our estimator is the calculation of the Voronoi diagram. The computationtime for Voronoi diagrams of k points is O ( k log k + k ⌊ d/ ⌋ ), see [5], which issomewhat slower than n -local algorithms for which the computation time for k data points is O ( k ). The computation time can be improved by ignoringinterior points as discussed in Section 5.3.The idea to base digital estimators for intrinsic volumes on an inversionof the Steiner formula as in (9) has occurred before in [16, 23]. In both24eferences, the authors define estimators for polyconvex sets which are notnecessarily of positive reach. This more ambitious aim leads to problemswith the convergence.In [16], the authors use a version of the Steiner formula for polycon-vex sets given in terms of the Schneider index, see [26]. Since its definitionis, however, n -local in nature, the authors choose an n -local algorithm toestimate it. As already mentioned, such algorithms are not multigrid con-vergent.In [23], it is used that the intrinsic volumes of a polyconvex set can, onthe one hand, be approximated by those of a parallel set with small parallelradius, and on the other hand, the closed complement of this parallel set haspositive reach, so that its intrinsic volumes can be computed via the Steinerformula. The authors employ a discretization of the parallel volumes ofdigital images, but without showing that the convergence is preserved.It is likely that the ideas of the present paper combined with the onesof [23] could be used to construct multigrid convergent digital algorithmsfor polyconvex sets. The price for this is that the notion of convergence in[23] is slightly artificial for practical purposes, requiring very small parallelradii in order to get good approximations and at the same time large radiicompared to resolution.In [33], n -local algorithms based on grey-valued images are suggested.They are shown to converge to the true value when the resolution tends toinfinity. However, they only apply to surface and certain mean curvaturetensors. Moreover, they are hard to apply in practice, since they require de-tailed information about the underlying point spread function which specifiesthe representation of the object as grey-value image. If grey-value imagesare given, the algorithm of the present paper could be applied to thresh-olded images, but there may be more efficient ways to exploit the additionalinformation of the grey-values. acknowledgements We wish to thank the referees for carefully reading the paper and mak-ing helpful suggestions for improvements. The first author was supportedin part by DFG grants FOR 1548 and HU 1874/4-2. The third authorwas supported by a grant from the Carlsberg Foundation. The second andthird authors were supported by the Centre for Stochastic Geometry andAdvanced Bioimaging, funded by the Villum Foundation.25 eferences [1] Arns, C.A., Knackstedt, M.A., Mecke, K.R.: Characterising the mor-phology of disordered materials. In: Mecke, K.R., Stoyan, D. (eds.):Morphology of Condensed Matter. Lecture Note in Physics vol 600, pp.37–74. Springer, Berlin (2002)[2] Barvinok, A.: A Course in Convexity. American Mathematical Society,Providence, RI (2002)[3] Beisbart, C., Barbosa, M.S., Wagner, H., da F. Costa, L.: Extendedmorphometric analysis of neuronal cells with Minkowski valuations.Eur. Phys. J. B 52, 531–546 (2006)[4] Chazal, F., Cohen-Steiner, D., M´erigot, Q.: Boundary measures forgeometric inference. Found. Comput. Math. 10, 221–240 (2010)[5] Chazelle, B.: An optimal convex hull algorithm in any fixed dimension.Discrete Comput. Geom. 10, 377–409 (1993)[6] Clarenz, U., Rumpf, M., Telea, A.: Robust feature detection and localclassification for surfaces based on moment analysis. IEEE Trans. Vis.Comp. Graphics 10, 516–524 (2004)[7] Fremlin, D.H.: Skeletons and central sets. Proc. London Math. Soc. 74,701–720 (1997)[8] Federer, H.: Curvature measures. Trans. Amer. Math. Soc. 93, 418–491(1959)[9] Federer, H.: Geometric measure theory. Springer, Berlin (1969)[10] H¨orrmann, J., Kousholt, A.: Reconstruction of convex bodies frommoments. arXiv:1605.06362v1 (2016)[11] Hug, D., Schneider, R.: Local tensor valuations. Geom. Funct. Anal.24, 1516–1564 (2014)[12] Hug, D., Last, G., Weil, W.: A local Steiner-type formula for generalclosed sets and applications. Math. Z. 246, 237–272 (2004)[13] Hug, D., Schneider, R., Schuster, R.: Integral geometry of tensor valu-ations. Adv. in Appl. Math. 41, 482–509 (2008)[14] Kapfer, S.C., Mickel, W., Schaller, F.M., Spanner, M., Goll, C., No-gawa, T., Ito, N., Mecke, K., Schr¨oder-Turk, G.E.: Local anisotropy offluids using Minkowski tensors. J. Stat. Mech: Theory and Exp. 2010P11010 (2010) 2615] Kiderlen, M., Jensen, E.B.V. (eds.): Tensor Valuations and their Ap-plications in Stochastic Geometry and Imaging. Lecture Notes in Math-emaics 2177. Springer (2017).[16] Klenk, S., Schmidt, V., Spodarev, E.: A new algorithmic approach tothe computation of Minkowski functionals of polyconvex sets. Comput.Geom. 34, 127–148 (2006)[17] Klette, R., Rosenfeld, A.: Digital Geometry. Elsevier, San Francisco(2004)[18] Kousholt, A., Kiderlen, M.: Reconstruction of convex bodies from sur-face tensors. Adv. in Appl. Math. 76, 1–33 (2016)[19] Lindblad, J.: Surface area estimation of digitized 3D objects usingweighted local configurations. Image Vis. Comput. 23, 111–122 (2005)[20] McMullen, P.: Isometry covariant valuations on convex bodies. Rend.Circ. Mat. Palermo (2), Suppl. 50, 259–271 (1997)[21] M´erigot, Q., Ovsjanikov, M., Guibas, L.: Voronoi-based curvature andfeature estimation from point clouds. IEEE Trans. Vis. Comp. Graph-ics. 17, 743–756 (2010)[22] Miles, R.E., Serra, J. (eds.): Geometrical Probability and BiologicalStructures: Buffons 200th Anniversary (Proceedings Paris, 1977), Lec-ture Notes in Biomath. 23. Springer, Berlin (1978)[23] Mrkviˇcka, T., Rataj, J.: On the estimation of intrinsic volume densitiesof stationary random closed sets. Stochastic Process. Appl. 118, 213–231 (2008)[24] Ohser, J., M¨ucklich, F.: Statistical Analysis of Microstructures. JohnWiley & Sons, Ltd, Chichester (2000)[25] Ohser, J., Schladitz, K.: 3D Images of Materials Structures: Processingand Analysis. Wiley-VCH, Weinheim (2009)[26] Schneider, R.: Convex bodies: The Brunn–Minkowski Theory. 2nd edn.Cambridge University Press, Cambridge (2014)[27] Schneider, R., Schuster, R.: Tensor valuations on convex bodies andintegral geometry, II. Rend. Circ. Mat. Palermo (2), Suppl. 70, 295–314 (2002)[28] Schr¨oder-Turk, G.E., Kapfer, S.C., Breidenbach, B., Beisbart, C.,Mecke, K.: Tensorial Minkowski functionals and anisotropy measuresfor planar patterns. J. Microsc. 238, 57–74 (2008)2729] Schr¨oder-Turk, G.E., Mickel, W., Kapfer, S.C., Klatt, M.A., Schaller,F. M., Hoffmann, M.J., Kleppmann, N., Armstrong, P., Inayat, A.,Hug, D., Reichelsdorfer, M., Peukert, W., Schwieger, W., Mecke, K.:Minkowski tensor shape analysis of cellular, granular and porous struc-tures. Adv. Mater. 23, 2535–2553 (2011)[30] Schr¨oder-Turk, G.E., Mickel, W., Kapfer, S.C., Schaller, F.M., Breiden-bach, B., Hug, D., Mecke, K.: Minkowski tensors of anisotropic spatialstructure. New J. Phys. 15, 083028 (2013)[31] Schr¨oder-Turk, G.E., Mickel, W., Schr¨oter, M., Delaney, G.W., Saadat-far, M., Senden, T.J., Mecke, K., Aste, T.: Disordered spherical beadpacks are anisotropic. Europhys. Lett. 90, 34001 (2010)[32] Svane, A.M.: On multigrid convergence of local algorithms for intrinsicvolumes. J. Math. Imaging Vis. 49, 352–376 (2014)[33] Svane, A.M.: Estimation of Minkowski tensors from digital grey-scaleimages. Image Anal. Stereol. 34, 51–61 (2015)[34] Svane, A.M.: Local digital algorithms for estimating the integratedmean curvature of rr