Category Archives: Optics and waves

Diffractionless Tomography

Classical tomography is based on the assumption that the illuminating light travels in straight lines through the object. While this assumption is only valid for waves of very short wavelength such as, for example, x-rays, it has also been extensively used as an approximation for the light in the visible spectrum.

Radon Transform

The Radon transform describes the relationship between a D-dimensional function and its projection onto a D-1 dimensional hyperplane. In practice, we are often interested in the cases when D = 2 or D = 3, which correspond to the projection of a two-dimensional (2D) surface onto a line or of a three-dimensional (3D) volume onto a plane, respectively. Generally, it is sufficient to consider the problem in 2D, since, the 3D case can be solved by decomposing it into a succession of 2D slices.

Consider the tomographic problem illustrated in the following figure.


Left: Schematic representation of the tomographic scenario. Right: Sinogram of the Shepp-Logan phantom
Left: Schematic representation of the tomographic scenario. Right: Sinogram of the Shepp-Logan phantom

Let \Omega \in \mathbb{R}^2 be a bounded, 2D domain which contains an object characterized by a function f(\boldsymbol{x}) with \boldsymbol{x} = (x, y). We assume that f varies within \Omega, but is zero outside of it. We probe the object by illuminating it with parallel, indefinitely thin beams of light that travel in the direction making an angle \theta \in [0, \pi] with respect to the y-axis. We collect measurements along the t-axis, orthogonal to the direction of light propagation. Our objective is to reconstruct the object f given the measurements p_\theta(t) for \theta \in [0, \pi] and t \in \mathbb{R}. This measurement problem can be described mathematically as

p_\theta(t) = \mathcal{R}\{f\}(t, \theta) =\int_{\mathbb{R}^2} f(\boldsymbol{x}) \, \delta (\langle \boldsymbol{x}, \boldsymbol{\theta} \rangle - t ) \mathrm{d} \boldsymbol{x},

where \delta is the Dirac distribution, \boldsymbol{\theta} = (\cos \theta, \sin \theta) is a unit vector.  We can rewrite the equation by using the following variable change

\begin{pmatrix}x\\y\end{pmatrix} =\begin{bmatrix}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{bmatrix}\begin{pmatrix}t\\s\end{pmatrix}


p_\theta(t) =\int_\mathbb{R} f(t \cos \theta - s \sin \theta, t \sin \theta + s \cos \theta) \mathrm{d} s.

The operator \mathcal{R}: L_2(\mathbb{R}^2) \rightarrow L_2(\mathbb{R} \times [0, \pi]) is called the Radon transform after Johann Radon (1887-1956) who introduced it in 1917. The data \{p_\theta(t)\}_{\theta \in [0, \pi], t \in \mathbb{R}} generated by the transform is often called the sinogram, since it maps a point into a sinusoid. To see this consider a point \boldsymbol{x}_0 = (x_0, y_0), whose projection onto the t-axis will be t_0(\theta) = x_0 \cos \theta + y_0 \sin \theta, which can be represented in polar coordinates (r, \phi), with x_0 = r \cos \phi and y_0 = r \sin \phi, as t_0(\theta) = r \cos (\theta-\phi). Hence, the trajectory of the point \boldsymbol{x}_0 is a cosine of radius r and shift \phi.

Fourier Slice Theorem

Fourier slice theorem is one of the most important results from classical tomography, which establishes the connection between the data p_\theta(t) and the function f(\boldsymbol{x}) in Fourier space.

Theorem: The 1D Fourier transform of the projected data \hat{p}_\theta(\omega) = \mathcal{F}_{\text{1D}}\{p_\theta(t)\} corresponds to a slice of the 2D Fourier transform of the function \hat{f}(\boldsymbol{\omega}) = \mathcal{F}_{\text{2D}}\{f(\boldsymbol{x})\}, with \boldsymbol{\omega} = (\omega_x, \omega_y), taken along the angle \theta through the origin, or more concisely

\hat{p}_\theta(\omega) = \hat{f}(\omega \cos \theta, \omega \sin \theta).

To establish this result, we consider the Fourier transform

