Wicksell's Corpuscle Problem Overview

• Wicksell’s corpuscle problem is the inverse stereological challenge of reconstructing 3D size distributions of spherical particles from 2D cross-sectional data, pivotal in fields like materials science and astronomy.
• It involves an ill-posed Abel-type integral inversion that requires regularization methods such as isotonic projection to enforce monotonicity and stability in estimation.
• Recent advancements integrate nonparametric and Bayesian approaches to handle the unique rate irregularities and shape constraints inherent in the inversion process.

Searching arXiv for recent and relevant papers on Wicksell's corpuscle problem. Wicksell’s corpuscle problem is the inverse stereological problem of recovering the size distribution of three-dimensional particles from the distribution of two-dimensional planar sections. In the classical setting, the unobserved corpuscles are spheres in $\mathbb{R}^3$, and the observations are circular cross-sections produced by a random plane. The problem arises in astronomy, materials science, microscopy, and petrography, because only projected or sectioned data are available while the scientific target is a latent three-dimensional distribution. The forward map from the latent distribution to the observed one is an Abel-type first-kind integral transform, so direct inversion is ill-posed without additional structure or regularization (Gili et al., 21 Feb 2025, Gili et al., 2023, Kiseľák et al., 2019).

1. Geometric and probabilistic formulation

In one formulation, unseen spheres have radii $X$ with distribution function $F(x)$, and slicing by a plane yields observed circular radii $Y$ with distribution function $G(y)$. The relation is

$G(y)=\frac{1}{E[X]}\int_y^\infty (x-y)\,dF(x),$

and, in density form,

$f_G(y)=\frac{1}{E[X]}\int_y^\infty f_F(x)\,dx.$

These formulas define the inverse problem: given $G$ or $f_G$, recover $F$ or $X$0 (Gili et al., 21 Feb 2025).

A second formulation, used in recent nonparametric work, is expressed in terms of squared radii. Here $X$1 denotes the squared sphere radius, $X$2 the squared section radius, and

$X$3

If $X$4 denotes the density of $X$5, then

$X$6

Defining

$X$7

one has

$X$8

Thus estimation of $X$9 reduces to estimation of $F(x)$0 at $F(x)$1 and at $F(x)$2 (Gili et al., 2023).

The original stereological content is already visible in Wicksell’s classical radii formulation: a random plane intersects only particles large enough to be hit, and larger particles are sampled with higher probability. This size-biasing is intrinsic to the forward operator and is one source of the nontrivial inversion geometry emphasized in later work (Kiseľák et al., 2019).

2. Abel-type inversion and ill-posedness

Wicksell’s problem is an inverse problem because the map $F(x)$3 is a compact integral operator of the first kind. In the Bayesian treatment based on a prior for the observable law, this compactness is stated explicitly: direct inversion is ill-posed without further constraints (Gili et al., 21 Feb 2025).

The classical Abel-type structure appears in several equivalent forms. In the squared-radii setting, the inversion is encoded by

$F(x)$4

or, with density notation,

$F(x)$5

In the Euclidean radii formulation in dimension three, an Abel-type inversion is

$F(x)$6

which is recovered as the zero-curvature limit of later curved-space formulas (Spanos et al., 11 Aug 2025).

The analytic difficulty is that the inverse functional is irregular. Recent frequentist work characterizes the fixed-$F(x)$7 fluctuations through the non-standard rate $F(x)$8 rather than the usual parametric $F(x)$9 rate under generic smoothness. This non-standardity is tied to the Abel boundary behavior and persists even after monotonicity regularization (Gili et al., 2023). A plausible implication is that successful procedures must exploit structural information—most prominently monotonicity, concavity of primitives, or explicit parametric modeling—rather than rely on naive plug-in inversion.

3. Isotonic inverse estimation in the classical nonparametric model

Given $Y$0 sampled from $Y$1, let $Y$2 denote the empirical cdf. The naive plug-in estimator of $Y$3 is

$Y$4

The corresponding plug-in cdf estimator is

$Y$5

As a function of $Y$6, this estimator is highly irregular and not monotone (Gili et al., 2023).

The key regularization is isotonic projection. Define the primitive

$Y$7

The true $Y$8 is concave and increasing, but $Y$9 is not. Let $G(y)$0 be the least concave majorant of $G(y)$1 on $G(y)$2, and let

$G(y)$3

Then $G(y)$4 is a nonincreasing estimator of $G(y)$5, and the isotonized plug-in estimator of $G(y)$6 is

$G(y)$7

Equivalently, $G(y)$8 is the $G(y)$9-projection of $G(y)=\frac{1}{E[X]}\int_y^\infty (x-y)\,dF(x),$0 onto the space of nondecreasing right-continuous functions on $G(y)=\frac{1}{E[X]}\int_y^\infty (x-y)\,dF(x),$1 evaluated at the observed sample points (Gili et al., 2023).

