diff options
Diffstat (limited to 'Wrappers')
| -rw-r--r-- | Wrappers/Python/ccpi/astra/astra_ops.py | 6 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/framework.py | 278 | ||||
| -rwxr-xr-x | Wrappers/Python/ccpi/processors.py | 357 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/reconstruction/algs.py | 89 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/reconstruction/ops.py | 85 | ||||
| -rwxr-xr-x | Wrappers/Python/wip/simple_demo_astra.py | 189 | ||||
| -rwxr-xr-x | Wrappers/Python/wip/simple_demo_ccpi.py | 270 |
7 files changed, 1183 insertions, 91 deletions
diff --git a/Wrappers/Python/ccpi/astra/astra_ops.py b/Wrappers/Python/ccpi/astra/astra_ops.py index 5f25eec..d8c2f4b 100644 --- a/Wrappers/Python/ccpi/astra/astra_ops.py +++ b/Wrappers/Python/ccpi/astra/astra_ops.py @@ -18,7 +18,7 @@ from ccpi.reconstruction.ops import Operator import numpy import astra -from ccpi.framework import AcquisitionData, ImageData +from ccpi.framework import AcquisitionData, ImageData, DataContainer from ccpi.reconstruction.ops import PowerMethodNonsquare from ccpi.astra.astra_processors import AstraForwardProjector, AstraBackProjector, \ AstraForwardProjectorMC, AstraBackProjectorMC @@ -68,6 +68,10 @@ class AstraProjectorSimple(Operator): 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""" diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py index 877becf..15dbe30 100644 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework.py @@ -83,6 +83,20 @@ class ImageGeometry: else: return 0 + def clone(self): + '''returns a copy of ImageGeometry''' + return ImageGeometry( + self.voxel_num_x, + self.voxel_num_y, + self.voxel_num_z, + self.voxel_size_x, + self.voxel_size_y, + self.voxel_size_z, + self.center_x, + self.center_y, + self.center_z, + self.channels) + class AcquisitionGeometry: @@ -137,6 +151,31 @@ class AcquisitionGeometry: self.pixel_size_v = pixel_size_v self.channels = channels + + def clone(self): + '''returns a copy of the AcquisitionGeometry''' + return AcquisitionGeometry(self.geom_type, + self.dimension, + self.angles, + self.pixel_num_h, + self.pixel_size_h, + self.pixel_num_v, + self.pixel_size_v, + self.dist_source_center, + self.dist_center_detector, + self.channels) + + def __str__ (self): + repres = "" + repres += "Number of dimensions: {0}\n".format(self.dimension) + repres += "angles: {0}\n".format(len(self.angles)) + repres += "voxel_num : h{0},v{1}\n".format(self.pixel_num_h, self.pixel_num_v) + repres += "voxel size: h{0},v{1}\n".format(self.pixel_size_h, self.pixel_size_v) + repres += "geometry type: {0}\n".format(self.geom_type) + repres += "distance source-detector: {0}\n".format(self.dist_source_center) + repres += "distance center-detector: {0}\n".format(self.dist_source_center) + repres += "number of channels: {0}\n".format(self.channels) + return repres @@ -178,6 +217,25 @@ class DataContainer(object): # assume it is parallel beam pass + def get_dimension_size(self, dimension_label): + if dimension_label in self.dimension_labels.values(): + acq_size = -1 + for k,v in self.dimension_labels.items(): + if v == dimension_label: + acq_size = self.shape[k] + return acq_size + else: + raise ValueError('Unknown dimension {0}. Should be one of'.format(dimension_label, + self.dimension_labels)) + def get_dimension_axis(self, dimension_label): + if dimension_label in self.dimension_labels.values(): + for k,v in self.dimension_labels.items(): + if v == dimension_label: + return k + else: + raise ValueError('Unknown dimension {0}. Should be one of'.format(dimension_label, + self.dimension_labels.values())) + def as_array(self, dimensions=None): '''Returns the DataContainer as Numpy Array @@ -194,7 +252,15 @@ class DataContainer(object): '''Creates a DataContainer containing a subset of self according to the labels in dimensions''' if dimensions is None: - return self.array.copy() + if kw == {}: + return self.array.copy() + else: + reduced_dims = [v for k,v in self.dimension_labels.items()] + for dim_l, dim_v in kw.items(): + for k,v in self.dimension_labels.items(): + if v == dim_l: + reduced_dims.pop(k) + return self.subset(dimensions=reduced_dims, **kw) else: # check that all the requested dimensions are in the array # this is done by checking the dimension_labels @@ -258,14 +324,31 @@ class DataContainer(object): return type(self)(cleaned , True, dimensions) - def fill(self, array): + def fill(self, array, **dimension): '''fills the internal numpy array with the one provided''' - if numpy.shape(array) != numpy.shape(self.array): - raise ValueError('Cannot fill with the provided array.' + \ - 'Expecting {0} got {1}'.format( - numpy.shape(self.array), - numpy.shape(array))) - self.array = array[:] + if dimension == {}: + if numpy.shape(array) != numpy.shape(self.array): + raise ValueError('Cannot fill with the provided array.' + \ + 'Expecting {0} got {1}'.format( + numpy.shape(self.array), + numpy.shape(array))) + self.array = array[:] + else: + + command = 'self.array[' + i = 0 + for k,v in self.dimension_labels.items(): + for dim_label, dim_value in dimension.items(): + if dim_label == v: + command = command + str(dim_value) + else: + command = command + ":" + if i < self.number_of_dimensions -1: + command = command + ',' + i += 1 + command = command + "] = array[:]" + exec(command) + def check_dimensions(self, other): return self.shape == other.shape @@ -301,7 +384,7 @@ class DataContainer(object): dimension_labels=self.dimension_labels, geometry=self.geometry) else: - raise ValueError('Wrong shape: {0} and {1}'.format(self.shape, + raise ValueError('__sub__ Wrong shape: {0} and {1}'.format(self.shape, other.shape)) elif isinstance(other, (int, float, complex)): return type(self)(self.as_array() - other, @@ -325,7 +408,7 @@ class DataContainer(object): dimension_labels=self.dimension_labels, geometry=self.geometry) else: - raise ValueError('Wrong shape: {0} and {1}'.format(self.shape, + raise ValueError('__div__ Wrong shape: {0} and {1}'.format(self.shape, other.shape)) elif isinstance(other, (int, float, complex)): return type(self)(self.as_array() / other, @@ -346,7 +429,7 @@ class DataContainer(object): dimension_labels=self.dimension_labels, geometry=self.geometry) else: - raise ValueError('Wrong shape: {0} and {1}'.format(self.shape, + raise ValueError('__pow__ Wrong shape: {0} and {1}'.format(self.shape, other.shape)) elif isinstance(other, (int, float, complex)): return type(self)(self.as_array() ** other, @@ -354,7 +437,7 @@ class DataContainer(object): dimension_labels=self.dimension_labels, geometry=self.geometry) else: - raise TypeError('Cannot {0} DataContainer with {1}'.format("power" , + raise TypeError('pow: Cannot {0} DataContainer with {1}'.format("power" , type(other))) # __pow__ @@ -367,7 +450,7 @@ class DataContainer(object): dimension_labels=self.dimension_labels, geometry=self.geometry) else: - raise ValueError('Wrong shape: {0} and {1}'.format(self.shape, + raise ValueError('*:Wrong shape: {0} and {1}'.format(self.shape, other.shape)) elif isinstance(other, (int, float, complex)): return type(self)(self.as_array() * other, @@ -464,12 +547,13 @@ class DataContainer(object): return self / other # __idiv__ - def __str__ (self): + def __str__ (self, representation=False): repres = "" 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: \n{0}\n".format(self.array) + if representation: + repres += "Representation: \n{0}\n".format(self.array) return repres def clone(self): @@ -479,6 +563,40 @@ class DataContainer(object): dimension_labels=self.dimension_labels, deep_copy=True, geometry=self.geometry ) + + def get_data_axes_order(self,new_order=None): + '''returns the axes label of self as a list + + if new_order is None returns the labels of the axes as a sorted-by-key list + if new_order is a list of length number_of_dimensions, returns a list + with the indices of the axes in new_order with respect to those in + self.dimension_labels: i.e. + self.dimension_labels = {0:'horizontal',1:'vertical'} + new_order = ['vertical','horizontal'] + returns [1,0] + ''' + if new_order is None: + + axes_order = [i for i in range(len(self.shape))] + for k,v in self.dimension_labels.items(): + axes_order[k] = v + return axes_order + else: + if len(new_order) == self.number_of_dimensions: + axes_order = [i for i in range(self.number_of_dimensions)] + + for i in range(len(self.shape)): + found = False + for k,v in self.dimension_labels.items(): + if new_order[i] == v: + axes_order[i] = k + found = True + if not found: + raise ValueError('Axis label {0} not found.'.format(new_order[i])) + return axes_order + else: + raise ValueError('Expecting {0} axes, got {2}'\ + .format(len(self.shape),len(new_order))) @@ -501,29 +619,45 @@ class ImageData(DataContainer): 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 channels > 1: - if vert > 1: - shape = (channels, vert, horiz_y, horiz_x) - dim_labels = ['channel' ,'vertical' , 'horizontal_y' , - 'horizontal_x'] + if dimension_labels is None: + 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: - shape = (channels , horiz_y, horiz_x) - dim_labels = ['channel' , 'horizontal_y' , - 'horizontal_x'] + 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'] + dimension_labels = dim_labels else: - 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'] - + shape = [] + for dim in dimension_labels: + if dim == 'channel': + shape.append(channels) + elif dim == 'horizontal_y': + shape.append(horiz_y) + elif dim == 'vertical': + shape.append(vert) + elif dim == 'horizontal_x': + shape.append(horiz_x) + if len(shape) != len(dimension_labels): + raise ValueError('Missing {0} axes'.format( + len(dimension_labels) - len(shape))) + shape = tuple(shape) + array = numpy.zeros( shape , dtype=numpy.float32) super(ImageData, self).__init__(array, deep_copy, - dim_labels, **kwargs) + dimension_labels, **kwargs) else: raise ValueError('Please pass either a DataContainer, ' +\ @@ -571,6 +705,12 @@ class ImageData(DataContainer): if key == 'spacing' : self.spacing = value + def subset(self, dimensions=None, **kw): + out = super(ImageData, self).subset(dimensions, **kw) + #out.geometry = self.recalculate_geometry(dimensions , **kw) + out.geometry = self.geometry + return out + class AcquisitionData(DataContainer): '''DataContainer for holding 2D or 3D sinogram''' @@ -582,37 +722,54 @@ class AcquisitionData(DataContainer): self.geometry = None if array is None: if 'geometry' in kwargs.keys(): - geometry = kwargs['geometry'] + geometry = kwargs['geometry'] self.geometry = geometry - channels = geometry.channels - horiz = geometry.pixel_num_h - vert = geometry.pixel_num_v - angles = geometry.angles + 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'] + if dimension_labels is None: + 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: - shape = (channels , num_of_angles, horiz) - dim_labels = ['channel' , 'angle' , - 'horizontal'] + if vert > 1: + shape = (num_of_angles, vert, horiz) + dim_labels = ['angle' , 'vertical' , + 'horizontal'] + else: + shape = (num_of_angles, horiz) + dim_labels = ['angle' , + 'horizontal'] + + dimension_labels = dim_labels else: - if vert > 1: - shape = (num_of_angles, vert, horiz) - dim_labels = ['angles' , 'vertical' , - 'horizontal'] - else: - shape = (num_of_angles, horiz) - dim_labels = ['angles' , - 'horizontal'] - + shape = [] + for dim in dimension_labels: + if dim == 'channel': + shape.append(channels) + elif dim == 'angle': + shape.append(num_of_angles) + elif dim == 'vertical': + shape.append(vert) + elif dim == 'horizontal': + shape.append(horiz) + if len(shape) != len(dimension_labels): + raise ValueError('Missing {0} axes'.format( + len(dimension_labels) - len(shape))) + shape = tuple(shape) + array = numpy.zeros( shape , dtype=numpy.float32) super(AcquisitionData, self).__init__(array, deep_copy, - dim_labels, **kwargs) + dimension_labels, **kwargs) else: if type(array) == DataContainer: @@ -732,6 +889,9 @@ class DataSetProcessor(object): def process(self): raise NotImplementedError('process must be implemented') + + + class DataSetProcessor23D(DataSetProcessor): '''Regularizers DataSetProcessor @@ -945,7 +1105,7 @@ if __name__ == '__main__': # create VolumeData from geometry vgeometry = ImageGeometry(voxel_num_x=2, voxel_num_y=3, channels=2) - vol = VolumeData(geometry=vgeometry) + vol = ImageData(geometry=vgeometry) sgeometry = AcquisitionGeometry(dimension=2, angles=numpy.linspace(0, 180, num=20), geom_type='parallel', pixel_num_v=3, diff --git a/Wrappers/Python/ccpi/processors.py b/Wrappers/Python/ccpi/processors.py index 3010053..d98ef12 100755 --- a/Wrappers/Python/ccpi/processors.py +++ b/Wrappers/Python/ccpi/processors.py @@ -17,12 +17,16 @@ # See the License for the specific language governing permissions and
# limitations under the License
-from ccpi.framework import DataSetProcessor, DataSet, AcquisitionData,\
- AcquisitionGeometry
+from ccpi.framework import DataSetProcessor, DataContainer, AcquisitionData,\
+ AcquisitionGeometry, ImageGeometry, ImageData
+from ccpi.reconstruction.parallelbeam import alg as pbalg
import numpy
import h5py
from scipy import ndimage
+import matplotlib.pyplot as plt
+
+
class Normalizer(DataSetProcessor):
'''Normalization based on flat and dark
@@ -130,7 +134,8 @@ class CenterOfRotationFinder(DataSetProcessor): if dataset.geometry.geom_type == 'parallel':
return True
else:
- raise ValueError('This algorithm is suitable only for parallel beam geometry')
+ raise ValueError('{0} is suitable only for parallel beam geometry'\
+ .format(self.__class__.__name__))
else:
raise ValueError("Expected input dimensions is 3, got {0}"\
.format(dataset.number_of_dimensions))
@@ -385,6 +390,275 @@ class CenterOfRotationFinder(DataSetProcessor): return cor
+class CCPiForwardProjector(DataSetProcessor):
+ '''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':
+ #int msize = ndarray_volume.shape(0) > ndarray_volume.shape(1) ? ndarray_volume.shape(0) : ndarray_volume.shape(1);
+ #int detector_width = msize;
+ # detector_width is the max between the shape[0] and shape[1]
+
+
+ #double rotation_center = (double)detector_width/2.;
+ #int detector_height = ndarray_volume.shape(2);
+
+ #int number_of_projections = ndarray_angles.shape(0);
+
+ ##numpy_3d pixels(reinterpret_cast<float*>(ndarray_volume.get_data()),
+ #boost::extents[number_of_projections][detector_height][detector_width]);
+
+ 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(DataSetProcessor):
+ '''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:
+ #number_of_projections][detector_height][detector_width
+
+ 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
+ #print (center_of_rotation)
+ 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(DataSetProcessor):
+ '''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 * cor
+
+ 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
+
+#class FiniteDifferentiator(DataSetProcessor):
+# def __init__(self):
+# kwargs = {
+#
+# }
+#
+# super(FiniteDifferentiator, self).__init__(**kwargs)
+#
+# def check_input(self, dataset):
+# return True
+#
+# def process(self):
+# axis = 0
+# d1 = numpy.diff(x,n=1,axis=axis)
+# d1 = numpy.resize(d1, numpy.shape(x))
+
+
+
def loadNexus(filename):
'''Load a dataset stored in a NeXuS file (HDF5)'''
###########################################################################
@@ -437,35 +711,94 @@ if __name__ == '__main__': pixel_num_h=numpy.shape(proj)[2],
pixel_num_v=numpy.shape(proj)[1],
)
- sino = AcquisitionData( proj , geometry=parallelbeam)
+
+ dim_labels = ['angles' , 'vertical' , 'horizontal']
+ sino = AcquisitionData( proj , geometry=parallelbeam,
+ dimension_labels=dim_labels)
normalizer = Normalizer()
normalizer.set_input(sino)
normalizer.set_flat_field(flat)
normalizer.set_dark_field(dark)
norm = normalizer.get_output()
- print ("Processor min {0} max {1}".format(norm.as_array().min(), norm.as_array().max()))
+ #print ("Processor min {0} max {1}".format(norm.as_array().min(), norm.as_array().max()))
- norm1 = numpy.asarray(
- [Normalizer.normalize_projection( p, flat, dark, 1e-5 )
- for p in proj]
- )
+ #norm1 = numpy.asarray(
+ # [Normalizer.normalize_projection( p, flat, dark, 1e-5 )
+ # for p in proj]
+ # )
- print ("Numpy min {0} max {1}".format(norm1.min(), norm1.max()))
+ #print ("Numpy min {0} max {1}".format(norm1.min(), norm1.max()))
cor_finder = CenterOfRotationFinder()
cor_finder.set_input(sino)
cor = cor_finder.get_output()
print ("center of rotation {0} == 86.25?".format(cor))
+
+ padder = AcquisitionDataPadder(center_of_rotation=cor,
+ pad_value=0.75,
+ acquisition_geometry=normalizer.get_output().geometry)
+ #padder.set_input(normalizer.get_output())
+ padder.set_input_processor(normalizer)
+
+ #print ("padder ", padder.get_output())
+
+ volume_geometry = ImageGeometry()
+
+ back = CCPiBackwardProjector(acquisition_geometry=parallelbeam,
+ image_geometry=volume_geometry)
+ back.set_input_processor(padder)
+ #back.set_input(padder.get_output())
+ #print (back.image_geometry)
+ forw = CCPiForwardProjector(acquisition_geometry=parallelbeam,
+ image_geometry=volume_geometry)
+ forw.set_input_processor(back)
+
+
+ out = padder.get_output()
+ fig = plt.figure()
+ # projections row
+ a = fig.add_subplot(1,4,1)
+ a.imshow(norm.array[80])
+ a.set_title('orig')
+ a = fig.add_subplot(1,4,2)
+ a.imshow(out.array[80])
+ a.set_title('padded')
+
+ a = fig.add_subplot(1,4,3)
+ a.imshow(back.get_output().as_array()[15])
+ a.set_title('back')
+
+ a = fig.add_subplot(1,4,4)
+ a.imshow(forw.get_output().as_array()[80])
+ a.set_title('forw')
+
+
+ plt.show()
+
+# back = CCPiBackwardProjector(acquisition_geometry=parallelbeam,
+# image_geometry=volume_geometry)
+
+# int msize = ndarray_volume.shape(0) > ndarray_volume.shape(1) ? ndarray_volume.shape(0) : ndarray_volume.shape(1);#
+# int detector_width = msize;
+
+# double rotation_center = (double)detector_width/2.;
+# int detector_height = ndarray_volume.shape(2);
+# int number_of_projections = ndarray_angles.shape(0);
+# std::cout << "pb_forward_project rotation_center " << rotation_center << std::endl;
+# boost::extents[number_of_projections][detector_height][detector_width]);
+
+
conebeam = AcquisitionGeometry('cone', '3D' ,
angles=angles,
pixel_num_h=numpy.shape(proj)[2],
pixel_num_v=numpy.shape(proj)[1],
)
- sino = AcquisitionData( proj , geometry=conebeam)
+
try:
- cor_finder.set_input(sino)
+ sino2 = AcquisitionData( proj , geometry=conebeam)
+ cor_finder.set_input(sino2)
cor = cor_finder.get_output()
except ValueError as err:
print (err)
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/reconstruction/algs.py b/Wrappers/Python/ccpi/reconstruction/algs.py index ec52fee..46865a8 100644 --- a/Wrappers/Python/ccpi/reconstruction/algs.py +++ b/Wrappers/Python/ccpi/reconstruction/algs.py @@ -23,14 +23,36 @@ import time from ccpi.reconstruction.funcs import BaseFunction def FISTA(x_init, f=None, g=None, opt=None): - + '''Fast Iterative Shrinkage-Thresholding Algorithm + + Beck, A. and Teboulle, M., 2009. A fast iterative shrinkage-thresholding + algorithm for linear inverse problems. + SIAM journal on imaging sciences,2(1), pp.183-202. + + Parameters: + x_init: initial guess + f: data fidelity + g: regularizer + h: + opt: additional algorithm + ''' # default inputs if f is None: f = BaseFunction() if g is None: g = BaseFunction() - if opt is None: opt = {'tol': 1e-4, 'iter': 1000} - + # algorithmic parameters - tol = opt['tol'] + if opt is None: + opt = {'tol': 1e-4, 'iter': 1000} + else: + try: + max_iter = opt['iter'] + except KeyError as ke: + opt[ke] = 1000 + try: + opt['tol'] = 1000 + except KeyError as ke: + opt[ke] = 1e-4 + tol = opt['tol'] max_iter = opt['iter'] # initialization @@ -76,15 +98,33 @@ def FISTA(x_init, f=None, g=None, opt=None): return x, it, timing, criter def FBPD(x_init, f=None, g=None, h=None, opt=None): - + '''FBPD Algorithm + + Parameters: + x_init: initial guess + f: constraint + g: data fidelity + h: regularizer + opt: additional algorithm + ''' # default inputs if f is None: f = BaseFunction() if g is None: g = BaseFunction() if h is None: h = BaseFunction() - if opt is None: opt = {'tol': 1e-4, 'iter': 1000} - + # algorithmic parameters - tol = opt['tol'] + if opt is None: + opt = {'tol': 1e-4, 'iter': 1000} + else: + try: + max_iter = opt['iter'] + except KeyError as ke: + opt[ke] = 1000 + try: + opt['tol'] = 1000 + except KeyError as ke: + opt[ke] = 1e-4 + tol = opt['tol'] max_iter = opt['iter'] # step-sizes @@ -127,13 +167,34 @@ def FBPD(x_init, f=None, g=None, h=None, opt=None): return x, it, timing, criter -def CGLS(A,b,max_iter,x_init): - '''Conjugate Gradient Least Squares algorithm''' +def CGLS(x_init, operator , data , opt=None): + '''Conjugate Gradient Least Squares algorithm + + Parameters: + x_init: initial guess + operator: operator for forward/backward projections + data: data to operate on + opt: additional algorithm + ''' + + if opt is None: + opt = {'tol': 1e-4, 'iter': 1000} + else: + try: + max_iter = opt['iter'] + except KeyError as ke: + opt[ke] = 1000 + try: + opt['tol'] = 1000 + except KeyError as ke: + opt[ke] = 1e-4 + tol = opt['tol'] + max_iter = opt['iter'] - r = b.clone() + r = data.clone() x = x_init.clone() - d = A.adjoint(r) + d = operator.adjoint(r) normr2 = (d**2).sum() @@ -145,11 +206,11 @@ def CGLS(A,b,max_iter,x_init): t = time.time() - Ad = A.direct(d) + Ad = operator.direct(d) alpha = normr2/( (Ad**2).sum() ) x = x + alpha*d r = r - alpha*Ad - s = A.adjoint(r) + s = operator.adjoint(r) normr2_new = (s**2).sum() beta = normr2_new/normr2 diff --git a/Wrappers/Python/ccpi/reconstruction/ops.py b/Wrappers/Python/ccpi/reconstruction/ops.py index d6f31eb..fa30b0c 100644 --- a/Wrappers/Python/ccpi/reconstruction/ops.py +++ b/Wrappers/Python/ccpi/reconstruction/ops.py @@ -19,7 +19,8 @@ import numpy from scipy.sparse.linalg import svds -from ccpi.framework import DataContainer +from ccpi.framework import DataContainer, ImageGeometry , ImageData +from ccpi.processors import CCPiBackwardProjector, CCPiForwardProjector # Maybe operators need to know what types they take as inputs/outputs # to not just use generic DataContainer @@ -110,7 +111,13 @@ class FiniteDiff2D(Operator): def adjoint(self,x): '''Backward differences, Newumann BC.''' - Nrows, Ncols, Nchannels = x.as_array().shape + #Nrows, Ncols, Nchannels = x.as_array().shape + print (x) + Nrows = x.get_dimension_size('horizontal_x') + Ncols = x.get_dimension_size('horizontal_x') + Nchannels = 1 + if len(x.shape) == 4: + Nchannels = x.get_dimension_size('channel') zer = numpy.zeros((Nrows,1)) xxx = x.as_array()[:,:-1,0] h = numpy.concatenate((zer,xxx), 1) - numpy.concatenate((xxx,zer), 1) @@ -129,13 +136,27 @@ class FiniteDiff2D(Operator): def PowerMethodNonsquare(op,numiters): # Initialise random - inputsize = op.size()[1] - x0 = DataContainer(numpy.random.randn(inputsize[0],inputsize[1])) + # Jakob's + #inputsize = op.size()[1] + #x0 = ImageContainer(numpy.random.randn(*inputsize) + # Edo's + #vg = ImageGeometry(voxel_num_x=inputsize[0], + # voxel_num_y=inputsize[1], + # voxel_num_z=inputsize[2]) + # + #x0 = ImageData(geometry = vg, dimension_labels=['vertical','horizontal_y','horizontal_x']) + #print (x0) + #x0.fill(numpy.random.randn(*x0.shape)) + + x0 = op.create_image_data() + s = numpy.zeros(numiters) # Loop for it in numpy.arange(numiters): x1 = op.adjoint(op.direct(x0)) x1norm = numpy.sqrt((x1**2).sum()) + #print ("x0 **********" ,x0) + #print ("x1 **********" ,x1) s[it] = (x1*x0).sum() / (x0*x0).sum() x0 = (1.0/x1norm)*x1 return numpy.sqrt(s[-1]), numpy.sqrt(s), x0 @@ -150,4 +171,58 @@ def PowerMethodNonsquare(op,numiters): # x1norm = np.sqrt(np.sum(np.dot(x1,x1))) # s[it] = np.dot(x1,x0) / np.dot(x1,x0) # x0 = (1.0/x1norm)*x1 -# return s, x0
\ No newline at end of file +# return s, x0 + +class CCPiProjectorSimple(Operator): + """ASTRA projector modified to use DataSet and geometry.""" + def __init__(self, geomv, geomp): + super(CCPiProjectorSimple, self).__init__() + + # Store volume and sinogram geometries. + self.acquisition_geometry = geomp + self.volume_geometry = geomv + + self.fp = CCPiForwardProjector(image_geometry=geomv, + acquisition_geometry=geomp, + output_axes_order=['angle','vertical','horizontal']) + + self.bp = CCPiBackwardProjector(image_geometry=geomv, + acquisition_geometry=geomp, + output_axes_order=['horizontal_x','horizontal_y','vertical']) + + # Initialise empty for singular value. + self.s1 = None + + def direct(self, image_data): + self.fp.set_input(image_data) + out = self.fp.get_output() + return out + + def adjoint(self, acquisition_data): + self.bp.set_input(acquisition_data) + out = self.bp.get_output() + return out + + #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']) + print (x0) + x0.fill(numpy.random.randn(*x0.shape)) + return x0 diff --git a/Wrappers/Python/wip/simple_demo_astra.py b/Wrappers/Python/wip/simple_demo_astra.py new file mode 100755 index 0000000..3365504 --- /dev/null +++ b/Wrappers/Python/wip/simple_demo_astra.py @@ -0,0 +1,189 @@ +#import sys +#sys.path.append("..") + +from ccpi.framework import ImageData , ImageGeometry, AcquisitionGeometry +from ccpi.reconstruction.algs import FISTA, FBPD, CGLS +from ccpi.reconstruction.funcs import Norm2sq, Norm1 , TV2D +from ccpi.astra.astra_ops import AstraProjectorSimple + +import numpy as np +import matplotlib.pyplot as plt + +test_case = 1 # 1=parallel2D, 2=cone2D + +# Set up phantom +N = 128 + + +vg = ImageGeometry(voxel_num_x=N,voxel_num_y=N) +Phantom = ImageData(geometry=vg) + +x = Phantom.as_array() +x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 +x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 + +plt.imshow(x) +plt.show() + +# 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 = AcquisitionGeometry('parallel', + '2D', + angles, + det_num,det_w) +elif test_case==2: + pg = AcquisitionGeometry('cone', + '2D', + angles, + det_num, + det_w, + dist_source_center=SourceOrig, + dist_center_detector=OrigDetec) + +# ASTRA operator using volume and sinogram geometries +Aop = AstraProjectorSimple(vg, pg, 'cpu') + +# Unused old astra projector without geometry +# Aop_old = AstraProjector(det_w, det_num, SourceOrig, +# OrigDetec, angles, +# N,'fanbeam','gpu') + +# Try forward and backprojection +b = Aop.direct(Phantom) +out2 = Aop.adjoint(b) + +plt.imshow(b.array) +plt.show() + +plt.imshow(out2.array) +plt.show() + +# Create least squares object instance with projector and data. +f = Norm2sq(Aop,b,c=0.5) + +# Initial guess +x_init = ImageData(np.zeros(x.shape),geometry=vg) + +# Run FISTA for least squares without regularization +x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None) + +plt.imshow(x_fista0.array) +plt.title('FISTA0') +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) + +plt.imshow(x_fista1.array) +plt.title('FISTA1') +plt.show() + +plt.semilogy(criter1) +plt.show() + +# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm +opt = {'tol': 1e-4, 'iter': 100} +x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt) + +plt.imshow(x_fbpd1.array) +plt.title('FBPD1') +plt.show() + +plt.semilogy(criter_fbpd1) +plt.show() + +# Now FBPD for least squares plus TV +#lamtv = 1 +#gtv = TV2D(lamtv) + +#x_fbpdtv, it_fbpdtv, timing_fbpdtv, criter_fbpdtv = FBPD(x_init,None,f,gtv,opt=opt) + +#plt.imshow(x_fbpdtv.array) +#plt.show() + +#plt.semilogy(criter_fbpdtv) +#plt.show() + + +# Run CGLS, which should agree with the FISTA0 +x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt ) + +plt.imshow(x_CGLS.array) +plt.title('CGLS') +#plt.title('CGLS recon, compare FISTA0') +plt.show() + +plt.semilogy(criter_CGLS) +plt.title('CGLS criterion') +plt.show() + + +#%% + +clims = (0,1) +cols = 3 +rows = 2 +current = 1 +fig = plt.figure() +# projections row +a=fig.add_subplot(rows,cols,current) +a.set_title('phantom {0}'.format(np.shape(Phantom.as_array()))) + +imgplot = plt.imshow(Phantom.as_array(),vmin=clims[0],vmax=clims[1]) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('FISTA0') +imgplot = plt.imshow(x_fista0.as_array(),vmin=clims[0],vmax=clims[1]) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('FISTA1') +imgplot = plt.imshow(x_fista1.as_array(),vmin=clims[0],vmax=clims[1]) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('FBPD1') +imgplot = plt.imshow(x_fbpd1.as_array(),vmin=clims[0],vmax=clims[1]) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('CGLS') +imgplot = plt.imshow(x_CGLS.as_array(),vmin=clims[0],vmax=clims[1]) + +#current = current + 1 +#a=fig.add_subplot(rows,cols,current) +#a.set_title('FBPD TV') +#imgplot = plt.imshow(x_fbpdtv.as_array(),vmin=clims[0],vmax=clims[1]) + +fig = plt.figure() +# projections row +b=fig.add_subplot(1,1,1) +b.set_title('criteria') +imgplot = plt.loglog(criter0 , label='FISTA0') +imgplot = plt.loglog(criter1 , label='FISTA1') +imgplot = plt.loglog(criter_fbpd1, label='FBPD1') +imgplot = plt.loglog(criter_CGLS, label='CGLS') +#imgplot = plt.loglog(criter_fbpdtv, label='FBPD TV') +b.legend(loc='right') +plt.show() +#%%
\ No newline at end of file diff --git a/Wrappers/Python/wip/simple_demo_ccpi.py b/Wrappers/Python/wip/simple_demo_ccpi.py new file mode 100755 index 0000000..4eeaaae --- /dev/null +++ b/Wrappers/Python/wip/simple_demo_ccpi.py @@ -0,0 +1,270 @@ +#import sys +#sys.path.append("..") + +from ccpi.framework import ImageData , AcquisitionData, ImageGeometry, AcquisitionGeometry +from ccpi.reconstruction.algs import FISTA, FBPD, CGLS +from ccpi.reconstruction.funcs import Norm2sq, Norm1 , TV2D +from ccpi.astra.astra_ops import AstraProjectorSimple +from ccpi.reconstruction.ops import CCPiProjectorSimple +from ccpi.reconstruction.parallelbeam import alg as pbalg +from ccpi.processors import CCPiForwardProjector, CCPiBackwardProjector + +import numpy as np +import matplotlib.pyplot as plt + +test_case = 1 # 1=parallel2D, 2=cone2D, 3=parallel3D + +# Set up phantom +N = 128 +vert = 1 +# Set up measurement geometry +angles_num = 20; # angles number +det_w = 1.0 +det_num = N +SourceOrig = 200 +OrigDetec = 0 + +if test_case==1: + angles = np.linspace(0,np.pi,angles_num,endpoint=False,dtype=np.float32)*180/np.pi + #nangles = angles_num + #angles = np.linspace(0,360, nangles, dtype=np.float32) + +elif test_case==2: + angles = np.linspace(0,2*np.pi,angles_num,endpoint=False) +elif test_case == 3: + angles = np.linspace(0,np.pi,angles_num,endpoint=False) +else: + NotImplemented + + +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'],det_w, + geoms['n_v'], det_w #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 ) + + #print (counter) + counter+=1 + #print (geoms , geoms_i) + 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: + print ("return geoms {0}".format(geoms)) + return geoms + else: + print ("return geoms_i {0}".format(geoms_i)) + return geoms_i + +geoms = setupCCPiGeometries(N,N,vert,angles,0) +#%% +#geoms = {'n_v': 12, 'output_volume_y': 128, 'n_angles': 20, +# 'output_volume_x': 128, 'output_volume_z': 12, 'n_h': 128} +vg = ImageGeometry(voxel_num_x=geoms['output_volume_x'], + voxel_num_y=geoms['output_volume_y'], + voxel_num_z=geoms['output_volume_z']) +Phantom = ImageData(geometry=vg,dimension_labels=['horizontal_x','horizontal_y','vertical']) + 0.1 + + +#x = Phantom.as_array() +i = 0 +while i < geoms['n_v']: + #x = Phantom.subset(vertical=i, dimensions=['horizontal_x','horizontal_y']).array + x = Phantom.subset(vertical=i).array + x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 + x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 0.98 + Phantom.fill(x, vertical=i) + i += 1 + +plt.imshow(Phantom.subset(vertical=0).as_array()) +plt.show() + + + +# Parallelbeam geometry test +if test_case==1: + #Phantom_ccpi = Phantom.subset(['horizontal_x','horizontal_y','vertical']) + #Phantom_ccpi.geometry = vg.clone() + center_of_rotation = Phantom.get_dimension_size('horizontal_x') / 2 + + pg = AcquisitionGeometry('parallel', + '3D', + angles, + geoms['n_h'],det_w, + geoms['n_v'], det_w #2D in 3D is a slice 1 pixel thick + ) +elif test_case==2: + raise NotImplemented('cone beam projector not yet available') + pg = AcquisitionGeometry('cone', + '2D', + angles, + det_num, + det_w, + vert, det_w, #2D in 3D is a slice 1 pixel thick + dist_source_center=SourceOrig, + dist_center_detector=OrigDetec) + +# ASTRA operator using volume and sinogram geometries +#Aop = AstraProjectorSimple(vg, pg, 'cpu') +Cop = CCPiProjectorSimple(vg, pg) + +# Try forward and backprojection +b = Cop.direct(Phantom) +out2 = Cop.adjoint(b) + +#%% +for i in range(b.get_dimension_size('vertical')): + plt.imshow(b.subset(vertical=i).array) + #plt.imshow(Phantom.subset( vertical=i).array) + #plt.imshow(b.array[:,i,:]) + plt.show() +#%% + +plt.imshow(out2.subset( vertical=0).array) +plt.show() + +# Create least squares object instance with projector and data. +f = Norm2sq(Cop,b,c=0.5) + +# Initial guess +x_init = ImageData(geometry=vg, dimension_labels=['horizontal_x','horizontal_y','vertical']) +#invL = 0.5 +#g = f.grad(x_init) +#print (g) +#u = x_init - invL*f.grad(x_init) + +#%% +# Run FISTA for least squares without regularization +opt = {'tol': 1e-4, 'iter': 100} +x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt=opt) + +plt.imshow(x_fista0.subset(vertical=0).array) +plt.title('FISTA0') +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=opt) + +plt.imshow(x_fista0.subset(vertical=0).array) +plt.title('FISTA1') +plt.show() + +plt.semilogy(criter1) +plt.show() + +# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm +x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt) + +plt.imshow(x_fbpd1.subset(vertical=0).array) +plt.title('FBPD1') +plt.show() + +plt.semilogy(criter_fbpd1) +plt.show() + +# Now FBPD for least squares plus TV +#lamtv = 1 +#gtv = TV2D(lamtv) + +#x_fbpdtv, it_fbpdtv, timing_fbpdtv, criter_fbpdtv = FBPD(x_init,None,f,gtv,opt=opt) + +#plt.imshow(x_fbpdtv.subset(vertical=0).array) +#plt.show() + +#plt.semilogy(criter_fbpdtv) +#plt.show() + + +# Run CGLS, which should agree with the FISTA0 +x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Cop, b, opt=opt) + +plt.imshow(x_CGLS.subset(vertical=0).array) +plt.title('CGLS') +plt.title('CGLS recon, compare FISTA0') +plt.show() + +plt.semilogy(criter_CGLS) +plt.title('CGLS criterion') +plt.show() + + +#%% + +clims = (0,1) +cols = 3 +rows = 2 +current = 1 +fig = plt.figure() +# projections row +a=fig.add_subplot(rows,cols,current) +a.set_title('phantom {0}'.format(np.shape(Phantom.as_array()))) + +imgplot = plt.imshow(Phantom.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1]) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('FISTA0') +imgplot = plt.imshow(x_fista0.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1]) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('FISTA1') +imgplot = plt.imshow(x_fista1.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1]) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('FBPD1') +imgplot = plt.imshow(x_fbpd1.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1]) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('CGLS') +imgplot = plt.imshow(x_CGLS.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1]) + +plt.show() +#%% +#current = current + 1 +#a=fig.add_subplot(rows,cols,current) +#a.set_title('FBPD TV') +#imgplot = plt.imshow(x_fbpdtv.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1]) + +fig = plt.figure() +# projections row +b=fig.add_subplot(1,1,1) +b.set_title('criteria') +imgplot = plt.loglog(criter0 , label='FISTA0') +imgplot = plt.loglog(criter1 , label='FISTA1') +imgplot = plt.loglog(criter_fbpd1, label='FBPD1') +imgplot = plt.loglog(criter_CGLS, label='CGLS') +#imgplot = plt.loglog(criter_fbpdtv, label='FBPD TV') +b.legend(loc='right') +plt.show() +#%%
\ No newline at end of file |