Geometric identities in stereological particle analysis

We review recent ﬁndings about geometric identities in integral geometry and geometric tomography, and their statistical application to stereological particle analysis. Open questions are discussed.


INTRODUCTION
Many identities in integral geometry and geometric tomography (Gardner, 1995) take the form where X is a subset of some Euclidean space, α and β are real valued functionals, T ranges over a class of sets such as the k-dimensional planes, and typically dT represents integration with respect to an appropriate 'uniform' measure.The simplest example is the representation of the volume α(X) = V (X) of a three-dimensional solid X as the integral of the areas β (X ∩ T ) = A(X ∩ T ) of its horizontal plane sections X ∩ T .
Such identities may be given a stochastic or statistical interpretation, which is the basis of the science of stereology (Baddeley and Jensen, 2005;Weibel, 1979;1980).For example, the representation of volume in terms of the areas of plane sections was used by 19th century geologist A.E. Delesse to obtain a practical method for determining the composition of rocks.Namely, the volume percentage (fraction of volume of rock) occupied by a particular mineral of interest can be estimated from the fraction of area occupied by the same mineral in a single plane section of the rock.Note that this is different from applications to computed tomography, which require information from all plane sections of the object.
Until the 1980's, stereological techniques were based mainly on the classical section formulae of integral geometry, taking the form (Eq. 1), where α and β are the intrinsic volumes or quermassintegrals in R d for d = 1, 2, 3, and T is a k-dimensional plane, k < d.These formulae, and the stochastic interpretations that were placed on them, make it possible to statistically estimate volumes, areas and lengths of three-dimensional structures using information obtained from random plane sections of the structure (Weibel, 1979;1980).

