v19.10 docs (#467)

updated docstrings and documentation

37 files changed, 1022 insertions, 187 deletions

diff --git a/README.md b/README.md
index cb4f856..7b52ba5 100644
--- a/README.md
+++ b/README.md
@@ -11,11 +11,11 @@ This package provides a common framework, hence the name, for the analysis of da
 
 ## Installation
 
-Binary installation of the CCPi Framework can be done with `conda`. Install a new environment as in the following. Some packages are optional (cvxpy and ccpi-astra)
+Binary installation of the CCPi Framework can be done with `conda`. Install a new environment as in the following. Some packages are optional (`cvxpy`, `ccpi-plugins`, `ccpi-astra`) but very useful.
 
 ```bash
 
-conda create -y --name cil python=3.6.7 cvxpy=1.0.15 lapack numpy=1.14 ccpi-framework ccpi-astra tomophantom ccpi-plugins -c oxfordcontrol -c conda-forge -c astra-toolbox -c ccpi/label/dev
+conda create -y --name cil python=3 cvxpy ccpi-framework=19.10 ccpi-astra=19.10 tomophantom ccpi-plugins=19.10 -c oxfordcontrol -c conda-forge -c astra-toolbox -c ccpi
 ```
 
 ### Components
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py
index 408a904..48d109e 100755
--- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py
@@ -30,14 +30,15 @@ class Algorithm(object):
     '''Base class for iterative algorithms
 
       provides the minimal infrastructure.
+
       Algorithms are iterables so can be easily run in a for loop. They will
       stop as soon as the stop cryterion is met.
-      The user is required to implement the set_up, __init__, update and
-      and update_objective methods
+      The user is required to implement the :code:`set_up`, :code:`__init__`, :code:`update` and
+      and :code:`update_objective` methods
       
-      A courtesy method run is available to run n iterations. The method accepts
-      a callback function that receives the current iteration number and the actual objective
-      value and can be used to trigger print to screens and other user interactions. The run
+      A courtesy method :code:`run` is available to run :code:`n` iterations. The method accepts
+      a :code:`callback` function that receives the current iteration number and the actual objective
+      value and can be used to trigger print to screens and other user interactions. The :code:`run`
       method will stop when the stopping cryterion is met. 
    '''
 
@@ -45,14 +46,15 @@ class Algorithm(object):
         '''Constructor
         
         Set the minimal number of parameters:
-            iteration: current iteration number
-            max_iteration: maximum number of iterations
-            memopt: whether to use memory optimisation ()
-            timing: list to hold the times it took to run each iteration
-            update_objectice_interval: the interval every which we would save the current
-                                       objective. 1 means every iteration, 2 every 2 iteration
-                                       and so forth. This is by default 1 and should be increased
+        
+        
+        :param max_iteration: maximum number of iterations
+        :type max_iteration: int, optional, default 0
+        :param update_objectice_interval: the interval every which we would save the current\
+                                       objective. 1 means every iteration, 2 every 2 iteration\
+                                       and so forth. This is by default 1 and should be increased\
                                        when evaluating the objective is computationally expensive.
+        :type update_objective_interval: int, optional, default 1
         '''
         self.iteration = 0
         self.__max_iteration = kwargs.get('max_iteration', 0)
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py
index 57292df..53804c5 100755
--- a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py
@@ -64,9 +64,9 @@ class CGLS(Algorithm):
     def set_up(self, x_init, operator, data, tolerance=1e-6):
         '''initialisation of the algorithm
 
-        :param operator : Linear operator for the inverse problem
-        :param x_init : Initial guess ( Default x_init = 0)
-        :param data : Acquired data to reconstruct       
+        :param operator: Linear operator for the inverse problem
+        :param x_init: Initial guess ( Default x_init = 0)
+        :param data: Acquired data to reconstruct       
         :param tolerance: Tolerance/ Stopping Criterion to end CGLS algorithm
         '''
         print("{} setting up".format(self.__class__.__name__, ))
@@ -94,6 +94,7 @@ class CGLS(Algorithm):
 
         
     def update(self):
+        '''single iteration'''
         
         self.q = self.operator.direct(self.p)
         delta = self.q.squared_norm()
@@ -121,9 +122,11 @@ class CGLS(Algorithm):
         self.loss.append(a)
         
     def should_stop(self):
+        '''stopping criterion'''
         return self.flag() or self.max_iteration_stop_cryterion()
  
     def flag(self):
+        '''returns whether the tolerance has been reached'''
         flag  = (self.norms <= self.norms0 * self.tolerance) or (self.normx * self.tolerance >= 1)
 
         if flag:
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
index 8c485b7..15a289d 100755
--- a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
@@ -40,9 +40,9 @@ class FISTA(Algorithm):
     
     Parameters :
         
-      :parameter x_init : Initial guess ( Default x_init = 0)
-      :parameter f : Differentiable function
-      :parameter g : Convex function with " simple " proximal operator
+      :param x_init: Initial guess ( Default x_init = 0)
+      :param f: Differentiable function
+      :param g: Convex function with " simple " proximal operator
 
 
     Reference:
@@ -60,9 +60,11 @@ class FISTA(Algorithm):
         initialisation can be done at creation time if all 
         proper variables are passed or later with set_up
         
-        :param x_init : Initial guess ( Default x_init = 0)
-        :param f : Differentiable function
-        :param g : Convex function with " simple " proximal operator'''
+        Optional parameters:
+
+        :param x_init: Initial guess ( Default x_init = 0)
+        :param f: Differentiable function
+        :param g: Convex function with " simple " proximal operator'''
         
         super(FISTA, self).__init__(**kwargs)
         
@@ -72,9 +74,9 @@ class FISTA(Algorithm):
     def set_up(self, x_init, f, g=ZeroFunction()):
         '''initialisation of the algorithm
 
-        :param x_init : Initial guess ( Default x_init = 0)
-        :param f : Differentiable function
-        :param g : Convex function with " simple " proximal operator'''
+        :param x_init: Initial guess ( Default x_init = 0)
+        :param f: Differentiable function
+        :param g: Convex function with " simple " proximal operator'''
 
         print("{} setting up".format(self.__class__.__name__, ))
         
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
index db1b8dc..8776875 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
@@ -34,22 +34,21 @@ class PDHG(Algorithm):
     .. math::
     
       \min_{x} f(Kx) + g(x)
-    |
-
-    Parameters : 
         
-        :parameter operator : Linear Operator = K
-        :parameter f : Convex function with "simple" proximal of its conjugate. 
-        :parameter g : Convex function with "simple" proximal 
-        :parameter sigma : Step size parameter for Primal problem
-        :parameter tau : Step size parameter for Dual problem
+    :param operator: Linear Operator = K
+    :param f: Convex function with "simple" proximal of its conjugate. 
+    :param g: Convex function with "simple" proximal 
+    :param sigma: Step size parameter for Primal problem
+    :param tau: Step size parameter for Dual problem
         
-        Remark: Convergence is guaranted provided that
+    Remark: Convergence is guaranted provided that
         
-        .. math:: \tau \sigma \|K\|^{2} <1
+    .. math:: 
+    
+      \tau \sigma \|K\|^{2} <1
         
             
-    Reference :
+    Reference:
         
         
         (a) A. Chambolle and T. Pock (2011), "A first-order primal–dual algorithm for convex
@@ -64,11 +63,14 @@ class PDHG(Algorithm):
     def __init__(self, f=None, g=None, operator=None, tau=None, sigma=1.,**kwargs):
         '''PDHG algorithm creator
 
-        :param operator : Linear Operator = K
-        :param f : Convex function with "simple" proximal of its conjugate. 
-        :param g : Convex function with "simple" proximal 
-        :param sigma : Step size parameter for Primal problem
-        :param tau : Step size parameter for Dual problem'''
+        Optional parameters
+
+        :param operator: a Linear Operator
+        :param f: Convex function with "simple" proximal of its conjugate. 
+        :param g: Convex function with "simple" proximal 
+        :param sigma: Step size parameter for Primal problem
+        :param tau: Step size parameter for Dual problem
+        '''
         super(PDHG, self).__init__(**kwargs)
         
 
@@ -78,11 +80,11 @@ class PDHG(Algorithm):
     def set_up(self, f, g, operator, tau=None, sigma=1.):
         '''initialisation of the algorithm
 
-        :param operator : Linear Operator = K
-        :param f : Convex function with "simple" proximal of its conjugate. 
-        :param g : Convex function with "simple" proximal 
-        :param sigma : Step size parameter for Primal problem
-        :param tau : Step size parameter for Dual problem'''
+        :param operator: a Linear Operator
+        :param f: Convex function with "simple" proximal of its conjugate. 
+        :param g: Convex function with "simple" proximal 
+        :param sigma: Step size parameter for Primal problem
+        :param tau: Step size parameter for Dual problem'''
 
         print("{} setting up".format(self.__class__.__name__, ))
         
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py
index 50398f4..a59ce5f 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py
@@ -35,27 +35,26 @@ class SIRT(Algorithm):
     
     .. math::  
     
-    A x = b
-    |
-
-    Parameters:
-        
-      :parameter operator : Linear operator for the inverse problem
-      :parameter x_init : Initial guess
-      :parameter data : Acquired data to reconstruct       
-      :parameter constraint : Function proximal method
-                   e.g.  x\in[0, 1], IndicatorBox to enforce box constraints
-                         Default is None).
+      A x = b
+    
+    :param x_init: Initial guess
+    :param operator: Linear operator for the inverse problem
+    :param data: Acquired data to reconstruct       
+    :param constraint: Function proximal method
+                e.g.  :math:`x\in[0, 1]`, :code:`IndicatorBox` to enforce box constraints
+                        Default is :code:`None`).
     '''
     def __init__(self, x_init=None, operator=None, data=None, constraint=None, **kwargs):
         '''SIRT algorithm creator
 
-        :param x_init : Initial guess
-        :param operator : Linear operator for the inverse problem
-        :param data : Acquired data to reconstruct       
-        :param constraint : Function proximal method
-                   e.g.  x\in[0, 1], IndicatorBox to enforce box constraints
-                         Default is None).'''
+       Optional parameters:
+
+      :param x_init: Initial guess
+      :param operator: Linear operator for the inverse problem
+      :param data: Acquired data to reconstruct       
+      :param constraint: Function proximal method
+                   e.g.  :math:`x\in[0, 1]`, :code:`IndicatorBox` to enforce box constraints
+                         Default is :code:`None`).'''
         super(SIRT, self).__init__(**kwargs)
 
         if x_init is not None and operator is not None and data is not None:
@@ -64,12 +63,12 @@ class SIRT(Algorithm):
     def set_up(self, x_init, operator, data, constraint=None):
         '''initialisation of the algorithm
 