Computation is explicit. One computes $G(y)=\frac{1}{E[X]}\int_y^\infty (x-y)\,dF(x),$2 at the sorted sample values, sets $G(y)=\frac{1}{E[X]}\int_y^\infty (x-y)\,dF(x),$3, applies the pooled-adjacent-violators algorithm to obtain the nearest nondecreasing sequence, then extends the fit flat between sample points; the complexity is $G(y)=\frac{1}{E[X]}\int_y^\infty (x-y)\,dF(x),$4 after sorting (Gili et al., 2023). This yields a tuning-free procedure: unlike kernel or spline estimators, the method uses no smoothing parameter, and the curvature of the least concave majorant automatically balances bias and variance locally (Gili et al., 2023).

Its asymptotic theory is non-standard but sharp. Under local Hölder-type behavior of order $G(y)=\frac{1}{E[X]}\int_y^\infty (x-y)\,dF(x),$5 at $G(y)=\frac{1}{E[X]}\int_y^\infty (x-y)\,dF(x),$6 and $G(y)=\frac{1}{E[X]}\int_y^\infty (x-y)\,dF(x),$7 at $G(y)=\frac{1}{E[X]}\int_y^\infty (x-y)\,dF(x),$8,

$G(y)=\frac{1}{E[X]}\int_y^\infty (x-y)\,dF(x),$9

with

$f_G(y)=\frac{1}{E[X]}\int_y^\infty f_F(x)\,dx.$0

A local asymptotic minimax lower bound shows that the isotonic inverse estimator attains the local minimax bound, and the efficiency proof is organized through a two-dimensional LAN perturbation path and the Hájek–Le Cam convolution theorem (Gili et al., 2023).

4. Flat regions, informed projections, and the limits of adaptivity

When $f_G(y)=\frac{1}{E[X]}\int_y^\infty f_F(x)\,dx.$1 or, equivalently, $f_G(y)=\frac{1}{E[X]}\int_y^\infty f_F(x)\,dx.$2 is exactly constant on an interval containing $f_G(y)=\frac{1}{E[X]}\int_y^\infty f_F(x)\,dx.$3, the isotonic inverse estimator adapts to the higher rate $f_G(y)=\frac{1}{E[X]}\int_y^\infty f_F(x)\,dx.$4, but its limit distribution is not normal. The limit is the non-Gaussian law of a slope at $f_G(y)=\frac{1}{E[X]}\int_y^\infty f_F(x)\,dx.$5 of the least concave majorant of a certain Gaussian process (Gili et al., 2024). This establishes an important distinction between adaptive rate improvement and asymptotic efficiency.

If the interval of constancy $f_G(y)=\frac{1}{E[X]}\int_y^\infty f_F(x)\,dx.$6 is known a priori, three informed projection-type estimators can be constructed on

$f_G(y)=\frac{1}{E[X]}\int_y^\infty f_F(x)\,dx.$7

They are: a direct naive-projection estimator, a two-stage projection that flattens the isotonic inverse estimator on $f_G(y)=\frac{1}{E[X]}\int_y^\infty f_F(x)\,dx.$8, and a direct projection of $f_G(y)=\frac{1}{E[X]}\int_y^\infty f_F(x)\,dx.$9 onto piecewise-constant functions on a fine grid. All three agree to first order on the constant interval (Gili et al., 2024).

Under the moment condition $G$0 and exact constancy of $G$1 on $G$2,

$G$3

with

$G$4

A local asymptotic minimax lower bound shows that no estimator can do better in asymptotic risk for smooth, symmetric losses, so these informed estimators are asymptotically efficient (Gili et al., 2024).

The same paper proves that the isotonic inverse estimator is not efficient on constant zones when such local information is available. Its limit law satisfies a convolution decomposition,

$G$5

with $G$6 independent of the normal part. A simulation with a model having a flat plateau on $G$7 found

$G$8

close to the theoretical $G$9; the finite-sample differences were therefore small even though the asymptotic limit laws differ (Gili et al., 2024).

5. Bayesian inference via the isotonized inverse posterior

A recent Bayesian approach departs from the classical strategy of placing a nonparametric prior on the latent distribution $f_G$0. Instead, it places a Dirichlet-process prior directly on the observable distribution:

$f_G$1

where $f_G$2 is a finite base measure on $f_G$3. Given an i.i.d. sample $f_G$4 from the true observable law $f_G$5, conjugacy gives

$f_G$6

with $f_G$7 (Gili et al., 21 Feb 2025).

The naive Bayes posterior for the inverse functional is obtained by plugging posterior draws of $f_G$8 into

$f_G$9

However, the true inverse functional $F$0 is nonincreasing and right-continuous, whereas naive plug-in draws need not satisfy this shape constraint. The isotonized inverse posterior is therefore defined by projecting each posterior draw onto the cone of nonincreasing, right-continuous functions. Writing

$F$1

one computes the least concave majorant of $F$2 and takes its right-derivative:

$F$3

Equivalently, $F$4 is the $F$5-projection of $F$6 onto the relevant convex cone, and the projection can be implemented by PAVA (Gili et al., 21 Feb 2025).

Under a local Hölder condition at $F$7 with exponent $F$8,

