From 737a7aa27e0f55718cbcef4909ff9d64457a7ef1 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 21 Mar 2018 13:17:07 +0000 Subject: Ccpi forward backprojector (#77) * added simple test for SinogramData * initial version of processors This contains a version of a simple normalization algorithm. * added default pixel size = 1 (mm) * Added FindCenter of rotation for parallel beam geometry closes #5 * WIP for ccpi forward/backprojectors * WIP * Added Forward and Backprojectors from ccpi added processor to pad acquisition data to center on center of rotation. requires update of ccpi-reconstruction. * Added simple_demo_ccpi.py Currently there is a problem with the order of the axis after a projection/backprojection with the ccpi operator. In the calculation of the PowerMethod there is a wrong order of the axis. Possibly the creation of the data should be handled by the operator. * fixed calculation of power method * further work * fixed the x_init for the FISTA example * fixed ccpi FISTA demo * Fix missing subset call for display * Added create_image_data method to AstraOperatorSimple * Stack direction 0, not 2, in TV * Added finite different forward processor * wip * add io to setup.py * wip * comment out FiniteDifferentiator * uniformed algorithm interface Added some docstring. Added parsing of opt. * removed ccpi.io * little changes * fix setupCCPiGeometries * added simple astra demo * use minimal size for example --- Wrappers/Python/ccpi/astra/astra_ops.py | 6 +- Wrappers/Python/ccpi/framework.py | 278 +++++++++++++++++----- Wrappers/Python/ccpi/processors.py | 357 +++++++++++++++++++++++++++- Wrappers/Python/ccpi/reconstruction/algs.py | 89 +++++-- Wrappers/Python/ccpi/reconstruction/ops.py | 85 ++++++- Wrappers/Python/wip/simple_demo_astra.py | 189 +++++++++++++++ Wrappers/Python/wip/simple_demo_ccpi.py | 270 +++++++++++++++++++++ 7 files changed, 1183 insertions(+), 91 deletions(-) create mode 100755 Wrappers/Python/wip/simple_demo_astra.py create mode 100755 Wrappers/Python/wip/simple_demo_ccpi.py (limited to 'Wrappers/Python') 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(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 -- cgit v1.2.3