toupy.registration package

Submodules

toupy.registration.cone_registration module

Cone-beam projection alignment and geometric calibration.

Adapts the parallel-beam registration tools from toupy.registration.registration to the cone-beam case. The key differences are:

  1. Centre-of-rotation (CoR) scaling — a detector-space horizontal shift of Δu pixels corresponds to an object-space shift of Δu * SOD / SDD voxels (inverse magnification factor).

  2. Fan-angle correction for vertical alignment — the apparent vertical position of the rotation axis varies with horizontal detector position due to the divergent beam. Tilt estimation must account for this fan-angle-dependent gradient.

  3. Geometric calibration — the sinusoidal trajectory of a ball-bearing marker on the detector encodes all misalignment parameters (SOD, SDD, detector tilt, CoR offset). ConeBeamRegistration.estimate_geometry() fits this trajectory.

  4. Fan-to-parallel rebinning — the central detector row can be exactly rebinned from fan-beam to parallel-beam coordinates (rebin_fan_to_parallel()), after which the full parallel-beam registration toolbox in toupy.registration.registration applies without modification.

All methods follow the NumPy docstring convention used throughout toupy.

class toupy.registration.cone_registration.ConeBeamRegistration(geometry: ConeBeamGeometry)[source]

Bases: object

Projection alignment and geometric calibration for cone-beam CT.

Wraps a ConeBeamGeometry and exposes methods that mirror the parallel-beam functions alignprojections_vertical(), alignprojections_horizontal(), and estimate_rot_axis(), with cone-beam magnification scaling applied where appropriate.

Parameters:

geometry (ConeBeamGeometry) – Validated acquisition geometry. The SOD, SDD, and det_pixel_size fields are used in every shift-scaling calculation.

geometry

The geometry passed at construction (stored as-is, not copied).

Type:

ConeBeamGeometry

Examples

>>> import numpy as np
>>> from toupy.tomo.geometry import ConeBeamGeometry
>>> from toupy.registration.cone_registration import ConeBeamRegistration
>>> 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()
>>> reg = ConeBeamRegistration(geom)
>>> cor_shift = reg.estimate_cor(sinogram)
align_horizontal(projections: ndarray, shiftstack: ndarray, **params) → ndarray[source]

Align projections horizontally using cone-beam tomographic consistency.

Mirrors alignprojections_horizontal() but replaces the parallel-beam FBP inner loop with an FDK reconstruction + cone-beam re-projection loop, and scales the CoR formula by the inverse magnification SOD / SDD when converting detector shifts to object-space corrections.

The cross-correlation inner loop is unchanged because it operates entirely in detector space (detector-space shifts are magnification-agnostic).

Parameters:

  • projections (ndarray, shape (n_angles, n_v, n_u)) – Stack of flat-field-corrected cone-beam projections.

  • shiftstack (ndarray, shape (2, n_angles)) – Current shift estimates [vertical_shifts, horizontal_shifts] in detector pixels. The horizontal row (index 1) is updated and returned; the vertical row is left unchanged.

  • **params –

    Algorithm control parameters. Required keys:

    maxitint

    Maximum number of outer alignment iterations.

    pixtolfloat

    Convergence tolerance in detector pixels.

    freqcutofffloat

    Low-pass cutoff applied to the sinogram before computing the cost function.

    shiftmethstr

    Interpolation method for sub-pixel shifting ('linear', 'fourier', or 'spline').

    filtertypestr

    Ramp-filter type forwarded to the FDK step.

    circlebool

    Apply a cylindrical mask to the FDK reconstruction before re-projection.

Returns:

shiftstack – Updated shift array with optimised horizontal detector shifts in row 1 (units: detector pixels).

Return type:

ndarray, shape (2, n_angles)

Notes

The synthetic sinogram re-projection uses cone_project() instead of projector() to ensure the forward model is consistent with cone-beam divergent-ray geometry.

align_horizontal_from_parallel(projections: ndarray, shiftstack: ndarray, row: int | None = None, **params) → ndarray[source]

Horizontal alignment using fan-to-parallel rebinning + FBP consistency.

Converts the selected detector row to a parallel-beam sinogram via rebin_fan_to_parallel(), runs the full parallel-beam tomographic-consistency alignment (alignprojections_horizontal()) on the rebinned sinogram, then maps the resulting per-angle shifts back to fan-beam detector-pixel shifts.

This approach is more accurate than the direct align_horizontal() method because:

  • The FBP inner loop in alignprojections_horizontal() expects sinusoidal feature traces, which are only present in a parallel-beam sinogram. After rebinning, every object feature follows an exact sinusoid s(φ) = x₀ cos φ + y₀ sin φ.

  • Per-angle shifts (not just a global CoR estimate) are computed.

  • Frequency-cutoff schedules and spatial multiresolution warm-starts are available via params.

Unit conversion summary

Let ds be the parallel-sinogram pixel pitch (mm) and p be the fan-beam detector pixel size (mm).

