From 1f9e7441198a6f088dfed311a067e8acf05d4d3b Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Tue, 17 Apr 2018 18:06:41 +0100
Subject: removed LinearOperatorMatrix and added check of geometries in ccpi
operator
---
Wrappers/Python/ccpi/plugins/ops.py | 72 +++++++++++++++++++++----------------
1 file changed, 42 insertions(+), 30 deletions(-)
(limited to 'Wrappers/Python')
diff --git a/Wrappers/Python/ccpi/plugins/ops.py b/Wrappers/Python/ccpi/plugins/ops.py
index 0028cf1..aeb51af 100755
--- a/Wrappers/Python/ccpi/plugins/ops.py
+++ b/Wrappers/Python/ccpi/plugins/ops.py
@@ -19,49 +19,61 @@
import numpy
from ccpi.optimisation.ops import Operator, PowerMethodNonsquare
-from ccpi.framework import ImageData, DataContainer
-from ccpi.plugins.processors import CCPiBackwardProjector, CCPiForwardProjector
-
-class LinearOperatorMatrix(Operator):
- def __init__(self,A):
- self.A = A
- self.s1 = None # Largest singular value, initially unknown
- super(LinearOperatorMatrix, self).__init__()
-
- def direct(self,x):
- return DataContainer(numpy.dot(self.A,x.as_array()))
-
- def adjoint(self,x):
- return DataContainer(numpy.dot(self.A.transpose(),x.as_array()))
-
- def size(self):
- return self.A.shape
-
- def get_max_sing_val(self):
- # If unknown, compute and store. If known, simply return it.
- if self.s1 is None:
- self.s1 = svds(self.A,1,return_singular_vectors=False)[0]
- return self.s1
- else:
- return self.s1
+from ccpi.framework import ImageData, DataContainer , \
+ ImageGeometry, AcquisitionGeometry
+from ccpi.plugins.processors import CCPiBackwardProjector, \
+ CCPiForwardProjector , setupCCPiGeometries
class CCPiProjectorSimple(Operator):
"""ASTRA projector modified to use DataSet and geometry."""
- def __init__(self, geomv, geomp):
+ def __init__(self, geomv, geomp, default=False):
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,
+ if geomp.geom_type == "cone":
+ raise TypeError('Can only handle parallel beam')
+
+ # set-up the geometries if compatible
+ geoms = setupCCPiGeometries(geomv.voxel_num_x,geomv.voxel_num_y,
+ geomv.voxel_num_z, geomp.angles, 0)
+
+
+ vg = ImageGeometry(voxel_num_x=geoms['output_volume_x'],
+ voxel_num_y=geoms['output_volume_y'],
+ voxel_num_z=geoms['output_volume_z'])
+
+ pg = AcquisitionGeometry('parallel',
+ '3D',
+ geomp.angles,
+ geoms['n_h'], geomp.pixel_size_h,
+ geoms['n_v'], geomp.pixel_size_v #2D in 3D is a slice 1 pixel thick
+ )
+ if not default:
+ # check if geometry is the same (on the voxels)
+ if not ( vg.voxel_num_x == geomv.voxel_num_x and \
+ vg.voxel_num_y == geomv.voxel_num_y and \
+ vg.voxel_num_z == geomv.voxel_num_z ):
+ msg = 'The required volume geometry will not work\nThe following would\n'
+ msg += vg.__str__()
+ raise ValueError(msg)
+ if not (pg.pixel_num_h == geomp.pixel_num_h and \
+ pg.pixel_num_v == geomp.pixel_num_v and \
+ len( pg.angles ) == len( geomp.angles ) ) :
+ msg = 'The required acquisition geometry will not work\nThe following would\n'
+ msg += pg.__str__()
+ raise ValueError(msg)
+
+ self.fp = CCPiForwardProjector(image_geometry=vg,
+ acquisition_geometry=pg,
output_axes_order=['angle','vertical','horizontal'])
- self.bp = CCPiBackwardProjector(image_geometry=geomv,
- acquisition_geometry=geomp,
+ self.bp = CCPiBackwardProjector(image_geometry=vg,
+ acquisition_geometry=pg,
output_axes_order=['horizontal_x','horizontal_y','vertical'])
# Initialise empty for singular value.
--
cgit v1.2.3
From 1d7313fe6f66d4b87efab95c88cd510ff44e26e1 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Tue, 17 Apr 2018 18:07:25 +0100
Subject: added setupCCPiGeometries
---
Wrappers/Python/ccpi/plugins/processors.py | 45 ++++++++++++++++++++++++++++++
1 file changed, 45 insertions(+)
(limited to 'Wrappers/Python')
diff --git a/Wrappers/Python/ccpi/plugins/processors.py b/Wrappers/Python/ccpi/plugins/processors.py
index e05c6ca..df580e0 100755
--- a/Wrappers/Python/ccpi/plugins/processors.py
+++ b/Wrappers/Python/ccpi/plugins/processors.py
@@ -27,6 +27,51 @@ from scipy import ndimage
import matplotlib.pyplot as plt
+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'], 1.0,
+ geoms['n_v'], 1.0 #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
+
+
class CCPiForwardProjector(DataProcessor):
'''Normalization based on flat and dark
--
cgit v1.2.3
From 0fc626543be8c0827ca1696ecae5f57ff05b0fbc Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Tue, 17 Apr 2018 18:08:22 +0100
Subject: started cleaning demo
---
Wrappers/Python/wip/simple_demo_ccpi.py | 78 +++++++--------------------------
1 file changed, 17 insertions(+), 61 deletions(-)
(limited to 'Wrappers/Python')
diff --git a/Wrappers/Python/wip/simple_demo_ccpi.py b/Wrappers/Python/wip/simple_demo_ccpi.py
index 10bb75f..3fdc2d4 100755
--- a/Wrappers/Python/wip/simple_demo_ccpi.py
+++ b/Wrappers/Python/wip/simple_demo_ccpi.py
@@ -8,7 +8,6 @@ from ccpi.optimisation.funcs import Norm2sq, Norm1 , TV2D
from ccpi.plugins.ops import CCPiProjectorSimple
from ccpi.plugins.processors import CCPiForwardProjector, CCPiBackwardProjector
-
from ccpi.reconstruction.parallelbeam import alg as pbalg
import numpy as np
@@ -18,7 +17,7 @@ test_case = 1 # 1=parallel2D, 2=cone2D, 3=parallel3D
# Set up phantom
N = 128
-vert = 1
+vert = 4
# Set up measurement geometry
angles_num = 20; # angles number
det_w = 1.0
@@ -38,71 +37,28 @@ elif test_case == 3:
else:
NotImplemented
+vg = ImageGeometry(voxel_num_x=N,
+ voxel_num_y=N,
+ voxel_num_z=vert)
-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
+while i < vert:
+ if vert > 1:
+ x = Phantom.subset(vertical=i).array
+ else:
+ x = Phantom.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)
+ if vert > 1 :
+ Phantom.fill(x, vertical=i)
i += 1
-plt.imshow(Phantom.subset(vertical=0).as_array())
+if vert > 1:
+ plt.imshow(Phantom.subset(vertical=0).as_array())
+else:
+ plt.imshow(Phantom.as_array())
plt.show()
@@ -116,8 +72,8 @@ if test_case==1:
pg = AcquisitionGeometry('parallel',
'3D',
angles,
- geoms['n_h'],det_w,
- geoms['n_v'], det_w #2D in 3D is a slice 1 pixel thick
+ N , det_w,
+ vert , det_w #2D in 3D is a slice 1 pixel thick
)
elif test_case==2:
raise NotImplemented('cone beam projector not yet available')
--
cgit v1.2.3