Radon transform

Definition

Let ${\displaystyle f(\mathbf{\text{x}})=f(x,y)}$ be a function that satisfies the three regularity conditions:

1. ${\displaystyle f(\mathbf{\text{x}})}$ is continuous;
2. the double integral ${\displaystyle \iint {\displaystyle \frac{|f(\mathbf{\text{x}})|}{\sqrt{x^2+y^2}}}dxdy}$, extending over the whole plane, converges;
3. for any arbitrary point ${\displaystyle (x,y)}$ on the plane it holds that ${\displaystyle \underset{r\to \mathrm{\infty }}{lim}{\int }_0^{2\pi }f(x+r\cos \phi ,y+r\sin \phi )d\phi =0.}$

The Radon transform, ${\displaystyle Rf}$, is a function defined on the space of straight lines ${\displaystyle L\subset {\mathbb{R}}^2}$ by the line integral along each such line as: $${\displaystyle Rf(L)={\int }_{L}f(\mathbf{x})|d\mathbf{x}|.}$$Concretely, the parametrization of any straight line ${\displaystyle L}$ with respect to arc length ${\displaystyle z}$ can always be written:$${\displaystyle (x(z),y(z))=((z\sin \alpha +s\cos \alpha ),(-z\cos \alpha +s\sin \alpha ))}$$where ${\displaystyle s}$ is the distance of ${\displaystyle L}$ from the origin and ${\displaystyle \alpha }$ is the angle the normal vector to ${\displaystyle L}$ makes with the ${\displaystyle X}$-axis. It follows that the quantities ${\displaystyle (\alpha ,s)}$ can be considered as coordinates on the space of all lines in ${\displaystyle {\mathbb{R}}^2}$, and the Radon transform can be expressed in these coordinates by: $${\displaystyle \begin{array}{cc}Rf(\alpha ,s)& ={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}f(x(z),y(z))dz\\ & ={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}f((z\sin \alpha +s\cos \alpha ),(-z\cos \alpha +s\sin \alpha ))dz.\end{array}}$$More generally, in the ${\displaystyle n}$-dimensional Euclidean space ${\displaystyle {\mathbb{R}}^n}$, the Radon transform of a function ${\displaystyle f}$ satisfying the regularity conditions is a function ${\displaystyle Rf}$ on the space ${\displaystyle {\mathrm{\Sigma }}_n}$ of all hyperplanes in ${\displaystyle {\mathbb{R}}^n}$. It is defined by:

$${\displaystyle Rf(\xi )={\int }_{\xi }f(\mathbf{x})d\sigma (\mathbf{x}),\mathrm{\forall }\xi \in {\mathrm{\Sigma }}_n}$$where the integral is taken with respect to the natural hypersurface measure, ${\displaystyle d\sigma }$ (generalizing the ${\displaystyle |d\mathbf{x}|}$ term from the ${\displaystyle 2}$-dimensional case). Observe that any element of ${\displaystyle {\mathrm{\Sigma }}_n}$ is characterized as the solution locus of an equation ${\displaystyle \mathbf{x}\cdot \alpha =s}$, where ${\displaystyle \alpha \in S^{n-1}}$ is a unit vector and ${\displaystyle s\in \mathbb{R}}$. Thus the ${\displaystyle n}$-dimensional Radon transform may be rewritten as a function on ${\displaystyle S^{n-1}\times \mathbb{R}}$ via: $${\displaystyle Rf(\alpha ,s)={\int }_{\mathbf{x}\cdot \alpha =s}f(\mathbf{x})d\sigma (\mathbf{x}).}$$It is also possible to generalize the Radon transform still further by integrating instead over ${\displaystyle k}$-dimensional affine subspaces of ${\displaystyle {\mathbb{R}}^n}$. The X-ray transform is the most widely used special case of this construction, and is obtained by integrating over straight lines.

Relationship with the Fourier transform

The Radon transform is closely related to the Fourier transform. We define the univariate Fourier transform here as: $${\displaystyle \hat{f}(\omega )={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}f(x)e^{-2\pi ix\omega }dx.}$$ For a function of a ${\displaystyle 2}$-vector ${\displaystyle \mathbf{x}=(x,y)}$, the univariate Fourier transform is: $${\displaystyle \hat{f}(\mathbf{w})={\iint }_{{\mathbb{R}}^2}f(\mathbf{x})e^{-2\pi i\mathbf{x}\cdot \mathbf{w}}dxdy.}$$ For convenience, denote ${\displaystyle {\mathcal{R}}_{\alpha }[f](s)=\mathcal{R}[f](\alpha ,s)}$. The Fourier slice theorem then states: $${\displaystyle \widehat{{\mathcal{R}}_{\alpha }[f]}(\sigma )=\hat{f}(\sigma \mathbf{n}(\alpha ))}$$ where ${\displaystyle \mathbf{n}(\alpha )=(\cos \alpha ,\sin \alpha ).}$