-        :param operator : Linear operator for the inverse problem
-        :param x_init : Initial guess
-        :param data : Acquired data to reconstruct       
-        :param constraint : Function proximal method
-                   e.g.  x\in[0, 1], IndicatorBox to enforce box constraints
-                         Default is None).'''
+        :param x_init: Initial guess
+        :param operator: Linear operator for the inverse problem
+        :param data: Acquired data to reconstruct       
+        :param constraint: Function proximal method
+                   e.g.  :math:`x\in[0, 1]`, :code:`IndicatorBox` to enforce box constraints
+                         Default is :code:`None`).'''
         print("{} setting up".format(self.__class__.__name__, ))
         
         self.x = x_init.copy()
diff --git a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py
index ee3ad78..57592cd 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py
@@ -52,7 +52,7 @@ class BlockFunction(Function):
         
             :param: x (BlockDataContainer): must have as many rows as self.length
 
-            returns ..math:: sum(f_i(x_i))
+            returns ..math:: \sum(f_i(x_i))
             
         '''
         
@@ -67,7 +67,7 @@ class BlockFunction(Function):
         
         r'''Convex conjugate of BlockFunction at x            
         
-            .. math:: returns sum(f_i^{*}(x_i))
+            .. math:: returns \sum(f_i^{*}(x_i))
         
         '''       
         t = 0                
@@ -80,7 +80,7 @@ class BlockFunction(Function):
         
         r'''Proximal operator of BlockFunction at x: 
                  