Stereologists use the generic term 'particles'
to describe solid objects that can be scientifically interpreted as discrete individuals (such as biological cells, mineral grains, or enclosed holes).A 'particle population' is a finite or countable collection of disjoint particles.
The 'particle problem' in stereology is the problem of inferring the number, average size, and distribution of sizes, of a population of particles, from a plane section of the population.For example, given a microscope image of a plane section of brain tissue, we may wish to estimate the number of neurons in the brain, their average size, and the distribution of their sizes.
A general solution to the particle problem was not found in classical stereology.The classical identities of integral geometry (and their classical stochastic interpretation) only allow us to estimate 'aggregate' quantities α(X) where X = i Y i is the union of all the particles Y i in the population, and not 'individual' particle properties such as the distribution of particle sizes α(Y i ) .
Until the 1980's, progress on the particle problem was achieved mainly by assuming that the particles have a common, simple, known shape such as spheres.Swedish mathematical statistician S.D. Wicksell first formulated the particle problem, and treated the case of estimating the size distribution of spheres in three dimensions from two-dimensional plane sections (Wicksell, 1925).In R n , the relation between the distribution function F of the radii of the spheres and the distribution function F (n,k) of the radii of the spheres sectioned with a k-dimensional plane is given by, cf.Schneider and Weil (2000), where M n−k is the (n − k)th moment of sphere radii assumed positive and finite.Wicksell's problem has attracted considerable attention among mathematical statisticians, one reason being that it is an ill-posed problem implying that a small perturbation of the data in the section may lead to a large change in the estimated sphere size distribution.
Geometrical assumptions about particles -for example that the particles are all spheres -are overly simplistic for most applications.It is therefore a major advance that with the introduction of new sampling and measurement techniques such fragile geometrical assumptions can be relaxed.An example of these new techniques is 'local' sectioning of each particle through a reference point in the particle.Stereological methods based on such local sections constitute the new field of local stereology which is closely related to geometric tomography, cf.Gardner et al. (2003).Moments of particle size (including mean and variance) can be estimated stereologically without specific shape assumptions, using information on such local sections.(As a side remark, Wicksell's problem becomes trivial with such information at hand since the radius of a sphere can be observed directly in a section through the centre of the sphere.) Prompted by the need to utilize the new type of data, new geometric identities have been discovered or rediscovered during the last couple of decades and have been used to renew stereological particle analysis.In the present paper we review these geometric identities in integral geometry and geometric tomography (Gardner, 1995) and point to some missing ones.
The first section describes the key geometrical identities.The second section discusses how the geometric identities can be used in the stereological analysis of particle populations.In section three we discuss inference about the particle size distribution.The fourth section describes approaches to modelling shapes.Finally in the last section we discuss some open problems.
An earlier review of stereology for geometers can be found in Weil (1983).Subsequently the development of three-dimensional microscopy, confocal microscopy and other imaging modes has led to a revolution in stereological sampling and measuring techniques.This involves both the development of new geometric identities, and alternative interpretations of existing geometric identities.A comprehensive treatment of stereology for statisticians has recently appeared in Baddeley and Jensen (2005), but for readers interested in the geometric background there seems to be no up-to-date review.

GEOMETRIC IDENTITIES THE BLASCHKE-PETKANTSCHIN FORMULAE
The Blaschke-Petkantschin formulae play a fundamental role in modern stereological particle analysis.They provide the mathematical tool for estimating moments of particle volume from information on local p-dimensional sections through the particles.

The
Blaschke-Petkantschin formulae are geometric measure decompositions of the q-fold product of Hausdorff measure in R n .In the special case of decomposition of one copy of n-dimensional Lebesgue measure V n , the Blaschke-Petkantschin formula takes the following form for X ⊆ R n compact and p > r: The normalizing constant is given by c is the surface area of the unit ball in R n .Note that by conventioyn c(n, 0) ≡ 1.The outer integration in (Eq. 3) is over the set L n p(r) of all p-dimensional linear subspaces L p (called for short p-subspaces in the following) of R n , containing a fixed r-subspace L r , say.Note that L n p(0) = L n p , the set of p-subspaces in R n .The differentials dx p and dL n p(r) are defined as dx p = H n p (dx) and dL n p(r) = µ n p(r) (dL p ), where H n p is p-dimensional Hausdorff measure in R n and µ n p(r) is the unique measure on More generally, suppose that X 1 , X 2 , . . ., X q are compact subsets of R n .For p ≥ q + r, we have where ∇ q (x 1 , . . ., x q ) = q!H n q (conv{O, x 1 , . . ., x q }) and conv{O, x 1 , . . ., x q } denotes the convex hull of the set {O, x 1 , . . ., x q }.Here, the point Early versions of the Blaschke-Petkantschin formulae for decomposition of Lebesgue measures are due to Blaschke (1935a;b) and Petkantschin (1936).The one presented above has been derived in Miles (1979).The general decomposition of the q-fold product of d-dimensional Hausdorff measure in R n has been established in Zähle (1990) and Jensen and Kiêu (1992).For d = n the decomposition reduces to the one presented above while for d < n it involves the so-called G-factors which contain information about the angle between the boundary of the sets X i and the intersecting p-subspace.
The formula for decomposition of Lebesgue measures has been proved by invariant measure theory (Møller, 1987) while the general formula for ddimensional Hausdorff measures has been established using geometric measure theory, see Zähle (1990) and Jensen and Kiêu (1992).A simplified proof, also utilizing induction in the dimension of the spaces involved, can be found in Jensen (1998, Ch. 5.4).

A LOCAL SLICE FORMULA
For d < n, the Blaschke-Petkantschin formulae depend on angles defined in R n which cannot be determined only from information on lowerdimensional sections X ∩ L p .Also in the case d = 0 the decomposition is trivial so it is not possible to get information about number using the Blaschke-Petkantschin formulae.This has been the motivation for deriving geometric identities for p-dimensional slices centred at a reference point.Such slices are psubspaces with some thickness.

The local slice formula gives a geometric decomposition of d-dimensional
we denote the set of pslices L p + B n (O,t) for which L p ⊇ L r , where L r is a fixed r-subspace.Similar as in the subsection before the differential dT n p(r) is defined as Let X be a d-dimensional differentiable and bounded manifold in R n .Then, for arbitrary p, r satisfying 0 ≤ r < p < n, the local slice formula takes the form where An early version of this geometric identity can be found in Jensen and Kiêu (1994).The proof of the slice identity uses the following reasoning that holds for any x ∈ R n : ) .At the last equality sign it is used that the length of the isotropic projection of a line segment is Beta distributed.The more general issue of isotropic projections of simplices has been studied in Nielsen (1999).Here the distribution becomes that of a product of independent Beta-distributed random variables.