Fan → parallel (warm-start):

\[\delta s_{\text{par}}[j] = \delta u_{\text{fan}}[i] \;\cdot\; \frac{\mathrm{SOD}}{\mathrm{SDD}} \;\cdot\; \frac{p}{ds}\]

Parallel → fan (final conversion):

\[\delta u_{\text{fan}} = \frac{\mathrm{SDD}}{p} \tan\!\left( \arcsin\!\left( \frac{\delta s_{\text{par}} \cdot ds}{\mathrm{SOD}} \right) \right)\]
Parameters:

  • projections (ndarray, shape (n_angles, n_v, n_u)) – Flat-field-corrected cone-beam projection stack.

  • shiftstack (ndarray, shape (2, n_angles)) – Current shift estimates [vertical_shifts, horizontal_shifts] in fan-beam detector pixels. Row 1 (horizontal) is updated in-place and returned; row 0 (vertical) is unchanged.

  • row (int or None, optional) – Detector row to use for the horizontal sinogram. Defaults to the central row n_v // 2.

  • **params –

    All keyword arguments are forwarded to alignprojections_horizontal(). Required keys:

    maxitint

    Maximum iterations per frequency stage.

    pixtolfloat

    Convergence tolerance in pixels.

    freqcutofffloat

    Low-pass sinogram cutoff (fraction of Nyquist).

    shiftmethstr

    Sub-pixel shift method ('linear', 'fourier', 'spline').

    circlebool

    Apply circular mask to the FBP reconstruction.

    derivativesbool

    Set True for phase-gradient projections.

    calc_derivativesbool

    Compute derivatives of the synthetic sinogram.

Returns:

shiftstack – Updated shift array with optimised fan-beam horizontal shifts in row 1 (units: detector pixels).

Return type:

ndarray, shape (2, n_angles)

See also

align_horizontal

Direct CoR-based horizontal alignment (no rebinning).

rebin_to_parallel

Standalone rebinning wrapper.

toupy.registration.registration.alignprojections_horizontal

Parallel-beam alignment routine used internally.

Examples

>>> shiftstack = np.zeros((2, len(geom.angles)))
>>> shiftstack = reg.align_horizontal_from_parallel(
...     projections, shiftstack,
...     maxit=20, pixtol=0.01,
...     freqcutoff=0.8, shiftmeth='fourier',
...     filtertype='ram-lak', circle=True,
...     derivatives=False, calc_derivatives=False,
... )
align_vertical(projections: ndarray, shiftstack: ndarray, **params) → ndarray[source]

Align projections vertically, accounting for cone-beam fan angle.

Mirrors alignprojections_vertical() with a cone-beam correction for the fan-angle-dependent vertical gradient. A detector tilt of angle α produces a spurious vertical shift that grows linearly with horizontal detector position:

\[\Delta v_{\text{fan}}(u) \approx u \cdot \frac{v_{\text{tilt}}}{\mathrm{SDD}}\]

This systematic term is estimated and subtracted from the vertical mass-fluctuation signal before the optimiser runs, so that only residual mechanical drift remains in the cost function.

Parameters:

  • projections (ndarray, shape (n_angles, n_v, n_u)) – Stack of flat-field-corrected cone-beam projections.

  • shiftstack (ndarray, shape (2, n_angles)) – Current shift estimates. Only the vertical row (index 0) is updated.

  • **params –

    Algorithm parameters. Required keys:

    polyorderint

    Polynomial order for baseline removal from the vertical fluctuation signal.

    limsylist of int or None

    Explicit [row_start, row_end] detector row limits. None uses the full detector height.

    deltaxint

    Horizontal margin in detector pixels to exclude from the mass-fluctuation ROI (avoids edge artefacts).

    maxitint

    Maximum number of outer iterations.

    pixtolfloat

    Convergence tolerance in detector pixels.

Returns:

shiftstack – Updated shift array with optimised vertical shifts in row 0 (units: detector pixels).

Return type:

ndarray, shape (2, n_angles)

See also

toupy.registration.registration.alignprojections_vertical

Parallel-beam analogue (no fan-angle correction).

align_vertical_cone(projections: ndarray, shiftstack: ndarray, **params) → tuple[ndarray, float][source]

Vertical alignment with per-column centroid and axis-tilt estimation.

Extends align_vertical() with two cone-beam-specific improvements:

1. Per-column centroid decomposition

Instead of collapsing the entire projection into a single vertical mass centroid, the detector is divided into n_col_bins vertical strips. For each strip at horizontal position u_k a centroid v_c(θ, u_k) is computed and the column-dependence is modelled as

\[v_c(\theta, u_k) = a(\theta) + b(\theta)\,(u_k - u_0)\]

where u_0 is the detector centre column. A least-squares fit over the column bins separates:

  • a(θ) — global vertical position of the centre of mass (mechanical drift + DC offset).

  • b(θ) — slope due to rotation-axis tilt in the detector plane. A tilt angle α gives b ≈ tan α / magnification.