-                 .. math:: prox_{tau*f}(x) = sum_{i} prox_{tau*f_{i}}(x_{i}) 
+                 .. math:: prox_{tau*f}(x) = \sum_{i} prox_{tau*f_{i}}(x_{i}) 
         
         
         '''
@@ -110,7 +110,7 @@ class BlockFunction(Function):
         
         r'''Proximal operator of the convex conjugate of BlockFunction at x:
         
-            .. math:: prox_{tau*f^{*}}(x) = sum_{i} prox_{tau*f^{*}_{i}}(x_{i}) 
+            .. math:: prox_{tau*f^{*}}(x) = \sum_{i} prox_{tau*f^{*}_{i}}(x_{i}) 
         '''
         
         if out is None:
diff --git a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py
index 4162134..ed5c1b1 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py
@@ -26,15 +26,24 @@ from ccpi.optimisation.functions import ScaledFunction
 
 class FunctionOperatorComposition(Function):
     
-    '''Function composition with Operator: (f o A)(x) = f(Ax)
+    r'''Function composition with Operator: :math:`(f \otimes A)(x) = f(Ax)`
     
-            : parameter A: operator
-            : parameter f: function
+    :param A: operator
+    :type A: :code:`Operator`
+    :param f: function
+    :type f: :code:`Function`
     
     '''
     
     def __init__(self, function, operator):
