Framework¶
The goal of the CCPi Framework is to allow the user to simply create iterative reconstruction methods which go beyond the standard filter back projection technique and which better suit the data characteristics. The framework comprises:
ccpi.framework
module which allows to simply translate real world CT systems into software.ccpi.optimisation
module allows the user to quickly create iterative methods to reconstruct acquisition data applying different types of regularisation, which better suit the data characteristics.ccpi.io
module which provides a number of loaders for real CT machines, e.g. Nikon. It also provides reader and writer to save to NeXuS file format.
CT Geometry¶
Please refer to this notebook on the CILDemos repository for full description.
In conventional CT systems, an object is placed between a source emitting Xrays and a detector array measuring the Xray transmission images of the incident Xrays. Typically, either the object is placed on a rotating sample stage and rotates with respect to the sourcedetector assembly, or the sourcedetector gantry rotates with respect to the stationary object. This arrangement results in socalled circular scanning trajectory. Depending on source and detector types, there are three conventional data acquisition geometries:
 parallel geometry (2D or 3D),
 fanbeam geometry, and
 conebeam geometry.
Parallel geometry¶
Parallel beams of Xrays are emitted onto 1D (single pixel row) or 2D detector array. This geometry is common for synchrotron sources. 2D parrallel geometry is illustrated below.
Fanbeam geometry¶
 A single pointlike Xray source emits a cone beam onto 1D detector pixel row. Conebeam is typically
 collimated to imaging field of view. Collimation allows greatly reduce amount of scatter radiation reaching the detector. Fanbeam geometry is used when scattering has significant influence on image quality or singleslice reconstruction is sufficient.
Conebeam geometry¶
A single pointlike Xray source emits a cone beam onto 2D detector array. Conebeam geometry is mainly used in labbased CT instruments. Depending on where the sample is placed between the source and the detector one can achieve a different magnification factor \(F\):
where \(r_1\) and \(r_2\) are the distance from the source to the center of the sample and the distance from the center of the sample to the detector, respectively.
AcquisitonGeometry and AcquisitionData¶
In the Framework, we implemented AcquisitionGeometry
class to hold acquisition parameters and
ImageGeometry
to hold geometry of a reconstructed volume. Corresponding data arrays are wrapped
as AcquisitionData
and ImageData
classes, respectively.
The simplest (of course from image processing point of view, not from physical implementation) geometry is the parallel geometry. Geometrical parameters for parallel geometry are depicted below:
In the Framework, we define AcquisitionGeometry
as follows.
# imports
from ccpi.framework import AcquisitionGeometry
import numpy as np
# acquisition angles
n_angles = 90
angles = np.linspace(0, np.pi, n_angles, dtype=np.float32)
# number of pixels in detector row
N = 256
# pixel size
pixel_size_h = 1
# # create AcquisitionGeometry
ag_par = AcquisitionGeometry(geom_type='parallel',
dimension='2D',
angles=angles,
pixel_num_h=N,
pixel_size_h=pixel_size_h)
AcquisitionGeometry
contains only metadata, the actual data is wrapped in AcquisitionData
class. AcquisitionGeometry
class also holds information about arrangement of the actual
acquisition data array. We use attribute dimension_labels
to label axis. The expected dimension labels are shown below.
The default order of dimensions for AcquisitionData
is [angle, horizontal]
, meaning that
the number of elements along 0 and 1 axes in the acquisition data array is expected to be n_angles
and N
, respectively.
To have consistent AcquisitionData
and AcquisitionGeometry
, we recommend to allocate
AcquisitionData
using allocate
method of the AcquisitionGeometry
instance:
# allocate AcquisitionData
ad_par = ag_par.allocate()
ImageGeometry and ImageData¶
To store reconstruction results, we implemented two classes: ImageGeometry
and ImageData
classes.
Similar to AcquisitionData
and AcquisitionGeometry
, we first define 2D ImageGeometry
and then allocate ImageData
.
# imports
from ccpi.framework import ImageData, ImageGeometry
# define 2D ImageGeometry
# given AcquisitionGeometry ag_par, default parameters for corresponding ImageData
ig_par = ImageGeometry(voxel_num_y=ag_par.pixel_num_h,
voxel_size_x=ag_par.pixel_size_h,
voxel_num_x=ag_par.pixel_num_h,
voxel_size_y=ag_par.pixel_size_h)
# allocate ImageData filled with 0 values with the specific geometry
im_data1 = ig_par.allocate()
# allocate ImageData filled with random values with the specific geometry
im_data2 = ig_par.allocate('random', seed=5)
Fanbeam, conebeam and 3D (multislice) parallel geometry can be setup similar to 2D parallel geometry.
3D parallel geometry¶
3D parallel beam AcquisitionGeometry
and default ImageGeometry
parameters can be set up
as follows:
# setup 3D parallel beam AcquisitionGeometry
# physical pixel size
pixel_size_h = 1
ag_par_3d = AcquisitionGeometry(geom_type='parallel',
dimension='3D',
angles=angles,
pixel_num_h=N,
pixel_size_h=pixel_size_h,
pixel_num_v=N,
pixel_size_v=pixel_size_h)
# setup 3D parallel beam ImageGeometry
ig_par_3d = ImageGeometry(voxel_num_x=ag_par_3d.pixel_num_h,
voxel_size_x=ag_par_3d.pixel_size_h,
voxel_num_y=ag_par_3d.pixel_num_h,
voxel_size_y=ag_par_3d.pixel_size_h,
voxel_num_z=ag_par_3d.pixel_num_v,
voxel_size_z=ag_par_3d.pixel_size_v)
Fanbeam geometry¶
Fanbeam AcquisitionGeometry
and
default ImageGeometry
parameters can be set up as follows:
# setup fanbeam AcquisitionGeometry
# distance from source to center of rotation
dist_source_center = 200.0
# distance from center of rotation to detector
dist_center_detector = 300.0
# physical pixel size
pixel_size_h = 2
ag_fan = AcquisitionGeometry(geom_type='cone',
dimension='2D',
angles=angles,
pixel_num_h=N,
pixel_size_h=pixel_size_h,
dist_source_center=dist_source_center,
dist_center_detector=dist_center_detector)
# calculate geometrical magnification
mag = (ag_fan.dist_source_center + ag_fan.dist_center_detector) / ag_fan.dist_source_center
ig_fan = ImageGeometry(voxel_num_x=ag_fan.pixel_num_h,
voxel_size_x=ag_fan.pixel_size_h / mag,
voxel_num_y=ag_fan.pixel_num_h,
voxel_size_y=ag_fan.pixel_size_h / mag)