2. Sinusoidal modulation removal

For a sample whose centre of mass is at (x_c, y_c, z_c) in object space, the expected centroid is

\[a_{\text{theory}}(\theta) \approx \frac{\mathrm{SDD}}{\mathrm{SOD}}\,z_c \left( 1 - \frac{x_c\sin\theta - y_c\cos\theta}{\mathrm{SOD}} \right)\]

The sinusoidal term B\sin\theta + C\cos\theta (whose amplitude grows with off-axis horizontal position) is fitted and removed from a(θ) before the drift estimate is computed. This prevents a laterally displaced sample from being misidentified as having per-angle vertical drift.

Parameters:

  • projections (ndarray, shape (n_angles, n_v, n_u)) – Flat-field-corrected cone-beam projection stack.

  • shiftstack (ndarray, shape (2, n_angles)) – Current shift estimates [vertical_shifts, horizontal_shifts] in detector pixels. Row 0 (vertical) is updated; row 1 is used to pre-shift projections but is not changed.

  • **params –

    Algorithm control parameters. Required keys:

    maxitint

    Maximum outer iterations (default 20).

    pixtolfloat

    Convergence tolerance in detector pixels (default 0.01).

    polyorderint

    Polynomial order for the baseline model of the centroid signal a(θ) (default 0).

    • 0 — baseline = mean of a(θ). The correction becomes det_v_centre − a(θ), which removes every trend (DC offset, linear ramp, sinusoidal drift). Recommended for circular-orbit scans.

    • 1 — baseline = best-fit line. Only the mean and fast jitter are corrected; slow linear trends are treated as natural and left in the data. Use this only if you deliberately want to preserve a slow linear vertical motion.

    deltaxint

    Horizontal margin in detector pixels to exclude from the column-strip ROI (default 0).

    limsylist of int or None

    Explicit [row_start, row_end] detector row limits. None uses the full detector height.

    n_col_binsint

    Number of equal-width vertical strips for the per-column centroid estimate (default 8). More bins give a more robust tilt estimate but noisier per-strip centroids. Must be ≥ 2.

    remove_sinusoidalbool

    If True, fit and remove the A + B·sinθ + C·cosθ component from a(θ) before drift estimation (default False).

    This is designed to remove the tiny 1/U(θ) modulation that arises when the sample’s horizontal centre of mass is far from the rotation axis. Do not enable this if the scan has genuine sinusoidal drift, because the fit cannot distinguish the two: both have the same frequency. Only enable if you have independently confirmed that the horizontal CoR offset is large (e.g. from estimate_cor) and that per-angle drift is negligible. Requires a full 360° scan.

    apply_tilt_correctionbool

    If True, apply the estimated tilt as a per-column vertical warp to each projection after drift correction (default False). Setting this to True modifies projections in-place because the tilt correction cannot be represented as a scalar shift per projection.

Returns:

  • shiftstack (ndarray, shape (2, n_angles)) – Updated array with optimised per-angle vertical shifts in row 0 (detector pixels).

  • tilt_deg (float) – Estimated rotation-axis tilt angle in degrees (in the detector plane, i.e. the angle between the rotation axis and the detector v-axis). Positive tilt means the top of the rotation axis leans toward positive u. Pass this value to estimate_geometry() or use it as a diagnostic.

Notes

The tilt correction (apply_tilt_correction=True) applies a shear warp v → v − b·(u − u_0) to every projection using scipy.ndimage.map_coordinates. This is an in-place modification of the projections array and cannot be undone without re-loading the raw data. Only enable it when the tilt estimate has converged and you are satisfied with the value.

See also

align_vertical

Simpler single-centroid vertical alignment.

align_horizontal_from_parallel

Horizontal alignment via rebinning.

Examples

>>> shiftstack = np.zeros((2, len(geom.angles)))
>>> shiftstack, tilt_deg = reg.align_vertical_cone(
...     projections, shiftstack,
...     maxit=20, pixtol=0.01,
...     polyorder=0, deltax=5,
...     n_col_bins=8,
...     remove_sinusoidal=False,
...     apply_tilt_correction=False,
... )
>>> print("Axis tilt: {:.3f} deg".format(tilt_deg))
estimate_cor(sinogram: ndarray) → float[source]

Estimate the centre-of-rotation offset from a horizontal sinogram.

Adapts the mass-centroid approach of center_of_mass_stack() to the cone-beam case. A detector-space shift Δu (in detector pixels) corresponds to an object-space CoR offset of:

\[\Delta x_{\mathrm{obj}} = \Delta u \cdot \frac{\mathrm{SOD}}{\mathrm{SDD}}\]

Algorithm:

  1. Compute the intensity-weighted horizontal centroid of each projection (one centroid per angle).

  2. Average the centroids over all angles.

  3. Compare the average centroid with the expected detector centre (n_u - 1) / 2 to obtain the detector-space offset Δu.

  4. Convert to object-space voxels via the inverse magnification.

