3D inversion

The first Born approximation for a 3D 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}) = \iiint \!\! d^3r' G(\mathbf{r-r'}) f(\mathbf{r'}) u_0(\mathbf{r'})\]

The Green’s function in 3D can be written as:

\[G(\mathbf{r-r'}) = \frac{ik_\mathrm{m}}{8\pi^2} \iint \!\! dpdq \frac{1}{M} \exp\! \left \lbrace i k_\mathrm{m} \left[ p(x-x') + q(y-y') + M(z-z') \right] \right \rbrace\]

with

\[M = \sqrt{1-p^2-q^2}\]

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

\[\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}, k_\mathrm{Dy}) \exp \! \left(-i k_\mathrm{m} M l_\mathrm{D} \right)\]

where \(\widehat{F}(k_\mathrm{x}, k_\mathrm{y}, k_\mathrm{z})\) is the Fourier transformed object function and \(\widehat{U}_{\mathrm{B}, \phi_0}(k_\mathrm{Dx}, k_\mathrm{Dy})\) 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} &= (p, q, M)\\\mathbf{s_0} &= (p_0, q_0, M_0) = (-\sin\phi_0, \, 0, \, \cos\phi_0)\\\mathbf{t_\perp} &= \left(\cos\phi_0, \, \frac{k_\mathrm{Dy}}{k_\mathrm{Dx}}, \, \sin\phi_0 \right)^\top\end{aligned}\end{align} \]

Method summary

backpropagate_3d(uSin, angles, res, nm[, …]) 3D backpropagation
backpropagate_3d_tilted(uSin, angles, res, nm) 3D backpropagation with a tilted axis of rotation

Backpropagation

odtbrain.backpropagate_3d(uSin, angles, res, nm, lD=0, coords=None, weight_angles=True, onlyreal=False, padding=(True, True), padfac=1.75, padval='edge', intp_order=2, dtype=None, num_cores=2, save_memory=False, copy=True, count=None, max_count=None, verbose=0)

3D backpropagation

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

This method implements the 3D backpropagation algorithm [Mueller2015arxiv].

