summaryrefslogtreecommitdiffstats
path: root/Wrappers
diff options
context:
space:
mode:
Diffstat (limited to 'Wrappers')
-rwxr-xr-xWrappers/Python/ccpi/plugins/ops.py1
-rwxr-xr-xWrappers/Python/ccpi/plugins/processors.py31
-rwxr-xr-xWrappers/Python/wip/demo_ccpi_simple.py222
-rwxr-xr-xWrappers/Python/wip/simple_demo_ccpi.py228
4 files changed, 227 insertions, 255 deletions
diff --git a/Wrappers/Python/ccpi/plugins/ops.py b/Wrappers/Python/ccpi/plugins/ops.py
index aeb51af..75c5db9 100755
--- a/Wrappers/Python/ccpi/plugins/ops.py
+++ b/Wrappers/Python/ccpi/plugins/ops.py
@@ -109,6 +109,5 @@ class CCPiProjectorSimple(Operator):
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 \ No newline at end of file
diff --git a/Wrappers/Python/ccpi/plugins/processors.py b/Wrappers/Python/ccpi/plugins/processors.py
index df580e0..9938b9e 100755
--- a/Wrappers/Python/ccpi/plugins/processors.py
+++ b/Wrappers/Python/ccpi/plugins/processors.py
@@ -17,15 +17,10 @@
# See the License for the specific language governing permissions and
# limitations under the License
-from ccpi.framework import DataProcessor, DataContainer, AcquisitionData,\
- AcquisitionGeometry, ImageGeometry, ImageData
+from ccpi.framework import DataProcessor, 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
-
def setupCCPiGeometries(voxel_num_x, voxel_num_y, voxel_num_z, angles, counter):
@@ -53,10 +48,9 @@ def setupCCPiGeometries(voxel_num_x, voxel_num_y, voxel_num_z, angles, counter):
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))
@@ -65,10 +59,8 @@ def setupCCPiGeometries(voxel_num_x, voxel_num_y, voxel_num_z, angles, counter):
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
@@ -119,18 +111,6 @@ class CCPiForwardProjector(DataProcessor):
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,
@@ -179,7 +159,6 @@ class CCPiBackwardProjector(DataProcessor):
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:
@@ -198,7 +177,7 @@ class CCPiBackwardProjector(DataProcessor):
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(),
diff --git a/Wrappers/Python/wip/demo_ccpi_simple.py b/Wrappers/Python/wip/demo_ccpi_simple.py
new file mode 100755
index 0000000..a8265ce
--- /dev/null
+++ b/Wrappers/Python/wip/demo_ccpi_simple.py
@@ -0,0 +1,222 @@
+
+# This demo illustrates how CCPi 2D parallel-beam projectors can be used with
+# the modular optimisation framework. The demo sets up a small 4-slice 3D test
+# case and demonstrates reconstruction using CGLS, as well as FISTA for least
+# squares and 1-norm regularisation and FBPD for 1-norm regularisation.
+
+# First make all imports
+from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry
+
+from ccpi.optimisation.algs import FISTA, FBPD, CGLS
+from ccpi.optimisation.funcs import Norm2sq, Norm1
+
+from ccpi.plugins.ops import CCPiProjectorSimple
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Set up phantom size N x N x vert by creating ImageGeometry, initialising the
+# ImageData object with this geometry and empty array and finally put some
+# data into its array, and display one slice as image.
+
+# Image parameters
+N = 128
+vert = 4
+
+# Set up image geometry
+ig = ImageGeometry(voxel_num_x=N,
+ voxel_num_y=N,
+ voxel_num_z=vert)
+
+# Set up empty image data
+Phantom = ImageData(geometry=ig,
+ dimension_labels=['horizontal_x',
+ 'horizontal_y',
+ 'vertical'])
+
+# Populate image data by looping over and filling slices
+i = 0
+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
+ if vert > 1 :
+ Phantom.fill(x, vertical=i)
+ i += 1
+
+# Display slice of phantom
+if vert > 1:
+ plt.imshow(Phantom.subset(vertical=0).as_array())
+else:
+ plt.imshow(Phantom.as_array())
+plt.show()
+
+
+# Set up AcquisitionGeometry object to hold the parameters of the measurement
+# setup geometry: # Number of angles, the actual angles from 0 to
+# pi for parallel beam, set the width of a detector
+# pixel relative to an object pixe and the number of detector pixels.
+angles_num = 20
+det_w = 1.0
+det_num = N
+
+angles = np.linspace(0,np.pi,angles_num,endpoint=False,dtype=np.float32)*\
+ 180/np.pi
+
+# Inputs: Geometry, 2D or 3D, angles, horz detector pixel count,
+# horz detector pixel size, vert detector pixel count,
+# vert detector pixel size.
+ag = AcquisitionGeometry('parallel',
+ '3D',
+ angles,
+ N,
+ det_w,
+ vert,
+ det_w)
+
+# Set up Operator object combining the ImageGeometry and AcquisitionGeometry
+# wrapping calls to CCPi projector.
+Cop = CCPiProjectorSimple(ig, ag)
+
+# Forward and backprojection are available as methods direct and adjoint. Here
+# generate test data b and do simple backprojection to obtain z. Display all
+# data slices as images, and a single backprojected slice.
+b = Cop.direct(Phantom)
+z = Cop.adjoint(b)
+
+for i in range(b.get_dimension_size('vertical')):
+ plt.imshow(b.subset(vertical=i).array)
+ plt.show()
+
+plt.imshow(z.subset(vertical=0).array)
+plt.title('Backprojected data')
+plt.show()
+
+# Using the test data b, different reconstruction methods can now be set up as
+# demonstrated in the rest of this file. In general all methods need an initial
+# guess and some algorithm options to be set. Note that 100 iterations for
+# some of the methods is a very low number and 1000 or 10000 iterations may be
+# needed if one wants to obtain a converged solution.
+x_init = ImageData(geometry=ig,
+ dimension_labels=['horizontal_x','horizontal_y','vertical'])
+opt = {'tol': 1e-4, 'iter': 100}
+
+# First a CGLS reconstruction can be done:
+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.show()
+
+plt.semilogy(criter_CGLS)
+plt.title('CGLS criterion')
+plt.show()
+
+# CGLS solves the simple least-squares problem. The same problem can be solved
+# by FISTA by setting up explicitly a least squares function object and using
+# no regularisation:
+
+# Create least squares object instance with projector, test data and a constant
+# coefficient of 0.5:
+f = Norm2sq(Cop,b,c=0.5)
+
+# Run FISTA for least squares without regularization
+x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt=opt)
+
+plt.imshow(x_fista0.subset(vertical=0).array)
+plt.title('FISTA Least squares')
+plt.show()
+
+plt.semilogy(criter0)
+plt.title('FISTA Least squares criterion')
+plt.show()
+
+# FISTA can also solve regularised forms by specifying a second function object
+# such as 1-norm regularisation with choice of regularisation parameter lam:
+
+# Create 1-norm function object
+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)
+
+plt.imshow(x_fista1.subset(vertical=0).array)
+plt.title('FISTA Least squares plus 1-norm regularisation')
+plt.show()
+
+plt.semilogy(criter1)
+plt.title('FISTA Least squares plus 1-norm regularisation criterion')
+plt.show()
+
+# The least squares plus 1-norm regularisation problem can also be solved by
+# other algorithms such as the Forward Backward Primal Dual algorithm. This
+# algorithm minimises the sum of three functions and the least squares and
+# 1-norm functions should be given as the second and third function inputs.
+# In this test case, this algorithm requires more iterations to converge, so
+# new options are specified.
+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('FBPD for least squares plus 1-norm regularisation')
+plt.show()
+
+plt.semilogy(criter_fbpd1)
+plt.title('FBPD for least squares plus 1-norm regularisation criterion')
+plt.show()
+
+
+# Compare all reconstruction and criteria
+
+clims = (0,1)
+cols = 3
+rows = 2
+current = 1
+
+fig = plt.figure()
+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])
+plt.axis('off')
+
+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.axis('off')
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('FISTA LS')
+imgplot = plt.imshow(x_fista0.subset(vertical=0).as_array(),
+ vmin=clims[0],vmax=clims[1])
+plt.axis('off')
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('FISTA LS+1')
+imgplot = plt.imshow(x_fista1.subset(vertical=0).as_array(),
+ vmin=clims[0],vmax=clims[1])
+plt.axis('off')
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('FBPD LS+1')
+imgplot = plt.imshow(x_fbpd1.subset(vertical=0).as_array(),
+ vmin=clims[0],vmax=clims[1])
+plt.axis('off')
+
+fig = plt.figure()
+b=fig.add_subplot(1,1,1)
+b.set_title('criteria')
+imgplot = plt.loglog(criter_CGLS, label='CGLS')
+imgplot = plt.loglog(criter0 , label='FISTA LS')
+imgplot = plt.loglog(criter1 , label='FISTA LS+1')
+imgplot = plt.loglog(criter_fbpd1, label='FBPD LS+1')
+b.legend(loc='lower left')
+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
deleted file mode 100755
index 3fdc2d4..0000000
--- a/Wrappers/Python/wip/simple_demo_ccpi.py
+++ /dev/null
@@ -1,228 +0,0 @@
-#import sys
-#sys.path.append("..")
-
-from ccpi.framework import ImageData , AcquisitionData, ImageGeometry, AcquisitionGeometry
-
-from ccpi.optimisation.algs import FISTA, FBPD, CGLS
-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
-import matplotlib.pyplot as plt
-
-test_case = 1 # 1=parallel2D, 2=cone2D, 3=parallel3D
-
-# Set up phantom
-N = 128
-vert = 4
-# 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
-
-vg = ImageGeometry(voxel_num_x=N,
- voxel_num_y=N,
- voxel_num_z=vert)
-
-Phantom = ImageData(geometry=vg,dimension_labels=['horizontal_x','horizontal_y','vertical']) + 0.1
-
-i = 0
-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
- if vert > 1 :
- Phantom.fill(x, vertical=i)
- i += 1
-
-if vert > 1:
- plt.imshow(Phantom.subset(vertical=0).as_array())
-else:
- plt.imshow(Phantom.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,
- 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')
- 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
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-05 10:37:19 +0000