-        
+        '''creator
+
+    :param A: operator
+    :type A: :code:`Operator`
+    :param f: function
+    :type f: :code:`Function`
+    '''
+
         super(FunctionOperatorComposition, self).__init__()
         
         self.function = function     
diff --git a/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py
index ac8978a..9e9e55c 100755
--- a/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py
@@ -30,20 +30,25 @@ class IndicatorBox(Function):
     r'''Indicator function for box constraint
             
       .. math:: 
-         f(x) = \mathbb{I}_{[a, b]} = \begin{cases}
-            
-                                            0, if x\in[a, b]
-                                            \infty, otherwise                            
-                                    \end{cases}
+         
+         f(x) = \mathbb{I}_{[a, b]} = \begin{cases}  
+                                            0, \text{ if } x \in [a, b] \\
+                                            \infty, \text{otherwise}
+                                     \end{cases}
     
     '''
     
     def __init__(self,lower=-numpy.inf,upper=numpy.inf):
+        '''creator
 
+        :param lower: lower bound
+        :type lower: float, default = :code:`-numpy.inf`
+        :param upper: upper bound
+        :type upper: float, optional, default = :code:`numpy.inf`
         super(IndicatorBox, self).__init__()
         self.lower = lower
         self.upper = upper
-        
+        '''
 
     def __call__(self,x):
         
diff --git a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py
index 1c2c43f..1fcfcca 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py
@@ -32,13 +32,21 @@ class L1Norm(Function):
     r'''L1Norm function: 
             
             Cases considered (with/without data):            
-                a) .. math:: f(x) = ||x||_{1}
-                b) .. math:: f(x) = ||x - b||_{1}
+                a) :math:`f(x) = ||x||_{1}`
+                b) :math:`f(x) = ||x - b||_{1}`
                                 
     '''   
            
     def __init__(self, **kwargs):
-        
+        '''creator
+
+        Cases considered (with/without data):            
+        a) :math:`f(x) = ||x||_{1}`
+        b) :math:`f(x) = ||x - b||_{1}`
+
+        :param b: translation of the function
+        :type b: :code:`DataContainer`, optional
+        '''
         super(L1Norm, self).__init__()
         self.b = kwargs.get('b',None)
         self.shinkage_operator = ShrinkageOperator()
diff --git a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py
index f5108ba..ef7c698 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py
@@ -37,7 +37,15 @@ class L2NormSquared(Function):
     '''    
     
     def __init__(self, **kwargs):
-                                
+        '''creator
+
+        Cases considered (with/without data):            
+                a) .. math:: f(x) = \|x\|^{2}_{2} 
+                b) .. math:: f(x) = \|\|x - b\|\|^{2}_{2}
+
+        :param b:  translation of the function
+        :type b: :code:`DataContainer`, optional
+        '''                        
         super(L2NormSquared, self).__init__()
         self.b = kwargs.get('b',None) 
         
diff --git a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
index 55e6e53..1af0e77 100755
--- a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
@@ -31,11 +31,14 @@ class MixedL21Norm(Function):
     
     
     '''