Parameters:

sinogram (ndarray, shape (n_u, n_angles)) – Horizontal sinogram — one row of the projection stack transposed so that detector columns are rows and angles are columns. Matches the layout used by alignprojections_horizontal().

Returns:

cor_shift – Centre-of-rotation offset in object-space voxels. A positive value means the rotation axis is displaced to the right of the detector centre. Pass the negated value as a horizontal correction to align_horizontal().

Return type:

float

See also

toupy.registration.registration.center_of_mass_stack

Parallel-beam analogue (no magnification correction).

toupy.registration.registration.estimate_rot_axis

Interactive parallel-beam CoR estimation from sinogram symmetry.

estimate_geometry(projections: ndarray, angles: ndarray, init_geometry: ConeBeamGeometry) → ConeBeamGeometry[source]

Calibrate the cone-beam geometry from a ball-bearing sinogram.

A small dense marker (e.g. a tungsten or steel ball-bearing) placed off-axis traces a sinusoidal path on the detector as the sample rotates. The trajectory parameters encode the true SOD, SDD, detector tilt, and horizontal CoR offset.

Algorithm

  1. Marker detection — locate the marker centre in each projection via thresholding and blob centroiding.

  2. Horizontal fit — the horizontal detector trajectory follows:

    \[u_d(\theta) = \frac{\mathrm{SDD}}{\mathrm{SOD}} \left( r\cos(\theta + \phi_0) + x_0 \right)\]

    where r is the marker orbit radius, φ₀ its initial phase, and x₀ the horizontal CoR offset on the detector.

  3. Vertical fit — the vertical detector trajectory is approximately:

    \[v_d(\theta) \approx v_0 + v_{\text{tilt}}\sin(\theta + \phi_1)\]

    The sinusoidal term arises from a small tilt of the detector or rotation axis.

  4. Geometry extraction — solve for SOD and SDD from the fitted amplitude of the horizontal sinusoid using the known marker orbit radius r (measured independently or estimated from the projection extent).

param projections:

Projections containing the calibration marker. Should be flat-field corrected but not log-normalised (the marker appears as a dense spot in transmission).

type projections:

ndarray, shape (n_angles, n_v, n_u)

param angles:

Projection angles in degrees, same ordering as the first axis of projections.

type angles:

ndarray, shape (n_angles,)

param init_geometry:

Initial geometry estimate used as the starting point for the non-linear sinusoidal fit and to supply det_pixel_size, n_u, and n_v.

type init_geometry:

ConeBeamGeometry

returns:

refined_geometry – A new ConeBeamGeometry with updated SOD, SDD, and angles fields. det_pixel_size, n_u, and n_v are inherited from init_geometry.

rtype:

ConeBeamGeometry

raises ValueError:

If fewer than 3 projection angles are provided (under-determined sinusoidal fit).

raises RuntimeError:

If the non-linear sinusoidal fit fails to converge within the allowed number of iterations.

See also

toupy.registration.registration.estimate_rot_axis

Interactive parallel-beam CoR estimation (no SOD/SDD fitting).

scipy.optimize.curve_fit

Underlying non-linear least-squares solver used for the fit.

rebin_to_parallel(projections_central: ndarray, n_s: int | None = None, n_phi: int | None = None, average_redundant: bool = True) → tuple[ndarray, ndarray, ndarray][source]

Rebin the central fan-beam row to a parallel-beam sinogram.

Thin wrapper around rebin_fan_to_parallel() that reads the geometry from self.geometry. After rebinning, the returned sinogram can be passed directly to the parallel-beam registration functions in toupy.registration.registration.

Parameters:

  • projections_central (ndarray, shape (n_angles, n_u)) – Central detector row (v = 0) of the cone-beam projection stack.

  • n_s (int or None, optional) – Number of output parallel-beam detector bins. Defaults to geometry.n_u.

  • n_phi (int or None, optional) – Number of output parallel-beam angles in [0°, 180°). Defaults to len(geometry.angles).

  • average_redundant (bool, optional) – If True (default), average the two redundant views available in a full 360° scan (each ray is measured from both sides). Set to False to use only the primary view.

Returns:

  • sino_parallel (ndarray, shape (n_s, n_phi)) – Parallel-beam sinogram compatible with mod_iradon() and estimate_rot_axis().

  • s_coords (ndarray, shape (n_s,)) – Centred physical detector coordinates of the parallel-beam grid in mm.

  • phi_angles (ndarray, shape (n_phi,)) – Projection angles in degrees, uniformly distributed in [0°, 180°).

See also

rebin_fan_to_parallel

Standalone version of this function.

toupy.registration.registration.estimate_rot_axis

Interactive parallel-beam CoR estimation (works directly on the returned sinogram).

toupy.registration.registration.alignprojections_horizontal

Parallel-beam horizontal alignment (works on the returned sinogram after extracting shiftstack).

