CIL Plugins#

CCPi Regularisation#

This plugin allows the use of regularisation functions from the CCPi Regularisation toolkit (10.1016/j.softx.2019.04.003, a set of CPU/GPU optimised regularisation modules for iterative image reconstruction and other image processing tasks.

Total variation#

class cil.plugins.ccpi_regularisation.functions.FGP_TV(alpha=1, max_iteration=100, tolerance=0, isotropic=True, nonnegativity=None, device='cpu', strong_convexity_constant=0)[source]#

Fast Gradient Projection Total Variation (FGP_TV)

The FGP_TV computes the proximal operator of the Total variation regulariser

\[\mathrm{prox}_{\tau (\alpha TV)}(x) = \underset{z}{\mathrm{argmin}} \,\alpha\,\mathrm{TV}(z) + \frac{1}{2}\|z - x\|^{2} .\]

The algorithm used for the proximal operator of TV is the Fast Gradient Projection algorithm applied to the _dual problem_ of the above problem, see [1], [2].

Note

In CIL Version 24.1.0 we change the default value of nonnegativity to False. This means non-negativity is not enforced by default.

Parameters:

alpha (Number (positive), default = 1.0 .) – Total variation regularisation parameter.
max_iteration (int. Default = 100 .) – Maximum number of iterations for the Fast Gradient Projection algorithm.
isotropic (boolean. Default = True .) –
Isotropic or Anisotropic definition of the Total variation regulariser.

\[|x|_{2} = \sqrt{x_{1}^{2} + x_{2}^{2}},\, (\mbox{isotropic})\]

\[|x|_{1} = |x_{1}| + |x_{2}|\, (\mbox{anisotropic})\]
nonnegativity (boolean. Default = False .) – Non-negativity constraint for the solution of the FGP algorithm.
tolerance (float, Default = 0 .) –
Stopping criterion for the FGP algorithm.

\[\|x^{k+1} - x^{k}\|_{2} < \mathrm{tolerance}\]
device (str, Default = ‘cpu’ .) – FGP_TV algorithm runs on cpu or gpu.
strong_convexity_constant (float, default = 0) –
A strongly convex term weighted by the strong_convexity_constant (\(\gamma\)) parameter is added to the Total variation. Now the TotalVariation function is \(\gamma\) - strongly convex and the proximal operator is

\[\underset{u}{\mathrm{argmin}} \frac{1}{2\tau}\|u - b\|^{2} + \mathrm{TV}(u) + \frac{\gamma}{2}\|u\|^{2} \Leftrightarrow\]

\[\underset{u}{\mathrm{argmin}} \frac{1}{2\frac{\tau}{1+\gamma\tau}}\|u - \frac{b}{1+\gamma\tau}\|^{2} + \mathrm{TV}(u)\]

Examples

\[\underset{u\qeq0}{\mathrm{argmin}} \frac{1}{2}\|u - b\|^{2} + \alpha TV(u)\]

>>> G = alpha * FGP_TV(max_iteration=100, device='gpu')
>>> sol = G.proximal(b)

Note

The FGP_TV regularisation does not incorparate information on the ImageGeometry, i.e., pixel/voxel size. Therefore a rescaled parameter should be used to match the same solution computed using TotalVariation.

>>> G1 = (alpha/ig.voxel_size_x) * FGP_TV(max_iteration=100, device='gpu')
>>> G2 = alpha * TotalVariation(max_iteration=100, lower=0.)

Other regularisation functions#

class cil.plugins.ccpi_regularisation.functions.TGV(alpha=1, gamma=1, max_iteration=100, tolerance=0, device='cpu', **kwargs)[source]#

__init__(alpha=1, gamma=1, max_iteration=100, tolerance=0, device='cpu', **kwargs)[source]#

Creator of Total Generalised Variation Function

Parameters:

alpha (number, default 1) – regularisation parameter
gamma (number, default 1, can range between 1 and 2) – ratio of TGV terms
max_iteration (integer, default 100) – max number of sub iterations. The algorithm will iterate up to this number of iteration or up to when the tolerance has been reached
tolerance (float, default 0) – minimum difference between previous iteration of the algorithm that determines the stop of the iteration earlier than max_iteration. If set to 0 only the max_iteration will be used as stop criterion.
device (string, default 'cpu', can be 'gpu' if GPU is installed) – determines if the code runs on CPU or GPU

__call__(x)[source]#: Call self as a function.

convex_conjugate(x)[source]#

Evaluation of the function F* at x, where F* is the convex conjugate of function F,

\[F^{*}(x^{*}) = \underset{x}{\sup} \langle x^{*}, x \rangle - F(x)\]

Parameters:: x (DataContainer)
Return type:: The value of the convex conjugate of the function at x.

__rmul__(scalar)[source]#

Define the multiplication with a scalar

this changes the regularisation parameter in the plugin

class cil.plugins.ccpi_regularisation.functions.FGP_dTV(reference, alpha=1, max_iteration=100, tolerance=0, eta=0.01, isotropic=True, nonnegativity=True, device='cpu')[source]#

Creator of FGP_dTV Function

Parameters:

reference (ImageData) – reference image
alpha (number, default 1) – regularisation parameter
max_iteration (integer, default 100) – max number of sub iterations. The algorithm will iterate up to this number of iteration or up to when the tolerance has been reached
tolerance (float, default 0) – minimum difference between previous iteration of the algorithm that determines the stop of the iteration earlier than max_iteration. If set to 0 only the max_iteration will be used as stop criterion.
eta (number, default 0.01) – smoothing constant to calculate gradient of the reference
isotropic (boolean, default True, can range between 1 and 2) – Whether it uses L2 (isotropic) or L1 (anisotropic) norm
nonnegativity (boolean, default True) – Whether to add the non-negativity constraint
device (string, default 'cpu', can be 'gpu' if GPU is installed) – determines if the code runs on CPU or GPU

__init__(reference, alpha=1, max_iteration=100, tolerance=0, eta=0.01, isotropic=True, nonnegativity=True, device='cpu')[source]#

__call__(x)[source]#: Call self as a function.

convex_conjugate(x)[source]#

Evaluation of the function F* at x, where F* is the convex conjugate of function F,

\[F^{*}(x^{*}) = \underset{x}{\sup} \langle x^{*}, x \rangle - F(x)\]

Parameters:: x (DataContainer)
Return type:: The value of the convex conjugate of the function at x.

__rmul__(scalar)[source]#

Define the multiplication with a scalar

this changes the regularisation parameter in the plugin

class cil.plugins.ccpi_regularisation.functions.TNV(alpha=1, max_iteration=100, tolerance=0)[source]#

__init__(alpha=1, max_iteration=100, tolerance=0)[source]#

Creator of TNV Function

Parameters:

alpha (number, default 1) – regularisation parameter
max_iteration (integer, default 100) – max number of sub iterations. The algorithm will iterate up to this number of iteration or up to when the tolerance has been reached
tolerance (float, default 0) – minimum difference between previous iteration of the algorithm that determines the stop of the iteration earlier than max_iteration. If set to 0 only the max_iteration will be used as stop criterion.

__call__(x)[source]#: Call self as a function.

convex_conjugate(x)[source]#

Evaluation of the function F* at x, where F* is the convex conjugate of function F,

\[F^{*}(x^{*}) = \underset{x}{\sup} \langle x^{*}, x \rangle - F(x)\]

Parameters:: x (DataContainer)
Return type:: The value of the convex conjugate of the function at x.

__rmul__(scalar)[source]#

Define the multiplication with a scalar

this changes the regularisation parameter in the plugin

check_input(input)[source]#: TNV requires 2D+channel data with the first dimension as the channel dimension

TomoPhantom#

This plugin allows the use of part of TomoPhantom (10.1016/j.softx.2018.05.003, a toolbox written in C language to generate customisable 2D-4D phantoms (with a temporal capability).

cil.plugins.TomoPhantom.get_ImageData(num_model, geometry)[source]#

Returns an ImageData relative to geometry with the model num_model from tomophantom

Parameters:

num_model (int) – model number
geometry (ImageGeometry) – geometrical info that describes the phantom

Example usage:

ndim = 2
N=128
angles = np.linspace(0, 360, 50, True, dtype=np.float32)
offset = 0.4
channels = 3

if ndim == 2:
    ag = AcquisitionGeometry.create_Cone2D((offset,-100), (offset,100))
    ag.set_panel(N)

else:
    ag = AcquisitionGeometry.create_Cone3D((offset,-100, 0), (offset,100,0))
    ag.set_panel((N,N-2))

ag.set_channels(channels)
ag.set_angles(angles, angle_unit=AngleUnit.DEGREE)

ig = ag.get_ImageGeometry()
num_model = 1
phantom = TomoPhantom.get_ImageData(num_model=num_model, geometry=ig)

TIGRE#

This plugin allows the use of TIGRE (10.1088/2057-1976/2/5/055010 for forward and back projections and filter back projection reconstruction.

FBP#

This reconstructs with FBP for parallel-beam data, and with FDK weights for cone-beam data

class cil.plugins.tigre.FBP(image_geometry=None, acquisition_geometry=None, **kwargs)[source]#

FBP Filtered Back Projection is a reconstructor for 2D and 3D parallel and cone-beam geometries. It is able to back-project circular trajectories with 2 PI angular range and equally spaced angular steps.

This uses the ram-lak filter This is provided for simple and offset parallel-beam geometries only

acquisition_geometryAcquisitionGeometry: A description of the acquisition data
image_geometryImageGeometry, default used if None: A description of the area/volume to reconstruct

Example

>>> from cil.plugins.tigre import FBP
>>> fbp = FBP(image_geometry, data.geometry)
>>> fbp.set_input(data)
>>> reconstruction = fbp.get_output()

get_output(out=None)#

Runs the configured processor and returns the processed data

Parameters:: out (DataContainer, optional) – Fills the referenced DataContainer with the processed data
Returns:: The processed data
Return type:: DataContainer

set_input(dataset)#

Set the input data to the processor

Parameters:: input (DataContainer) – The input DataContainer

Projection Operator#

class cil.plugins.tigre.ProjectionOperator(image_geometry=None, acquisition_geometry=None, direct_method='interpolated', adjoint_weights='matched', **kwargs)[source]#

ProjectionOperator configures and calls TIGRE Projectors for your dataset.

Please refer to the TIGRE documentation for further descriptions CERN/TIGRE https://iopscience.iop.org/article/10.1088/2057-1976/2/5/055010

Parameters:

image_geometry (ImageGeometry, default used if None) – A description of the area/volume to reconstruct
acquisition_geometry (AcquisitionGeometry, BlockGeometry) – A description of the acquisition data. If passed a BlockGeometry it will return a BlockOperator.
direct_method (str, default 'interpolated') – The method used by the forward projector, ‘Siddon’ for ray-voxel intersection, ‘interpolated’ for interpolated projection
adjoint_weights (str, default 'matched') – The weighting method used by the cone-beam backward projector, ‘matched’ for weights to approximately match the ‘interpolated’ forward projector, ‘FDK’ for FDK weights

Example

>>> from cil.plugins.tigre import ProjectionOperator
>>> PO = ProjectionOperator(image.geometry, data.geometry)
>>> forward_projection = PO.direct(image)
>>> backward_projection = PO.adjoint(data)

ASTRA#

This plugin allows the use of ASTRA-toolbox (10.1364/OE.24.025129) for forward and back projections and filter back projection reconstruction.

FBP#

This reconstructs with FBP for parallel-beam data, and with FDK weights for cone-beam data

class cil.plugins.astra.FBP(image_geometry=None, acquisition_geometry=None, device='gpu')[source]#

FBP configures and calls an appropriate ASTRA FBP or FDK algorithm for your dataset.

The best results will be on data with circular trajectories of a 2PI angular range and equally spaced small angular steps.

Parameters:

image_geometry (ImageGeometry, default used if None) – A description of the area/volume to reconstruct
acquisition_geometry (AcquisitionGeometry) – A description of the acquisition data
device (string, default='gpu') – ‘gpu’ will run on a compatible CUDA capable device using the ASTRA FDK_CUDA algorithm ‘cpu’ will run on CPU using the ASTRA FBP algorithm - see Notes for limitations

Example

>>> from cil.plugins.astra import FBP
>>> fbp = FBP(image_geometry, data.geometry)
>>> fbp.set_input(data)
>>> reconstruction = fbp.get_output()

Notes

A CPU version is provided for simple 2D parallel-beam geometries only, any offsets and rotations in the acquisition geometry will be ignored.

This uses the ram-lak filter only.

set_input(dataset)[source]#

Set the input data to the processor

Parameters:: input (DataContainer) – The input DataContainer

get_output(out=None)[source]#

Runs the configured processor and returns the processed data

Parameters:: out (DataContainer, optional) – Fills the referenced DataContainer with the processed data
Returns:: The processed data
Return type:: DataContainer

Projection Operator#

class cil.plugins.astra.ProjectionOperator(image_geometry=None, acquisition_geometry=None, device='gpu', **kwargs)[source]#

ProjectionOperator configures and calls appropriate ASTRA Projectors for your dataset.

Parameters:

image_geometry (ImageGeometry, default used if None) – A description of the area/volume to reconstruct
acquisition_geometry (AcquisitionGeometry, BlockGeometry) – A description of the acquisition data. If passed a BlockGeometry it will return a BlockOperator.
device (string, default='gpu') – ‘gpu’ will run on a compatible CUDA capable device using the ASTRA 3D CUDA Projectors, ‘cpu’ will run on CPU using the ASTRA 2D CPU Projectors

Example

>>> from cil.plugins.astra import ProjectionOperator
>>> PO = ProjectionOperator(image.geometry, data.geometry)
>>> forward_projection = PO.direct(image)
>>> backward_projection = PO.adjoint(data)

Notes

For multichannel data the ProjectionOperator will broadcast across all channels.

Return Home

Table of Contents

CIL Plugins#

CCPi Regularisation#

Total variation#

Other regularisation functions#

TomoPhantom#

TIGRE#

FBP#

Projection Operator#

ASTRA#

FBP#

Projection Operator#