FRACTIONAL TREND OF THE VARIANCE IN CAVALIERI SAMPLING

Cavalieri sampling is often used to estimate the volume of an object with systematic sections a constant distance T apart. The variance of the corresponding estimator can be expressed as the sum of the extension term (which gives the overall trend of the variance and is used to estimate it), the ’Zitterbewegung’ (which oscillates about zero) and higher order terms. The extension term is of order T for small T, where m is the order of the first non-continuous derivative of the measurement function f, (namely of the area function if the target is the volume). A key condition is that the jumps of the mth derivative f (m) of f are finite. When this is not the case, then the variance exhibits a fractional trend, and the current theory fails. Indeed, in practice the mentioned trend is often of order T, typically with 0 <q <1. We obtain a general representation of the variance, and thereby of the extension term, by means of a new Euler-MacLaurin formula involving fractional derivatives of f. We also present a new and general estimator of the variance, see Eq. 26a, b, and apply it to real data (white matter of a human brain).


INTRODUCTION
Cavalieri sampling is widely used in design stereology, notably to estimate the volume of a bounded object from systematic sections.To predict the variance of the corresponding estimator is therefore a relevant problem.In spite of recent, important advances, the current theory is failing to suit many real situations (Roberts et al., 1994;McNulty et al., 2000;García-Fiñana and Cruz-Orive, 1999).We generalize the theory and present new variance predictors which improve on the existing ones.
The target parameter is often the volume Q of a bounded connected and non-random subset X ⊂ ! 3 , which may be expressed as follows, where f (x) denotes the area function, namely the area of the intersection between X and a plane normal to a convenient axis at the point of abscissa x.In general f : !→ !+ is a bounded non-random function, called the measurement function, which is integrable on a finite support [a, b] ⊂ ! and vanishes outside [a, b].Recent theory (Souchet, 1995;Kiêu, 1997;Kiêu et al., 1999;Gundersen et al., 1999;García-Fiñana and Cruz-Orive, 1998;1999), applies when f is (m, p)-smooth, a concept we recall briefly.Given a function h : !→ !+ , the amplitude of the jump of the kth derivative h (k) of h at the abscissa is expressed as follows, ( ) ( ) ( ) ( ) ( ) ( ) Sh x h y h y , x , k , ,... , provided that the limits exist.Further, is the set of points where h (k) is non-continuous; in this paper we extend Dh (k) to include the points where Sh (k)  may not exist.We say that h is (m, p)-smooth, with smoothness constant m = 0,1,... and p = 1,2,..., both known and fixed, if the support of h is bounded and if a. Dh (k) = ∅, (k = 0,1,..., m -1), and m is the order of the first non-continuous derivative of h, and b. h (k) , (k = m, m + 1,..., m+p), may have only a finite number of jumps, and they have to be finite.
An unbiased estimator of Q is: where x 1 is a uniform random variable in [0,T ), T the distance between consecutive planes and f 1 , f 2 , ..., f n the observations of f at the sampling points which lie in [a, b].The mean value of n is (b -a)/T, (denoted also by n throughout, for short).
Suppose that f is (m, p ≥ 1)-smooth.Then, where ) the extension term and it explains the overall trend of the variance; Z(T) is the 'Zitterbewegung' term and it is made of functions which oscillate about zero.The extension term is proportional to even powers of T (i.e., to even powers of 1/n), and it can be expressed as follows, (Souchet, 1995;Kiêu, 1997), where P 2l (0) = B 2l (0)/(2l)!and B 2l (0) is a Bernoulli number, (Abramowitz and Stegun, 1965).
The preceding results, however, are valid if f is (m, p ≥ 1)-smooth, where m, p are integers and the conditions (a), (b) above hold.In the following section we give examples of objects whose area function violates condition (b).In Section 3 we present a new variance representation, more general than (5), which allows for a fractional smoothness constant.The general validity of the approach is examined in Section 4. In Section 5 the new variance predictors are presented and illustrated with real data.For reasons of space, the mathematical proofs are only sketched; complete proofs will be given in a forthcoming paper.

