toupy.tomo package

Submodules

toupy.tomo.cone_projector module

Cone-beam forward projector.

Computes the cone-beam Radon transform (line integrals along divergent rays) of a 3-D volume for a circular source orbit.

Primary use-cases:

  • Consistency checks — compare measured projections with the re-projection of the FDK reconstruction to quantify residual reconstruction artefacts.

  • Iterative reconstruction — the forward model is the mandatory companion to the FDK back-projector in a SART-style cone-beam update step.

The function signature mirrors projector() for the parallel-beam case.

GPU pathstub analogous to radon_cuda();

will delegate to a CuPy CUDA ray-casting kernel.

CPU path : NumPy reference using trilinear interpolation.

toupy.tomo.cone_projector.cone_project(volume: ndarray, geometry: ConeBeamGeometry, cuda: bool = False) → ndarray[source]

Compute the cone-beam forward projection of a 3-D volume.

For each projection angle θ in geometry.angles and each detector pixel (u_i, v_j) the function evaluates the line integral along the cone ray from the point source through the voxel grid to the detector.

The perspective-projection mapping of a voxel at object-space coordinates (x, y, z) onto the detector at angle θ is:

\[U(\theta) = \mathrm{SOD} + x\sin\theta - y\cos\theta\]
\[u_d = \frac{\mathrm{SDD}}{U}\,(x\cos\theta + y\sin\theta), \qquad v_d = \frac{\mathrm{SDD}}{U}\, z\]

The voxel grid is assumed isotropic with pitch geometry.effective_pixel_size, centred at the origin (the rotation axis).

Parameters:

  • volume (ndarray, shape (n_v, N, N)) – The 3-D volume to project. Axis order: (z, y, x) where z is the rotation (axial) direction and (y, x) is the transaxial plane. N is the transaxial voxel count and may differ from geometry.n_u; the coordinate mapping is scaled accordingly.

  • geometry (ConeBeamGeometry) – Validated acquisition geometry. Provides SOD, SDD, det_pixel_size, and the set of projection angles.

  • cuda (bool, optional) – If True, attempt to use the CuPy GPU path. Falls back to CPU with a warnings warning when CuPy is unavailable. Default False.

Returns:

projections – Simulated cone-beam projections. Axis order matches the input convention expected by fdk_weight().

Return type:

ndarray, shape (n_angles, n_v, n_u)

Raises:

  • ValueError – If volume.ndim != 3.

  • ValueError – If volume.shape[0] != geometry.n_v (axial size mismatch between the volume and the detector).

See also

toupy.tomo.radon.projector

Parallel-beam forward projector (2-D).

toupy.tomo.fdk.fdk_reconstruct

Paired cone-beam reconstruction.

Notes

This forward projector is the exact adjoint of fdk_backproject() (without the (SOD/U)² weighting applied in the back-projector). For a consistent SART-style iterative solver, the weighting must be applied symmetrically in both the forward and back steps, or the normal equations must be formed explicitly.

toupy.tomo.fdk module

Feldkamp–Davis–Kress (FDK) cone-beam reconstruction.

Implements the three-step FDK pipeline from:

Feldkamp, L. A., Davis, L. C., & Kress, J. W. (1984). Practical cone-beam algorithm. Journal of the Optical Society of America A, 1(6), 612–619. https://doi.org/10.1364/JOSAA.1.000612

Step 1 — fdk_weight():

Multiply each projection pixel (u, v) by the cosine factor SOD / sqrt(SOD² + u² + v²), where u and v are centred physical detector coordinates.

Step 2 — fdk_filter():

Apply a 1-D ramp (or windowed-ramp) filter row-by-row along the horizontal detector direction u. Reuses compute_filter().

Step 3 — fdk_backproject():

3-D divergent-ray back-projection. For each output voxel (x, y, z) and each projection angle θ, compute the detector coordinates at which the cone ray through the voxel intersects the detector, interpolate the filtered projection, and accumulate weighted by (SOD / U)² where U is the perspective-projection denominator (see function docstring).