A GEOMETRIC IDENTITY INVOLVING VERTICAL SECTIONS
Vertical sections are hyperplanes parallel to a fixed axis.Such sections are often used in microscopy to reveal structural information or for practical convenience.In the 1980s, Baddeley (1984;1985;1987) and collaborators (Baddeley et al., 1986) showed how a rotation invariant line in R n can be generated via a vertical section.The important consequence is that results based on rotation invariant lines can be transformed into results for lines generated on vertical sections.We will in this subsection present this key result and show how it can be used to derive an alternative local geometric identity for volume.The main application of Baddeley's result has, however, been in the development of a stereological estimator for surface area.This estimator has also relevance for stereological particle analysis, as shown in section two.Another important issue are vertical probes, which have been developed to vertical projections by Gokhale (1990).We do not consider vertical probes in this paper but for the interested reader we refer to a survey in Beneš and Rataj (2004).
Then, the key result proved by Baddeley can be stated as follows: where In Beneš and Rataj (2004, Ch. 4.1.3),(Eq.4) has been proved by first establishing the following identity, using the coarea formula and spherical coordinates, where S n−2 (v ⊥ ) is the unit ball in the hyperplane v ⊥ through the origin, perpendicular to v, span{l, v} is the linear subspace of R n , spanned by l and v, S 1 See Fig. 1 for an illustration in 3D.Note that (Eq.4) is indeed implied by (Eq.5) since Fig. 1.Illustration relevant for the identity (Eq.5), relating to vertical sectioning.
As announced in the beginning of this subsection, Eq. 4 has been used to derive a local identity for volume.Combining with the Blaschke-Petkantschin formula (Eq. 3) with p = 1 and r = 0, we get where

A GEOMETRIC IDENTITY FOR SURFACE AREA
The Blaschke-Petkantschin formulae depend for d < n on angles in R n which cannot be determined on a lower dimensional section.As an example, for an (n − 1)-dimensional manifold X we have where | and α d,p is a generalized angle between the subspaces L d and L p .Furthermore, Tan[X, x] is the tangent space of X at x ∈ X.In Jensen (1998, Ch. 5.6), it is shown that the geometric identity can be modified such that only lower-dimensional Image Anal Stereol 2006;25:63-74 angular information is used.The modified geometric identity is of the following form: where is the angle between L 1 and the unit normal vector of the tangent space of X ∩ L p at x, and h can be expressed as a hypergeometric function Very recently, an alternative approach to estimate surface area has been suggested by Cruz-Orive (2005).
The method is based on a new principle to generate isotropic uniform random test lines hitting a particle.

STEREOLOGICAL PARTICLE ANALYSIS
We will now discuss how the geometric identities can be used in the stereological analysis of particle populations.