class
ccpi.framework.
ImageGeometry
(voxel_num_x=0, voxel_num_y=0, voxel_num_z=0, voxel_size_x=1, voxel_size_y=1, voxel_size_z=1, center_x=0, center_y=0, center_z=0, channels=1, **kwargs)[source]¶

class
ccpi.framework.
AcquisitionGeometry
(geom_type, dimension=None, angles=None, pixel_num_h=0, pixel_size_h=1, pixel_num_v=0, pixel_size_v=1, dist_source_center=None, dist_center_detector=None, channels=1, **kwargs)[source]¶ 
allocate
(value=0, dimension_labels=None, **kwargs)[source]¶ allocates an AcquisitionData according to the size expressed in the instance
Parameters:  value (number or string, default None allocates empty memory block) – accepts numbers to allocate an uniform array, or a string as ‘random’ or ‘random_int’ to create a random array or None.
 dimension_labels – labels for the dimension axis
 dtype (numpy type, default numpy.float32) – numerical type to allocate


class
ccpi.framework.
VectorGeometry
(length)[source]¶ Geometry describing VectorData to contain 1D array
DataContainer
and subclasses AcquisitionData
and ImageData
are
meant to contain data and metadata in AcquisitionGeometry
and
ImageGeometry
respectively.
DataContainer and subclasses¶
AcquisiionData
and ImageData
inherit from the same parent DataContainer
class,
therefore they largely behave the same way.
There are algebraic operations defined for both AcquisitionData
and ImageData
.
Following operations are defined:
 binary operations (between two DataContainers or scalar and DataContainer)
+
addition
subtraction/
division*
multiplication**
powermaximum
minimum
 inplace operations
+=
=
*=
**=
/=
 unary operations
abs
sqrt
sign
conjugate
 reductions
sum
norm
dot
product
AcquisitionData
and ImageData
provide a simple method to transpose the data and to
produce a subset of itself based on the axis we would like to have. This method is based on the label of
the axes of the data rather than the way it is stored. We think that the user should describe what she
wants and not bother with knowing the actual layout of the data in the memory.
# transpose data using subset method
data_transposed = data.subset(['horizontal_y', 'horizontal_x'])
# extract single row
data_profile = data_subset.subset(horizontal_y=100)

class
ccpi.framework.
DataContainer
(array, deep_copy=True, dimension_labels=None, **kwargs)[source]¶ Generic class to hold data
Data is currently held in a numpy arrays

