Image Anal Stereol 2010;29:45-52
Original Research Paper
VARIANCES OF LENGTH AND SURFACE AREA ESTIMATES BY SPATIAL GRIDS: PRELIMINARY STUDY
jiiil JanaAccek and Lucie KubInovaA
Institute of Physiology, Academy of Sciences of the Czech Republic, Videiiska 1083, CZ-14220 Praha, Czech Republic
e-mail: janacek@biomed.cas.cz, kubinova@biomed.cas.cz
(Accepted December 4, 2009)
ABSTRACT
Periodic spatial grids can be used for unbiased estimation of length and surface area of objects by counting or measuring intersections of the objects with the grids. The estimators are theoretically based on discrete approximation of well established integral geometric formulas. The variance of the estimates depends on properties of both the grid and the measured objects. Main results of the theory of variance of the isotropic uniform random (lUR) volume estimation by spatial grids, especially a formula relating the variance of the volume estimator with the object surface area and the grid constant, are recapitulated. To identify main features of length and surface area lUR estimates the variance due to rotation and simple asymptotic formulas for the residual variance of estimates of selected model objects is calculated. Surface area estimates by multiple grids of parallel lines in 3D and of the variance of length estimates by periodic grids of planes or spheres in d -dimensional space are studied.
Keywords: stereology, variance, 3D.
INTRODUCTION
Basic geometrical characteristics of 2D or 3D objects can be estimated by counting or measuring intersections of the objects with randomized periodic patterns of manifolds (Barbier, 1860; Santal(3, 1976; Cruz-Orive, 1997). Estimation of geometrical characteristics by measuring the objects intersections with the randomized probes can be regarded as a generalization of the classical topic of stereology - estimation of the characteristics measuring lower dimensional randomized sections of the objects. The formulas for unbiased estimators using intersections with randomized spatial grids follow from well established general formulas of integral geometry (Baddeley and Vedel-Jensen, 2005). On the other hand the variance of the estimates may depend on properties of both the grid and the measured objects in a complicated way.
The variance of volume estimation by measuring intersection of the object with arbitrary periodic grid in random position is already well understood. A useful approximate formula (Eq. 1) for the variance using only simple characteristic of the measured object is available. The formula for variance of d—dimensional volume estimate by point lattice was proved in (Hlawka, 1950) for strictly convex sets with 6 d—times differentiable support function, the proof for convex sets or sets with C1-5 smooth boundary can be found in (Janacek, 2008) and it can be easily generalized to an arbitrary periodic grid. The volume of a bounded
object K in d—dimensional space can be estimated by a point counting method using regular point lattice T = AZd in Rd, where Z is the set of integers and A is a regular matrix. The number of points of the randomly shifted point lattice inside the object multiplied by the specific volume of the point is an unbiased estimate of the volume. Under mild regularity conditions on the object K the variance of the estimator with the randomly rotated lattice scaled by the factor u > 0 is
CT^d—1 (dK) F (u—1) ud+1
(Janacek, 2008), where CT is a constant of the lattice, Hd—1 (dK) is the surface measure of the object, F > 0 is a function which ergodic average tends to 1:
1 r
li^ - F (t) dt = 1 -
xj 0
Th(e consta)nt of the lattice can be expressed as
Ct = (2n2dKd)—'|X|—d—1 where T* = A-1Zd
and Kd is the volume of the unit ball in Rd. The value of the coefficient is approximately 0.072 837 for the quadratic lattice and 0.071701 for the hexagonal lattice normalized to unit intensity in plane, 0.066649 for the cubic lattice, 0.064350 for the normalized body centered cubic lattice and 0.064390 for the normalized face centered cubic lattice in space.
The result can be easily generalized from the counting measure of lattice T to an arbitrary T— periodic measure m with constant of the measure C^ =
{2K2dKd
)—1 yX=0 I Z^T*
MX ^ IX l d 1 where jn^ are Fourier
coefficients of the measure (Janaček, 2006). Software package pgs for calculation of coefficients of various grids is available (Kieu and Mora, 2006).
For the purpose of the variance estimation the function F oscillating around 1 is neglected and the variance can be expressed approximately:
Var (estV) CsHd-1 (dK) NV
d+1 " d
(1)
where Nv = |^|-1 is the grid points intensity and the lattice S is the homothetic copy of the point lattice T with unit intensity. Similar simple approximations for estimators of length or surface area are not yet known. The aim of this paper is to study the behavior of variance of selected estimators of length or surface area, design of the expressions for the estimators variance similar to Eq. 1 and evaluation of the coefficients in the expressions.
The length of a fibre-like structure in Rd can be measured with a randomly rotated and shifted (lUR) periodic grid of surfaces. The surface area of objects in Rd can be estimated by counting intersections with a grid of lines. The variances of the estimators can be decomposed into the variance due to the orientation of the grid and to the residual component of variance:
Var (est Q) =
VargeSOd (E (est Q|g)) + Eg^SOd (Var (est Q|g)) ,
where Q is either length (L) or surface area (S) and SOd is the group of rotations in Rd.
The variances of length or surface area estimates in Rd using spatial grids due to orientation are calculated in the same way as the variance of the surface area estimators using combined projection, treated in the following section.
VARIANCE OF SURFACE AREA ESTIMATE IN Rd USING TOTAL PROJECTION
Crofton formula relates the suface area of a convex body to the mean area of its isotropic projection. Surface area of a convex body in Rd is average area of its projection multiplied by dKd/Kd-1 where Kd is the volume of the unit ball in Rd and using this equality we can construct an estimator of the surface area from the object total projection area (area weighted by Euler characteristics of projection inverse, see Serra, 1982):
est S = -^^^Atot kd-1
Variance of the projection area depends on the object anisotropy, i.e., on the distribution of the angles between the normals to the body. It is well known that the variance of the surface area estimators can be decreased by using average of areas of projections into a set of directions (with fixed mutual position, e.g., perpendicular) instead of one direction. The decrease of the variance is due to a negative covariances of the projections areas (Mattfeldt etal., 1985).
The formula for coefficient of variation (CV) of surface estimator from multiple weighted projection areas of a body (Janacek, 1999) is:
CV2 =
0 J0
Kd(y,c) dF(y)dG(c) , (2)
where F is the distribution function of angles between normals to the object surface, G is the distribution function of angles between projections directions and the kernel Kd is defined as
Kd (y, c) =
dkd
\2kd-u J SOd
l(gy, x)(gu,v)l dg-1
Kd(y,c) can be evaluated either by direct integration or using convolution formula and Parseval equality for harmonic analysis on Sd-1. For details and examples see Appendix A.
LENGTH ESTIMATION BY THE GRID OF SURFACES
The length of a curve in Rd can be estimated from the number of intersection of the curve with a spatial grid ofplanes
est L =
dkd N 2Kd-1 Sv
where N is the number of intersections and SV is the surface area of the grid per unit volume. The variance of the estimator of the length due to orientation can be calculated by Eq. 2 substituting for F the distribution function of angles between the tangents to the curve and for G the distribution function of angles between the normals to the grid surfaces. The simplest of grids of surfaces is the grid composed of parallel planes. The grids with different directions of normals can be combined obtaining decreased variance due to orientation. The variance is decreased due to negative covariances as in the case of surface estimation by projections. The residual component of variance, EgeSOd (Var (est L|g)), is the mean over all orientations of variance of an estimate by equidistant point samples
of an integral of total projection of the space curve to a line of given orientation. The total projection of the space curve to a line is a step function and the variance of integral estimate of such step function is well known (Moran, 1950). An approximation of the residual component of variance of estimate of length L of the object K composed of space curves by applying grids of parallel hyperplanes in then follows from differential geometry:
EgesOd (Var (estL|g)) - CLkabs(K) S-2 ,
kABS (K)= |k|d5 + P (Nends + Noddbranchings) ,
JK 4 rL W dkd ^
CG =
points (Fig. 1), and so the grid of touching spheres shall be effective for the length estimation in 3D.
3^ 2kd_ 1
where Sv is the grid surface intensity, k is curvature of the curve (i.e., the reciprocal of the radius of the osculating circle), Nends is number of endings of curves, Aoddbranchings is number of branchings where an odd number of curves is joining together and kd is the volume of the unit ball in Rd.
Another grid suitable for length estimation is composed of identical spheres centered in points forming periodic lattice T. The grid is for given lattice type defined by two parameters, the point lattice density and by the spheres diameter. Let us define shape parameter p of the grid equal to the radius of the spheres of the grid homothetically transformed in such a way, that the point lattice T is transformed to the unit lattice S. If the grid is randomly rotated the variance due to orientation of the grid, Varg£SOd (E (est QIg)), is obviously zero.
The residual component of variance of the estimate of length of a line segment can be approximated by the variance of the estimate of volume of a cylinder with a diameter equal to the diameter of the spheres in the original grid by a point lattice T. The covariogram of the cylinder can be approximated by a direct product of covariogram of (d — 1) — dimensional ball and a line segment. An approximate expression for the residual component of variance of estimate of length L of the object K by the grid of spheres in Rd is:
Eg,,SOd (Var (estL|g)) - C'g(p)L(K) S-1 ,
CL = CG(p) =
d — 1 kd-1
X=0
d/2—1 (2nXp) ,
where J is the Bessel function of the first kind, follows from the limit transition Sv ^ ¥.
The function CGg((p) has a local minimum close to the half of the distance between the nearest lattice
cL p
Fig. 1. Values of CG(p) for grid of spheres S with centres in face centered cubic lattice.
In estimates of the curve length the grids of flat planes yield quite different order of asymptotics than grids of curved surfaces, due to increasing curvatures during homothetic downscaling of the curved grids. The variance is also controlled by different properties of the measured object, namely by the total absolute curvature in estimations using parallel surfaces and by the length in estimation by application of spheres.
SURFACE AREA ESTIMATION BY GRIDS OF STRAIGHT LINES
The surface area of an object in Rd can be estimated from the number of intersection of the object surface with a spatial grid of lines
est S =
dkd N 2kd-1 Lv
where N is the number of intersections and Lv is the length of the grid per unit volume. The grids composed from sets of parallel lines can be used for estimation of the surface area of spatial objects. The estimators can be optimized by finding optimal mutual position of the sets (Kubinova and Janacek, 1998). Parametric expressions of two simple grids and three multiple grids in basic and optimized shifted versions (Fig. 2), are listed below, where i, j G Z are discrete parameters, u G R is a continuous parameter.
3
0
p
S»
Fig. 2. Spat^ial gr^ids of lines in R3 used for sur^face ar^ea est^imat^ion. (a) g^ri^d w^th quadr^at^ic cross-sect^ion; (b) gr^d w^it^h t^r^iang^ular cross-sect^ion; (c) 3-fold basic g^r^id; (d) 3-^old shii^ted g^r^id; (e) 4-^old basic g^r^id; (f) 4-^old shi^^ed ^r^id; (g) 7-^old basic ^r^id; (h^) 7-^old shi^^ed g^r^id.
Unit grid with quadratic cross-section: (i, j, u)
Grid with triangular cross-section: (i + u, j + u, u)
Unit grid with triangular cross-section - the above grid scaled by the factor
•	3-fold basic grid: (u, i, j), (i, u, j), (i, j, u),
•	3-fold shifted grid: (u, i, j + 2), (i + 2, u, j),
(i, j + 1. u .
•	Unit 3-fold grids - the above grids scaled by the
factor \/3.
•	4-fold basic grid: (i + u, j + u, u), (i + u, j + u, -u)
(i + u, j - u, u), (i + u, j - u, -u).
•	4-fold shifted grid: (i + u, j + u, u),
{i + u, j + u + 2, -u), (i + u + 1, j - u, ^ , (i + u, j - u,-u + 2).
•	Unit 4-fold grids - the above grids scaled by the factor 2Ay3.
•	7-fold grids are combinations of above 4-fold grids with 1 homothetic images of 3-fold grids.
•	7-fold ba)ic grid: (u, 2i, 2 j, (2 i, u, 2 j,
(2i, 1 j, u), (i +u, j +u, u), (i +u, j +u,-u),
(i + u, j - u, u), (i + u, j - u, -u).
•	7-fold shifted grid: {u, 1 i, 1 j + 4), {2 i + 4, u, 2 j),
{1 i, 2 j + 4, ^, (i' + u, u),
{i + u, j + u + 2, -u), (i + u + 2, j - u, ^, (i + u, j - u,-u + 2).
•	Unit 7-fold grids - the above grids scaled by the factor ^ V3 + V3.
The operation of estimation of the area of a surface by the grid of parallel straight lines can be decomposed into parallel projection to the plane perpendicular to the lines, and estimation of the total projection area by a point lattice. The variance of the estimator can be decomposed into the variance due to rotation of the grid and to the residual component of variance.
THE VARIANCE OF SURFACE AREA ESTIMATES DUE TO ROTATION.
Coefficients of variation of estimate with the grids of lines due to orientation S-^Vargesod (E (est S|g)) can be calculated by Eq. 2 substituting for F the distribution function of angles between normals to surface and for G the distribution function of angles between tangents to grid lines. The variances due to orientation for a flat object, that are also upper bounds for the variances of the surface area estimates for arbitrary objects, are shown in Table 1. and the results for a cylindrical surface are presented in Table 2. The gain in effectivity with respect to independent measurements with a simple grid is also shown.
Table 1. Coefficient of variation of projection of a flat object into n directions offene gr^ids in R3 calculated using Eq. 4 from Appendix A.
n	CV	Cy2/nCy2
1	0.57735	1
3	0.10163	10.8
4	0.07523	14.7
7	0.03822	32.6
Table 2. Coefficient of variation of projection of a cy^lindrical surface into n directions of line grids in R3 calculated using Eq. 5 from Appendix A.
n	CV	Cy2/nCy2
1	0.28418	1
3	0.03706	19.6
4	0.02668	28.4
7	0.01196	80.6
Table 3. Coefficients CG (K) in Eq. 3 for residual component of variance of estimates of surface area w^ith unit (Ly = 1) line gr^ids G in R3; namely sample (1) grids with square point lattice cross-section and w^ith triangular point lattice cross-section (K is ar^bitr^ary) and multiple (3, 4, 7) basic and shifted gr^ids (K is ball or ^isc).
n	G	CG (K)	G	C5 (K)
1	square	3.661	triangular	3.604
3	basic, ball	9.717	basic, disc	8.490
3	shifted, ball	3.810	shifted, disc	4.730
4	basic, ball	11.920	basic, disc	10.376
4	shifted, ball	3.675	shifted, disc	4.832
7	basic, ball	16.594	basic, disc	14.265
7	shifted, ball	4.022	shifted, disc	5.769
THE RESIDUAL COMPONENT OF VARIANCE OF SURFACE AREA ESTIMATES.
The residual component of variance of surface area (5) of object K estimate by grids of lines in R3 can be expressed approximately similarly to Eq. 1 as
Ege503 (Var (est5|g)) = CG{K)HAB5(K) L-2 , (3)
where Ly is the grid length intensity, HAB5 (K) is the absolute width of K defined using the absolute mean curvature by equality 2nHAB5 = Kf (Baddeley, 1980). The values of the coefficients in Eq. 3 for our grids (Table 3) and for ball and disc are calculated in Appendix B.
The coefficients of simple grids are independent of the body. In the multiple grids it is not true as demonstrated by the examples of the estimates of the disc and ball surface areas.
DISCUSSION
The study of variances of estimators with periodic grids is important for the design of efficient measurement methods. Approximate asymptotic results for homothetic images of grids with spatial density increasing to infinity, as those applied in the paper to variance of volume estimates or to residual components of variances of the length and surface area estimators, use simple properties of the measured objects and yield useful results for model objects under study. We suppose that the formulas are generally valid, and the coefficients of ball and disc in formula for the residual component of variance of surface area estimates represent extreme values of coefficients for arbitrary objects. Using the knowledge on the behavior of the estimators we were able to design efficient estimators of surface area using spatial grids of shifted line probes and estimators of length of fibre like objects using spatial grids of touching spheres. The asymptotic results are established rigorously in the case of volume estimates only, the theory of variance in more difficult cases - the surface area and length estimates - is far from completness yet, however the presented study with model objects exhibits some important features of the variance of length and surface area estimators.
ACKNOWLEDGEMENT
This study was supported by the Academy of Sciences of the Czech Republic, grants No. A100110502 and AV0Z 50110509.
REFERENCES
Baddeley AJ (1980). Absolute curvatures in integral geometry. MathProc Camb Phil Soc 88: 45-58.
Baddeley A, Vedel-Jensen EB (2005). Stereology for statisticians. Boca Raton: Chapman & Hall/CRC.
Barbier JE (1860). Note sur probleme de l'aiguille et le jeu du joint couvert. J Math Pure Appl 5:273-87.
Cruz-Orive LM (1997). Stereology of single objects. J Microsc 186:93-107.
Hlawka E (1950). Über Integrale auf konvexen Körpern I. Monatsh Math 54:1-36.
Janacek J (1999). Errors of spatial grids estimators of volume and surface area. Acta Stereol 18:389-96.
Janjiček J (2006). Variance of periodic measure of bounded set with random position. Comment Math Univ Carolinae 47:473-82.
Jan^Cek J (2008). An asymptotics of variance of the lattice points count. CzechMath J 58:751-75.
Kieu K, Mora M (2006). Precision of stereological planar area predictors. J Microsc 222:201-11.
Kubinova L, Jan^cek J (1998). Estimating surface area by the isotropic fakir method from thick slices cut in an arbitrary direction. JMicrosc 191:201-11.
Matheron G (1965). Les variables regionalisees et leur estimations. Paris: Masson et Cie.
Mattfeldt T, Mobius H-J, Mall G (1985). Orthogonal triplet probes: an efficient method for unbiased estimation of length and surface of objects with unknown orientation in space. J Microsc 139:279-89.
Moran PAP (1950). Numerical integration by systematic sampling. Proc Camb Phil Soc 46:111-5.
Santalo LA (1976). Integral Geometry and Geometric Probability. Reading: Addison-Wesley.
Serra J (1982). Image Analysis and Mathematical Morphology. London: Academic Press.
APPENDIX A
Let us start with auxiliary computations. Let x = x,),a = («1,...a,), = l^!1...|^dd, |a| = Sd^jtti^then
S^ 1
l^a dS (x) =
nd=1 r
ai+1
2 rf
where dS is surface measure. The equality can be proved by calculating fRd |^a exp	dx in two
ways: by Fubini theorem and in spherical coordinates.
The indentity f^ ^ 1 dS(x) = dk, can serve as an
example where Kd = p 2 ^ 2 + ^ 1 is the volume of the unit ball B, (1) in Rd.
As f^ ^ l^il dS(x) = 2kd-1, the mean of the projection length of a unit segment to a line (a simple projection) is 2kd-1/dkd. As f^ ^ |^1|2dS(x) = Kd, the squared coefficient of variation of the projection is
f dkd \ '1 \2kd-1j
CV2 =
1
The values of the the CV2 for d = 1,2,3,4,... are 0,p2/8 - 1,1/3,9p2/64 - 1,... = 0,0.23,0.33,0.39... (the limit is p/2 - 1= 0.57).
EgeSO^d (L1L2), the mixed second order moment of lengths of projections of a unit segment to two lines spanning angle y in R2 (a double projection) is P/p=0sin — y| sin f df =
p (sin y+ (p - y cos y. In corresponding value
Rd we obtain the
p
, , sin y+ — - y cos y dp 2
(4)
multiplying the value for R2 by the value fSd-1 (x? + dS(x) = 2/d. EgeSOd (L1L2), the mixed second order moment of lengths of projections of unit segments spanning angle c to corresponding lines spanning angle y is
|(g7,x)(^u, v)l dg,
SOd
where x, j, u, v G Rd, l^i = = | u| = | ^ = 1, (y, u) = cos y, (x, v) = cosc, SOd is the group of rotations equipped with the probabilistic invariant measure. From this formula the kernel Kd can be calculated. Special values of Kd can be calculated from preceding formulas. Thus CV2 of the simple projection of unit segment gives Kd (0,0) and CV2 of the double projection of a unit segment gives Kd (y, 0). It is easy to see that K2 (y,c) = K2 (y- c, 0).
By Parseval equality for harmonic analysis on Sd-1 we obtain the following expression:
2
pJy=0 S¥=0 (4n + 1)
k3 (y, c) dy =
dkd \2
2kd-1y (2n - 3)!!
/(9n- 3)!^ ^(2n- 1)!!
2n!!
vm
m
V(2n + 2)!!; V (2n-2l- 1)!! (2n + 2l- 1)!! (2n + 2l)!! (2n - 2l)!!
cos (2lc) - 1 ,
(5)
which enables us to calculate the variance of multiple projections of a cylindrical surface in R3. It was obtained from Fourier series in spherical harmonics of absolute value of sine of lattitude, of 1-dimensional measure supported by equator and of sum of Dirac measures located c radians apart from each other.
APPENDIX B
The variance of estimates of surface area of a convex body by periodic grid oflines can be calculated using generalization of isotropic (point) covariogram of the body (Matheron, 1965): the value of generalized covariogram for pair of lines is the measure of set
2
x
2
x
of such euclidean motions of the lines, that the both lines intersect the body, multiplied by the sine of angle between the lines, or the limit for parallel lines. The value for a pair of parallel lines is equal to mean isotropic point covariogram of the planar projection of the body. For nonparallel lines we obtained so far a simple result for ball or disc only. For ball we obtained montee of covariogram of the circle (Matheron, 1965). For disc we obtained montee of mean covariogram of the disc planar projection multiplied by (sin y + (Pp - y cos y where y is angle between the lines. Using Poisson summation formula yields the following formulas for coefficients in Eq. 3 of the ball (as in Janacek, 1999) and the disc (Table 3).
The Fourier coefficient of a aZ3 - periodic measure m with index X G a-1Z3, a > 0, is
mx =
3
CS (ball) = -2
p2
{0,a)3 X=0
I
= a-1Z3
exp(-2nixX)dm(x) , (6)
mx I2 IXI
3
Xea
Fourier coeficients of unit grid with quadratic cross-section are ]lx = 1 for X = (i, j, 0), i, j G Z, mix = 0 otherwise.
Grid with triangular cross-section:
•	-mx = for X = (i, j,-i- j), i, j G Z
•	]xx = 0 otherwise.
Coefficients of unit grid with triangular cross-section are multiples of above coefficients with indices x
scaled by factor . 3-fold basic grid:
•	mx = 1 for x = (1, m, 0), x = (m, 0,1) and x =
(0,1, m).
•	mx = 2 for x = (1,0,0), x = (0,1,0) and x = (0,0,1).
•	Jix = 3 for x = (0,0,0).
•	Jix = 0 otherwise. 3-fold shifted grid:
•	mx = (-1)m for x = (0,1, m), x = (m, 0,1) and x = (1, m, 0).
•	mx = 2 for x = (21,0,0), x = (0,21,0) and x = (0,0,21).
•	mx = 3 for x = (0,0,0).
•	mix = 0 otherwise.
Coefficients of unit 3-fold grids are ^ multiples of above coefficients with indices x scaled by factor --=.
4-fold basic grid:
•	mix ^ \/3 for x = (1, m, n) where 1, m, n G Z \ {0}, 1 = ±m = ±n = 1, 1 + m + n = 0, 1 + m - n = 0, 1 - m + n = 0or 1 - m - n = 0.
•	Jlx = 2-3 for x = (1, -1,0), x = (1,0, -1), x = (0,1, -1), x = (0,1,1), x = (1,0,1) and x = (1,1,0).
•	mix = 4-3 for x = (0,0,0). 4-fold shifted grid:
•	Ji^x ^ ^^ for x = (1, m, n) where 1, m, n G Z \ {0}, 1 = ±m = ±n = 1, 1 + m + n = 0.
•	Jix =(-1)^—3 for x = (1, m, n) where 1+m - n = 0.
•	Jlx = (-1)^—3 for x = (1,m,n) where 1 - m+n = 0.
•	Jix = (-1)^—3 for x = (1, m, n) where 1 - m - n = 0.
•	mx = 2-3 for x = (21, -21,0), x = (21,0, -21), x = (0,21, -21), x = (0,21,21), x = (21,0,21) and
x = (21,21, 0).
•	mx = 4-3 for x = (0,0,0).
Coefficients of unit 4-fold grids are multiples of above coefficients with indices x scaled by factor
7-fold basic grid:
•	jix = 4 for x = (21,2m,0), x = (2m, 0,21) and x = (0,21,2m) where 1, m G Z \ {0}, 1 = m.
•	Jix = 8 for x = (21,0,0), x = (0,21,0) and x = (0,0,21).
•	Jix ^ —3 for x = (1, m, n) where 1, m, n G Z \ {0}, 1 = ±m = ±n = 1 and 1+m + n = 0,1+m - n = 0, 1 - m + n = 0or 1 - m - n = 0.
•	Jix = 4 + 2-3 for x = (21, -21, 0), x = (21, 0, -21), x = (0,21, -21), x = (0,21,21), x = (21,0,21) and
x = (21,21, 0).
•	mix = 2-3 for x = (2i + 1, -2i + 1,0), x =
(2i + 1,0, -2i + 1), x = (0,2i + 1, -2i + 1), x =
(0,2i + 1,2i + 1), x = (2i + 1,0,2i + 1) and x = (2i + 1,2i + 1,0) where i G Z.
•	Jix = 12 + 4-3 for x =(0,0,0). 7-fold shifted grid:
•	= (-l)m4 for X = (2l,2m,0), X = (2m,0,2l) and X = (0,2l, 2m) where l, m G Z \ {0}, l = m.
•	= 8 for X = (4l, 0,0), X = (0,4l, 0) and X = (0,0,4l).
•	ilX = for X = (l, m, n) where l, m, n G Z \ {0}, l = ±m = ±n = l, l+m + n = 0.
•	mi^X = (-l)^V^for X = (l, m, n) where l+m - n = 0.
•	JlX = (-l)^V^for X = (l,m,n) where l — m + n = 0.
•	mi-X = (—X = (l,m,n) where l—m — n = 0.
•	ßX = (—l)l4 + 2^3 for X = (2l, —2l, 0), X =
(2l, 0, —2l), X = (0,2l, —2l), X = (0,2l,2l), X =
(2l, 0,2l) and X =(21, 2l, 0).
•	mX = 12 + 4^3 for X = (0,0,0).
Coefficients of unit 7-fold grids are (4 (3 ^ v^^ multiples of above coefficients with
indices X scaled by factor ^2 V3+V3
4 ■ 34 Z(^2,3) + 6Z(I2,3) +
We can see that the shifting caused vanishing of the coefficients with smallest indices, which explains the superior performance of the shifted grids (Table 3).
Properly grouping the terms in expression for the coefficients we can use the Epstein zeta functions Z (A, 3) defined as
n=0
Z (A, s)= X |An|—s ,
neZd
with matrix A equal to Id , d = 1 or 2, the identity matrix in Rd and
A V2 ( ./I 0 \
A2=2 l)
The values of coefficients Cq(ball) are
•	z(I2,3) for simple grid with quadratic crossection,
•	pj Z (A2, 3) for simple grid with triangular crossection,
•	(3z(I2,3) + 6Z(I1,3)) for basic 3-fold grid and ball,
•	(3z(I2,3) — 2Z(Ii,3)) for shifted 3-fold grid and ball,
•	p42i3i (4■ 3—iZ(A2,3) + 12■ 2—3Z(Ii,3^ for basic 4-fold grid and ball,
•	■p2^23'4 (4 ■ 3—i Z(A2,3) — 9 ■ 2—2 Z (Ii, 3^ for shifted 4-fold grid and ball,
p2 2V3W3
+ 3 (4 + ^ V2	Z (Ii, 3^ for basic sevenfold
grid and ball,
4 ■ 34 Z (A2,3) + 6Z (I2,3) —
p2 2V^V3
— 9 (4 + ^ V2 ^V^ Z (Ii, 3^ for shifted sevenfold grid and ball.
The Epstein zeta functions Z(I1,3), Z(I2,3) and Z(A2, 3) can be calculated using the following identities valid between the Epstein zeta functions, the Riemann zeta function Z and Dirichlet function Lp
n
Z (s)= X n—s , Lp (s) = X (p|n)
n=1	n=0
where (p|n) is the Kronecker symbol:
Z(I1, s) = 2C (s),
ss
Z (I2, s)= 2)L-4( 2 Z (A2, s) = 21+2 31-4 ^ |)L—
The covariances of the simple estimates of the disc surface area must be mult)plied t)y additional factor f (y) = p (siny+ (p — y cosy where y is angle spanned by tangents to parallel lines of the simple subgrids.
Thus we obtain the values of coefficients CQ(disc):
• p4^^ (3Z(I2,3) + f (f) 6Z(I1,3)) for basic 3-fold grid and disc,
(3Z(I2,3) — f (f) 9Z(I1,3)) for shifted 3-
fc)ld grid and disc,
Li3if4 ■ 3—3 Z (A2,3) +
A1 f22
+ f (arccosl) ■ 12 ■ 2—2 Z(I1,3 for basic 4-fold grid and disc.
442i3i (4■ 3—iZ(A2,3) —
— f (arccosl) ■ 9 ■ 2—3Z(I1,3) for shifted 4-fold grid and disc.
2731+73 )4 ■ Z 3)+ 3)+(
+ 3 (4f (f) + 372f (arccos3) + 76f (arccos713 Z(I1, 3)) for basic sevenfold grid and disc.
f22 (S^^
4 ■ 34 Z (A2' 3)+ 6Z (/2,3) —
4 (4f(p) + f (arccosD ^76f (arccos7li Z(I1, 3)) for shifted sevenfold grid and disc.
4