Tidy up demo, remove prints and old code.

3 files changed, 94 insertions, 114 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/simple_demo_ccpi.py b/Wrappers/Python/wip/simple_demo_ccpi.py
index 6344c06..a8265ce 100755
--- a/Wrappers/Python/wip/simple_demo_ccpi.py
+++ b/Wrappers/Python/wip/simple_demo_ccpi.py
@@ -5,16 +5,12 @@
 # squares and 1-norm regularisation and FBPD for 1-norm regularisation.
 
 # First make all imports
-from ccpi.framework import ImageData, AcquisitionData, ImageGeometry, \
-                           AcquisitionGeometry
+from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry
 
 from ccpi.optimisation.algs import FISTA, FBPD, CGLS
-from ccpi.optimisation.funcs import Norm2sq, Norm1, TV2D
+from ccpi.optimisation.funcs import Norm2sq, Norm1
 
 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
@@ -63,14 +59,16 @@ plt.show()
 # 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; # angles number
+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
-
-#center_of_rotation = Phantom.get_dimension_size('horizontal_x') / 2
+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,
@@ -79,142 +77,146 @@ ag = AcquisitionGeometry('parallel',
                          vert,
                          det_w)
 
-# CCPi operator using image and acquisition geometries
+# Set up Operator object combining the ImageGeometry and AcquisitionGeometry
+# wrapping calls to CCPi projector.
 Cop = CCPiProjectorSimple(ig, ag)
 
-# Try forward and backprojection
+# 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)
-out2 = Cop.adjoint(b)
+z = 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.imshow(z.subset(vertical=0).array)
+plt.title('Backprojected data')
 plt.show()
 
-# Create least squares object instance with projector and data.
+# 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)
 
-# 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.title('FISTA Least squares')
 plt.show()
 
-# Now least squares plus 1-norm regularization
+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=opt)
+x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt)
 
-plt.imshow(x_fista0.subset(vertical=0).array)
-plt.title('FISTA1')
+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()
 
-# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm
+# 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('FBPD1')
+plt.title('FBPD for least squares plus 1-norm regularisation')
 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.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()
-# 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])
+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('FISTA0')
-imgplot = plt.imshow(x_fista0.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1])
+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('FISTA1')
-imgplot = plt.imshow(x_fista1.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1])
+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('FBPD1')
-imgplot = plt.imshow(x_fbpd1.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1])
+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('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])
+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()
-# 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
+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