THE TOMATO SALAD PROBLEM IN SPHERICAL STEREOLOGY^{1}

R. Gorenflo

Abstract

In 1925 S.D. Wicksell studied the problem of determining the probability density of the sizes of spherical particles embedded in a solid opaque medium from that of the sizes of their intersections with a random plane cut. Mathematical modelling leads to a first kind Abel type integral equation. Beginning in 1958 G. Bach published a series of papers on a modified problem, replacing the cut by a slice of positive thickness bounded by parallel planes and measuring the radii of circles obtained by orthogonal projection of the cut-out parts of spheres on a plane of observation lying parallel to the slice. P.L. Goldsmith in 1967 presented a solution formula to this ”tomato salad problem”.

We here give a description of the problem and another derivation of a solution formula, using the formalism of fractional calculus, hoping this new derivation to be more transparent (at least to people familiar with fractional calculus) than Goldsmith’s. Furthermore, we discuss aspects that, to a large extent, have been neglected by previous authors, namely the transition from the well-posed tomato salad problem to Wicksell’s ill-posed problem (as the thickness of the slice tends to zero), and non-negative measures instead of classical functions as densities.

1991 Mathematics Subject Classification: 45E10, 26A33, 65R30.

Key Words and Phrases: Applications of Abel integral equations in stereology of spherical particles, Weyl transform, transition from well-posedness to ill-posedness.

1Lecture presented at the Second International Workshop on Transform Methods & Special Functions, Varna/Bulgaria, 23-30 August 1996. To be published in the Proceedings of this workshop via Science Culture Technology Publishing (SCTP) Singapore.

1 Introduction

An important problem of stereology is the determination of the probability density of the sizes of randomly distributed particles in a solid opaque medium from that of the sizes of their intersections with a random plane cut.

Applications are, among others, the study of corpuscles embedded in the (solidified) tissue of a biological organ, in sedimentary petrography, in crystallography, in investigation of carbide particles in steel, in statistics of sizes of bubbles in Swiss cheese. See references given in [9] and literature on stochastic geometry, e.g. [3], [12], [16], [17]. In spherical stereology the case of spherical particles is considered. The pioneering paper in this field is that by Wicksell [18], published in 1925. As one of the most lucid papers we recommend [15].

The tomato salad problem has been treated by Bach in a series of pioneering papers of which we quote [1], [2]. A solution formula has been derived by Goldsmith [5]. A random slice with positive thickness t, bounded by two parallel planes, is cut out from the solid opaque medium. The thickness t should be so small that the slice (apart from parts of spheres contained in it) is still transparent. The cut-out parts of spheres are orthogonally projected on another plane (lying parallel to the slice), where they produce circles as shadows whose radii are observable and are called ”apparent” radii in distinction from the ”actual” radii of the spheres. A simplifying assumption is that the slice is so thin or the spatial density of spheres is so small that overlapping of shadows is negligible. The problem is to determine the probability distribution of the actual radii from the observed one of the apparent radii.

Naturally, t > 0, but we formally require t ≥ 0 so that Wicksell’s classical problem of spherical stereology (described in the Abstract) is a special case. Whereas this classical problem is modelled by an Abel integral equation of first kind (ill-posed in customary function spaces) the case t > 0 leads to a well-posed Abel integral equation of second kind.

We imagine our slice cut out from the solid opaque medium as horizontal. This means that we can (and henceforth shall) speak of the ”upper bounding plane” and the ”lower bounding plane” (in case t > 0). The third plane, also horizontal, is the ”observation plane”. We imagine it as lying below the lower bounding plane. The image of a cut-out part of a sphere, produced by the orthogonal projection on the plane of observation, is called the ”shadow” of this sphere. The reader is advised to make sketches of the geometrical configurations described in words. See Fig. 1.

We concentrate in this paper our attention to analytical aspects. The reader who is interested in numerical methods may look into [5], [6], [7], [8], [10] and [13].

[figure: r; 2x; t]

Fig. 1: General configuration

2 Notation and generalities

Let us first introduce some notations. By λ we denote the spatial density of spheres with centers in a unit volume of the medium, i.e. the expected number of spheres with centers in a unit volume. Now, let

r_{2} > r_{1} ≥ 0, x_{2} > x_{1} ≥ 0, τ_{2} > τ_{1} ≥ 0, σ_{2} > σ_{1} ≥ 0,

and understand ∫_{α}^{β} for β > α to mean ∫_{α−}^{β−}. Then we define non-negative generalized functions f̃, f, g̃, g, p̃, p, q̃, q by the following properties.

∫_{r_{1}}^{r_{2}} f̃(r)dr = expected number of spheres, with centers in a unit volume,

