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/source | |
| 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/source')
| -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.pngBinary files differ new 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.pngBinary files differ new 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.pngBinary files differ new 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.pngBinary files differ new 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.pngBinary files differ new 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.pngBinary files differ new 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.pngBinary files differ new 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.pngBinary files differ new 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.pngBinary files differ new 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.pngBinary files differ new 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  ============ | 
