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