diff options
author | Edoardo Pasca <edo.paskino@gmail.com> | 2020-01-06 16:51:02 +0000 |
---|---|---|
committer | GitHub <noreply@github.com> | 2020-01-06 16:51:02 +0000 |
commit | f959a1c7f903fb31b40105f48701aadce2bd7b4c (patch) | |
tree | 6b596f6be0c36b301662b4fbb713bfe0c1e3d88a /docs | |
parent | 3d3a0958fad475c6b0493ad85459e1c04ba4ba62 (diff) | |
download | framework-f959a1c7f903fb31b40105f48701aadce2bd7b4c.tar.gz framework-f959a1c7f903fb31b40105f48701aadce2bd7b4c.tar.bz2 framework-f959a1c7f903fb31b40105f48701aadce2bd7b4c.tar.xz framework-f959a1c7f903fb31b40105f48701aadce2bd7b4c.zip |
v19.10 docs (#467)
updated docstrings and documentation
Diffstat (limited to 'docs')
-rw-r--r-- | docs/source/astra.rst | 3 | ||||
-rwxr-xr-x | docs/source/conf.py | 6 | ||||
-rw-r--r-- | docs/source/contrib.rst | 1 | ||||
-rw-r--r-- | docs/source/framework.rst | 474 | ||||
-rwxr-xr-x | docs/source/images/cone.png | bin | 0 -> 127928 bytes | |||
-rwxr-xr-x | docs/source/images/fan.png | bin | 0 -> 86375 bytes | |||
-rwxr-xr-x | docs/source/images/fan_data.png | bin | 0 -> 87766 bytes | |||
-rwxr-xr-x | docs/source/images/fan_geometry.png | bin | 0 -> 136263 bytes | |||
-rwxr-xr-x | docs/source/images/parallel.png | bin | 0 -> 29796 bytes | |||
-rwxr-xr-x | docs/source/images/parallel3d.png | bin | 0 -> 375145 bytes | |||
-rwxr-xr-x | docs/source/images/parallel3d_data.png | bin | 0 -> 371872 bytes | |||
-rwxr-xr-x | docs/source/images/parallel3d_geometry.png | bin | 0 -> 423629 bytes | |||
-rwxr-xr-x | docs/source/images/parallel_data.png | bin | 0 -> 21843 bytes | |||
-rwxr-xr-x | docs/source/images/parallel_geometry.png | bin | 0 -> 79825 bytes | |||
-rwxr-xr-x | docs/source/index.rst | 41 | ||||
-rw-r--r-- | docs/source/io.rst | 29 | ||||
-rw-r--r-- | docs/source/optimisation.rst | 282 | ||||
-rw-r--r-- | docs/source/plugins.rst | 2 |
18 files changed, 788 insertions, 50 deletions
diff --git a/docs/source/astra.rst b/docs/source/astra.rst index b80d2a4..a8759fd 100644 --- a/docs/source/astra.rst +++ b/docs/source/astra.rst @@ -22,7 +22,6 @@ Processors .. autoclass:: ccpi.astra.processors.AstraForwardProjectorMC :members: :special-members: -| Operators ========= @@ -35,7 +34,5 @@ Operators .. autoclass:: ccpi.astra.operators.AstraProjectorMC :members: :special-members: -| - :ref:`Return Home <mastertoc>` diff --git a/docs/source/conf.py b/docs/source/conf.py index 62790cc..b3084fa 100755 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -24,9 +24,9 @@ copyright = '2019, Edoardo Pasca' author = 'Edoardo Pasca'
# The short X.Y version
-version = '19.07'
+version = '19.10'
# The full version, including alpha/beta/rc tags
-release = '19.07'
+release = '19.10'
# -- General configuration ---------------------------------------------------
@@ -80,7 +80,7 @@ pygments_style = None # The theme to use for HTML and HTML Help pages. See the documentation for
# a list of builtin themes.
#
-html_theme = 'classic'
+html_theme = 'sphinx_rtd_theme'
# Theme options are theme-specific and customize the look and feel of a theme
# further. For a list of options available for each theme, see the
diff --git a/docs/source/contrib.rst b/docs/source/contrib.rst index 336097e..eaccc64 100644 --- a/docs/source/contrib.rst +++ b/docs/source/contrib.rst @@ -9,7 +9,6 @@ Contributed by Dr. Matthias Ehrhardt. .. autoclass:: ccpi.contrib.optimisation.algorithms.spdhg.spdhg :members: :special-members: -| :ref:`Return Home <mastertoc>` diff --git a/docs/source/framework.rst b/docs/source/framework.rst index 2b8ebf0..35d68fb 100644 --- a/docs/source/framework.rst +++ b/docs/source/framework.rst @@ -1,9 +1,339 @@ 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 <https://github.com/vais-ral/CIL-Demos/blob/v19.10.1/Notebooks/00_building_blocks.ipynb>`_ 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: @@ -15,37 +345,153 @@ DataContainer and subclasses .. autoclass:: ccpi.framework.VectorData :members: -.. autoclass:: ccpi.framework.ImageGeometry - :members: -.. autoclass:: ccpi.framework.AcquisitionGeometry - :members: -.. autoclass:: ccpi.framework.VectorGeometry - :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 <https://en.wikipedia.org/wiki/Tikhonov_regularization>`_: + +.. 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.CenterOfRotationFinder - :members: -.. autoclass:: ccpi.processors.Normalizer - :members: .. 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 <https://doi.org/10.1364/OE.22.019078>`_. + +.. autoclass:: ccpi.processors.CenterOfRotationFinder + :members: + :private-members: + :special-members: + :ref:`Return Home <mastertoc>` + +.. _optimisation problems: optimisation.html diff --git a/docs/source/images/cone.png b/docs/source/images/cone.png Binary files differnew file mode 100755 index 0000000..bd8896f --- /dev/null +++ b/docs/source/images/cone.png diff --git a/docs/source/images/fan.png b/docs/source/images/fan.png Binary files differnew file mode 100755 index 0000000..4f20da4 --- /dev/null +++ b/docs/source/images/fan.png diff --git a/docs/source/images/fan_data.png b/docs/source/images/fan_data.png Binary files differnew file mode 100755 index 0000000..cc7a470 --- /dev/null +++ b/docs/source/images/fan_data.png diff --git a/docs/source/images/fan_geometry.png b/docs/source/images/fan_geometry.png Binary files differnew file mode 100755 index 0000000..c1b6a50 --- /dev/null +++ b/docs/source/images/fan_geometry.png diff --git a/docs/source/images/parallel.png b/docs/source/images/parallel.png Binary files differnew file mode 100755 index 0000000..a58f79e --- /dev/null +++ b/docs/source/images/parallel.png diff --git a/docs/source/images/parallel3d.png b/docs/source/images/parallel3d.png Binary files differnew file mode 100755 index 0000000..f5dc76f --- /dev/null +++ b/docs/source/images/parallel3d.png diff --git a/docs/source/images/parallel3d_data.png b/docs/source/images/parallel3d_data.png Binary files differnew file mode 100755 index 0000000..2b5536a --- /dev/null +++ b/docs/source/images/parallel3d_data.png diff --git a/docs/source/images/parallel3d_geometry.png b/docs/source/images/parallel3d_geometry.png Binary files differnew file mode 100755 index 0000000..fdcff6f --- /dev/null +++ b/docs/source/images/parallel3d_geometry.png diff --git a/docs/source/images/parallel_data.png b/docs/source/images/parallel_data.png Binary files differnew file mode 100755 index 0000000..7adea39 --- /dev/null +++ b/docs/source/images/parallel_data.png diff --git a/docs/source/images/parallel_geometry.png b/docs/source/images/parallel_geometry.png Binary files differnew file mode 100755 index 0000000..2880b55 --- /dev/null +++ b/docs/source/images/parallel_geometry.png diff --git a/docs/source/index.rst b/docs/source/index.rst index 654a083..266a03a 100755 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -6,6 +6,24 @@ Welcome to CCPi-Framework's documentation!
==========================================
+The aim of this package is to enable rapid prototyping of optimisation-based
+reconstruction problems, i.e. defining and solving different optimization problems to enforce different properties on the reconstructed image, while being
+powerful enough to be employed on real scale problems.
+
+Firstly, it provides a framework to handle acquisition and reconstruction
+data and metadata; it also provides a basic input/output package to read data
+from different sources, e.g. Nikon X-Radia, NeXus.
+
+Secondly, it provides an object-oriented framework for defining mathematical
+operators and functions as well a collection of useful example operators and
+functions. Both smooth and non-smooth functions can be used.
+
+Further, it provides a number of high-level generic implementations of
+optimisation algorithms to solve genericlly formulated optimisation problems
+constructed from operator and function objects.
+
+A number of demos can be found on the `CIL-Demos`_ repository.
+
.. toctree::
:maxdepth: 2
:caption: Contents:
@@ -13,15 +31,26 @@ Welcome to CCPi-Framework's documentation! framework
- optimisation
io
+ optimisation
plugins
astra
contrib
-Indices and tables
-==================
+.. Indices and tables
+.. ==================
+
+.. * :ref:`genindex`
+.. * :ref:`modindex`
+.. * :ref:`search`
+
+Contacts
+========
+
+Please refer to the main `CCPi website`_ for up-to-date information.
+
+The CCPi developers may be contacted joining the `devel mailing list`_
-* :ref:`genindex`
-* :ref:`modindex`
-* :ref:`search`
+.. _devel mailing list: https://www.jiscmail.ac.uk/cgi-bin/webadmin?A0=CCPI-DEVEL
+.. _CCPi website: https://www.ccpi.ac.uk
+.. _CIL-Demos: https://github.com/vais-ral/CIL-Demos
diff --git a/docs/source/io.rst b/docs/source/io.rst index fb24a3a..9ac78a4 100644 --- a/docs/source/io.rst +++ b/docs/source/io.rst @@ -1,9 +1,34 @@ -Input/Output -************ +Read/ write AcquisitionData and ImageData +***************************************** + NeXus ===== +The CCPi Framework provides classes to read and write :code:`AcquisitionData` and :code:`ImageData` +as NeXuS files. + +.. code:: python + + # imports + from ccpi.io import NEXUSDataWriter, NEXUSDataReader + + # initialise NEXUS Writer + writer = NEXUSDataWriter() + writer.set_up(file_name='tmp_nexus.nxs', + data_container=my_data) + # write data + writer.write_file() + + # read data + # initialize NEXUS reader + reader = NEXUSDataReader() + reader.set_up(nexus_file='tmp_nexus.nxs') + # load data + ad1 = reader.load_data() + # get AcquisiionGeometry + ag1 = reader.get_geometry() + .. autoclass:: ccpi.io.NEXUSDataReader :members: :special-members: diff --git a/docs/source/optimisation.rst b/docs/source/optimisation.rst index eec54e1..59f3dd3 100644 --- a/docs/source/optimisation.rst +++ b/docs/source/optimisation.rst @@ -8,9 +8,112 @@ Further, it provides a number of high-level generic implementations of optimisat The fundamental components are: -+ Operator: A class specifying a (currently linear) operator -+ Function: A class specifying mathematical functions such as a least squares data fidelity. -+ Algorithm: Implementation of an iterative optimisation algorithm to solve a particular generic optimisation problem. Algorithms are iterable Python object which can be run in a for loop. Can be stopped and warm restarted. ++ :code:`Operator`: A class specifying a (currently linear) operator ++ :code:`Function`: A class specifying mathematical functions such as a least squares data fidelity. ++ :code:`Algorithm`: Implementation of an iterative optimisation algorithm to solve a particular generic optimisation problem. Algorithms are iterable Python object which can be run in a for loop. Can be stopped and warm restarted. + +To be able to express more advanced optimisation problems we developed the +`Block Framework`_, which provides a generic strategy to treat variational +problems in the following form: + +.. math:: + \min \text{Regulariser} + \text{Fidelity} + +The block framework consists of: + ++ BlockDataContainer ++ BlockFunction ++ BlockOperator + +`BlockDataContainer`_ holds `DataContainer` as column vector. It is possible to +do basic algebra between ``BlockDataContainer`` s and with numbers, list or numpy arrays. + +`BlockFunction`_ acts on ``BlockDataContainer`` as a separable sum function: + + .. math:: + + f = [f_1,...,f_n] \newline + + f([x_1,...,x_n]) = f_1(x_1) + .... + f_n(x_n) + +`BlockOperator`_ represent a block matrix with operators + +.. math:: + K = \begin{bmatrix} + A_{1} & A_{2} \\ + A_{3} & A_{4} \\ + A_{5} & A_{6} + \end{bmatrix}_{(3,2)} * \quad \underbrace{\begin{bmatrix} + x_{1} \\ + x_{2} + \end{bmatrix}_{(2,1)}}_{\textbf{x}} = \begin{bmatrix} + A_{1}x_{1} + A_{2}x_{2}\\ + A_{3}x_{1} + A_{4}x_{2}\\ + A_{5}x_{1} + A_{6}x_{2}\\ + \end{bmatrix}_{(3,1)} = \begin{bmatrix} + y_{1}\\ + y_{2}\\ + y_{3} + \end{bmatrix}_{(3,1)} = \textbf{y} + +Column: Share the same domains :math:`X_{1}, X_{2}` + +Rows: Share the same ranges :math:`Y_{1}, Y_{2}, Y_{3}` + +.. math:: + K : (X_{1}\times X_{2}) \rightarrow (Y_{1}\times Y_{2} \times Y_{3}) + +:math:`A_{1}, A_{3}, A_{5}`: share the same domain :math:`X_{1}` and +:math:`A_{2}, A_{4}, A_{6}`: share the same domain :math:`X_{2}` + +.. math:: + + A_{1}: X_{1} \rightarrow Y_{1} \\ + A_{3}: X_{1} \rightarrow Y_{2} \\ + A_{5}: X_{1} \rightarrow Y_{3} \\ + A_{2}: X_{2} \rightarrow Y_{1} \\ + A_{4}: X_{2} \rightarrow Y_{2} \\ + A_{6}: X_{2} \rightarrow Y_{3} + +For instance with these ingredients one may write the following objective +function, + +.. math:: + \alpha ||\nabla u||_{2,1} + ||u - g||_2^2 + +where :math:`g` represent the measured values, :math:`u` the solution +:math:`\nabla` is the gradient operator, :math:`|| ~~ ||_{2,1}` is a norm for +the output of the gradient operator and :math:`|| x-g ||^2_2` is +least squares fidelity function as + +.. math:: + K = \begin{bmatrix} + \nabla \\ + \mathbb{1} + \end{bmatrix} + + F(x) = \Big[ \alpha \lVert ~x~ \rVert_{2,1} ~~ , ~~ || x - g||_2^2 \Big] + + w = [ u ] + +Then we have rewritten the problem as + +.. math:: + F(Kw) = \alpha \left\lVert \nabla u \right\rVert_{2,1} + ||u-g||^2_2 + +Which in Python would be like + +.. code-block:: python + + op1 = Gradient(ig, correlation=Gradient.CORRELATION_SPACE) + op2 = Identity(ig, ag) + + # Create BlockOperator + K = BlockOperator(op1, op2, shape=(2,1) ) + + # Create functions + F = BlockFunction(alpha * MixedL21Norm(), 0.5 * L2NormSquared(b=noisy_data)) + Algorithm ========= @@ -22,12 +125,13 @@ Gradient (PDHG) and Fast Iterative Shrinkage Thresholding Algorithm (FISTA). An algorithm is designed for a particular generic optimisation problem accepts and number of -Functions and/or Operators as input to define a specific instance of +:code:`Function`s and/or :code:`Operator`s as input to define a specific instance of the generic optimisation problem to be solved. They are iterable objects which can be run in a for loop. The user can provide a stopping criterion different than the default max_iteration. -New algorithms can be easily created by extending the Algorithm class. The user is required to implement only 4 methods: set_up, __init__, update and update_objective. +New algorithms can be easily created by extending the :code:`Algorithm` class. +The user is required to implement only 4 methods: set_up, __init__, update and update_objective. + :code:`set_up` and :code:`__init__` are used to configure the algorithm + :code:`update` is the actual iteration updating the solution @@ -43,7 +147,9 @@ algorithm to minimise a Function will only be: def update_objective(self): self.loss.append(self.objective_function(self.x)) -The :code:`Algorithm` provides the infrastructure to continue iteration, to access the values of the objective function in subsequent iterations, the time for each iteration. +The :code:`Algorithm` provides the infrastructure to continue iteration, to access the values of the +objective function in subsequent iterations, the time for each iteration, and to provide a nice +print to screen of the status of the optimisation. .. autoclass:: ccpi.optimisation.algorithms.Algorithm :members: @@ -55,6 +161,7 @@ The :code:`Algorithm` provides the infrastructure to continue iteration, to acce :members: .. autoclass:: ccpi.optimisation.algorithms.FISTA :members: + :special-members: .. autoclass:: ccpi.optimisation.algorithms.PDHG :members: .. autoclass:: ccpi.optimisation.algorithms.SIRT @@ -69,6 +176,14 @@ The output is another :code:`DataContainer` object or subclass hereof. An important special case is to represent the tomographic forward and backprojection operations. + +Operator base classes +--------------------- + +All operators extend the :code:`Operator` class. A special class is the :code:`LinearOperator` +which represents an operator for which the :code:`adjoint` operation is defined. +A :code:`ScaledOperator` represents the multiplication of any operator with a scalar. + .. autoclass:: ccpi.optimisation.operators.Operator :members: :special-members: @@ -78,35 +193,57 @@ forward and backprojection operations. .. autoclass:: ccpi.optimisation.operators.ScaledOperator :members: :special-members: -.. autoclass:: ccpi.optimisation.operators.GradientOperator - :members: - :special-members: + +Trivial operators +----------------- + +Trivial operators are the following. + .. autoclass:: ccpi.optimisation.operators.Identity :members: :special-members: -.. autoclass:: ccpi.optimisation.operators.LinearOperatorMatrix + +.. autoclass:: ccpi.optimisation.operators.ZeroOperator :members: :special-members: -.. autoclass:: ccpi.optimisation.operators.ShrinkageOperator + +.. autoclass:: ccpi.optimisation.operators.LinearOperatorMatrix :members: :special-members: -.. autoclass:: ccpi.optimisation.operators.SparseFiniteDiff + + +Gradient +----------------- + +In the following the required classes for the implementation of the :code:`Gradient` operator. + +.. autoclass:: ccpi.optimisation.operators.Gradient :members: :special-members: -.. autoclass:: ccpi.optimisation.operators.SymmetrizedGradientOperator + +.. autoclass:: ccpi.optimisation.operators.FiniteDiff :members: :special-members: -.. autoclass:: ccpi.optimisation.operators.ZeroOperator + +.. autoclass:: ccpi.optimisation.operators.SparseFiniteDiff :members: :special-members: -.. autoclass:: ccpi.optimisation.operators.BlockOperator + +.. autoclass:: ccpi.optimisation.operators.SymmetrizedGradient :members: :special-members: -.. autoclass:: ccpi.optimisation.operators.BlockScaledOperator + + +Shrinkage operator +------------------ + +.. autoclass:: ccpi.optimisation.operators.ShrinkageOperator :members: :special-members: + + Function ======== @@ -124,36 +261,143 @@ input point. The function value is evaluated by calling the function itself, e.g. :code:`f(x)` for a :code:`Function f` and input point :code:`x`. +Base classes +------------ + .. autoclass:: ccpi.optimisation.functions.Function :members: :special-members: + +.. autoclass:: ccpi.optimisation.functions.ScaledFunction + :members: + :special-members: + +Composition of operator and a function +-------------------------------------- + +This class allows the user to write a function which does the following: + +.. math:: + + F ( x ) = G ( Ax ) + +where :math:`A` is an operator. For instance the least squares function l2norm_ :code:`Norm2Sq` can +be expressed as + +.. math:: + + F(x) = || Ax - b ||^2_2 + +.. code::python + + F1 = Norm2Sq(A, b) + # or equivalently + F2 = FunctionOperatorComposition(L2NormSquared(b=b), A) + + .. autoclass:: ccpi.optimisation.functions.FunctionOperatorComposition :members: :special-members: + +Indicator box +------------- + .. autoclass:: ccpi.optimisation.functions.IndicatorBox :members: :special-members: + + +KullbackLeibler +--------------- + .. autoclass:: ccpi.optimisation.functions.KullbackLeibler :members: :special-members: + +L1 Norm +------- + .. autoclass:: ccpi.optimisation.functions.L1Norm :members: :special-members: + +Squared L2 norm +--------------- +.. l2norm: + .. autoclass:: ccpi.optimisation.functions.L2NormSquared :members: :special-members: + +And a least squares function: + +.. autoclass:: ccpi.optimisation.functions.Norm2Sq + :members: + :special-members: + +Mixed L21 norm +-------------- + .. autoclass:: ccpi.optimisation.functions.MixedL21Norm :members: :special-members: -.. autoclass:: ccpi.optimisation.functions.Norm2Sq + +.. autoclass:: ccpi.optimisation.functions.ZeroFunction :members: :special-members: -.. autoclass:: ccpi.optimisation.functions.ScaledFunction + + +Block Framework +*************** + +Block Operator +============== + + +.. autoclass:: ccpi.optimisation.operators.BlockOperator :members: :special-members: -.. autoclass:: ccpi.optimisation.functions.ZeroFunction +.. autoclass:: ccpi.optimisation.operators.BlockScaledOperator + :members: + :special-members: + + +Block Function +============== +A Block vector of functions, Size of vector coincides with the rows of :math:`K`: + +.. math:: + + Kx = \begin{bmatrix} + y_{1}\\ + y_{2}\\ + y_{3}\\ + \end{bmatrix}, \quad f = [ f_{1}, f_{2}, f_{3} ] + + f(Kx) : = f_{1}(y_{1}) + f_{2}(y_{2}) + f_{3}(y_{3}) + + +.. autoclass:: ccpi.optimisation.functions.BlockFunction :members: :special-members: +Block DataContainer +============== + +.. math:: + + x = [x_{1}, x_{2} ]\in (X_{1}\times X_{2}) + + y = [y_{1}, y_{2}, y_{3} ]\in(Y_{1}\times Y_{2} \times Y_{3}) + + +.. autoclass:: ccpi.framework.BlockDataContainer + :members: + :special-members: + :ref:`Return Home <mastertoc>` + +.. _BlockDataContainer: framework.html#ccpi.framework.BlockDataContainer +.. _BlockFunction: optimisation.html#ccpi.optimisation.functions.BlockFunction +.. _BlockOperator: optimisation.html#ccpi.optimisation.operators.BlockOperators diff --git a/docs/source/plugins.rst b/docs/source/plugins.rst index 948980c..4348f62 100644 --- a/docs/source/plugins.rst +++ b/docs/source/plugins.rst @@ -7,7 +7,6 @@ Operators .. autoclass:: ccpi.plugins.operators.CCPiProjectorSimple :members: :special-members: -| Processors ========== @@ -23,7 +22,6 @@ Processors .. autoclass:: ccpi.plugins.processors.setupCCPiGeometries :members: :special-members: -| Regularisers ============ |