__weakref__
¶ list of weak references to the object (if defined)

as_array
(dimensions=None)[source]¶ Returns the DataContainer as Numpy Array
Returns the pointer to the array if dimensions is not set. If dimensions is set, it first creates a new DataContainer with the subset and then it returns the pointer to the array

axpby
(a, b, y, out, dtype=<class 'numpy.float32'>, num_threads=2)[source]¶ performs axpby with cilacc C library
Does the operation .. math:: a*x+b*y and stores the result in out, where x is self
Parameters:  a (float) – scalar
 b (float) – scalar
 y – DataContainer
 out – DataContainer instance to store the result
 dtype (numpy type, optional, default numpy.float32) – data type of the DataContainers
 num_threads (int, optional, default 1/2 CPU of the system) – number of threads to run on

dot
(other, *args, **kwargs)[source]¶ return the inner product of 2 DataContainers viewed as vectors
applies to real and complex data. In such case the dot method returns
a.dot(b.conjugate())

dtype
¶ Returns the type of the data array

fill
(array, **dimension)[source]¶ fills the internal numpy array with the one provided
Parameters:  array (DataContainer, numpy array or number) – numpy array to copy into the DataContainer
 dimension – dictionary, optional

get_data_axes_order
(new_order=None)[source]¶ returns the axes label of self as a list
if new_order is None returns the labels of the axes as a sortedbykey list if new_order is a list of length number_of_dimensions, returns a list with the indices of the axes in new_order with respect to those in self.dimension_labels: i.e.
self.dimension_labels = {0:’horizontal’,1:’vertical’} new_order = [‘vertical’,’horizontal’] returns [1,0]

size
¶ Returns the number of elements of the DataContainer


class
ccpi.framework.
ImageData
(array=None, deep_copy=False, dimension_labels=None, **kwargs)[source]¶ DataContainer for holding 2D or 3D DataContainer

class
ccpi.framework.
AcquisitionData
(array=None, deep_copy=True, dimension_labels=None, **kwargs)[source]¶ DataContainer for holding 2D or 3D sinogram
Multi channel data¶
Both AcquisitionGeometry
, AcquisitionData
and ImageGeometry
, ImageData
can be defined for multichannel (spectral) CT data using channels
attribute.
# multichannel fanbeam geometry
ag_fan_mc = AcquisitionGeometry(geom_type='cone',
dimension='2D',
angles=angles,
pixel_num_h=N,
pixel_size_h=1,
dist_source_center=200,
dist_center_detector=300,
channels=10)
# define multichannel 2D ImageGeometry
ig3 = ImageGeometry(voxel_num_y=5,
voxel_num_x=4,
channels=2)
Block Framework¶
The block framework allows writing more advanced optimisation problems. Consider the typical Tikhonov regularisation:
where,
 \(A\) is the projection operator
 \(b\) is the acquired data
 \(u\) is the unknown image to be solved for
 \(\alpha\) is the regularisation parameter
 \(L\) is a regularisation operator
The first term measures the fidelity of the solution to the data. The second term meausures the fidelity to the prior knowledge we have imposed on the system, operator \(L\).
This can be rewritten equivalently in the block matrix form:
With the definitions:
 \(\tilde{A} = \binom{A}{\alpha L}\)
 \(\tilde{b} = \binom{b}{0}\)
this can now be recognised as a least squares problem which can be solved by any algorithm in the ccpi.optimisation
which can solve least squares problem, e.g. CGLS.
To be able to express our optimisation problems in the matrix form above, we developed the socalled,
Block Framework comprising 4 main actors: BlockGeometry
, BlockDataContainer
,
BlockFunction
and BlockOperator
.
A BlockDataContainer
can be instantiated from a number of DataContainer
and subclasses
represents a column vector of :code:`DataContainer`s.
bdc = BlockDataContainer(DataContainer0, DataContainer1)
. These
classes are required for it to work. They provide a base class that will
behave as normal DataContainer
.

class
ccpi.framework.
BlockDataContainer
(*args, **kwargs)[source]¶ Class to hold DataContainers as column vector
Provides basic algebra between BlockDataContainer’s, DataContainer’s and subclasses and Numbers
 algebra between `BlockDataContainer`s will be elementwise, only if the shape of the 2 `BlockDataContainer`s is the same, otherwise it will fail
 algebra between BlockDataContainer`s and `list or numpy array will work as long as the number of rows and element of the arrays match, indipendently on the fact that the BlockDataContainer could be nested
 algebra between BlockDataContainer and one DataContainer is possible. It will require that all the DataContainers in the block to be compatible with the DataContainer we want to algebra with. Should we require that the DataContainer is the same type? Like ImageData or AcquisitionData?
 algebra between BlockDataContainer and a Number is possible and it will be done with each element of the BlockDataContainer even if nested
