Image Anal Stereol 2015;34:147-159
Original Research Paper
doi:10.5566/ias.1324
THE INVARIATOR DESIGN: AN UPDATE
Luis M. Cruz-Orive1 and Ximo Gual-Arnau^'2
1 Department of Mathematics, Statistics and Computation, Faculty of Sciences, University of Cantabria, E-39005 Santander, Spain; 2Department of Mathematics-INIT, University Jaume I, E-12071 Castello, Spain. e-mail: luis.cruz@unican.es, gual@uji.es
(Received April 30, 2015; revised July 23, 2015; accepted September 21, 2015) ABSTRACT
The invariator is a method to generate a test line within an isotropically oriented plane through a fixed point, in such a way that the test line is effectively motion invariant in three dimensional space. Generalizations exist for non Euclidean spaces. The invariator design is convenient to estimate surface area and volume simultaneously. In recent years a number of new results have appeared which call for an updated survey. We include two new estimators, namely the a posteriori weighting estimator for surface area and volume, and the peak-and-valley formula for surface area.
Keywords: invariator, peak-and-valley formula, stereology, surface area, test line weighting, volume.
INTRODUCTION
Since the publication of the invariator design (Cruz-Orive, 2005, although the name 'invariator' was first coined in Cruz-Orive, 2009b), a number of related papers have appeared with additional results - hence we thought that an update would now be opportune.
Briefly, the invariator design consists of two stages. In the first stage a randomly oriented plane, called the pivotal plane, is taken through a fixed point O, called the pivotal point. In the second stage a randomly oriented test line is sampled in the pivotal plane with a weight proportional to its distance r from O. This test line is effectively equipped with the motion invariant density in space, whereby surface area and volume, for instance, can be estimated unbiasedly by design. The advantage of the invariator is that the necessary observations can be made solely in a pivotal plane which has to be randomly oriented, but not randomly located.
In the original paper the aforementioned r-weighting was made a priori, namely it was implicit in the design: each test line was drawn through a point from a uR grid on the pivotal plane, with a direction perpendicular to the axis joining the point with O. Here a new, more natural procedure, based on a posteriori r-weighting, is proposed which requires only a system of parallel test lines (Section Test lines with a posteriori weighting).
In Cruz-Orive (2005) a method was also proposed to estimate the surface area of a convex object by measuring Feret rays emanating from O in a pivotal section. Recently, Thorisdottir and Kiderlen (2014) have generalized the result to non convex objects, see
also Thorisdottir et al. (2014). Following a different approach, Gual-Arnau and Cruz-Orive (2015) have obtained a simplified version of the relevant estimator, see Section Case of a general object: the peak-and-valley formula.
To make the survey widely accessible, an informal description of concepts like 'invariant density' is given in the next section. For completeness, elementary proofs of relevant stereological results are given in the appendixes.
PREREQUISITES AND NOTATION
Here we introduce basic tools of integral geometry (Santalo, 1976; De-Lin, 1994) for geometric sampling. We concentrate on points, straight lines and planes to be used as test probes, namely on geometric elements equipped with a well defined mechanism of randomness, used to probe (namely to hit, or intersect, with a sampling purpose) a target subset in Euclidean space.
A fundamental prerequisite is to find the measure of a set of geometric objects satisfying a given condition - for instance of test points belonging to a given region, of intersections between a surface and all possible test lines, etc. - which is invariant under the group of rigid motions (namely of rotations and translations). In other words, the measure should not depend on the choice of the reference frame. Consider for instance a bounded interval Y of the real line R1. The problem is to find a translation invariant measure for sets of points of abscissa x, namely a measure ^ such that an integral of the form
147
Cruz-Orive LM et al: The invariator design
M (Y) = j M (dx)	(1)
does not depend on the choice of the origin of abscissas. The preceding integral may be interpreted as the measure of points in Y. Here we consider a measure element m (dx) associated with the length element dx which may be expressed as m (dx) = w(x)dx, where w(x) is a non negative weight function. The problem is to find w(x) such that, for any translation vector z G R1, the identity
w(x) dx = w(x) dx	(2)
Jy	JY+z
holds, where Y + z represents the translate of the interval Y by the vector z. in the left hand side integral make the change of variable x = y — z. Then,
w(x) dx =/ w(y — z) dy .	(3)
JY	JY+z
The right hand sides of the preceding two identities must coincide for all z G R1. Thus, for each x G R1 we must have that w(x) = w(x — z) for all z G R1, (up to a set of points of zero length), which implies that w(x) = constant. The constant is a scale factor that may be taken to be equal to 1, and therefore the translation invariant element of measure for points in R1 is the length element,
M (dx) = dx .	(4)
in integral geometry, a motion invariant measure element such as the preceding one is usually called a motion invariant density. Any such density is always taken in absolute value because it must be non negative.
Similarly, the translation invariant density of a point in Rd is the d-dimensional volume element (namely the Lebesgue measure element) in Rd , see Fig. 1a and Fig. 2c. Hence, in this case M(Y) is the volume of Y, which does clearly not depend on the location and orientation of Y.
A test probe equipped with the motion invariant density is called an invariant probe, for short. Henceforth the pertinent invariant densities are given without proof, which can be found for instance in the aforementioned books.
For coherence with the sequel a point of abscissa p on an arbitrary axis (i.e. on R1) may be denoted by L0 = L0(p). Its translation invariant density is
dL0 = dp ,	(5)
as we have seen. From a statistical viewpoint, the probability that a point endowed with the preceding
invariant density belongs to an interval of length dx, given that the point x belongs to an interval of length h > 0, is equal to the ratio of the corresponding measures, namely dx/h. This is the probability element of a uniform random (uR) variable in an interval of length h > 0. Therefore, a point endowed with the invariant density, and belonging to any bounded interval of the real axis, is uR in that interval.
Fix a rectangular frame with origin O in the plane R2. An axis L2[0] of direction 0 G [0, n) is a unoriented straight line through O (hence the subscript '[0]') making an angle 0 with the positive half axis of abscissas. This is equivalent to joining O with a point 0 of the unit semicircle S+. The rotation invariant density of an axis is,
dL2[0] = d0 ,	(6)
namely the length element in S+, see Fig. 1b. Thus, 0 is uniform random (uR) in any interval of the semicircle. In geometric sampling, an axis equipped with the rotation invariant density is said to be isotropic random (IR).
A straight line l2 = l2(p, 0) in R2 is normal to an axis l2[0] of direction 0, called the orthogonal complement of the line, at a signed distance p G (—<x, <x) from O. The pair (p, 0) are called the normal coordinates of L2, see Fig. 1c. The motion invariant density of Lj, namely the density invariant with respect to rotations and translations, is,
dL? = d p d0.	(7)
It is equivalent to take p G [0, and 0 G [0,2n).
ju fr
o	o
® ®
Fig. 1. (a) The translation invariant density for a point p on the real axis is the length element dp. (b) The rotation invariant density for a point 0 on the unit semicircle, or equivalently for an axial direction in the plane, is the arc element d0. (c) A straight line in the plane with its normal coordinates (p, 0).
148
Image Anal Stereol 2015;34:199-208
Because the orientation of an invariant line is IR, and its translation parameter p is UR in any bounded interval of the orthogonal complement of the line, an invariant test line hitting a bounded subset of the plane is said to be isotropic uniform random (IUR) hitting the subset (Miles and Davy, 1976). The latter term applies to any invariant probe hitting a target subset in any dimension.
Fix a rectangular trihedron with origin O in space R3. An axis L^ = ¿1(0, u) of direction u = u(0, d) e
S+ is a unoriented straight line joining O with a point u of the unit hemisphere S+, see Fig. 2a. The angles (0, d) are the spherical polar coordinates of u, namely the longitude 0 e [0,2n), and the colatitude d e [0, n/2]. The rotation invariant density of an IR axis is,
dL3[0] = du = sin d d0 d6 ,	(8)
namely the area element in S+, see Fig. 2b. Thus, u is UR in any region of the hemisphere.
The rotation invariant density of an IR plane l2[0] = L2(0, u) through O is equal to du, namely the same as the invariant density of its normal axis ¿3(0,u).
The plane ¿2 = ¿3(p, u) is parallel to ¿3(0, u) at a distance p e (—from O, see Fig. 3a. Equivalently, ¿2(p, u) is a translate of the plane ¿3(0, u) by a distance p along the orthogonal complement ¿3(0, u). The motion invariant density of ¿2^
dZ| = dp du .	(9)
A straight line Lj = Lj (z, u) in R3 is a translate of the axis Lj(0, u) to a point z in the orthogonal complement of the line, namely in the perpendicular plane L3(0, u), see Fig. 3b. The motion invariant density of Lj is,
dLi = dz du .
(10)
In turn, if (p, q) denote the Cartesian coordinates of z in l2(0, u), then,
dz = d p dq
(11)
(Fig. 2c), namely the motion invariant density for points in L2(0, u). Thus the point z is UR in any bounded region of the latter plane.
Fig. 2. (a) Axial direction through a point u of the unit hemisphere. (b) The rotation invariant density for a direction u in three dimensional space is the area element du on the unit hemisphere. (c) The motion invariant density for a point z in the plane is the area element dz.
Fig. 3. (a) Parametrization of a plane. (b) Parametrization of a straight line.
THE INVARIATOR DESIGN
CLASSICAL CROFTON INTERSECTION
FORMULAE FOR SURFACE AREA AND
VOLUME
The object of interest is a fixed, bounded, nonvoid subset Y of three dimensional Euclidean space R3, (often said 'in 3D'), with piecewise smooth boundary d Y. The target parameters are the volume V of Y, and the surface area S of d Y.
The geometric probe adopted here to estimate S and V is an IUR test line L3 = L3 (z, u) hitting Y. The classical representations of S and V in terms of the intersection measure of an IUR test line hitting Y are given by the following Crofton formulae (for a cursory derivation see Appendix I):
S = I
n
n l3) dL3
(12)
149
Cruz-Orive LM et al: The invariator design
V = 2nJL(YnL?) dL? , (13)
where the density dL? is given by Eq. 10, and /(•), L(^) denote intersection number and intercept length, respectively, with /(0) = L(0) = 0. The integrals are extended to the following set,
{(z,u) : z e ¿3(0,u), u e S+}, (14)
where L?(0, u) denotes the orthogonal complement of the test line L? (z, u).
The preceding identities underlie the classical isotropic sampling designs. An invariant line L\ is normal to an IR plane L^j = L|(0, u) at an invariant point z in that plane. This principle leads to the fakir probe, see for instance Fig. 1a from Cruz-Orive et al. (2010). It can be shown that an equivalent procedure is to take an invariant plane L2 (not just an IR plane L2[0]), and then an invariant test line
within L2. This leads to the isotropic Cavalieri design with test lines, see for instance Fig. 3 from the latter paper. Both designs therefore require translating test lines in 3D. The invariator principle, however, establishes that an invariant test line L13 can be effectively generated within an isotropic plane l2[0]. This requires translating test lines in 2D, which constitutes a practical advantage.
THE INVARIATOR PRINCIPLE
The ensuing approach, and the mathematical derivations given in the appendixes, are somewhat informal - for a more rigorous treatment see Gual-Arnau and Cruz-Orive (2009), Auneau and Jensen (2010), Gual-Arnau et al. (2010), and Thorisdottir and Kiderlen (2014). The idea was first developed in Cruz-Orive (2005), and it was inspired by a result of Varga (1935, Eq. 10). It is noteworthy that the key density decomposition - see Eq. 16 below - appeared in Jensen and Gundersen (1989, Eq. 5.35), but its potential was not exploited further. Schneider and Weil (2008, p. 285) indicate that the mentioned decomposition is related with a general result of Petkantschin (1936). In any case, the relevance of the invariator lies mainly in the practical ramifications of a simple principle.
Recall that the invariant density dL^ of a test line L\ in 3D is given by Eq. 10 as the product of the invariant density dz of a point z in the plane L2[0j orthogonal to the line, times the invariant density du of that plane. In the invariator context the IR plane L^ is called the 'pivotal plane', because it is a plane through a fixed point O and free to rotate or 'pivot' in space around O.
The density dz may be expressed in polar coordinates (r, m) as follows,
dz = rdrdm r e [0, m e [0,2n) . (15)
Thus, dL31 may be decomposed as follows,
dL3 = r dr dm du = r dL^du,	(16)
where dL12 = drdm is the invariant density of a straight line in the plane (Eq. 7). The preceding decomposition is valid provided that the orthogonality of the geometric elements involved is preserved, that is, the distance r must be measured along an axis perpendicular to L\ within the pivotal plane L|[0j (Fig. 4). Moreover, the angle m is measured in the pivotal plane, and therefore the line L12 must also be contained in the pivotal plane at a distance r from O. The preceding construction is the invariator principle.
Fig. 4. Scheme of the invariator principle.
With the invariator principle the classical Eqs. 12, 13 become the invariator equations,
S = ^JdL3[0]/r•1 {(dYnl2[0]) nL1 j dL2,
(17)
V = hjdL2[0] /r• L{ (Ynl3[0]) nl2} dL1,
(18)
respectively. The outer integral is extended over the unit hemisphere S+, whereas, for each direction u of the pivotal plane L|[0] = L^(0, u), the inner integral is extended to the parameter space
{(r,m) : r e [0,«>), m e [0,2n)} (19)
of the test line Lj = Lj (r, m) in the pivotal plane. Thus, instead of hitting the target object Y with an invariant
150
Image Anal Stereol 2015;34:199-208
test line directly in 3D, as it was the case for Eqs. 12, 13, with the invariator design Y is first hit by a pivotal plane, and then the pivotal section is in turn hit by an invariant test line in the pivotal plane. Note also that the distance factor 'r' appears in the integrand.
Eq. 18 may be written,
V = 2n / a (Y n l3[0]) dL3[0]' (20)
where
a(Y n 4[0]) = / r ■ ¿{(y n Z^]) n L2}dL2 (21)
is a functional depending on the pivotal section Y n l2[0] only - and similarly for Eq. 17.
As explained in the following sections, conditional on a given IR direction u e S+ of the pivotal plane, Eqs. 17, 18 give rise to a variety of different unbiased estimators of S and V. Translated to the present special cases, however, Conjecture 4.1 from Gual-Arnau et a/. (2010) says that, conditional on a given pivotal section, the averages of all possible estimators of V will always boil down to Eq. 21 up to a constant factor, and analogously for S. In other words, if yi(u), y2(u) are the averages of any two unbiased estimators of V defined on a pivotal plane l2(0, u), namely if
2n(u)du = ¿/V2(u)du = V , (22)
then the conjecture implies that (u) = y2(u) for all u e S+, and similarly for S. (Note that, in general, if two integrals coincide, the corresponding integrands do not need to coincide).
In Cruz-Orive (2012) the preceding conjecture was proved on the one hand for the different invariator estimators given below for the volume V of an arbitrarily shaped object Y, and on the other hand for the surfactor estimator of the surface area S of a convex body Y, (Jensen and Gundersen 1987, 1989). (Recall that Y is convex if the straight line segment joining any two points of Y, lies entirely in Y). The general conjecture remains open, however, partly because the examined estimators need not be the only possible ones.
It should be stressed that for a finite number of linear probes generated in a pivotal section, the volume estimators given below are different, even though their means will always coincide for any given pivotal section, and similarly for the surface area estimators. Consequently, the estimators will generally have different variances under similar sampling intensities - for an illustration of this see Cruz-Orive (2008).
INVARIATOR ESTIMATORS OF SURFACE AREA AND VOLUME USING TEST LINES
The estimators given in the next subsection emerge directly from Eqs. 17, 18. A pivotal section is hit by an invariant test line with density dL1 in the pivotal plane, and the pertinent measures (/(•), or L(-), respectively) are weighted a posteriori for each test line, namely multiplied by the distance r from the pivotal point O to the test line.
On the other hand, the estimators given in the second subsection below use the fact that the factor rdL2 in the integrand of the mentioned equations is the invariant density dz of a point in the pivotal plane (Eq. 15). Thus, a UR grid of points may be generated in the pivotal plane and then a test line is drawn through each grid point perpendicular to the axis joining the pivotal point O with the grid point. In this way the factor r is implicit in the UR generation of the grid, which means that each test line is weighted automatically by the sampling design. Such test line, which we denote by l2[ ] = l2[ ] (z) may be called a pivotal test line and, as indicated above, its invariant density
dL2[.](z)= rdL2(r, m)= dz , (23)
is that of the point z. The corresponding grid may be called a pivotal grid, see Cruz-Orive (2009a). Thus, for the a priori weighted design the invariator Eqs. 17, 18 take the following form,
S = 1 f dL2[0] f l{ (dYnl3,]) nL1[.]} dL?[.] ,
(24)
V = / dL3P]/L UY n l3P]) n l2[]} dL1[]-
(25)
The preceding design was the one originally described in Cruz-Orive (2005), and it was applied in Cruz-Orive et a/. (2010), where a computer routine was described to generate a pivotal grid in a disk. A fresh pivotal grid has to be generated for each pivotal section, however, which is tedious to do by hand. This was pointed out in Cruz-Orive (2009b), where the question was raised of whether the pivotal grid could be replaced with a simpler one. The answer is provided by the a posteriori weighting design (next subsection) which, in spite of being the more natural, has apparently escaped attention so far.
151
Cruz-Orive LM et al: The invariator design
TEST LINES WITH A POSTERIORI WEIGHTING
On a pivotal plane Z|(0, u) generate a UR grid of parallel test lines a distance T > 0 apart with IR orientation a e [0, n), namely,
{L2(rk, a), k e Z} ,
(26)
where
rk = (Ux + k)T, a = U2n ,	(27)
denote the signed distance from the pivotal point O to the kth test line of the grid, and the orientation of the normal to the test line, respectively, whereas Ui, U2 are two independent UR numbers in the interval [0, 1) , (Fig. 5). Set,
Ik = I{(dYnl2(0,u)) nL2(rk,a)} , Lk = L{(Ynl2(0, u)) nL2(rk, a)} .
Then, (Appendix II),
y = 2nT £ |rk| Ik ,
keZ
y = nT £ |rk|Lk ,
keZ
are unbiased estimators (UE) of S and V, respectively.
To increase precision two mutually perpendicular IUR grids of parallel test lines may be generated on the pivotal plane, whereby the summations in the right hand sides of the preceding two equations should be replaced with their corresponding averages.
(28)
(29)
(30)
Fig. 5. Illustration of the a posteriori weighting method to estimate S and V using Eqs. 29 and 30 respectively. The parallel test lines are IUR hitting the pivotal section - here they are shown horizontal for convenience. The relevant intersections are marked with white dots, the relevant intercepts with thick line segments. Here I-i = 2, I0 = 4, Ii = 2. The intercept lengths and their distances from a parallel axis through the pivotal point O have to be measured at the specimen scale.
TEST LINES WITH A PRIORI WEIGHTING
On the pivotal plane l2(0, u) generate a UR square grid of test points of size T > 0, namely,
{zij = ((Ui + i)T, (U2 + j)T), i, j e Z} . (31)
Let l2[ ] (zij) represent a pivotal line, namely a test line through the point zi j in the pivotal plane, and normal to the axis joining the pivotal point O with zij, see Fig. 6. Then
L2^] (zij), i, j e z}
(32)
a pivotal grid. Set,
I = ££ i{ (d Y n l2 (0, u)) n L2ih](zij)
ieZ jeZ
L = ££l{(ynL2(0,u)) nLuzj)
ieZ jeZ
(33)
namely the total number of intersections and the total intercept length, respectively, determined by the pivotal grid in the pivotal section.
Fig. 6. Illustration of the a priori weighting method to estimate S and V using Eqs. 34 and 35 respectively. The relevant probe is a pivotal grid of test lines superimposed on a pivotal section. Each test line of the grid is drawn through a vertex of a UR grid of points, and it is perpendicular to the axis joining the pivotal point O and the vertex. The relevant intersections are marked with white dots, the relevant intercepts with thick line segments. Only the latter have been drawn, not the complete test lines. Here I = i2, whereas L is the total intercept length at the specimen scale.
i52
Image Anal Stereol 2015;34:199-208
Then (Appendix III),
5 = 2a/,	(34)
V = aL ,	(35)
where a = T2, are unbiased estimators (UE) of S and V, respectively (Cruz-Orive, 2005).
Recall that the invariant density of a pivotal line is that of a point. Therefore a pivotal line is not an invariant test line in the pivotal plane, but it is so in 3D by virtue of the invariator principle. Thus, not surprisingly the preceding estimators are identical to those corresponding to the isotropic fakir probe (Cruz-Orive, 2013, Section 5).
INVARIATOR ESTIMATORS OF VOLUME ONLY
The estimators given below are not recent - for instance the integrated nucleator was already described in (Jensen, 1991; 1998) and applied in Hansen et al. (2011), whereas the nucleator was described in Gundersen (1988). They are included here, however, because they also emerge from Eq. 18, (for general accounts see Jensen and Rataj, 2008; Auneau and Jensen, 2010; Gual-Arnau et al., 2010).
The ordinary version of the integrated nucleator is defined on a pivotal section, and it is based on the following representation,
V = i/du/ 3 p (z, w) dz, (36)
n JS+ J y nL3(0,w)
(Appendix IV), where p (z, w) denotes the absolute distance from the pivotal point O to a point z of the pivotal section. It is assumed that O e Y .In the next subsection, the corresponding volume estimator based on a sample of points in the pivotal section is called the discretized nucleator.
The initial idea underlying the nucleator (Cruz-Orive, 1987, Appendix B), was based on a ray emanating from a fixed point O directly in 3D space, namely in a direction w e S2. If such ray hits the sampled object Y, then the corresponding intersection will in general consist of say m(w) > 1 separate intercept segments. The distances of the end points of these intercepts from O, arranged in increasing order of magnitude, may be denoted as follows,
{/;_ (w), //+(«); i = 1, 2,..., m(w)} . (37)
Then the direct 3D nucleator estimator stems from the following integral,
V = 3/, L (li+(u) -/3-(«)) du . (38) 3 J s2 ¿=1
If O e Y, then Zi_(u) = 0 for all u. If the object Y is moreover star convex with respect to O e Y, namely if the ray joining O with any point of the boundary d Y is always simply connected, then
V = 1j 1+ (u) du .	(39)
where 1+(u) = 11+(u) is the intercept length determined by a ray in the direction u.
The nucleator version stemming from the invariator Eq. 18, however, is a two stage one. In the first stage an IR pivotal plane is generated through the pivotal point O and then, in the second stage an IR ray emanating from O is generated within the pivotal plane making an IR angle p in [0,2n) with a fixed half axis in the pivotal plane. Thus - assuming for the clarity of exposition that Y is star convex with respect to O e Y - the relevant integral now is,
V = -1- /du/2n 1+ (p; u) dp , (40) 3n Js+ Jo +
(Appendix IV). As opposed to the direct nucleator (Eq. 39), the preceding version may be called the pivotal nucleator.
The two nucleator estimators stemming directly from Eqs. 39 and 40 will be clearly unbiased, but not identical, and they will generally not share the same precision for a fixed sample size of n rays, say. Note that in the direct nucleator the n rays may be sampled according to a pseudosystematic design directly on the unit sphere and they will not be coplanar in general (Gual-Arnau and Cruz-Orive, 2002), whereas in the pivotal nucleator the n rays will be sampled on the unit circle within a pivotal section. In practice, the nucleator has generally been applied as a two stage procedure, whereby a pivotal section (either isotropic, or vertical) is sampled first, and then the rays are sampled in that section (Gundersen, 1988).
THE DISCRETIZED NUCLEATOR
On the pivotal plane L|(0, u) generate a UR square grid of test points of size T > 0 exactly as the one represented by Eq. 31. Let n denote the points of the grid hitting the pivotal section Y n L^(0, u), and renumber these points (in any convenient way) as
153
Cruz-Orive LM et al: The invariator design
{zi,Z2,...,zn}• Let di > 0 denote the distance of the hitting grid point zi from O (Fig. 7a). Then,
V = 2a £ di ,
(41)
i=i
where a = T2, is an UE of V.
THE PIVOTAL SECTION NUCLEATOR
On the pivotal plane ¿3(0, u) generate n systematic rays an angle 2n/n apart, emanating from the pivotal point O. Thus the angles made by the n rays with a fixed half axis in the pivotal plane are,
. 2n
yk = (U + k) •—, k = 1,2,...,n n
where U is a UR number in the interval [0,1).
(42)
Fig. 7. (a) Illustration of the discretized nucleator, Eq. 41, (n = 4 distances are generated). (b) Idem of the pivotal nucleator with a single ray, (n = 1, mi = 2), Eq. 44.
If the k-th ray hits the sampled object Y, then the corresponding intersection will consist of say mk > 1 separate intercepts. The distances of the end points of these intercepts from O, in increasing order of magnitude, are denoted as follows,
{h,i-, 4,i+; i =1,2,...,m}■
(43)
(Fig. 7b). If O e Y, then lk>1_ = 0 for all k. The nucleator estimator defined on a pivotal section Y n J|(0, u) reads
n mk
V
n (&+-1
3 ï
k,i-j
(44)
k=1 i=1
and it is unbiased for V.
If Y is star convex with respect to a pivotal point O £ Y, namely if the intersection between Y and any ray emanating from O is always simply connected for any pivotal section, then
- 4n " 3
V = Yn £ l¡+ 3n k=1
(45)
where lk+ = lk,1+ is the intercept length determined by the kth ray in the pivotal section.
INVARIATOR ESTIMATORS OF SURFACE AREA ONLY
CASE OF A CONVEX OBJECT: THE FLOWER FORMULA
In Cruz-Orive (2005) an invariator estimator of surface area, thereafter known as the flower estimator (Cruz-Orive, 2011; Dvorak and Jensen, 2013), was obtained from Eq. 24 for the special case in which the object Y is a convex body. Clearly, in this case
the intersection number I
d Y n l3[O]
C\L

is
equal to 2 almost surely (namely with probability 1) whenever the pivotal test line ¿2^ (z) hits the pivotal section, and equal to 0 otherwise. Equivalently, the hitting event takes place whenever the point z belongs
to a set HYnJ3 called the support set or 'flower' of the
Y J2[0]
pivotal section with respect to the pivotal point O. The boundary d HK of the flower HK of a convex set K with respect to a fixed point O, is the geometric locus of the intersection between a variable tangent to dK and the normal to this tangent from O, see Fig. 8. In other words, dHk is the pedal curve of dK with respect to O, see for instance Lockwood (1961). Thus, in this case Eq. 24 becomes the flower formula:
S = 4f£[ dz
YnL;
■2[0]
4E u{A[hY nL3[0]
(46)
namely S is 4 times the mean value (over all possible isotropic orientations of the pivotal plane) of the area of the flower of the pivotal section.
154
Image Anal Stereol 2015;34:199-208
Fig. 8. (a) Pivotal section (red boundary) of a convex body Y. (b) Flower of the pivotal section.
Let h(a), (0 < a < 2n) denote the radial function of the flower with respect to an interior pivotal point O E Y. Then the flower formula may be written,
S = 4EA 2f h2(m)dmj .
(47)
On the pivotal plane generate n systematic rays an angle 2n/n apart emanating from O, as in the pivotal nucleator design. Measure the corresponding radial lengths {hk+, k = 1,2,..., n} of the flower, or equivalently the 'Feret rays' of the pivotal section. Then an UE of the surface area of a convex body is,
S
4n vh2
n L hk+ n k=i
(48)
which may be regarded as a companion of Eq. 45 when Y is convex. In Cruz-Orive (2005) the choice n = 4 was recommended. Later, Dvorak and Jensen (2013) discovered certain optimality properties of this choice as far as estimation precision is concerned.
CASE OF A GENERAL OBJECT: THE PEAK-AND-VALLEY FORMULA
For a non convex object Y an adaptation of the flower formula is possible starting from Eq. 17. Let C = d Y n l2(0, u) represent the pivotal trace curve. The main task is to identify the number of intersections I(C n L\(r, m)) in terms of (r, m). For each direction angle m e [0,2n) consider a line Lj^r, m) sweeping the pivotal plane parallel to itself from r = down to r = 0, in search of critical points of the height function restricted to the curve C in the given direction. In general, a critical point may be a local maximum, or minimum, with a tangent, or a local supremum, or infimum, without a tangent. To include all cases a local maximum, or supremum, will be called a peak,
whereas a local minimum, or infimum, will be called a valley. If C has points above the axis L2(0, a), then we assign an index £k = +1 if the kth critical point encountered is a peak, and £k = —1 if it is a valley. The first critical point is necessarily a peak, the second is also a peak if C is not convex. Thus, £1 = £2 = +1. Thereafter the critical point may be a peak, or a valley. Immediately after a peak is met, two new intersections appear whereas, after a valley is met, two intersections are lost. Suppose that m critical points are met successively, and let h1 > h2 > . .. hm > hm+1 = 0 denote their corresponding distances from the axis l2(0, a). Then,
I(C n l2 (r, m )) = <
0,
2 = 2e1, 4 = 2(ei + £2),
2 ¿1 £i,
r > h1, r g (h2, h1] , r g (h3, h2] ,
r e (hm+1, hm] .
(49)
Naturally m and {(ek, hk), k = 1,2,..., m} depend on (u, m) in general. Now Eq. 17 yields (Gual-Arnau and Cruz-Orive, 2015),
1 f	¡-2% m phk ( k \
S = - du / dm Y	2 Y eAr dr
n A+ J0 k~í-/hk+^ i=1 )
1	2n m	k
= 1/du/ dm V (h2 - h2+1) V ei (50)
n Js+ J0
i=1
k=1
1 ç	r2n m
- 2dU dm L £kh2 .
n A+ Jo =
From the preceding result a UE of S - called the peak-and-valley estimator - follows, namely,
5 = 4n £ £kh2
(51)
k=1
(with S = 0 if m = 0) which is based on a single IR direction angle a e [0,2n) of the sweeping line in a IR pivotal plane of direction u E S+. To increase precision a number n of systematic orientations may be sampled on the pivotal plane, whereby the right hand side of Eq. 51 should be replaced with the corresponding average. If the target object Y is convex and O E Y, then the preceding estimator reduces to Eq. 48. Note, however, that the general estimator is also valid for
OE Y.
A basically equivalent result was obtained by Thorisdottir and Kiderlen (2014) using a different route. The corresponding estimator, given by Eq. 8 from Thorisdottir et al. (2014), can be shown to simplify into the direct estimator given by Eq. 51.
155
Cruz-Orive LM et al: The invariator design
Example. Fig. 9a represents a pivotal trace C = dY n L2(0, u) determined in the boundary of the target object Y by a pivotal plane l3(0, u) through a fixed pivotal point O previously identifiable in Y, (e.g. a nucleolus of a neuron). The axis Lj(0, m) in the pivotal plane has been conveniently oriented as horizontal, but it is supposed to be isotropically oriented about O. In Fig. 9b the section is the same but, for the sake of illustration, the pivotal point has a different location relative to Y. A sweeping line L2(r, m) moves parallel to Lj (0, m) from top to bottom. Instead of considering that r varies from to 0 for each m e [0,2n), which was convenient to derive the result (50), in practice it is convenient to consider that r varies from to —~ for each m e [0, n). It then suffices to replace the factor 4 with 2 in the right hand side of Eq. 51. Thus the line L2(r, m) sweeps the trace C entirely from top to bottom for each sampled orientation angle m e [0, n). Below the axis L2(0, m), however, it is useful to imagine the pivotal trace rotated by an angle of 180o and use the same criterion to identify peaks and valleys. In this manner m = 4 critical points are met in each case. In Fig. 9a the third critical point is a valley, whereby £3 = — 1. The remaining three critical points are peaks, hence £1 = £2 = £4 = +1. Thus, in this case an unbiased estimate of S would be,
S = 2n(h? + ¿2 - ¿2 + ¿2) .
(52)
In Fig. 9b, however, the four critical points are all peaks, whereby, in this case,
S = 2n(h? + ¿2 + ¿2 + ¿2) .
(53)
Fig. 9. Illustration of the peak-and-valley method to estimate S from a pivotal section. See text.
DISCUSSION
As mentioned in the Section The invariator design, the various estimators described here for volume
coincide in the mean for each pivotal section. For a finite sample of test lines, however, the estimators will have different variances in general. Apart from the empirical study carried out in Cruz-Orive et al. (2010) and the synthetic ones of Cruz-Orive (2008) and Dvorak and Jensen (2013), little is known about the variance of the different estimators. Moreover, variance prediction formulae from a single sample will be precluded in practice whenever a single pivotal section is available. There is therefore scope for further research in this area. The same can be said for the invariator estimator of surface area, which coincides in the mean (for each pivotal section) with the surfactor at least for convex objects (Cruz-Orive, 2012).
The distinction between a priori and a posteriori weighting (Sections The invariator design, Invariator estimators of surface area and volume using test lines) deserves further comments.
(a)	An a priori weighted test line ¿2^ in a pivotal plane is effectively a motion invariant test line in 3D. On the contrary, a posteriori weighting uses a test line L2 which is motion invariant in the pivotal plane, but not in 3D. For this reason, a factor 'r' has to be inserted in the integrand of Eqs. 17, 18 - this means "a posteriori weighting".
(b)	Because a weighted test line L2^ is effectively motion invariant in 3D, it will sample intercept lengths directly from the corresponding intercept length distribution of a given object. Since the latter distribution depends only on the object, the variance of the estimator given by Eq. 35 will not depend on the location of the pivotal point (Cruz-Orive, 2008). On the contrary the variance of the nucleator, for instance, will depend on that location. For an illustration of this fact see Fig. 7 from the latter paper.
Regarding the peak-and-valley formula to estimate surface area, it is noteworthy that the critical points of the pivotal trace do not need to have a tangent. Since we have supposed that d Y is a piecewise smooth surface, the pivotal trace C might be a piecewise smooth curve, and therefore the critical points of its height function might not be defined for every direction. For instance, the third critical point in Fig. 9a is a valley with no tangent, but this is of no consequence for the estimator given by Eq. 51.
In the applications it should be borne in mind that the invariator is a local isotropic design - it is based on an isotropic plane through a fixed point. Among the tools described here, only the nucleator (Eq. 38) can be implemented on local vertical sections.
156
Image Anal Stereol 2015;34:199-208
ACKNOWLEDGEMENTS
The second author acknowledges financial
support from the UJI project P11B2012-24 and the
PROMETEOII/2014/062 project.
REFERENCES
Auneau J, Jensen EBV (2010). Expressing intrinsic volumes as rotational integrals. Adv Appl Math 45:1-11.
Cruz-Orive LM (1987). Particle number can be estimated using a disector of unknown thickness: the selector. J Microsc 145:121-142.
Cruz-Orive LM (2002). Stereology: meeting point of integral geometry, probability, and statistics, Mathematicae Notae. In memory of Professor Luis A. Santalo (1911-2001). 41:49-98.
Cruz-Orive LM (2005). A new stereological principle for test lines in 3D. J Microsc 219:18-28.
Cruz-Orive LM (2008). Comparative precision of the pivotal estimators of particle size. Image Anal Stereol 27:1722.
Cruz-Orive LM (2009a). The pivotal tessellation. Image Anal Stereol 28:101-05.
Cruz-Orive LM (2009b). Stereology: old and new, In: Capasso V, Aletti G, Micheletti A, eds. Proc 10th Eur Congr Stereol Image Anal. Bologna: Esculapio.
Cruz-Orive LM (2011). Flowers and wedges for the stereology of particles. J Microsc 243:86-102.
Cruz-Orive LM (2012). Uniqueness properties of the invariator, leading to simple computations. Image Anal Stereol 31:87-96.
Cruz-Orive LM (2013). Variance predictors for isotropic geometric sampling, with applications in forestry. Stat Methods Appl 22:3-31.
Cruz-Orive LM, Ramos-Herrera ML, Artacho-Perula E (2010). Stereology of isolated objects with the invariator. J Microsc 240:94-110.
De-Lin R (1994). Topics in integral geometry. Singapore: World Scientific.
Dvorak J, Jensen EBV (2013). On semiautomatic estimation of surface area. J Microsc 250:142-57.
Gual-Arnau X, Cruz-Orive LM (2002). Variance prediction for pseudosystematic sampling on the sphere. Adv Appl Prob 34:469-83.
Gual-Arnau X, Cruz-Orive LM (2009). A new expression for the density of totally geodesic submanifolds in space forms, with stereological applications. Diff Geom Appl 27:124-28.
Gual-Arnau X, Cruz-Orive LM (2015). New rotational integrals in space forms, with an application to surface area estimation. Appl Math-Czech. To appear.
Gual-Arnau X, Cruz-Orive LM, Nuno-Ballesteros JJ (2010). A new rotational integral formula for intrinsic volumes in space forms. Adv Appl Math 44:298-308.
Gundersen HJG (1988). The nucleator. J Microsc 151:3-21.
Hansen LV, Nyengaard JR, Andersen JB, Jensen EBV (2011). The semi-automatic nucleator. J Microsc 242:206-15.
Jensen EBV (1991). Recent developments in the stereological analysis of particles. Ann Inst Statist Math 43:455-68.
Jensen EBV(1998). Local stereology. Singapore: World Scientific.
Jensen EBV, Gundersen HJG (1987). Stereological estimation of surface area of arbitrary particles. Acta Stereol 6 (Suppl. III):25-30.
Jensen EBV, Gundersen HJG (1989). Fundamental stereological formulae based on isotropically orientated probes through fixed points with applications to particle analysis. J Microsc 153:249-67.
Jensen EBV, Rataj J (2008). A rotational integral formula for intrinsic volumes. Adv Appl Math 41:530-60.
Lockwood EH (1961). A book of curves. Cambridge: Cambridge University Press.
Miles RE, Davy PJ (1976). Precise and general conditions for the validity of a comprehensive set of stereological fundamental formulae. J Microsc 107:211-26.
Petkantschin P (1936). Zusammenhange zwischen den Dichten der linearen Unterräume im n-dimensionalen Raum. Abh Math Sem Hamburg 11:249-310.
Santalo LA (1976). Integral geometry and geometric probability. Addison-Wesley, Reading, Massachusetts.
Schneider R, Weil W (2008). Stochastic and integral geometry. Berlin: Springer-Verlag.
Thorisdottir O, Kiderlen M (2014). The invariator principle in convex geometry. Adv Appl Math 58:63-87.
Thorisdottir O, Rafati AH, Kiderlen M (2014). Estimating the surface area of nonconvex particles from central planar sections. J Microsc 255:49-64.
Varga O (1935). Integralgeometrie 3. Croftons Formeln fur den Raum. Math Z 40:387-405.
APPENDIX I: CROFTON
INTERSECTION FORMULAE
The volume element determined in the object Y by a test line L\(z,u) of fixed orientation u e S+ is dv(z, u) = L(Y n Lj(z, u)) dz. Integration over the domain given by Eq. 14 yields Eq. 13.
To obtain Eq. 12 consider a test line L\ (z, u) hitting
157
Cruz-Orive LM et al: The invariator design
an essentially planar surface element dy c d Y of area ds. Thus, I(dy n L?) = 1 if dy n L? = 0, and I(dy n L? ) = 0 if dy n L? = 0. Because the density dLl is motion invariant, the polar axis Ox? may be taken along the normal to the surface element, whereby dz = ds| cos d |, see Fig. 10, and
dLl = dz du = ds | cos d | du
(54)
Now,
JI(dY nLl)dLl = i	I(dy n Ll) dL?
Id Y JdynLl=0 i i
Id YJ d ynL3=0
dL?
Jd Y
n S,
r 2n pn/2
ds / d<j) (
J0 J0
cos d sin d d0
(55)
which verifies Eq. 12.
Fig. 10. Geometric elements to obtain the motion invariant density of a straight line hitting a surface element.
APPENDIX II: ESTIMATORS OF S, V USING A POSTERIORI WEIGHTED TEST LINES
Write I(r, rn, u) = I{(dYnL?(0, u)) nL?(r,«)}, for short. Eq. 17 may be expressed as follows,
r	j-n	/•+»
nS = du l dm |r| I(r, a, u) dr
J S+ J0
du f dm ^
c(k+1)T
10
ke Z'
'kT
| r | I(r, a, u) dr
du I dm i ^ |r + kT| I(r + kT, a, u) dr.
keZ
(56)
According to the design described in the Subsection Test lines with a posteriori weighting, the variables (u, ®,r) are independent UR in S+, [0,2n), [0, T), respectively. Therefore, their joint probability element (namely their joint probability density function times dud®dr) reads,
du da dr
P(du, dm, dr) = —— - .
2n n 1
(57)
Returning to the notation introduced in the Subsection Test lines with a posteriori weighting, it follows that
nS = 2n2T ■ Eu,®,r £ |r*|Ik ,	(58)
ke Z
which verifies the unbiasedness of the estimator given by Eq. 29. The proof for Eq. 30 is analogous.
APPENDIX III: ESTIMATORS OF S, V USING A PRIORI WEIGHTED TEST LINES
The identities in Eq. 56 were obtained by partitioning the real axis into congruent, non overlapping segments of length T, as R = UkeZ [kT, (k + 1)T). In the present case it is convenient to partition the plane into non overlapping tiles congruent with a fundamental tile J0, namely,
R2 = U Jk, Jk = J0 + Tk, Jk n Jl = 0
keZ
if k = l, k, l e Z. (59)
Every tile Jk can be brought to coincide with the fundamental tile J0 by means of a translation — Tk which leaves the partition unchanged. Here J0 c l3(0, u) is adopted to be a square of area a. Write
I(z, u)= I {(d Y n l2(0, u)) n L21h] (z^, for short. Now Eq. 24 may be written as follows,
n S
IS+
du I(z, u) dz
J R2
du ^ I(z, u) dz
keZJ Jk
du ^ I(z + Tk, u) dz
Jj0 keZ
(60)
158
Image Anal Stereol 2015;34:199-208
The variables (u, z) are independent UR in S+, Jo, respectively, so that,
_ , , . du dz
P(du, dz) = ---.
2n a
(61)
Returning to the notation introduced in the Subsection
Test lines with a priori weighting,
n S = 2n a ■ Eu,z (I )
(62)
which verifies the unbiasedness of the estimator given by Eq. 34. The proof for Eq. 35 is analogous.
The technique used in this, and in the preceding subsection, is a direct application of a theorem given in Ch.8 of Santalo (1976); see also Cruz-Orive (2002), (Eq. 4.6).
APPENDIX IV: THE INTEGRATED NUCLEATOR, AND THE PIVOTAL NUCLEATOR, AS DIRECT CONSEQUENCES OF THE INVARIATOR EQUATION FOR VOLUME
Write L(r, a, u) = L {(Y n L|(0, u)) n L?(r,a)}, for short. Eq. 18 may be written as follows,
2nV
j- j-n /•+»
/ du/ d® / |r| L(r, ®, u) dr . (63)
J s2 Jo J-™
The length of the chord (Y n L3,(0, u)) n L?(r, a) = C, say, (Fig. 11a), may be expressed as the integral of a length element dl along the extent of C, so that,
/+"	rr
|r|L(r, a, u) dr = |r| dr dl. (64)
J—^	JC
The length elements dr and dl are orthogonal, and therefore dr dl = dz is the area element at a point of the pivotal section, namely of z e Y n L\(0, u), see Fig. 11b. Let (p, <p) denote the polar coordinates of the point z. Then |r| = p | cos(a — <p) |, whereby,
2nV
Í du í d® Í p | cos(® - ç)| dz
J S+ Jo J y n¿3(o,u)
du
S+ J Y nL|(0,u)
p(z,u)dz J |cos(® - ç)| d®
= 2/ du/ 3 p (z, u) dz ,
J S+ JYnL3(o,u)
which is the identity Eq. 36.
(65)
Fig. 11. Geometric elements involved in Appendix IV, see text.
The verification of the unbiasedness of the discretized nucleator Eq. 41 is immediate using Santalo's theorem as in the preceding appendix.
To obtain the pivotal nucleator identity (40) for a star convex set with respect to O, replace dz with its polar coordinate version dz = p dp da in the last Eq. 65, whereby,
2nV = 2
r 2n r l+(ç ,u)
du I dç
Jo Jo
P 2 dp
2	r	r 2n
ö/2du/ l+ (ç,u)dç .
3	Js+ Jo +
(66)
APPENDIX V: EQUIVALENCE BETWEEN THE DIRECT AND THE PIVOTAL NUCLEATOR REPRESENTATIONS
Let (p, u), p > 0, u e S2 denote the spherical polar coordinates of a point x e R3. The corresponding volume element is dx = p2dpdu. The volume of a star convex set Y with respect to an interior origin O e Y may then be expressed as follows,
r rl+(u) 0 V = 2 du p2 dp
J s2 Jo
= - l+ (u) du ,
3 ys2+v ; '
(67)
which is the volume representation corresponding to the direct nucleator (Eq. 39). Therefore the right hand sides of Eqs. 39 and 40 must coincide. As pointed out in the Section Invariator estimators of volume only, the corresponding estimators are unbiased but different from one another in general.
159