---
title: Wicksell's Corpuscle Problem Overview
url: https://www.emergentmind.com/topics/wicksell-s-corpuscle-problem
type: topic
---

# Wicksell's Corpuscle Problem Overview

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 [2502.15352], [2310.05463], [1912.01663].

## 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 $f_F$ [2502.15352].

A second formulation, used in recent nonparametric work, is expressed in terms of squared radii. Here $X$ denotes the squared sphere radius, $Y$ the squared section radius, and
$$
Y=(1-U^2)X,\qquad U\sim \mathrm{Uniform}[0,1].
$$
If $g$ denotes the density of $Y$, then
$$
g(y)=\frac{1}{2m_0}\int_y^\infty \frac{dF(x)}{\sqrt{x-y}},\qquad m_0=E[\sqrt X]<\infty.
$$
Defining
$$
V(x)=\int_x^\infty \frac{g(z)}{\sqrt{z-x}}\,dz,
$$
one has
$$
F(x)=1-\frac{V(x)}{V(0)}.
$$
Thus estimation of $F$ reduces to estimation of $V$ at $x$ and at $0$ [2310.05463].

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 [1912.01663].

## 2. Abel-type inversion and ill-posedness

Wicksell’s problem is an inverse problem because the map $F\mapsto G$ 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 [2502.15352].

The classical Abel-type structure appears in several equivalent forms. In the squared-radii setting, the inversion is encoded by
$$
V(x)=\int_x^\infty (z-x)^{-1/2}\,dG(z),
$$
or, with density notation,
$$
V(x)=\int_x^\infty \frac{g(z)}{\sqrt{z-x}}\,dz.
$$
In the Euclidean radii formulation in dimension three, an Abel-type inversion is
$$
f(R)= -\frac{1}{\pi R}\frac{d}{dR}\int_R^\infty \frac{g(r)}{\sqrt{r^2-R^2}}\,dr,
$$
which is recovered as the zero-curvature limit of later curved-space formulas [2508.07762].

The analytic difficulty is that the inverse functional is irregular. Recent frequentist work characterizes the fixed-$x$ fluctuations through the non-standard rate $\sqrt{\log n/n}$ rather than the usual parametric $\sqrt{n}^{-1}$ rate under generic smoothness. This non-standardity is tied to the Abel boundary behavior and persists even after monotonicity regularization [2310.05463]. 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_1,\dots,Y_n$ sampled from $g$, let $G_n$ denote the empirical cdf. The naive plug-in estimator of $V$ is
$$
V_n(x)=\int_x^\infty \frac{dG_n(z)}{\sqrt{z-x}}
      =\frac1n\sum_{i=1}^n (Y_i-x)_+^{-1/2}.
$$
The corresponding plug-in cdf estimator is
$$
F_n^{\mathrm{plug}}(x)=1-\frac{V_n(x)}{V_n(0)}.
$$
As a function of $x$, this estimator is highly irregular and not monotone [2310.05463].

The key regularization is isotonic projection. Define the primitive
$$
U_n(x)=\int_0^x V_n(y)\,dy.
$$
The true $U$ is concave and increasing, but $U_n$ is not. Let $\widetilde U_n$ be the least concave majorant of $U_n$ on $[0,\infty)$, and let
$$
\hat V_n(x)=d^+\widetilde U_n(x)/dx.
$$
Then $\hat V_n$ is a nonincreasing estimator of $V$, and the isotonized plug-in estimator of $F$ is
$$
\hat F_n(x)=1-\frac{\hat V_n(x)}{\hat V_n(0)}.
$$
Equivalently, $\hat F_n$ is the $L_2$-projection of $F_n^{\mathrm{plug}}$ onto the space of nondecreasing right-continuous functions on $\mathbb{R}_+$ evaluated at the observed sample points [2310.05463].

Computation is explicit. One computes $V_n(Y_{(i)})$ at the sorted sample values, sets $W_i=V_n(Y_{(i)})/V_n(0)$, applies the pooled-adjacent-violators algorithm to obtain the nearest nondecreasing sequence, then extends the fit flat between sample points; the complexity is $O(n)$ after sorting [2310.05463]. 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 [2310.05463].

Its asymptotic theory is non-standard but sharp. Under local Hölder-type behavior of order $\gamma_x>1/2$ at $x$ and $\gamma_0$ at $0$,
$$
\sqrt{n/\log n}\,\bigl(\hat F_n(x)-F(x)\bigr)\to \mathcal N(0,\sigma^2),
$$
with
$$
\sigma^2=\frac{4m_0^2}{\pi^2}\left[\frac{g(x)}{2\gamma_x}+(1-F(x))^2\frac{g(0)}{2\gamma_0}\right].
$$
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 [2310.05463].

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