having radii in the semi-open interval [r_{1}, r_{2}).

∫_{x_{1}}^{x_{2}} g̃(x)dx = expected number of spheres meeting the slice

whose centers are projected on a unit area of the observation plane and whose apparent radii are in [x_{1}, x_{2}).

∫_{r_{1}}^{r_{2}} f(r)dr = probability of actual radius in [r_{1}, r_{2}).
∫_{x_{1}}^{x_{2}} g(x)dx = probability of a sphere that meets the slice

having apparent radius in [x_{1}, x_{2}).

m = ∫_{0}^{∞} rf(r)dr = mean actual radius of a sphere.  (2.1)

Remark: We shall always assume m to be finite.

∫_{τ_{1}}^{τ_{2}} p̃(τ)dτ = expected number of spheres, with centers in a unit volume,

having central cross-section area in [τ_{1}, τ_{2}).

Explanation: A sphere with radius r has central cross-section area πr^{2}.

∫_{σ_{1}}^{σ_{2}} q̃(σ)dσ = expected number of spheres meeting the slice

whose centers are projected on a unit area of the observation plane and whose shadows have an area in [σ_{1}, σ_{2}).

∫_{τ_{1}}^{τ_{2}} p(τ)dτ = probability of a sphere having central cross-section area

in [τ_{1}, τ_{2}).

∫_{σ_{1}}^{σ_{2}} q(σ)dσ = probability of a sphere that meets the slice

having shadow area (of its cut-out part) in [σ_{1}, σ_{2}).

Remark: By admitting generalized functions (actually measures) we can treat by the same formalism continuous as well as more general probability distributions. By our understanding of definite integrals the primitives

F (r) = ∫_{0}^{r} f(s)ds, G(x) = ∫_{0}^{x} g(s)ds, P (τ) = ∫_{0}^{τ} p(s)ds, Q(σ) = ∫_{0}^{σ} q(s)ds

for nonnegative arguments are nonnegative non-decreasing left-continuous functions tending to 1 for their respective argument tending to infinity. If, for example, f(r) = δ(r − a) with a ≥ 0, then F (r) = 0 for r ≤ a, F (r) = 1 for r > a. A reader who does not like integration of generalized functions may rewrite this paper replacing f(τ)dτ by dF (τ) etc.

We note the basic relations between our functions and their arguments. With r as actual and x as apparent radius of a sphere, τ as central cross-section area and σ as shadow area we have

τ = πr^{2}, σ = πx^{2}, r = \sqrt{τ/π}, x = \sqrt{σ/π},  (2.2)

and the equalities of differentials

f(r)dr = p(τ)dτ, g(x)dx = q(σ)dσ,
f̃(r)dr = p̃(τ)dτ, g̃(x)dx = q̃(σ)dσ,
dτ = 2πrdr, dσ = 2πxdx

imply

\begin{matrix}p(τ) = f(r)/(2πr), & q(σ) = g(x)/(2πx), \\ p̃(τ) = f̃(r)/(2πr), & q̃(σ) = g̃(x)/(2πx). \end{matrix} \}  (2.3)

Further, we observe the relations

f̃(r) = λf(r),  (2.4)
g̃(x) = λ(t + 2m)g(x).  (2.5)

The validity of (2.4) is obvious, and (2.5) follows from the fact that the total expected number ∫_{0}^{∞} g̃(x)dx of spheres that meet the slice and whose centers are projected on a unit area of the observation plane is λ(m + t + m), with m as mean radius of a sphere. Here t comes from the sphere centers inside the slice whereas m + m comes from centers above the upper and below the lower bounding plane. The requirement ∫_{0}^{∞} g(x)dx = 1 and the natural scaling relation

g(x) = (∫_{0}^{∞} g̃(s)ds)^{−1} g̃(x)  (2.6)

also yield (2.5). Note that scaling relations analogous to (2.6) hold for f and f̃, q and q̃, p and p̃.

Another derivation of (2.5) will be given in the next section in which we shall derive two integral equations: Bach’s equation for determining f̃from g̃, Goldsmith’s equation for determining f from g.

We shall use the notion of characteristic function of a set. For M ⊂ R we define

χ_{M}(s) = 1 if s ∈ M, = 0 if s ∉ M.  (2.7)

In Sections 4 and 5 an essential role will be played by the Mittag-Leffler function with index 1/2. On general properties of Mittag-Leffler functions the reader can find orientation in Chapter 18 of [4]. In particular,

E_{1/2}(z) = ∑_{n=0}^{∞} \frac{z^{n}}{Γ (\frac{n}{2} + 1)} for z ∈ C.

We need this function for non-positive real argument. It is known that