STOCHASTIC MODEL
A 'particle' is a compact, nonempty subset of R n .For instance, the cells in a biological tissue or the inclusions in steel may be modelled as a collection of particles in R 3 .To study the particle problem we must consider the intersection between a collection of particles and a 'probe' T such as a plane or a line.Either T is randomly positioned (a "design-based" approach) or the collection of particles is random (a "model-based" approach).In this review we focus on the model-based approach.D. Stoyan, J. Mecke and collaborators (Hanisch and Stoyan, 1981;Mecke and Stoyan, 1980;Stoyan, 1979;1982;1990;1996;Stoyan and Mecke, 1983) formulated the particle problem using the theory of point processes (Daley and Vere-Jones, 1988).First consider the very simple case of a collection of spheres B(x i , r i ) in R 3 with random centres x i ∈ R 3 and random radii r i > 0. The collection of pairs (x i , r i ) can be treated as a point process in R 3 × R + .Equivalently this is a 'marked point process' in R 3 with marks in R + .
A marked point process in R n with marks in some space M can be defined formally as a point process in R n × M satisfying certain conditions.It is interpreted as a point process in R n where each point x i ∈ R n of the process carries additional information in the form of a mark m i ∈ M .For details see e.g., Stoyan et al. (1987).
The class K of all particles is a complete separable metric space, so a random particle (a random element of K ) can be defined in a straightforward manner (Matheron, 1975).A random collection of particles is most conveniently modelled as a marked point process in R n with marks in K .
The particles are regarded as a realization of a marked point process Ψ = {[x i ; Ξ i ]} in R n , where the x i 's are points in R n and the marks Ξ i are compact subsets of R n .The ith particle of the process is represented by the set X i = x i + Ξ i .In this framework, x i is called the nucleus of the ith particle and Ξ i the 'primary' or 'centred' particle.The corresponding process of nuclei is denoted by The marked point process Ψ is assumed to be stationary, i.e., for all x ∈ R n we have Ψ x ∼ Ψ, where Since Ψ is stationary, the intensity measure of Ψ can be written as is the mean number of particles per unit volume with marks in K.The mark distribution is defined as where N V ∈ (0, ∞) is the intensity of Ψ n .Using Eq. 6 and Eq. 7 it can be shown for any measurable nonnegative function h defined on R n × K , In what follows, we let Ξ 0 be a random compact set with distribution P m .If Ψ is isotropic, then AΞ 0 ∼ Ξ 0 for all rotations A.

GEOMETRIC SAMPLING EFFECTS
When a population of particles is sectioned by a plane, two 'sampling effects' occur.First there is sizedependent sampling bias: the probability that a given particle is hit by the section plane depends on the size of the particle.Second, there is random reduction in size: a plane section of a particle is smaller than the particle itself.
These effects were first explained by Wicksell (1925) in the case of spherical particles.Consider a stationary process of spheres in R 3 , intersected by a fixed two-dimensional plane.The plane section is a stationary process of circles in R 2 .Size-dependent sampling bias occurs because, roughly speaking, the probability that a given sphere is intersected by the section plane is proportional to the sphere's radius.Let N V denote the expected number of sphere centres per unit volume in the population, and N A the expected number of circle centres per unit area in the plane section.Then N V and N A are not equal, but instead are related by where ED is the mean sphere diameter in the original particle process.If the sphere radii have cumulative distribution function F, then the spheres selected by a random section plane have radii with distribution function the size-biased counterpart of F. Here, ER denotes the mean sphere radius in the original particle process.
For particles of more general shape, it was established by DeHoff and Rhines (1961) and others, then in complete generality by Stoyan (1979), that where N V is the intensity of the particle process (expected number of particle centres x i per unit volume), N A is the intensity of the process of section profiles (the expected number of particle profiles per unit area in the plane section) and EH is the mean particle height.The height H of a particle is the length of its projection onto the subspace normal to the section plane.Relation (Eq.11) is known as the 'Rhines-DeHoff equation'.See Fig. 2 for an illustration in 3D.