GPU pathstub for a CuPy/CUDA implementation, falls back to CPU with a

warning (same pattern as mod_iradon_cuda()).

CPU path : NumPy reference implementation (loops over angles).

Public API

fdk_weight, fdk_filter, fdk_backproject, fdk_reconstruct

toupy.tomo.fdk.fdk_backproject(projections_filtered: ndarray, geometry: ConeBeamGeometry, output_size: int | None = None) → ndarray[source]

3-D cone-beam back-projection (Step 3 of FDK).

For each output voxel at object-space coordinates (x, y, z) and for each projection angle θ, the perspective-projection denominator is:

\[U(\theta) = \mathrm{SOD} + x\sin\theta - y\cos\theta\]

The detector coordinates of the cone ray through the voxel are:

\[u_d = \frac{\mathrm{SDD}}{U}\, (x\cos\theta + y\sin\theta)\]
\[v_d = \frac{\mathrm{SDD}}{U}\, z\]

The filtered projection value at (u_d, v_d) is obtained by bilinear interpolation and accumulated into the voxel, weighted by the solid-angle correction factor (SOD / U)². The final volume is scaled by π / n_angles:

\[f(x, y, z) = \frac{\pi}{N_\theta} \sum_{k=1}^{N_\theta} \left(\frac{\mathrm{SOD}}{U_k}\right)^2 p_f^{\theta_k}(u_d, v_d)\]
Parameters:

  • projections_filtered (ndarray, shape (n_angles, n_v, n_u)) – Output of fdk_filter().

  • geometry (ConeBeamGeometry) – Acquisition geometry. geometry.angles supplies θ values in degrees; geometry.u_coords() and geometry.v_coords() supply the centred detector coordinates for interpolation.

  • output_size (int, optional) – Side length N (in voxels) of the square transaxial output grid. Defaults to geometry.n_u (one detector column per voxel, at effective_pixel_size sampling).

Returns:

volume – Reconstructed 3-D volume. Axis order: (z, y, x) where z is the rotation axis. The rotation axis voxel is at index (n_v//2, N//2, N//2).

Return type:

ndarray, shape (n_v, N, N)

Notes

The (SOD/U)² weighting is the cone-beam analogue of the cos²(γ) factor in 2-D fan-beam backprojection. It corrects for the solid-angle element of each voxel as seen from the source.

The CPU reference uses bilinear interpolation via scipy.interpolate.RegularGridInterpolator. A future CUDA kernel (stub in fdk_reconstruct()) will use texture-memory interpolation for speed.

toupy.tomo.fdk.fdk_filter(projections_weighted: ndarray, filter_type: str = 'ram-lak', freqcutoff: float = 1.0) → ndarray[source]

Apply a 1-D ramp filter row-by-row to cosine-weighted projections.

This is Step 2 of the FDK algorithm. Each detector row (fixed v, varying u) is independently filtered along the u direction using an FFT-based convolution with the chosen ramp kernel. The filter is computed by compute_filter(), so all filter types supported by the parallel-beam FBP are available here.

Parameters:

  • projections_weighted (ndarray, shape (n_angles, n_v, n_u)) – Output of fdk_weight().

  • filter_type (str, optional) – Frequency-domain filter name. One of 'ram-lak' (default), 'shepp-logan', 'cosine', 'hamming', 'hann', or None (no filtering, useful for debugging).

  • freqcutoff (float, optional) – Normalised frequency cutoff in (0, 1]. Default 1.0 (no cutoff beyond Nyquist).

Returns:

projections_filtered – Ramp-filtered projections, same shape and dtype as input.

Return type:

ndarray, shape (n_angles, n_v, n_u)

Notes

The filter is applied in the Fourier domain row-by-row:

\[p_f(u) = \mathcal{F}^{-1}\!\left[ \mathcal{F}[p_w](\xi)\, H(\xi) \right]\]

where \(H(\xi)\) is the chosen window function. This step is mathematically identical to the FBP filtering step; only the geometry of the subsequent back-projection differs.

toupy.tomo.fdk.fdk_reconstruct(projections: ndarray, geometry: ConeBeamGeometry, filter_type: str = 'ram-lak', freqcutoff: float = 1.0, output_size: int | None = None, cuda: bool = False) → ndarray[source]

Full FDK cone-beam reconstruction pipeline.

Chains fdk_weight() → fdk_filter() → fdk_backproject() in sequence, mirroring the interface of mod_iradon() for the parallel-beam case.

An optional GPU path (cuda=True) is stubbed out; when CuPy is unavailable it falls back to the CPU path with a warnings message, exactly as mod_iradon_cuda() does.

Parameters:

  • projections (ndarray, shape (n_angles, n_v, n_u)) – Raw cone-beam projections (Beer-Lambert log-normalised line integrals).

  • geometry (ConeBeamGeometry) – Acquisition geometry. validate() is called internally; no need to call it again before passing in.

  • filter_type (str, optional) – Ramp-filter window forwarded to fdk_filter(). Default 'ram-lak'.

  • freqcutoff (float, optional) – Normalised frequency cutoff forwarded to fdk_filter(). Default 1.0.

  • output_size (int, optional) – Transaxial reconstruction grid side length N, forwarded to fdk_backproject(). Defaults to geometry.n_u.

  • cuda (bool, optional) – If True, attempt to use the CuPy GPU back-projector. Falls back to CPU with a warnings warning when CuPy is unavailable. Default False.

Returns:

volume – Reconstructed 3-D volume.

Return type:

ndarray, shape (n_v, N, N)

Raises:

ValueError – If projections.ndim != 3 or its trailing dimensions are inconsistent with geometry.n_v and geometry.n_u.

See also

toupy.tomo.iradon.mod_iradon

Parallel-beam FBP (2-D, single slice).

toupy.tomo.iradon.mod_iradon_cuda

GPU parallel-beam FBP.

toupy.tomo.tomorecons.fdk_tomo_recons

High-level wrapper with interactive display and optional saving.

Notes

The π / N_θ scale factor is absorbed into fdk_backproject(). It differs from the π / (2 N_θ) factor in parallel-beam FBP because FDK sums over a full 360° orbit while the parallel-beam formula assumes 180°.

toupy.tomo.fdk.fdk_weight(projections: ndarray, geometry: ConeBeamGeometry) → ndarray[source]

Apply FDK cosine weights to a stack of cone-beam projections.

This is Step 1 of the FDK algorithm. Each pixel at detector position (u, v) is multiplied by the angle-independent cosine factor:

\[w(u, v) = \frac{\mathrm{SOD}}{\sqrt{\mathrm{SOD}^2 + u^2 + v^2}}\]

The weight map is computed once and broadcast over the angle axis because the geometry is fixed across all views.

Parameters:

  • projections (ndarray, shape (n_angles, n_v, n_u)) – Raw cone-beam projections. Axis order: angle first, detector row (v) second, detector column (u) last. Values are typically Beer-Lambert log-normalised line integrals.

  • geometry (ConeBeamGeometry) – Fully validated acquisition geometry. geometry.u_coords() and geometry.v_coords() supply the centred physical detector coordinates used to evaluate w(u, v).

Returns:

projections_weighted – Cosine-weighted projections, same shape and dtype as projections.

Return type:

ndarray, shape (n_angles, n_v, n_u)

Raises:

ValueError – If projections.shape[1:] does not match (geometry.n_v, geometry.n_u).

Notes

The weight approaches 1 on the central beam axis (u = v = 0) and decreases toward the detector periphery, compensating for the longer effective path length of oblique cone rays. This is the exact extension of the cos(γ) factor in the fan-beam Parker weighting to the 2-D detector case.

toupy.tomo.geometry module

Cone-beam acquisition geometry descriptor.

ConeBeamGeometry is a plain dataclass that encodes every metric parameter of a circular-orbit cone-beam CT setup. It is intentionally kept free of reconstruction imports so that the registration sub-package can import it without pulling in NumPy FFT or CuPy.

Typical usage:

from toupy.tomo.geometry import ConeBeamGeometry
import numpy as np

geom = ConeBeamGeometry(
    SOD=500.0, SDD=1000.0,
    det_pixel_size=0.1,
    n_u=1024, n_v=512,
    angles=np.linspace(0, 360, 720, endpoint=False),
)
geom.validate()
u = geom.u_coords()   # shape (1024,)
v = geom.v_coords()   # shape (512,)
class toupy.tomo.geometry.ConeBeamGeometry(SOD: float, SDD: float, det_pixel_size: float, n_u: int, n_v: int, angles: ndarray = <factory>)[source]

Bases: object

Circular-orbit cone-beam acquisition geometry.

All distances are expected in consistent physical units (e.g. mm or µm). Angles are in degrees.

Parameters:

  • SOD (float) – Source-to-object distance (source to rotation axis), in physical length units. Must satisfy SOD < SDD.

  • SDD (float) – Source-to-detector distance, in the same units as SOD.

  • det_pixel_size (float) – Physical pitch of one detector pixel (square pixels assumed), in the same units as SOD/SDD. Must be strictly positive.

  • n_u (int) – Number of detector columns (horizontal / transaxial direction). Must be a positive integer.

  • n_v (int) – Number of detector rows (vertical / axial direction). Must be a positive integer.

  • angles (ndarray, shape (n_angles,)) – Projection angles in degrees. Arbitrary spacing is allowed; the FDK weighting step uses them only to iterate over views.

magnification

Geometric magnification SDD / SOD. A voxel at the rotation axis appears enlarged by this factor on the detector.

Type:

float

effective_pixel_size

Physical size of one object-space pixel at the rotation axis, det_pixel_size / magnification. Use this value to set the reconstruction voxel size.

Type:

float

Examples

>>> import numpy as np
>>> from toupy.tomo.geometry import ConeBeamGeometry
>>> geom = ConeBeamGeometry(
...     SOD=500.0, SDD=1000.0, det_pixel_size=0.055,
...     n_u=2048, n_v=2048,
...     angles=np.linspace(0, 360, 1800, endpoint=False),
... )
>>> geom.validate()
>>> print(geom.magnification)          # 2.0
>>> print(geom.effective_pixel_size)   # 0.0275
SDD: float
SOD: float
angles: ndarray
det_pixel_size: float
property effective_pixel_size: float

Object-space pixel size at the rotation axis.

Defined as det_pixel_size / magnification, i.e. the physical extent of one detector pixel back-projected to the object plane. Use this as the reconstruction voxel pitch.

Returns:

Effective voxel pitch in the same units as det_pixel_size.

Return type:

float

property magnification: float

Geometric magnification factor SDD / SOD.

Returns:

The ratio of source-to-detector distance to source-to-object distance. Always > 1 for a valid geometry (SDD > SOD).

Return type:

float

Raises:

ZeroDivisionError – If SOD is zero (caught earlier by validate()).

n_u: int
n_v: int
u_coords() → ndarray[source]

Centred horizontal (transaxial) detector coordinate array.

Returns a 1-D array of length n_u containing the physical detector-plane u-coordinates of each column centre, measured from the beam axis. The coordinate of column i is:

u_i = (i - (n_u - 1) / 2) * det_pixel_size
Returns:

u – Detector u-coordinates in physical units, centred on zero.

Return type:

ndarray, shape (n_u,)

v_coords() → ndarray[source]

Centred vertical (axial) detector coordinate array.

Analogous to u_coords() but along the vertical detector axis. The coordinate of row j is:

v_j = (j - (n_v - 1) / 2) * det_pixel_size
Returns:

v – Detector v-coordinates in physical units, centred on zero.

Return type:

ndarray, shape (n_v,)

validate() → None[source]

Check self-consistency of all geometry fields.

Should be called once after construction, before passing the geometry object to any reconstruction or registration function.

Raises:

ValueError – With a descriptive message for each invalid condition: * SOD <= 0 — source-to-object distance must be positive. * SDD <= 0 — source-to-detector distance must be positive. * SOD >= SDD — the detector must lie beyond the object; equivalently, magnification > 1 must hold. * det_pixel_size <= 0 — pixel pitch must be positive. * n_u <= 0 or n_v <= 0 — detector dimensions must be positive integers. * len(angles) == 0 — at least one projection angle is required.

toupy.tomo.iradon module

Filtered Back-Projection and SART reconstruction.

GPU path : CuPy + a custom CUDA kernel CPU path : pure NumPy/SciPy, identical to the original mod_iradon.

The public API is unchanged:

compute_angle_weights, compute_filter, mod_iradon, mod_iradon_cuda, backprojector, reconsSART

toupy.tomo.iradon.backprojector(sinogram, theta, **params)[source]

Wrapper that dispatches to the GPU or CPU FBP implementation.

Parameters:

  • sinogram (ndarray shape (nbins, nangles))

  • theta (ndarray)

  • params (dict) – params["opencl"] — repurposed as params["cuda"]: set True to use the CUDA path. The key opencl is still accepted as an alias for backwards compatibility. All other keys are forwarded to the chosen implementation.

Returns:

recons

Return type:

ndarray shape (output_size, output_size)

toupy.tomo.iradon.compute_angle_weights(theta)[source]

Compute per-angle weights for non-equally-spaced angular distributions.

Parameters:

theta (ndarray) – Angles in degrees.

Returns:

weights – Weight for each angle, proportional to the angular distance to its neighbours. Forked from odtbrain.util.compute_angle_weights_1d.

Return type:

ndarray

toupy.tomo.iradon.compute_filter(nbins, filter_type='ram-lak', derivatives=False, freqcutoff=1)[source]

Compute the frequency-domain filter for FBP reconstruction.

Parameters:

  • nbins (int) – Number of detector bins (sinogram rows).

  • filter_type (str, optional) – One of ram-lak, shepp-logan, cosine, hamming, hann. Default is ram-lak.

  • derivatives (bool, optional) – Use a Hilbert filter suited for derivative projections. Default False.

  • freqcutoff (float, optional) – Normalised frequency cutoff in [0, 1]. Default 1 (no cutoff).

Returns:

fourier_filter

Return type:

ndarray shape (projection_size_padded, 1) complex64

toupy.tomo.iradon.mod_iradon(radon_image, theta=None, output_size=None, filter_type='ram-lak', derivatives=False, interpolation='linear', circle=False, freqcutoff=1)[source]

Inverse Radon transform — CPU path (pure NumPy/SciPy).

Reconstruct an image from its Radon transform using filtered back-projection.

Parameters:

  • radon_image (ndarray shape (nbins, nangles)) – Sinogram. Each column corresponds to one projection angle.

  • theta (ndarray, optional) – Projection angles in degrees. Defaults to nangles evenly-spaced values in [0, 180).

  • output_size (int, optional) – Side length of the square output image.

  • filter_type (str, optional) – See compute_filter().

  • derivatives (bool, optional) – Set True if the sinogram contains projection derivatives.

  • interpolation (str, optional) – linear (default), nearest, or cubic.

  • circle (bool, optional) – Zero pixels outside the inscribed circle.

  • freqcutoff (float, optional) – Normalised frequency cutoff.

Returns:

reconstructed

Return type:

ndarray shape (output_size, output_size)

toupy.tomo.iradon.mod_iradon_cuda(radon_image, theta=None, output_size=None, filter_type='ram-lak', derivatives=False, circle=False, freqcutoff=1)[source]

Inverse Radon transform — CUDA/GPU path.

Falls back to mod_iradon() automatically if CuPy is unavailable.

Parameters:

  • radon_image (ndarray shape (nbins, nangles)) – Sinogram on CPU. Converted to float32 internally.

  • theta (ndarray, optional) – Projection angles in degrees.

  • output_size (int, optional) – Side length of the square output image.

  • filter_type (str, optional) – See compute_filter().

  • derivatives (bool, optional) – Set True if the sinogram contains projection derivatives.

  • circle (bool, optional) – Zero pixels outside the inscribed circle.

  • freqcutoff (float, optional) – Normalised frequency cutoff.

Returns:

reconstructed

Return type:

ndarray shape (output_size, output_size) float32

toupy.tomo.iradon.reconsSART(sinogram, theta, num_iter=2, FBPinitial_guess=True, relaxation_params=0.15, **params)[source]

Tomographic reconstruction with the SART algorithm.

Parameters:

  • sinogram (ndarray shape (nbins, nangles))

  • theta (ndarray)

  • num_iter (int, optional) – Number of SART iterations. Default 2.

  • FBPinitial_guess (bool, optional) – Seed SART with an FBP reconstruction. Default True.

  • relaxation_params (float, optional) – SART relaxation parameter. Default 0.15.

Returns:

recons_sart

Return type:

ndarray

toupy.tomo.radon module

Forward Radon transform (projection).

GPU path : CuPy + the CUDA kernel defined in iradon.py CPU path : skimage.transform.radon (unchanged fallback).

Public API is unchanged: projector()

toupy.tomo.radon.projector(recons, theta, **params)[source]

Wrapper to choose between CUDA and CPU forward Radon transform.

Parameters:

  • recons (ndarray shape (N, N)) – Tomographic slice to project.

  • theta (ndarray) – Projection angles in degrees.

  • params (dict) – params["cuda"] — True → use GPU path. params["opencl"] — accepted as alias for backwards compatibility.

Returns:

sinogram

Return type:

ndarray shape (N, nangles)

toupy.tomo.radon.radon_cuda(recons, theta)[source]

Forward Radon transform — CUDA/GPU path.

Falls back to skimage.transform.radon automatically if CuPy is unavailable.

Parameters:

  • recons (ndarray shape (N, N)) – The tomographic slice to project. Converted to float32 internally.

  • theta (ndarray) – Projection angles in degrees.

Returns:

sinogram – The computed sinogram, in the same layout as skimage’s radon output (each column = one projection).

Return type:

ndarray shape (N, nangles)

toupy.tomo.tomorecons module

toupy.tomo.tomorecons.fdk_tomo_recons(projections, geometry: ConeBeamGeometry, **params)[source]

High-level FDK cone-beam reconstruction wrapper.

Mirrors tomo_recons() for cone-beam geometry. Validates the geometry, applies optional pre-processing, calls the three-step FDK pipeline, and optionally displays the central reconstructed slice.

Parameters:

  • projections (ndarray, shape (n_angles, n_v, n_u)) – Flat-field-corrected, Beer-Lambert log-normalised cone-beam projections.

  • geometry (ConeBeamGeometry) – Validated acquisition geometry.

  • **params –

    Algorithm parameters. Recognised keys:

    filtertypestr

    Ramp-filter window (forwarded to fdk_filter()). Default 'ram-lak'.

    freqcutofffloat

    Normalised frequency cutoff in (0, 1]. Default 1.

    output_sizeint or None

    Transaxial reconstruction grid side length. None uses geometry.n_u.

    cudabool

    Use CuPy GPU path. Default False.

    circlebool

    Multiply the reconstructed volume by a cylindrical mask to suppress border artefacts. Default False.

    showreconsbool

    Display the central axial slice after reconstruction. Default False.

    colormapstr

    Colormap for display. Default 'bone'.

    vmin_plotfloat or None

    Display window minimum.

    vmax_plotfloat or None

    Display window maximum.

    autosavebool

    Save the volume to disk without prompting. Default False.

Returns:

volume – Reconstructed 3-D volume.

Return type:

ndarray, shape (n_v, N, N)

Raises:

ValueError – If projections.ndim != 3 or shape is inconsistent with geometry.

See also

tomo_recons

Parallel-beam (2-D) slice reconstruction wrapper.

full_tomo_recons

Slice-by-slice parallel-beam full reconstruction.

toupy.tomo.fdk.fdk_reconstruct

Bare FDK pipeline (no display/save).

toupy.tomo.tomorecons.full_tomo_recons(input_stack, theta, **params)[source]

Full tomographic reconstruction

Parameters:

  • input_stack (ndarray) – A 3-dimensional array containing the stack of projections. The order should be [projection_num, row, column]

  • theta (ndarray) – A 1-dimensional array of thetas

  • params (dict) – Dictionary containing additional parameters

  • params["algorithm"] (str) – Choice of algorithm. Two algorithm implemented: “FBP” and “SART”

  • params["slicenum"] (int) – Slice number

  • params["filtertype"] (str) – Filter to use for FBP

  • params["freqcutoff"] (float) – Frequency cutoff (between 0 and 1)

  • params["circle"] (bool) – Multiply the reconstructed slice by a circle to remove borders

  • params["derivatives"] (bool) – If the projections are derivatives. Only for FBP.

  • params["calc_derivatives"] (bool) – Calculate derivatives of the sinogram if not done yet.

  • params["opencl"] (bool) – Implement the tomographic reconstruction in CUDA

  • params["autosave"] (bool) – Save the data at the end without asking

  • params["vmin_plot"] (float) – Minimum value for the gray level at each display

  • params["vmax_plot"] (float) – Maximum value for the gray level at each display

  • params["colormap"] (str) – Colormap

  • params["showrecons"] (bool) – If to show the reconstructed slices

Returns:

tomogram – A 3-dimensional array containing the full reconstructed tomogram.

Return type:

ndarray

toupy.tomo.tomorecons.tomo_recons(sinogram, theta, **params)[source]

Wrapper to select tomographic algorithm

sinogramndarray

A 2-dimensional array containing the sinogram

thetandarray

A 1-dimensional array of thetas

paramsdict

Dictionary containing additional parameters

params[“algorithm”]str

Choice of algorithm. Two algorithm implemented: “FBP” and “SART”

params[“slicenum”]int

Slice number

params[“filtertype”]str

Name of the filter to be applied in frequency domain filtering. The options are: ram-lak, shepp-logan, cosine, hamming, hann. Assign None to use no filter.

params[“freqcutoff”]float

Frequency cutoff (between 0 and 1)

params[“circle”]bool

Multiply the reconstructed slice by a circle to remove borders

params[“weight_angles”]bool

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

\[\Delta \phi_0 \longmapsto \Delta \phi_j =\]

rac{phi_{j+1} - phi_{j-1}}{2}

params[“derivatives”]bool

If the projections are derivatives. Only for FBP.

params[“calc_derivatives”]bool

Calculate derivatives of the sinogram if not done yet.

params[“opencl”]bool

Implement the tomographic reconstruction in opencl as implemented in Silx

params[“autosave”]bool

Save the data at the end without asking

params[“vmin_plot”]float

Minimum value for the gray level at each display

params[“vmax_plot”]float

Maximum value for the gray level at each display

params[“colormap”]str

Colormap

params[“showrecons”]bool

If to show the reconstructed slices

reconsndarray

A 2-dimensional array containing the reconstructed slice.

toupy.tomo.tv_recons module

Total-Variation (TV) regularised tomographic reconstruction.

Implements the Chambolle-Pock primal-dual algorithm (CP) to solve:

min_{x >= 0} (1/2) ||A x - b||^2 + lambda * TV(x)

where A is the Radon forward operator, b is the (possibly noisy) sinogram, and TV is the isotropic discrete total-variation.

The algorithm is an exact first-order method and handles the non-smooth TV term via proximal splitting. It is well-suited to ptychographic tomography, where the phase sinogram is noisy or has limited angular sampling.

Reference:

Chambolle, A., & Pock, T. (2011). A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision, 40(1), 120–145. https://doi.org/10.1007/s10107-010-0436-y

Sidky, E. Y., & Pan, X. (2008). Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization. Physics in Medicine & Biology, 53(17), 4777. https://doi.org/10.1088/0031-9155/53/17/014

GPU notes

Both public functions accept cuda=True. In GPU mode all primal and dual variables reside on the device as CuPy float32 arrays throughout the iteration; only the final reconstruction is transferred back to the host. The Radon forward and backward operators reuse the CUDA kernels compiled in toupy.tomo.iradon; no additional compilation is required.

The step-size estimate is also computed entirely on the GPU via a five-step power iteration, so the first call is slightly slower than subsequent calls (kernel reuse).

toupy.tomo.tv_recons.chambolle_pock_tv(sinogram, theta, lam=0.01, n_iter=100, output_size=None, cuda=False, verbose=False)[source]

Chambolle-Pock TV-regularised reconstruction.

Solves the primal-dual saddle-point problem:

\[\min_{x \ge 0}\; \max_{p,q}\; \langle K x,\, (p, q) \rangle - G^*(p, q) + F(x)\]

where \(K = (A,\, \nabla)\) stacks the Radon operator and the gradient, \(G^*\) is the convex conjugate of the data and TV terms, and \(F\) is the non-negativity indicator.

Parameters:

  • sinogram (ndarray, shape (n_pixels, n_angles)) – Measured sinogram with axes (detector pixel, angle).

  • theta (array_like, shape (n_angles,)) – Projection angles [degrees].

  • lam (float, optional) – TV regularisation weight. Default 1e-2.

  • n_iter (int, optional) – Number of Chambolle-Pock iterations. Default 100.

  • output_size (int or None, optional) – Side length of the square reconstruction grid. When None the value is inferred from the sinogram as round(n_pixels / sqrt(2)), matching skimage.transform.radon with circle=False.

  • cuda (bool, optional) – Run on GPU via CuPy. All primal/dual variables reside on the device throughout the iteration; only the final image is transferred back. Falls back to CPU with a warning when CuPy is unavailable. Default False.

  • verbose (bool, optional) – Print cost every 10 iterations. Default False.

Returns:

x – Reconstructed image. Always returned as a CPU numpy array.

Return type:

ndarray, shape (output_size, output_size)

See also

tv_reconstruction

High-level wrapper with FBP warm-start.

toupy.tomo.tv_recons.tv_reconstruction(sinogram, theta, lam=0.01, n_iter=100, output_size=None, init='fbp', cuda=False, verbose=False)[source]

TV-regularised tomographic reconstruction.

High-level wrapper around chambolle_pock_tv() with an optional FBP warm-start and the same sinogram-axis convention used throughout toupy (detector pixels along axis 0).

Parameters:

  • sinogram (ndarray, shape (n_pixels, n_angles)) – Sinogram — same axis convention as mod_iradon().

  • theta (array_like, shape (n_angles,)) – Projection angles in degrees.

  • lam (float, optional) – TV regularisation weight. Default 1e-2.

  • n_iter (int, optional) – Number of Chambolle-Pock iterations. Default 100.

  • output_size (int or None, optional) – Side length of the square reconstruction grid. None infers from the sinogram (matches skimage.transform.radon).

  • init ({'fbp', 'zeros'}, optional) – Initialisation strategy. 'fbp' uses a Ram-Lak FBP as the starting point (recommended; faster convergence). 'zeros' starts from the zero image. Default 'fbp'.

  • cuda (bool, optional) – Run on GPU via CuPy. Default False.

  • verbose (bool, optional) – Print cost every 10 iterations. Default False.

Returns:

recons – Reconstructed image.

Return type:

ndarray, shape (n, n)

Notes

The TV weight lam should be tuned to the noise level. A cross-validation or L-curve approach is recommended for quantitative work.

Examples

>>> import numpy as np
>>> from skimage.transform import radon
>>> from toupy.simulation import phantom
>>> from toupy.tomo import tv_reconstruction
>>> img = phantom(256)
>>> theta = np.linspace(0, 180, 180, endpoint=False)
>>> sino = radon(img, theta=theta, circle=False)
>>> recons = tv_reconstruction(sino, theta, lam=0.01, n_iter=50)