A = [ [B,C] , D] A * 3 = [ 3 * [B,C] , 3* D] = [ [ 3*B, 3*C] , 3*D ]

__radd__
(other)[source]¶ Reverse addition
to make sure that this method is called rather than the __mul__ of a numpy array the class constant __array_priority__ must be set > 0 https://docs.scipy.org/doc/numpy1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__

__rdiv__
(other)[source]¶ Reverse division
to make sure that this method is called rather than the __mul__ of a numpy array the class constant __array_priority__ must be set > 0 https://docs.scipy.org/doc/numpy1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__

__rmul__
(other)[source]¶ Reverse multiplication
to make sure that this method is called rather than the __mul__ of a numpy array the class constant __array_priority__ must be set > 0 https://docs.scipy.org/doc/numpy1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__

__rpow__
(other)[source]¶ Reverse power
to make sure that this method is called rather than the __mul__ of a numpy array the class constant __array_priority__ must be set > 0 https://docs.scipy.org/doc/numpy1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__

__rsub__
(other)[source]¶ Reverse subtraction
to make sure that this method is called rather than the __mul__ of a numpy array the class constant __array_priority__ must be set > 0 https://docs.scipy.org/doc/numpy1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__

__rtruediv__
(other)[source]¶ Reverse truedivision
to make sure that this method is called rather than the __mul__ of a numpy array the class constant __array_priority__ must be set > 0 https://docs.scipy.org/doc/numpy1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__

__weakref__
¶ list of weak references to the object (if defined)

