ImageAnalStereol2009;28:165-177 OriginalResearchPaper ERRORBOUNDSFOR SURFACE AREAESTIMATORSBASED ON CROFTON'SFORMULA MARKUSKIDERLEN1ANDDANIELMESCHENMOSER2 1DepartmentofMathematicalSciences,UniversityofAarhus ,8000Aarhus,Denmark;2Insitute ofStochastics, Ulm University,89069Ulm, Germany e-mail: kiderlen@imf.au.dk,daniel.meschenmoser@uni-u lm.de (AcceptedSeptember1,2009) ABSTRACT According to Crofton’s formula, the surface area S(A)of a sufficiently regular compact set AinRdis proportional to the mean of all total projections pA(u)on a linear hyperplane with normal u, uniformly averaged over all unit vectors u. In applications, pA(u)is only measured in kdirections and the mean is approximated by a finite weighted sum /hatwideS(A)of the total projections in these directions. The choice of t he weightsdependson the selected quadraturerule. We define an associated zonotope Z(dependingonlyon the projection directions and the quadrature rule ), and show that the relative error /hatwideS(A)/S(A)is bounded from below by the inradius of Zand from above by the circumradiusof Z. Applying a strengthened isoperimetric inequalityduetoBonnesen,we showthattherectangularqua dratureruledoesnotgivethebest possibleerror boundsfor d=2. Inaddition,wederiveasymptoticbehavioroftheerror(w ithincreasing k)intheplanarcase. Thepaperconcludeswith applicationsto surfacearea estim ationin design-baseddigital stereologywherewe showthat the weightsdueto Bonnesen’sinequalityare bette r thanthe usual weightsbased onthe rectangular rule and almost optimal in the sense that the relative error o f the surface area estimator is very close to the minimalerror. Keywords: associated zonotope, Crofton formula, digitiza tion, isoperimetric inequality, minimal annulus, perimeter,surfacearea. INTRODUCTION Onecommonapproachtoapproximatethesurface areaS(A)of an unknown set A⊂Rdfrom its digitization is based on a discretization of Crofton’s formula. We discuss the worst case error introduced by the discretization of the rotational integral in dependence of the quadrature rule chosen. As the methods apply generally to surface area estimators based on Crofton’s formula, we describe them in a general framework and return to its application to digital images in the third section, entitled “Error Boundsfor DigitalSurfaceAreaEstimators” . Throughout the paper a direction is a vector on the unit sphere Sd−1inRd. Ifuis a direction, u⊥ denotesthe linear hyperplane with normal u, ander,u isthestraightlinewith direction uthroughr∈u⊥.Let A⊂Rdbe afull-dimensional compact set in the class UPR; (definitions can be found in the next section ). A special case of Crofton’s formula (Rother and Z¨ ahle, 1990)expressesthe surfacearea S(A)ofAin terms of the Eulercharacteristic χoflinearsections S(A) =2 γd/integraldisplay Sd−1/integraldisplay u⊥χ(A∩er,u)drµ(du).(1) Hereγd= (2κd−1)/(dκd), whereκdis thevolumeof thed-dimensional unit ball, and µis thenormalizedHaar measure on the unit sphere Sd−1; see,e.g., SchneiderandWeil(1992,p.18),butnotethedifferent normalization. ForsetsAthatarenotfull-dimensional, Eq. 1 still holds if its left hand side S(A)is defined in such a way that lower dimensional parts of Aare countedtwice. Theinnerintegralof Eq.1, pu=/integraldisplay u⊥χ(A∩er,u)dr, (2) is called total projection of A in direction u , as it is obtained by measuring the (d−1)-volume of the orthogonal projection of Aonu⊥with multiplicities. In Eq. 2 the integration is understood with respect to the Lebesgue measure on u⊥.If total projections can be determined exactly for finitely many directions u1,...,uk∈Sd−1,say,ak-pointquadraturerulecanbe used to discretize the outer integral in Eq.1 and one obtainstheapproximation /hatwideSd k(A) =2 γdk ∑ i=1cipui, (3) whichdependsonthechoiceofweights c1,...,ck≥0. To assure that the quadrature rule is exact whenever Ais a ball, we assume throughout that the weights sum up to 1. If, for example, the rectangular rule is chosen in the planar case, then the weights are proportional to the arc-lengths of the corresponding 165 KIDERLEN MET AL:Errorboundsforsurfaceareaestimators spherical Voronoi cells generated by {u1,...,uk}on S1. This geometric interpretation generalizes readily to higher dimensions. If {P1,...,Pk}is the spherical Voronoi tessellation of Sd−1generated by the set of projection directions {u1,...,uk}withui∈Pithen the weights cV i=µ(Pi),i=1,...,k, will be called Voronoi weights associated to {u1,...,uk}. These weights are commonly used in applications for d=2,3.In the special case where d=2andu1,...,ukareequidistant,theweightsforthe rectangularquadraturerule(Voronoiweights)coincide with thoseforthe trapezoidalquadraturerule. The discretization of the spherical integral introducesabias,whichtypicallydependsontheset A. We are interested in the worst case behavior. Already Steinhaus(1930)treatedthe specialcasewhere d=2, kis even, and {u1,...,uk}forms an equidistant set of points in S1. For/hatwideS2 k(A), given by Eq. 3 with Voronoi weightscV 1,...,cV k, he derived sharp bounds for the relativeerror: π kcos(π/k) sin(π/k)≤/hatwideS2 k(A) S(A)≤π k1 sin(π/k).(4) The left hand side and the right hand side of Eq.4 arethe endpointsofthe intervalofallpossiblerelative errors, as Avaries.Such aninterval canbe established without the assumption of equidistant directions and in all dimensions. We refer to this interval as error intervalin thefollowing. Using a translative Crofton formula in Section “Error Bounds for /hatwideSd k(A)”, we will define an origin- symmetric convex body Z⊂Rdassociated to the discretization, only depending on the projection directions and the quadrature rule. We will show in Lemma 1 that the relative error /hatwideSd k(A)/S(A)is in a sharp way bounded from below by (a multiple of) the inradius of Zand from above by (a multiple of) the outer radius of Z. Thus, the thickness of the minimal annulus of Zis proportional to the length of the error interval and describes the quality of the estimator. Givenkprojection directions, the quadrature rule (in other words, the values of the associatedweights) that minimizes the minimal thickness of Zcan typically only be determined numerically. In the planar case, we suggest to bound the thickness of the minimal annulus of Zfrom above by an isoperimetric deficit using a strengthened isoperimetric inequality due to Bonnesen.Thisisoperimetricdeficitcanbeminimized withrespecttoallquadraturerulesinclosedform.The weights minimizing the isoperimetric deficit will becalledBonnesen weights and are proportional to the lengths of the edges of a polygon circumscribing the unitdiskandtouchingitexactlyatthepoints u1,...,uk. WewillshowthatVoronoiweightsarenotminimizing the length of the error interval by giving an example wheretheBonnesenweightsyieldbettererrorbounds. We will determine the asymptotic behavior (as k→ ∞) of the relative error for the Bonnesen weights in Theorem 4. At the end of the second section we will consider the case where the directions u1,...,uk∈S1 areobtainedusingsystematicrandomsamplingon S1. We will show that the coefficient of error of /hatwideS2 k(A) can be bounded from above by a geometric quantity involving Z,namely a multiple of the L2-distance betweenZandits Steinerball. Inthe subsequentsection we discuss error bounds for digital surface area estimators . The digitization of Aon a randomly translated, rectangular grid will be considered.Asymptotic boundsforthe expectedvalue oftheestimatorfor S(A)inthegridwillbeestablished. ”Asymptotic” relates here to increasing resolution of the grid. The vectors u1,...,ukare chosen as grid directions, i.e.,normalizedvectorsconnectingtwogrid points. These two grid points are usually neighbours andwewillconsiderthe 4-,8-and16-neighborhoodin 2Dandthe6-and26-neighborhoodin3D.Forallthese settings the Voronoi and Bonnesen weights together with the corresponding in- and circumradii randR, respectively, will be computed analytically, exceptfor the 26 directions in 3D where numerical methods will be used. We will compare the relative errors with the minimal error achievedby numerically optimizing the weights and show that the Bonnesen weights lead, at least in 2D, to smaller errors than the widely used Voronoi weights. We then restrict to quadratic grids in the plane and consider boundary length estimators based on pairs of grid points that are contained in some (n−1)×(n−1)square of grid cells, n≥2. This generalizes the case n=2, which corresponds to the boundary length estimator based on 8-neighbours. Theorem7considerssuchestimatorsforgeneral n≥2 and shows that the asymptotic mean relative error of /hatwideSd k(A)for Bonnesenweightsdecreasesas n−2. The application of Bonnesen’s improved isoperimetric inequality restricts many of the above arguments to the two-dimensional case. In the last section we discuss the possibility of extensions to higherdimensions. ERRORBOUNDS FOR /hatwideSd k(A) A setA⊂Rdisfull-dimensional , if its tangent cone ataspans Rdfor almost all a∈Awith respect 166 ImageAnalStereol2009;28:165-177 to(d−1)-dimensional Hausdorff measure. If Ais topologically regular ( Ais the closure of its interior) and convex, it is also full-dimensional. The set Ahas positive reach if there is anε>0 such that each point in theε-neighborhood of Ahas a unique closest point inA. Throughout the following we assume that Ais an element of the family UPRof all sets in Rdwhich can be written as a finite union of compact sets A1,...,Am with positive reach such that any intersection/intersectiontext i∈IAi withI⊂ {1,...,m}is either empty or a set of positive reach, as well. In particular, convex bodies (nonempty compact convex subsets of Rd) andpolyconvex sets (finite unions of convex bodies) are elements of UPR. ForA∈UPRthe surface area measure S(A,·)of orderd−1 is defined, and Eq. 1 holds with S(A) = S(A,Sd−1). IfAis in addition full-dimensional, S(A) coincides with the usual surface area of A.In view of theapplicationsindigitalstereology,itisconvenientto extendthetotalprojectionmapping v/ma√sto→pvtoRd\{o} by positive homogeneity of degree 0. The translative Crofton formula pv=1 2/bardblv/bardbl/integraldisplay Sd−1|/a\}bracketle{tu,v/a\}bracketri}ht|S(A,du)(5) holds for almost all v∈Rd; see Rataj (2002, Theorem 2.1 and Theorem 2.3) , where /a\}bracketle{tu,v/a\}bracketri}htis the usual inner product of uandv. IfAis polyconvex, Eq.5 holds for all v∈Rd\{o}. Ifv1,...,vk∈Rd\ {o}are such that Eq.5 holds withv=vifor alli=1,...,k, then the definition of /hatwideSd k(A)in combinationwith Eq.5yields /hatwideSd k(A) =1 γd/integraldisplay Sd−1h(u)S(A,du),(6) with h:=k ∑ i=1ci /bardblvi/bardbl|/a\}bracketle{tvi,·/a\}bracketri}ht|. (7) The key observation is that the integrand his the supportfunctionofaconvexbody.Wereferthereader to Schneider (1993) for relevant notions and concepts inconvexgeometryandonlyrecallthemostimportant facts here. The support function hKof a convex body Kis givenby hK(u) =max{/a\}bracketle{tx,u/a\}bracketri}ht:x∈K},u∈Sd−1. Here and in the following, we consider the support function as a function on the unit sphere. For convex bodiesKandMandscalars α,β≥0, wehave αhK+βhM=hαK⊕βM, (8)where the Minkowski addition ⊕of sets and the multiplication of a set with a scalar are understood pointwise. We will repeatedly use the monotonicity property, K⊂M⇐⇒hK≤hM, (9) and the fact that the support function hBdof the Euclidean unit ball BdinRdis the constant 1. Eq. 9 implies in particular that any convex body is uniquely determined by its support function. Consequently, the definition, δ2 2(K,M):=/integraldisplay Sd−1(hK(u)−hM(u))2du, for convex bodies KandM, givesrise to the so-called L2-metricδ2(·,·)on the family of convex bodies. The support function of the line segment [−x,x]with endpoints −xandx∈Rdis|/a\}bracketle{tx,·/a\}bracketri}ht|. Due to Eq.8, the function hinEq.7 is the support function of a finite sum Zof line segments. Such sets are called zonotopes and play a prominent role in functional analysis, convex and stochastic geometry (see, e.g., Goodey and Weil, 1993 and the references therein). Explicitly, wehave Z=c1[−u1,u1]⊕...⊕ck[−uk,uk],(10) with the unit vectors ui=vi//bardblvi/bardblfori=1,...,k. In view of Eq.6 the approximation /hatwideSd k(A)can be expressedin termsofthe associatedzonotope Z,as /hatwideSd k(A) =1 γd/integraldisplay Sd−1hZ(u)S(A,du).(11) To obtain lower and upper bounds of /hatwideSd k(A), we have to find maxima and minima of hZ. Due to Eq.9, r≥0 is the minimum of hZonSd−1if and only ifrBdis the largest ball contained in Z. Similarly, R≥0 is the maximum of hZ, if and only if RBdis the smallest ball containing Z. With these optimal values of 0≤r≤R, the setRBd\rBdis called the minimal annulus of Z . The difference R−ris called the width oftheminimalannulus anddenotedby T(Z).Forlater reference we summarize this geometric interpretation forpolyconvexsets(forwhich Eq.5holdsforarbitrary v/\e}atio\slash=o). As formulations for UPR-sets are obtainedin a straightforward manner,wewillrestrict topolyconvex setsfrom nowon. Lemma1 Let A ⊂Rdbe a polyconvex set with positive surface area, and fix k ≥2and v1,...,vk∈ Rd\ {o}. Letγd= (2κd−1)/(dκd). If/hatwideSd k(A)is given byEq.3, thenthesharpbounds r γd≤/hatwideSd k(A) S(A)≤R γd, (12) for the relative estimation error hold, where r is the smaller and R is larger radius of the minimal annulus ofthezonotopeZ givenby Eq.10. 167 KIDERLEN MET AL:Errorboundsforsurfaceareaestimators ProofAsRis the maximum of hZonSd−1, it follows from Eq.11,that /hatwideSd k(A)≤R γdS(A), which yields the upper bound. The lower bound follows analoguously from Eq. 11 and the fact that r is theminimum of hZonSd−1. That the bounds in the above Lemma are sharp follows from thenextexample. Example2 Fix k≥2, u1,...,uk∈Sd−1and weights c1,...,ck≥0for a quadrature rule. Define Z according to Eq.10. Due to symmetry, the ball rBd touches the boundary of Z in at least two antipodal points rw and −rw, w ∈Sd−1. Let A be a ball of (d−1)-volume1/2in the hyperplane w⊥. (As A is lower dimensional, the proper interpretation of S (A) istwicethe(d−1)-dimensional Hausdorff measure, so S(A) =1.) The surface area measure of A is concentratedonthepointswand −w,andh Zcoincides in bothofthesedirectionswith r,so Eq.11implies /hatwideSd k(A) =1 γd/integraldisplay Sd−1rS(A,du) =r γdS(A), andequalityholdsontheleft handsideof Eq.12. Fig.1illustratesthisford =2,k=2,u1= (1,0)⊤, u2= (0,1)⊤and c1=c2=1/2. Obviously A is not topologically regular, but it can be approximated by a sequence of topologically regular convex bodies (Am)in such a way that limm→∞/hatwideSd k(Am)/S(Am) = /hatwideSd k(A)/S(A). This implies that the left hand side ofEq.12 cannot be improved, even if we restrict considerations to topologically regular sets. To show thatthe secondinequalityin Eq.12is sharp,asimilar argument can be used, if ±w are directions for which hZbecomesmaximal,andthuscoincideswith R. Fig.1.ApossiblesetAforthecasewhereZ istheunit cube;seeExample2.THETWO-DIMENSIONALCASE In the following, we will restrict to the case d= 2, although some of the concepts can be transferred to higher dimensions. As the aim is to minimize the length of the error interval of /hatwideS2 k(A)in Eq.6, the difference R−rshould be as small as possible. This can be achieved by an appropriate choice of the weightsc1,...,ck. To obtain an exact value for the integralin Eq.6inthecasewhere Aisadisk,wemust assumethattheweightssumuptoone. Itfollowsfrom Schneider (1993) , thatc1+...+ck=1 is equivalent to the condition that the zonotope Zgiven by Eq.10 has perimeter 4. Let Zbe the family of all zonotopes that can be written as sum of line-segmentsparallel to given unit vectors u1,...,uk. LetZ4be the family of thoseZ∈Zthat have perimeter 4. We are therefore facedwith the problem of findinga zonotope Z∗∈Z4 thatsatisfies T(Z∗) =min{T(Z):Z∈Z4}.(13) If Z∗=c∗ 1[−u1,u1]⊕...⊕c∗ k[−uk,uk], thenc∗ 1,...,c∗ k≥0 are the bestweights in Eq. 3, in the sense that among all weights summing up to one they yield the shortest interval of possible relative errors. A solution Z∗of the optimization problem in Eq.13 alwaysexistsduetoacompactnessargumentbasedon theBlaschkeselectiontheorem. For asymptotic results, it is enough to replace the objective function in Eq.13 by a simpler one. Bonnesen (1929) improved the isoperimetric inequality for an arbitrary planar convex body K, statingthat S2(K) 4π−V(K)≥π 4T2(K), (14) whereS(K)andV(K)are perimeter and area of K, respectively.For K=Z∈Z4wehaveS(Z)=4andthe lefthandsideof Eq.14isminimalforthezonotope /tildewideZ∈ Z4that has the greatest area. According to a classical result of Lindel¨ of (1869), /tildewideZis characterized among all zonotopes in Z4by the fact that it circumscribes a circle. Due to origin-symmetry, this circle is the incircleof/tildewideZ,centeredattheorigin,andwithradius /tildewiderk. Thisallowsanexplicitconstructionof /tildewideZ.Uptoscaling with the factor 1 //tildewiderkthe zonotope/tildewideZcoincideswith the polytope /tildewideP:=k/intersectiondisplay i=1/braceleftbig x∈R2:|/a\}bracketle{tui,x/a\}bracketri}ht| ≤1/bracerightbig ,(15) 168 ImageAnalStereol2009;28:165-177 obtained by intersecting all supporting half-planes of the unitdiskwith outernormalin {±u1,...,±uk}.We now assumewithoutloss ofgenerality thatthe vectors u1,...,ukall are located on the positive half-sphere {(cosϕ,sinϕ):0≤ϕ<π}andareorderedaccording toincreasing angles. We write <) (u,v)∈[0,π]for the (smaller)anglebetweentheunitvectors uandv.Letαi betheouterangleofthevertexbetweenthe ithandthe (i+1)st edge (i.e.,αiis the angle of the normal cone at this vertex, in other words π−αiis the usual inner angle;seeFig.2). Explicitly, wehave αi=/braceleftBigg <) (ui,ui+1),i=1,...,k−1, π−<) (u1,uk),i=k,(16) asthe normalofthe (k+1)stedgeis −u1. Fig. 2.Construction of the angles αiand the polytope /tildewideP givenby Eq.15with k =4. As the length of the ith edge of/tildewidePis tan(αi/2)+ tan(αi−1/2)andS(/tildewideZ) =4,we obtain /tildewideZ=/tildewidec1[−u1,u1]⊕...⊕/tildewideck[−uk,uk], with /tildewideci=/tildewiderktan(αi/2)+tan(αi−1/2) 2,(17) and /tildewiderk=/parenleftBigg k ∑ i=1tan(αi/2)/parenrightBigg−1 , (18) where we have put α0=αk. This and V(/tildewideZ) =2/tildewiderk wasalsoderivedbyKnebelman(1941).Theweights /tildewideci are based on the application of Bonnesen’s inequality Eq.14 and will be called Bonnesen weights in the following.The ithweight/tildewideciistherelativelengthofthe ith edge of the polygon with facet normals ±uiwhich circumscribesacircleofradius /tildewiderk.Theouterradius /tildewideRkof/tildewideZis the largest distance of a vertex of /tildewideZfrom the origin, andthisis /tildewideRk=/tildewiderkkmax i=11 cos(αi/2). (19) Summarizing, we have shown the following. If /hatwideS2 k(A) is an estimator of S(A)>0 given by Eq.6 with the Bonnesen weights ci=/tildewideci,i=1,...,k, fromEq.17, thenthe relativeerrorsobey π 2/tildewiderk≤/hatwideS2 k(A) S(A)≤π 2/tildewideRk. (20) Theseerrorboundsaresharp;seeExample2. It should be noted that Eq.14 is always a strict inequality unless Kis a disk. As Eq.14 is used for zonotopeshere,thisapproachwillnotnecessarilylead to the optimal choice of the weights. However, the choice may be better than choices for the weights motivatedbyusualquadraturerules. Example3 We consider the integrand in Eq.6 with k=3and the directions u i=/parenleftbig cos(ϕi),sin(ϕi)/parenrightbig , i= 1,2,3, where ϕ1=0,ϕ2=π 16,ϕ3=π 8. As mentioned before, the rectangular quadrature rule leads to the Voronoi weights cV i: The ith weight is the normalizedlengthoftheVoronoiarccorrespondingto ui(arcin S1ofallpointscloserto u ithanto anyother pointin {±u1,±u2,±u3}). Thisgives cV i=ϕi+1−ϕi−1 2π, where we assumed π-periodicity. For the present example,we obtain cV=/parenleftbigg15 32,1 16,15 32/parenrightbigg , and the corresponding zonotope Z Vhas inradius rV≈0.18290and circumradius R V≈0.98199; see Fig. 3. Thus, the width of the minimal annulus is approximately 0.79909. Using instead the Bonnesen weights yields approximately /tildewidec= (0.49057 ,0.01885 ,0.49057 ), leadingtotheinradius /tildewider≈0.191412andcircumradius /tildewideR≈0.981147, respectively. The width of the minimal annulus is now approximately 0.789735, which is an improvementofabout 1%. 169 KIDERLEN MET AL:Errorboundsforsurfaceareaestimators ±1±0.50.51 ±1 ±0.5 0.5 1 Fig. 3.The zonoid from Example 3 associated to Voronoiweightstogetherwith its minimalannulus. To formulate an asymptotic result, we have to specify how close the set of the directions u1,...,uk isto asetofequidistantdirections.FollowingGardner et al.(2006), we introduce the symmetrizedspread Δ∗ k ofu1,...,ukby Δ∗ k=max u∈S1min 1≤i≤kmin{/bardblu−ui/bardbl,/bardblu−(−ui)/bardbl}. Geometrically, Δ∗ kis the maximal distance of a unit vector from the set {±u1,...,±uk}. In particular, {±u1,...,±uk}is aΔ∗ k-net in S1. Forαidefined by Eq.16,we have 2sinαi 4≤Δ∗ k,i=1,...,k.(21) The following theorem shows that the choice ci= /tildewidecileads to a relative error of /hatwideS2 k(A)that depends quadratically on Δ∗ k. Here we only consider sampling sets{u1,...,uk}such that every closed sub-arc of S1of length π/2 contains at least one point of {±u1,...,±uk}. Equivalently, Δ∗ k≤/radicalbig 2−√ 2. Theorem4 Let A ⊂R2be a polyconvex set with positive perimeter. Let k ≥2and{v1,...,vk} ⊂R2\ {0}such that the symmetrized spread of the vectors ui=vi//bardblvi/bardbl,i=1,...,k,isΔ∗ k≤/radicalbig 2−√ 2.If/hatwideS2 k(A)in Eq.3iscalculatedusingtheBonnesenweightsc i=/tildewideci, i=1,...,k,fromEq.17,thenthe relativeerrorobeys /vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/hatwideS2 k(A)−S(A) S(A)/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle≤π2 3(Δ∗ k)2.(22) ProofFromEq.20we get π/tildewiderk−2≤2/parenleftBigg/hatwideS2 k(A)−S(A) S(A)/parenrightBigg ≤π/tildewideRk−2.(23)We estimate the left hand side of Eq.23. In view of Eq.21, we have αi/2≤2arcsin/parenleftbig Δ∗ k/2/parenrightbig ≤π/4 for all i=1,...,k.Taylor’stheoremimplies tan(αi/2)≤αi/2+c′(αi/2)3,foralli=1,...,k, wherec′=8/3is the third derivative of tan (x)/3! evaluatedat π/4.Relations18,21andarcsin (x)≤π 2x, 0≤x≤1, imply that π/tildewiderk≥π ∑k i=1/parenleftBig αi 2/parenleftBig 1+c′/parenleftbigαi 2/parenrightbig2/parenrightBig/parenrightBig≥2 1+c′π2 4/parenleftbig Δ∗ k/parenrightbig2, andthisgives π/tildewiderk−2≥ −2π2 3(Δ∗ k)2. The right hand side of Eq.23 can be estimated in an even easier way using the fact that the perimeter of the incircle of /tildewideZis bounded by S(/tildewideZ) =4 and hence /tildewiderk≤2/π. Togetherwith Eq.19this gives π/tildewideRk−2≤2 1−/parenleftbig Δ∗ k/parenrightbig2(Δ∗ k)2≤2√ 2−1(Δ∗ k)2, asΔ∗ k≤/radicalbig 2−√ 2.Putting thingstogetherwearriveat /vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/hatwideS2 k(A)−S(A) S(A)/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle≤c(Δ∗ k)2 withc=max/braceleftbig π2/3,1/(√ 2−1)/bracerightbig =π2/3. Using the fact that Voronoi weights deviate only slightly from Bonnesen weights as kincreases, it can be shown that the same order of convergence also holds for the relative error of /hatwideS2 k(A)if the estimator is basedonVoronoiweights. Theexampleofequidistant sampling shows that quadratic behavior is the best possible. Example5 Considerthespecialcasewhereu 1,...,uk are equidistant on the upper half circle, meaning that ui= (cos(iπ/k),sin(iπ/k)), i=1,...,k.Hence Δ∗ k=2sinπ 4k∼π 2k,k→∞. By symmetry arguments, the weights leading to the minimal width of the corresponding minimal annulus must all be equal and thus c 1=...=ck=1/k and ci=/tildewidecifor i=1,...,k. According to Eqs.18 and 19, the inner and outer radii of the associated zonotope are /tildewiderk=/parenleftBig ktan/parenleftBigπ 2k/parenrightBig/parenrightBig−1 , 170 ImageAnalStereol2009;28:165-177 and /tildewideRk=/tildewiderk/parenleftBig cos/parenleftBigπ 2k/parenrightBig/parenrightBig−1 =/parenleftBig ksin/parenleftBigπ 2k/parenrightBig/parenrightBig−1 , cf.Eq.4. Therefore, the width of the minimal annulus is /tildewideRk−/tildewiderk=π 4k−2+O/parenleftbig k−4/parenrightbig ,k→∞. This shows that /tildewideRk−/tildewiderk, and thus the relative worst caseerror,areoforder/parenleftbig Δ∗ k/parenrightbig2. Insteadofusing theabovegeometricargumentsto obtain asymptotics for the worst case error, one might also use methods from optimum quantization (see Gruber, 2004). Among other important applications, this theory yields asymptotic minimum errors of numericalintegrationforclassesofH¨ oldercontinuous functions. As the function gu:v/ma√sto→ |/a\}bracketle{tv,u/a\}bracketri}ht|, and hence the function pvinEq.5 are Lipschitz continuous, optimum quantization gives an upper bound for the worst case error depending linearly on Δk. This suboptimal rate is due to the fact that the class of H¨ older continuous functions with H¨ older exponent 1 is considerably larger than its subspace spanned by/braceleftbig gu:u∈S1/bracerightbig . ASEMI-RANDOMIZED APPROACH The associated zonotope for quadrature rules can also be used in the context of a semi- randomized approach, which generalizes systematic random sampling designs. The idea of this design based approach is to evaluate the total projections of the randomly rotated set ϑAinkdirections. In other words, given kvectorsv1,...,vk∈Sd−1and weights c1,...,ck, theestimatorfor S(A)is definedby /hatwideSd k(ϑA) =2 γdk ∑ i=1cipϑ−1vi, (24) whereϑis a random rotation whose distribution is the normalized Haar measure on the compact group SOdof proper rotations. Clearly, Eq.24 defines a randomvariable and Crofton’s formula implies that this variable is an unbiased estimator for S(A). In particular, if d=2 and the set {±v1,...,±vk}is equidistant in S1, the estimator /hatwideSd k(ϑA)inEq.24 is the one obtained from systematic random sampling. Moran (1966) considered this special case and gave worst case bounds for the variance of /hatwideSd k(ϑA). His approach allows a geometric interpretation which is notrestricted to the planarsetting:Let Zbe,again,thezonotope associated to a fixed quadrature rule. From Eq.11andthe unbiasednessofthe estimator, weget Var/parenleftBig /hatwideSd k(ϑA)/parenrightBig =Eϑ/parenleftBig /hatwideSd k(ϑA)−S(A)/parenrightBig2 =γ−2 dEϑ/parenleftbigg/integraldisplay Sd−1(hZ(ϑu)−γd)S(A,du)/parenrightbigg2 . Hence Var/parenleftBig /hatwideSd k(ϑA)/parenrightBig = γ−2 d/integraldisplay Sd−1/integraldisplay Sd−1J(u,v)S(A,du)S(A,dv),(25) where J(u,v) =Eϑ((hZ(ϑu)−γd)(hZ(ϑv)−γd)). H¨ older’sinequalityimplies that J(u,v)≤/parenleftBig Eϑ(hZ(ϑu)−γd)2/parenrightBig1/2 ×/parenleftBig Eϑ(hZ(ϑv)−γd)2/parenrightBig1/2 =/integraldisplay Sd−1(hZ(u)−γd)2µ(du). Hence Var/parenleftBig /hatwideSd k(ϑA)/parenrightBig ≤S(A)2 γ2 dϖdδ2 2/parenleftBig Z,γdBd/parenrightBig ,(26) whereδ2(·,·)denotes the L2-metric defined earlier. This inequality is sharp, as equality holds here for example whenever Ais asetof dimension d− 1. The quadrature rule is exact when A=Bdand thus/hatwideSd k/parenleftbig Bd/parenrightbig =dκd, which implies that 2 γdis the mean width b(Z)ofZandγdBdis theSteiner ballofZ; see Schneider (1993, p. 353). Hence, the coefficient of variation/radicalbigg Var/parenleftBig /hatwideSd k(ϑA)/parenrightBig /S(A) is bounded from above by a multiple of the L2- distance of Zto its Steiner ball. Again, a geometric inequality (Groemer, 1990) would allow to give an upper bound of the right hand side of Ineq. 26. However,a more direct evaluationis possibleand was carried out by Jan´ aˇ cek (2001). In particular he found optimal weights c1,...,ckto minimize the variance in the case of a “totally anisotropic object”, i.e., a setAcontained in a hyperplane. The optimal weights are obtained by inverting the covariance matrix of (pϑ−1u1,...,pϑ−1uk). 171 KIDERLEN MET AL:Errorboundsforsurfaceareaestimators ERRORBOUNDSFOR DIGITAL SURFACEAREAESTIMATORS In this section we consider digitizations of a topologically regular polyconvex setA⊂Rdon rectangular grids and assume that the directions uiare givenas normalized difference vectors of grid points, whichusuallyareneigbours.Forexample,weconsider thecommon 4- and 8-neighborhoods in the plane. We compute Voronoi- and Bonnesen-weightstogether with the associated in- and circumradii explicitly and give asymptotic error bounds for /hatwideS(A)for increasing resolutionofthe digitization. To digitize Awe consider the rectangular point grid G=ζ1Z×ζ2Z×...×ζd−1Z×Z, where ζ1,...,ζd−1>0arethegriddistancesinthedirections of the axes, and we used the griddistance in the last coordinate direction as unit. A grid cell is any d-dimensional cuboid z+ [0,ζ1]×...×[0,ζd−1]× [0,1],z∈G. Motivated by design based stereology, we consider digitizations of Aby a randomly translatedgrid.With therandomvariable ξ,uniformly distributed in an arbitrarily chosen grid cell, the random grid ξ+Gis a stationary random closed set and is called a stationary grid in Kiderlen and Rataj (2006). In order to increase resolution, we scale the grid byafactor t>0anddenotethedigitizationof Ainthe scaled grid t(ξ+G)byΔt(A). LetQbe a non-empty compact set, called the sampling element . We assume thateachgridpoint x∈t(ξ+G)isthecenterofasmall sampling window x+tQ, which can be thought of as a pixel or voxel. The pixel digitization consists of all grid points xfor which this sampling window x+tQ hits the set A. HenceΔt(A) =/parenleftbig A⊕tˇQ/parenrightbig ∩t(ξ+G), whereˇQis the reflection of Qat the origin. For Q={o}the pixel digitization reduces to the Gauss digitization (sometimescalled hit-or-missdigitization ) containing all points of the scaled grid in A. All results on error bounds in this section will be stated for the pixel digitization and therefore also hold for the Gauss digitization. In image processing, the term digitization often denotes the union of pixels (grid cells) which are centered at grid points in Δt(A). As such pixel unions are in one-to-one correspondence with the unions of their centers, one may equivalently consider digitizations as subsets of the grid, and we will dosothrougoutthe paper. We fix a set Aand estimate its surface area S(A) fromtheinformationavailableinitsdigitization Δt(A) using a discretized Crofton formula; cf. Ohser and M¨ ucklich (2000). We will focus on the asymptotic error of this estimator when the grid spacing getsfiner,i.e., when the digitized set Δt(A)becomes a betterapproximationoftheoriginalset A.Thefunction pvgiven by Eq.2 can easily be estimated from the digitized set Δt(A)by comparing the values of neighboring points, if vis a gridpoint. For any vector v∈G\{o}suchanestimatorisgivenby /hatwidepv(t) = td−1 /bardblv/bardbl#{x∈t(ξ+G):x∈Δt(A),x+tv/∈Δt(A)}. This estimator counts the number of points xin the digitized set Δt(A)such that x+tvdoes not lie in the digitization of A. From Kiderlen and Rataj (2006, Theorem 5) it follows that this estimator is asymptoticallyunbiased, i.e., lim tց0E/hatwidepv(t) =pv. (27) Havingchosen kvectorsv1,...,vk∈Gandscalars c1,...,ck≥0,onecandefine /hatwideSd k(A;t) =2 γdk ∑ i=1ci/hatwidepvi(t), (28) and it follows from Eq.27 that/hatwideSd k(A;t)is an asymptoticallyunbiasedestimatorfor /hatwideSd k(A),astց0. Note that the estimator given by Eq.28 can be calculated from the knowledge of the digitization Δt(A)ofAalone./hatwideSd k(A;t)behavesapproximately like a discretized Crofton integral, when tis small. Thus, the methods and results of the previous section can be appliedto obtain asymptoticerrorbounds. To illustrate this approach we discretize the Crofton integral using only directions parallel to the coordinate axes: in Rdwe choose the 2 dgrid points vi=−vd+i=ζiei,i=1,...,d−1andvd=−v2d=ed whereeidenotesthe ith unit vector. Due to symmetry, the weights leading to a minimal error interval are all equal,ci=1/(2d),i=1,...,2d, and coincide with the Voronoi weights. The zonotope Zdefined in Eq. 10 is given by Z= [−1/d,1/d]dand it has inradius r=1/dand circumradius R=1/√ d. Lemma 1 and theasymptoticunbiasednessof /hatwideSd 2d(A;t)imply 1 dγd≤lim tց0E/hatwideSd 2d(A;t) S(A)≤1√ dγd. In theplanarcasewehave γ2=2/πand 0.785≈π 4≤lim tց0E/hatwideS2 4(A;t) S(A)≤π 4√ 2≈1.111, 172 ImageAnalStereol2009;28:165-177 which means the asymptotic relative error is 21 .5% in the worst case. In three dimensions, the asymptotic relative erroris atmost33 .3%asγ3=1/2and 0.667≈2 3≤lim tց0E/hatwideS3 6(A;t) S(A)≤2 3√ 3≈1.155. Due to Stirling’s formula we have√ dγd→/radicalbig 2/πas d→∞. As√ dγdis decreasingin d,wehave 0≤lim tց0E/hatwideSd 2d(A;t) S(A)≤/radicalbigg π 2≈1.253, for alld, where/radicalbig π/2 is the best upper bound that holdsfor arbitrary dimension d. We do not obtain a non-trivial uniform lower bound, as there are d∈N andsetsAinRdwithS(A) =1,butsuchthat /hatwideSd 2d(A)is arbitrarily closeto 0. Due to the large worst caseerror eveninlowdimensions,theabovechoiceof gridpoints isnotrecommendedforpracticalapplications.Instead, a larger number of grid pointsshould be used. On the otherhand,theestimatorof S(A)inEq.28isbasedon asymptotic considerations and becomes less reliable whenthelengthsofthevectors viarelarge.Onewould therefore restrict to vectors with bounded length, i.e., vectors contained in a small Euclidean disk around the origin. For computational reasons it is easier to replacetheEuclideandiskbyadiskwithrespecttothe maximumnorm. Tobespecific,wechoose {v1,...,vk} in the set V(d) n:= [−ζ1(n−1),ζ1(n−1)]×... ×[−ζd−1(n−1),ζd−1(n−1)]×[−(n−1),n−1]\{o}. Forn=2, the most common choice in applications,we have V(d) 2={−ζ1,0,ζ1}×... ×{−ζd−1,0,ζd−1}×{−1,0,1}\{o},(29) andk=#V(d) 2=3d−1. We determine asymptotic worstcaseerrorsfor n=2andn=3intheplanarcase andforn=2in dimension d=3. THETWO-DIMENSIONALCASE In the planar case for n=2, thek=32−1= 8grid points of the set in Eq. 29 are just the 8-neighbours of the origin, and the corresponding directions are parallel to the edges and diagonals of a grid cell. Both the Voronoi weights and the Bonnesen weights can be computed analytically. It turns out that the Voronoi weights and the Bonnesen weights do not coincide, and the corresponding in-and circumradii are different. In the following let β= arccos (ζ1/√ ζ2 1+1). For the Voronoi case, the inradius rVor(ζ1)andthecircumradius RVor(ζ1)aregivenby rVor(ζ1) =  1 πβ+ζ1 2√ ζ2 1+1,ifζ1≤1, 1 2−1 πβ+1 2√ ζ2 1+1,otherwise, and RVor(ζ1) =  1 π/bracketleftBigg β2+/parenleftbigg π 2−β+π 2√ ζ2 1+1/parenrightbigg2/bracketrightBigg1/2 , ifζ1≤1, 1 π/bracketleftBigg/parenleftbigg β+π 2ζ1√ ζ2 1+1/parenrightbigg2 +/parenleftbigπ 2−β/parenrightbig2/bracketrightBigg1/2 , otherwise, respectively. In the Bonnesen case the inradius rBon(ζ1)is givenby rBon(ζ1) =ζ1/parenleftbigg 2/radicalBig ζ2 1+1(ζ1+1)−2/parenleftbig ζ2 1+1/parenrightbig/parenrightbigg−1 andthecircumradius RBon(ζ1)is givenby RBon(ζ1) =  rBon(ζ1)√ 2/parenleftbigg 1+/parenleftBig 1 ζ2 1+1/parenrightBig−1/2/parenrightbigg−1/2 , ifζ1≤1 rBon(ζ1)√ 2/parenleftBig 1+/parenleftbig ζ2 1+1/parenrightbig−1/2/parenrightBig−1/2 , otherwise. Neither the Voronoi nor the Bonnesenweights are optimal in the sense that they minimize the length of the asymptotic error interval. The optimal weights were found by numerically solving a constrained minimizationproblemusingM ATLAB.Itturnsoutthat the Bonnesen weights are very close to the optimum as can be seen from Fig. 4 which shows the thickness of the minimal annulus for Voronoi, Bonnesen, and optimized weights, respectively,dependingon ζ1>0. Interchanging the x- and the y-axis, if necessary, we mayrestrict to ζ1≤1. 173 KIDERLEN MET AL:Errorboundsforsurfaceareaestimators Fig. 4.Thickness of the minimal annulus for eight directional vectors on a rectangular grid in 2D dependingon ζ1forVoronoi,Bonnesen,andoptimized weights,respectively(seethe textfordetails). Fig. 5.Thickness of the minimal annulus for 16 directional vectors on a rectangular grid in 2D dependingon ζ1forVoronoi,Bonnesen,andoptimized weights,respectively(seethe textfordetails). It is easy to see that in the case of a square grid (ζ1=1) the width R−rof the minimal annulus is minimal if and only if all weights are equal, c1= ...=c8=1/8, a choice which coincides with the VoronoiandtheBonnesenweights.Thecorresponding zonotope Zis a regular octagon with side length 1/2 and facets parallel to the vectors v1,...,v8.Z has circumradius R=/radicalbig 4+2√ 2/4 and inradius r=/parenleftBig 1+√ 2/parenrightBig /4.Withγ2=2/πweget 0.948≈π 8/parenleftBig 1+√ 2/parenrightBig ≤lim tց0E/hatwideS2 8(A;t) S(A) ≤π 8/radicalBig 4+2√ 2≈1.026. ThiswasalsoobtainedinKiderlen andJensen(2003).We consider also the case n=3 in the plane, i.e., we use all 16 directions obtained by normalizing the grid points in the cuboid [−2ζ1,2ζ1]×[−2,2]\ {0}. Now even in the special case when ζ1=1 the 16 weights leading to a shortest asymptotic error interval are not equal. We refrain from explicitly stating the formulas for the in- and circumradii because of their complexity. The results are qualitatively similar to the casen=2. The use of the Bonnesen weights yields a smaller thickness of the minimal annulus than the use of the Voronoi weights. The minimal thickness achievedwithoptimizedweightsisonlyslightlybetter thanforBonnesenweights(Fig. 5). THETHREE-DIMENSIONAL CASE In three dimensions, we restrict considerations to the cubic grid ( ζ1=ζ2=1) andn=2.The directions associated to the k=33−1=26 grid points viof the set in Eq. 29 consist of the 6 directions along the edges ( i=1,...,6), 12 diagonals of the faces (i=7,...,18) and 8 spatial diagonals ( i=19,...,26) of the unit cube [0,1]3. Hence, the surface area estimator in Eq. 28 is based on comparison of pixels with neighbors in the so-called 26-neighborhood. The Voronoiweights canbe calculatedasthe relative sizes of the Voronoi cells on the unit sphere generated by{v1//bardblv1/bardbl,...,vk//bardblvk/bardbl}.They can be derived by computing the area of the Voronoi cells on the unit sphere analytically with the help of spherical trigonometry,andaregivenby ci=  1 2−2 πarccos/parenleftbigg√ 2+√ 3√ 2√ 3−√ 3sin/parenleftbigπ 8/parenrightbig/parenrightbigg ≈0.0457779 , ifi=1,...,6, 1 2−1 πarccos/parenleftbigg√ 6−2 2√ 3−√ 6sin/parenleftbigπ 8/parenrightbig/parenrightbigg ≈0.0369806 , ifi=7,...,18, 1 2−3 2πarccos/parenleftbigg (2−√ 3)(2−√ 6)+2 4√ 3−√ 3√ 3−√ 6/parenrightbigg ≈0.0351956 , ifi=19,...,26. The zonotope Zis the convex hull of all points of the form∑26 i=1εivi//bardblvi/bardbl, where (ε1,...,ε26)runs through all vectors in {−1,0,1}26. Thequickhull-algorithm (Barberet al., 1996) was used to find this convex hull.Zhas 96 vertices, inradius rVor≈0.463312, circumradius RVor≈0.511386, and thickness of minimal annulus TVor≈0.0480748. As γ3=1/2 we obtainthebounds 0.927≤lim tց0E/hatwideS3 26(A;t) S(A)≤1.023. Although the definition of the Bonnesen weights was restricted to the planar case, it can naturally 174 ImageAnalStereol2009;28:165-177 be generalized to higher dimensions. We will do so for comparison with the established Voronoi weights. Recall the construction for the Bonnesen weights in the plane for given vectorsv1,...,vk∈R2\ {o}. We definedui=vi//bardblvi/bardbl,i=1,...,k, and constructed the polygon/tildewidePinEq.15 with outer unit vectors ±u1,...,±ukcircumscribing the unit ball. We then chose/tildewideciproportionaltothelengthoftheedgeof /tildewidePwith outer unit normal ui. In higher dimensions, for given v1,...,vk∈Rd\{o}, we setui=vi//bardblvi/bardbl,i=1,...,k, and /tildewideP:=k/intersectiondisplay i=1/braceleftBig x∈Rd:|/a\}bracketle{tui,x/a\}bracketri}ht| ≤1/bracerightBig , incompleteanalogyto Eq.15.Hence/tildewidePisthepolytope circumscribing the unit ball with facet normals in {±u1,...,±uk}. We then choose /tildewideciproportional to the(d−1)-dimensional volume of the facet of /tildewideP which has outer unit normal ui. In the present three- dimensional example (with ζ1=ζ2=1and allgrid pointsinV(3) 2) the Bonnesen weights were computed withquickhull andare givenby /tildewideci≈  0.0465894 ,ifi=1,...,6, 0.0367439 ,ifi=7,...,18, 0.0349421 ,ifi=19,...26. The associated inradius is rBon≈0.462424, the circumradius is RBon≈0.511243, and thickness of minimal annulus is TBon≈0.0488187. This shows that Bonnesen weights are not better than Voronoi weights for the 3d−1 directions in dimension d=3. But this is not surprising because they are based on Bonnesen’s inequality Eq.14 which is valid only for two-dimensional convex bodies. In the last section we discuss how the approach developed in this paper couldbeextendedto higherdimensions. Finally,weshowhowtheasymptoticrelativeerror bounds in the case of a quadratic ( ζ1=1) planar grid depend on the choice of n. To do so, the symmetrized spread of the normalized vectors of V(2) nmust be determined. Lemma 6 For n ≥2the symmetrized spread of/braceleftBig x//bardblx/bardbl:x∈V(2) n/bracerightBig is equalto dn=/radicaltp/radicalvertex/radicalvertex/radicalvertex/radicalvertex/radicalbt2−/radicaltp/radicalvertex/radicalvertex/radicalvertex/radicalbt2 1+n−1/radicalBig 1+(n−1)2 ≤1 2(n−1). (30)ProofLetΔbe the symmetrized spread of D:=/braceleftBig x//bardblx/bardbl:x∈V(2) n/bracerightBig and setm:=n−1. Clearly, the arcC⊂S1in the first quadrant with endpoints (1,0)⊤,/parenleftbig m2+1/parenrightbig−1/2(m,1)⊤∈Ddoesnotcontainany other points in Dand thus the spread Δis at least the distance of the midpoint of Cto one of its endpoints. Hence Δ≥/radicalbigg 2/parenleftBig 1−cosϕ 2/parenrightBig =/radicaltp/radicalvertex/radicalvertex/radicalbt2/parenleftBigg 1−/radicalbigg 1+cosϕ 2/parenrightBigg =dn, whereϕ=arccos/parenleftBig m/√ 1+m2/parenrightBig is the length of C. This interpretation of dnalso shows the inequality in Eq.30, asdncannot be larger then half the distance of(1,0)⊤from (1,1/m)⊤. To show that Δ≤dn, letv,v′∈Dtwo points such that the sub-arc of S1 connecting them does not contain any other points ofD. Using reflections and translations leaving Z2 invariant, we may assume that v=x//bardblx/bardblandv′= x′//bardblx′/bardblwherex,x′∈V(2) nand the angles they form with the x-axis are at most π/4. We refer to Fig. 6 for a sketch of the situation. The cone between the rays spanned by xandx′cannot contain any other points ofV(2) nin its interior. Let y(y′) denote the point in V(2) n∩ {(m,t):t≥0}with largest (smallest) second coordinate below (above) this cone. Then the length of the arc in S1with endpoints vandv′is at most the length of the arc C′⊂S1with endpoints y//bardbly/bardbland y′//bardbly′/bardbl.Thesegment [y,y′]doesnotcontainanypoints ofV(2) n, soyandy′are distance one apart, and the length of C′is bounded from above by the length of C. Hereweusethefactthatamongallunitintervalsin/braceleftBig (m,t)⊤:t≥0/bracerightBig , the interval/bracketleftBig (m,0)⊤,(m,1)⊤/bracketrightBig has the largest gnomonicprojection. This gives Δ≤dn, as vandv′wherearbitrary in D. InviewofTheorem4,Lemma6yieldsanestimate for the asymptotic relative error using the Bonnesen- weightsci=/tildewideciwhenever n≥2. Theorem7 Let A ⊂R2be a topologically regular polyconvex set with positive perimeter, n ≥2, and ζ1=1. Let/hatwideS2 k(A;t)be defined by Eq.28, where {v1,...,vk}=V(2) nand ci=/tildewideci, i=1,...,k, are the Bonnesen weights. Then the asymptotic relative mean errorobeys lim tց0/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingleE/hatwideS2 k(A;t)−S(A) S(A)/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle≤π2 12(n−1)−2. 175 KIDERLEN MET AL:Errorboundsforsurfaceareaestimators Fig. 6.Construction to determine the symmetrized spread in Lemma 6: the points x and x′inV(2) n are contained in the cube [0,n−1]2and their normalizationsv ,v′inS1. EXTENSIONTO HIGHER DIMENSIONS Large parts of the present worst case analysis for quadrature rules, including the use of an associated zonotope Z, are not restricted to the two-dimensional setting. In order to find an easily accessible upper bound for the width of the minimal annulus, a joint extension of Bonnesen’s refined isoperimetric inequalities Eq.14and the geometric inequality of Groemer (1990) to higher dimensions (which is known) is not suitable, as it involves the surface area ofZ. Instead, a strengthened version of Uhrysohn’s inequalityis appropriate.Itreads /parenleftbiggb(Z) b(Bd)/parenrightbiggd −Vd(Z) Vd(Bd)≥cd(Z)δ2 2/parenleftBig Z,γdBd/parenrightBig .(31) Herecdis an explicitly known constant, depending ond, the mean width b(Z)ofZ, and the second intrinsic volume of Z. This inequality is a specialcase ofawholefamily ofgeometricinequalitiesderivedby Groemer and Schneider (1991), who also showed that it implies theBonnesentypeinequality /parenleftbiggb(Z) b(Bd)/parenrightbiggd −Vd(Z) Vd(Bd)≥c′ d(Z)(R−r)(d+3)/2.(32) The constant c′ d, again, depends on d,b(Z)and the second intrinsic volume of Z. Ford=2Ineq.32 is of the same form as Ineq.14, but with a weaker exponent. This apparently suboptimal exponent, and the problem to determine the zonotope with given mean width which minimizes the left hand side of Eq. 32ford>2,limits theusefulnessoftheseinequalities forapplications.ACKNOWLEDGMENTS We want to thank Luis Cruz-Orive, Maria Hernandez Cifre, Ilya Molchanov, and the two anonymous referees for their helpful comments on an earlier version of this paper. The first author’s work was partially supported by the Carlsberg foundation and the Danish Council for Strategic Research. The second author’s work was partially supported by a Marie Curie Fellowship of the European Community Programme. REFERENCES Barber C, Dobkin D, Huhdanpaa H (1996). The Quickhull algorithm for convex hulls. ACM Trans Math Soft 22:469–83. Bonnesen T (1929). Probl` emes des Isop´ erim` etres et des Is´ epiphanes.Paris: Gauthier-Villars. Gardner R, Kiderlen M, Milanfar P (2006). Convergence of algorithms for reconstructing convex bodies and directionalmeasures.AnnStat34:1331–74. GoodeyP, Weil W (1993).Zonoidsand generalizations.In: GruberP,WillsJ,eds.,HandbookofConvexGeometry, vol.B. Amsterdam:Elsevier,1297–326. Groemer H (1990). Stability properties of geometric inequalities.AmMathMon97:382–94. Groemer H, Schneider R (1991). Stability estimates for some geometric inequalities. Bull Lond Math Soc 23:67–74. GruberP(2004).Optimumquantizationanditsapplications . AdvMath186:456–97. Jan´ aˇ cek J (2001). Estimating length and surface area by systematic projections. Image Anal Stereol 20 (Suppl. 1):89–94. Kiderlen M, Jensen E (2003). Estimation of the directional measure of planar random sets by digitization. Adv ApplProbab35:583–602. Kiderlen M, Rataj J (2006). On infinitesimal increase of volumes of morphological transforms. Mathematika 53:103–27. Knebelman M (1941). Two isoperimetric problems. Am MathMon48:623–7. Lindel¨ of L (1869). Propri´ et´ es g´ en´ erales des poly` edr es qui, sous une ´ etendue superficielle donn´ ee, renferment le plusgrandvolume.BullAcadSciStPetersburg14:258– 69.MathAnn(1870)2:150-159. Moran P (1966). Measuring the length of a curve. Biometrika53:359–64. 176 ImageAnalStereol2009;28:165-177 Ohser J, M¨ ucklich F (2000). Statistical Analysis of Microstructures in Material Science. Chichester: John Wiley& Sons. Rataj J (2002). Determination of spherical area measure by meansofdilationvolumes.MathNachr235:143–62. Rother W, Z¨ ahle M (1990). A short proof of a principal kinematic formula and extensions. Trans Amer Math Soc321:547–58.SchneiderR(1993).ConvexBodies:TheBrunn-Minkowski Theory,vol.44ofEncyclopediaofMathematicsandits Applications.Cambridge:CambridgeUniversityPress. Schneider R, Weil W (1992). Integralgeometrie. Stuttgart: Teubner. Steinhaus H (1930). Praxis der Rektifikation und zum L¨ angenbegriff. Berich Sachs Akad Wiss Leipzig 82:120–30. 177