diff options
Diffstat (limited to 'Wrappers')
| -rw-r--r-- | Wrappers/Python/ccpi/astra/astra_processors.py | 273 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/astra/astra_utils.py | 61 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/framework.py | 30 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/reconstruction/astra_ops.py | 161 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/reconstruction/geoms.py | 45 | ||||
| -rw-r--r-- | Wrappers/Python/test/simple_mc_demo.py | 134 |
6 files changed, 582 insertions, 122 deletions
diff --git a/Wrappers/Python/ccpi/astra/astra_processors.py b/Wrappers/Python/ccpi/astra/astra_processors.py new file mode 100644 index 0000000..650e11b --- /dev/null +++ b/Wrappers/Python/ccpi/astra/astra_processors.py @@ -0,0 +1,273 @@ +from ccpi.framework import DataSetProcessor, DataSet, VolumeData, SinogramData +from ccpi.astra.astra_utils import convert_geometry_to_astra +import astra +import numpy + + +class AstraForwardProjector(DataSetProcessor): + '''AstraForwardProjector + + Forward project VolumeDataSet to SinogramDataSet using ASTRA proj_id. + + Input: VolumeDataSet + Parameter: proj_id + Output: SinogramDataSet + ''' + + 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.setVolumeGeometry(volume_geometry) + self.setSinogramGeometry(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' is only for parallel (2D) and only one option + self.setProjector(astra.create_projector('line', proj_geom, vol_geom) ) + elif device == 'gpu': + self.setProjector(astra.create_projector('cuda', proj_geom, vol_geom) ) + else: + NotImplemented + + def checkInput(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 setProjector(self, proj_id): + self.proj_id = proj_id + + def setVolumeGeometry(self, volume_geometry): + self.volume_geometry = volume_geometry + + def setSinogramGeometry(self, sinogram_geometry): + self.sinogram_geometry = sinogram_geometry + + def process(self): + IM = self.getInput() + sinogram_id, DATA = astra.create_sino(IM.as_array(), self.proj_id) + astra.data2d.delete(sinogram_id) + return SinogramData(DATA,geometry=self.sinogram_geometry) + +class AstraBackProjector(DataSetProcessor): + '''AstraBackProjector + + Back project SinogramDataSet to VolumeDataSet using ASTRA proj_id. + + Input: SinogramDataSet + Parameter: proj_id + Output: VolumeDataSet + ''' + + 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.setVolumeGeometry(volume_geometry) + self.setSinogramGeometry(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' is only for parallel (2D) and only one option + self.setProjector(astra.create_projector('line', proj_geom, vol_geom) ) + elif device == 'gpu': + self.setProjector(astra.create_projector('cuda', proj_geom, vol_geom) ) + else: + NotImplemented + + def checkInput(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 setProjector(self, proj_id): + self.proj_id = proj_id + + def setVolumeGeometry(self, volume_geometry): + self.volume_geometry = volume_geometry + + def setSinogramGeometry(self, sinogram_geometry): + self.sinogram_geometry = sinogram_geometry + + def process(self): + DATA = self.getInput() + rec_id, IM = astra.create_backprojection(DATA.as_array(), self.proj_id) + astra.data2d.delete(rec_id) + return VolumeData(IM,geometry=self.volume_geometry) + +class AstraForwardProjectorMC(DataSetProcessor): + '''AstraForwardProjector + + Forward project VolumeDataSet to SinogramDataSet using ASTRA proj_id. + + Input: VolumeDataSet + Parameter: proj_id + Output: SinogramDataSet + ''' + + 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(AstraForwardProjectorMC, self).__init__(**kwargs) + + self.setVolumeGeometry(volume_geometry) + self.setSinogramGeometry(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' is only for parallel (2D) and only one option + self.setProjector(astra.create_projector('line', proj_geom, vol_geom) ) + elif device == 'gpu': + self.setProjector(astra.create_projector('cuda', proj_geom, vol_geom) ) + else: + NotImplemented + + def checkInput(self, dataset): + if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2: + return True + else: + return True + #raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ + # .format(dataset.number_of_dimensions)) + + def setProjector(self, proj_id): + self.proj_id = proj_id + + def setVolumeGeometry(self, volume_geometry): + self.volume_geometry = volume_geometry + + def setSinogramGeometry(self, sinogram_geometry): + self.sinogram_geometry = sinogram_geometry + + def process(self): + IM = self.getInput() + DATA = numpy.zeros((self.sinogram_geometry.angles.size, + self.sinogram_geometry.pixel_num_h, + 1, + self.sinogram_geometry.channels), + 'float32') + for k in range(self.sinogram_geometry.channels): + sinogram_id, DATA[:,:,0,k] = astra.create_sino(IM.as_array()[:,:,0,k], + self.proj_id) + astra.data2d.delete(sinogram_id) + return SinogramData(DATA,geometry=self.sinogram_geometry) + +class AstraBackProjectorMC(DataSetProcessor): + '''AstraBackProjector + + Back project SinogramDataSet to VolumeDataSet using ASTRA proj_id. + + Input: SinogramDataSet + Parameter: proj_id + Output: VolumeDataSet + ''' + + 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(AstraBackProjectorMC, self).__init__(**kwargs) + + self.setVolumeGeometry(volume_geometry) + self.setSinogramGeometry(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' is only for parallel (2D) and only one option + self.setProjector(astra.create_projector('line', proj_geom, vol_geom) ) + elif device == 'gpu': + self.setProjector(astra.create_projector('cuda', proj_geom, vol_geom) ) + else: + NotImplemented + + def checkInput(self, dataset): + if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2: + return True + else: + return True + #raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ + # .format(dataset.number_of_dimensions)) + + def setProjector(self, proj_id): + self.proj_id = proj_id + + def setVolumeGeometry(self, volume_geometry): + self.volume_geometry = volume_geometry + + def setSinogramGeometry(self, sinogram_geometry): + self.sinogram_geometry = sinogram_geometry + + def process(self): + DATA = self.getInput() + IM = numpy.zeros((self.volume_geometry.voxel_num_x, + self.volume_geometry.voxel_num_y, + 1, + self.volume_geometry.channels), + 'float32') + for k in range(self.volume_geometry.channels): + rec_id, IM[:,:,0,k] = astra.create_backprojection(DATA.as_array()[:,:,0,k], + self.proj_id) + astra.data2d.delete(rec_id) + return VolumeData(IM,geometry=self.volume_geometry)
\ 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..8b5043a --- /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.getMinX(), + volume_geometry.getMaxX(), + volume_geometry.getMinY(), + volume_geometry.getMaxY()) + + 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.getMinX(), + volume_geometry.getMaxX(), + volume_geometry.getMinY(), + volume_geometry.getMaxY(), + volume_geometry.getMinZ(), + volume_geometry.getMaxZ()) + + if sinogram_geometry.proj_geom == '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 diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py index 754e0ee..82ad65c 100644 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework.py @@ -437,20 +437,20 @@ class VolumeData(DataSet): if type(array) == DataSet: # if the array is a DataSet get the info from there - if not ( array.number_of_dimensions == 2 or \ - array.number_of_dimensions == 3 ): - raise ValueError('Number of dimensions are not 2 or 3: {0}'\ - .format(array.number_of_dimensions)) + #if not ( array.number_of_dimensions == 2 or \ + # array.number_of_dimensions == 3 ): + # raise ValueError('Number of dimensions are not 2 or 3: {0}'\ + # .format(array.number_of_dimensions)) #DataSet.__init__(self, array.as_array(), deep_copy, # array.dimension_labels, **kwargs) super(VolumeData, self).__init__(array.as_array(), deep_copy, array.dimension_labels, **kwargs) elif type(array) == numpy.ndarray: - if not ( array.ndim == 3 or array.ndim == 2 ): - raise ValueError( - 'Number of dimensions are not 3 or 2 : {0}'\ - .format(array.ndim)) + #if not ( array.ndim == 3 or array.ndim == 2 ): + # raise ValueError( + # 'Number of dimensions are not 3 or 2 : {0}'\ + # .format(array.ndim)) if dimension_labels is None: if array.ndim == 3: @@ -485,17 +485,17 @@ class SinogramData(DataSet): if type(array) == DataSet: # if the array is a DataSet get the info from there - if not ( array.number_of_dimensions == 2 or \ - array.number_of_dimensions == 3 ): - raise ValueError('Number of dimensions are not 2 or 3: {0}'\ - .format(array.number_of_dimensions)) + #if not ( array.number_of_dimensions == 2 or \ + # array.number_of_dimensions == 3 ): + # raise ValueError('Number of dimensions are not 2 or 3: {0}'\ + # .format(array.number_of_dimensions)) DataSet.__init__(self, array.as_array(), deep_copy, array.dimension_labels, **kwargs) elif type(array) == numpy.ndarray: - if not ( array.ndim == 3 or array.ndim == 2 ): - raise ValueError('Number of dimensions are != 3: {0}'\ - .format(array.ndim)) + #if not ( array.ndim == 3 or array.ndim == 2 ): + # raise ValueError('Number of dimensions are != 3: {0}'\ + # .format(array.ndim)) if dimension_labels is None: if array.ndim == 3: dimension_labels = ['angle' , diff --git a/Wrappers/Python/ccpi/reconstruction/astra_ops.py b/Wrappers/Python/ccpi/reconstruction/astra_ops.py index 452c86a..1346f22 100644 --- a/Wrappers/Python/ccpi/reconstruction/astra_ops.py +++ b/Wrappers/Python/ccpi/reconstruction/astra_ops.py @@ -16,118 +16,105 @@ # along with this program. If not, see <http://www.gnu.org/licenses/>. from ccpi.reconstruction.ops import Operator -import astra import numpy from ccpi.framework import SinogramData, VolumeData from ccpi.reconstruction.ops import PowerMethodNonsquare +from ccpi.astra.astra_processors import * class AstraProjectorSimple(Operator): """ASTRA projector modified to use DataSet and geometry.""" def __init__(self, geomv, geomp, device): super(AstraProjectorSimple, self).__init__() - # Store our volume and sinogram geometries. Redundant with also - # storing in ASTRA format below but needed to assign to - # SinogramData in "direct" method and VolumeData in "adjoint" method + # Store volume and sinogram geometries. self.sinogram_geometry = geomp self.volume_geometry = geomv - # ASTRA Volume geometry - if geomp.dimension == '2D': - self.vol_geom = astra.create_vol_geom(geomv.voxel_num_x, - geomv.voxel_num_y, - geomv.getMinX(), - geomv.getMaxX(), - geomv.getMinY(), - geomv.getMaxY()) - elif geomp.dimension == '3D': - self.vol_geom = astra.create_vol_geom(geomv.voxel_num_x, - geomv.voxel_num_y, - geomv.voxel_num_z, - geomv.getMinX(), - geomv.getMaxX(), - geomv.getMinY(), - geomv.getMaxY(), - geomv.getMinZ(), - geomv.getMaxZ()) - else: - NotImplemented - - - # ASTRA Projections geometry - if geomp.dimension == '2D': - if geomp.geom_type == 'parallel': - self.proj_geom = astra.create_proj_geom('parallel', - geomp.pixel_size_h, - geomp.pixel_num_h, - geomp.angles) - elif geomp.geom_type == 'cone': - self.proj_geom = astra.create_proj_geom('fanflat', - geomp.pixel_size_h, - geomp.pixel_num_h, - geomp.angles, - geomp.dist_source_center, - geomp.dist_center_detector) - else: - NotImplemented - elif geomp.dimension == '3D': - if geomp.proj_geom == 'parallel': - self.proj_geom = astra.create_proj_geom('parallel3d', - geomp.pixel_size_h, - geomp.pixel_size_v, - geomp.pixel_num_v, - geomp.pixel_num_h, - geomp.angles) - elif geomp.geom_type == 'cone': - self.proj_geom = astra.create_proj_geom('cone', - geomp.pixel_size_h, - geomp.pixel_size_v, - geomp.pixel_num_v, - geomp.pixel_num_h, - geomp.angles, - geomp.dist_source_center, - geomp.dist_center_detector) - else: - NotImplemented - else: - NotImplemented - - # ASTRA projector - if device == 'cpu': - # Note that 'line' is only for parallel (2D) and only one option - self.proj_id = astra.create_projector('line', self.proj_geom, - self.vol_geom) # for CPU - elif device == 'gpu': - self.proj_id = astra.create_projector('cuda', self.proj_geom, - self.vol_geom) # for GPU - else: - NotImplemented + 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): - - sinogram_id, DATA = astra.create_sino(IM.as_array(), self.proj_id) - astra.data2d.delete(sinogram_id) - return SinogramData(DATA,geometry=self.sinogram_geometry) + self.fp.setInput(IM) + out = self.fp.getOutput() + return out def adjoint(self, DATA): - rec_id, IM = astra.create_backprojection(DATA.as_array(), self.proj_id) - astra.data2d.delete(rec_id) - return VolumeData(IM,geometry=self.volume_geometry) + self.bp.setInput(DATA) + out = self.bp.getOutput() + return out - def delete(self): - astra.data2d.delete(self.proj_id) + #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): - return ( (self.proj_geom['ProjectionAngles'].size, \ - self.proj_geom['DetectorCount']), \ - (self.vol_geom['GridColCount'], \ - self.vol_geom['GridRowCount']) ) + # 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) ) + +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 = 50 + + def direct(self, IM): + self.fp.setInput(IM) + out = self.fp.getOutput() + return out + + def adjoint(self, DATA): + self.bp.setInput(DATA) + out = self.bp.getOutput() + 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) ) class AstraProjector(Operator): diff --git a/Wrappers/Python/ccpi/reconstruction/geoms.py b/Wrappers/Python/ccpi/reconstruction/geoms.py index edce3b3..9f81fa0 100644 --- a/Wrappers/Python/ccpi/reconstruction/geoms.py +++ b/Wrappers/Python/ccpi/reconstruction/geoms.py @@ -1,16 +1,17 @@ class VolumeGeometry: - def __init__(self, \ - voxel_num_x=None, \ - voxel_num_y=None, \ - voxel_num_z=None, \ - voxel_size_x=None, \ - voxel_size_y=None, \ - voxel_size_z=None, \ - center_x=0, \ - center_y=0, \ - center_z=0): + def __init__(self, + voxel_num_x=None, + voxel_num_y=None, + voxel_num_z=None, + voxel_size_x=None, + voxel_size_y=None, + voxel_size_z=None, + center_x=0, + center_y=0, + center_z=0, + channels=1): self.voxel_num_x = voxel_num_x self.voxel_num_y = voxel_num_y @@ -21,6 +22,7 @@ class VolumeGeometry: self.center_x = center_x self.center_y = center_y self.center_z = center_z + self.channels = channels def getMinX(self): return self.center_x - 0.5*self.voxel_num_x*self.voxel_size_x @@ -43,16 +45,17 @@ class VolumeGeometry: class SinogramGeometry: - def __init__(self, \ - geom_type, \ - dimension, \ - angles, \ - pixel_num_h=None, \ - pixel_size_h=1, \ - pixel_num_v=None, \ - pixel_size_v=1, \ - dist_source_center=None, \ - dist_center_detector=None, \ + def __init__(self, + geom_type, + dimension, + angles, + pixel_num_h=None, + pixel_size_h=1, + pixel_num_v=None, + pixel_size_v=1, + dist_source_center=None, + dist_center_detector=None, + channels=1 ): """ General inputs for standard type projection geometries @@ -91,6 +94,8 @@ class SinogramGeometry: self.pixel_size_h = pixel_size_h self.pixel_num_v = pixel_num_v self.pixel_size_v = pixel_size_v + + self.channels = channels diff --git a/Wrappers/Python/test/simple_mc_demo.py b/Wrappers/Python/test/simple_mc_demo.py new file mode 100644 index 0000000..a6959f9 --- /dev/null +++ b/Wrappers/Python/test/simple_mc_demo.py @@ -0,0 +1,134 @@ + + +import sys + +sys.path.append("..") + +from ccpi.framework import * +from ccpi.reconstruction.algs import * +from ccpi.reconstruction.funcs import * +from ccpi.reconstruction.ops import * +from ccpi.reconstruction.astra_ops import * +from ccpi.reconstruction.geoms import * + +import numpy as np +import matplotlib.pyplot as plt + +test_case = 2 # 1=parallel2D, 2=cone2D + +# Set up phantom +N = 128 + +numchannels = 3 + +x = np.zeros((N,N,1,numchannels)) + +x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,0] = 1.0 +x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,0] = 2.0 + +x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,1] = 0.7 +x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,1] = 1.2 + +x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,2] = 1.5 +x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,2] = 2.2 + +f, axarr = plt.subplots(1,numchannels) +for k in numpy.arange(3): + axarr[k].imshow(x[:,:,0,k],vmin=0,vmax=2.5) +plt.show() + +vg = VolumeGeometry(N,N,None, 1,1,None,channels=numchannels) + + +Phantom = VolumeData(x,geometry=vg) + +# Set up measurement geometry +angles_num = 20; # angles number + +if test_case==1: + angles = np.linspace(0,np.pi,angles_num,endpoint=False) +elif test_case==2: + angles = np.linspace(0,2*np.pi,angles_num,endpoint=False) +else: + NotImplemented + +det_w = 1.0 +det_num = N +SourceOrig = 200 +OrigDetec = 0 + +# Parallelbeam geometry test +if test_case==1: + pg = SinogramGeometry('parallel', + '2D', + angles, + det_num, + det_w, + channels=numchannels) +elif test_case==2: + pg = SinogramGeometry('cone', + '2D', + angles, + det_num, + det_w, + dist_source_center=SourceOrig, + dist_center_detector=OrigDetec, + channels=numchannels) + +# ASTRA operator using volume and sinogram geometries +Aop = AstraProjectorMC(vg, pg, 'gpu') + + +# Try forward and backprojection +b = Aop.direct(Phantom) + +fb, axarrb = plt.subplots(1,numchannels) +for k in numpy.arange(3): + axarrb[k].imshow(b.as_array()[:,:,0,k],vmin=0,vmax=250) +plt.show() + +out2 = Aop.adjoint(b) + +fo, axarro = plt.subplots(1,numchannels) +for k in range(3): + axarro[k].imshow(out2.as_array()[:,:,0,k],vmin=0,vmax=3500) +plt.show() + +# Create least squares object instance with projector and data. +f = Norm2sq(Aop,b,c=0.5) + +# Initial guess +x_init = VolumeData(np.zeros(x.shape),geometry=vg) + +# FISTA options +opt = {'tol': 1e-4, 'iter': 200} + +# Run FISTA for least squares without regularization +x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt) + + +ff0, axarrf0 = plt.subplots(1,numchannels) +for k in numpy.arange(3): + axarrf0[k].imshow(x_fista0.as_array()[:,:,0,k],vmin=0,vmax=2.5) +plt.show() + +plt.semilogy(criter0) +plt.title('Criterion vs iterations, least squares') +plt.show() + +# Now least squares plus 1-norm regularization +lam = 0.1 +g0 = Norm1(lam) + + +# Run FISTA for least squares plus 1-norm function. +x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt) + +ff1, axarrf1 = plt.subplots(1,numchannels) +for k in numpy.arange(3): + axarrf1[k].imshow(x_fista1.as_array()[:,:,0,k],vmin=0,vmax=2.5) +plt.show() + +plt.semilogy(criter1) +plt.title('Criterion vs iterations, least squares plus 1-norm regu') +plt.show()
\ No newline at end of file |