# Function reference

## 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, geometry, img_shape, angles, n_detectors=None, angle_weights=None, detector_width=2.0, image_width=None, R=None, RE=None, detector_shift=0.0, midpoint_shift=[0.0, 0.0], reverse_detector=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 (int) – 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 (list[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:
create_sparse_matrix(dtype=dtype('float32'), order='F')[source]

Creates a sparse matrix representation of the associated forward projection.

Parameters:
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.

show_geometry(angle, figure=None, axes=None, show=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:
Returns:

Figure and axes in which the geometry visualization is plotted.

Return type:

## 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, projectionsetting, sino=None, wait_for=[])[source]

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

Parameters:
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, projectionsetting, img=None, wait_for=[])[source]

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

Parameters:
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, projectionsetting, number_iterations=100, w=1)[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:
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

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

Parameters:
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, projectionsetting, mu, number_iterations=1000, slice_thickness=1.0, stepsize_weighting=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 – Regularization parameter, the smaller the stronger the applied regularization.

• number_iterations (float, 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, number_iterations=50, dtype='float32', allocator=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:
Returns:

An estimate of the spectral norm for the projection operator.

Return type:

float

gratopy.weight_sinogram(sino, projectionsetting, sino_out=None, divide=False, wait_for=[])[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, N, modified=True, E=None, ret_E=False, dtype='double', allocator=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:

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.

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

Parameters:
Returns:

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

Return type:

pyopencl.Event

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

Parameters:
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, img_shape, angles, angle_weights, n_detectors=None, detector_width=2.0, image_width=2.0, midpoint_shift=[0, 0], detector_shift=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 (list[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, img, projectionsetting, wait_for=[])[source]

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

Parameters:
Returns:

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

Return type:

pyopencl.Event

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

Parameters:
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, img_shape, angles, detector_width, source_detector_dist, source_origin_dist, angle_weights, n_detectors=None, detector_shift=0.0, image_width=None, midpoint_shift=[0, 0], reverse_detector=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 (list[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()[source]

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

Returns:

The toolbox’s CL code.

Return type:

str