Basic Routines

Fast approximate discrete Radon transform for NumPy arrays.

The adrt package provides a fast implementation of an approximate discrete Radon transform and several related routines such as backprojection, inverses, and utilities.

This main module contains the basic functions of the package, designed to be simple to use. The module adrt.core contains lower-level routines which provide more control over the computation. Other utility routines can be found in adrt.utils.

Forward Transforms

These functions implement the fundamental operations of this package: the forward approximate discrete Radon transform (ADRT) and the associated backprojection operation.

adrt.adrt(a, /)

The Approximate Discrete Radon Transform (ADRT).

This is the fundamental routine of this package, computing the ADRT of the provided array. The array a must store square input images with sizes a power of two. The input may optionally include an additional leading batch dimension, so an array of either two or three dimensions.

If padding is needed for the input array, consider numpy.pad().

The returned array will have the shape of an ADRT output of size N. The output is divided into four quadrants, each one less than twice as tall as the input. The taller height dimension represents ADRT offsets, while the final dimension has the same size as the input and represents the ADRT angles.

For more information on the construction of the quadrants and the contents of this array see: Forward Transform.

Parameters:

a (numpy.ndarray of float) – Array for which the ADRT should be computed. This should be a square image with side length a power of two, and optionally a leading batch dimension.

Returns:

The ADRT of the provided data. For input images of size N, each member of the batch will have shape (4, 2*N-1, N).

Return type:

numpy.ndarray of float

Notes

The transform implemented here is an approximation to the Radon transform and approximates the sums along lines with carefully chosen angles. Each quadrant slices along a range of \(\pi/4\) radians. For a detailed description of the algorithm see Forward Transform and refer to the source papers [1] [2].

adrt.bdrt(a, /)

Backprojection operator for the ADRT.

The transform implemented in adrt() is a linear operation. This function computes a generalized transpose of this operation, producing an output with the same shape as its input.

To retrieve the entries of the transpose of the ADRT in the proper order, apply adrt.utils.truncate() to the array produced by this function. The truncation operation will remove the extended entries and rotate each quadrant into the same orientation as the original, square shape ADRT input. The quadrants can then be combined, if desired, potentially by numpy.mean().

Parameters:

a (numpy.ndarray of float) – An ADRT output array to backproject.

Returns:

Backprojection of a with the same shape.

Return type:

numpy.ndarray of float

Notes

For more details on the backprojection implemented here see the source paper [2].

Examples

As discussed above, this function can be used to compute the transpose of the operator applied by adrt() as follows:

def adrt_tranpose(a):
    return adrt.utils.truncate(adrt.bdrt(a)).mean(axis=-3)

Inverse Transforms

We provide two possible inverses to adrt(). The first is an exact (although ill-conditioned) inverse, while the second implements an approximate inverse by the full multigrid method.

adrt.iadrt(a, /)

An exact inverse to the ADRT.

Computes an exact inverse to the ADRT, but only works for exact ADRT outputs. The array a may have an optional batch dimension with the shape of an ADRT output.

The returned array has the same shape as a, but each quadrant should have only zeros below the square at the top of the array. However, this may not be the case due to conditioning or imprecision in the calculations performed by this routine.

The upper square of each quadrant can be extracted and rotated using adrt.utils.truncate().

Parameters:

a (numpy.ndarray of float) – An ADRT output for which to compute the inverse.

Returns:

The computed inverse with the same shape as a.

Return type:

numpy.ndarray of float

Warning

This inverse is ill-conditioned and will be exact only if a is an exact output of the forward ADRT and the floating point type provides sufficient precision. In other cases this inverse is not appropriate.

For an alternative, see the Iterative Inverse example.

Notes

For details of the algorithm see Inverse Transforms or the source paper [3].

adrt.iadrt_fmg(a, /, *, max_iters=None)

Approximate inverse to the ADRT by the full multigrid method.

Estimated inverses are computed and iteratively refined until the norm of the residual error fails to decrease from one iteration to the next. This iteration can also be terminated early if max_iters is specified.

Particularly with a limited iteration count this inverse can be relatively quick to compute, but may not achieve the best possible precision. If your output array a is exact and floating-point precision is sufficient you may consider iadrt() for an exact inverse. Otherwise, for an inverse that may perform more reliably for certain inputs, consider the iadrt_cg recipe proposed in the Iterative Inverse example.

See adrt.core.iadrt_fmg_iter() and adrt.core.iadrt_fmg_step() for more information on the iterative refinement applied internally and the full multigrid method.

Parameters:

• a (numpy.ndarray of float) – The array for which the inverse is to be computed. This array must have the shape of an ADRT output. This array must not have a batch dimension.

• max_iters (int, optional) – If None (default), the number of internal iterations is unbounded. The computation will terminate only when the norm of the residual error fails to decrease. Otherwise, this must be an integer argument at least 1, and provides an upper bound on the number of iterations performed internally.

Returns:

A square inverse computed by the full multigrid method with the lowest residual error observed so long as the decrease was monotonic, or the last computation output if iteration was terminated by max_iters.

Return type:

numpy.ndarray of float

Notes

Unlike iadrt() the output of this routine will be square—the same size as the original input to adrt() would have been.

This function does not support batch dimensions.

[1]

Martin L. Brady, Discrete Radon transform has an exact, fast inverse and generalizes to operations other than sums along lines, Proceedings of the National Academy of Sciences, 103 doi:10.1137/S0097539793256673.

[2] (1,2)

William H. Press, A Fast Discrete Approximation Algorithm for the Radon Transform Related Databases, SIAM Journal on Computing, 27 doi:10.1073/pnas.0609228103.

[3]

Donsub Rim, Exact and fast inversion of the approximate discrete Radon transform from partial data, Applied Mathematics Letters, 102 doi:10.1016/j.aml.2019.106159.