From 6d609d54f828882ec46e11af4d3e09fc83a20535 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 1 Mar 2019 17:42:53 +0000 Subject: Geometry allocation (#211) * initial revision * Removed class members of Algorithm class added update_objective * initial version. Fix inline __idiv__ * First implementation of CompositeOperator/DataContainer * removed __getitem__ added get_item added shape * added CGLS * working unit test, initial tomography test * added reverse multiplication of operator with number * added operators directory * fixed typo * added unittest for CompositeDataContainer * fix TomoIdentity with scalar * check numerical types from numpy * add default stop criterion and run method * add run method * first working implementation of CGLS with CompositeOperator/DataContainer notice problem with _rmul_ and _mul_ methods precedence with numpy. * removed line endings * removed dos line ending * Added allocate to ImageGeometry and AcquisitionGeometry * remove composite operator * Delete Algorithms.py * remove setup.py from PR * readded setup.py * added newline at the end of the file * added newline at EOF in setup.py --- Wrappers/Python/ccpi/framework.py | 25 ++- .../Python/ccpi/optimisation/algorithms/FBPD.py | 172 ++++++++++----------- .../ccpi/optimisation/algorithms/__init__.py | 58 +++---- Wrappers/Python/test/test_DataContainer.py | 43 ++++++ 4 files changed, 176 insertions(+), 122 deletions(-) mode change 100755 => 100644 Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py mode change 100755 => 100644 Wrappers/Python/ccpi/optimisation/algorithms/__init__.py (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py index 0c23628..71a9f3f 100644 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework.py @@ -111,8 +111,12 @@ class ImageGeometry(object): repres += "voxel_size : x{0},y{1},z{2}\n".format(self.voxel_size_x, self.voxel_size_y, self.voxel_size_z) repres += "center : x{0},y{1},z{2}\n".format(self.center_x, self.center_y, self.center_z) return repres - - + def allocate(self, value=0, dimension_labels=None): + '''allocates an ImageData according to the size expressed in the instance''' + out = ImageData(geometry=self, dimension_labels=dimension_labels) + if value != 0: + out += value + return out class AcquisitionGeometry(object): def __init__(self, @@ -192,9 +196,12 @@ class AcquisitionGeometry(object): repres += "distance center-detector: {0}\n".format(self.dist_source_center) repres += "number of channels: {0}\n".format(self.channels) return repres - - - + def allocate(self, value=0, dimension_labels=None): + '''allocates an AcquisitionData according to the size expressed in the instance''' + out = AcquisitionData(geometry=self, dimension_labels=dimension_labels) + if value != 0: + out += value + return out class DataContainer(object): '''Generic class to hold data @@ -741,6 +748,7 @@ class DataContainer(object): return numpy.sqrt(self.squared_norm()) + class ImageData(DataContainer): '''DataContainer for holding 2D or 3D DataContainer''' def __init__(self, @@ -903,8 +911,11 @@ class AcquisitionData(DataContainer): elif dim == 'horizontal': shape.append(horiz) if len(shape) != len(dimension_labels): - raise ValueError('Missing {0} axes'.format( - len(dimension_labels) - len(shape))) + raise ValueError('Missing {0} axes.\nExpected{1} got {2}}'\ + .format( + len(dimension_labels) - len(shape), + dimension_labels, shape) + ) shape = tuple(shape) array = numpy.zeros( shape , dtype=numpy.float32) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py b/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py old mode 100755 new mode 100644 index 322e9eb..798fb61 --- a/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py @@ -1,86 +1,86 @@ -# -*- coding: utf-8 -*- -# This work is part of the Core Imaging Library developed by -# Visual Analytics and Imaging System Group of the Science Technology -# Facilities Council, STFC - -# Copyright 2019 Edoardo Pasca - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -""" -Created on Thu Feb 21 11:09:03 2019 - -@author: ofn77899 -""" - -from ccpi.optimisation.algorithms import Algorithm -from ccpi.optimisation.funcs import ZeroFun - -class FBPD(Algorithm): - '''FBPD Algorithm - - Parameters: - x_init: initial guess - f: constraint - g: data fidelity - h: regularizer - opt: additional algorithm - ''' - constraint = None - data_fidelity = None - regulariser = None - def __init__(self, **kwargs): - pass - def set_up(self, x_init, operator=None, constraint=None, data_fidelity=None,\ - regulariser=None, opt=None): - - # default inputs - if constraint is None: - self.constraint = ZeroFun() - else: - self.constraint = constraint - if data_fidelity is None: - data_fidelity = ZeroFun() - else: - self.data_fidelity = data_fidelity - if regulariser is None: - self.regulariser = ZeroFun() - else: - self.regulariser = regulariser - - # algorithmic parameters - - - # step-sizes - self.tau = 2 / (self.data_fidelity.L + 2) - self.sigma = (1/self.tau - self.data_fidelity.L/2) / self.regulariser.L - - self.inv_sigma = 1/self.sigma - - # initialization - self.x = x_init - self.y = operator.direct(self.x) - - - def update(self): - - # primal forward-backward step - x_old = self.x - self.x = self.x - self.tau * ( self.data_fidelity.grad(self.x) + self.operator.adjoint(self.y) ) - self.x = self.constraint.prox(self.x, self.tau); - - # dual forward-backward step - self.y = self.y + self.sigma * self.operator.direct(2*self.x - x_old); - self.y = self.y - self.sigma * self.regulariser.prox(self.inv_sigma*self.y, self.inv_sigma); - - # time and criterion - self.loss = self.constraint(self.x) + self.data_fidelity(self.x) + self.regulariser(self.operator.direct(self.x)) +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2019 Edoardo Pasca + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +""" +Created on Thu Feb 21 11:09:03 2019 + +@author: ofn77899 +""" + +from ccpi.optimisation.algorithms import Algorithm +from ccpi.optimisation.funcs import ZeroFun + +class FBPD(Algorithm): + '''FBPD Algorithm + + Parameters: + x_init: initial guess + f: constraint + g: data fidelity + h: regularizer + opt: additional algorithm + ''' + constraint = None + data_fidelity = None + regulariser = None + def __init__(self, **kwargs): + pass + def set_up(self, x_init, operator=None, constraint=None, data_fidelity=None,\ + regulariser=None, opt=None): + + # default inputs + if constraint is None: + self.constraint = ZeroFun() + else: + self.constraint = constraint + if data_fidelity is None: + data_fidelity = ZeroFun() + else: + self.data_fidelity = data_fidelity + if regulariser is None: + self.regulariser = ZeroFun() + else: + self.regulariser = regulariser + + # algorithmic parameters + + + # step-sizes + self.tau = 2 / (self.data_fidelity.L + 2) + self.sigma = (1/self.tau - self.data_fidelity.L/2) / self.regulariser.L + + self.inv_sigma = 1/self.sigma + + # initialization + self.x = x_init + self.y = operator.direct(self.x) + + + def update(self): + + # primal forward-backward step + x_old = self.x + self.x = self.x - self.tau * ( self.data_fidelity.grad(self.x) + self.operator.adjoint(self.y) ) + self.x = self.constraint.prox(self.x, self.tau); + + # dual forward-backward step + self.y = self.y + self.sigma * self.operator.direct(2*self.x - x_old); + self.y = self.y - self.sigma * self.regulariser.prox(self.inv_sigma*self.y, self.inv_sigma); + + # time and criterion + self.loss = self.constraint(self.x) + self.data_fidelity(self.x) + self.regulariser(self.operator.direct(self.x)) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py old mode 100755 new mode 100644 index 52fe6d7..903bc30 --- a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py @@ -1,29 +1,29 @@ -# -*- coding: utf-8 -*- -# This work is part of the Core Imaging Library developed by -# Visual Analytics and Imaging System Group of the Science Technology -# Facilities Council, STFC - -# Copyright 2019 Edoardo Pasca - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -""" -Created on Thu Feb 21 11:03:13 2019 - -@author: ofn77899 -""" - -from .Algorithm import Algorithm -from .CGLS import CGLS -from .GradientDescent import GradientDescent -from .FISTA import FISTA -from .FBPD import FBPD \ No newline at end of file +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2019 Edoardo Pasca + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +""" +Created on Thu Feb 21 11:03:13 2019 + +@author: ofn77899 +""" + +from .Algorithm import Algorithm +from .CGLS import CGLS +from .GradientDescent import GradientDescent +from .FISTA import FISTA +from .FBPD import FBPD diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py index 05f3fe8..3def054 100755 --- a/Wrappers/Python/test/test_DataContainer.py +++ b/Wrappers/Python/test/test_DataContainer.py @@ -457,6 +457,49 @@ class TestDataContainer(unittest.TestCase): pixel_num_h=5, channels=2) sino = AcquisitionData(geometry=sgeometry) self.assertEqual(sino.shape, (2, 10, 3, 5)) + def test_ImageGeometry_allocate(self): + vgeometry = ImageGeometry(voxel_num_x=4, voxel_num_y=3, channels=2) + image = vgeometry.allocate() + self.assertEqual(0,image.as_array()[0][0][0]) + image = vgeometry.allocate(1) + self.assertEqual(1,image.as_array()[0][0][0]) + default_order = ['channel' , 'horizontal_y' , 'horizontal_x'] + self.assertEqual(default_order[0], image.dimension_labels[0]) + self.assertEqual(default_order[1], image.dimension_labels[1]) + self.assertEqual(default_order[2], image.dimension_labels[2]) + order = [ 'horizontal_x' , 'horizontal_y', 'channel' ] + image = vgeometry.allocate(0,dimension_labels=order) + self.assertEqual(order[0], image.dimension_labels[0]) + self.assertEqual(order[1], image.dimension_labels[1]) + self.assertEqual(order[2], image.dimension_labels[2]) + def test_AcquisitionGeometry_allocate(self): + ageometry = AcquisitionGeometry(dimension=2, angles=numpy.linspace(0, 180, num=10), + geom_type='parallel', pixel_num_v=3, + pixel_num_h=5, channels=2) + sino = ageometry.allocate() + shape = sino.shape + print ("shape", shape) + self.assertEqual(0,sino.as_array()[0][0][0][0]) + self.assertEqual(0,sino.as_array()[shape[0]-1][shape[1]-1][shape[2]-1][shape[3]-1]) + + sino = ageometry.allocate(1) + self.assertEqual(1,sino.as_array()[0][0][0][0]) + self.assertEqual(1,sino.as_array()[shape[0]-1][shape[1]-1][shape[2]-1][shape[3]-1]) + print (sino.dimension_labels, sino.shape, ageometry) + + default_order = ['channel' , ' angle' , + 'vertical' , 'horizontal'] + self.assertEqual(default_order[0], sino.dimension_labels[0]) + self.assertEqual(default_order[1], sino.dimension_labels[1]) + self.assertEqual(default_order[2], sino.dimension_labels[2]) + self.assertEqual(default_order[3], sino.dimension_labels[3]) + order = ['vertical' , 'horizontal', 'channel' , 'angle' ] + sino = ageometry.allocate(0,dimension_labels=order) + print (sino.dimension_labels, sino.shape, ageometry) + self.assertEqual(order[0], sino.dimension_labels[0]) + self.assertEqual(order[1], sino.dimension_labels[1]) + self.assertEqual(order[2], sino.dimension_labels[2]) + self.assertEqual(order[2], sino.dimension_labels[2]) def assertNumpyArrayEqual(self, first, second): res = True -- cgit v1.2.3