2D inversion

The first Born approximation for a 2D scattering problem with a plane wave \(u_0(\mathbf{r}) = a_0 \exp(-ik_\mathrm{m}\mathbf{s_0r})\) reads:

\[u_\mathrm{B}(\mathbf{r}) = \iint \!\! d^2r' G(\mathbf{r-r'}) f(\mathbf{r'}) u_0(\mathbf{r'})\]

The Green’s function in 2D is the zero-order Hankel function of the first kind:

\[G(\mathbf{r-r'}) = \frac{i}{4} H_0^\mathrm{(1)}(k_\mathrm{m} \left| \mathbf{r-r'} \right|)\]

Solving for \(f(\mathbf{r})\) yields the Fourier diffraction theorem in 2D

\[\widehat{F}(k_\mathrm{m}(\mathbf{s-s_0})) = - \sqrt{\frac{2}{\pi}} \frac{i k_\mathrm{m}}{a_0} M \widehat{U}_{\mathrm{B},\phi_0}(k_\mathrm{Dx}) \exp \! \left(-i k_\mathrm{m} M l_\mathrm{D} \right)\]

where \(\widehat{F}(k_\mathrm{x}, k_\mathrm{z})\) is the Fourier transformed object function and \(\widehat{U}_{\mathrm{B}, \phi_0}(k_\mathrm{Dx})\) is the Fourier transformed complex wave that travels along \(\mathbf{s_0}\) (in the direction of \(\phi_0\)) measured at the detector \(\mathbf{r_D}\).

The following identities are used:

\[ \begin{align}\begin{aligned}k_\mathrm{m} (\mathbf{s-s_0}) &= k_\mathrm{Dx} \, \mathbf{t_\perp} + k_\mathrm{m}(M - 1) \, \mathbf{s_0}\\\mathbf{s_0} &= \left(p_0 , \, M_0 \right) = (-\sin\phi_0, \, \cos\phi_0)\\\mathbf{t_\perp} &= \left(- M_0 , \, p_0 \right) = (\cos\phi_0, \, \sin\phi_0)\end{aligned}\end{align} \]

Method summary

backpropagate_2d(uSin, angles, res, nm[, …]) 2D backpropagation with the Fourier diffraction theorem
fourier_map_2d(uSin, angles, res, nm[, lD, …]) 2D Fourier mapping with the Fourier diffraction theorem
integrate_2d(uSin, angles, res, nm[, lD, …]) (slow) 2D reconstruction with the Fourier diffraction theorem

Backpropagation

odtbrain.backpropagate_2d(uSin, angles, res, nm, lD=0, coords=None, weight_angles=True, onlyreal=False, padding=True, padval=0, count=None, max_count=None, verbose=0)

2D backpropagation with the Fourier diffraction theorem

Two-dimensional diffraction tomography reconstruction algorithm for scattering of a plane wave \(u_0(\mathbf{r}) = u_0(x,z)\) by a dielectric object with refractive index \(n(x,z)\).

This method implements the 2D backpropagation algorithm [MSG15b].

\[f(\mathbf{r}) = -\frac{i k_\mathrm{m}}{2\pi} \sum_{j=1}^{N} \! \Delta \phi_0 D_{-\phi_j} \!\! \left \{ \text{FFT}^{-1}_{\mathrm{1D}} \left \{ \left| k_\mathrm{Dx} \right| \frac{\text{FFT}_{\mathrm{1D}} \left \{ u_{\mathrm{B},\phi_j}(x_\mathrm{D}) \right \} }{u_0(l_\mathrm{D})} \exp \! \left[i k_\mathrm{m}(M - 1) \cdot (z_{\phi_j}-l_\mathrm{D}) \right] \right \} \right \}\]

with the forward \(\text{FFT}_{\mathrm{1D}}\) and inverse \(\text{FFT}^{-1}_{\mathrm{1D}}\) 1D fast Fourier transform, the rotational operator \(D_{-\phi_j}\), the angular distance between the projections \(\Delta \phi_0\), the ramp filter in Fourier space \(|k_\mathrm{Dx}|\), and the propagation distance \((z_{\phi_j}-l_\mathrm{D})\).