COUNTEREXAMPLES TO THE CURRENT THEORY
Figs. 1b, c represent two solids of revolution generated by a curve of equation |x| a + |y| a = 1, (a > 0), rotating around the (horizontal) OX axis, for a = 3/2, 8 respectively.When a = 2 the solid is a sphere, and when a > 2, a 'superegg', (Gardner, 1978).The solid in Fig. 1d is a right circular cylinder with its axis parallel to the OY axis.
Suppose that the objects are swept by a plane normal to the OX axis, see Fig. 1a.For the solids in Figs.1b,  c the corresponding measurement function is (up to a factor π) of the form, see Fig. 1e.This measurement function is continuous, that is Df (0) = ∅.However: Thus, the current theory can be applied only when a = 1, 2, because then the jumps Sf (1) (c) are all finite.
For the right circular cylinder (Fig. 1d) the measurement function is (up to a constant factor), see Fig. 1e, red curve.This measurement function is continuous, that is Df (0) = ∅.However, Sf (1) (±1) = +∞, and therefore the current theory cannot be applied either.The wavy lines in Fig. 1f -d.The broken lines represent the corresponding extension terms (Eq.17 below), whose slopes are seen to be fractional.Thus, none of these trends can be explained by the extension term ( 6), and a new approach is needed.