\[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{2D}} \left \{ \left| k_\mathrm{Dx} \right| \frac{\text{FFT}_{\mathrm{2D}} \left \{ u_{\mathrm{B},\phi_j}(x_\mathrm{D}, y_\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{2D}}\) and inverse \(\text{FFT}^{-1}_{\mathrm{2D}}\) 2D 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, Ny, Nx) ndarray) – Three-dimensional sinogram of plane recordings \(u_{\mathrm{B}, \phi_j}(x_\mathrm{D}, y_\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 [(3, M) ndarray]) – Only compute 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 (tuple of 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. The default is padding in x and y: padding=(True, True). For padding only in x-direction (e.g. for cylindrical symmetries), set padding to (True, False). To turn off padding, set it to (False, False).
  • padfac (float) – Increase padding size of the input data. A value greater than one will trigger padding to the second-next power of two. For example, a value of 1.75 will lead to a padded size of 256 for an initial size of 144, whereas it will lead to a padded size of 512 for an initial size of 150. Values geater than 2 are allowed. This parameter may greatly increase memory usage!
  • padval (float or "edge") – The value used for padding. This is important for the Rytov approximation, where an approximat 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 “edge”, then the edge values are used for padding (see documentation of numpy.pad()). If padval is a float, then padding is done with a linear ramp.
  • intp_order (int between 0 and 5) – Order of the interpolation for rotation. See scipy.ndimage.interpolation.rotate() for details.
  • dtype (dtype object or argument for numpy.dtype()) – The data type that is used for calculations (float or double). Defaults to numpy.float_.
  • num_cores (int) – The number of cores to use for parallel operations. This value defaults to the number of cores on the system.
  • save_memory (bool) –

    Saves memory at the cost of longer computation time.

    New in version 0.1.5.

  • copy (bool) –

    Copy input sinogram uSin for data processing. If copy is set to False, then uSin will be overridden.

    New in version 0.1.5.

  • 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 (Nx, Ny, Nx), complex if onlyreal==False

See also

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).

Backpropagation with tilted axis of rotation

odtbrain.backpropagate_3d_tilted(uSin, angles, res, nm, lD=0, tilted_axis=[0, 1, 0], coords=None, weight_angles=True, onlyreal=False, padding=(True, True), padfac=1.75, padval='edge', intp_order=2, dtype=None, num_cores=2, save_memory=False, copy=True, count=None, max_count=None, verbose=0)

3D backpropagation with a tilted axis of rotation

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

This method implements the 3D backpropagation algorithm with a rotational axis that is tilted by \(\theta_\mathrm{tilt}\) w.r.t. the imaging plane [Mueller2015tilted].

\[f(\mathbf{r}) = -\frac{i k_\mathrm{m}}{2\pi} \sum_{j=1}^{N} \! \Delta \phi_0 D_{-\phi_j}^\mathrm{tilt} \!\! \left \{ \text{FFT}^{-1}_{\mathrm{2D}} \left \{ \left| k_\mathrm{Dx} \cdot \cos \theta_\mathrm{tilt}\right| \frac{\text{FFT}_{\mathrm{2D}} \left \{ u_{\mathrm{B},\phi_j}(x_\mathrm{D}, y_\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 a modified rotational operator \(D_{-\phi_j}^\mathrm{tilt}\) and a different filter in Fourier space \(|k_\mathrm{Dx} \cdot \cos \theta_\mathrm{tilt}|\) when compared to backpropagate_3d().

New in version 0.1.2.

Parameters:
  • uSin ((A, Ny, Nx) ndarray) – Three-dimensional sinogram of plane recordings \(u_{\mathrm{B}, \phi_j}(x_\mathrm{D}, y_\mathrm{D})\) divided by the incident plane wave \(u_0(l_\mathrm{D})\) measured at the detector.
  • angles (ndarray of shape (A,3) or 1D array of length A) – If the shape is (A,3), then angles consists of vectors on the unit sphere that correspond to the direction of illumination and acquisition (s₀). If the shape is (A,), then angles is a one-dimensional array of angles in radians that determines the angular position \(\phi_j\). In both cases, tilted_axis must be set according to the tilt of the rotational axis.
  • 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.
  • tilted_axis (list of floats) – The coordinates [x, y, z] on a unit sphere representing the tilted axis of rotation. The default is (0,1,0), which corresponds to a rotation about the y-axis and follows the behavior of odtbrain.backpropagate_3d().
  • coords (None [(3, M) ndarray]) – Only compute 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}\]

    This currently only works when angles has the shape (A,).

  • onlyreal (bool) – If True, only the real part of the reconstructed image will be returned. This saves computation time.
  • padding (tuple of 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. The default is padding in x and y: padding=(True, True). For padding only in x-direction (e.g. for cylindrical symmetries), set padding to (True, False). To turn off padding, set it to (False, False).
  • padfac (float) – Increase padding size of the input data. A value greater than one will trigger padding to the second-next power of two. For example, a value of 1.75 will lead to a padded size of 256 for an initial size of 144, whereas it will lead to a padded size of 512 for an initial size of 150. Values greater than 2 are allowed. This parameter may greatly increase memory usage!
  • padval (float or "edge") – 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 “edge”, then the edge values are used for padding (see documentation of numpy.pad()). If padval is a float, then padding is done with a linear ramp.
  • intp_order (int between 0 and 5) – Order of the interpolation for rotation. See scipy.ndimage.interpolation.affine_transform() for details.
  • dtype (dtype object or argument for numpy.dtype()) – The data type that is used for calculations (float or double). Defaults to numpy.float_.
  • num_cores (int) – The number of cores to use for parallel operations. This value defaults to the number of cores on the system.
  • save_memory (bool) –

    Saves memory at the cost of longer computation time.

    New in version 0.1.5.

  • copy (bool) –

    Copy input sinogram uSin for data processing. If copy is set to False, then uSin will be overridden.

    New in version 0.1.5.

  • 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 (Nx, Ny, Nx), complex if onlyreal==False

See also

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

Notes

This implementation can deal with projection angles that are not distributed along a circle about the rotational axis. If there are slight deviations from this circle, simply pass the 3D rotational positions instead of the 1D angles to the angles argument. In principle, this should improve the reconstruction. The general problem here is that the backpropagation algorithm requires a ramp filter in Fourier space that is oriented perpendicular to the rotational axis. If the sample does not rotate about a single axis, then a 1D parametric representation of this rotation must be found to correctly determine the filter in Fourier space. Such a parametric representation could e.g. be a spiral between the poles of the unit sphere (but this kind of rotation is probably difficult to implement experimentally).

If you have input images with rectangular shape, e.g. Nx!=Ny and the rotational axis deviates by approximately PI/2 from the axis (0,1,0), then data might get cropped in the reconstruction volume. You can avoid that by rotating your input data and the rotational axis by PI/2. For instance, change`tilted_axis` from [1,0,0] to [0,1,0] and np.rot90 the sinogram images.

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).