\hat{p}_\theta(\omega) &= \int_\mathbb{R} p_\theta(t) \mathrm{e}^{-\mathrm{j} \omega t} \mathrm{d} t \\
&= \int_\mathbb{R} \left [\int_{\mathbb{R}^2} f(\boldsymbol{x}) \delta(\langle \boldsymbol{x}, \boldsymbol{\theta} \rangle - t) \mathrm{d} \boldsymbol{x} \right] \mathrm{e}^{-\mathrm{j} \omega t} \mathrm{d} t \\
&= \int_{\mathbb{R}^2} f(\boldsymbol{x}) \mathrm{e}^{- \mathrm{j} \omega \langle \boldsymbol{x}, \boldsymbol{\theta} \rangle } \mathrm{d} \boldsymbol{x} \\
&= \hat{f}(\omega \cos \theta, \omega \sin \theta),

where in the second equality use used the definition of the Radon transform.

The remarkable aspect of the Fourier Slice Theorem is that it allows one to invert the Radon transform by (a) collecting measurements for all \theta \in [0, \pi], (b) evaluating the 1D Fourier transform of these measurements and using them to fill the 2D Fourier space, and (c) evaluating inverse 2D Fourier transform to form the image.


We now define the backprojection operator \mathcal{R}^\ast: L_2(\mathbb{R} \times [0, \pi]) \rightarrow L_2(\mathbb{R}^2)

b(\boldsymbol{x}) = \mathcal{R}^\ast\{p_\theta(t)\}(\boldsymbol{x}) = \int_0^\pi p_\theta (x \cos \theta + y \sin \theta) \mathrm{d} \theta

which is in fact the adjoint of the Radon transform, i.e., for all p \in L_2(\mathbb{R} \times [0, \pi]) and for all f \in L_2(\mathbb{R}^2), we have that

\langle p, \mathcal{R} f\rangle_{L_2(\mathbb{R} \times [0, \pi])} = \langle \mathcal{R}^\ast p, f\rangle_{L_2(\mathbb{R}^2)}.

Intuitively, for each point \boldsymbol{x}, backprojection collects all the projections that passed through that point.

Filtered Backprojection

Consider the following inverse Fourier transform

&= \frac{1}{(2\pi)^2} \int_{\mathbb{R}^2} \hat{f}(\boldsymbol{\omega}) \mathrm{e}^{\mathrm{j} \langle \boldsymbol{x}, \boldsymbol{\omega} \rangle} \mathrm{d} \boldsymbol{\omega} \\
&= \frac{1}{(2\pi)^2} \int_0^\pi \int_\mathbb{R} \hat{f}(\omega \cos \theta, \omega \sin \theta) \mathrm{e}^{\mathrm{j} \omega \langle \boldsymbol{x}, \boldsymbol{\theta} \rangle} |\omega| \mathrm{d} \omega \mathrm{d} \theta\\
&= \frac{1}{(2\pi)^2} \int_0^\pi \int_\mathbb{R} |\omega| \hat{p}_\theta(\omega) \mathrm{e}^{\mathrm{j} \omega \langle \boldsymbol{x}, \boldsymbol{\theta} \rangle} \mathrm{d} \omega \mathrm{d} \theta \\
&= \int_0^\pi \left[\frac{1}{2\pi}\int_\mathbb{R} \frac{|\omega|}{2\pi} \, \hat{p}_\theta(\omega) \mathrm{e}^{\mathrm{j} \omega \langle \boldsymbol{x}, \boldsymbol{\theta} \rangle} \mathrm{d} \omega\right] \mathrm{d} \theta \\
&= \int_0^\pi q_\theta(x \cos \theta + y \sin \theta) \mathrm{d} \theta
where we performed the variable change \omega_x = \omega \cos \theta and \omega_y = \omega \sin \theta, which gives the Jacobian J(\omega, \theta) = \omega. This implies that it is possible to reconstruct the function f by 1D filtering of the Radon data with

h(t) = \mathcal{F}^{-1}\left\{\frac{|\omega|}{2\pi}\right\}(t)

and then applying the backprojection operator \mathcal{R}^\ast to the result q_\theta(t) = (h \ast p_\theta)(t). To read the PDF version of this post click here.

From Maxwell to Helmholtz

An accurate physical model is one of the main contributing factors towards the quality of  images that can be formed with an optical tomographic microscope. Our model in [1] was based on the beam propagation method, which splits the sample into layers, and computationally propagates the light layer-by-layer from the source all the way to the detector. The starting point for our derivations was the scalar Helmholtz equation

\left(\Delta + k^2(\boldsymbol{x}) \mathrm{I}\right) u(\boldsymbol{x}) = 0,

