diff options
Diffstat (limited to 'Wrappers')
| -rw-r--r-- | Wrappers/Python/ccpi/astra/astra_processors.py | 75 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/framework.py | 244 | ||||
| -rwxr-xr-x | Wrappers/Python/ccpi/processors.py | 14 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/reconstruction/astra_ops.py | 51 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/reconstruction/geoms.py | 55 | ||||
| -rw-r--r-- | Wrappers/Python/setup.py | 2 | ||||
| -rw-r--r-- | Wrappers/Python/wip/DemoRecIP.py (renamed from Wrappers/Python/test/DemoRecIP.py) | 0 | ||||
| -rw-r--r-- | Wrappers/Python/wip/regularizers.py (renamed from Wrappers/Python/test/regularizers.py) | 0 | ||||
| -rw-r--r-- | Wrappers/Python/wip/simple_demo.py (renamed from Wrappers/Python/test/simple_demo.py) | 28 | ||||
| -rwxr-xr-x | Wrappers/Python/wip/simple_mc_demo.py | 143 |
10 files changed, 501 insertions, 111 deletions
diff --git a/Wrappers/Python/ccpi/astra/astra_processors.py b/Wrappers/Python/ccpi/astra/astra_processors.py index c3127eb..f51aec9 100644 --- a/Wrappers/Python/ccpi/astra/astra_processors.py +++ b/Wrappers/Python/ccpi/astra/astra_processors.py @@ -1,6 +1,7 @@ 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): @@ -50,7 +51,8 @@ class AstraForwardProjector(DataSetProcessor): NotImplemented def checkInput(self, dataset): - if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2: + 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}"\ @@ -67,9 +69,14 @@ class AstraForwardProjector(DataSetProcessor): def process(self): IM = self.getInput() - sinogram_id, DATA = astra.create_sino(IM.as_array(), self.proj_id) + DATA = SinogramData(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 SinogramData(DATA,geometry=self.sinogram_geometry) + #return SinogramData(array=DATA, geometry=self.sinogram_geometry) + return DATA class AstraBackProjector(DataSetProcessor): '''AstraBackProjector @@ -135,6 +142,64 @@ class AstraBackProjector(DataSetProcessor): def process(self): DATA = self.getInput() - rec_id, IM = astra.create_backprojection(DATA.as_array(), self.proj_id) + IM = VolumeData(geometry=self.volume_geometry) + rec_id, IM.array = astra.create_backprojection(DATA.as_array(), + self.proj_id) astra.data2d.delete(rec_id) - return VolumeData(IM,geometry=self.volume_geometry)
\ No newline at end of file + return IM + +class AstraForwardProjectorMC(AstraForwardProjector): + '''AstraForwardProjector Multi channel + + Forward project VolumeDataSet to SinogramDataSet using ASTRA proj_id. + + Input: VolumeDataSet + Parameter: proj_id + Output: SinogramDataSet + ''' + def checkInput(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.getInput() + #create the output Sinogram + DATA = SinogramData(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 SinogramDataSet to VolumeDataSet using ASTRA proj_id. + + Input: SinogramDataSet + Parameter: proj_id + Output: VolumeDataSet + ''' + def checkInput(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.getInput() + + IM = VolumeData(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/framework.py b/Wrappers/Python/ccpi/framework.py index 5ac17ee..5a507d9 100644 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework.py @@ -23,6 +23,7 @@ import numpy import sys from datetime import timedelta, datetime import warnings +from ccpi.reconstruction import geoms if sys.version_info[0] >= 3 and sys.version_info[1] >= 4: ABC = abc.ABC @@ -119,7 +120,7 @@ class DataSet(object): if deep_copy: self.array = array[:] else: - self.array = array + self.array = array else: raise TypeError('Array must be NumpyArray, passed {0}'\ .format(type(array))) @@ -141,8 +142,9 @@ class DataSet(object): if dimensions is not None: return self.subset(dimensions).as_array() return self.array - - def subset(self, dimensions=None): + + + def subset(self, dimensions=None, **kw): '''Creates a DataSet containing a subset of self according to the labels in dimensions''' if dimensions is None: @@ -164,7 +166,7 @@ class DataSet(object): else: axis_order.append(find_key(self.dimension_labels, dl)) if not proceed: - raise KeyError('Unknown key specified {0}'.format(dl)) + raise KeyError('Subset error: Unknown key specified {0}'.format(dl)) # slice away the unwanted data from the array unwanted_dimensions = self.dimension_labels.copy() @@ -176,13 +178,21 @@ class DataSet(object): #print ("left_dimensions {0}".format(left_dimensions)) #new_shape = [self.shape[ax] for ax in axis_order] #print ("new_shape {0}".format(new_shape)) - command = "self.array" + command = "self.array[" for i in range(self.number_of_dimensions): if self.dimension_labels[i] in unwanted_dimensions.values(): - command = command + "[0]" + value = 0 + for k,v in kw.items(): + if k == self.dimension_labels[i]: + value = v + + command = command + str(value) else: - command = command + "[:]" - #print ("command {0}".format(command)) + command = command + ":" + if i < self.number_of_dimensions -1: + command = command + ',' + command = command + ']' + cleaned = eval(command) # cleaned has collapsed dimensions in the same order of # self.array, but we want it in the order stated in the @@ -413,7 +423,7 @@ class DataSet(object): repres += "Number of dimensions: {0}\n".format(self.number_of_dimensions) repres += "Shape: {0}\n".format(self.shape) repres += "Axis labels: {0}\n".format(self.dimension_labels) - repres += "Representation: {0}\n".format(self.array) + repres += "Representation: \n{0}\n".format(self.array) return repres @@ -421,40 +431,81 @@ class DataSet(object): class VolumeData(DataSet): '''DataSet for holding 2D or 3D dataset''' def __init__(self, - array, + array = None, deep_copy=True, dimension_labels=None, **kwargs): - 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)) - - #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)) + self.geometry = None + if array is None: + if 'geometry' in kwargs.keys(): + geometry = kwargs['geometry'] + self.geometry = geometry + channels = geometry.channels + horiz_x = geometry.voxel_num_x + horiz_y = geometry.voxel_num_y + vert = 1 if geometry.voxel_num_z is None\ + else geometry.voxel_num_z # this should be 1 for 2D - if dimension_labels is None: - if array.ndim == 3: - dimension_labels = ['horizontal_x' , - 'horizontal_y' , - 'vertical'] + if channels > 1: + if vert > 1: + shape = (channels, vert, horiz_y, horiz_x) + dim_labels = ['channel' ,'vertical' , 'horizontal_y' , + 'horizontal_x'] + else: + shape = (channels , horiz_y, horiz_x) + dim_labels = ['channel' , 'horizontal_y' , + 'horizontal_x'] else: - dimension_labels = ['horizontal' , - 'vertical'] - - #DataSet.__init__(self, array, deep_copy, dimension_labels, **kwargs) - super(VolumeData, self).__init__(array, deep_copy, - dimension_labels, **kwargs) + if vert > 1: + shape = (vert, horiz_y, horiz_x) + dim_labels = ['vertical' , 'horizontal_y' , + 'horizontal_x'] + else: + shape = (horiz_y, horiz_x) + dim_labels = ['horizontal_y' , + 'horizontal_x'] + + array = numpy.zeros( shape , dtype=numpy.float32) + super(VolumeData, self).__init__(array, deep_copy, + dim_labels, **kwargs) + + else: + raise ValueError('Please pass either a DataSet, ' +\ + 'a numpy array or a geometry') + else: + 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 or \ + array.number_of_dimensions == 4): + raise ValueError('Number of dimensions are not 2 or 3 or 4: {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 == 2 or array.ndim == 3 or array.ndim == 4 ): + raise ValueError( + 'Number of dimensions are not 2 or 3 or 4 : {0}'\ + .format(array.ndim)) + + if dimension_labels is None: + if array.ndim == 4: + dimension_labels = ['channel' ,'vertical' , 'horizontal_y' , + 'horizontal_x'] + elif array.ndim == 3: + dimension_labels = ['vertical' , 'horizontal_y' , + 'horizontal_x'] + else: + dimension_labels = ['horizontal_y' , + 'horizontal_x'] + + #DataSet.__init__(self, array, deep_copy, dimension_labels, **kwargs) + super(VolumeData, self).__init__(array, deep_copy, + dimension_labels, **kwargs) # load metadata from kwargs if present for key, value in kwargs.items(): @@ -469,34 +520,79 @@ class VolumeData(DataSet): class SinogramData(DataSet): '''DataSet for holding 2D or 3D sinogram''' def __init__(self, - array, + array = None, deep_copy=True, dimension_labels=None, **kwargs): - - 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)) - - 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 dimension_labels is None: - if array.ndim == 3: - dimension_labels = ['angle' , - 'horizontal' , - 'vertical'] + self.geometry = None + if array is None: + if 'geometry' in kwargs.keys(): + geometry = kwargs['geometry'] + self.geometry = geometry + channels = geometry.channels + horiz = geometry.pixel_num_h + vert = geometry.pixel_num_v + angles = geometry.angles + num_of_angles = numpy.shape(angles)[0] + + + if channels > 1: + if vert > 1: + shape = (channels, num_of_angles , vert, horiz) + dim_labels = ['channel' , ' angle' , + 'vertical' , 'horizontal'] + else: + shape = (channels , num_of_angles, horiz) + dim_labels = ['channel' , 'angle' , + 'horizontal'] else: - dimension_labels = ['angle' , - 'horizontal'] - DataSet.__init__(self, array, deep_copy, dimension_labels, **kwargs) + if vert > 1: + shape = (num_of_angles, vert, horiz) + dim_labels = ['angles' , 'vertical' , + 'horizontal'] + else: + shape = (num_of_angles, horiz) + dim_labels = ['angles' , + 'horizontal'] + + array = numpy.zeros( shape , dtype=numpy.float32) + super(SinogramData, self).__init__(array, deep_copy, + dim_labels, **kwargs) + else: + 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 or \ + array.number_of_dimensions == 4): + raise ValueError('Number of dimensions are not 2 or 3 or 4: {0}'\ + .format(array.number_of_dimensions)) + + #DataSet.__init__(self, array.as_array(), deep_copy, + # array.dimension_labels, **kwargs) + super(SinogramData, self).__init__(array.as_array(), deep_copy, + array.dimension_labels, **kwargs) + elif type(array) == numpy.ndarray: + if not ( array.ndim == 2 or array.ndim == 3 or array.ndim == 4 ): + raise ValueError( + 'Number of dimensions are not 2 or 3 or 4 : {0}'\ + .format(array.ndim)) + + if dimension_labels is None: + if array.ndim == 4: + dimension_labels = ['channel' ,'angle' , 'vertical' , + 'horizontal'] + elif array.ndim == 3: + dimension_labels = ['angle' , 'vertical' , + 'horizontal'] + else: + dimension_labels = ['angle' , + 'horizontal'] + + #DataSet.__init__(self, array, deep_copy, dimension_labels, **kwargs) + super(SinogramData, self).__init__(array, deep_copy, + dimension_labels, **kwargs) + class DataSetProcessor(object): '''Defines a generic DataSet processor @@ -772,5 +868,31 @@ if __name__ == '__main__': s = numpy.reshape(numpy.asarray(s), (3,4,4)) sino = SinogramData( s ) - + shape = (4,3,2) + a = [i for i in range(2*3*4)] + a = numpy.asarray(a) + a = numpy.reshape(a, shape) + print (numpy.shape(a)) + ds = DataSet(a, True, ['X', 'Y','Z']) + # this means that I expect the X to be of length 2 , + # y of length 3 and z of length 4 + subset = ['Y' ,'Z'] + b0 = ds.subset( subset ) + print ("shape b 3,2? {0}".format(numpy.shape(b0.as_array()))) + # expectation on b is that it is + # 3x2 cut at z = 0 + + subset = ['X' ,'Y'] + b1 = ds.subset( subset , Z=1) + print ("shape b 2,3? {0}".format(numpy.shape(b1.as_array()))) + + + # create VolumeData from geometry + vgeometry = geoms.VolumeGeometry(voxel_num_x=2, voxel_num_y=3, channels=2) + vol = VolumeData(geometry=vgeometry) + + sgeometry = geoms.SinogramGeometry(dimension=2, angles=numpy.linspace(0, 180, num=20), + geom_type='parallel', pixel_num_v=3, + pixel_num_h=5 , channels=2) + sino = SinogramData(geometry=sgeometry)
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/processors.py b/Wrappers/Python/ccpi/processors.py index 09e7229..c509a4e 100755 --- a/Wrappers/Python/ccpi/processors.py +++ b/Wrappers/Python/ccpi/processors.py @@ -32,19 +32,23 @@ class Normalizer(DataSetProcessor): Input: SinogramDataSet
Parameter: 2D projection with flat field (or stack)
2D projection with dark field (or stack)
- Output: SinogramDataSet
+ Output: SinogramDataSetn
'''
- def __init__(self):
+ def __init__(self, flat_field = None, dark_field = None, tolerance = 1e-5):
kwargs = {
- 'flat_field' :None,
- 'dark_field' :None,
+ 'flat_field' : None,
+ 'dark_field' : None,
# very small number. Used when there is a division by zero
- 'tolerance' : 1e-5
+ 'tolerance' : tolerance
}
#DataSetProcessor.__init__(self, **kwargs)
super(Normalizer, self).__init__(**kwargs)
+ if not flat_field is None:
+ self.setFlatField(flat_field)
+ if not dark_field is None:
+ self.setDarkField(dark_field)
def checkInput(self, dataset):
if dataset.number_of_dimensions == 3:
diff --git a/Wrappers/Python/ccpi/reconstruction/astra_ops.py b/Wrappers/Python/ccpi/reconstruction/astra_ops.py index 6af7af7..1346f22 100644 --- a/Wrappers/Python/ccpi/reconstruction/astra_ops.py +++ b/Wrappers/Python/ccpi/reconstruction/astra_ops.py @@ -19,7 +19,7 @@ from ccpi.reconstruction.ops import Operator import numpy from ccpi.framework import SinogramData, VolumeData from ccpi.reconstruction.ops import PowerMethodNonsquare -from ccpi.astra.astra_processors import AstraForwardProjector, AstraBackProjector +from ccpi.astra.astra_processors import * class AstraProjectorSimple(Operator): """ASTRA projector modified to use DataSet and geometry.""" @@ -66,6 +66,55 @@ class AstraProjectorSimple(Operator): 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..a887c0c 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=0, + voxel_num_y=0, + voxel_num_z=0, + voxel_size_x=1, + voxel_size_y=1, + voxel_size_z=1, + 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 @@ -35,24 +37,31 @@ class VolumeGeometry: return self.center_y + 0.5*self.voxel_num_y*self.voxel_size_y def getMinZ(self): - return self.center_z - 0.5*self.voxel_num_z*self.voxel_size_z + if not voxel_num_z == 0: + return self.center_z - 0.5*self.voxel_num_z*self.voxel_size_z + else: + return 0 def getMaxZ(self): - return self.center_z + 0.5*self.voxel_num_z*self.voxel_size_z + if not voxel_num_z == 0: + return self.center_z + 0.5*self.voxel_num_z*self.voxel_size_z + else: + return 0 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=0, + pixel_size_h=1, + pixel_num_v=0, + pixel_size_v=1, + dist_source_center=None, + dist_center_detector=None, + channels=1 ): """ General inputs for standard type projection geometries @@ -91,6 +100,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/setup.py b/Wrappers/Python/setup.py index bbd210f..3f83081 100644 --- a/Wrappers/Python/setup.py +++ b/Wrappers/Python/setup.py @@ -35,7 +35,7 @@ if cil_version == '': setup( name="ccpi-common", version=cil_version, - packages=['ccpi' , 'ccpi.reconstruction'], + packages=['ccpi' , 'ccpi.reconstruction' , 'ccpi.astra'], # Project uses reStructuredText, so ensure that the docutils get # installed or upgraded on the target machine diff --git a/Wrappers/Python/test/DemoRecIP.py b/Wrappers/Python/wip/DemoRecIP.py index 483a202..483a202 100644 --- a/Wrappers/Python/test/DemoRecIP.py +++ b/Wrappers/Python/wip/DemoRecIP.py diff --git a/Wrappers/Python/test/regularizers.py b/Wrappers/Python/wip/regularizers.py index 04ac3aa..04ac3aa 100644 --- a/Wrappers/Python/test/regularizers.py +++ b/Wrappers/Python/wip/regularizers.py diff --git a/Wrappers/Python/test/simple_demo.py b/Wrappers/Python/wip/simple_demo.py index 0bbb687..7c5f43e 100644 --- a/Wrappers/Python/test/simple_demo.py +++ b/Wrappers/Python/wip/simple_demo.py @@ -1,34 +1,30 @@ +#import sys +#sys.path.append("..") -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 * +from ccpi.framework import VolumeData +from ccpi.reconstruction.algs import FISTA +from ccpi.reconstruction.funcs import Norm2sq, Norm1 +from ccpi.reconstruction.astra_ops import AstraProjectorSimple +from ccpi.reconstruction.geoms import VolumeGeometry, SinogramGeometry import numpy as np import matplotlib.pyplot as plt -test_case = 2 # 1=parallel2D, 2=cone2D +test_case = 1 # 1=parallel2D, 2=cone2D # Set up phantom N = 128 -x = np.zeros((N,N)) +vg = VolumeGeometry(voxel_num_x=N,voxel_num_y=N) +Phantom = VolumeData(geometry=vg) + +x = Phantom.as_array() x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 1.0 x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 2.0 plt.imshow(x) plt.show() -vg = VolumeGeometry(N,N,None, 1,1,None) - -Phantom = VolumeData(x,geometry=vg) - # Set up measurement geometry angles_num = 20; # angles number diff --git a/Wrappers/Python/wip/simple_mc_demo.py b/Wrappers/Python/wip/simple_mc_demo.py new file mode 100755 index 0000000..0d976d7 --- /dev/null +++ b/Wrappers/Python/wip/simple_mc_demo.py @@ -0,0 +1,143 @@ +#import sys +#sys.path.append("..") + +from ccpi.framework import VolumeData, SinogramData +from ccpi.reconstruction.algs import FISTA +from ccpi.reconstruction.funcs import Norm2sq, Norm1 +from ccpi.reconstruction.astra_ops import AstraProjectorMC +from ccpi.reconstruction.geoms import VolumeGeometry, SinogramGeometry + +import numpy +import matplotlib.pyplot as plt + +test_case = 1 # 1=parallel2D, 2=cone2D + +# Set up phantom +N = 128 +M = 100 +numchannels = 3 + +vg = VolumeGeometry(voxel_num_x=N,voxel_num_y=N,channels=numchannels) +Phantom = VolumeData(geometry=vg) + +x = Phantom.as_array() +x[0 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4) ] = 1.0 +x[0 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 2.0 + +x[1 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4) ] = 0.7 +x[1 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 1.2 + +x[2 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4) ] = 1.5 +x[2 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 2.2 + +#x = numpy.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[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 = numpy.linspace(0,numpy.pi,angles_num,endpoint=False) +elif test_case==2: + angles = numpy.linspace(0,2*numpy.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, 'cpu') + + +# 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()[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()[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(numpy.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()[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()[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 |