From dac9d2ce209487765f66332547ce59f92a6c1d78 Mon Sep 17 00:00:00 2001
From: Jakob Jorgensen <jakob.jorgensen@manchester.ac.uk>
Date: Tue, 24 Apr 2018 11:50:53 +0100
Subject: Partial cleanup demo
---
Wrappers/Python/wip/simple_demo_ccpi.py | 98 +++++++++++++++------------------
1 file changed, 45 insertions(+), 53 deletions(-)
diff --git a/Wrappers/Python/wip/simple_demo_ccpi.py b/Wrappers/Python/wip/simple_demo_ccpi.py
index 3fdc2d4..6344c06 100755
--- a/Wrappers/Python/wip/simple_demo_ccpi.py
+++ b/Wrappers/Python/wip/simple_demo_ccpi.py
@@ -1,48 +1,44 @@
-#import sys
-#sys.path.append("..")
-from ccpi.framework import ImageData , AcquisitionData, ImageGeometry, AcquisitionGeometry
+# 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, AcquisitionData, 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, 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 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.
-# Set up phantom
+# Image parameters
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,
+# Set up image geometry
+ig = 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
-
+# 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:
@@ -55,6 +51,7 @@ while i < vert:
Phantom.fill(x, vertical=i)
i += 1
+# Display slice of phantom
if vert > 1:
plt.imshow(Phantom.subset(vertical=0).as_array())
else:
@@ -62,33 +59,28 @@ else:
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; # angles number
+det_w = 1.0
+det_num = N
-# 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)
+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
+
+ag = AcquisitionGeometry('parallel',
+ '3D',
+ angles,
+ N,
+ det_w,
+ vert,
+ det_w)
+
+# CCPi operator using image and acquisition geometries
+Cop = CCPiProjectorSimple(ig, ag)
# Try forward and backprojection
b = Cop.direct(Phantom)
--
cgit v1.2.3