Parameters:
  • uSin ((A,N) ndarray) – Two-dimensional sinogram of line recordings \(u_{\mathrm{B}, \phi_j}(x_\mathrm{D})\) divided by the incident plane wave \(u_0(l_\mathrm{D})\) measured at the detector.
  • angles ((A,) ndarray) – Angular positions \(\phi_j\) of uSin in radians.
  • res (float) – Vacuum wavelength of the light \(\lambda\) in pixels.
  • nm (float) – Refractive index of the surrounding medium \(n_\mathrm{m}\).
  • lD (float) – Distance from center of rotation to detector plane \(l_\mathrm{D}\) in pixels.
  • coords (None [(2,M) ndarray]) – Computes only the output image at these coordinates. This keyword is reserved for future versions and is not implemented yet.
  • weight_angles (bool) –

    If True, weights each backpropagated projection with a factor proportional to the angular distance between the neighboring projections.

    \[\Delta \phi_0 \longmapsto \Delta \phi_j = \frac{\phi_{j+1} - \phi_{j-1}}{2}\]

    New in version 0.1.1.

  • onlyreal (bool) – If True, only the real part of the reconstructed image will be returned. This saves computation time.
  • padding (bool) – Pad the input data to the second next power of 2 before Fourier transforming. This reduces artifacts and speeds up the process for input image sizes that are not powers of 2.
  • padval (float) – The value used for padding. This is important for the Rytov approximation, where an approximate zero in the phase might translate to 2πi due to the unwrapping algorithm. In that case, this value should be a multiple of 2πi. If padval is None, then the edge values are used for padding (see documentation of numpy.pad()).
  • max_count (count,) – Can be used to monitor the progress of the algorithm. Initially, the value of max_count.value is incremented by the total number of steps. At each step, the value of count.value is incremented.
  • verbose (int) – Increment to increase verbosity.
Returns:

f – Reconstructed object function \(f(\mathbf{r})\) as defined by the Helmholtz equation. \(f(x,z) = k_m^2 \left(\left(\frac{n(x,z)}{n_m}\right)^2 -1\right)\)

Return type:

ndarray of shape (N,N), complex if onlyreal is False

See also

odt_to_ri()
conversion of the object function \(f(\mathbf{r})\) to refractive index \(n(\mathbf{r})\)
radontea.backproject()
backprojection based on the Fourier slice theorem

Notes

Do not use the parameter lD in combination with the Rytov approximation - the propagation is not correctly described. Instead, numerically refocus the sinogram prior to converting it to Rytov data (using e.g. odtbrain.sinogram_as_rytov()) with a numerical focusing algorithm (available in the Python package nrefocus).

Fourier mapping

odtbrain.fourier_map_2d(uSin, angles, res, nm, lD=0, semi_coverage=False, coords=None, count=None, max_count=None, verbose=0)

2D Fourier mapping with the Fourier diffraction theorem

Two-dimensional diffraction tomography reconstruction algorithm for scattering of a plane wave \(u_0(\mathbf{r}) = u_0(x,z)\) by a dielectric object with refractive index \(n(x,z)\).

This function implements the solution by interpolation in Fourier space.

Parameters:
  • uSin ((A,N) ndarray) – Two-dimensional sinogram of line recordings \(u_{\mathrm{B}, \phi_j}(x_\mathrm{D})\) divided by the incident plane wave \(u_0(l_\mathrm{D})\) measured at the detector.
  • angles ((A,) ndarray) – Angular positions \(\phi_j\) of uSin in radians.
  • res (float) – Vacuum wavelength of the light \(\lambda\) in pixels.
  • nm (float) – Refractive index of the surrounding medium \(n_\mathrm{m}\).
  • lD (float) – Distance from center of rotation to detector plane \(l_\mathrm{D}\) in pixels.
  • semi_coverage (bool) – If set to True, it is assumed that the sinogram does not necessarily cover the full angular range from 0 to 2π, but an equidistant coverage over 2π can be achieved by inferring point (anti)symmetry of the (imaginary) real parts of the Fourier transform of f. Valid for any set of angles {X} that result in a 2π coverage with the union set {X}U{X+π}.
  • coords (None [(2,M) ndarray]) – Computes only the output image at these coordinates. This keyword is reserved for future versions and is not implemented yet.
  • max_count (count,) – Can be used to monitor the progress of the algorithm. Initially, the value of max_count.value is incremented by the total number of steps. At each step, the value of count.value is incremented.
  • verbose (int) – Increment to increase verbosity.
