Image Anal Stereol 2000;19:71-79 Original Research Paper FRACTIONAL TREND OF THE VARIANCE IN CAVALIERI SAMPLING Marta García-Fiñana and Luis M Cruz-Orive Department of Mathematics, Statistics and Computation, Faculty of Sciences, University of Cantabria, E-39005 Santander (Accepted May 31, 2000) ABSTRACT 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 T2m+2 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 T2q+2, typically with 0 R+ is a bounded non-random function, called the measurement function, which is integrable on a finite support [a, b] c R 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/is (m, /?)-smooth, a concept we recall briefly. Given a function h : R -> R+, the amplitude of the jump of the Mh derivative h(k) of h at the abscissa is expressed as follows, 1S*/2(*)(x)=lim/2W(j)-lim/2W(>'), (xel, £ = 0,1,...), (2) y^>x+ y^>x provided that the limits exist. Further, Dh(k)={x1\Sh(k)(x1)^0,ieZ], (£ = 0,1,...), (3) 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, /?)-smooth, with smoothness constant m = 0,1,... and/? = 1,2,..., both known and fixed, if the support of h is bounded and if García-Fiñana et al: Fractional trend of the variance in Cavalieri sampling a. Dh(k) = 1)-smooth. Then, Var( ˆ )=Var£( ˆ ) + Z (r) + o(r2m+2), (5) where Var£ (q) is called 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, Var£( ˆ ) = ( -1f-r2m+2-i?2m+2 ( 0 ) £ (sy(m) ( c ) )2 , (6) (Souchet, 1995; Kiêu, 1997), where P2i(0) = B2l(0)/(2l)! and B2[(0) is a Bernoulli number, (Abramowitz and Stegun, 1965). The preceding results, however, are valid if/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. 2. 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 ti) of the form, f(x) 1-\x\") ", xe[-1,1l a>0, (7) 0, x<£ [-1,1] see Fig. 1e. This measurement function is continuous, that is Df (0) = 0. However: For 0 < a < 1, Df 1) = {0} and Sf )(0) is unbounded. For 1 < a < 2, Df 1) = 0, Df 2) = {-1, 0, 1}, 5^2)(0) is not defined and Sf 2)(c) is unbounded for c e {-1, 1}. For a > 2, Sf (1)(c) is unbounded for c e Df 1) = {-1, 1}. Thus, the current theory can be applied only when a = 1, 2, because then the jumps Sf )(c) are all finite. For the right circular cylinder (Fig. 1d) the measurement function is (up to a constant factor), 72 Image Anal Stereol 2000;19:71-79 fix) l-x' xe[-1,1], ä[-1,1], (8) see Fig. le, red curve. This measurement function is continuous, that is Df (0) = 0. However, Sf (V)(±\) = +00, and therefore the current theory cannot be applied either. The wavy lines in Fig. If represent the true Ce(q) as a function of the mean number of sections used for the Cavalieri estimator Qˆ of the volume Q for each of the solids in Figs. 1b-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. 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. 3. GENERAL VARIANCE REPRESENTATION 3.1. 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 = 0, (k = 0, 1, ..., q - 1), where q is the smallest integer > q. b. Dfirq) = {c„ i e N} is a finite set, and for each ci we can find two real constants d~ , di , and corresponding rationals a: ,a\, (not necessarily different), such that: 0, 73 García-Fiñana et al: Fractional trend of the variance in Cavalieri sampling lim (ix-ciK ¦f{q)(x)) = df, af e [0,1). We have that q is precisely q-a, with a = max la. , a. ), and q is called the fractional smoothness constant of f Note that if all af=0, then the jumps of f{q) at the points {ci} are all finite, and we recover the integral smoothness case (Section 1) with p = 0. The Riemann-Liouville definition of the q–Ha. fractional derivative of f (x), (a < x < b), reads, dqf(x) dq-1xq q+\ dx q+\ 1 T(l +q-q) •s x — u\ ¦f(u)-du (10) (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 T(-) is the Gamma function. With this definition the constants df in the right hand side of Eq. 9 can be connected with the fractional derivatives of f at the corresponding discontinuity points as follows, f q i ci ) = r(i-<) di (ii) 3.2 GENERAL EULER-MACLAURIN FORMULA The estimator Qˆ is a periodic function of x1 with period T, and it therefore admits the Fourier expansion Q(xi)= 2, akQx$(2nikxllT), where ak = \ f(x)-exp(-2nikx/T)dx (12) Bearing in mind that ak = a_k we obtain (ink Q{xl)-Q = 2jJ\ab f{x)cosMT (x-xi) ]dx. k=i (13) If is q–smooth, from the preceding expression we obtain a new, generalized version of the Euler-MacLaurin formula, namely Q ( x0-Q = - q rfq c Pq T x-c+Hq f^c-Pq Tc-xOTq+ofTq1) (14) where PyT (x) = P (x/T), and 1 cos 2njx-—ny P(x) = ^2-f----Lj?----L , 00 ) . (16) The extension term has the following expression P (O} 2 VaE(Q)= (-lq-T2q+2------2q+2 ; Y (Sf{q )( c )), (17) V ' cos(7v(q-q))ceDfqy ' where Sfq) is the ’fractional jump’ of f q), defined as follows {Sfiq) (c))2 := {f{q) (c+ ))2 + (f (c- ))2 - 2fq (c+ )f (c" )cos (* (q -q)). (18) To obtain Eq. 17 we have used the following identities, E[^r(^-c)] = Efer(c-x)] = P q+2(0)/cos(^), EPq TCx-ci-Pq Tfcj-x^-P qfcj-ci). (19) If q = m e N {0}, then cos (tt (q - q)) = 1 and f(m) (c+) - f(m) (c") = Sf{m) (c), whereby Eq. 17 reduces to Eq. 6. 3.4 EXAMPLES We illustrate the applicability of Eq. 17 to obtain the extension term for the measurement function (8). In virtue of Eq. 9, d~ = lim ( ( 1 - xf ¦ f{l) ( x )\ = —\=, with a2 =-. (20) Likewise, d? =l/>/2 with < =1/2, whereas d; =d+2 = 0, (a; =a2+=o), because f (x) = 0 for xz [-1,1]. Thus, q = 1 - max(0, 1/2, 1/2, 0) = 1/2, hence f is (l/2)-smooth. By Eq. 11, f(1/2) (-l+) = -f(1/2) (r ) = JOll whereas f(1/2) (-1") = f(1/2) (l+) = 0. Thus, by Eq. 18, Sf{ll2) (-1) = Sf{ll2) (l) = jnll, and finally Eq. 17 yields VaE Q=T3-----3±L- \- + -\=^1.T\ (21) v ; cos(tt/2) ^2 2) An where C(s) = X^ik~s , (s > *) is Riemann’s Zeta function (e. g. Abramowitz and Stegun, 1965). The graph of the extension term (21) is shown in Fig. If, broken red line, and it shows the trend of Var(Q) rather faithfully. For the measurement function (7) the same method yields VarE ( q ) = O ( T1413 ) for a = 3/2, (Fig. If, broken green line), and VarE ( q ) = O ( t512 ) for a = 8, (Fig. If, broken blue line). 4. GENERAL VALIDITY The special measurement functions (7) and (8) arise when the solids in Figs. lb, c and in Fig. Id, respectively, are cut at one particular orientation (Fig. la). 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. 75 García-Fiñana et al: Fractional trend of the variance in Cavalieri sampling 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 Var ( Q ) = O ( T2 ) in general, whereas for ? = p/2 we have seen that Var ( Q ) = O ( T3 ) , (Fig. 1). For 0 < ? < p/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. 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. 3. (a) Area function of the cylinder in Fig. 2 for ? = 7°. (b) Model to comply with current theory. (c) Model to comply with fractional theory. (d) True Ce(q), (wavy line), approximation using (b), (steeper broken line), approximation using (c), (red line) and step function case (less steep broken line). 76 Image Anal Stereol 2000;19:71-79 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 =\, integer, and fl) 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 fm 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. 5. VARIANCE ESTIMATION AND APPLICATIONS [5] 5.1 THE NEW VARIANCE PREDICTOR The classical representation of Var(Q) in terms of the covariogram g of f reads Var(Q) = T]T g ( kT ) -J"°° g ( h ) dh, g ( h ) = ff ( x ) f ( x + h ) dx,heR. (22) k=-oo 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 heR), g{h) = b0+b2q+l-\hfq+l +b2-h2, he[a-b,b-a], qe[0,\]. (23) Application of the general Euler-MacLaurin formula (14) to the first Eq. 22 with the model 23 yields the following extension term, NaxE[Q) =-2-b2q+l-T{2q + 2 )- P2q+2{ 3) see Eq. 4, the usual technique (e.g. Kieu, 1977; Cruz-Orive, 1999) yields the following estimator of b2q¥1, IC -AC +C ^j b q^ =-r5-------\------, Cj =Tfifi+„ (j = 0,1,...,n-1), (25) yz —tj-T i=\ whereby Eq. 24 yields the following fractional variance predictor, Var(Q) = a( q )-( 3C0-4C+C2 )-T2, (26a) a(q)= ' q+ ,2q+2' = ^ q+ ',\q+ >co^qn\ (26b) V ! 2-22q (2n) q \\-2lq-1 ) As particular cases we have the familiar values a(0) = 1/12, a(\) = 1/240, see Fig. 4a. Also, a(l/2) = C(3)/(87t2log2)~ 0.022. An estimator oiq is obtained similarly as in Souchet (1995), Kiêu (1997), and Cruz-Orive (1999), Eq. 2.19, namely, q = log *^C/-\ ^Ca I C A 1 1 log(4) 2 (27) Note that now q is not rounded off; then a(q) is calculated via Eq. 26b. If there are local measurement errors (’nugget’), then a(q) replaces a in Cruz-Orive (1999), §2.3. 77 García-Fiñana et al: Fractional trend of the variance in Cavalieri sampling Fig. 4. (a) Graph of 1/a(q), where a(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. 5.2 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 CE(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. 6. 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 CE(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 CE(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 Xth International Congress for Stereology, Melbourne, Australia, 1-4 November 1999. 78 Image Anal Stereol 2000;19:71-79 ACKNOWLEDGMENTS The authors are grateful to Ms Vicky McNulty and Dr Neil Roberts for the data used in Subsection 5.2. This research was supported by the DGES Project No. PM97–0043, and Biomed II (E.U.) No. BMH4–CT97–2437. REFERENCES Abramowitz M, Stegun IA (1965). Handbook of Mathematical Functions. New York: Dover. Cruz-Orive LM (1999). Precision of Cavalieri sections and slices with local errors. J Microsc 193:182-98. García-Fiñana M, Cruz-Orive LM (1998). Explanation of apparent paradoxes in Cavalieri sampling. Acta Stereol 17:293-302. García-Fiñana M, Cruz-Orive LM (1999). New approximations for the efficiency of Cavalieri sampling. Internal report 3/99. University of Cantabria: Department of Mathematics. J Microsc, in press. Gardner M (1978). Mathematical Carnival. Harmondsworth: Penguin. Gundersen HJG, Jensen EB, Kiêu K, Nielsen J (1999). The efficiency of systematic sampling in stereology - reconsidered. J Microsc 193:199-211. Kiêu K (1997). Three lectures on systematic geometric sampling. Memoirs 13/1997. University of Aarhus: Department of Theoretical Statistics. Kiêu K, Souchet S, Istas J (1999). Precision of systematic sampling and transitive methods. J Statist Planning Inf 77:263-79. Matheron G (1965). Les Variables Régionalisées et Leur Estimation. Paris: Masson et Cie. Matheron G (1971). The Theory of Regionalized Variables and Its Applications. Les Cahiers du Centre de Morphologie Mathématique de Fontainebleau, No 5. Fontainebleau: École Nationale Supérieure des Mines de Paris. McNulty V, Cruz-Orive LM, Roberts N, Holmes CJ, Gual-Arnau X (2000). Estimation of brain compartment volume from MR Cavalieri slices. J Comput Assist Tomogr, in press. Nishimoto K (1984). Fractional Calculus. Koriyama: Descartes Press CO. Roberts N, Garden AS, Cruz-Orive LM, Whitehouse GH, Edwards RHT (1994). Estimation of fetal volume by MRI and stereology. Br J Radiol 67:1067-77. Souchet S (1995). Précision de l’estimateur de Cavalieri. Rapport de stage, D.E.A. de Statistiques et Modéles Aléatoires appliqués à la Finance. Université Paris-VII. INRA-Versailles: Laboratoire de Biométrie. 79