CHARACTERIZATION OF THE MEASUREMENT FUNCTION
We say that the function f is q-smooth, where q ≥ 0, if its support is bounded and if: a. Df (k) = Ø, (k = 0, 1, ..., q -1), where q is the smallest integer ≥ q. b.
We have that q is precisely q α − , with , and q is called the fractional smoothness constant of f.Note that if all 0 i α ± = , then the jumps of ( ) q f at the points {c i } are all finite, and we recover the integral smoothness case (Section 1) with p = 0.
The Riemann-Liouville definition of the q-th fractional derivative of f (x), (a (Nishimoto, 1984), where (r, s) = (a, x) for the right, and (r, s) = (x, b) for the left derivative, respectively; q is the largest integer ≤ q and Γ(•) is the Gamma function.With this definition the constants i d ± in the right hand side of Eq. 9 can be connected with the fractional derivatives of f at the corresponding discontinuity points as follows, ( ) ( ) ( )

GENERAL EULER-MACLAURIN FORMULA
The estimator Q is a periodic function of x 1 with period T, and it therefore admits the Fourier expansion ( ) ( ) Bearing in mind that k k a a − = we obtain ( ) ( ) ( ) If is q-smooth, from the preceding expression we obtain a new, generalized version of the Euler-MacLaurin formula, namely ( ) where is a generalized Bernoulli polynomial (Matheron, 1965).When q ∈ #, Eq. 14 reduces to the refined Euler-MacLaurin formula of Kiêu (1997, p. 36).

THE NEW VARIANCE REPRESENTATION
If is q-smooth, then taking the expectation on the square of the right hand side of Eq. 14 with respect to x 1 we obtain the general expansion The extension term has the following expression where Sf (q) is the 'fractional jump' of f (q) , defined as follows To obtain Eq. 17 we have used the following identities, , whereby Eq. 17 reduces to Eq. 6.

GENERAL VALIDITY
The special measurement functions ( 7) and ( 8) arise when the solids in Figs.1b, c and in Fig. 1d, respectively, are cut at one particular orientation (Fig. 1a).One might therefore wonder whether, for any object, fractional theory will at best apply for isolated orientations of total measure zero.The following example shows that this does not need to be the case.Var Q O T = , (Fig. 1).For 0 < θ < π/2, f and f (1) are continuous, whereas f (2) is not.It would therefore seem that the smoothness constant would be q = m = 2, but this is in fact not the case because f (2) is unbounded at its points of discontinuity (Fig. 2): the fractional theory is therefore needed.When T exceeds the distance between points of discontinuity, see Figs. 3a, d, terms should be recruited from the Zitterbewegung into the right hand side of Eq. 17 to improve the extension term, as in current theory (García-Fiñana and Cruz-Orive, 1998;1999).To simplify the problem we have two main options.
Option 1.To model f so that q = m = 1, integer, and f (1) has finite jumps (Fig. 3b).The corresponding extension term is not satisfactory, see Fig. 3d, steeper broken line.The alternative choice m = 0 also fails (less steep broken line).
Option 2. To model f allowing f (1) to have infinite jumps at the ends of its support, see Fig. 3c, and to apply fractional theory.The result (Fig. 3d, red line) is far more satisfactory.

.1 THE NEW VARIANCE PREDICTOR
The classical representation of ( ) The extension term depends only on the smoothness properties of g at the origin, hence it is very effective to model g allowing for infinite jumps of its first non-continuous derivative at the origin.Consequently, we adopt the following model (also hinted in Matheron, 1971, for h ∈ ! ), Application of the general Euler-MacLaurin formula ( 14) to the first Eq.22 with the model 23 yields the following extension term, ( ) Note that when q = 0, 1 we obtain the well known expressions respectively.To estimate the extension term (24) we need to estimate q and b 2q+1 .Suppose that q is known.From a given sample {f 1 , f 2 , ..., f n }, (n ≥ 3) see Eq. 4, the usual technique (e.g.Kiêu, 1977;Cruz-Orive, 1999) yields the following estimator of b 2q+1 , ( ) whereby Eq. 24 yields the following fractional variance predictor, As particular cases we have the familiar values α(0) = 1/12, α(1) = 1/240, see Fig. 4a.Also, α(1/2) = ζ(3)/(8π 2 log2) ≈ 0.022.

APPLICATION
The material used was the white matter of a human brain, stemming from the study of McNulty et al. (2000), with kind permission from the authors.The basic data consisted of 177 consecutive coronal slice images of 1 mm thickness, obtained by non invasive magnetic resonance imaging techniques.The projected area of white matter in each slice was regarded as a planar section area, and measured by automatic pixel counting (pixel side = 1 mm).The observed section areas were therefore practically free from measurement errors, but the latter would not be a limitation of the method.The graph of the empirical ( ) ĈE Q , obtained by resampling from the 177 data points as indicated e.g. in McNulty et al. (2000), is represented in Fig. 4b, wavy line.The broken lines represent the predictions with q = m = 0 and 1, whereas the red line represents the fractional prediction computed via Eq.26a with q = 0.45.In turn, the latter was estimated from 88 sections 2 mm apart using Eq.27.

FINAL COMMENTS
Our results constitute an extension of the current theory, and they were motivated by practical and theoretical considerations.The geometrical object in Section 4 and the real one in Subsection 5.2 suggest that a fractional behaviour of the trend of ( ) ĈE Q might be a rule rather than an exception in practice.To use an integral smoothness constant (typically m = 1) with or without local measurement errors, seems to provide a poorer fit to ( ) ĈE Q than the fractional approach.The latter provides more freedom to model a measurement function respecting its main features, yielding a better extension term (Fig. 3).
A preliminary report of this paper was presented at the X th International Congress for Stereology, Melbourne, Australia, 1-4 November 1999.
the mean number of sections used for the Cavalieri estimator Q of the volume Q for each of the solids in Figs.1b

Fig. 1 .
Fig. 1.(a) Sweeping plane.(b-d) Solids described in the text, Section 2. (e) Corresponding measurement (area) functions.The red curve corresponds to the cylinder.(f) Graphs of the true coefficient of error of the Cavalieri estimator of the volume Q for each solid; the colours match those of (e).Note that the trends are fractional.

Fig. 2 ,
Fig. 2, left, shows a right circular cylinder and a plane whose normal makes an angle θ with the axis of the cylinder.For θ = 0 the measurement function f is a rectangular function and ( ) ( ) 2 Var Q O T = in general,

Fig. 2 .
Fig. 2. Left: Right circular cylinder and oblique sweeping plane.Right: Corresponding measurement (area) functions for two specific sectioning angles; the circles indicate the points at which f (2) exhibits infinite jumps.See text, Section 4.

Fig. 4 .
Fig. 4. (a) Graph of 1/α(q), where α(q) is given by Eq. 26b.(b) The fractional predictor of the CE, red line, obtained via Eq.26a, represents a clear improvement over the integral predictors used so far, see Subsection 5.2.