E_{1/2}(−s) = \frac{2}{\sqrt{π}} exp(s^{2}) ∫_{s}^{∞} exp(−x^{2})dx for s ∈ R.  (2.8)

For completeness we present a short nice proof (adapted from [5]) of (2.8). Indeed, for all s ∈ C, we have E_{1/2}(−s) = u(s) − v(s) with

u(s) = ∑_{k=0}^{∞} \frac{s^{2k}}{Γ(k + 1)} = exp(s^{2}), v(s) = ∑_{k=0}^{∞} \frac{s^{2k+1}}{Γ (k + \frac{3}{2})}.

With y = \sqrt{2}s we get, using Γ(1/2) = \sqrt{π} and the functional equation of the gamma function, v(s) = \sqrt{2/π}w(y) with

w(y) = y + \frac{1}{3}y^{3} + \frac{1}{3 × 5}y^{5} + \frac{1}{3 × 5 × 7}y^{7} + …,

and readily we see that

w′(y) = yw(y) + 1 for y ∈ C, w(0) = 0.

This initial value problem is uniquely solved by

w(y) = exp (\frac{y^{2}}{2} ) ∫_{0}^{y} exp (−\frac{z^{2}}{2} ) dz.

Back substitution x = z/\sqrt{2} gives

v(s) = \frac{2}{\sqrt{π}} exp(s^{2}) ∫_{0}^{s} exp(−x^{2})dx,

E_{1/2}(−s) = exp(s^{2}) \{1 − \frac{2}{\sqrt{π}} ∫_{0}^{s} exp(−x^{2})dx\} .  (2.9)

Using the well-known identity 2 ∫_{0}^{∞} exp(−x^{2})dx = \sqrt{π} we now deduce (2.8) from (2.9).

It is not difficult to show that E_{1/2}(−s) is strictly decreasing from 1 towards 0 as s runs from 0 towards ∞ (through real s).

3 Derivation of integral equations

In this section we assume t > 0 if not explicitly stated otherwise. First we consider the special case of all spheres having the same actual radius ρ > 0. Then, with Dirac’s δ-function, f(r) = δ(r − ρ) =: f(r, ρ). Let us use the notations g(x, ρ) and g̃(x, ρ) for the corresponding functions g(x) and g̃(x), respectively.

The expected number of spheres with centers in the slice and projected center lying in a unit observational area is tλ. For shortness of expression we shall henceforth say ”unit observational area” for ”unit area of the observation plane”. The expected number of spheres with centers outside the slice but meeting it and having projected center in a unit observational area is 2ρλ, the factor 2ρ stemming from the fact that intersection happens if the center of the sphere is at most distant by ρ from the slice but can lie above the upper bounding plane or below the lower bounding plane. See Fig. 2.

[figure: t; ρ; x; w]

Fig. 2: A sphere immersed from above

Denoting by w the depth of immersion below the upper bounding plane (if the center is above this plane) or above the lower bounding plane (if the center is below this plane) we see, by a theorem of elementary geometry, that for the apparent radius x we have

x^{2} = w(2ρ − w).  (3.1)

The random variable w being uniformly distributed in the interval [0, ρ] the probability for immersion depth lying in the interval [w, w + dw) is \frac{dw}{ρ} for 0 ≤ w < ρ and positive differential dw. From (3.1) we get 2xdx = (2ρ − 2w)dw, hence dw = \frac{xdx}{ρ−w}, and we conclude that the expected number of such spheres that produce apparent radii in [x, x + dx) is

2λdw = \frac{2λxdx}{ρ − w} for 0 ≤ w < ρ,

the factor 2 resulting from contributions from above and below the slice. From (3.1) or by invoking Pythagoras we obtain ρ − w = (ρ^{2} − x^{2})^{1/2}, hence

2λdw = \frac{2λxdx}{\sqrt{ρ^{2} − x^{2}}} for 0 ≤ x < ρ.

Remembering now in addition the expected number λt of sphere centers inside the slice and projected on a unit observational area (these spheres produce an apparent radius equal to the actual radius ρ) we find (using the notation (2.7)), for x ≥ 0,

g̃(x, ρ) = λtδ(x − ρ) + \frac{2λx}{\sqrt{ρ^{2} − x^{2}}} χ_{[0,ρ)}(x).  (3.2)

Using the scaling relation (2.6) and taking into account the expected number

λ(t + 2ρ) = ∫_{0}^{∞} g̃(x, ρ)dx

of sphere centers inside the slice or at distance at most ρ from it that are projected on a unit observational area we arrive at

g(x, ρ) = \frac{1}{t + 2ρ} \{tδ(x − ρ) + \frac{2x}{\sqrt{ρ^{2} − x^{2}}} χ_{[0,ρ)}(x)\} .  (3.3)

