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: * :code:`ccpi.framework` module which allows to simply translate real world CT systems into software. * :code:`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. * :code:`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 CIL-Demos repository for full description. In conventional CT systems, an object is placed between a source emitting X-rays and a detector array measuring the X-ray transmission images of the incident X-rays. Typically, either the object is placed on a rotating sample stage and rotates with respect to the source-detector assembly, or the source-detector gantry rotates with respect to the stationary object. This arrangement results in so-called circular scanning trajectory. Depending on source and detector types, there are three conventional data acquisition geometries: * parallel geometry (2D or 3D), * fan-beam geometry, and * cone-beam geometry. Parallel geometry ----------------- Parallel beams of X-rays are emitted onto 1D (single pixel row) or 2D detector array. This geometry is common for synchrotron sources. 2D parrallel geometry is illustrated below. .. figure:: images/parallel.png :align: center :alt: alternate text :figclass: align-center 2D Parallel geometry .. figure:: images/parallel3d.png :align: center :alt: alternate text :figclass: align-center 3D Parallel geometry Fan-beam geometry ----------------- A single point-like X-ray source emits a cone beam onto 1D detector pixel row. Cone-beam is typically collimated to imaging field of view. Collimation allows greatly reduce amount of scatter radiation reaching the detector. Fan-beam geometry is used when scattering has significant influence on image quality or single-slice reconstruction is sufficient. .. figure:: images/fan.png :align: center :alt: alternate text :figclass: align-center Fan beam geometry Cone-beam geometry ------------------ A single point-like X-ray source emits a cone beam onto 2D detector array. Cone-beam geometry is mainly used in lab-based CT instruments. Depending on where the sample is placed between the source and the detector one can achieve a different magnification factor :math:`F`: .. math:: F = \frac{r_1 + r_2}{r_1} where :math:`r_1` and :math:`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. .. figure:: images/cone.png :align: center :alt: alternate text :figclass: align-center Cone beam geometry AcquisitonGeometry and AcquisitionData ====================================== In the Framework, we implemented :code:`AcquisitionGeometry` class to hold acquisition parameters and :code:`ImageGeometry` to hold geometry of a reconstructed volume. Corresponding data arrays are wrapped as :code:`AcquisitionData` and :code:`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: .. figure:: images/parallel_geometry.png :align: center :alt: alternate text :figclass: align-center Parallel geometry In the Framework, we define :code:`AcquisitionGeometry` as follows. .. code:: python # 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) :code:`AcquisitionGeometry` contains only metadata, the actual data is wrapped in :code:`AcquisitionData` class. :code:`AcquisitionGeometry` class also holds information about arrangement of the actual acquisition data array. \ We use attribute :code:`dimension_labels` to label axis. The expected dimension labels are shown below. The default order of dimensions for :code:`AcquisitionData` is :code:`[angle, horizontal]`, meaning that the number of elements along 0 and 1 axes in the acquisition data array is expected to be :code:`n_angles` and :code:`N`, respectively. .. figure:: images/parallel_data.png :align: center :alt: alternate text :figclass: align-center Parallel data To have consistent :code:`AcquisitionData` and :code:`AcquisitionGeometry`, we recommend to allocate :code:`AcquisitionData` using :code:`allocate` method of the :code:`AcquisitionGeometry` instance: .. code:: python # allocate AcquisitionData ad_par = ag_par.allocate() ImageGeometry and ImageData =========================== To store reconstruction results, we implemented two classes: :code:`ImageGeometry` and :code:`ImageData` classes. Similar to :code:`AcquisitionData` and :code:`AcquisitionGeometry`, we first define 2D :code:`ImageGeometry` and then allocate :code:`ImageData`. .. code:: python # 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) 3D parallel, fan-beam and cone-beam geometries ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Fan-beam, cone-beam and 3D (multi-slice) parallel geometry can be set-up similar to 2D parallel geometry. 3D parallel geometry -------------------- .. figure:: images/parallel3d_geometry.png :align: center :alt: alternate text :figclass: align-center Geometrical parameters and dimension labels for 3D parallel beam geometry 3D parallel beam :code:`AcquisitionGeometry` and default :code:`ImageGeometry` parameters can be set up as follows: .. code:: python # set-up 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) # set-up 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) Fan-beam geometry ----------------- .. figure:: images/fan_geometry.png :align: center :alt: alternate text :figclass: align-center Geometrical parameters and dimension labels for fan-beam geometry. .. figure:: images/fan_data.png :align: center :alt: alternate text :figclass: align-center Geometrical parameters and dimension labels for fan-beam data. Fan-beam :code:`AcquisitionGeometry` and default :code:`ImageGeometry` parameters can be set up as follows: .. code :: python # set-up fan-beam 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) .. autoclass:: ccpi.framework.ImageGeometry :members: .. autoclass:: ccpi.framework.AcquisitionGeometry :members: .. autoclass:: ccpi.framework.VectorGeometry :members: ======= ``DataContainer`` and subclasses ``AcquisitionData`` and ``ImageData`` are meant to contain data and meta-data in ``AcquisitionGeometry`` and ``ImageGeometry`` respectively. DataContainer and subclasses ============================ :code:`AcquisiionData` and :code:`ImageData` inherit from the same parent :code:`DataContainer` class, therefore they largely behave the same way. There are algebraic operations defined for both :code:`AcquisitionData` and :code:`ImageData`. Following operations are defined: * binary operations (between two DataContainers or scalar and DataContainer) * :code:`+` addition * :code:`-` subtraction * :code:`/` division * :code:`*` multiplication * :code:`**` power * :code:`maximum` * :code:`minimum` * in-place operations * :code:`+=` * :code:`-=` * :code:`*=` * :code:`**=` * :code:`/=` * unary operations * :code:`abs` * :code:`sqrt` * :code:`sign` * :code:`conjugate` * reductions * :code:`sum` * :code:`norm` * :code:`dot` product :code:`AcquisitionData` and :code:`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. .. code:: python # transpose data using subset method data_transposed = data.subset(['horizontal_y', 'horizontal_x']) # extract single row data_profile = data_subset.subset(horizontal_y=100) .. autoclass:: ccpi.framework.DataContainer :members: :private-members: :special-members: .. autoclass:: ccpi.framework.ImageData :members: .. autoclass:: ccpi.framework.AcquisitionData :members: .. autoclass:: ccpi.framework.VectorData :members: Multi channel data ------------------ Both :code:`AcquisitionGeometry`, :code:`AcquisitionData` and :code:`ImageGeometry`, :code:`ImageData` can be defined for multi-channel (spectral) CT data using :code:`channels` attribute. .. code:: python # multi-channel fan-beam 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 multi-channel 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 `_: .. math:: \underset{u}{\mathrm{argmin}}\begin{Vmatrix}A u - b \end{Vmatrix}^2_2 + \alpha^2\|Lu\|^2_2 where, * :math:`A` is the projection operator * :math:`b` is the acquired data * :math:`u` is the unknown image to be solved for * :math:`\alpha` is the regularisation parameter * :math:`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 :math:`L`. This can be re-written equivalently in the block matrix form: .. math:: \underset{u}{\mathrm{argmin}}\begin{Vmatrix}\binom{A}{\alpha L} u - \binom{b}{0}\end{Vmatrix}^2_2 With the definitions: * :math:`\tilde{A} = \binom{A}{\alpha L}` * :math:`\tilde{b} = \binom{b}{0}` this can now be recognised as a least squares problem which can be solved by any algorithm in the :code:`ccpi.optimisation` which can solve least squares problem, e.g. CGLS. .. math:: \underset{u}{\mathrm{argmin}}\begin{Vmatrix}\tilde{A} u - \tilde{b}\end{Vmatrix}^2_2 To be able to express our optimisation problems in the matrix form above, we developed the so-called, Block Framework comprising 4 main actors: :code:`BlockGeometry`, :code:`BlockDataContainer`, :code:`BlockFunction` and :code:`BlockOperator`. A :code:`BlockDataContainer` can be instantiated from a number of :code:`DataContainer` and subclasses represents a column vector of :code:`DataContainer`s. .. code:: python bdc = BlockDataContainer(DataContainer0, DataContainer1) . These classes are required for it to work. They provide a base class that will behave as normal ``DataContainer``. .. autoclass:: ccpi.framework.BlockDataContainer :members: :private-members: :special-members: .. autoclass:: ccpi.framework.BlockGeometry :members: :private-members: :special-members: DataProcessor ============= A :code:`DataProcessor` takes as an input a :code:`DataContainer` or subclass and returns either another :code:`DataContainer` or some number. The aim of this class is to simplify the writing of processing pipelines. .. autoclass:: ccpi.framework.DataProcessor :members: :private-members: :special-members: Resizer ------- Quite often we need either crop or downsample data; the :code:`Resizer` provides a convenient way to perform these operations for both :code:`ImageData` and :code:`AcquisitionData`. .. code:: python # 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 .. autoclass:: ccpi.processors.Resizer :members: :private-members: :special-members: 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 (double-borders). :code:`CenterOfRotationFinder` allows to estimate offset of center of rotation from theoretical. In the current release :code:`CenterOfRotationFinder` supports only parallel geometry. :code:`CenterOfRotationFinder` is based on Nghia Vo's `method `_. .. autoclass:: ccpi.processors.CenterOfRotationFinder :members: :private-members: :special-members: :ref:`Return Home ` .. _optimisation problems: optimisation.html