Code reference¶
Parallel beam geometry¶
radon_parallel (arr, angles[, count, max_count]) |
Compute the Radon transform (sinogram) of a circular image. |
backproject (sinogram, angles[, filtering, …]) |
2D backprojection algorithm |
fourier_map (sinogram, angles[, intp_method, …]) |
2D Fourier mapping with the Fourier slice theorem |
integrate (sinogram, angles[, count, max_count]) |
2D sum-reconstruction with the Fourier slice theorem |
art (sinogram, angles[, initial, iterations, …]) |
Algebraic Reconstruction Technique |
sart (sinogram, angles[, initial, …]) |
Simultaneous Algebraic Reconstruction Technique |
Radon transform¶
-
radontea.
radon_parallel
(arr, angles, count=None, max_count=None)¶ Compute the Radon transform (sinogram) of a circular image.
The
scipy
Radon transform performs this operation on the entire image, whereas this implementation requires an input image that has gray-scale values of 0 outside of a circle with diameter equal to the image size.Parameters: - arr (ndarray, shape (N,N)) – the input image.
- angles (ndarray, length A) – angles or projections in radians
- 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.
Returns: outarr – Sinogram of the input image. The i’th row contains the projection data of the i’th angle.
Return type: ndarray of floats, shape (A, N)
See also
scipy.ndimage.interpolation.rotate()
- The interpolator used to rotate the image.
Non-iterative reconstruction¶
Computes the inverse Radon transform with non-iterative techniques.
The linear system of equations that describes the forward process can be
inverted with several algorithms, most notably the backprojection
algorithm radontea.backproject()
. The reconstruction is based on
the Fourier slice theorem. A Fourier-based interpolation algorithm is
implemented in radontea.fourier_map()
.
Backprojection¶
-
radontea.
backproject
(sinogram, angles, filtering='ramp', weight_angles=True, padding=True, padval=0, count=None, max_count=None, verbose=0)¶ 2D backprojection algorithm
Computes the inverse of the Radon transform using filtered backprojection.
Parameters: - sinogram (ndarray, shape (A,N)) – Two-dimensional sinogram of line recordings.
- angles ((A,) ndarray) – Angular positions of the sinogram in radians from zero to PI.
- filtering ({'ramp', 'shepp-logan', 'cosine', 'hamming', 'hann'}, optional) –
Specifies the Fourier filter. Either of
- ”ramp”: mathematically correct reconstruction
- ”shepp-logan”
- ”cosine”
- ”hamming”
- ”hann”
- 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.9.
- padding (bool, optional) – 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. 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: out – The reconstructed image.
Return type: ndarray
Examples
>>> import numpy as np >>> from radontea import backproject >>> N = 5 >>> A = 7 >>> x = np.linspace(-N/2, N/2, N) >>> projection = np.exp(-x**2) >>> sinogram = np.tile(projection, A).reshape(A,N) >>> angles = np.linspace(0, np.pi, A, endpoint=False) >>> backproject(sinogram, angles) array([[ 0.03813283, -0.01972347, -0.02221885, 0.03822303, 0.01903376], [ 0.04033526, -0.01591647, 0.02262173, 0.04203002, 0.02123619], [ 0.02658797, 0.11375576, 0.41525594, 0.17170226, 0.0074889 ], [ 0.04008672, 0.11139209, 0.27650193, 0.16933859, 0.02098764], [ 0.02140445, 0.0334597 , 0.07691547, 0.0914062 , 0.00230537]])
Fourier mapping¶
-
radontea.
fourier_map
(sinogram, angles, intp_method='cubic', count=None, max_count=None)¶ 2D Fourier mapping with the Fourier slice theorem
Computes the inverse of the Radon transform using Fourier interpolation. Warning: This is the naive reconstruction that assumes that the image is rotated through the upper left pixel nearest to the actual center of the image. We do not have this problem for odd images, only for even images.
Parameters: - sinogram ((A,N) ndarray) – Two-dimensional sinogram of line recordings.
- angles ((A,) ndarray) – Angular positions of the sinogram in radians equally distributed from zero to PI.
- intp_method ({'cubic', 'nearest', 'linear'}, optional) –
Method of interpolation. For more information see scipy.interpolate.griddata. One of
- ”nearest”: instead of interpolating, use the points closest to the input data
- ”linear”: bilinear interpolation between data points
- ”cubic”: interpolate using a two-dimensional poolynimial surface
- 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.
Returns: out – The reconstructed image.
Return type: ndarray
See also
scipy.interpolate.griddata()
- interpolation method used
Slow integration¶
-
radontea.
integrate
(sinogram, angles, count=None, max_count=None)¶ 2D sum-reconstruction with the Fourier slice theorem
Computes the inverse of the Radon transform by computing the integral in real space.
Parameters: - sinogram ((A,N) ndarray) – Two-dimensional sinogram of line recordings.
- angles ((A,) ndarray) – Angular positions of the sinogram in radians equally distributed from zero to PI.
- 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.
Returns: out – The reconstructed image.
Return type: ndarray
Iterative reconstruction¶
Inversion of Radon-based tomography methods using iterative algorithms. The convegence of these algorithms might be slow. The implementation is not optimized.
ART¶
-
radontea.
art
(sinogram, angles, initial=None, iterations=1, count=None, max_count=None)¶ Algebraic Reconstruction Technique
The Algebraic Reconstruction Technique (ART) iteratively computes the inverse of the Radon transform in two dimensions. The reconstruction technique uses rays of the diameter of one pixel to iteratively solve the system of linear equations that describe the projection process. The binary weighting factors are
- 1, if the center of the a pixel is within the ray
- 0, else
Parameters: - sinogram (ndarrayy, shape (A,N)) – Two-dimensional sinogram of line recordings.
- angles (ndarray, length A) – Angular positions of the sinogram in radians. The angles
at which the sinogram slices were recorded do not have to be
distributed equidistantly as in
backproject()
. The angles are internaly converted to modulo PI. - initial (ndarray, shape (N,N), optional) – The initial guess for the solution.
- iterations (int) – Number of iterations to perform.
- 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.
Notes
For theoretical backround, see Kak, A. C., & Slaney, M.. Principles of Computerized Tomographic Imaging, SIAM, (2001)
Sec. 7.2: “ART reconstrutions usually suffer from salt and pepper noise, which is caused by the inconsitencies introuced in the set of equations by the approximations commonly used for \(w_{ik}\) ‘s.”
SART¶
-
radontea.
sart
(sinogram, angles, initial=None, iterations=1, count=None, max_count=None)¶ Simultaneous Algebraic Reconstruction Technique
SART computes an inverse of the Radon transform in two dimensions. The reconstruction technique uses “rays” of the diameter of one pixel to iteratively solve the system of linear equations that describe the image. The weighting factors are bilinear elements. At the beginning and end of each ray, only partial weights are used. The pixel values of the image are updated only after each iteration is complete.
Parameters: - sinogram (ndarrayy, shape (A,N)) – Two-dimensional sinogram of line recordings.
- angles (ndarray, length A) – Angular positions of the sinogram in radians. The angles at which the sinogram slices were recorded do not have to be distributed equidistantly as in backprojection techniques. The angles are internaly converted to modulo PI.
- initial (ndarray, shape (N,N), optional) – the initial guess for the solution.
- iterations (integer) – Number of iterations to perform.
- 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.
Notes
- Algebraic reconstruction technique (ART) (see art):
- Iterations are performed over each ray of each projection. Weighting factors are binary (1 if center of pixel is within ray, 0 else). This leads to salt and pepper noise.
- Simultaneous iterative reconstruction technique (SIRT):
- Same idea as ART, but for each iteration, the change of the image f is computed for all rays and projections separately and the weights are applied simultaneously after each iteration. The result is a slower convergence but the final image is also less noisy.
This implementation does NOT use a hamming window to filter the data and to emphasize points at the center of the recon- struction region.
For theoretical backround, see Kak, A. C., & Slaney, M.. Principles of Computerized Tomographic Imaging, SIAM, (2001)
Sec 7.4: “[SART] seems to combine the best of ART and SIRT. […] Here are the main features if SART: First, […] the traditional pixel basis is abandonded in favor of bilinear elements [e.g. interpolation]. Also, for a circular reconstruction region, only partial weights are assigned to the first and last picture elements on the individual rays. To further reduce the noise […], the correction terms are simultaneously applied for all the rays in one projection […].”
Fan beam geometry¶
fan.radon_fan (arr, det_size[, det_spacing, …]) |
Compute the Radon transform for a fan beam geometry |
fan.get_det_coords (size, spacing) |
Compute pixel-center positions for 2D detector |
fan.get_fan_coords (size, spacing, distance, …) |
Compute equispaced angular coordinates for 2d detector |
fan.lino2sino (linogram, lds[, stepsize, …]) |
Convert linogram to sinogram for an equispaced detector. |
fan.fan_rec (linogram, lds, method[, …]) |
2D synthetic aperture reconstruction |
Coordinate transforms¶
-
radontea.fan.
get_det_coords
(size, spacing)¶ Compute pixel-center positions for 2D detector
The centers of the pixels of a detector are usually not aligned to a pixel grid. If we center the detector at the origin, odd images will have a pixel at zero, whereas even images will have two pixels next to the actual center.
Parameters: Returns: arr – Pixel coordinates.
Return type: 1D ndarray
-
radontea.fan.
get_fan_coords
(size, spacing, distance, numang)¶ Compute equispaced angular coordinates for 2d detector
The centers of the pixels of a detector are usually not aligned to a pixel grid. If we center the detector at the origin, odd images will have a pixel at zero, whereas even images will have two pixels next to the actual center. We want to compute an equispaced array of numang angles that go between the centers of the outmost far pixels of the detector.
Parameters: Returns: angles, latpos – Angles and pixel coordinates at the detector.
Return type: two 1D ndarrays
Notes
Actually one would not need to define spacing and distance, but for convenience, these parameters are separated and an arbitrary uint ‘pixel’ is defined.
-
radontea.fan.
radon_fan
(arr, det_size, det_spacing=1, shift_size=1, lS=1, lD=None, return_ang=False, count=None, max_count=None)¶ Compute the Radon transform for a fan beam geometry
In contrast to
radon_parallel()
, this function uses (1) a fan-beam geometry (the integral is taken along rays that meet at one point), and (2) translates the object laterally instead of rotating it. The result is sometimes referred to as ‘linogram’.Problem sketch:
x ^ | ----> z source object detector / . (+det_size/2, lD) (0,-lS) ./ (0,0) . \ . \ . (-det_size/2, lD)
The algorithm computes all angular projections for discrete movements of the object. The position of the object is changed such that its lower boundary starts at (det_size/2, 0) and its upper boundary ends at (-det_size/2, 0) at increments of shift_size.
Parameters: - arr (ndarray, shape (N,N)) – the input image.
- det_size (int) – The total detector size in pixels. The detector centered to the source. The axial position of the detector is the center of the pixels on the far right of the object.
- det_spacing (float) – Distance between detection points in pixels.
- shift_size (int) – The amount of pixels that the object is shifted between projections.
- lS (multiples of 0.5) – Source position relative to the center of arr. lS >= 1.
- lD (int) – Detector position relative to the center arr. Default is N/2.
- return_ang (bool) – Also return the angles corresponding to the detector pixels.
- 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.
Returns: outarr – Linogram of the input image. Where N+det_size determines the lateral position of the sample.
Return type: ndarray of floats, shape (N+det_size,det_size)
See also
scipy.ndimage.interpolation.rotate()
- The interpolator used to rotate the image.
-
radontea.fan.
lino2sino
(linogram, lds, stepsize=1, det_spacing=1, numang=None, retang=False, count=None, max_count=None)¶ Convert linogram to sinogram for an equispaced detector.
Parameters: - linogram (real 2d ndarray of shape (D, A*)) – Linogram from synthetic aperture measurements.
- lds (float) – Distance from point source to detector in au.
- stepsize (float) – Translational increment of object in au (stepsize in D).
- det_spacing (float) – Distance between detector positions in au.
- numang (int) – Number of equispaced angles, defaults to linogram.shape[1]
- retang (bool) – Return the corresponding angles for the sinogram.
- 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.
Returns: sinogram – The distortion-corrected sinogram. If retang is True, then the equispaced angles are returned as well.
Return type: 2d ndarray of shape (D, A)
Notes
This function can be used to convert a linogram obtained with fan-beam tomography to a sinogram, which then can be reconstructed with the backprojection or fourier mapping algorithms.
Reconstruction¶
-
radontea.fan.
fan_rec
(linogram, lds, method, stepsize=1, det_spacing=1, numang=None, count=None, max_count=None, **kwargs)¶ 2D synthetic aperture reconstruction
Computes the inverse of the fan-beam Radon transform using interpolation of the linogram and one of the inverse algorithms for tomography with the Fourier slice theorem.
Parameters: - linogram (2d ndarray of shape (D, A)) – Input linogram from the synthetic aperture measurement.
- lds (float) – Distance in pixels between source and detector.
- method (callable) – Reconstruction method, e.g. radontea.backproject.
- numang (int) – Number of angles to be used for the sinogram. A higher number increases quality, but interpolation takes longer. By default numang = linogram.shape[1].
- 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.
- **kwargs (dict) – Keyword arguments for method.
Volumetric reconstruction¶
For a slice-wise 3D reconstruction, radontea can use multiprocessing to parallelize the reconstruction process.
Convenience functions¶
-
radontea.
backproject_3d
(sinogram, angles, filtering='ramp', weight_angles=True, padding=True, padval=0, count=None, max_count=None, ncpus=None)¶ Convenience wrapper for 3D backprojection reconstruction
See
backproject()
for parameter definitions. The additional parameter ncpus sets the number of CPUs used.Returns: out – The reconstructed volume. Changed in version 0.4.0: Output indexing now follows the ODTbrain convention. For the the old behavior, use
out.transpose(1, 0, 2)
.Return type: ndarray, shape (N,M,N)
-
radontea.
fourier_map_3d
(sinogram, angles, intp_method='cubic', count=None, max_count=None, ncpus=None)¶ Convenience wrapper for 3D Fourier mapping reconstruction
See
fourier_map()
for parameter definitions. The additional parameter ncpus sets the number of CPUs used.Returns: out – The reconstructed volume. Changed in version 0.4.0: Output indexing now follows the ODTbrain convention. For the the old behavior, use
out.transpose(1, 0, 2)
.Return type: ndarray, shape (N,M,N)
General 3D reconstruction¶
-
radontea.
volume_recon
(func2d, sinogram=None, angles=None, count=None, max_count=None, ncpus=None, **kwargs)¶ Slice-wise 3D inversion of the Radon transform
Computes the slice-wise 3D inverse of the Radon transform using multiprocessing.
Parameters: - func2d (callable) – A method for the slice-wise reconstruction
(e.g.
backproject()
). - sinogram (ndarray, shape (A,M,N)) – Three-dimensional sinogram of line recordings. The axis 1 iterates through the M slices. The rotation takes place through axis 1.
- angles ((A,) ndarray) – Angular positions of the sinogram in radians in the interval [0, PI).
- max_count (count,) – Can be used to monitor the progress of this function. The value of max_count.value is set initially and the value of count.value is incremented until it reaches the end of the algorithm (max_count.value).
- **kwargs (dict) – Additional keyword arguments to func2d.
Returns: out – The reconstructed volume.
Changed in version 0.4.0: Output indexing now follows the ODTbrain convention. For the the old behavior, use
out.transpose(1, 0, 2)
.Return type: ndarray, shape (N,M,N)
- func2d (callable) – A method for the slice-wise reconstruction
(e.g.