Remark: It is known that very small spheres are difficult to detect in Wicksell’s classical arrangement (t = 0). Let us look what happens with equation (3.2) in the limiting cases t = 0, ρ = 0. Observe that we have defined (3.2) and (3.3) under the assumption (i) t > 0, ρ > 0. Consider the three remaining cases (ii) t > 0, ρ = 0, (iii) t = 0, ρ > 0, (iv) t = 0, ρ = 0. A critical re-view shows that (3.2) becomes

g̃(x, 0) = λtδ(x) in case (ii),

g̃(x, ρ) = \frac{2λx}{\sqrt{ρ^{2} − x^{2}}} χ_{[0,ρ)}(x) in case (iii),

g̃(x, 0) ≡ 0 in case (iv).

This last relation reveals that point-like particles are practically undetectable in Wicksell’s arrangement (they are hit by the plane cut with probability zero). We shall return to this question in Section 6.

Now we have the situation that there is only one possible actual radius. Generalizing our treatment and working with positive differentials dx and dr we see (a) and (b).

(a) The expected number of spheres centered inside the slice above a unit observational area, producing apparent radii in [x, x + dx) is λtf(x)dx (for such spheres actual and apparent radii coincide).

(b) The expected number of spheres meeting the slice but having their centers outside of it above a unit observational area and having actual radius in [r, r + dr) that produce apparent radius in [x, x + dx) is

\frac{2λx}{\sqrt{r^{2} − x^{2}}} χ_{[0,r)}(x)f(r)drdx.

Integrating this from r = 0 to r = ∞ gives the expected number of spheres meeting the slice but centered outside of it above a unit observational area as

2λx ∫_{r=x}^{∞} \frac{f(r)dr}{\sqrt{r^{2} − x^{2}}} dx.

Adding the contributions (a) and (b) and taking account of (2.4) f̃ = λf gives us Bach’s integral equation connecting the densities g̃ and f̃ of expected numbers

g̃(x) = tf̃(x) + 2x ∫_{r=x}^{∞} \frac{f̃(r)dr}{\sqrt{r^{2} − x^{2}}} .  (3.4)

For given g̃ this is an Abel type integral equation of second kind for determining f̃ from g̃.

To find an integral relation combining the probability densities f and g we remember the scaling condition (2.6) g(x) = g̃(x)/μ with μ = ∫_{0}^{∞} g̃(x)dx. Using ∫_{0}^{∞} f(r)dr = 1 and again (2.4) we obtain by integrating (3.4) from x = 0 to x = ∞ after interchanging the order of integrations μ = λ(t + 2m) and thus have another proof of (2.5). Dividing (3.4) by λ we obtain Goldsmith’s integral equation

(t + 2m)g(x) = tf(x) + 2x ∫_{r=x}^{∞} \frac{f(r)dr}{\sqrt{r^{2} − x^{2}}}  (3.5)

which for determination of f from given g is to be complemented by the scaling condition

∫_{0}^{∞} f(r)dr = 1.  (3.6)

We can solve the system {(3.5), (3.6)} for the unknown function f and the unknown mean radius m by the following procedure. Solve (3.5) for ϕ(r) := f(r)/(t + 2m). Then, with γ = ∫_{0}^{∞} ϕ(r)dr, we have f(r) = ϕ(r)/γ, m = (\frac{1}{γ} − t)/2.

4 Solution of Bach’s integral equation

The essential work to be carried out is to solve Bach’s integral equation (3.4), g̃ given, f̃ unknown, or the integral equation resulting from (3.5) by substitution f(r) = (t + 2m)ϕ(r) and division by (t + 2m). Both integral equations have the same form, and so we concentrate on the solution of (3.4),

g̃(x) = tf̃(x) + 2x ∫_{r=x}^{∞} \frac{f̃(r)dr}{\sqrt{r^{2} − x^{2}}} ,

treating the case t > 0 of positive thickness of the slice. At the end of this section we shall make a short remark on Wicksell’s case t = 0.

In analogy to [15] (where only the case t = 0 is treated) and deviating from Goldsmith [5] we take as independent variables the central cross-section area τ = πr^{2} and the shadow area σ = πx^{2} of a sphere instead of its actual and apparent radius r and x, respectively. See Section 2 for the notations. Bach’s equation (3.4) is by this substitution equivalent to the integral equation

q̃(σ) = tp̃(σ) + \frac{1}{\sqrt{π}} ∫_{τ=σ}^{∞} \frac{p̃(τ)dτ}{\sqrt{τ − σ}} ,  (4.1)

whereas Goldsmith’s equation (3.5) with the scaling condition (3.6) is equivalent to