-        f(x) = ||x||_{2,1} = \sum |x|_{2}                   
+        .. math::
+          
+          f(x) = ||x||_{2,1} = \sum |x|_{2}
     '''      
     
     def __init__(self, **kwargs):
-
+        '''creator        
+        '''
         super(MixedL21Norm, self).__init__()                      
         self.SymTensor = kwargs.get('SymTensor',False)
         
@@ -43,7 +46,7 @@ class MixedL21Norm(Function):
         
         ''' Evaluates L2,1Norm at point x
             
-            :param: x is a BlockDataContainer
+            :param x: is a BlockDataContainer
                                 
         '''
         if not isinstance(x, BlockDataContainer):
@@ -60,8 +63,9 @@ class MixedL21Norm(Function):
                             
     def convex_conjugate(self,x):
         
-        ''' This is the Indicator function of ||\cdot||_{2, \infty}
-            which is either 0 if ||x||_{2, \infty} or \infty        
+        ''' This is the Indicator function of :math:`||\cdot||_{2, \infty}` which is either 0 if :math:`||x||_{2, \infty}` or :math:`\infty`
+
+        Notice this returns 0
         '''
         
         return 0.0
diff --git a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
index 0dd7d4b..3c563fb 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
@@ -27,12 +27,14 @@ class FiniteDiff(LinearOperator):
     
     '''Finite Difference Operator:
             
-            Computes first-order forward/backward differences 
-                     on 2D, 3D, 4D ImageData
-                     under Neumann/Periodic boundary conditions
+        Computes first-order forward/backward differences 
+                    on 2D, 3D, 4D ImageData
+                    under Neumann/Periodic boundary conditions
 
         Order of the Gradient ( ImageGeometry may contain channels ):
-                            
+
+        .. code:: python          
+
             Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x']
             Grad_order = ['channels', 'direction_y', 'direction_x']
             Grad_order = ['direction_z', 'direction_y', 'direction_x']
@@ -43,7 +45,18 @@ class FiniteDiff(LinearOperator):
 
         
     def __init__(self, gm_domain, gm_range=None, direction=0, bnd_cond = 'Neumann'):
+        '''creator
 
+        :param gm_domain: domain of the operator
+        :type gm_domain: :code:`AcquisitionGeometry` or :code:`ImageGeometry`
+        :param gm_range: optional range of the operator
+        :type gm_range: :code:`AcquisitionGeometry` or :code:`ImageGeometry`, optional
+        :param direction: optional axis in the input :code:`DataContainer` along which to calculate the finite differences, default 0
+        :type direction: int, optional, default 0
+        :param bnd_cond: boundary condition, either :code:`Neumann` or :code:`Periodic`.
+        :type bnd_cond: str, default :code:`Neumann`
+        
+        '''
         super(FiniteDiff, self).__init__() 
 
         self.gm_domain = gm_domain
diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
index 8e07802..2ff0b20 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
@@ -35,7 +35,9 @@ CORRELATION_SPACECHANNEL = "SpaceChannels"
 
 class Gradient(LinearOperator):
 
-    """This is a class to compute the first-order forward/backward differences on ImageData
+
+    r'''Gradient Operator: Computes first-order forward/backward differences on 
+        2D, 3D, 4D ImageData under Neumann/Periodic boundary conditions
     
     :param gm_domain: Set up the domain of the function
     :type gm_domain: `ImageGeometry`
@@ -50,17 +52,17 @@ class Gradient(LinearOperator):
           'Space' or 'SpaceChannels', defaults to 'Space'
         * *backend* (``str``) --
           'c' or 'numpy', defaults to 'c' if correlation is 'SpaceChannels' or channels = 1