When $F$ or, equivalently, $V$ is exactly constant on an interval containing $x$, the isotonic inverse estimator adapts to the higher rate $\sqrt n$, but its limit distribution is not normal. The limit is the non-Gaussian law of a slope at $x$ of the least concave majorant of a certain Gaussian process [2410.14263]. This establishes an important distinction between adaptive rate improvement and asymptotic efficiency.

If the interval of constancy $[\underline x,\bar x]$ is known a priori, three informed projection-type estimators can be constructed on
$$
\mathcal V_{\underline x,\bar x}
=\{\,V\in\mathcal V:V\text{ is constant on }[\underline x,\bar x]\}.
$$
They are: a direct naive-projection estimator, a two-stage projection that flattens the isotonic inverse estimator on $[\underline x,\bar x]$, and a direct projection of $V_n$ onto piecewise-constant functions on a fine grid. All three agree to first order on the constant interval [2410.14263].

Under the moment condition $\int s^{3/2}dF(s)<\infty$ and exact constancy of $V$ on $[\underline x,\bar x]$,
$$
\sqrt n\bigl(V_n^{(\underline x,\bar x)}(x)-V(x)\bigr)
\rightsquigarrow N(0,\sigma^2_{\underline x,\bar x}),
$$
with
$$
\sigma^2_{\underline x,\bar x}
=\Var\!\Bigl[2(\bar x-\underline x)^{-1}
\bigl\{\sqrt{(Z-\underline x)_+}-\sqrt{(Z-\bar x)_+}\bigr\}\Bigr].
$$
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 [2410.14263].

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,
$$
L_x=W+N(0,\sigma^2_{\underline x,\bar x}),
$$
with $W$ independent of the normal part. A simulation with a model having a flat plateau on $[2,3]$ found
$$
\SD[\sqrt n\,(V_n^{(2,3)}-V)(2.5)]\approx 0.5464,\qquad
\SD[\sqrt n\,(\hat V_n-V)(2.5)]\approx 0.5697,
$$
close to the theoretical $\sigma_{2,3}=0.5466$; the finite-sample differences were therefore small even though the asymptotic limit laws differ [2410.14263].

## 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$. Instead, it places a Dirichlet-process prior directly on the observable distribution:
$$
G\sim \mathrm{DP}(\alpha),
$$
where $\alpha$ is a finite base measure on $[0,\infty)$. Given an i.i.d. sample $Z_1,\dots,Z_n$ from the true observable law $G_0$, conjugacy gives
$$
G\mid Z_1,\dots,Z_n \sim \mathrm{DP}(\alpha+nG_n),
$$
with $G_n=(1/n)\sum_{i=1}^n\delta_{Z_i}$ [2502.15352].

The naive Bayes posterior for the inverse functional is obtained by plugging posterior draws of $G$ into
$$
V_G(x)=\int_x^\infty (z-x)^{-1/2}\,dG(z).
$$
However, the true inverse functional $V_0(x)$ 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
$$
U_G(y)=\int_0^y V_G(t)\,dt,
$$
one computes the least concave majorant of $U_G$ and takes its right-derivative:
$$
\hat V_G(x)=d^+\,\mathrm{LCM}\{U_G(t)\}/dt.
$$
Equivalently, $\hat V_G$ is the $L^2$-projection of $V_G$ onto the relevant convex cone, and the projection can be implemented by PAVA [2502.15352].

Under a local Hölder condition at $x$ with exponent $\gamma>1/2$,
$$
\int_0^1 \bigl(V_0(x)-V_0(x+u\delta)\bigr)\,du \sim K\delta^\gamma
\qquad (\delta\to 0),
$$
the posterior satisfies a semiparametric Bernstein–von Mises theorem. If $\hat V_n(x)$ is the isotonized empirical plug-in estimator and
$$
\Delta_G(x)=\sqrt{n/\log n}\,\bigl(\hat V_G(x)-\hat V_n(x)\bigr),
$$
then, in $G_0$-probability,
$$
\Delta_G(x)\mid Z_1,\dots,Z_n \Rightarrow N\bigl(0,g_0(x)/(2\gamma)\bigr),
$$
where $g_0(x)=f_{G_0}(x)$ is the true density of the visible cross-section radii at $x$ [2502.15352].