(t + 2m)q(σ) = tp(σ) + \frac{1}{\sqrt{π}} ∫_{τ=σ}^{∞} \frac{p(τ)dτ}{\sqrt{τ − σ}}  (4.2)

with the scaling condition

∫_{0}^{∞} p(τ)dτ = 1.  (4.3)

The system { (4.2), (4.3)} can be treated in analogy to the procedure described at the end of the preceding Section 3. After having calculated p̃(τ) from (4.1) or p(τ) from {(4.2), (4.3)} we get by back-substitution, see (2.3),

f̃(r) = 2πrp̃(πr^{2}), f(r) = 2πrp(πr^{2}).  (4.4)

To solve (4.1) formally we apply the Weyl fractional calculus described in [14] and then discuss conditions for validity of the formal solution obtained. We use Weyl operators V^{α}, α > 0, according to

(V^{α}u)(y) = \frac{1}{Γ(α)} ∫_{s=y}^{∞}(s − y)^{α−1}u(s)ds, y ≥ 0.  (4.5)

This expression is meaningful if the generalized function u(s) behaves sufficiently well, e.g. if it is bounded and tends to zero faster than every negative power of s as s → ∞. Let us remark that Miller and Ross in [14] write W^{−α} instead of V^{α}. With V^{0} := I = identity operator the operators V^{α} form a semigroup:

V^{α+β} = V^{α}V^{β} for α ≥ 0, β ≥ 0.

Writing (4.1) in the form

q̃ = (tI + V^{1/2})p̃:= B_{t}p̃  (4.6)

we obtain formally

p̃ = (tI + V^{1/2})^{−1}q̃ = B_{t}^{−1}q̃

and by Neumann series expansion

B_{t}^{−1} = t^{−1}(I + t^{−1}V^{1/2})^{−1} = t^{−1} ∑_{n=0}^{∞}(−t)^{−n}V^{n/2},

p̃(τ) = t^{−1} \{q̃(τ) + ∑_{n=1}^{∞} \frac{(−t)^{−n}}{Γ(n/2)} ∫_{τ}^{∞}(σ − τ)^{\frac{n}{2} −1}q̃(σ)dσ\} .

If q̃(τ) is so well behaved that the series converges and ∑ and ∫ can be interchanged we obtain

p̃(τ) = \frac{1}{t} \{q̃(τ) + ∫_{τ}^{∞} w_{t}(σ − τ)q̃(σ)dσ\} = (B_{t}^{−1}q̃)(τ)  (4.7)

with

w_{t}(y) = ∑_{n=1}^{∞} \frac{(−t)^{−n}}{Γ(n/2)}y^{\frac{n}{2} −1} = \frac{d}{dy} E_{1/2} (−\frac{\sqrt{y}}{t} ) .  (4.8)

Here E_{1/2} is the Mittag-Leffler function of index 1/2, see Section 2 for details. We see that asymptotically

w_{t}(y) ∼ −\frac{t^{−1}y^{−1/2}}{\sqrt{π}} for y → 0,

and from formulas (1.17), (2.5), (3.4), (4.3) of [9] we deduce

w_{t}(y) ∼ −\frac{ty^{−3/2}}{2\sqrt{π}} for y → ∞.

Hence (4.7) is meaningful (the convolution integral converges) if q̃(τ) is integrable and bounded. Appropriately interpreted, (4.7) is valid as solution formula under much wider conditions.

In situations arising in practice the density p̃ may have the form

p̃(τ) = γ(τ) + ∑_{j∈M} c_{j}δ(τ − a_{j}), τ ≥ 0.

Here γ(τ) is non-negative, piecewise continuous, exponentially decreasing to zero for τ → ∞ (or having bounded support), M is a finite or denumerably infinite index set, all a_{j} > 0, all c_{j} > 0, and ∑_{j∈M} c_{j} < ∞. See [6], [7], [8] and [13].

By careful direct calculation it can be shown that our formalism, in particular formula (4.7), does work correctly if applied to a generalized function q̃corresponding by (4.1) to a generalized function p̃ of this special form.

Remark: In the special case t = 0 we have, with Du = u′, formally

B_{0}^{−1} = V_{1/2}^{−1} = −DV^{1/2}  (4.9)

(see [14]).

5 Noise amplification

Again, assume t > 0 if not stated explicitly otherwise.

In practice, the function q̃(σ) is given by statistical measuring and counting procedures, hence unavoidably numerically not precisely. It is therefore important to investigate the amplification of noise effected in applying the inverse operator B_{t}^{−1} to noisy data. In Banach space setting the worst case amplification factor is given by the operator norm ∥B_{t}^{−1}∥.