Thus the two-dimensional Fourier transform of the initial function along a line at the inclination angle ${\displaystyle \alpha }$ is the one variable Fourier transform of the Radon transform (acquired at angle ${\displaystyle \alpha }$) of that function. This fact can be used to compute both the Radon transform and its inverse. The result can be generalized into n dimensions: $${\displaystyle \hat{f}(r\alpha )={\int }_{\mathbb{R}}\mathcal{R}f(\alpha ,s)e^{-2\pi isr}ds.}$$

Dual transform

The dual Radon transform is a kind of adjoint to the Radon transform. Beginning with a function g on the space ${\displaystyle {\mathrm{\Sigma }}_n}$, the dual Radon transform is the function ${\displaystyle {\mathcal{R}}^{\ast }g}$ on Rⁿ defined by: $${\displaystyle {\mathcal{R}}^{\ast }g(\mathbf{x})={\int }_{\mathbf{x}\in \xi }g(\xi )d\mu (\xi ).}$$The integral here is taken over the set of all hyperplanes incident with the point ${\displaystyle \mathbf{\text{x}}\in {\mathbb{R}}^n}$, and the measure ${\displaystyle d\mu }$ is the unique probability measure on the set ${\displaystyle \{\xi |\mathbf{x}\in \xi \}}$ invariant under rotations about the point ${\displaystyle \mathbf{x}}$.

Concretely, for the two-dimensional Radon transform, the dual transform is given by: $${\displaystyle {\mathcal{R}}^{\ast }g(\mathbf{x})=\frac{1}{2\pi }{\int }_{\alpha =0}^{2\pi }g(\alpha ,\mathbf{n}(\alpha )\cdot \mathbf{x})d\alpha .}$$ In the context of image processing, the dual transform is commonly called back-projection as it takes a function defined on each line in the plane and 'smears' or projects it back over the line to produce an image.

Intertwining property

Let ${\displaystyle \mathrm{\Delta }}$ denote the Laplacian on ${\displaystyle {\mathbb{R}}^n}$ defined by:$${\displaystyle \mathrm{\Delta }=\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}_1^2}+\cdots +\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}_n^2}}$$This is a natural rotationally invariant second-order differential operator. On ${\displaystyle {\mathrm{\Sigma }}_n}$, the "radial" second derivative ${\displaystyle Lf(\alpha ,s)\equiv \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{s}^2}f(\alpha ,s)}$ is also rotationally invariant. The Radon transform and its dual are intertwining operators for these two differential operators in the sense that: $${\displaystyle \mathcal{R}(\mathrm{\Delta }f)=L(\mathcal{R}f),{\mathcal{R}}^{\ast }(Lg)=\mathrm{\Delta }({\mathcal{R}}^{\ast }g).}$$In analysing the solutions to the wave equation in multiple spatial dimensions, the intertwining property leads to the translational representation of Lax and Philips. In imaging and numerical analysis this is exploited to reduce multi-dimensional problems into single-dimensional ones, as a dimensional splitting method.

Reconstruction approaches

The process of reconstruction produces the image (or function ${\displaystyle f}$ in the previous section) from its projection data. Reconstruction is an inverse problem.

Radon inversion formula

In the two-dimensional case, the most commonly used analytical formula to recover ${\displaystyle f}$ from its Radon transform is the Filtered Back-projection Formula or Radon Inversion Formula: $${\displaystyle f(\mathbf{x})={\int }_0^{\pi }(\mathcal{R}f(\cdot ,\theta )\ast h)(⟨\mathbf{x},{\mathbf{n}}_{\theta }⟩)d\theta }$$where ${\displaystyle h}$ is such that ${\displaystyle \hat{h}(k)=|k|}$. The convolution kernel ${\displaystyle h}$ is referred to as Ramp filter in some literature.

Ill-posedness

Intuitively, in the filtered back-projection formula, by analogy with differentiation, for which $\left(\widehat{\frac{d}{dx}f}\right)(k)=ik\hat{f}(k)$, we see that the filter performs an operation similar to a derivative. Roughly speaking, then, the filter makes objects more singular. A quantitive statement of the ill-posedness of Radon inversion goes as follows:$${\displaystyle \widehat{{\mathcal{R}}^{\ast }\mathcal{R}[g]}(\mathbf{k})=\frac{1}{\Vert \mathbf{k}\Vert }\hat{g}(\mathbf{k})}$$ where ${\displaystyle {\mathcal{R}}^{\ast }}$ is the previously defined adjoint to the Radon transform. Thus for ${\displaystyle g(\mathbf{x})=e^{i⟨{\mathbf{k}}_0,\mathbf{x}⟩}}$, we have: $${\displaystyle {\mathcal{R}}^{\ast }\mathcal{R}[g](\mathbf{x})=\frac{1}{\Vert {\mathbf{k}}_{\mathbf{0}}\Vert }{e}^{i⟨{\mathbf{k}}_0,\mathbf{x}⟩}}$$ The complex exponential ${\displaystyle e^{i⟨{\mathbf{k}}_0,\mathbf{x}⟩}}$ is thus an eigenfunction of ${\displaystyle {\mathcal{R}}^{\ast }\mathcal{R}}$ with eigenvalue $\frac{1}{\Vert {\mathbf{k}}_0\Vert }$. Thus the singular values of ${\displaystyle \mathcal{R}}$ are $\frac{1}{\sqrt{\Vert \mathbf{k}\Vert }}$. Since these singular values tend to ${\displaystyle 0}$, ${\displaystyle {\mathcal{R}}^{-1}}$ is unbounded.