add
(other, *args, **kwargs)[source]¶ Algebra: add method of BlockDataContainer with number/DataContainer or BlockDataContainer
Param: other (number, DataContainer or subclasses or BlockDataContainer Param: out (optional): provides a placehold for the resul.

axpby
(a, b, y, out, dtype=<class 'numpy.float32'>, num_threads=2)[source]¶ performs axpby elementwise on the BlockDataContainer containers
Does the operation .. math:: a*x+b*y and stores the result in out, where x is self
Parameters:  a – scalar
 b – scalar
 y – compatible (Block)DataContainer
 out – (Block)DataContainer to store the result
 dtype – optional, data type of the DataContainers

binary_operations
(operation, other, *args, **kwargs)[source]¶ Algebra: generic method of algebric operation with BlockDataContainer with number/DataContainer or BlockDataContainer
Provides commutativity with DataContainer and subclasses, i.e. this class’s reverse algebric methods take precedence w.r.t. direct algebric methods of DataContainer and subclasses.
This method is not to be used directly

divide
(other, *args, **kwargs)[source]¶ Algebra: divide method of BlockDataContainer with number/DataContainer or BlockDataContainer
Param: other (number, DataContainer or subclasses or BlockDataContainer Param: out (optional): provides a placehold for the resul.

maximum
(other, *args, **kwargs)[source]¶ Algebra: power method of BlockDataContainer with number/DataContainer or BlockDataContainer
Param: other (number, DataContainer or subclasses or BlockDataContainer Param: out (optional): provides a placehold for the resul.

minimum
(other, *args, **kwargs)[source]¶ Algebra: power method of BlockDataContainer with number/DataContainer or BlockDataContainer
Param: other (number, DataContainer or subclasses or BlockDataContainer Param: out (optional): provides a placehold for the resul.

multiply
(other, *args, **kwargs)[source]¶ Algebra: multiply method of BlockDataContainer with number/DataContainer or BlockDataContainer
Param: other (number, DataContainer or subclasses or BlockDataContainer Param: out (optional): provides a placehold for the resul.

power
(other, *args, **kwargs)[source]¶ Algebra: power method of BlockDataContainer with number/DataContainer or BlockDataContainer
Param: other (number, DataContainer or subclasses or BlockDataContainer Param: out (optional): provides a placehold for the resul.
DataProcessor¶
A DataProcessor
takes as an input a DataContainer
or subclass and returns either
another DataContainer
or some number. The aim of this class is to simplify the writing of
processing pipelines.

class
ccpi.framework.
DataProcessor
(**attributes)[source]¶ Defines a generic DataContainer processor
accepts DataContainer as inputs and outputs DataContainer additional attributes can be defined with __setattr__

__weakref__
¶ list of weak references to the object (if defined)

Resizer¶
Quite often we need either crop or downsample data; the Resizer
provides a convenient way to
perform these operations for both ImageData
and AcquisitionData
.
# imports
from ccpi.processors import Resizer
# crop ImageData along 1st dimension
# initialise Resizer
resizer_crop = Resizer(binning = [1, 1], roi = [1, (20,180)])
# pass DataContainer
resizer_crop.input = data
data_cropped = resizer_crop.process()
# get new ImageGeometry
ig_data_cropped = data_cropped.geometry

class
ccpi.processors.
Resizer
(roi=1, binning=1)[source]¶ 
__init__
(roi=1, binning=1)[source]¶ Constructor
Input:
 roi regionofinterest to crop. If roi = 1 (default), then no crop.
Otherwise roi is given by a list with ndim elements, where each element is either 1 if no crop along this dimension or a tuple with beginning and end coodinates to crop to. Example:
to crop 4D array along 2nd dimension: roi = [1, 1, (100, 900), 1] binning number of pixels to bin (combine) along each dimension.
If binning = 1, then projections in original resolution are loaded. Otherwise, binning is given by a list with ndim integers. Example:
to rebin 3D array along 1st direction: binning = [1, 5, 1]

Calculation of Center of Rotation¶
In the ideal alignment of a CT instrument, orthogonal projection of an axis of rotation onto a
detector has to coincide with a vertical midline of the detector. This is barely feasible in practice
due to misalignment and/or kinematic errors in positioning of CT instrument components.
A slight offset of the center of rotation with respect to the theoretical position will contribute
to the loss of resolution; in more severe cases, it will cause severe artifacts in the reconstructed
volume (doubleborders). CenterOfRotationFinder
allows to estimate offset of center of rotation
from theoretical. In the current release CenterOfRotationFinder
supports only parallel geometry.
CenterOfRotationFinder
is based on Nghia Vo’s method.

class
ccpi.processors.
CenterOfRotationFinder
(smin=None, smax=None, srad=None, step=None, ratio=None, drop=None)[source]¶ Processor to find the center of rotation in a parallel beam experiment
This processor read in a AcquisitionDataSet and finds the center of rotation based on Nghia Vo’s method. https://doi.org/10.1364/OE.22.019078
Input: AcquisitionDataSet Set_slice: Slice index or ‘centre’ smin, smax : int, optional
Reference to the horizontal center of the sinogram. srad : float, optional
 Fine search radius.
 step : float, optional
 Step of fine searching.
 ratio : float, optional
 The ratio between the FOV of the camera and the size of object. It’s used to generate the mask.
 drop : int, optional
 Drop lines around vertical center of the mask.
Output: float. center of rotation in pixel coordinate

__init__
(smin=None, smax=None, srad=None, step=None, ratio=None, drop=None)[source]¶ Initialize self. See help(type(self)) for accurate signature.

static
_search_coarse
(sino, smin, smax, ratio, drop)[source]¶ Coarse search for finding the rotation center.

static
_search_fine
(sino, srad, step, init_cen, ratio, drop)[source]¶ Fine search for finding the rotation center.

check_input
(dataset)[source]¶ Checks parameters of the input DataContainer
Should raise an Error if the DataContainer does not match expectation, e.g. if the expected input DataContainer is 3D and the Processor expects 2D.

static
find_center_vo
(tomo, ind=None, smin=40, smax=40, srad=10, step=0.5, ratio=2.0, drop=20)[source]¶ Find rotation axis location using Nghia Vo’s method. :cite:`Vo:14`.
 tomo : ndarray
 3D tomographic data.
 ind : int, optional
 Index of the slice to be used for reconstruction.
 smin, smax : int, optional
 Reference to the horizontal center of the sinogram.
 srad : float, optional
 Fine search radius.
 step : float, optional
 Step of fine searching.
 ratio : float, optional
 The ratio between the FOV of the camera and the size of object. It’s used to generate the mask.
 drop : int, optional
 Drop lines around vertical center of the mask.
 float
 Rotation axis location.
The function may not yield a correct estimate, if:
 the sample size is bigger than the field of view of the camera.
In this case the
ratio
argument need to be set larger than the default of 2.0.  there is distortion in the imaging hardware. If there’s no correction applied, the center of the projection image may yield a better estimate.
 the sample contrast is weak. Paganin’s filter need to be applied to overcome this.
 the sample was changed during the scan.