If we work in the Lebesgue space L_{β} = L_{β}[0, ∞), β ∈ [1, ∞], of functions defined on [0, ∞) then we observe that the operator B_{t} of (4.6) may be applicable only in part of L_{β}. For example,

∫_{σ}^{∞}(τ − σ)^{−1/2} min(1, τ^{−α})dτ ≡ ∞, if 0 ≤ α ≤ 1/2 and σ ≥ 0,

but the function min(1, τ^{−α}) is in L_{β} if β > 1/α. However, B_{t} is applicable to all functions in L_{β} with bounded support.

We shall now show that the inverse operator B_{t}^{−1} of the solution formula (4.7) is a bounded linear operator from L_{β} into L_{β}, for all β ∈ [1, ∞]. Our argumentation will be analogous to one used in [11].

In fact, by appropriate changes of variables and extension by zero of the involved functions, the integral in (4.7) can be written as a convolution of functions defined on R, and using Young’s inequality for convolutions,

∥u ∗ v∥_{L_{β}(R)} ≤ ∥u∥_{L_{1}(R)}∥v∥_{L_{β}(R)}

we find, with ∥.∥_{β} = ∥.∥_{L_{β}[0, ∞)}, from (4.7)

∥p̃∥_{β} ≤ \frac{1}{t}∥q∥_{β}(1 + ∥w_{t}∥_{1}).

Because w_{t}(x) < 0 for x > 0 we have

∥w_{t}∥_{1} = ∫_{0}^{∞} |w_{t}(x)|dx = − ∫_{0}^{∞} \frac{d}{dx}E_{1/2} (−\frac{\sqrt{x}}{t}) dx = E_{1/2}(0) = 1.

Hence

∥p̃∥_{β} ≤ \frac{2}{t}∥q∥_{β}, ∥B_{t}^{−1}∥ ≤ \frac{2}{t}.  (5.1)

This bound tends to ∞ as t → 0, in accordance with the well-known fact that B_{0} = V^{1/2} is a smoothing operator whose inverse is unbounded.

When we consider B_{t}^{−1} as an operator acting only on functions with support in an interval [0, S], with S ∈ (0, ∞), an analogous argument teaches us that

∥B_{t}^{−1}∥_{L_{β}[0, S]} ≤ \frac{1}{t} (2 − E_{1/2} (−\frac{\sqrt{S}}{t} )) .

Let us illustrate by a simple example that noise amplification by a factor like 1/t is realistic. For 0 ≤ a < b < ∞ the function q̃(σ) = χ_{(a,b]}(σ) is in L_{β} for every β ∈ [1, ∞].

If t > 0 we find from (4.7)

p̃(τ) = \frac{1}{t} \{E_{1/2} (−\frac{\sqrt{b − τ}}{t} ) χ_{(0,b]}(τ) − E_{1/2} (−\frac{\sqrt{a − τ}}{t} ) χ_{(0,a]}(τ)\} ,

in particular

p̃(a−) = \frac{1}{t} \{E_{1/2} (−\frac{\sqrt{b − a}}{t}) − 1\} , E_{1/2}(b−) = 1/t.

In the special case t = 0 by (4.9)

p̃(τ) = \frac{1}{\sqrt{π}} \{(b − τ)^{−1/2}χ_{(0,b]}(τ) − (a − τ)^{−1/2}χ_{(0,a]}(τ)\} ,

and we see that p̃ is unbounded although q̃ is bounded. This illustrates ill-posedness. We also recognize the transition from well-posedness to ill-posedness in the factor 1/t tending to ∞ as t tends to 0.

6 Superposition of δ-functions as density of actual

radius

In the preceding section we have seen the advantage of the positive thickness slice method (t > 0) versus Wicksell’s method of a plane cut (t = 0) with respect to noise amplification (well-posedness versus ill-posedness). In what follows we shall exhibit another advantage (already hinted at in Section 3) of taking t > 0, namely: pointlike particles (spheres with radius 0, in practice spheres with extremely small radius) are undetectable if t = 0. They are hit by a plane cut with probability 0, however can be found with positive probability (if they exist) in a slice of positive thickness. This seems clear to naive intuition, but we shall make it precise.

Consider first the case of finitely many possible values of the actual radius r, equivalently the case of finitely many possible values a_{1}, a_{2}, a_{3}, …, a_{N} of the central cross-section area of a sphere. For definiteness assume N ∈ N,

0 ≤ a_{1} < a_{2} < … < a_{N} < ∞,

and, with all c_{j} > 0,

p̃(τ) = ∑_{j=1}^{N} c_{j}δ(τ − a_{j}).  (6.1)

Then by (4.1)