PARTICLE NUMBER
An important task in stereology is the estimation of N V , the mean number of particles per unit volume of the material of interest.Unfortunately, it is not possible to estimate N V , in complete generality, from a plane section alone.This can be seen from the Rhines-DeHoff relation (Eq.11).The population average height of the particles EH cannot generally be determined from a plane section.The problem can be circumvented only when we have extra information, for example when we have some knowledge about the particle shape, or when particle height measurements are available.
In order to construct an unbiased estimator for N V in a situation where the particles are of general shape, we need a different sampling scheme.Consider a sampling scheme with the property that for any fixed sampling window W ∈ K with positive volume and any particle X ∈ K , Then, using Eq. 8 and Eq. 12, we conclude that Therefore an unbiased estimator for N V is given by In practice Miles' associated point rule (Miles, 1974) and Gundersen's tiling rule (Gundersen, 1977;1978) satisfy the basic requirement (Eq.12).See Fig. 3 for an illustration in 2D of these two sampling rules.

MOMENTS OF PARTICLE SIZE
The geometric identities presented in the first section can be used to derive estimators of moments of particle size Eϕ(Ξ 0 ) q , q ≥ 1, where ϕ is a selected size parameter.Below we give two examples of such estimators.

Estimators based on the Blaschke-Petkantschin formula
Suppose that the distribution of Ξ 0 is invariant under the action of SO(n, L r ).Using Eq. 3, we get Since the distribution of Ξ 0 is invariant under rotations around L r , the inner integral of the right-hand side of (Eq.13) does not depend on L p and we get for a particular p-subspace L p0 ∈ L n p(r) , is an unbiased estimator of EV n (Ξ 0 ).
In practice, a sample of particles {x i + Ξ i : x i ∈ W } is collected in a sampling window W and a central section is determined through each particle.An illustration of this sampling scheme for n = 2 and p = 1 is shown in Fig. 4. For each sampled particle, we determine Using Eq. 8, we find that We can thereby construct an unbiased estimator of EV n (Ξ 0 ) if an estimate of N V can be constructed.

Estimators based on vertical sections
In the same way, other identities from the first section can be used to construct estimators of EV n (Ξ 0 ) q or ES n (Ξ 0 ).
The identity (Eq.4) involving vertical sections provides an alternative for estimating ES n (Ξ 0 ) if the distribution of Ξ 0 is invariant under rotations around the vertical axis v.The identity is used in combination with the classical Crofton formula.Let us suppose that the particles are non-overlapping and let Then, using (Eq.8) we find that Image Anal Stereol 2006;25:63-74 only three types of possible limiting distributions M (de Haan, 1975;Galambos, 1987), namely and M 1,α , M 2,α and M 3 are called Fréchet, Weibull and Gumbel distributions respectively.
For spherical particles, Eq. 2 has been used to show that A proof of these results for n − k = 1 can be found in Drees and Reiss (1992).The proof for n − k ≥ 2 is similar.
These stability properties can be used to estimate the normalizing constants of F (n,k) from data in the section plane L k and the aim is then to get back to the normalizing constants of F and predict its extremes.In a series of papers this was done in the spherical case, using a generalized gamma model for the distribution function F (Takahashi and Sibuya, 1996;1998;2001;2002).
A slight extension of this approach was recently obtained for spheroidal particles in R 3 (Hlubinka, 2003a;b).Suppose the particles are ellipsoidal with two equal major semiaxes whose length is X.Denote the minor semiaxis by W . Then a particle is completely determined (up to position) by the bivariate variable (X, T ), where is called a shape factor.Note that the spherical case is obtained for T = 0 and deviation from spherical shape increases with T .Under certain isotropy and independence assumptions the observed ellipses in the section plane can be described by a size-shape variable (Y, Z), where Y is the length of the major semiaxis of the ellipse and Z is a shape factor similarly defined as in (Eq.16).Assuming that (X, T ) has a density f (x,t), one can derive a relation between f (x,t) and the density g(y, z) of (Y, Z), see Cruz-Orive (1976).Using this relationship it is possible to prove stability properties as in the spherical case assuming a fixed size or fixed shape for the spheroidal particles.A simulation study where f (x,t) follows a continuous Farlie-Gumbel-Morgenstern type of distribution can be found in Beneš et al. (2003).
Instead of considering maxima one could also be interested in minima.Stability properties for minimal radii in the spherical Wicksell corpuscle problem has been proved recently in Kötzer and Molchanov (2006).