where \boldsymbol{x} = (x, y, z) denotes a spatial coordinate, u is what we called the total light-field at \boldsymbol{x}, \Delta \triangleq (\partial^2/\partial x^2 + \partial^2/\partial y^2 + \partial^2/\partial z^2) is the Laplacian, \mathrm{I} is the identity operator, and k(\boldsymbol{x}) is the wavenumber of the field at \boldsymbol{x}. The spatial dependence of the wavenumber k is due to variations of the speed of light c induced by the inhomogeneous nature of the sample under consideration. Once we have the Helmholtz equation, we can repeat the process described in [1], and obtain an algorithm for computationally forming an image from the measured data.

However, how does the Helmholtz equation above relate to the actual Maxwell's equations?

Maxwell's equations provide the full description of the optical waves, and there are four equations to consider. The first equation is often called the Faraday's law of induction

\nabla \times \mathbf{E}(\boldsymbol{x}, t) = -\frac{\partial \mathbf{B}(\boldsymbol{x},t)}{\partial t},

where \mathbf{E} and \mathbf{B} are electric and magnetic fields, respectively, and \nabla \times is the curl operator. The second equation is sometimes called the Ampère's circuital law

\nabla \times \mathbf{H}(\boldsymbol{x},t) = \mathbf{J}(\boldsymbol{x}, t) + \frac{\partial \mathbf{D}(\boldsymbol{x},t)}{\partial t},

where \mathbf{J} is the is the total current density, \mathbf{H} is a magnetic field, and \mathbf{D} is the electric displacement field. The total current density \mathbf{J} is related to the electric field \mathbf{E} as \mathbf{J} = \sigma \mathbf{E}, where \sigma(\boldsymbol{x}) is the distribution of condictivity in the sample. The magnetic field \mathbf{H} is related to \mathbf{B} as \mathbf{B} = \mu \mathbf{H}, where \mu(\boldsymbol{x}) is the distribution of permeability in space. Similarly, the electric displacement field \mathbf{D} is related to \mathbf{E} as \mathbf{D} = \epsilon \mathbf{E}, where \epsilon(\boldsymbol{x}) is the distribution of permittivity in space. The third equation is the Gauss's law

\nabla \cdot \mathbf{D}(\boldsymbol{x},t) = \rho(\boldsymbol{x}),

where \nabla \cdot is divergence operator and \rho is the distribution of charge in the sample. Finally, the fourth equation is the Gauss's law for magnetism

 \nabla \cdot \mathbf{B}(\boldsymbol{x},t) = 0.

To obtain the scalar Helmhotz equation, we assume that the sample is charge free and non-magnetic, i.e., \rho(\boldsymbol{x}) = 0 and \mu(\boldsymbol{x}) = \mu_0. We also limit ourselves to the electric fields that have transverse magnetic (TM) polarizations, i.e., \mathbf{E}(\boldsymbol{x},t) = (0, 0, E_z(x,y,t)) and \mathbf{H}(\boldsymbol{x},t) = (H_x(\boldsymbol{x},t), H_y(\boldsymbol{x},t), 0), and that are thus divergence free \nabla \cdot \mathbf{E} = 0. Then by applying the identity \nabla \times (\nabla \times \mathbf{E}) = \nabla (\nabla \cdot \mathbf{E}) - \Delta \mathbf{E}, and limiting ourselves to the scalar E_z, we obtain

\Delta E(\boldsymbol{x},t) = \mu_0\sigma(\boldsymbol{x}) \frac{\partial E(\boldsymbol{x},t)}{\partial t} + \mu_0\epsilon(\boldsymbol{x})\frac{\partial^2 E(\boldsymbol{x},t)}{\partial t^2}

The usual Helmholtz equation is then obtained by taking a temporal Fourier transform form both sides, ignoring \omega, and rearranging the terms

\Delta E(\boldsymbol{x}) + k^2(\boldsymbol{x}) E(\boldsymbol{x}) =0,

where k^2 = \omega^2 \mu_0 \epsilon(1 - \mathrm{j}\frac{\sigma}{\omega \epsilon}). Typically, it is more convenient to write k^2 = k_0^2 \epsilon, where k_0^2 = \omega^2 \epsilon_0 \mu_0 corresponds to the wave vector in the air, and \epsilon has been modified to incorporate \sigma, and is thus complex.



[1] U. S. Kamilov, I. N. Papadopoulos, M. H. Shoreh, A. Goy, C. Vonesch, M. Unser, and D. Psaltis, “Optical tomographic image reconstruction based on beam propagation and sparse regularization,” IEEE Trans. Comp. Imag., vol. 2, no. 1, 2016,