This yields automatic uncertainty quantification. Pointwise credible intervals of the form
$$
\hat V_n(x)\pm z_{1-\alpha/2}\sqrt{g_0(x)/(2\gamma\,n/\log n)}
$$
have asymptotically correct frequentist coverage, and no separate estimation of $\gamma$ 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 [2502.15352].

## 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 $10$–$200$ 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 [2008.09091].

Poliakova introduced a density approximation based on a “polygonal revolution” model. A regular $4m$-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
$$
a=\frac{\pi}{2m\sin(\pi/(2m))},
$$
the approximating profile-diameter density becomes
$$
g(y)=\frac{1}{aE[D]}
\sum_{i=1}^m
p_i\,
\frac{F(y/(ax_i))-F(y/(ax_{i-1}))}{x_{i-1}-x_i},
$$
with $x_i=\cos(i\pi/(2m))$ and
$$
p_i=\sin(i\pi/(2m))-\sin((i-1)\pi/(2m)).
$$
In practice $m=10$–$15$ already gives density values within $0.1\%$ of the exact Wicksell integral, and the approximation is substantially cheaper numerically [2008.09091].

For a parametric family $f(R;\theta)$, one observes independent profile diameters $y_1,\dots,y_n$ and forms
$$
L(\theta)=\prod_{j=1}^n g(y_j;\theta),\qquad
\ell(\theta)=\sum_{j=1}^n \log g(y_j;\theta).
$$
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 $3$–$5$, and the overall MLE search by a factor of $7$–$20$ [2008.09091].

The simulation findings are specific. For $n\ge 50$ profiles, maximum likelihood has bias and standard deviation typically $10$–$30\%$ 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 $n<200$; and $m=8$ already gives essentially the same standard deviation as $m=100$, while $m=15$ was used throughout [2008.09091]. 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 $-2\ell(\hat\theta)$ [2008.09091].

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 $E[D]$ were $3.48\,(2.82,4.16)$ for “3437-6 layer,” $8.82\,(5.96,11.51)$ for “3437-6 matrix,” $6.71\,(4.85,8.56)$ for “3434-7,” $4.09\,(3.74,4.44)$ for “3438-3,” and $3.71\,(2.84,4.56)$ for “3438-3 sub” [2008.09091].

## 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 $\mathbb{R}^3$, Santaló’s integral equations replace the sphere-specific Abel kernel. If $H(\lambda)d\lambda$ is the number-density of particles of scale $\lambda$, and $\phi$ 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 $(\mathrm{RP})$ and $(\mathrm{RL})$. 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
$$
\mathbb H_{pl}(\lambda)=\frac{1}{\alpha}\,\mathfrak M^{-1}\!\Bigl(\frac{h^*(s)}{\phi^*(s)}\Bigr)(\lambda),
$$
and the line-section analogue is
$$
\mathbb H_{li}(\lambda)=\frac{1}{\beta\,\lambda^2}\,
\mathfrak M^{-1}\!\Bigl(\frac{h^*(s)}{\phi^*(s)}\Bigr)(\lambda).
$$
Under the stated Mellin-integrability assumptions, these formulas yield partial existence and uniqueness results [1912.01663].

When the template body $K$ 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 $\phi^*(s)$ [1912.01663].

A second extension replaces Euclidean $\mathbb{R}^d$ by a constant-curvature space $M_k^d$. For a stationary process of random balls invariant under the full isometry group, and a fixed totally geodesic hypersurface $L$, the induced section-radius law $\nu_{d-1}$ is related to the original radius law $\nu_d$ by
$$
g_k(r)=\frac{N_d}{N_{d-1}}\int_{R=r}^{l_k} K_k(R,r)f_k(R)\,dR,
$$
where $K_k(R,r)$ is obtained by differentiating an equidistant-decomposition integral involving $\cos^{d-1}(\sqrt{k}\,h)$ for $k>0$, $\cosh^{d-1}(\sqrt{|k|}\,h)$ for $k<0$, and the Euclidean limit when $k=0$ [2508.07762].

The curved-space theory includes explicit inversion formulas for both negative and positive curvature and proves that, as $k\to 0$, the section and inversion formulas converge to the classical Euclidean results. The stated applications include biology, materials science, planetary geology, network science, and cosmology [2508.07762]. 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.

Source: https://www.emergentmind.com/topics/wicksell-s-corpuscle-problem