SHAPE MODELLING
Shape modelling of planar and spatial objects with no obvious landmarks has attracted much attention in the last ten years.One approach in this direction uses the normal deformation of a sphere, which is defined via the spherical-harmonic basis on the unit sphere, see e.g., Grenander and Miller (1994);Hobolth (2003); Hobolth and Jensen (2002).This shape model can easily be generalized to higher dimensions.Thus, let Ξ 0 ⊆ R n be star-shaped relative to z ∈ Ξ 0 .The boundary of Ξ 0 is determined by where r(ω) is the distance from z to the boundary of Ξ 0 in direction ω.We can express the radius-vector function r(ω) in terms of the spherical harmonics which constitute an orthonormal basis on S n−1 (Groemer, 1996;Mueller, 1966).Here N(n, k) denotes the number of linearly independent spherical harmonics of degree k in n variables.The Fourier-Legendre series expansion of the radius-vector function is given by For modelling shape variability, the coefficients a where ω 1 , ω 2 ∈ S n−1 and P (n) k (•) denotes the Legendre polynomial of degree k in n variables.¿From (Eq.18) we see that assumption (Eq.17) implies stationarity on the sphere, in the sense that the covariance between two points on the sphere depends only on the angle between these points.Since in (Eq.18) σ n , N(n, k) and P in the planar and spatial case (n = 2, 3) are discussed in Grenander and Miller (1994;1998); Hobolth and Jensen (2002).

OPEN PROBLEMS
The list of geometric identities presented in this paper is not complete in the sense that it would be of great practical interest in stereological particle analysis if identities of the type (Eq. 1) could be derived for functionals α other than powers of volume and surface area.Such an identity is, however, only of interest in applications if it is possible to determine β (X ∩ T ) from information inside the section T .In particular, it is important to derive identities of this type with α equal to squared surface area.Another class of identities still missing are those with α equal to intrinsic volumes in R n or β equal to intrinsic volumes on the sections.Some progress on the second problem has been made in Jensen and Rataj (2005) but the first problem remains open.In fact, the first problem is the one with most practical interest.
Prediction of extreme particle sizes in the case of spherical or spheroidal particles is the first step in creating a theory of "Stereology of extremes", but having applications in mind it is clearly of great importance to develop methods for particles with more complicated shape.A major problem is that only for simple particle shapes is it possible to derive relationships between the particle sizes of interest such as (Eq.2).For example, Cruz-Orive (1976) has shown that the size distribution of ellipsoidal particles in R 3 with three different semiaxes cannot be uniquely determined using information from plane sections.

Fig. 2 .
Fig. 2.An explanation that Q must depend on the mean particle height EH.The two boxes contain equal numbers of particles.The particles on the left are smaller than those on the right.A typical section plane cuts fewer particles in the left box than it does on the right.

Fig. 4 .
Fig. 4. Local sampling scheme.Particles x i + Ξ i satisfying x i ∈ W are sampled using a fixed line L 10 .
k,m ∼ N(0, λ(n) k,m ) , k ∈ N 0 , m = 1, . . ., N(n, k) .Assuming moreover that λ (n) k,m = λ (n) k , k ∈ N 0 , m = 1, . . ., N(n, k) , ) are known, the covariance is completely determined by the variances λ (n) k .Parametric models for the variances λ (n) k the subgroup of rotations in R n , keeping the r-subspace L r fixed.Furthermore, the symbol || • || in (Eq. 3) denotes the Euclidean norm in R n and π L ⊥ r is the orthogonal projection onto L ⊥ r , the orthogonal complement of L r .