$F$9

the posterior satisfies a semiparametric Bernstein–von Mises theorem. If $X$00 is the isotonized empirical plug-in estimator and

$X$01

then, in $X$02-probability,

$X$03

where $X$04 is the true density of the visible cross-section radii at $X$05 (Gili et al., 21 Feb 2025).

This yields automatic uncertainty quantification. Pointwise credible intervals of the form

$X$06

have asymptotically correct frequentist coverage, and no separate estimation of $X$07 is required. The result is stated as the first semiparametric Bernstein–von Mises theorem for projection-based posteriors with a Dirichlet-process prior in inverse problems (Gili et al., 21 Feb 2025).

6. Parametric likelihood methods and small-sample practice

A different strand of work studies Wicksell’s problem under parametric assumptions, motivated by microscopy applications with only a few $X$08–$X$09 profile measurements. In this setting, the particles are approximately spherical, parametric assumptions are regarded as reasonable, and the inferential target is the parameter vector of the three-dimensional size distribution rather than a fully nonparametric cdf (Poliakova, 2020).

Poliakova introduced a density approximation based on a “polygonal revolution” model. A regular $X$10-sided polygon is inscribed in the projected disk and revolved about one of its diagonals, producing a piecewise-conical solid whose profile density is easier to evaluate. After rescaling by

$X$11

the approximating profile-diameter density becomes

$X$12

with $X$13 and

$X$14

In practice $X$15–$X$16 already gives density values within $X$17 of the exact Wicksell integral, and the approximation is substantially cheaper numerically (Poliakova, 2020).

For a parametric family $X$18, one observes independent profile diameters $X$19 and forms

$X$20

All computations in the paper were done in R using optim() with method="BFGS" or "L-BFGS-B". Relative to brute-force trapezoidal integration of the Wicksell integral, the polygonal approximation speeds up each likelihood evaluation by a factor of $X$21–$X$22, and the overall MLE search by a factor of $X$23–$X$24 (Poliakova, 2020).

The simulation findings are specific. For $X$25 profiles, maximum likelihood has bias and standard deviation typically $X$26–$X$27 smaller than method of moments in log-normal models, with comparable performance in Weibull models; minimum-distance estimation is much slower and exhibits both larger bias and variance, especially for $X$28; and $X$29 already gives essentially the same standard deviation as $X$30, while $X$31 was used throughout (Poliakova, 2020). Confidence regions are based on Wilks’ theorem, and model selection is carried out by comparing maximized log-likelihoods or, since the candidate families have the same number of parameters, equivalently comparing $X$32 (Poliakova, 2020).

The applied scope includes glacier ice petrography. In the Vostok ice-core example, model choice by AIC favored Weibull in all but one smallest sample, and the reported ML-1 estimates of the mean diameter $X$33 were $X$34 for “3437-6 layer,” $X$35 for “3437-6 matrix,” $X$36 for “3434-7,” $X$37 for “3438-3,” and $X$38 for “3438-3 sub” (Poliakova, 2020).

7. Extensions beyond the classical Euclidean spherical model

The classical corpuscle problem has been generalized in two directions: from spheres to arbitrary convex similar bodies, and from Euclidean space to spaces of constant curvature.

For convex similar bodies in $X$39, Santaló’s integral equations replace the sphere-specific Abel kernel. If $X$40 is the number-density of particles of scale $X$41, and $X$42 describes the section statistics of a fixed template body, then the random-plane and random-line observation models become integral equations of the forms labeled $X$43 and $X$44. Kiseľák and Balúchová solve these equations by the Method of Model Solutions, implemented through the Mellin transform. The principal plane-section solution is

$X$45

and the line-section analogue is

$X$46

Under the stated Mellin-integrability assumptions, these formulas yield partial existence and uniqueness results (Kiseľák et al., 2019).

When the template body $X$47 is a sphere, the general Mellin-ratio inversion reproduces Wicksell’s Abel kernel and its classical inversion. In this sense, the spherical corpuscle problem is a special case of a broader stereological inversion theory in which the geometry of the particle class is encoded through the Mellin moments $X$48 (Kiseľák et al., 2019).

A second extension replaces Euclidean $X$49 by a constant-curvature space $X$50. For a stationary process of random balls invariant under the full isometry group, and a fixed totally geodesic hypersurface $X$51, the induced section-radius law $X$52 is related to the original radius law $X$53 by

$X$54

where $X$55 is obtained by differentiating an equidistant-decomposition integral involving $X$56 for $X$57, $X$58 for $X$59, and the Euclidean limit when $X$60 (Spanos et al., 11 Aug 2025).

The curved-space theory includes explicit inversion formulas for both negative and positive curvature and proves that, as $X$61, the section and inversion formulas converge to the classical Euclidean results. The stated applications include biology, materials science, planetary geology, network science, and cosmology (Spanos et al., 11 Aug 2025). This suggests that Wicksell’s problem is best understood not as a single Abel inversion formula, but as a family of geometrically structured inverse problems whose kernels depend on the ambient space and on the underlying particle class.