Iterative reconstruction methods

Compared with the Filtered Back-projection method, iterative reconstruction costs large computation time, limiting its practical use. However, due to the ill-posedness of Radon Inversion, the Filtered Back-projection method may be infeasible in the presence of discontinuity or noise. Iterative reconstruction methods (e.g. iterative Sparse Asymptotic Minimum Variance) could provide metal artefact reduction, noise and dose reduction for the reconstructed result that attract much research interest around the world.

Inversion formulas

Explicit and computationally efficient inversion formulas for the Radon transform and its dual are available. The Radon transform in ${\displaystyle n}$ dimensions can be inverted by the formula: $${\displaystyle c_{n}f=(-\mathrm{\Delta }{)}^{(n-1)/2}{R}^{\ast }Rf}$$where ${\displaystyle c_n=(4\pi {)}^{(n-1)/2}\frac{\mathrm{\Gamma }(n/2)}{\mathrm{\Gamma }(1/2)}}$, and the power of the Laplacian ${\displaystyle (-\mathrm{\Delta }{)}^{(n-1)/2}}$ is defined as a pseudo-differential operator if necessary by the Fourier transform: $${\displaystyle \left[\mathcal{F}(-\mathrm{\Delta }{)}^{(n-1)/2}\phi \right](\xi )=|2\pi \xi {|}^{n-1}(\mathcal{F}\phi )(\xi ).}$$For computational purposes, the power of the Laplacian is commuted with the dual transform ${\displaystyle R^{\ast }}$ to give: $${\displaystyle c_{n}f=\{\begin{array}{ll}{R}^{\ast }\frac{d^{n-1}}{d{s}^{n-1}}Rf& n\text{ odd}\\ R^{\ast }{\mathcal{H}}_s\frac{d^{n-1}}{d{s}^{n-1}}Rf& n\text{ even}\end{array}}$$where ${\displaystyle {\mathcal{H}}_s}$ is the Hilbert transform with respect to the s variable. In two dimensions, the operator ${\displaystyle {\mathcal{H}}_s\frac{d}{ds}}$ appears in image processing as a ramp filter. One can prove directly from the Fourier slice theorem and change of variables for integration that for a compactly supported continuous function ${\displaystyle f}$ of two variables: $${\displaystyle f=\frac{1}{2}{R}^{\ast }{\mathcal{H}}_s\frac{d}{ds}Rf.}$$Thus in an image processing context the original image ${\displaystyle f}$ can be recovered from the 'sinogram' data ${\displaystyle Rf}$ by applying a ramp filter (in the ${\displaystyle s}$ variable) and then back-projecting. As the filtering step can be performed efficiently (for example using digital signal processing techniques) and the back projection step is simply an accumulation of values in the pixels of the image, this results in a highly efficient, and hence widely used, algorithm.

Explicitly, the inversion formula obtained by the latter method is: $${\displaystyle f(x)=\{\begin{array}{ll}{\displaystyle -ı2\pi (2\pi {)}^{-n}(-1{)}^{n/2}{\int }_{S^{n-1}}\frac{{\mathrm{\partial }}^{n-1}}{2\mathrm{\partial }{s}^{n-1}}Rf(\alpha ,\alpha \cdot x)d\alpha }& n\text{ odd}\\ {\displaystyle (2\pi {)}^{-n}(-1{)}^{n/2}{\iint }_{\mathbb{R}\times S^{n-1}}\frac{{\mathrm{\partial }}^{n-1}}{q\mathrm{\partial }{s}^{n-1}}Rf(\alpha ,\alpha \cdot x+q)d\alpha dq}& n\text{ even}\end{array}}$$The dual transform can also be inverted by an analogous formula: $${\displaystyle c_{n}g=(-L{)}^{(n-1)/2}R(R^{\ast }g).}$$

Radon transform in algebraic geometry

In algebraic geometry, a Radon transform (also known as the Brylinski–Radon transform) is constructed as follows.

Write

${\displaystyle {\mathbf{P}}^d\stackrel{p_1}{\leftarrow }H\stackrel{p_2}{\to }{\mathbf{P}}^{\vee ,d}}$

for the universal hyperplane, i.e., H consists of pairs (x, h) where x is a point in d-dimensional projective space ${\displaystyle {\mathbf{P}}^d}$ and h is a point in the dual projective space (in other words, x is a line through the origin in (d+1)-dimensional affine space, and h is a hyperplane in that space) such that x is contained in h.

Then the Brylinksi–Radon transform is the functor between appropriate derived categories of étale sheaves

${\displaystyle \mathrm{Rad}:=R{p}_{2,\ast }{p}_1^{\ast }:D({\mathbf{P}}^d)\to D({\mathbf{P}}^{\vee ,d}).}$

The main theorem about this transform is that this transform induces an equivalence of the categories of perverse sheaves on the projective space and its dual projective space, up to constant sheaves.