q̃(σ) = t ∑_{j=1}^{N} c_{j}δ(σ − a_{j}) + \frac{1}{\sqrt{π}} ∑_{j=1}^{N} \frac{c_{j}}{\sqrt{a_{j} − σ}} χ_{(0,a_{j}]}(σ).  (6.2)

Now consider the special case N = 1, a_{1} = 0, c_{1} = c > 0. Then q̃(σ) = ctδ(σ), and for the expected number Q(t) = ∫_{0}^{∞} q̃(σ)dσ of particles producing a shadow in a unit observational area we get

Q(t) = tc  (6.3)

which is positive if t > 0, but 0 if t = 0. Hence pointlike particles can be expected to be seen if and only if t > 0.

If you do not trust manipulations with δ-functions consider the case of extremely small central cross-section areas (equivalent to extremely small actual radii). For h > 0 let p̃_{h}(τ) = \frac{c}{h}χ_{(0,h]}(τ). Then the spatial density of sphere-centers is λ = ∫_{0}^{∞} p̃_{h}(τ)dτ = c > 0, and by (4.1) we obtain

q̃_{h}(σ) = \frac{c}{h} \{t + \frac{2}{\sqrt{π}} (h − σ)^{1/2}\} χ_{(0,h]}(σ).

For the expected number Q_{h}(t) = ∫_{0}^{∞} q̃_{h}(σ)dσ of spheres producing a shadow center in a unit observational area we find

Q_{h}(t) = c (t + \frac{4}{3\sqrt{π}} h^{1/2}) ,

and Q_{h}(t) → tc as h → 0 in agreement with (6.3).

References

[1] G. B a c h, Über die Größenverteilung von Kugelschnitten in durchsichtigen Schnitten endlicher Dicke. Zeitschrift für angewandte Mathematik und Mechanik 38 (1958), 256-258.

[2] G. B a c h, Über die Bestimmung von charakteristischen Größen einer Kugelverteilung aus der Verteilung der Schnittkreise. Zeitschrift für wissenschaftliche Mikroskopie 65 (1963), 285-291.

[3] A. B a d d e l e y, Stochastic geometry, an introduction and reading list. International Statistical Review 50 (1982), 179-193.

[4] A. E r d é l y i e t a l i i, Higher Transcendental Functions, Vol. 3, McGraw Hill Book Company, New York 1955.

[5] P. L. G o l d s m i t h, The calculation of true particle size distributions from the sizes observed in a thin slice. Brit. J. Appl. Phys. 18 (1967), 813-830.

[6] R. G o r e n f l o, Approximation of discrete probability distributions in spherical stereology. International Series of Numerical Mathematics 77 (1986), 103-115.

[7] R. G o r e n f l o, Abel integral equations: applications-motivated solution concepts. Methoden und Verfahren der Mathematischen Physik 34 (1987), 151-174.

[8] R. G o r e n f l o, Computation of rough solutions of Abel integral equations. In: H.W. Engl and C.W. Groetsch (editors): Inverse and Ill-Posed Problems, 195-210, Academic Press, Boston 1987.