toupy.registration.cone_registration.rebin_fan_to_parallel(projections_central: ndarray, geometry: ConeBeamGeometry, n_s: int | None = None, n_phi: int | None = None, average_redundant: bool = True) → tuple[ndarray, ndarray, ndarray][source]

Exactly rebin a central-row fan-beam sinogram to parallel-beam.

For a flat-panel detector the fan-to-parallel mapping for the central detector row (v = 0) is an exact geometric identity (not an approximation). Each cone ray from the source through detector column u at source angle θ corresponds to a unique parallel-beam ray characterised by:

\[s = \frac{\mathrm{SOD}\,u}{\sqrt{\mathrm{SDD}^2 + u^2}} = \mathrm{SOD}\sin\gamma\]
\[\phi = \left(\theta + \gamma - \frac{\pi}{2}\right) \bmod \pi\]

where \(\gamma = \arctan(u / \mathrm{SDD})\) is the fan angle, \(s\) is the perpendicular distance from the rotation axis to the ray, and \(\phi\) is the parallel-beam projection angle in the convention used by mod_iradon().

The inverse mapping (used for output-driven interpolation) is:

\[\gamma = \arcsin(s / \mathrm{SOD}),\quad u = \mathrm{SDD}\tan\gamma,\quad \theta = \phi - \gamma - \frac{\pi}{2} \pmod{2\pi}\]

For a full 360° acquisition each parallel ray is measured twice (from opposite source positions, θ and θ + π). When average_redundant=True these two measurements are averaged, improving SNR and suppressing ring artefacts from detector inhomogeneities.

Parameters:

  • projections_central (ndarray, shape (n_angles, n_u)) – Central detector row (v = 0) at each projection angle. projections_central[k, i] is the value at source angle geometry.angles[k] and detector column i.

  • geometry (ConeBeamGeometry) – Validated acquisition geometry.

  • n_s (int or None, optional) – Number of output parallel-beam detector bins. Defaults to geometry.n_u.

  • n_phi (int or None, optional) – Number of output parallel-beam angles uniformly distributed in [0°, 180°). Defaults to len(geometry.angles).

  • average_redundant (bool, optional) – Average the two views of each parallel ray available in a 360° scan. Default True.

