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:
The Green’s function in 3D can be written as:
with
Solving for \(f(\mathbf{r})\) yields the Fourier diffraction theorem in 3D
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:
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 packagenrefocus
).
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 packagenrefocus
).