[9] R. G o r e n f l o and F. M a i n a r d i, Fractional oscillations and Mittag-Leffler functions. Preprint No. A-14/96, Freie Universität Berlin, Serie A Mathematik [http://www.math.fu-berlin.de/publ/index.html].

[10] R. G o r e n f l o and S. V e s s e l l a, Abel Integral Equations: Analysis and Applications. Lecture Notes in Mathematics 1461, Springer-Verlag, Berlin 1991.

[11] R. G o r e n f l o and M. Y a m a m o t o, On regularized inversion of Abel integral operators. In: N.H. Anh, D.M. Duc, B.D. Khanh, T.D. Van (editors): Analysis and Mechanics of Continuous Media. Proceedings of An International Conference HoChiMinh City 27-29 December 1995 in honor of D.D. Ang. Publications of the HoChiMinh City Mathematical Society Vol. 3, 1995, 162-182.

[12] M. G. K e n d a l l and P. A. M o r a n, Geometric Probability. Charles Griffin & Company, London 1963.

[13] C. K u t s c h e, Produktquadraturverfahren für nicht-glatte Lösungen Abelscher Integralgleichungen erster und zweiter Art. Dissertation Freie Universität Berlin 1994.

[14] K. S. M i l l e r and B. R o s s, An Introduction to the Fractional Calculus and Fractional Differential Equations. John Wiley & Sons, New York 1993.

[15] W. R e i d, Distribution of sizes of spheres in a solid from a study of slices of the solid. Journal of Mathematics and Physics 34 (1955), 95-102.

[16] A. S a n t a l ó, Sobre la distribucion de los tamaños de corpusculos contenidos en un cuerpo a partir de la distribucion en sus secciones o proyecciones. Trabajos de Estadistica 6 (1955), 181-196.

[17] D. S t o y a n and J. M e c k e, Stochastische Geometrie. Akademie-Verlag, Berlin 1983.

[18] S. D. W i c k s e l l, The corpuscle problem. A mathematical study of a biometric problem. Biometrica 17 (1925), 84-99.

Author’s address Prof. Dr. Rudolf Gorenflo First Mathematical Institute Research Group ”Convolutions” Free University of Berlin Arnimallee 2-6 D-14195 Berlin Germany e-mail: gorenflo@math.fu-berlin.de

Supplement (not contained in the paper to be published)

For the interested reader’s convenience we add a few more references in which problems of spherical stereology are treated. This list illustrates the wide range of applicability.

[19] R. S. A n d e r s s o n, The linear functional strategy for improperly posed problems. International Series of Numerical Mathematics 77 (1986), 11-30. (Title of this volume: Inverse Problems, editors J.R.Cannon and U.Hornung).

[20] R. S. A n d e r s s e n and A. J. J a k e m a n, Product integration for functionals of particle size distributions. Utilitas Mathematica 8 (1975), 111-126.

[21] R. S. A n d e r s s e n and A. J. J a k e m a n, On computational Stereology. 6-th Australian Computer Conference Sydney, May 20-24, 1974. Proceedings Vol. II, 353-362.

[22] G. B a c h, Über die Größenverteilung von Kugelschnitten in durchsichtigen Scheiben endlicher Dicke. Zeitschrift für wissenschaftliche Mikroskopie 64 (1959), 265-270.

[23] G. B a c h, Über die Bestimmung von charakteristischen Größen einer Kugelverteilung aus der unvollständigen Verteilung der Schnittkreise. Metrika 9 (1965), 228-233.

[24] R. L. F u l l m a n, Measurement of particle sizes in opaque bodies. Transactions AIME, Journal of Metals March 1953, 447-452.

[25] E. D. H y a m and J. N u t t i n g, The tempering of plain carbon steels. Journal of the Iron and Steel Institute October 1956, 148-165.

[26] A. J a k e m a n and R. S. A n d e r s s o n, Abel type integral equations in stereology. General discussion. Journal of Microscopy 105 (1975), 121-133.

[27] K e n - i c h i K a n a t a n i, Error analysis for the stereological estimation of sphere size distribution: Abel type integral equation. Journal of Computational Physics 57 (1985), 229-250.

[28] W. C. K r u m b e i n, Thin-section mechanical analysis of indurated sediments. Journal of Geology 43 (1935), 482-496.

[29] J. L a s c h and R. K ü h n a u, Evaluation of enzyme distribution in spherical supports by computed fluorescence tomography. Enzyme Microb. Technol. 8 (1986), 115-119.

[30] J. M e c k e and D. S t o y a n, Stereological Problems for Spherical Particles. Mathematische Nachrichten 96 (1980), 311-317.

[31] J. M e i s n e r, Estimation of the distribution of diameters of spherical particles from a given grouped distribution of diameters of observed circles formed by a plane section. Statistica Neerlandica 21 (1967), 11-30.

[32] W. L. N i c h o l s o n, Estimation of linear properties of particle size distributions. Biometrika 57 (1970), 273-297.

[33] D. N y c h k a, G. W a h b a, S. G o l d f a r b, and Th. P u g h, Cross-validated spline methods for the estimation of three-dimensional tumor size distributions from observations on two-dimensional cross-sections. Journal of the American Statistical Society 79 (1984), 832-846.

[34] L. A. S a n t a l ó, Integral Geometry and Geometric Probability. Encyclopedia of Mathematics and its Applications. Addison Wesley Publishing Company, 1976.

[35] E. S c h e i l, Die Berechnung der Anzahl und Größenverteilung kugelförmiger Kristalle in undurchsichtigen Körpern mit Hilfe der durch einen ebenen Schnitt erhaltenen Schnittkreise. Zeitschrift für anorganische und allgemeine Chemie 201 (1931), 259-264.

[36] D. S t o y a n, W. S. K e n d a l l, J. M e c k e, Stochastic Geometry and its Applications. Wiley, 1986.

[37] W. M. V i s s c h e r and A. S. G o l d m a n, Estimating the distribution of spherical particles from plane sections: an optimal algorithm for solution of the Abel integral equation. SIAM J. Sci. Stat. Comput. 4 (1983), 280-290.

[38] S. D. W i c k s e l l, The corpuscle problem. Second memoir. Case of ellipsoidal corpuscles. Biometrica 18 (1926), 151-172.