Reference manual

For the experimental operator-based interface, see Operator syntax.

Definition of geometry

A cornerstone in applying projection methods is to define for which geometry the projection has to be computed. Thus, the first step in using gratopy is always creating an instance of gratopy.ProjectionSettings defining the geometry, and thus internally precomputing relevant quantities.

class gratopy.ProjectionSettings(queue: CommandQueue, geometry: GeometryType, img_shape: tuple[int, int], angles: int | list[float] | ndarray | list[tuple[int | list[float] | ndarray, float, float]], n_detectors: int | None = None, angle_weights: float | list[float] | ndarray | None = None, detector_width: float = 2.0, image_width: float | None = None, R: float | None = None, RE: float | None = None, detector_shift: float = 0.0, midpoint_shift: tuple[float, float] = (0.0, 0.0), reverse_detector: bool = False)[source]

Creates and stores all relevant information concerning the projection geometry. Serves as a parameter for virtually all gratopy’s functions.

Parameters:

  • queue (pyopencl.CommandQueue) – The OpenCL command queue with which the computations are associated.

  • geometry (GeometryType) – Determines whether parallel beam (gratopy.RADON) or fanbeam geometry (gratopy.FANBEAM) is considered.

  • img_shape (tuple \((N_x,N_y)\)) – The number of pixels of the image in x- and y-direction respectively, i.e., the image dimension. It is assumed that by default, the center of rotation is in the middle of the grid of quadratic pixels. The midpoint can, however, be shifted, see the midpoint_shift parameter.

  • angles (int, list[float] / numpy.ndarray, list[tuple(int/list[float]/numpy.ndarray, float, float)]) –

    Determines which angles are considered for the projection. An integer is interpreted as the number \(N_a\) of uniformly distributed angles in the angular range \([0,\pi[\), \([0,2\pi[\) for Radon and fanbeam transform, respectively, where for negative integers the same angles according to its modulus but with reversed order are generated. Alternatively, the angles can be given explicitly as a list or numpy.ndarray. These two options also imply a full angle setting (as opposed to limited angle setting).

    A limited angle setting can be specified in two ways. First, a list of angular range sections can be passed as input. An angular range section is a tuple with either an integer or a list/array of angles (first element) together with a pair specifying the lower and upper bound of the angular range interval (second and third element), i.e., of type tuple(int, float, float) or tuple(list[float], float, float). If the first element is an integer, the angular interval will be uniformly partitioned into the modulus number of angles (note that the first and last angles are not the lower/upper bounds to ensure uniform angle weights) again in increasing or decreasing order, depending on the sign. Otherwise, a list or array specifying the individual angles is expected. In particular, multiple angular sections can be specified, by passing a list of angular range sections.

    Alternatively, one can use a list of angles and set angle_weigths (see below) manually to suitable values by passing a scalar, a list or an array.

  • n_detectors (int, default None) – The number \(N_s\) of (equi-spaced) detector pixels considered. When None, \(N_s\) will be chosen as \(\sqrt{N_x^2+N_y^2}\).

  • angle_weights (None, float, list[float] or numpy.ndarray, default None) – The weights \((\Delta_a)_a\) associated with the angles, which influences the weighting of the rays for the backprojection. See Adjointness in gratopy for a more detailed description. If None (by default), the weights are computed automatically based on the angles parameter. In the full angle setting, this automatism considers a partition of the half circle for parallel beam and the full circle for fanbeam geometry based on the given angles and sets the angle weight to the average of the distances from of an angle to its two neighbors (in the sense of a circle). Similarly, in the limited angle case, each angle section is partitioned by the angles associated with this section and the weights are chosen taking additionally the boundary of the section into account. In case of a scalar input, this scalar will be used as the (constant) angle weight for all angles. Further, all angle weights can directly be set by passing an input of type list[float] or numpy.ndarray of suitable length.

  • detector_width (float, default 2.0) – Physical length of the detector. For standard Radon transformation this can usually remain fixed at the default value (together with image_width).

  • image_width (float, default None) – Physical size of the image indicated by the length of the longer side of the rectangular image domain. For parallel beam geometry, when None, image_width is chosen as 2.0. For fanbeam geometry, when None, image_width is chosen such that the projections exactly capture the image domain. To illustrate, choosing image_width = detector_width results in the standard Radon transform with each projection touching the entire object, while img_width = 2 detector_width results in each projection capturing only half of the image.

  • R (float, must be set for fanbeam geometry) – Physical (orthogonal) distance from source to detector line. Has no impact for parallel beam geometry.

  • RE (float, must be set for fanbeam geometry) – Physical distance from source to origin (center of rotation). Has no impact for parallel beam geometry.

  • detector_shift (float, default 0.0) – Physical shift of all detector pixels along the detector line. Defaults to the application of no shift, i.e., the detector pixels span the range [-detector_width/2, detector_width/2].

  • midpoint_shift (list, default [0.0, 0.0]) – Two-dimensional vector representing the shift of the image away from center of rotation. Defaults to the application of no shift, i.e., the center of rotation is also the center of the image.

  • reverse_detector (bool, default False) – When True, the detector direction is flipped in case of fanbeam geometry, i.e., the positive and negative detector positions are swapped. This parameter has no effect for parallel geometry. When activated together with swapping the sign of the angles, this has the same effect for projection as mirroring the image.

These input parameters create attributes of the same name in an instance of ProjectionSettings, though the corresponding values might be slightly restructured by internal processes. Further useful attributes are listed below. It is advised not to set these attributes directly but rather to choose suitable input parameters for the initialization.

Variables:

  • is_parallel (bool) – True if the geometry is for parallel beams, False otherwise.

  • is_fan (bool) – True if the geometry is for fanbeam geometry, False otherwise.

  • angles (numpy.ndarray) – Angles from which projections are considered.

  • n_angles (int) – Number of all angles \(N_a\).

  • sinogram_shape (tuple \((N_s,N_a)\)) – Represents the number of considered detectors (n_detectors) and angles (n_angles).

  • delta_x (float) – Physical width and height \(\delta_x\) of the image pixels.

  • delta_s (float) – Physical width \(\delta_s\) of a detector pixel.

  • delta_ratio (float) – Ratio \({\delta_s}/{\delta_x}\), i.e. the detector pixel width relative to unit image pixels.

  • angle_weights (numpy.ndarray) – Represents the angular discretization width for each angle which are used to weight the projections, see parameter angle_weights above. When none was given as input, the angle_weights chosen by the automatism will be written to this variable.

  • prg (gratopy.Program) – OpenCL program containing the gratopy OpenCL kernels. For the corresponding code, see gratopy.create_code

  • struct (dict see radon_struct() and fanbeam_struct() returns) – Data used in the projection operator. Contains in particular dictionaries of numpy.ndarray associated to precision single and double with the angular information necessary for computations.

backprojection(img: Array, sino: Array, projectionsetting: ProjectionSettings, wait_for: list[Event] = [])[source]

Calls the appropriate back projection given the configured geometry.

create_sparse_matrix(dtype: DTypeLike = dtype('float32'), order: Literal['C', 'F'] = 'F') → coo_matrix[source]

Creates a sparse matrix representation of the associated forward projection.

Parameters:

  • dtype (numpy.dtype, default numpy.float32) – Precision to compute the sparse representation in.

  • order (str, default F) – Contiguity of the image and sinogram array to the transform, can be F or C.

Returns:

Sparse matrix corresponding to the forward projection.

Return type:

scipy.sparse.coo_matrix

Note that for high resolution projection operators, this may require infeasibly much time and memory.

forwardprojection(sino: Array, img: Array, projectionsetting: ProjectionSettings, wait_for: list[Event] = []) → None[source]

Calls the appropriate forward projection given the configured geometry.

property geometry: GeometryType

Returns the configured geometry.

get_projection_kernel(sinogram: Array, image: Array, adjoint: bool = False) → Kernel[source]

Returns the compiled projection function for the given sinogram and image. If adjoint is True, the adjoint projection function is returned.

Parameters:
Returns:

The projection function.

Return type:

callable

property is_fan: bool

Returns whether the projection settings are for fanbeam geometry.

property is_parallel: bool

Returns whether the projection settings are for parallel beam geometry.

show_geometry(angle: float, figure: Figure | None = None, axes: Axes | None = None, show: bool = True)[source]

Visualize the geometry associated with the projection settings. This can be useful in checking that indeed, the correct input for the desired geometry was given.

Parameters:

  • angle (float) – The angle for which the projection is considered.

  • figure (matplotlib.figure.Figure, default None) – Figure in which to plot. If neither figure nor axes are given, a new figure (figure(0)) will be created.

  • axes (matplotlib.axes.Axes, default None) – Axes to plot into. If None, a new axes inside the figure is created.

  • show (bool, default True) – Determines whether the resulting plot is immediately shown (True). If False, matplotlib.pyplot.show() can be used at a later point to show the figure.

Returns:

Figure and axes in which the geometry visualization is plotted.

Return type:

tuple(matplotlib.figure.Figure, matplotlib.axes.Axes)

Transforms

The functions forwardprojection() and backprojection() perform the projection operations based on the geometry defined in projectionsetting. The images img and the sinograms sino need to be interpreted and behave as described in Getting started.

gratopy.forwardprojection(img: Array, projectionsetting: ProjectionSettings, sino: Array | None = None, wait_for: list[Event] = [])[source]

Performs the forward projection (either for the Radon or the fanbeam transform) of a given image using the given projection settings.

Parameters:

  • img (pyopencl.array.Array with compatible dimensions) – The image to be transformed.

  • projectionsetting (gratopy.ProjectionSettings) – The geometry settings for which the forward transform is computed.

  • sino (pyopencl.array.Array with compatible dimensions, default None) – The array in which the result of transformation is saved. If None (per default) is given, a new array will be created and returned.

  • wait_for (list[pyopencl.Event], default []) – The events to wait for before performing the computation in order to avoid, e.g., race conditions, see pyopencl.Event. This program will always wait for img.events and sino.events (so you need not add them to wait_for).

Returns:

The sinogram associated with the projection of the image. If the sino is not None, the same pyopencl array is returned with the values in its data overwritten.

Return type:

pyopencl.array.Array

The forward projection can be performed for single or double precision arrays. The dtype (precision) of img and sino (if given) have to coincide and the output will be of the same precision. It respects any combination of C and F contiguous arrays where output will be of the same contiguity as img if no sino is given. The OpenCL events associated with the transform will be added to the output’s events. In case the output array is created, it will use the allocator of img. If the image and sinogram have a third dimension (z-direction) the operator is applied slicewise.

gratopy.backprojection(sino: Array, projectionsetting: ProjectionSettings, img: Array | None = None, wait_for: list[Event] = [])[source]

Performs the backprojection (either for the Radon or the fanbeam transform) of a given sinogram using the given projection settings.

Parameters:

  • sino (pyopencl.array.Array with compatible dimensions) – Sinogram to be backprojected.

  • projectionsetting (gratopy.ProjectionSettings) – The geometry settings for which the forward transform is computed.

  • img (pyopencl.array.Array with compatible dimensions, default None) – The array in which the result of backprojection is saved. If None is given, a new array will be created and returned.

  • wait_for (list[pyopencl.Event], default []) – The events to wait for before performing the computation in order to avoid, e.g., race conditions, see pyopencl.Event. This program will always wait for img.events and sino.events (so you need not add them to wait_for).

Returns:

The image associated with the backprojected sinogram, coinciding with the img if not None, with the values in its data overwritten.

Return type:

pyopencl.array.Array

The backprojection can be performed for single or double precision arrays. The dtype (precision) of img and sino have to coincide. If no img is given, the output precision coincides with sino’s. The operation respects any combination of C and F contiguous arrays, where if img is None, the result’s contiguity coincides with sino’s. The OpenCL events associated with the transform will be added to the output’s events. In case the output array is created, it will use the allocator of sino. If the sinogram and image have a third dimension (z-direction), the operator is applied slicewise.

Solvers

Based on these forward and backward operators, one can implement a variety of reconstruction algorithms, where the toolbox’s focus is on iterative methods (as those in particular are dependent on efficient implementation). The following constitute a few easy-to-use examples which also serve as illustration on how gratopy can be included in custom pyopencl implementations.

gratopy.landweber(sino: Array, projectionsetting: ProjectionSettings, number_iterations: int = 100, w: float = 1.0) → Array[source]

Performs a Landweber iteration [L1951] to approximate a solution to the image reconstruction problem associated with a projection and sinogram. This method is also known as SIRT.

Parameters:

  • sino (pyopencl.array.Array) – Sinogram data to reconstruct from.

  • projectionsetting (gratopy.ProjectionSettings) – The geometry settings for which the projection is considered.

  • number_iterations (int, default 100) – Number of iteration steps to be performed.

  • w (float, default 1) – Relaxation parameter weighted by the norm of the projection operator (w<1 guarantees convergence).

Returns:

Reconstruction from given sinogram gained via Landweber iteration.

Return type:

pyopencl.array.Array

[L1951]

Landweber, L. “An iteration formula for Fredholm integral equations of the first kind.” Amer. J. Math. 73, 615–624 (1951). https://doi.org/10.2307/2372313

gratopy.conjugate_gradients(sino: Array, projectionsetting: ProjectionSettings, number_iterations: int = 20, epsilon: float = 0.0, x0: Array | None = None)[source]

Performs a conjugate gradients iteration [HS1952] to approximate a solution to the image reconstruction problem associated with a projection and sinogram.

Parameters:

  • sino (pyopencl.array.Array) – Sinogram data to invert.

  • projectionsetting (gratopy.ProjectionSettings) – The geometry settings for which the projection is considered.

  • number_iterations (float, default 20) – Maximal number of iteration steps to be performed.

  • x0 (pyopencl.array.Array, default None) – Initial guess for iteration (defaults to zeros if None).

  • epsilon (float, default 0.00) – Tolerance parameter, the iteration stops if relative residual<epsilon.

Returns:

Reconstruction gained via conjugate gradients iteration.

Return type:

pyopencl.array.Array

[HS1952]

Hestenes, M. R., Stiefel, E. “Methods of Conjugate Gradients for Solving Linear Systems.” Journal of Research of the National Bureau of Standards, 49:409–436 (1952). https://doi.org/10.6028/jres.049.044

gratopy.total_variation(sino: Array, projectionsetting: ProjectionSettings, mu: float, number_iterations: int = 1000, slice_thickness: float = 1.0, stepsize_weighting: float = 10.0)[source]

Performs a primal-dual algorithm [CP2011] to solve a total-variation regularized reconstruction problem associated with a given projection operator and sinogram. This corresponds to the approximate solution of \(\min_{u} {\frac\mu2}\|\mathcal{P}u-f\|_{L^2}^2+\mathrm{TV}(u)\) for \(\mathcal{P}\) the projection operator, \(f\) the sinogram and \(\mu\) a positive regularization parameter (i.e., an \(L^2-\mathrm{TV}\) reconstruction approach).

Parameters:

  • sino (pyopencl.array.Array) – Sinogram data to invert.

  • projectionsetting (gratopy.ProjectionSettings) – The geometry settings for which the projection is considered.

  • mu (float) – Regularization parameter, the smaller the stronger the applied regularization.

  • number_iterations (int, default 1000) – Number of iterations to be performed.

  • slice_thickness (float, default 1.0, i.e., isotropic voxels) – When 3-dimensional data sets are considered, regularization is also applied across slices. This parameter represents the ratio of the slice thickness to the length of one pixel within a slice. The choice slice_thickness =0 results in no coupling across slices.

  • stepsize_weighting (float, default 10.0) – Allows to weight the primal-dual algorithm’s step sizes \(\sigma\) (stepsize for dual update) and \(\tau\) (stepsize for primal update) (with \(\sigma\tau\|\mathcal{P}\|^2\leq 1\)) by multiplication and division, respectively, with the given value.

Returns:

Reconstruction gained via primal-dual iteration for the total-variation regularized reconstruction problem.

Return type:

pyopencl.array.Array

[CP2011]

Chambolle, A., Pock, T. “A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging.” J Math Imaging Vis 40, 120–145 (2011). https://doi.org/10.1007/s10851-010-0251-1

gratopy.normest(projectionsetting: ProjectionSettings, number_iterations: int = 50, dtype: DTypeLike = 'float32', allocator: AllocatorBase = None)[source]

Estimate the spectral norm of the projection operator via power iteration, i.e., the operator norm with respect to the norms discussed in section concerning adjointness. Useful for iterative methods that require such an estimate, e.g., landweber() or total_variation().

Parameters:

  • projectionsetting (gratopy.ProjectionSettings) – The geometry settings for which the projection is considered.

  • number_iterations (int, default 50) – The number of iterations to terminate after.

  • dtype (numpy.dtype, default numpy.float32) – Precision for which to apply the projection operator (which is not supposed to impact the estimate significantly).

Returns:

An estimate of the spectral norm for the projection operator.

Return type:

float

gratopy.weight_sinogram(sino: Array, projectionsetting: ProjectionSettings, sino_out: Array | None = None, divide: bool = False, wait_for: list[Event] = []) → Array[source]

Performs an angular rescaling of a given sinogram via multiplication (or division) with the projection’s angle weights (size of projections in angle dimension, see attributes of ProjectionSettings) to the respective projections. This can be useful, e.g., for computing norms or dual pairings in the appropriate Hilbert space.

Parameters:

  • sino (pyopencl.array.Array) – The sinogram whose rescaling is computed. This array itself remains unchanged unless the same array is given as sino_out.

  • projectionsetting (gratopy.ProjectionSettings) – The geometry settings for which the rescaling is computed.

  • sino_out (pyopencl.array.Array default None) – The array in which the result of rescaling is saved. If None (per default) is given, a new array will be created and returned. When giving the same array as sino, the values in sino will be overwritten.

  • divide (bool, default False) – Determines whether the sinogram is divided (instead of multiplied) by the angular weights. If True, a division is performed, otherwise, the weights are multiplied.

  • wait_for (list[pyopencl.Event], default []) – The events to wait for before performing the computation in order to avoid, e.g., race conditions, see pyopencl.Event. This program will always wait for img.events and sino.events (so you need not add them to wait_for).

Returns:

The weighted sinogram. If sino_out is not None, it is returned with the values in its data overwritten. In particular, giving the same array for sino and sino_out will overwrite this array.

Return type:

pyopencl.array.Array

Data generation

For convenient testing, a phantom generator is included which creates a modified two-dimensional phantom of arbitrary size.

gratopy.phantom(queue: CommandQueue, N: int | tuple[int, int] | ndarray, modified: bool = True, E: ndarray | None = None, ret_E: bool = False, dtype: str | DTypeLike = 'double', allocator: AllocatorBase | None = None)

Generate an OpenCL Shepp-Logan phantom of size (N, N).

Parameters:

  • queue (pyopencl.CommandQueue) – The OpenCL command queue.

  • N (int or array_like) – Matrix size, (N, N) or (M, N).

  • modified (bool) – Use original gray-scale values as given in [SL1974]. Most implementations use modified values for better contrast (for example, see [2] and [3]).

  • E (array_like or None) –

    \(e \times 6\) numeric matrix defining \(e\) ellipses. The six columns of E are:

    • Gray value of the ellipse (in [0, 1])

    • Length of the horizontal semi-axis of the ellipse

    • Length of the vertical semi-axis of the ellipse

    • x-coordinate of the center of the ellipse (in [-1, 1])

    • y-coordinate of the center of the ellipse (in [-1, 1])

    • Angle between the horizontal semi-axis of the ellipse and the x-axis of the image (in rad)

  • ret_E (bool) – Return the matrix E used to generate the phantom.

  • dtype (str or numpy.dtype) – The pyopencl data type in which the phantom is created.

  • allocator (An implementation of pyopencl.tools.AllocatorInterface or None) – The pyopencl allocator used for memory allocation.

Returns:

Phantom/parameter pair (ph [, E]).

Variables:

  • ph (pyopencl.array.Array) – The Shepp-Logan phantom.

  • E (array_like, optional) – The ellipse parameters used to generate ph.

This much abused phantom is due to [SL1974]. The tabulated values in the paper are reproduced in the Wikipedia entry [1]. The original values do not produce great contrast, so modified values are used by default (see Table B.1 in [TS1996] or implementations [2] and [3]).

[SL1974] (1,2)

Shepp, Lawrence A., and Benjamin F. Logan. “The Fourier reconstruction of a head section.” IEEE Transactions on nuclear science 21.3 (1974): 21-43.

[TS1996]

Toft, Peter Aundal, and John Aasted Sørensen. “The Radon transform-theory and implementation.” (1996).

Internal functions

The following contains the documentation for a set of internal functions which could be of interest for developers. Note that these might be subject to change in the future.

gratopy.radon(sino: Array, img: Array, projectionsetting: ProjectionSettings, wait_for: list[Event] = [])[source]

Performs the Radon transform of a given image using the given projectionsetting.

Parameters:

  • sino (pyopencl.array.Array) – The array in which the resulting sinogram is written.

  • img (pyopencl.array.Array) – The image to transform.

  • projectionsetting (gratopy.ProjectionSettings) – The geometry settings for which the Radon transform is performed.

  • wait_for (list[pyopencl.Event], default []) – The events to wait for before performing the computation in order to avoid, e.g., race conditions, see pyopencl.Event. This program will always wait for img.events and sino.events (so you need not add them to wait_for).

Returns:

Event associated with the computation of the Radon transform (which is also added to the events of sino).

Return type:

pyopencl.Event

gratopy.radon_ad(img: Array, sino: Array, projectionsetting: ProjectionSettings, wait_for: list[Event] = [])[source]

Performs the Radon backprojection of a given sinogram using the given projectionsetting.

Parameters:

  • img (pyopencl.array.Array) – The array in which the resulting backprojection is written.

  • sino (pyopencl.array.Array) – The sinogram to transform.

  • projectionsetting (gratopy.ProjectionSettings) – The geometry settings for which the Radon backprojection is performed.

  • wait_for (list[pyopencl.Event], default []) – The events to wait for before performing the computation in order to avoid, e.g., race conditions, see pyopencl.Event. This program will always wait for img.events and sino.events (so you need not add them to wait_for).

Returns:

Event associated with the computation of the Radon backprojection (which is also added to the events of img).

Return type:

pyopencl.Event

gratopy.radon_struct(queue: CommandQueue, img_shape: tuple[int, int], angles: ndarray, angle_weights: ndarray, n_detectors: int | None = None, detector_width: float = 2.0, image_width: float = 2.0, midpoint_shift: tuple[float, float] = (0.0, 0.0), detector_shift: float = 0.0)[source]

Creates the structure storing geometry information required for the Radon transform and its adjoint.

Parameters:

  • queue (pyopencl.CommandQueue) – OpenCL command queue in which context the computations are to be performed.

  • img_shape (tuple \((N_x,N_y)\)) – The number of pixels of the image in x- and y-direction respectively, i.e., the image size. It is assumed that by default, the center of rotation is in the middle of the grid of quadratic pixels. The midpoint can, however, be shifted, see the midpoint_shift parameter.

  • angles (numpy.ndarray) – Determines which angles are considered for the projection.

  • angle_weights (numpy.ndarray) – The weights associated to the angles, e.g., how much of the angular range is covered by this angle. This impacts the weighting of rays for the backprojection.

  • n_detectors (int, default None) – The number \(N_s\) of considered (equi-spaced) detectors. If None, \(N_s\) will be chosen as \(\sqrt{N_x^2+N_y^2}\).

  • detector_width (float, default 2.0) – Physical length of the detector line.

  • image_width (float, default 2.0) – Physical size of the image indicated by the length of the longer side of the rectangular image domain. Choosing image_width = detector_width results in the standard Radon transform with each projection touching the entire object, while img_width = 2 detector_width results in each projection capturing only half of the image.

  • midpoint_shift (list[float], default [0.0, 0.0]) – Two-dimensional vector representing the shift of the image away from center of rotation. Defaults to the application of no shift.

  • detector_shift (float, default 0.0) – Physical shift of the detector along the detector line in detector pixel offsets. Defaults to the application of no shift, i.e., the detector reaches from [- detector_width/2, detector_width/2].

Returns:

Struct dictionary with the following variables as entries, where the keys are strings of the same names:

Return type:

dict

Variables:

  • ofs_dict (dict{numpy.dtype: numpy.ndarray}) –

    Dictionary containing the relevant angular information as numpy.ndarray for the data types numpy.float32 and numpy.float64. The arrays have dimension \((8, N_a)\) with columns:

    0

    weighted cosine

    1

    weighted sine

    2

    detector offset

    3

    inverse of cosine/sine

    4

    angular weight

    5

    flipped

    The remaining columns are unused. The value flipped indicates whether the x and y axis are flipped (1) or not (0), which is done for reasons of numerical stability. The 4th entry contains the inverse of sine if the axes are flipped and the inverse of cosine otherwise.

  • shape (tuple) – Tuple of integers \((N_x,N_y)\) representing the size of the image.

  • sinogram_shape (tuple) – Tuple of integers \((N_s,N_a)\) representing the size of the sinogram.

  • geo_dict (dict{numpy.dtype: numpy.ndarray}) – Dictionary mapping the allowed data types to an array containing the values [\(\delta_x, \delta_s, N_x, N_y, N_s, N_a\)].

  • angles_diff_buf – Dictionary containing the same values as in ofs_dict [4] representing the weights associated with the angles (i.e., the length of sinogram pixels in the angular direction).

gratopy.fanbeam(sino: Array, img: Array, projectionsetting: ProjectionSettings, wait_for: list[Event] = [])[source]

Performs the fanbeam transform of a given image using the given projectionsetting.

Parameters:

  • sino (pyopencl.array.Array) – The array in which the resulting sinogram is written.

  • img (pyopencl.array.Array) – The image to transform.

  • projectionsetting (gratopy.ProjectionSettings) – The geometry settings for which the fanbeam transform is performed.

  • wait_for (list[pyopencl.Event], default []) – The events to wait for before performing the computation in order to avoid, e.g., race conditions, see pyopencl.Event. This program will always wait for img.events and sino.events (so you need not add them to wait_for).

Returns:

Event associated with the computation of the fanbeam transform (which is also added to the events of sino).

Return type:

pyopencl.Event

gratopy.fanbeam_ad(img: Array, sino: Array, projectionsetting: ProjectionSettings, wait_for: list[Event] = [])[source]

Performs the fanbeam backprojection of a given sinogram using the given projectionsetting.

Parameters:

  • img (pyopencl.array.Array) – The array in which the resulting backprojection is written.

  • sino (pyopencl.array.Array) – The sinogram to transform.

  • projectionsetting (gratopy.ProjectionSettings) – The geometry settings for which the fanbeam backprojection is performed.

  • wait_for (list[pyopencl.Event], default []) – The events to wait for before performing the computation in order to avoid, e.g., race conditions, see pyopencl.Event. This program will always wait for img.events and sino.events (so you need not add them to wait_for).

Returns:

Event associated with the computation of the fanbeam backprojection (which is also added to the events of img).

Return type:

pyopencl.Event

gratopy.fanbeam_struct(queue: CommandQueue, img_shape: tuple[int, int], angles: ndarray, detector_width: float, source_detector_dist: float, source_origin_dist: float, angle_weights: ndarray, n_detectors: int | None = None, detector_shift: float = 0.0, image_width: float | None = None, midpoint_shift: tuple[float, float] = (0.0, 0.0), reverse_detector: bool = False)[source]

Creates the structure storing geometry information required for the fanbeam transform and its adjoint.

Parameters:

  • queue (pyopencl.CommandQueue) – OpenCL command queue in which context the computations are to be performed.

  • img_shape (tuple \((N_x,N_y)\)) – The number of pixels of the image in x- and y-direction respectively, i.e., the image size. It is assumed that by default, the center of rotation is in the middle of the grid of quadratic pixels. The midpoint can, however, be shifted, see the midpoint_shift parameter.

  • angles (numpy.ndarray) – Determines which angles are considered for the projection.

  • detector_width (float, default 2.0) – Physical length of the detector line.

  • source_detector_dist (float) – Physical (orthogonal) distance R from the source to the detector line.

  • source_origin_dist (float) – Physical distance RE from the source to the origin (center of rotation).

  • angle_weights (numpy.ndarray) – The weights associated to the angles, e.g., how much of the angular range is covered by this angle. This impacts the weighting of rays for the backprojection.

  • n_detectors (int or None, default None) – The number \(N_s\) of considered (equi-spaced) detectors. If None, \(N_s\) will be chosen as \(\sqrt{N_x^2+N_y^2}\).

  • detector_shift (float, default 0.0) – Physical shift of the detector along the detector line in detector pixel offsets. Defaults to the application of no shift, i.e., the detector reaches from [- detector_width/2, detector_width/2].

  • image_width (float, default None) – Physical size of the image indicated by the length of the longer side of the rectangular image domain. If None, image_width is chosen to capture just all rays.

  • midpoint_shift (list[float], default [0.0, 0.0]) – Two-dimensional vector representing the shift of the image away from center of rotation. Defaults to the application of no shift.

  • reverse_detector (bool, default False) – When True, the detector direction is flipped.

Returns:

Struct dictionary with the following variables as entries, where the keys are strings of the same names:

Return type:

dict

Variables:

  • img_shape (tuple) – Tuple of integers \((N_x,N_y)\) representing the size of the image.

  • sinogram_shape (tuple) – Tuple of integers \((N_s,N_a)\) representing the size of the sinogram.

  • ofs_dict (dict{numpy.dtype: numpy.ndarray}) –

    Dictionary containing the relevant angular information as numpy.ndarray for the data types numpy.float32 and numpy.float64. The arrays have dimension \((8, N_a)\) with columns:

    0 1

    vector of length \(\delta_s\) pointing in positive detector direction

    2 3

    vector connecting source and center of rotation

    4 5

    vector connection the origin and its projection onto the detector line

    6

    angular weight

    The remaining column is unused.

  • sdpd_dict (dict{numpy.dtype: numpy.ndarray}) – Dictionary mapping numpy.float32 and numpy.float64 to a numpy.ndarray representing the values \(\sqrt{(s^2+R^2)}\) for the weighting in the fanbeam transform (weighted by delta_ratio, i.e., \(\delta_s/\delta_x\)).

  • image_width – Physical size of the image. Equal to the input parameter if given, or to the determined image size if image_width is None (see parameter image_width).

  • geo_dict (dict{numpy.dtype: numpy.ndarray}) – Dictionary mapping the allowed data types to an array containing the values [source detector distance/\(\delta_x\), source origin distance/\(\delta_x\), width of a detector_pixel relative to width of image_pixels i.e. \(\delta_s\)/\(\delta_x\), image midpoint x-coordinate (in pixels), image midpoint y-coordinate (in pixels), detector line midpoint (in detector-pixels), \(N_x\), \(N_y\), \(N_s\), \(N_a\), width of an image pixel (\(\delta_x\))].

  • angles_diff (dict{numpy.dtype: numpy.ndarray}) – Dictionary containing the same values as in ofs_dict [6] representing the weights associated with the angles (i.e., the length of sinogram pixels in the angular direction).

gratopy.create_code(cl_context: Context | None = None)[source]

Reads and creates CL code containing all OpenCL kernels of the gratopy toolbox.

Parameters:

cl_context – The OpenCL context in which the code is to be compiled. Used for checking support for double precision floating point numbers.

Returns:

The toolbox’s CL code.

Return type:

str

Operator package

The operator package provides an experimental operator-oriented interface to a subset of gratopy’s functionality.

Operator algebra

Generic implementation of operators including basic arithmetic.

class gratopy.operator.base.Operator(name: str | None = None, scalar: float | int | int64 | float32 | float64 = 1, state: dict[str, Any] | None = None, arithmetic_operation: OperatorArithmeticOperation | None = None, operands: list[Operator] | None = None, input_shape: tuple[int, ...] | None = None, output_shape: tuple[int, ...] | None = None)[source]

Bases: object

Base class for all operators.

apply_to(argument: Any, output: Any | None = None, **kwargs: Any) → Any[source]

Application of this operator to some given argument.

For composite products, output is passed only to the final child operator. Intermediate results are computed normally while the final result can reuse the caller-provided array.

is_composite() → bool[source]

Check if the operator is composite.

property scalar: float | int | int64 | float32 | float64
class gratopy.operator.base.OperatorArithmeticOperation(*values)[source]

Bases: Enum

Enum for operations that can be performed on operators.

ADDITION = 'sum'
MULTIPLICATION = 'prod'

Projection operators

Gratopy projection operators.

This module contains concrete operator implementations for gratopy’s experimental operator interface, most notably Radon.

The focus is currently on a compositional, operator-style interface for Radon transforms and their adjoints. Fanbeam support is not yet implemented at the same level.

Examples

>>> import numpy as np
>>> import pyopencl as cl
>>> import gratopy
>>> ctx = cl.create_some_context(interactive=False)
>>> queue = cl.CommandQueue(ctx)
>>> Nx = 300
>>> img = np.zeros((Nx, Nx), dtype=np.float32)
>>> R = gratopy.operator.Radon(image_domain=Nx, angles=180)
>>> sinogram = R.apply_to(img, queue=queue)
>>> backprojection = R.T.apply_to(sinogram)

The same forward and adjoint applications can also be written via operator syntax:

>>> sinogram = R * img
>>> backprojection = R.T * sinogram
class gratopy.operator.projection.Fanbeam(source_distances: float | tuple[float, float], image_domain: int | tuple[int, int] | ImageDomain, angles: Angles | int, detectors: Detectors | int | None = None, adjoint: bool = False)[source]

Bases: Operator

class gratopy.operator.projection.Radon(image_domain: int | tuple[int, int] | ImageDomain, angles: Angles | int, detectors: Detectors | int | None = None, adjoint: bool = False, kernel_spec: OpenCLKernelSpec | None = None)[source]

Bases: _OpenCLOperator

Parallel-beam Radon transform operator.

This class provides the main entry point to gratopy’s experimental operator-based projection interface. A Radon object represents either the forward projection operator or, when created through T, its adjoint.

Parameters

image_domain:

Description of the image grid and its physical extent. This can be given as:

Plain integer and tuple inputs are converted to ImageDomain with a default extent of 2.0.

angles:

Angular sampling of the operator. This can be given either as an integer or as an explicit gratopy.utilities.Angles object.

If an integer is given, uniformly weighted angles are created via gratopy.utilities.Angles.uniform() with the half_circle parameter set to True.

detectors:

Detector configuration. This can be given as:

  • None to use the default detector count ceil(hypot(Nx, Ny)),

  • an integer specifying the number of detector pixels,

  • or an explicit gratopy.utilities.Detectors object.

Plain integer inputs are converted to Detectors with default extent handling.

adjoint:

If False (the default), construct the forward Radon transform. If True, construct the adjoint operator directly. In practice, the adjoint is typically accessed via T.

kernel_spec:

Optional gratopy.operator.opencl.OpenCLKernelSpec describing which OpenCL kernel bundle to use. When omitted, the operator uses the Radon kernels shipped with gratopy.

Notes

The operator accepts both pyopencl.array.Array inputs and NumPy arrays. When applying the operator to a NumPy array, a queue must be available so that the data can be transferred to the device.

The operator supports 2D as well as slicewise 3D data. For example, a forward operator with image shape (Nx, Ny) maps:

  • (Nx, Ny) to (Ns, Na),

  • (Nx, Ny, Nz) to (Ns, Na, Nz).

Correspondingly, the adjoint maps:

  • (Ns, Na) to (Nx, Ny),

  • (Ns, Na, Nz) to (Nx, Ny, Nz).

Extent placeholders are supported experimentally for Radon geometry when exactly one extent is fixed numerically and the other is given as gratopy.utilities.ExtentPlaceholder. The operator can resolve an image extent from a fixed detector extent, or a detector extent from a fixed image extent. Passing placeholders for both extents at once remains unsupported and raises NotImplementedError; geometries for which no valid placeholder resolution exists raise ValueError.

Examples

>>> import numpy as np
>>> import pyopencl as cl
>>> import gratopy
>>> ctx = cl.create_some_context(interactive=False)
>>> queue = cl.CommandQueue(ctx)
>>> img = np.zeros((128, 128), dtype=np.float32)
>>> R = gratopy.operator.Radon(image_domain=128, angles=180)
>>> sino = R.apply_to(img, queue=queue)
>>> backproj = R.T.apply_to(sino)

Operator algebra is also supported:

>>> G = R.T * R
>>> gram_img = G.apply_to(img, queue=queue)
property T: Radon
property adjoint: bool
property angles: Angles
apply_to(argument: Any, output: Any | None = None, queue: CommandQueue | None = None, return_event: bool = False, **kwargs: Any) → Array | tuple[Array, list[Event]][source]

Standard OpenCL-backed operator execution pipeline.

Passing output lets callers provide a preallocated device array and avoid allocating a new OpenCL array for the result.

property detectors: Detectors
property image_domain: ImageDomain
substitute_placeholder() → None[source]

Resolve extent placeholders to concrete numeric values.

OpenCL-backed operator helpers

OpenCL helpers for gratopy operators.

This module contains small building blocks for OpenCL-backed operators:

The classes in this module are intentionally lightweight. They are meant to support concrete operators such as gratopy.operator.Radon without pulling OpenCL-specific concerns into gratopy.operator.base.Operator.

class gratopy.operator.opencl.OpenCLKernelSpec(paths: tuple[str, ...], base_name: str, build_options: tuple[str, ...] = ())[source]

Bases: object

Description of an OpenCL kernel bundle.

This class describes the OpenCL source files used by an experimental gratopy operator as well as the kernel base name used for lookup.

Parameters

paths:

One or more kernel source files. Paths are normalized to absolute paths during initialization.

base_name:

Base name used for generated kernel lookup. With the current naming convention, forward and adjoint kernels are expected to follow the pattern

{base_name}_{dtype}_{out_order}{in_order}

and

{base_name}_ad_{dtype}_{out_order}{in_order},

respectively.

build_options:

Optional extra compiler flags forwarded to cl.Program.build().

Notes

The current operator implementation expects gratopy-style template placeholders in the source files, most notably \my_variable_type, \order1, and \order2.

The signature property is content-based: it depends on the file contents, their order, the base name, and the build options. This is used as part of the OpenCL program cache key so that changing a kernel file on disk automatically leads to recompilation when needed.

Examples Use the default shipped Radon kernels implicitly via gratopy.operator.projection.Radon, or provide an explicit custom kernel bundle:

>>> from gratopy.operator import OpenCLKernelSpec
>>> spec = OpenCLKernelSpec.from_path("scratch/my_radon.cl", base_name="radon")

The resulting spec can then be passed to a custom or built-in operator.

base_name: str
build_options: tuple[str, ...] = ()
classmethod from_path(path: str | Path, base_name: str, build_options: tuple[str, ...] = ()) → OpenCLKernelSpec[source]

Create a kernel spec from a single source file.

classmethod from_paths(paths: list[str | Path] | tuple[str | Path, ...], base_name: str, build_options: tuple[str, ...] = ()) → OpenCLKernelSpec[source]

Create a kernel spec from multiple source files.

paths: tuple[str, ...]
read_sources() → tuple[str, ...][source]

Return the source code of all configured kernel files.

property signature: str

Return a content-based signature used for program cacheing.

gratopy.operator.opencl.invalidate_kernel_cache() → None[source]

Clear the internal operator OpenCL program cache.

class gratopy.operator.opencl._OpenCLOperator(*, name: str | None = None, kernel_spec: OpenCLKernelSpec | None = None, **operator_kwargs: Any)[source]

Bases: Operator

Internal base class for OpenCL-backed operators.

This class intentionally only implements shared execution helpers. Concrete subclasses are still responsible for geometry handling and for deciding which kernels to load.

static _allocate_output(queue: CommandQueue, shape: tuple[int, ...], dtype: DTypeLike, order: Literal['C', 'F'] = 'F', allocator: Any = None) → Array[source]

Allocate an output array on the device.

static _coerce_argument(argument: ArrayLike | Array, queue: CommandQueue) → Array[source]

Coerce a NumPy/array-like input to a device array.

_default_kernel_spec() → OpenCLKernelSpec[source]

Return the default kernel spec for this operator.

Concrete subclasses are expected to override this if they want to rely on _OpenCLOperator construction directly.

_expected_output_shape(argument_shape: tuple[int, ...]) → tuple[int, ...][source]

Return the expected output shape for a given input shape.

_get_kernel(output: Array, argument: Array, queue: CommandQueue) → Kernel[source]

Return the kernel used for standard OpenCL operator execution.

_global_shape(output: Array, argument: Array) → tuple[int, ...][source]

Return the global shape used for kernel execution.

_infer_queue(argument: ArrayLike | Array | None = None, output: Array | None = None, queue: CommandQueue | None = None) → CommandQueue[source]

Infer the queue to use for computation.

_kernel_arguments(output: Array, argument: Array, queue: CommandQueue) → tuple[Any, ...][source]

Return additional kernel arguments beyond output/input buffers.

_validate_argument(argument: Array) → None[source]

Validate an input argument before kernel execution.

_validate_output(output: Array, argument: Array) → Array[source]

Validate and normalize an output array before kernel execution.

apply_to(argument: Any, output: Any | None = None, queue: CommandQueue | None = None, return_event: bool = False, **kwargs: Any) → Array | tuple[Array, list[Event]][source]

Standard OpenCL-backed operator execution pipeline.

Passing output lets callers provide a preallocated device array and avoid allocating a new OpenCL array for the result.

Utility geometry classes

The experimental operator API uses a small collection of helper classes for image domains, detector geometry, angular sampling, and extent placeholders.

Miscellaneous utility functions for gratopy.

class gratopy.utilities.Angles(angles: ArrayLike, weights: ArrayLike, half_circle: bool = False)[source]

Angular sampling together with associated weights.

An Angles object stores the projection angles in radians and the corresponding angle weights used, for instance, in the adjoint operator. It serves as the main explicit representation of angular geometry in the experimental operator API.

Parameters

angles:

One-dimensional array-like object containing angles in radians.

weights:

One-dimensional array-like object containing the weights associated with the given angles.

half_circle:

Whether the angular sampling should be interpreted with half-circle (\(\pi\)) periodicity instead of full-circle (\(2\pi\)) periodicity. For constructors that generate a complete equispaced sampling, this also determines whether angles are generated in \([0, \pi)\) or \([0, 2\pi)\). For limited-angle or explicit samplings, it is stored as metadata for downstream code.

Notes

The static constructors on this class provide a few common ways of constructing angular samplings:

static from_list(angles: ArrayLike, half_circle: bool = False) → Angles[source]

Automatically compute natural weights for a given list of angles.

Parameters:

  • angles – List of angles.

  • half_circle – If False (the default), angles are in \([0, 2\pi)\). Otherwise, angles are in \([0, \pi)\).

Returns:

Angles object with unitary weights.

static intervals(number_list: list[int], start_list: list[float], end_list: list[float], half_circle: bool = False) → Angles[source]

Create angles uniformly distributed across multiple intervals.

Parameters:

  • number_list – List of numbers of angles for each interval.

  • start_list – List of start points for each interval.

  • end_list – List of end points for each interval.

  • half_circle – Angular periodicity metadata. If True, the sampling is marked as using \(\pi\) periodicity; otherwise, it is marked as using \(2\pi\) periodicity.

Returns:

Angles object with angles and weights across all intervals.

static sparse(number: int, half_circle: bool = False) → Angles[source]

Store a list of angles with all weights set to 1.

Parameters:

  • number – Number of angles.

  • half_circle – If False (the default), angles are in \([0, 2\pi)\). Otherwise, angles are in \([0, \pi)\).

Returns:

Angles object with unitary weights.

static uniform(number: int, half_circle: bool = False) → Angles[source]

Generate uniformly distributed and weighted angles.

Parameters:

  • number – Number of angles.

  • half_circle – If False (the default), angles are in \([0, 2\pi)\). Otherwise, angles are in \([0, \pi)\).

Returns:

Angles object with uniform angles and weights.

static uniform_interval(start: float, end: float, number: int, half_circle: bool = False) → Angles[source]

Generate uniformly distributed and weighted angles in a specified interval.

Parameters:

  • start – Start of the interval.

  • end – End of the interval.

  • number – Number of angles.

  • half_circle – Angular periodicity metadata. If True, the sampling is marked as using \(\pi\) periodicity; otherwise, it is marked as using \(2\pi\) periodicity.

Returns:

Angles object with uniform angles and weights.

class gratopy.utilities.Detectors(number: int, extent: float | ExtentPlaceholder = ExtentPlaceholder.FULL, center: float = 0.0, reversed: bool | None = None)[source]

Detector discretization and physical placement.

This class describes the detector line used by an operator:

  • number is the number of detector pixels,

  • extent is the physical detector width or an extent placeholder,

  • center shifts the detector along its detector axis,

  • reversed flips the detector orientation.

Parameters

number:

Number of detector pixels. Negative values are accepted as a shorthand for reversed=True with abs(number) detector pixels.

extent:

Physical detector width or an ExtentPlaceholder.

center:

Physical shift of the detector center along the detector axis.

reversed:

Whether the detector orientation is reversed. If omitted, the sign of number is used to infer the orientation.

Notes

In the experimental Radon operator, ExtentPlaceholder.FULL and ExtentPlaceholder.VALID can be used here to resolve the detector extent from a fixed image extent. Passing placeholders for both detector and image extents at once is unsupported and raises NotImplementedError.

center: float = 0.0
extent: float | ExtentPlaceholder = 'full'
number: int
reversed: bool = False
class gratopy.utilities.ExtentPlaceholder(*values)[source]

Placeholder values for image and detector extents.

These placeholders are primarily intended for the experimental operator API, where image and detector extents can be specified independently. They express geometric intent rather than an immediate numerical value.

The placeholder mechanism in the experimental operator API is still evolving. The experimental Radon operator can resolve one placeholder extent at a time: either the image extent from a fixed detector extent, or the detector extent from a fixed image extent. Passing placeholders for both extents at once is unsupported and raises NotImplementedError.

FULL = 'full'

From the perspective of the detector: smallest detector extent such that each ray hitting the image domain also hits the detector. From the perspective of the image domain: largest image extent such that each ray hitting the image domain also hits the detector.

VALID = 'valid'

From the perspective of the detector: the largest extent such that each ray hitting the detector also hits the image domain. From the perspective of the image domain: the smallest extent such that each ray hitting the detector also hits the image domain.

class gratopy.utilities.ImageDomain(size: int | tuple[int, int], extent: float | ExtentPlaceholder = 2.0, center: tuple[float, float] = (0.0, 0.0))[source]

Image grid and physical image extent.

This class bundles the discrete image shape together with the physical extent and center shift used by an operator.

Parameters

size:

Image grid size. An integer is interpreted as a square domain of shape (N, N); a tuple specifies (Nx, Ny) directly.

extent:

Physical extent of the image domain or an ExtentPlaceholder. In gratopy’s conventions, this corresponds to the longer side length of the rectangular image domain.

center:

Physical shift of the image center relative to the rotation center.

Notes

In the experimental Radon operator, ExtentPlaceholder.FULL and ExtentPlaceholder.VALID can be used here to resolve the image extent from a fixed detector extent. Passing placeholders for both image and detector extents at once is unsupported and raises NotImplementedError.

center: tuple[float, float] = (0.0, 0.0)
extent: float | ExtentPlaceholder = 2.0
size: tuple[int, int]