-    """
+   
+                                                        
+        Example (2D): 
 
-    r'''Gradient Operator: .. math:: \nabla : X -> Y           
-            
-            Computes first-order forward/backward differences 
-                     on 2D, 3D, 4D ImageData
-                     under Neumann/Periodic boundary conditions
-                                                             
-                Example (2D): u\in X, \nabla(u) = [\partial_{y} u, \partial_{x} u]
-                              u^{*}\in Y, \nabla^{*}(u^{*}) = \partial_{y} v1 + \partial_{x} v2
+        .. math::
         
+          \nabla : X -> Y \\
+          u\in X, \nabla(u) = [\partial_{y} u, \partial_{x} u] \\
+          u^{*}\in Y, \nabla^{*}(u^{*}) = \partial_{y} v1 + \partial_{x} v2
+            
+
     '''
 
     #kept here for backwards compatability
@@ -126,7 +128,15 @@ class Gradient(LinearOperator):
 class Gradient_numpy(LinearOperator):
     
     def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs):
-        
+        '''creator
+        
+        :param gm_domain: domain of the operator
+        :type gm_domain: :code:`AcquisitionGeometry` or :code:`ImageGeometry`
+        :param bnd_cond: boundary condition, either :code:`Neumann` or :code:`Periodic`.
+        :type bnd_cond: str, optional, default :code:`Neumann`
+        :param correlation: optional, :code:`SpaceChannel` or :code:`Space`
+        :type correlation: str, optional, default :code:`Space`
+        '''
         super(Gradient_numpy, self).__init__() 
                 
         self.gm_domain = gm_domain # Domain of Grad Operator
diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
index f4d97b8..fb09819 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
@@ -40,7 +40,15 @@ class LinearOperator(Operator):
     

     @staticmethod

     def PowerMethod(operator, iterations, x_init=None):

-        '''Power method to calculate iteratively the Lipschitz constant'''

+        '''Power method to calculate iteratively the Lipschitz constant

+        

+        :param operator: input operator

+        :type operator: :code:`LinearOperator`

+        :param iterations: number of iterations to run

+        :type iteration: int

+        :param x_init: starting point for the iteration in the operator domain

+        :returns: tuple with: L, list of L at each iteration, the data the iteration worked on.

+        '''

         

         # Initialise random

         if x_init is None:

@@ -73,11 +81,11 @@ class LinearOperator(Operator):
 

     @staticmethod

     def dot_test(operator, domain_init=None, range_init=None, verbose=False):

-        '''Does a dot linearity test on the operator

+        r'''Does a dot linearity test on the operator

         

         Evaluates if the following equivalence holds

         

-        :math: ..

+        .. math::

         

           Ax\times y = y \times A^Tx

         

diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py
index 7d18ea1..a84ea94 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py
@@ -30,6 +30,10 @@ class LinearOperatorMatrix(LinearOperator):
     '''Matrix wrapped into a LinearOperator'''
     
     def __init__(self,A):
+        '''creator
+
+        :param A: numpy ndarray representing a matrix
+        '''
         self.A = A
         M_A, N_A = self.A.shape
         self.gm_domain = VectorGeometry(N_A)
diff --git a/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py
index c5db47d..d1ad07c 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py
@@ -28,9 +28,8 @@ class ScaledOperator(object):
     of the result of direct and adjoint of the operator with the scalar.
     For the rest it behaves like the operator it holds.
     
-    Args:
-       :param operator (Operator): a Operator or LinearOperator
-       :param scalar (Number): a scalar multiplier
+    :param operator: a Operator or LinearOperator
+    :param scalar: a scalar multiplier
     
     Example:
        The scaled operator behaves like the following:
@@ -47,18 +46,25 @@ class ScaledOperator(object):
     '''
     
     def __init__(self, operator, scalar):
+        '''creator
+
+        :param operator: a Operator or LinearOperator
+        :param scalar: a scalar multiplier
+        :type scalar: float'''
         super(ScaledOperator, self).__init__()
         if not isinstance (scalar, Number):
             raise TypeError('expected scalar: got {}'.format(type(scalar)))
         self.scalar = scalar
         self.operator = operator
     def direct(self, x, out=None):
+        '''direct method'''
         if out is None:
             return self.scalar * self.operator.direct(x, out=out)
         else:
             self.operator.direct(x, out=out)
             out *= self.scalar
     def adjoint(self, x, out=None):
+        '''adjoint method'''
         if self.operator.is_linear():
             if out is None:
                 return self.scalar * self.operator.adjoint(x, out=out)
@@ -68,11 +74,17 @@ class ScaledOperator(object):
         else:
             raise TypeError('Operator is not linear')
     def norm(self, **kwargs):
+        '''norm of the operator'''
         return numpy.abs(self.scalar) * self.operator.norm(**kwargs)
     def range_geometry(self):
+        '''range of the operator'''
         return self.operator.range_geometry()
     def domain_geometry(self):
+        '''domain of the operator'''
         return self.operator.domain_geometry()
     def is_linear(self):
+        '''returns whether the operator is linear
+        
+        :returns: boolean '''
         return self.operator.is_linear()
 
diff --git a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
index c85abfa..d82c5c0 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
@@ -30,27 +30,35 @@ class SymmetrizedGradient(LinearOperator):
     
     r'''Symmetrized Gradient Operator:  E: V -> W
         
-            V : range of the Gradient Operator
-            W : range of the Symmetrized Gradient          
+        V : range of the Gradient Operator
+        W : range of the Symmetrized Gradient          
+
+        Example (2D): 
+        
+        .. math::
+            v = (v1, v2) \\
         
-            Example (2D): 
-            .. math::
-               v = (v1, v2), 
+            Ev = 0.5 * ( \nabla\cdot v + (\nabla\cdot c)^{T} ) \\
             
-                           Ev = 0.5 * ( \nabla\cdot v + (\nabla\cdot c)^{T} )
-                           
-                           \begin{matrix} 
-                               \partial_{y} v1 & 0.5 * (\partial_{x} v1 + \partial_{y} v2) \\
-                               0.5 * (\partial_{x} v1 + \partial_{y} v2) & \partial_{x} v2 
-                           \end{matrix}
-              |                                                      
+            \begin{matrix} 
+                \partial_{y} v1 & 0.5 * (\partial_{x} v1 + \partial_{y} v2) \\
+                0.5 * (\partial_{x} v1 + \partial_{y} v2) & \partial_{x} v2 
+            \end{matrix}
+                                                                  
     '''
     
     CORRELATION_SPACE = "Space"
     CORRELATION_SPACECHANNEL = "SpaceChannels"
     
     def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs):
+        '''creator
         
+        :param gm_domain: domain of the operator
+        :param bnd_cond: boundary condition, either :code:`Neumann` or :code:`Periodic`.
+        :type bnd_cond: str, optional, default :code:`Neumann`
+        :param correlation: :code:`SpaceChannel` or :code:`Channel`
+        :type correlation: str, optional, default :code:`Channel`
+        '''
         super(SymmetrizedGradient, self).__init__() 
                 
         self.gm_domain = gm_domain
diff --git a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py
index c37e15e..f677dc2 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py
@@ -27,19 +27,19 @@ from ccpi.optimisation.operators import LinearOperator
 
 class ZeroOperator(LinearOperator):
     
-    r'''ZeroOperator:  O: X -> Y,  maps any element of x\in X into the zero element in Y
-                       O(x) = O_{Y}
+    r'''ZeroOperator:  O: X -> Y,  maps any element of :math:`x\in X` into the zero element :math:`\in Y,  O(x) = O_{Y}`
                        
-                       X : gm_domain
-                       Y : gm_range ( Default: Y = X )
-                       
-                       
-                       Note: 
-                       .. math::
+        :param gm_domain: domain of the operator 
+        :param gm_range: range of the operator, default: same as domain
+        
+        
+        Note: 
+        
+        .. math::
 
-                              O^{*}: Y^{*} -> X^{*} (Adjoint)
-                       
-                              < O(x), y > = < x, O^{*}(y) >
+                O^{*}: Y^{*} -> X^{*} \text{(Adjoint)}
+        
+                < O(x), y > = < x, O^{*}(y) >
                        
      '''
     
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
 ============