Returns:

  • sino_parallel (ndarray, shape (n_s, n_phi)) – Rebinned parallel-beam sinogram, directly compatible with mod_iradon() and the registration functions in toupy.registration.registration.

  • s_coords (ndarray, shape (n_s,)) – Physical detector coordinates of the parallel-beam grid in mm. The centre of the array (index n_s // 2) corresponds to s = 0 (the rotation axis).

  • phi_angles (ndarray, shape (n_phi,)) – Parallel-beam projection angles in degrees in [0°, 180°).

Raises:

ValueError – If projections_central shape does not match geometry.

Notes

Why rebinning enables parallel-beam tools

Registration algorithms such as estimate_rot_axis() rely on the fact that a point object at (x₀, y₀) traces a sinusoid s(φ) = x₀ cos φ + y₀ sin φ in the parallel-beam sinogram. In a fan-beam sinogram the trace is a distorted sinusoid, so the parallel-beam tools do not apply directly. After rebinning, the sinusoidal structure is restored and the full parallel-beam registration toolbox becomes valid.

Jacobian note

The rebinned sinogram is not multiplied by the Jacobian \(\partial(u, \theta) / \partial(s, \phi)\) of the coordinate change. This is intentional: the Jacobian is only required when the rebinned data will be used for quantitative reconstruction (the Jacobian is 1 to first order near the central ray and is commonly omitted). For registration, only the positions of features matter, not their amplitudes.

Examples

Extract the central row, rebin to parallel, then run parallel-beam CoR estimation and display:

>>> import numpy as np
>>> from toupy.registration.cone_registration import (
...     ConeBeamRegistration, rebin_fan_to_parallel)
>>> from toupy.tomo.geometry import ConeBeamGeometry
>>>
>>> geom = ConeBeamGeometry(
...     SOD=500.0, SDD=1000.0, det_pixel_size=0.5,
...     n_u=128, n_v=96,
...     angles=np.linspace(0, 360, 180, endpoint=False),
... )
>>> geom.validate()
>>>
>>> # projections shape (n_angles, n_v, n_u)
>>> central_row = projections[:, geom.n_v // 2, :]   # (n_angles, n_u)
>>>
>>> sino_par, s_coords, phi_deg = rebin_fan_to_parallel(central_row, geom)
>>> # sino_par shape (n_u, n_angles), s_coords in mm, phi_deg in [0, 180)
>>>
>>> # Use directly with mod_iradon for the central slice:
>>> from toupy.tomo.iradon import mod_iradon
>>> recon_central = mod_iradon(sino_par, phi_deg)
>>>
>>> # Or use with parallel-beam registration to estimate CoR:
>>> # from toupy.registration.registration import estimate_rot_axis
>>> # estimate_rot_axis(sino_par, ...)

See also

ConeBeamRegistration.rebin_to_parallel

OOP wrapper.

ConeBeamRegistration.estimate_cor

CoR estimation directly in detector space (no rebinning required).

toupy.tomo.iradon.mod_iradon

Parallel-beam FBP reconstruction.

toupy.registration.registration.estimate_rot_axis

Interactive sinogram-based CoR estimation.

toupy.registration.registration module

registration_gd.py

Tomographic data alignment using gradient descent optimisation.

Vertical alignment : Adam optimiser on the L2 mass-fluctuation cost. Horizontal alignment: Adam optimiser on the L2 sinogram-consistency cost.

Both replace the original discrete line-search / parabolic-fit approach (_search_vshift_direction / _search_hshift_direction) with a proper gradient-based update rule. All other helper functions are unchanged.

toupy.registration.registration.alignprojections_horizontal(sinogram, theta, shiftstack, **params)[source]

Horizontal alignment of projections by tomographic consistency.

Shifts are optimised with a Newton gradient-descent algorithm that minimises the L2 distance between the measured sinogram and a synthetic sinogram computed from a filtered-back-projection reconstruction. Supports a multi-stage frequency-cutoff schedule and an optional spatial multiresolution warm-start.

Parameters:

  • sinogram (ndarray, shape (nr, nc)) – Measured sinogram (rows = detector pixels, columns = projections).

  • theta (ndarray, shape (nc,)) – Projection angles in radians.

  • shiftstack (ndarray, shape (2, n)) – Current shift estimates [vertical_shifts, horizontal_shifts]. Only the horizontal row (index 1) is updated.

  • **params –

    Algorithm parameters. Required keys:

    maxitint

    Maximum number of outer iterations per stage.

    pixtolfloat

    Convergence tolerance in pixels.

    freqcutofffloat

    Fraction of Nyquist used as the sinogram low-pass cut-off.

    shiftmethstr

    Interpolation method ('linear', 'fourier', 'spline').

    circlebool

    Apply a circular mask to the reconstruction.

    cliplowfloat or None

    Lower clipping threshold for the reconstruction values.

    cliphighfloat or None

    Upper clipping threshold for the reconstruction values.

    derivativesbool

    Whether projections are phase-gradient images.

    calc_derivativesbool

    Whether to compute derivatives on the synthetic sinogram.

    freqcutoff_schedulelist of float, optional

    Sequence of freqcutoff values (coarsest first). Each stage warm-starts the next. Default is [params['freqcutoff']] (single stage).

    multiresolutionbool, optional

    Spatial multiresolution warm-start at stage 0. Default False.

    mr_factorint, optional

    Down-sampling factor for the spatial warm-start. Default 2.

    n_coarse_iterint, optional

    Number of coarse iterations in the warm-start. Default 5.

    rtolfloat, optional

    Relative improvement threshold for early stopping. Default 0.

    silentbool, optional

    If True, suppress all matplotlib output (required for sub-process workers). Default False.

Returns:

shiftstack – Updated shift array with optimised horizontal shifts in row 1.

Return type:

ndarray, shape (2, n)

toupy.registration.registration.alignprojections_vertical(input_stack, shiftstack, **params)[source]

Vertical alignment of projections using the mass-fluctuation approach.

Shifts are optimised with a Newton gradient-descent algorithm that minimises the L2 distance between each projection’s vertical mass-fluctuation signal and the stack mean.

Parameters:

  • input_stack (ndarray, shape (n, nr, nc)) – Stack of projection images to be aligned.

  • shiftstack (ndarray, shape (2, n)) – Initial shift estimates [vertical_shifts, horizontal_shifts]. Modified in-place and also returned.

  • **params –

    Algorithm parameters. Required keys:

    maxitint

    Maximum number of outer iterations.

    pixtolfloat

    Convergence tolerance in pixels.

    deltaxint

    Horizontal margin (pixels) to exclude from the ROI.

    limsylist of int or None

    Explicit [row_start, row_end] row limits. None uses the full image height.

    shiftmethstr

    Interpolation method ('linear', 'fourier', 'spline').

    polyorderint

    Polynomial order for baseline removal in the fluctuation signal.

    alignxbool, optional

    If True, estimate horizontal shifts from the centre-of-mass before the vertical loop. Default False.

    subpixelbool, optional

    Ignored (always runs at sub-pixel precision). Default True.

Returns:

  • shiftstack (ndarray, shape (2, n)) – Optimised shift array.

  • output_stack (ndarray, shape (n, nr, nc)) – Vertically aligned projection stack.

toupy.registration.registration.center_of_mass_stack(input_stack, lims, shiftstack, shift_method='fourier')[source]

Compute the centre-of-mass position for each projection in the stack.

Parameters:

  • input_stack (ndarray, shape (n, nr, nc)) – Stack of projection images.

  • lims (tuple of array_like) – (limrow, limcol) — row and column index arrays defining the region of interest.

  • shiftstack (ndarray, shape (2, n)) – Current shift estimates used to pre-align each projection before computing the centre of mass.

  • shift_method (str, optional) – Interpolation method for the shift operation. Default 'fourier'.

Returns:

Array [centerx, centery] where centerx[i] and centery[i] are the horizontal and vertical centre-of-mass offsets (in pixels) for projection i.

Return type:

ndarray, shape (2, n)

toupy.registration.registration.compute_aligned_horizontal(input_stack, shiftstack, shift_method='linear', **kwargs)[source]

Horizontal-only alignment (copy variant).

**kwargs absorbs any extra keys when called as compute_aligned_horizontal(..., **params). shift_method defaults to params["shiftmeth"] if that key is present and shift_method was not supplied explicitly.

toupy.registration.registration.compute_aligned_sino(input_sino, shiftslice, shift_method='linear')[source]

Compute the aligned sinogram given the correction for object positions.

Parameters:

  • input_sino (array_like) – Sinogram to be shifted.

  • shiftslice (array_like) – Per-projection horizontal shifts (1, n).

  • shift_method (str) – ‘linear’, ‘fourier’, or ‘spline’.

Returns:

output_sino – Aligned sinogram.

Return type:

array_like

toupy.registration.registration.compute_aligned_stack(input_stack, shiftstack, shift_method='linear')[source]

Compute the aligned stack given the correction for object positions.

Parameters:

  • input_stack (array_like) – Stack of images to be shifted.

  • shiftstack (array_like) – Array of object motion corrections (2, n).

  • shift_method (str) – ‘linear’, ‘fourier’, or ‘spline’.

Returns:

output_stack – Aligned image stack.

Return type:

array_like

toupy.registration.registration.estimate_rot_axis(input_array, theta, **params)[source]

Interactively estimate the rotation-axis offset.

Opens a GUI figure showing the reconstructed slice and sinogram for the current offset (params["rot_axis_offset"]). Enter a new integer value in the text box and click Update to recompute. Click Confirm when satisfied.

Parameters:

  • input_array (ndarray, shape (n, nr, nc)) – Stack of derivative projections.

  • theta (ndarray) – Projection angles (degrees).

  • **params – Must contain: slicenum, rot_axis_offset, colormap, cliplow, cliphigh, sinolow, sinohigh, filtertype, freqcutoff, circle, algorithm, derivatives, calc_derivatives.

Returns:

  • Terminal mode – rot_axis_offset (int) — the confirmed offset value.

  • Jupyter mode (two-cell workflow) – picker (_RotAxisPicker) — access picker.rot_axis_offset in the next cell after clicking Confirm.

toupy.registration.registration.oneslicefordisplay(sinogram, theta, **params)[source]

Reconstruct and display one tomographic slice.

Pass updated values via params (e.g. params["freqcutoff"] = 0.5) and re-call to change reconstruction settings.

toupy.registration.registration.refine_horizontalalignment(input_stack, theta, shiftstack, **params)[source]

Run one refinement pass of horizontal alignment.

Refinement uses the current freqcutoff only — freqcutoff_schedule and multiresolution are ignored because the input shifts are already a good warm-start. Use params["rtol"] (e.g. 1e-3) to stop early when the improvement per iteration becomes negligible.

To run multiple refinement passes or change parameters between passes, update params in your script/notebook and call this function again.

Parameters:

  • input_stack (ndarray, shape (n, nr, nc)) – Stack of derivative projections.

  • theta (ndarray) – Projection angles.

  • shiftstack (ndarray, shape (2, n)) – Current shift array (modified in-place and returned).

  • **params – Alignment parameters forwarded to alignprojections_horizontal().

Returns:

  • shiftstack (ndarray)

  • params (dict)

toupy.registration.registration.register_2Darrays(image1, image2, subpixel=False)[source]

Image registration using phase cross-correlations.

Parameters:

  • image1 (array_like) – Reference image.

  • image2 (array_like) – Image to be shifted relative to image1.

  • subpixel (bool, optional) – If False (default) use pixel-precision registration. If True use sub-pixel precision (upsample factor = 100).

Returns:

  • shift (list of floats) – [row_shift, col_shift].

  • diffphase (float) – Phase difference between the two images.

  • offset_image2 (array_like) – Shifted image2 aligned to image1.

toupy.registration.registration.tomoconsistency_multiple(input_stack, theta, shiftstack, **params)[source]

Tomographic consistency alignment over multiple sinogram slices.

Each slice is aligned independently (using the same shiftstack warm-start) and the resulting per-slice horizontal shifts are averaged to produce a robust final estimate. Slices can be processed in parallel via ProcessPoolExecutor.

Parameters:

  • input_stack (ndarray, shape (n, nr, nc)) – Full 3-D projection stack.

  • theta (ndarray, shape (n,)) – Projection angles in radians.

  • shiftstack (ndarray, shape (2, n)) – Current shift estimates used as a warm-start for every slice. The horizontal row (index 1) is updated with the averaged result.

  • **params –

    Algorithm parameters forwarded to alignprojections_horizontal(). Extra keys:

    slicenumint

    Central slice index.

    n_slices_tcint, optional

    Number of slices to process (centred on slicenum). Default 10.

    n_workers_tcint, optional

    Number of parallel worker processes. Default max(1, cpu_count // 2). Pass 1 for sequential.

Returns:

shiftstack – Updated shift array with the averaged horizontal shifts in row 1 (or the original shifts if the user declines the result).

Return type:

ndarray, shape (2, n)

toupy.registration.registration.vertical_fluctuations(input_stack, lims, shiftstack, shift_method='fourier', polyorder=2)[source]

Compute the vertical mass-fluctuation signal for each projection.

The signal is obtained by integrating each (shifted) projection over the horizontal axis within the region of interest, then subtracting a polynomial baseline fit.

Parameters:

  • input_stack (ndarray, shape (n, nr, nc)) – Stack of projection images.

  • lims (tuple of array_like) – (rows, cols) — row and column index arrays for the ROI.

  • shiftstack (ndarray, shape (2, n)) – Current vertical and horizontal shift estimates.

  • shift_method (str, optional) – Interpolation method for the shift operation. Default 'fourier'.

  • polyorder (int, optional) – Polynomial order used for baseline removal. Default 2.

Returns:

vert_fluct – Array of vertical fluctuation signals, one row per projection.

Return type:

ndarray, shape (n, n_rows_roi)

toupy.registration.registration.vertical_shift(input_array, lims, vstep, maxshift, shift_method='linear', polyorder=2)[source]

Compute the vertical mass-fluctuation signal for a single projection shifted by a given amount.

Parameters:

  • input_array (ndarray, shape (nr, nc)) – Single projection image.

  • lims (tuple of array_like) – (rows, cols) — row and column index arrays for the ROI.

  • vstep (float) – Vertical shift to apply (pixels).

  • maxshift (int) – Safety margin (pixels) around the ROI to absorb border effects.

  • shift_method (str, optional) – Interpolation method for the shift operation. Default 'linear'.

  • polyorder (int, optional) – Polynomial order used for baseline removal. Default 2.

Returns:

shift_calc – Vertical fluctuation signal after shifting and polynomial baseline subtraction.

Return type:

ndarray, shape (n_rows_roi,)

toupy.registration.shift module

class toupy.registration.shift.ShiftFunc(**params)[source]

Bases: Variables

Collection of shift functions for 1-D and 2-D arrays.

The desired method is selected at construction time via the shiftmeth parameter and exposed through __call__ so that the same interface works regardless of the underlying algorithm.

Parameters:

**params –

Must contain:

shiftmethstr

Shift algorithm to use. One of 'linear' (bilinear interpolation), 'fourier' (FFT-based sub-pixel shift via pyFFTW), or 'spline' (spline interpolation via scipy.ndimage.interpolation.shift()).

__call__(*args)[source]

Apply the selected shift method to an array.

Parameters:

*args –

args[0]array_like

Input 1-D or 2-D array.

args[1]float or tuple of float

Shift amplitude in pixels. For 1-D arrays pass a scalar; for 2-D arrays pass (row_shift, col_shift).

args[2]str, optional

Padding mode (overrides the instance default).

args[3]bool, optional

True to return a complex array; False (default) to return a real array.

Returns:

Shifted array with the same shape as args[0]. Returns the input unchanged when the shift is zero.

Return type:

ndarray

shift_fft(input_array, shift)[source]

Performs pixel and subpixel shift (with wraping) using pyFFTW.

Since FFTW has efficient functions for array sizes which can be decompose in prime factor, the input_array is padded to the next fast size given by pyFFTW.next_fast_len. The padding is done in mode = ‘reflect’ by default to reduce border artifacts.

Parameters:

  • input_array (array_like) – Input image to calculate the shifts.

  • shift (int or tuple) – Number of pixels to shift. For 1D, use a integer value. For 2D, use a tuple of integers where the first value corresponds to shifts in the rows and the second value corresponds to shifts in the columns.

Returns:

output_array – Shifted image

Return type:

array_like

shift_linear(input_array, shift)[source]

Shifts an image with wrap around and bilinear interpolation

Parameters:

  • input_array (array_like) – Input image to calculate the shifts.

  • shift (int or tuple) – Number of pixels to shift. For 1D, use a integer value. For 2D, use a tuple of integers where the first value corresponds to shifts in the rows and the second value corresponds to shifts in the columns.

Returns:

output_array – Shifted image

Return type:

array_like

shift_spline_wrap(input_array, shift)[source]

Performs pixel and subpixel shift (with wraping) using splines

Parameters:

  • input_array (array_like) – Input image to calculate the shifts.

  • shift (int or tuple) – Number of pixels to shift. For 1D, use a integer value. For 2D, use a tuple of integers where the first value corresponds to shifts in the rows and the second value corresponds to shifts in the columns.

Returns:

output_array – Shifted image

Return type:

array_like