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
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.png
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.png
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.png
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.png
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.png
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.png
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.png
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.png
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.png
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
 ============