Returns:

f – Reconstructed object function \(f(\mathbf{r})\) as defined by the Helmholtz equation. \(f(x,z) = k_m^2 \left(\left(\frac{n(x,z)}{n_m}\right)^2 -1\right)\)

Return type:

ndarray of shape (N,N), complex if onlyreal is False

See also

backpropagate_2d()
implementation by backpropagation
odt_to_ri()
conversion of the object function \(f(\mathbf{r})\) to refractive index \(n(\mathbf{r})\)

Notes

Do not use the parameter lD in combination with the Rytov approximation - the propagation is not correctly described. Instead, numerically refocus the sinogram prior to converting it to Rytov data (using e.g. odtbrain.sinogram_as_rytov()) with a numerical focusing algorithm (available in the Python package nrefocus).

The interpolation in Fourier space (which is done with scipy.interpolate.griddata()) may be unstable and lead to artifacts if the data to interpolate contains sharp spikes. This issue is not handled at all by this method (in fact, a test has been removed in version 0.2.6 because griddata gave different results on Windows and Linux).

Direct sum

odtbrain.integrate_2d(uSin, angles, res, nm, lD=0, coords=None, count=None, max_count=None, verbose=0)

(slow) 2D reconstruction with the Fourier diffraction theorem

Two-dimensional diffraction tomography reconstruction algorithm for scattering of a plane wave \(u_0(\mathbf{r}) = u_0(x,z)\) by a dielectric object with refractive index \(n(x,z)\).

This function implements the solution by summation in real space, which is extremely slow.

Parameters:
  • uSin ((A,N) ndarray) – Two-dimensional sinogram of line recordings \(u_{\mathrm{B}, \phi_j}(x_\mathrm{D})\) divided by the incident plane wave \(u_0(l_\mathrm{D})\) measured at the detector.
  • angles ((A,) ndarray) – Angular positions \(\phi_j\) of uSin in radians.
  • res (float) – Vacuum wavelength of the light \(\lambda\) in pixels.
  • nm (float) – Refractive index of the surrounding medium \(n_\mathrm{m}\).
  • lD (float) – Distance from center of rotation to detector plane \(l_\mathrm{D}\) in pixels.
  • coords (None or (2,M) ndarray]) – Computes only the output image at these coordinates. This keyword is reserved for future versions and is not implemented yet.
  • max_count (count,) – Can be used to monitor the progress of the algorithm. Initially, the value of max_count.value is incremented by the total number of steps. At each step, the value of count.value is incremented.
  • verbose (int) – Increment to increase verbosity.
Returns:

f – Reconstructed object function \(f(\mathbf{r})\) as defined by the Helmholtz equation. \(f(x,z) = k_m^2 \left(\left(\frac{n(x,z)}{n_m}\right)^2 -1\right)\)

Return type:

ndarray of shape (N,N), complex if onlyreal is False

See also

backpropagate_2d()
implementation by backprojection
fourier_map_2d()
implementation by Fourier interpolation
odt_to_ri()
conversion of the object function \(f(\mathbf{r})\) to refractive index \(n(\mathbf{r})\)

Notes

This method is not meant for production use. The computation time is very long and the reconstruction quality is bad. This function is included in the package, because of its educational value, exemplifying the backpropagation algorithm.

Do not use the parameter lD in combination with the Rytov approximation - the propagation is not correctly described. Instead, numerically refocus the sinogram prior to converting it to Rytov data (using e.g. odtbrain.sinogram_as_rytov()) with a numerical focusing algorithm (available in the Python package nrefocus).