diff options
Diffstat (limited to 'Wrappers/Python/ccpi')
| -rw-r--r-- | Wrappers/Python/ccpi/astra/__init__.py | 0 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/astra/astra_ops.py | 132 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/astra/astra_processors.py | 205 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/astra/astra_utils.py | 61 |
4 files changed, 398 insertions, 0 deletions
diff --git a/Wrappers/Python/ccpi/astra/__init__.py b/Wrappers/Python/ccpi/astra/__init__.py new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/Wrappers/Python/ccpi/astra/__init__.py diff --git a/Wrappers/Python/ccpi/astra/astra_ops.py b/Wrappers/Python/ccpi/astra/astra_ops.py new file mode 100644 index 0000000..45b8238 --- /dev/null +++ b/Wrappers/Python/ccpi/astra/astra_ops.py @@ -0,0 +1,132 @@ +# -*- coding: utf-8 -*- +# This work is independent part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see <http://www.gnu.org/licenses/>. + +from ccpi.reconstruction.ops import Operator +import numpy +import astra +from ccpi.framework import AcquisitionData, ImageData, DataContainer +from ccpi.reconstruction.ops import PowerMethodNonsquare +from ccpi.astra.astra_processors import AstraForwardProjector, AstraBackProjector, \ + AstraForwardProjectorMC, AstraBackProjectorMC + +class AstraProjectorSimple(Operator): + """ASTRA projector modified to use DataSet and geometry.""" + def __init__(self, geomv, geomp, device): + super(AstraProjectorSimple, self).__init__() + + # Store volume and sinogram geometries. + self.sinogram_geometry = geomp + self.volume_geometry = geomv + + self.fp = AstraForwardProjector(volume_geometry=geomv, + sinogram_geometry=geomp, + proj_id=None, + device=device) + + self.bp = AstraBackProjector(volume_geometry=geomv, + sinogram_geometry=geomp, + proj_id=None, + device=device) + + # Initialise empty for singular value. + self.s1 = None + + def direct(self, IM): + self.fp.set_input(IM) + out = self.fp.get_output() + return out + + def adjoint(self, DATA): + self.bp.set_input(DATA) + out = self.bp.get_output() + return out + + #def delete(self): + # astra.data2d.delete(self.proj_id) + + def get_max_sing_val(self): + self.s1, sall, svec = PowerMethodNonsquare(self,10) + return self.s1 + + def size(self): + # Only implemented for 2D + return ( (self.sinogram_geometry.angles.size, \ + self.sinogram_geometry.pixel_num_h), \ + (self.volume_geometry.voxel_num_x, \ + self.volume_geometry.voxel_num_y) ) + + def create_image_data(self): + inputsize = self.size()[1] + return DataContainer(numpy.random.randn(inputsize[0], + inputsize[1])) + + +class AstraProjectorMC(Operator): + """ASTRA Multichannel projector""" + def __init__(self, geomv, geomp, device): + super(AstraProjectorMC, self).__init__() + + # Store volume and sinogram geometries. + self.sinogram_geometry = geomp + self.volume_geometry = geomv + + self.fp = AstraForwardProjectorMC(volume_geometry=geomv, + sinogram_geometry=geomp, + proj_id=None, + device=device) + + self.bp = AstraBackProjectorMC(volume_geometry=geomv, + sinogram_geometry=geomp, + proj_id=None, + device=device) + + # Initialise empty for singular value. + self.s1 = None + + def direct(self, IM): + self.fp.set_input(IM) + out = self.fp.get_output() + return out + + def adjoint(self, DATA): + self.bp.set_input(DATA) + out = self.bp.get_output() + return out + + #def delete(self): + # astra.data2d.delete(self.proj_id) + + def get_max_sing_val(self): + if self.s1 is None: + self.s1, sall, svec = PowerMethodNonsquare(self,10) + return self.s1 + else: + return self.s1 + + def size(self): + # Only implemented for 2D + return ( (self.sinogram_geometry.angles.size, \ + self.sinogram_geometry.pixel_num_h), \ + (self.volume_geometry.voxel_num_x, \ + self.volume_geometry.voxel_num_y) ) + + def create_image_data(self): + inputsize = self.size()[1] + return DataContainer(numpy.random.randn(self.volume_geometry.channels, + inputsize[0], + inputsize[1])) + diff --git a/Wrappers/Python/ccpi/astra/astra_processors.py b/Wrappers/Python/ccpi/astra/astra_processors.py new file mode 100644 index 0000000..3de43bf --- /dev/null +++ b/Wrappers/Python/ccpi/astra/astra_processors.py @@ -0,0 +1,205 @@ +from ccpi.framework import DataSetProcessor, ImageData, AcquisitionData +from ccpi.astra.astra_utils import convert_geometry_to_astra +import astra + + +class AstraForwardProjector(DataSetProcessor): + '''AstraForwardProjector + + Forward project ImageData to AcquisitionData using ASTRA proj_id. + + Input: ImageData + Parameter: proj_id + Output: AcquisitionData + ''' + + def __init__(self, + volume_geometry=None, + sinogram_geometry=None, + proj_id=None, + device='cpu'): + kwargs = { + 'volume_geometry' : volume_geometry, + 'sinogram_geometry' : sinogram_geometry, + 'proj_id' : proj_id, + 'device' : device + } + + #DataSetProcessor.__init__(self, **kwargs) + super(AstraForwardProjector, self).__init__(**kwargs) + + self.set_ImageGeometry(volume_geometry) + self.set_AcquisitionGeometry(sinogram_geometry) + + # Set up ASTRA Volume and projection geometry, not to be stored in self + vol_geom, proj_geom = convert_geometry_to_astra(self.volume_geometry, + self.sinogram_geometry) + + # ASTRA projector, to be stored + if device == 'cpu': + # Note that 'line' only one option + if self.sinogram_geometry.geom_type == 'parallel': + self.set_projector(astra.create_projector('line', proj_geom, vol_geom) ) + elif self.sinogram_geometry.geom_type == 'cone': + self.set_projector(astra.create_projector('line_fanflat', proj_geom, vol_geom) ) + else: + NotImplemented + elif device == 'gpu': + self.set_projector(astra.create_projector('cuda', proj_geom, vol_geom) ) + else: + NotImplemented + + 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 set_projector(self, proj_id): + self.proj_id = proj_id + + def set_ImageGeometry(self, volume_geometry): + self.volume_geometry = volume_geometry + + def set_AcquisitionGeometry(self, sinogram_geometry): + self.sinogram_geometry = sinogram_geometry + + def process(self): + IM = self.get_input() + DATA = AcquisitionData(geometry=self.sinogram_geometry) + #sinogram_id, DATA = astra.create_sino( IM.as_array(), + # self.proj_id) + sinogram_id, DATA.array = astra.create_sino(IM.as_array(), + self.proj_id) + astra.data2d.delete(sinogram_id) + #return AcquisitionData(array=DATA, geometry=self.sinogram_geometry) + return DATA + +class AstraBackProjector(DataSetProcessor): + '''AstraBackProjector + + Back project AcquisitionData to ImageData using ASTRA proj_id. + + Input: AcquisitionData + Parameter: proj_id + Output: ImageData + ''' + + def __init__(self, + volume_geometry=None, + sinogram_geometry=None, + proj_id=None, + device='cpu'): + kwargs = { + 'volume_geometry' : volume_geometry, + 'sinogram_geometry' : sinogram_geometry, + 'proj_id' : proj_id, + 'device' : device + } + + #DataSetProcessor.__init__(self, **kwargs) + super(AstraBackProjector, self).__init__(**kwargs) + + self.set_ImageGeometry(volume_geometry) + self.set_AcquisitionGeometry(sinogram_geometry) + + # Set up ASTRA Volume and projection geometry, not to be stored in self + vol_geom, proj_geom = convert_geometry_to_astra(self.volume_geometry, + self.sinogram_geometry) + + # ASTRA projector, to be stored + if device == 'cpu': + # Note that 'line' only one option + if self.sinogram_geometry.geom_type == 'parallel': + self.set_projector(astra.create_projector('line', proj_geom, vol_geom) ) + elif self.sinogram_geometry.geom_type == 'cone': + self.set_projector(astra.create_projector('line_fanflat', proj_geom, vol_geom) ) + else: + NotImplemented + elif device == 'gpu': + self.set_projector(astra.create_projector('cuda', proj_geom, vol_geom) ) + else: + NotImplemented + + 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 set_projector(self, proj_id): + self.proj_id = proj_id + + def set_ImageGeometry(self, volume_geometry): + self.volume_geometry = volume_geometry + + def set_AcquisitionGeometry(self, sinogram_geometry): + self.sinogram_geometry = sinogram_geometry + + def process(self): + DATA = self.get_input() + IM = ImageData(geometry=self.volume_geometry) + rec_id, IM.array = astra.create_backprojection(DATA.as_array(), + self.proj_id) + astra.data2d.delete(rec_id) + return IM + +class AstraForwardProjectorMC(AstraForwardProjector): + '''AstraForwardProjector Multi channel + + Forward project ImageData to AcquisitionDataSet using ASTRA proj_id. + + Input: ImageDataSet + Parameter: proj_id + Output: AcquisitionData + ''' + def check_input(self, dataset): + if dataset.number_of_dimensions == 2 or \ + dataset.number_of_dimensions == 3 or \ + dataset.number_of_dimensions == 4: + return True + else: + raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ + .format(dataset.number_of_dimensions)) + def process(self): + IM = self.get_input() + #create the output AcquisitionData + DATA = AcquisitionData(geometry=self.sinogram_geometry) + + for k in range(DATA.geometry.channels): + sinogram_id, DATA.as_array()[k] = astra.create_sino(IM.as_array()[k], + self.proj_id) + astra.data2d.delete(sinogram_id) + return DATA + +class AstraBackProjectorMC(AstraBackProjector): + '''AstraBackProjector Multi channel + + Back project AcquisitionData to ImageData using ASTRA proj_id. + + Input: AcquisitionData + Parameter: proj_id + Output: ImageData + ''' + def check_input(self, dataset): + if dataset.number_of_dimensions == 2 or \ + dataset.number_of_dimensions == 3 or \ + dataset.number_of_dimensions == 4: + return True + else: + raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ + .format(dataset.number_of_dimensions)) + def process(self): + DATA = self.get_input() + + IM = ImageData(geometry=self.volume_geometry) + + for k in range(IM.geometry.channels): + rec_id, IM.as_array()[k] = astra.create_backprojection( + DATA.as_array()[k], + self.proj_id) + astra.data2d.delete(rec_id) + return IM
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/astra/astra_utils.py b/Wrappers/Python/ccpi/astra/astra_utils.py new file mode 100644 index 0000000..9f8fe46 --- /dev/null +++ b/Wrappers/Python/ccpi/astra/astra_utils.py @@ -0,0 +1,61 @@ +import astra + +def convert_geometry_to_astra(volume_geometry, sinogram_geometry): + # Set up ASTRA Volume and projection geometry, not stored + if sinogram_geometry.dimension == '2D': + vol_geom = astra.create_vol_geom(volume_geometry.voxel_num_x, + volume_geometry.voxel_num_y, + volume_geometry.get_min_x(), + volume_geometry.get_max_x(), + volume_geometry.get_min_y(), + volume_geometry.get_max_y()) + + if sinogram_geometry.geom_type == 'parallel': + proj_geom = astra.create_proj_geom('parallel', + sinogram_geometry.pixel_size_h, + sinogram_geometry.pixel_num_h, + sinogram_geometry.angles) + elif sinogram_geometry.geom_type == 'cone': + proj_geom = astra.create_proj_geom('fanflat', + sinogram_geometry.pixel_size_h, + sinogram_geometry.pixel_num_h, + sinogram_geometry.angles, + sinogram_geometry.dist_source_center, + sinogram_geometry.dist_center_detector) + else: + NotImplemented + + elif sinogram_geometry.dimension == '3D': + vol_geom = astra.create_vol_geom(volume_geometry.voxel_num_x, + volume_geometry.voxel_num_y, + volume_geometry.voxel_num_z, + volume_geometry.get_min_x(), + volume_geometry.get_max_x(), + volume_geometry.get_min_y(), + volume_geometry.get_max_y(), + volume_geometry.get_min_z(), + volume_geometry.get_max_z()) + + if sinogram_geometry.geom_type == 'parallel': + proj_geom = astra.create_proj_geom('parallel3d', + sinogram_geometry.pixel_size_h, + sinogram_geometry.pixel_size_v, + sinogram_geometry.pixel_num_v, + sinogram_geometry.pixel_num_h, + sinogram_geometry.angles) + elif sinogram_geometry.geom_type == 'cone': + proj_geom = astra.create_proj_geom('cone', + sinogram_geometry.pixel_size_h, + sinogram_geometry.pixel_size_v, + sinogram_geometry.pixel_num_v, + sinogram_geometry.pixel_num_h, + sinogram_geometry.angles, + sinogram_geometry.dist_source_center, + sinogram_geometry.dist_center_detector) + else: + NotImplemented + + else: + NotImplemented + + return vol_geom, proj_geom
\ No newline at end of file |