From fc4d8d28da1e8148f781a494338f6fcef48ef1fe Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 26 Apr 2019 15:34:46 +0100 Subject: Block Framework (#23) * fix for new block framework * refactoring closes #22 * updated setup.py --- .../ccpi/plugins/operators/CCPiProjectorSimple.py | 160 +++++++++++ Wrappers/Python/ccpi/plugins/operators/__init__.py | 1 + Wrappers/Python/ccpi/plugins/ops.py | 151 ---------- Wrappers/Python/ccpi/plugins/processors.py | 303 -------------------- .../Python/ccpi/plugins/processors/__init__.py | 4 + .../Python/ccpi/plugins/processors/processors.py | 309 +++++++++++++++++++++ Wrappers/Python/setup.py | 3 +- 7 files changed, 476 insertions(+), 455 deletions(-) create mode 100755 Wrappers/Python/ccpi/plugins/operators/CCPiProjectorSimple.py create mode 100644 Wrappers/Python/ccpi/plugins/operators/__init__.py delete mode 100755 Wrappers/Python/ccpi/plugins/ops.py delete mode 100755 Wrappers/Python/ccpi/plugins/processors.py create mode 100644 Wrappers/Python/ccpi/plugins/processors/__init__.py create mode 100755 Wrappers/Python/ccpi/plugins/processors/processors.py diff --git a/Wrappers/Python/ccpi/plugins/operators/CCPiProjectorSimple.py b/Wrappers/Python/ccpi/plugins/operators/CCPiProjectorSimple.py new file mode 100755 index 0000000..eedcac3 --- /dev/null +++ b/Wrappers/Python/ccpi/plugins/operators/CCPiProjectorSimple.py @@ -0,0 +1,160 @@ +# -*- 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 2018 Jakob Jorgensen, Daniil Kazantsev and 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. + +import numpy +from ccpi.optimisation.operators import Operator, LinearOperator +from ccpi.framework import ImageData, DataContainer , \ + ImageGeometry, AcquisitionGeometry +from ccpi.plugins.processors import CCPiBackwardProjector, \ + CCPiForwardProjector , setupCCPiGeometries + + + +class CCPiProjectorSimple(Operator): + """ASTRA projector modified to use DataSet and geometry.""" + def __init__(self, geomv, geomp, default=False): + super(CCPiProjectorSimple, self).__init__() + + # Store volume and sinogram geometries. + self.acquisition_geometry = geomp + self.volume_geometry = geomv + + if geomp.geom_type == "cone": + raise TypeError('Can only handle parallel beam') + + # set-up the geometries if compatible + geoms = setupCCPiGeometries(geomv, geomp, 0) + + + vg = ImageGeometry(voxel_num_x=geoms['output_volume_x'], + voxel_num_y=geoms['output_volume_y'], + voxel_num_z=geoms['output_volume_z']) + + pg = AcquisitionGeometry('parallel', + '3D', + geomp.angles, + geoms['n_h'], geomp.pixel_size_h, + geoms['n_v'], geomp.pixel_size_v #2D in 3D is a slice 1 pixel thick + ) + if not default: + # check if geometry is the same (on the voxels) + if not ( vg.voxel_num_x == geomv.voxel_num_x and \ + vg.voxel_num_y == geomv.voxel_num_y and \ + vg.voxel_num_z == geomv.voxel_num_z ): + msg = 'The required volume geometry will not work\nThe following would\n' + msg += vg.__str__() + raise ValueError(msg) + if not (pg.pixel_num_h == geomp.pixel_num_h and \ + pg.pixel_num_v == geomp.pixel_num_v and \ + len( pg.angles ) == len( geomp.angles ) ) : + msg = 'The required acquisition geometry will not work\nThe following would\n' + msg += pg.__str__() + raise ValueError(msg) + + self.fp = CCPiForwardProjector(image_geometry=vg, + acquisition_geometry=pg, + output_axes_order=['angle','vertical','horizontal']) + + self.bp = CCPiBackwardProjector(image_geometry=vg, + acquisition_geometry=pg, + output_axes_order=[ ImageGeometry.HORIZONTAL_X, + ImageGeometry.HORIZONTAL_Y, ImageGeometry.VERTICAL]) + + # Initialise empty for singular value. + self.s1 = None + self.ag = pg + self.vg = vg + + def is_linear(self): + return True + + def direct(self, image_data, out=None): + self.fp.set_input(image_data.subset(dimensions=['horizontal_x','horizontal_y','vertical'])) + if out is None: + out = self.fp.get_output() + #return out.subset(dimensions=[AcquisitionGeometry.ANGLE , AcquisitionGeometry.VERTICAL , + # AcquisitionGeometry.HORIZONTAL]) + return out + else: + out.fill(self.fp.get_output()) + + def adjoint(self, acquisition_data, out=None): + self.bp.set_input(acquisition_data.subset(dimensions=['angle','vertical','horizontal'])) + if out is None: + out = self.bp.get_output() + #return out.subset(dimensions=[ImageGeometry.VERTICAL, ImageGeometry.HORIZONTAL_Y, + # ImageGeometry.HORIZONTAL_X]) + return out + else: + out.fill(self.bp.get_output()) + + #def delete(self): + # astra.data2d.delete(self.proj_id) + + def norm(self): + x0 = self.domain_geometry().allocate(ImageGeometry.RANDOM, + #dimension_labels=['horizontal_x', 'horizontal_y','vertical'] + ) + print (x0.dimension_labels) + + a = LinearOperator.PowerMethod(self,10, x0) + self.s1 = a[0] + return self.s1 + + def size(self): + # Only implemented for 3D + return ( (self.acquisition_geometry.angles.size, \ + self.acquisition_geometry.pixel_num_v, + self.acquisition_geometry.pixel_num_h), \ + (self.volume_geometry.voxel_num_x, \ + self.volume_geometry.voxel_num_y, + self.volume_geometry.voxel_num_z) ) + def create_image_data(self): + x0 = ImageData(geometry = self.volume_geometry, + dimension_labels=self.bp.output_axes_order)#\ + #.subset(['horizontal_x','horizontal_y','vertical']) + x0.fill(numpy.random.randn(*x0.shape)) + return x0 + def domain_geometry(self): + return ImageGeometry( + self.vg.voxel_num_x, + self.vg.voxel_num_y, + self.vg.voxel_num_z, + self.vg.voxel_size_x, + self.vg.voxel_size_y, + self.vg.voxel_size_z, + self.vg.center_x, + self.vg.center_y, + self.vg.center_z, + self.vg.channels + ) + + def range_geometry(self): + return AcquisitionGeometry(self.ag.geom_type, + self.ag.dimension, + self.ag.angles, + self.ag.pixel_num_h, + self.ag.pixel_size_h, + self.ag.pixel_num_v, + self.ag.pixel_size_v, + self.ag.dist_source_center, + self.ag.dist_center_detector, + self.ag.channels, + ) + diff --git a/Wrappers/Python/ccpi/plugins/operators/__init__.py b/Wrappers/Python/ccpi/plugins/operators/__init__.py new file mode 100644 index 0000000..f5ffd0c --- /dev/null +++ b/Wrappers/Python/ccpi/plugins/operators/__init__.py @@ -0,0 +1 @@ +from .CCPiProjectorSimple import CCPiProjectorSimple \ No newline at end of file diff --git a/Wrappers/Python/ccpi/plugins/ops.py b/Wrappers/Python/ccpi/plugins/ops.py deleted file mode 100755 index f2e0dd3..0000000 --- a/Wrappers/Python/ccpi/plugins/ops.py +++ /dev/null @@ -1,151 +0,0 @@ -# -*- 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 2018 Jakob Jorgensen, Daniil Kazantsev and 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. - -import numpy -from ccpi.optimisation.ops import Operator, PowerMethodNonsquare -from ccpi.framework import ImageData, DataContainer , \ - ImageGeometry, AcquisitionGeometry -from ccpi.plugins.processors import CCPiBackwardProjector, \ - CCPiForwardProjector , setupCCPiGeometries - - - -class CCPiProjectorSimple(Operator): - """ASTRA projector modified to use DataSet and geometry.""" - def __init__(self, geomv, geomp, default=False): - super(CCPiProjectorSimple, self).__init__() - - # Store volume and sinogram geometries. - self.acquisition_geometry = geomp - self.volume_geometry = geomv - - if geomp.geom_type == "cone": - raise TypeError('Can only handle parallel beam') - - # set-up the geometries if compatible - geoms = setupCCPiGeometries(geomv.voxel_num_x,geomv.voxel_num_y, - geomv.voxel_num_z, geomp.angles, 0) - - - vg = ImageGeometry(voxel_num_x=geoms['output_volume_x'], - voxel_num_y=geoms['output_volume_y'], - voxel_num_z=geoms['output_volume_z']) - - pg = AcquisitionGeometry('parallel', - '3D', - geomp.angles, - geoms['n_h'], geomp.pixel_size_h, - geoms['n_v'], geomp.pixel_size_v #2D in 3D is a slice 1 pixel thick - ) - if not default: - # check if geometry is the same (on the voxels) - if not ( vg.voxel_num_x == geomv.voxel_num_x and \ - vg.voxel_num_y == geomv.voxel_num_y and \ - vg.voxel_num_z == geomv.voxel_num_z ): - msg = 'The required volume geometry will not work\nThe following would\n' - msg += vg.__str__() - raise ValueError(msg) - if not (pg.pixel_num_h == geomp.pixel_num_h and \ - pg.pixel_num_v == geomp.pixel_num_v and \ - len( pg.angles ) == len( geomp.angles ) ) : - msg = 'The required acquisition geometry will not work\nThe following would\n' - msg += pg.__str__() - raise ValueError(msg) - - self.fp = CCPiForwardProjector(image_geometry=vg, - acquisition_geometry=pg, - output_axes_order=['angle','vertical','horizontal']) - - self.bp = CCPiBackwardProjector(image_geometry=vg, - acquisition_geometry=pg, - output_axes_order=['horizontal_x','horizontal_y','vertical']) - - # Initialise empty for singular value. - self.s1 = None - self.ag = pg - self.vg = vg - - def is_linear(self): - return True - - def direct(self, image_data, out=None): - self.fp.set_input(image_data) - if out is None: - out = self.fp.get_output() - return out - else: - out.fill(self.fp.get_output()) - - def adjoint(self, acquisition_data, out=None): - self.bp.set_input(acquisition_data) - if out is None: - out = self.bp.get_output() - return out - else: - out.fill(self.bp.get_output()) - - #def delete(self): - # astra.data2d.delete(self.proj_id) - - def get_max_sing_val(self): - a = PowerMethodNonsquare(self,10) - self.s1 = a[0] - return self.s1 - - def size(self): - # Only implemented for 3D - return ( (self.acquisition_geometry.angles.size, \ - self.acquisition_geometry.pixel_num_v, - self.acquisition_geometry.pixel_num_h), \ - (self.volume_geometry.voxel_num_x, \ - self.volume_geometry.voxel_num_y, - self.volume_geometry.voxel_num_z) ) - def create_image_data(self): - x0 = ImageData(geometry = self.volume_geometry, - dimension_labels=self.bp.output_axes_order)#\ - #.subset(['horizontal_x','horizontal_y','vertical']) - x0.fill(numpy.random.randn(*x0.shape)) - return x0 - def domain_geometry(self): - return ImageGeometry( - self.vg.voxel_num_x, - self.vg.voxel_num_y, - self.vg.voxel_num_z, - self.vg.voxel_size_x, - self.vg.voxel_size_y, - self.vg.voxel_size_z, - self.vg.center_x, - self.vg.center_y, - self.vg.center_z, - self.vg.channels, - ['horizontal_x','horizontal_y','vertical'] ) - - def range_geometry(self): - return AcquisitionGeometry(self.ag.geom_type, - self.ag.dimension, - self.ag.angles, - self.ag.pixel_num_h, - self.ag.pixel_size_h, - self.ag.pixel_num_v, - self.ag.pixel_size_v, - self.ag.dist_source_center, - self.ag.dist_center_detector, - self.ag.channels, - ['angle','vertical','horizontal']) - diff --git a/Wrappers/Python/ccpi/plugins/processors.py b/Wrappers/Python/ccpi/plugins/processors.py deleted file mode 100755 index d31c03b..0000000 --- a/Wrappers/Python/ccpi/plugins/processors.py +++ /dev/null @@ -1,303 +0,0 @@ -# -*- 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 2018 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 - -from ccpi.framework import DataProcessor, AcquisitionData,\ - AcquisitionGeometry, ImageGeometry, ImageData -from ccpi.reconstruction.parallelbeam import alg as pbalg -import numpy - -def setupCCPiGeometries(voxel_num_x, voxel_num_y, voxel_num_z, angles, counter): - - vg = ImageGeometry(voxel_num_x=voxel_num_x,voxel_num_y=voxel_num_y, voxel_num_z=voxel_num_z) - Phantom_ccpi = ImageData(geometry=vg, - dimension_labels=['horizontal_x','horizontal_y','vertical']) - #.subset(['horizontal_x','horizontal_y','vertical']) - # ask the ccpi code what dimensions it would like - - voxel_per_pixel = 1 - geoms = pbalg.pb_setup_geometry_from_image(Phantom_ccpi.as_array(), - angles, - voxel_per_pixel ) - - pg = AcquisitionGeometry('parallel', - '3D', - angles, - geoms['n_h'], 1.0, - geoms['n_v'], 1.0 #2D in 3D is a slice 1 pixel thick - ) - - center_of_rotation = Phantom_ccpi.get_dimension_size('horizontal_x') / 2 - ad = AcquisitionData(geometry=pg,dimension_labels=['angle','vertical','horizontal']) - geoms_i = pbalg.pb_setup_geometry_from_acquisition(ad.as_array(), - angles, - center_of_rotation, - voxel_per_pixel ) - - counter+=1 - - if counter < 4: - if (not ( geoms_i == geoms )): - print ("not equal and {0}".format(counter)) - X = max(geoms['output_volume_x'], geoms_i['output_volume_x']) - Y = max(geoms['output_volume_y'], geoms_i['output_volume_y']) - Z = max(geoms['output_volume_z'], geoms_i['output_volume_z']) - return setupCCPiGeometries(X,Y,Z,angles, counter) - else: - return geoms - else: - return geoms_i - - -class CCPiForwardProjector(DataProcessor): - '''Normalization based on flat and dark - - This processor read in a AcquisitionData and normalises it based on - the instrument reading with and without incident photons or neutrons. - - Input: AcquisitionData - Parameter: 2D projection with flat field (or stack) - 2D projection with dark field (or stack) - Output: AcquisitionDataSetn - ''' - - def __init__(self, - image_geometry = None, - acquisition_geometry = None, - output_axes_order = None): - if output_axes_order is None: - # default ccpi projector image storing order - output_axes_order = ['angle','vertical','horizontal'] - - kwargs = { - 'image_geometry' : image_geometry, - 'acquisition_geometry' : acquisition_geometry, - 'output_axes_order' : output_axes_order, - 'default_image_axes_order' : ['horizontal_x','horizontal_y','vertical'], - 'default_acquisition_axes_order' : ['angle','vertical','horizontal'] - } - - super(CCPiForwardProjector, self).__init__(**kwargs) - - def check_input(self, dataset): - if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2: - # sort in the order that this projector needs it - return True - else: - raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ - .format(dataset.number_of_dimensions)) - - def process(self): - - volume = self.get_input() - volume_axes = volume.get_data_axes_order(new_order=self.default_image_axes_order) - if not volume_axes == [0,1,2]: - volume.array = numpy.transpose(volume.array, volume_axes) - pixel_per_voxel = 1 # should be estimated from image_geometry and - # acquisition_geometry - if self.acquisition_geometry.geom_type == 'parallel': - - pixels = pbalg.pb_forward_project(volume.as_array(), - self.acquisition_geometry.angles, - pixel_per_voxel) - out = AcquisitionData(geometry=self.acquisition_geometry, - label_dimensions=self.default_acquisition_axes_order) - out.fill(pixels) - out_axes = out.get_data_axes_order(new_order=self.output_axes_order) - if not out_axes == [0,1,2]: - out.array = numpy.transpose(out.array, out_axes) - return out - else: - raise ValueError('Cannot process cone beam') - -class CCPiBackwardProjector(DataProcessor): - '''Backward projector - - This processor reads in a AcquisitionData and performs a backward projection, - i.e. project to reconstruction space. - Notice that it assumes that the center of rotation is in the middle - of the horizontal axis: in case when that's not the case it can be chained - with the AcquisitionDataPadder. - - Input: AcquisitionData - Parameter: 2D projection with flat field (or stack) - 2D projection with dark field (or stack) - Output: AcquisitionDataSetn - ''' - - def __init__(self, - image_geometry = None, - acquisition_geometry = None, - output_axes_order=None): - if output_axes_order is None: - # default ccpi projector image storing order - output_axes_order = ['horizontal_x','horizontal_y','vertical'] - kwargs = { - 'image_geometry' : image_geometry, - 'acquisition_geometry' : acquisition_geometry, - 'output_axes_order' : output_axes_order, - 'default_image_axes_order' : ['horizontal_x','horizontal_y','vertical'], - 'default_acquisition_axes_order' : ['angle','vertical','horizontal'] - } - - super(CCPiBackwardProjector, self).__init__(**kwargs) - - def check_input(self, dataset): - if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2: - - return True - else: - raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ - .format(dataset.number_of_dimensions)) - - def process(self): - projections = self.get_input() - projections_axes = projections.get_data_axes_order(new_order=self.default_acquisition_axes_order) - if not projections_axes == [0,1,2]: - projections.array = numpy.transpose(projections.array, projections_axes) - - pixel_per_voxel = 1 # should be estimated from image_geometry and acquisition_geometry - image_geometry = ImageGeometry(voxel_num_x = self.acquisition_geometry.pixel_num_h, - voxel_num_y = self.acquisition_geometry.pixel_num_h, - voxel_num_z = self.acquisition_geometry.pixel_num_v) - # input centered/padded acquisitiondata - center_of_rotation = projections.get_dimension_size('horizontal') / 2 - - if self.acquisition_geometry.geom_type == 'parallel': - back = pbalg.pb_backward_project( - projections.as_array(), - self.acquisition_geometry.angles, - center_of_rotation, - pixel_per_voxel - ) - out = ImageData(geometry=self.image_geometry, - dimension_labels=self.default_image_axes_order) - - out_axes = out.get_data_axes_order(new_order=self.output_axes_order) - if not out_axes == [0,1,2]: - back = numpy.transpose(back, out_axes) - out.fill(back) - - return out - - else: - raise ValueError('Cannot process cone beam') - -class AcquisitionDataPadder(DataProcessor): - '''Normalization based on flat and dark - - This processor read in a AcquisitionData and normalises it based on - the instrument reading with and without incident photons or neutrons. - - Input: AcquisitionData - Parameter: 2D projection with flat field (or stack) - 2D projection with dark field (or stack) - Output: AcquisitionDataSetn - ''' - - def __init__(self, - center_of_rotation = None, - acquisition_geometry = None, - pad_value = 1e-5): - kwargs = { - 'acquisition_geometry' : acquisition_geometry, - 'center_of_rotation' : center_of_rotation, - 'pad_value' : pad_value - } - - super(AcquisitionDataPadder, self).__init__(**kwargs) - - def check_input(self, dataset): - if self.acquisition_geometry is None: - self.acquisition_geometry = dataset.geometry - if dataset.number_of_dimensions == 3: - return True - else: - raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ - .format(dataset.number_of_dimensions)) - - def process(self): - projections = self.get_input() - w = projections.get_dimension_size('horizontal') - delta = w - 2 * self.center_of_rotation - - padded_width = int ( - numpy.ceil(abs(delta)) + w - ) - delta_pix = padded_width - w - - voxel_per_pixel = 1 - geom = pbalg.pb_setup_geometry_from_acquisition(projections.as_array(), - self.acquisition_geometry.angles, - self.center_of_rotation, - voxel_per_pixel ) - - padded_geometry = self.acquisition_geometry.clone() - - padded_geometry.pixel_num_h = geom['n_h'] - padded_geometry.pixel_num_v = geom['n_v'] - - delta_pix_h = padded_geometry.pixel_num_h - self.acquisition_geometry.pixel_num_h - delta_pix_v = padded_geometry.pixel_num_v - self.acquisition_geometry.pixel_num_v - - if delta_pix_h == 0: - delta_pix_h = delta_pix - padded_geometry.pixel_num_h = padded_width - #initialize a new AcquisitionData with values close to 0 - out = AcquisitionData(geometry=padded_geometry) - out = out + self.pad_value - - - #pad in the horizontal-vertical plane -> slice on angles - if delta > 0: - #pad left of middle - command = "out.array[" - for i in range(out.number_of_dimensions): - if out.dimension_labels[i] == 'horizontal': - value = '{0}:{1}'.format(delta_pix_h, delta_pix_h+w) - command = command + str(value) - else: - if out.dimension_labels[i] == 'vertical' : - value = '{0}:'.format(delta_pix_v) - command = command + str(value) - else: - command = command + ":" - if i < out.number_of_dimensions -1: - command = command + ',' - command = command + '] = projections.array' - #print (command) - else: - #pad right of middle - command = "out.array[" - for i in range(out.number_of_dimensions): - if out.dimension_labels[i] == 'horizontal': - value = '{0}:{1}'.format(0, w) - command = command + str(value) - else: - if out.dimension_labels[i] == 'vertical' : - value = '{0}:'.format(delta_pix_v) - command = command + str(value) - else: - command = command + ":" - if i < out.number_of_dimensions -1: - command = command + ',' - command = command + '] = projections.array' - #print (command) - #cleaned = eval(command) - exec(command) - return out \ No newline at end of file diff --git a/Wrappers/Python/ccpi/plugins/processors/__init__.py b/Wrappers/Python/ccpi/plugins/processors/__init__.py new file mode 100644 index 0000000..e688671 --- /dev/null +++ b/Wrappers/Python/ccpi/plugins/processors/__init__.py @@ -0,0 +1,4 @@ +from .processors import setupCCPiGeometries +from .processors import CCPiForwardProjector +from .processors import CCPiBackwardProjector +from .processors import AcquisitionDataPadder diff --git a/Wrappers/Python/ccpi/plugins/processors/processors.py b/Wrappers/Python/ccpi/plugins/processors/processors.py new file mode 100755 index 0000000..84de79f --- /dev/null +++ b/Wrappers/Python/ccpi/plugins/processors/processors.py @@ -0,0 +1,309 @@ +# -*- 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 2018 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 + +from ccpi.framework import DataProcessor, AcquisitionData,\ + AcquisitionGeometry, ImageGeometry, ImageData +from ccpi.reconstruction.parallelbeam import alg as pbalg +import numpy + +def setupCCPiGeometries(ig, ag, counter): + Phantom_ccpi = ig.allocate(dimension_labels=[ImageGeometry.HORIZONTAL_X, + ImageGeometry.HORIZONTAL_Y, + ImageGeometry.VERTICAL]) + + voxel_per_pixel = 1 + angles = ag.angles + geoms = pbalg.pb_setup_geometry_from_image(Phantom_ccpi.as_array(), + angles, + voxel_per_pixel ) + + pg = AcquisitionGeometry('parallel', + '3D', + angles, + geoms['n_h'], 1.0, + geoms['n_v'], 1.0 #2D in 3D is a slice 1 pixel thick + ) + + center_of_rotation = Phantom_ccpi.get_dimension_size('horizontal_x') / 2 + #ad = AcquisitionData(geometry=pg,dimension_labels=['angle','vertical','horizontal']) + ad = pg.allocate(dimension_labels=[AcquisitionGeometry.ANGLE, + AcquisitionGeometry.VERTICAL, + AcquisitionGeometry.HORIZONTAL]) + geoms_i = pbalg.pb_setup_geometry_from_acquisition(ad.as_array(), + angles, + center_of_rotation, + voxel_per_pixel ) + + counter+=1 + + if counter < 4: + print (geoms, geoms_i) + if (not ( geoms_i == geoms )): + print ("not equal and {} {} {}".format(counter, geoms['output_volume_z'], geoms_i['output_volume_z'])) + X = max(geoms['output_volume_x'], geoms_i['output_volume_x']) + Y = max(geoms['output_volume_y'], geoms_i['output_volume_y']) + Z = max(geoms['output_volume_z'], geoms_i['output_volume_z']) + return setupCCPiGeometries(X,Y,Z,angles, counter) + else: + print ("happy now {} {} {}".format(counter, geoms['output_volume_z'], geoms_i['output_volume_z'])) + + return geoms + else: + return geoms_i + + +class CCPiForwardProjector(DataProcessor): + '''Normalization based on flat and dark + + This processor read in a AcquisitionData and normalises it based on + the instrument reading with and without incident photons or neutrons. + + Input: AcquisitionData + Parameter: 2D projection with flat field (or stack) + 2D projection with dark field (or stack) + Output: AcquisitionDataSetn + ''' + + def __init__(self, + image_geometry = None, + acquisition_geometry = None, + output_axes_order = None): + if output_axes_order is None: + # default ccpi projector image storing order + output_axes_order = ['angle','vertical','horizontal'] + + kwargs = { + 'image_geometry' : image_geometry, + 'acquisition_geometry' : acquisition_geometry, + 'output_axes_order' : output_axes_order, + 'default_image_axes_order' : ['horizontal_x','horizontal_y','vertical'], + 'default_acquisition_axes_order' : ['angle','vertical','horizontal'] + } + + super(CCPiForwardProjector, self).__init__(**kwargs) + + def check_input(self, dataset): + if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2: + # sort in the order that this projector needs it + return True + else: + raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ + .format(dataset.number_of_dimensions)) + + def process(self, out=None): + + volume = self.get_input() + volume_axes = volume.get_data_axes_order(new_order=self.default_image_axes_order) + if not volume_axes == [0,1,2]: + volume.array = numpy.transpose(volume.array, volume_axes) + pixel_per_voxel = 1 # should be estimated from image_geometry and + # acquisition_geometry + if self.acquisition_geometry.geom_type == 'parallel': + + pixels = pbalg.pb_forward_project(volume.as_array(), + self.acquisition_geometry.angles, + pixel_per_voxel) + + out = self.acquisition_geometry.allocate( + dimension_labels=self.output_axes_order) + out_axes = out.get_data_axes_order(new_order=self.output_axes_order) + if not out_axes == [0,1,2]: + print ("transpose") + pixels = numpy.transpose(pixels, out_axes) + out.fill(pixels) + + return out + else: + raise ValueError('Cannot process cone beam') + +class CCPiBackwardProjector(DataProcessor): + '''Backward projector + + This processor reads in a AcquisitionData and performs a backward projection, + i.e. project to reconstruction space. + Notice that it assumes that the center of rotation is in the middle + of the horizontal axis: in case when that's not the case it can be chained + with the AcquisitionDataPadder. + + Input: AcquisitionData + Parameter: 2D projection with flat field (or stack) + 2D projection with dark field (or stack) + Output: AcquisitionDataSetn + ''' + + def __init__(self, + image_geometry = None, + acquisition_geometry = None, + output_axes_order=None): + if output_axes_order is None: + # default ccpi projector image storing order + #output_axes_order = ['horizontal_x','horizontal_y','vertical'] + output_axes_order = ['vertical', 'horizontal_y','horizontal_x',] + kwargs = { + 'image_geometry' : image_geometry, + 'acquisition_geometry' : acquisition_geometry, + 'output_axes_order' : output_axes_order, + 'default_image_axes_order' : ['horizontal_x','horizontal_y','vertical'], + 'default_acquisition_axes_order' : ['angle','vertical','horizontal'] + } + + super(CCPiBackwardProjector, self).__init__(**kwargs) + + def check_input(self, dataset): + if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2: + + return True + else: + raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ + .format(dataset.number_of_dimensions)) + + def process(self, out=None): + projections = self.get_input() + projections_axes = projections.get_data_axes_order(new_order=self.default_acquisition_axes_order) + if not projections_axes == [0,1,2]: + projections.array = numpy.transpose(projections.array, projections_axes) + + pixel_per_voxel = 1 # should be estimated from image_geometry and acquisition_geometry + image_geometry = ImageGeometry(voxel_num_x = self.acquisition_geometry.pixel_num_h, + voxel_num_y = self.acquisition_geometry.pixel_num_h, + voxel_num_z = self.acquisition_geometry.pixel_num_v) + # input centered/padded acquisitiondata + center_of_rotation = projections.get_dimension_size('horizontal') / 2 + + if self.acquisition_geometry.geom_type == 'parallel': + back = pbalg.pb_backward_project( + projections.as_array(), + self.acquisition_geometry.angles, + center_of_rotation, + pixel_per_voxel + ) + out = self.image_geometry.allocate() + out_axes = out.get_data_axes_order(new_order=self.output_axes_order) + if not out_axes == [0,1,2]: + back = numpy.transpose(back, out_axes) + out.fill(back) + + return out + + else: + raise ValueError('Cannot process cone beam') + +class AcquisitionDataPadder(DataProcessor): + '''Normalization based on flat and dark + + This processor read in a AcquisitionData and normalises it based on + the instrument reading with and without incident photons or neutrons. + + Input: AcquisitionData + Parameter: 2D projection with flat field (or stack) + 2D projection with dark field (or stack) + Output: AcquisitionDataSetn + ''' + + def __init__(self, + center_of_rotation = None, + acquisition_geometry = None, + pad_value = 1e-5): + kwargs = { + 'acquisition_geometry' : acquisition_geometry, + 'center_of_rotation' : center_of_rotation, + 'pad_value' : pad_value + } + + super(AcquisitionDataPadder, self).__init__(**kwargs) + + def check_input(self, dataset): + if self.acquisition_geometry is None: + self.acquisition_geometry = dataset.geometry + if dataset.number_of_dimensions == 3: + return True + else: + raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ + .format(dataset.number_of_dimensions)) + + def process(self, out=None): + projections = self.get_input() + w = projections.get_dimension_size('horizontal') + delta = w - 2 * self.center_of_rotation + + padded_width = int ( + numpy.ceil(abs(delta)) + w + ) + delta_pix = padded_width - w + + voxel_per_pixel = 1 + geom = pbalg.pb_setup_geometry_from_acquisition(projections.as_array(), + self.acquisition_geometry.angles, + self.center_of_rotation, + voxel_per_pixel ) + + padded_geometry = self.acquisition_geometry.clone() + + padded_geometry.pixel_num_h = geom['n_h'] + padded_geometry.pixel_num_v = geom['n_v'] + + delta_pix_h = padded_geometry.pixel_num_h - self.acquisition_geometry.pixel_num_h + delta_pix_v = padded_geometry.pixel_num_v - self.acquisition_geometry.pixel_num_v + + if delta_pix_h == 0: + delta_pix_h = delta_pix + padded_geometry.pixel_num_h = padded_width + #initialize a new AcquisitionData with values close to 0 + out = AcquisitionData(geometry=padded_geometry) + out = out + self.pad_value + + + #pad in the horizontal-vertical plane -> slice on angles + if delta > 0: + #pad left of middle + command = "out.array[" + for i in range(out.number_of_dimensions): + if out.dimension_labels[i] == 'horizontal': + value = '{0}:{1}'.format(delta_pix_h, delta_pix_h+w) + command = command + str(value) + else: + if out.dimension_labels[i] == 'vertical' : + value = '{0}:'.format(delta_pix_v) + command = command + str(value) + else: + command = command + ":" + if i < out.number_of_dimensions -1: + command = command + ',' + command = command + '] = projections.array' + #print (command) + else: + #pad right of middle + command = "out.array[" + for i in range(out.number_of_dimensions): + if out.dimension_labels[i] == 'horizontal': + value = '{0}:{1}'.format(0, w) + command = command + str(value) + else: + if out.dimension_labels[i] == 'vertical' : + value = '{0}:'.format(delta_pix_v) + command = command + str(value) + else: + command = command + ":" + if i < out.number_of_dimensions -1: + command = command + ',' + command = command + '] = projections.array' + #print (command) + #cleaned = eval(command) + exec(command) + return out \ No newline at end of file diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py index f4c42e8..ad12d1c 100644 --- a/Wrappers/Python/setup.py +++ b/Wrappers/Python/setup.py @@ -32,7 +32,8 @@ if cil_version == '': setup( name="ccpi-plugins", version=cil_version, - packages=['ccpi' , 'ccpi.plugins'], + packages=['ccpi' , 'ccpi.plugins', 'ccpi.plugins.operators', + 'ccpi.plugins.processors'], # Project uses reStructuredText, so ensure that the docutils get # installed or upgraded on the target machine -- cgit v1.2.3