diff options
| -rw-r--r-- | Wrappers/Python/ccpi/framework.py | 63 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/reconstruction/astra_ops.py | 10 | ||||
| -rw-r--r-- | Wrappers/Python/test/simple_demo.py | 20 |
3 files changed, 55 insertions, 38 deletions
diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py index 2b0ba76..5ac17ee 100644 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework.py @@ -105,7 +105,8 @@ class DataSet(object): self.shape = numpy.shape(array) self.number_of_dimensions = len (self.shape) self.dimension_labels = {} - + self.geometry = None # Only relevant for SinogramData and VolumeData + if dimension_labels is not None and \ len (dimension_labels) == self.number_of_dimensions: for i in range(self.number_of_dimensions): @@ -123,6 +124,13 @@ class DataSet(object): raise TypeError('Array must be NumpyArray, passed {0}'\ .format(type(array))) + # finally copy the geometry + if 'geometry' in kwargs.keys(): + self.geometry = kwargs['geometry'] + else: + # assume it is parallel beam + pass + def as_array(self, dimensions=None): '''Returns the DataSet as Numpy Array @@ -212,7 +220,8 @@ class DataSet(object): out = self.as_array() + other.as_array() return type(self)(out, deep_copy=True, - dimension_labels=self.dimension_labels) + dimension_labels=self.dimension_labels, + geometry=self.geometry) else: raise ValueError('Wrong shape: {0} and {1}'.format(self.shape, other.shape)) @@ -220,7 +229,8 @@ class DataSet(object): return type(self)( self.as_array() + other, deep_copy=True, - dimension_labels=self.dimension_labels) + dimension_labels=self.dimension_labels, + geometry=self.geometry) else: raise TypeError('Cannot {0} DataSet with {1}'.format("add" , type(other))) @@ -232,14 +242,16 @@ class DataSet(object): out = self.as_array() - other.as_array() return type(self)(out, deep_copy=True, - dimension_labels=self.dimension_labels) + dimension_labels=self.dimension_labels, + geometry=self.geometry) else: 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, deep_copy=True, - dimension_labels=self.dimension_labels) + dimension_labels=self.dimension_labels, + geometry=self.geometry) else: raise TypeError('Cannot {0} DataSet with {1}'.format("subtract" , type(other))) @@ -254,14 +266,16 @@ class DataSet(object): out = self.as_array() / other.as_array() return type(self)(out, deep_copy=True, - dimension_labels=self.dimension_labels) + dimension_labels=self.dimension_labels, + geometry=self.geometry) else: 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, deep_copy=True, - dimension_labels=self.dimension_labels) + dimension_labels=self.dimension_labels, + geometry=self.geometry) else: raise TypeError('Cannot {0} DataSet with {1}'.format("divide" , type(other))) @@ -273,14 +287,16 @@ class DataSet(object): out = self.as_array() ** other.as_array() return type(self)(out, deep_copy=True, - dimension_labels=self.dimension_labels) + dimension_labels=self.dimension_labels, + geometry=self.geometry) else: 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, deep_copy=True, - dimension_labels=self.dimension_labels) + dimension_labels=self.dimension_labels, + geometry=self.geometry) else: raise TypeError('Cannot {0} DataSet with {1}'.format("power" , type(other))) @@ -292,14 +308,16 @@ class DataSet(object): out = self.as_array() * other.as_array() return type(self)(out, deep_copy=True, - dimension_labels=self.dimension_labels) + dimension_labels=self.dimension_labels, + geometry=self.geometry) else: 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, deep_copy=True, - dimension_labels=self.dimension_labels) + dimension_labels=self.dimension_labels, + geometry=self.geometry) else: raise TypeError('Cannot {0} DataSet with {1}'.format("multiply" , type(other))) @@ -315,19 +333,22 @@ class DataSet(object): out = numpy.abs(self.as_array() ) return type(self)(out, deep_copy=True, - dimension_labels=self.dimension_labels) + dimension_labels=self.dimension_labels, + geometry=self.geometry) def maximum(self,otherscalar): out = numpy.maximum(self.as_array(),otherscalar) return type(self)(out, deep_copy=True, - dimension_labels=self.dimension_labels) + dimension_labels=self.dimension_labels, + geometry=self.geometry) def sign(self): out = numpy.sign(self.as_array() ) return type(self)(out, deep_copy=True, - dimension_labels=self.dimension_labels) + dimension_labels=self.dimension_labels, + geometry=self.geometry) # reverse operand def __radd__(self, other): @@ -353,11 +374,13 @@ class DataSet(object): if isinstance(other, (int, float)) : fother = numpy.ones(numpy.shape(self.array)) * other return type(self)(fother ** self.array , - dimension_labels=self.dimension_labels) + dimension_labels=self.dimension_labels, + geometry=self.geometry) elif issubclass(other, DataSet): if self.checkDimensions(other): return type(self)(other.as_array() ** self.array , - dimension_labels=self.dimension_labels) + dimension_labels=self.dimension_labels, + geometry=self.geometry) else: raise ValueError('Dimensions do not match') # __rpow__ @@ -473,13 +496,7 @@ class SinogramData(DataSet): dimension_labels = ['angle' , 'horizontal'] DataSet.__init__(self, array, deep_copy, dimension_labels, **kwargs) - - # finally copy the instrument geometry - if 'instrument_geometry' in kwargs.keys(): - self.instrument_geometry = kwargs['instrument_geometry'] - else: - # assume it is parallel beam - pass + class DataSetProcessor(object): '''Defines a generic DataSet processor diff --git a/Wrappers/Python/ccpi/reconstruction/astra_ops.py b/Wrappers/Python/ccpi/reconstruction/astra_ops.py index bc6f5e1..6d527ba 100644 --- a/Wrappers/Python/ccpi/reconstruction/astra_ops.py +++ b/Wrappers/Python/ccpi/reconstruction/astra_ops.py @@ -26,6 +26,12 @@ class AstraProjectorSimple(Operator): def __init__(self, geomv, geomp, device): super(AstraProjectorSimple, self).__init__() + # Store our volume and sinogram geometries. Redundant with also + # storing in ASTRA format below but needed to assign to + # SinogramData in "direct" method and VolumeData in "adjoint" method + self.sinogram_geometry = geomp + self.volume_geometry = geomv + # ASTRA Volume geometry self.vol_geom = astra.create_vol_geom(geomv.voxel_num_x, \ geomv.voxel_num_y, \ @@ -63,12 +69,12 @@ class AstraProjectorSimple(Operator): sinogram_id, DATA = astra.create_sino(IM.as_array(), self.proj_id) astra.data2d.delete(sinogram_id) - return SinogramData(DATA) + return SinogramData(DATA,geometry=self.sinogram_geometry) def adjoint(self, DATA): rec_id, IM = astra.create_backprojection(DATA.as_array(), self.proj_id) astra.data2d.delete(rec_id) - return VolumeData(IM) + return VolumeData(IM,geometry=self.volume_geometry) def delete(self): astra.data2d.delete(self.proj_id) diff --git a/Wrappers/Python/test/simple_demo.py b/Wrappers/Python/test/simple_demo.py index 7a6b049..7a28ae9 100644 --- a/Wrappers/Python/test/simple_demo.py +++ b/Wrappers/Python/test/simple_demo.py @@ -25,32 +25,26 @@ plt.show() vg = VolumeGeometry(N,N,None, 1,1,None) -Phantom = VolumeData(x) -#Phantom = VolumeData(x) +Phantom = VolumeData(x,geometry=vg) # Set up measurement geometry angles_num = 20; # angles number angles = np.linspace(0,np.pi,angles_num,endpoint=False) - det_w = 1.0 det_num = N - SourceOrig = 500 OrigDetec = 0 # Parallelbeam geometry test pg = SinogramGeometry('parallel','2D',angles,det_num,det_w) -# Set up ASTRA projector -#Aop = AstraProjector(vg, angles, N,'gpu') -#Aop = AstraProjector(det_w, det_num, angles, projtype='parallel',device='gpu') - -Aop_old = AstraProjector(det_w, det_num, SourceOrig, - OrigDetec, angles, - N,'fanbeam','gpu') - +# ASTRA operator using volume and sinogram geometries Aop = AstraProjectorSimple(vg, pg, 'gpu') +# 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) @@ -66,7 +60,7 @@ plt.show() f = Norm2sq(Aop,b,c=0.5) # Initial guess -x_init = VolumeData(np.zeros(x.shape)) +x_init = VolumeData(np.zeros(x.shape),geometry=vg) # Run FISTA for least squares without regularization x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None) |