From f959a1c7f903fb31b40105f48701aadce2bd7b4c Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Mon, 6 Jan 2020 16:51:02 +0000 Subject: v19.10 docs (#467) updated docstrings and documentation --- README.md | 4 +- .../ccpi/optimisation/algorithms/Algorithm.py | 26 +- .../Python/ccpi/optimisation/algorithms/CGLS.py | 9 +- .../Python/ccpi/optimisation/algorithms/FISTA.py | 20 +- .../Python/ccpi/optimisation/algorithms/PDHG.py | 44 +- .../Python/ccpi/optimisation/algorithms/SIRT.py | 45 +- .../ccpi/optimisation/functions/BlockFunction.py | 8 +- .../functions/FunctionOperatorComposition.py | 17 +- .../ccpi/optimisation/functions/IndicatorBox.py | 17 +- .../Python/ccpi/optimisation/functions/L1Norm.py | 14 +- .../ccpi/optimisation/functions/L2NormSquared.py | 10 +- .../ccpi/optimisation/functions/MixedL21Norm.py | 14 +- .../operators/FiniteDifferenceOperator.py | 21 +- .../optimisation/operators/GradientOperator.py | 32 +- .../ccpi/optimisation/operators/LinearOperator.py | 14 +- .../optimisation/operators/LinearOperatorMatrix.py | 4 + .../ccpi/optimisation/operators/ScaledOperator.py | 18 +- .../operators/SymmetrizedGradientOperator.py | 32 +- .../ccpi/optimisation/operators/ZeroOperator.py | 22 +- docs/source/astra.rst | 3 - docs/source/conf.py | 6 +- docs/source/contrib.rst | 1 - docs/source/framework.rst | 474 ++++++++++++++++++++- docs/source/images/cone.png | Bin 0 -> 127928 bytes docs/source/images/fan.png | Bin 0 -> 86375 bytes docs/source/images/fan_data.png | Bin 0 -> 87766 bytes docs/source/images/fan_geometry.png | Bin 0 -> 136263 bytes docs/source/images/parallel.png | Bin 0 -> 29796 bytes docs/source/images/parallel3d.png | Bin 0 -> 375145 bytes docs/source/images/parallel3d_data.png | Bin 0 -> 371872 bytes docs/source/images/parallel3d_geometry.png | Bin 0 -> 423629 bytes docs/source/images/parallel_data.png | Bin 0 -> 21843 bytes docs/source/images/parallel_geometry.png | Bin 0 -> 79825 bytes docs/source/index.rst | 41 +- docs/source/io.rst | 29 +- docs/source/optimisation.rst | 282 +++++++++++- docs/source/plugins.rst | 2 - 37 files changed, 1022 insertions(+), 187 deletions(-) create mode 100755 docs/source/images/cone.png create mode 100755 docs/source/images/fan.png create mode 100755 docs/source/images/fan_data.png create mode 100755 docs/source/images/fan_geometry.png create mode 100755 docs/source/images/parallel.png create mode 100755 docs/source/images/parallel3d.png create mode 100755 docs/source/images/parallel3d_data.png create mode 100755 docs/source/images/parallel3d_geometry.png create mode 100755 docs/source/images/parallel_data.png create mode 100755 docs/source/images/parallel_geometry.png 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 ` 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 ` 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 `_ 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 `_: + +.. 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 `_. + +.. autoclass:: ccpi.processors.CenterOfRotationFinder + :members: + :private-members: + :special-members: + :ref:`Return Home ` + +.. _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 Binary files /dev/null and b/docs/source/images/cone.png differ diff --git a/docs/source/images/fan.png b/docs/source/images/fan.png new file mode 100755 index 0000000..4f20da4 Binary files /dev/null and b/docs/source/images/fan.png differ diff --git a/docs/source/images/fan_data.png b/docs/source/images/fan_data.png new file mode 100755 index 0000000..cc7a470 Binary files /dev/null and b/docs/source/images/fan_data.png differ diff --git a/docs/source/images/fan_geometry.png b/docs/source/images/fan_geometry.png new file mode 100755 index 0000000..c1b6a50 Binary files /dev/null and b/docs/source/images/fan_geometry.png differ diff --git a/docs/source/images/parallel.png b/docs/source/images/parallel.png new file mode 100755 index 0000000..a58f79e Binary files /dev/null and b/docs/source/images/parallel.png differ diff --git a/docs/source/images/parallel3d.png b/docs/source/images/parallel3d.png new file mode 100755 index 0000000..f5dc76f Binary files /dev/null and b/docs/source/images/parallel3d.png differ diff --git a/docs/source/images/parallel3d_data.png b/docs/source/images/parallel3d_data.png new file mode 100755 index 0000000..2b5536a Binary files /dev/null and b/docs/source/images/parallel3d_data.png differ diff --git a/docs/source/images/parallel3d_geometry.png b/docs/source/images/parallel3d_geometry.png new file mode 100755 index 0000000..fdcff6f Binary files /dev/null and b/docs/source/images/parallel3d_geometry.png differ diff --git a/docs/source/images/parallel_data.png b/docs/source/images/parallel_data.png new file mode 100755 index 0000000..7adea39 Binary files /dev/null and b/docs/source/images/parallel_data.png differ diff --git a/docs/source/images/parallel_geometry.png b/docs/source/images/parallel_geometry.png new file mode 100755 index 0000000..2880b55 Binary files /dev/null and b/docs/source/images/parallel_geometry.png differ 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 ` + +.. _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 ============ -- cgit v1.2.3