summaryrefslogtreecommitdiffstats
path: root/Wrappers
diff options
context:
space:
mode:
authorjakobsj <jakobsj@users.noreply.github.com>2018-05-11 13:13:43 +0100
committerEdoardo Pasca <edo.paskino@gmail.com>2018-05-11 14:13:43 +0200
commit992146ad44767f9f34515393b608ec2ca0304cd1 (patch)
treee1f2b14746bab724d76e39c0e14b51ee3bd87c51 /Wrappers
parent00b6766211f6fa77a599d6925b6674b0339170fd (diff)
downloadframework-plugins-992146ad44767f9f34515393b608ec2ca0304cd1.tar.gz
framework-plugins-992146ad44767f9f34515393b608ec2ca0304cd1.tar.bz2
framework-plugins-992146ad44767f9f34515393b608ec2ca0304cd1.tar.xz
framework-plugins-992146ad44767f9f34515393b608ec2ca0304cd1.zip
Simplify to TV only and comment RGLTK recon demo (#10)
Diffstat (limited to 'Wrappers')
-rw-r--r--Wrappers/Python/wip/demo_simple_RGLTK.py247
1 files changed, 105 insertions, 142 deletions
diff --git a/Wrappers/Python/wip/demo_simple_RGLTK.py b/Wrappers/Python/wip/demo_simple_RGLTK.py
index 9f0a4c3..3831603 100644
--- a/Wrappers/Python/wip/demo_simple_RGLTK.py
+++ b/Wrappers/Python/wip/demo_simple_RGLTK.py
@@ -1,214 +1,177 @@
+# This demo illustrates how the CCPi Regularisation Toolkit can be used
+# as TV regularisation for use with the FISTA algorithm of the modular
+# optimisation framework and compares with the FBPD TV implementation.
+
+# All own imports
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.astra.ops import AstraProjectorSimple
from ccpi.plugins.regularisers import _ROF_TV_, _FGP_TV_
+# All external imports
import numpy as np
import matplotlib.pyplot as plt
-test_case = 1 # 1=parallel2D, 2=cone2D
+# Choose either a parallel-beam (1=parallel2D) or fan-beam (2=cone2D) test case
+test_case = 1
-# Set up phantom
+# Set up phantom size NxN by creating ImageGeometry, initialising the
+# ImageData object with this geometry and empty array and finally put some
+# data into its array, and display as image.
N = 128
-
-
-vg = ImageGeometry(voxel_num_x=N,voxel_num_y=N)
-Phantom = ImageData(geometry=vg)
+ig = ImageGeometry(voxel_num_x=N,voxel_num_y=N)
+Phantom = ImageData(geometry=ig)
x = Phantom.as_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)] = 1
-#plt.imshow(x)
-#plt.show()
+plt.imshow(x)
+plt.title('Phantom image')
+plt.show()
-# Set up measurement geometry
-angles_num = 20; # angles number
+# 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 and 0 to 2pi for fanbeam, set the width of a detector
+# pixel relative to an object pixel, the number of detector pixels, and the
+# source-origin and origin-detector distance (here the origin-detector distance
+# set to 0 to simulate a "virtual detector" with same detector pixel size as
+# object pixel size).
+angles_num = 20
+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)
+ ag = AcquisitionGeometry('parallel',
+ '2D',
+ angles,
+ det_num,det_w)
elif test_case==2:
angles = np.linspace(0,2*np.pi,angles_num,endpoint=False)
+ ag = AcquisitionGeometry('cone',
+ '2D',
+ angles,
+ det_num,
+ det_w,
+ dist_source_center=SourceOrig,
+ dist_center_detector=OrigDetec)
else:
NotImplemented
-det_w = 1.0
-det_num = N
-SourceOrig = 200
-OrigDetec = 0
+# Set up Operator object combining the ImageGeometry and AcquisitionGeometry
+# wrapping calls to ASTRA as well as specifying whether to use CPU or GPU.
+Aop = AstraProjectorSimple(ig, ag, 'cpu')
-# Parallelbeam geometry test
-if test_case==1:
- pg = AcquisitionGeometry('parallel',
- '2D',
- angles,
- det_num,det_w)
-elif test_case==2:
- pg = AcquisitionGeometry('cone',
- '2D',
- angles,
- det_num,
- det_w,
- dist_source_center=SourceOrig,
- dist_center_detector=OrigDetec)
-
-# ASTRA operator using volume and sinogram geometries
-Aop = AstraProjectorSimple(vg, pg, 'cpu')
-
-# Unused old astra projector without geometry
-# Aop_old = AstraProjector(det_w, det_num, SourceOrig,
-# OrigDetec, angles,
-# N,'fanbeam','gpu')
-
-# 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.
b = Aop.direct(Phantom)
-out2 = Aop.adjoint(b)
+z = Aop.adjoint(b)
-#plt.imshow(b.array)
-#plt.show()
+plt.imshow(b.array)
+plt.title('Simulated data')
+plt.show()
-#plt.imshow(out2.array)
-#plt.show()
+plt.imshow(z.array)
+plt.title('Backprojected data')
+plt.show()
# Create least squares object instance with projector and data.
f = Norm2sq(Aop,b,c=0.5)
# Initial guess
-x_init = ImageData(np.zeros(x.shape),geometry=vg)
-#%%
-# FISTA with ROF-TV regularisation
-g_rof = _ROF_TV_(lambdaReg = 10.0,iterationsTV=50,tolerance=1e-5,time_marchstep=0.01,device='cpu')
+x_init = ImageData(np.zeros(x.shape),geometry=ig)
-opt = {'tol': 1e-4, 'iter': 100}
-
-x_fista_rof, it1, timing1, criter_rof = FISTA(x_init, f, g_rof,opt)
+# Set up FBPD algorithm for TV reconstruction and solve
+opt_FBPD = {'tol': 1e-4, 'iter': 10000}
-plt.figure()
-plt.subplot(121)
-plt.imshow(x_fista_rof.array,cmap="BuPu")
-plt.title('FISTA-ROF-TV')
-plt.subplot(122)
-plt.semilogy(criter_rof)
-plt.show()
-#%%
-# FISTA with FGP-TV regularisation
-g_fgp = _FGP_TV_(lambdaReg = 10.0,iterationsTV=50,tolerance=1e-5,methodTV=0,nonnegativity=0,printing=0,device='cpu')
+lamtv = 1.0
+gtv = TV2D(lamtv)
-x_fista_fgp, it1, timing1, criter_fgp = FISTA(x_init, f, g_fgp,opt)
+x_fbpdtv, it_fbpdtv, timing_fbpdtv, criter_fbpdtv = FBPD(x_init,
+ None,
+ f,
+ gtv,
+ opt=opt_FBPD)
plt.figure()
plt.subplot(121)
-plt.imshow(x_fista_fgp.array,cmap="BuPu")
-plt.title('FISTA-FGP-TV')
+plt.imshow(x_fbpdtv.array)
+plt.title('FBPD TV')
plt.subplot(122)
-plt.semilogy(criter_fgp)
-plt.show()
-#%%
-# Run FISTA for least squares without regularization
-x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt)
-
-plt.imshow(x_fista0.array)
-plt.title('FISTA0')
+plt.semilogy(criter_fbpdtv)
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)
+# Set up the ROF variant of TV from the CCPi Regularisation Toolkit and run
+# TV-reconstruction using FISTA
+g_rof = _ROF_TV_(lambdaReg = lamtv,
+ iterationsTV=50,
+ tolerance=1e-5,
+ time_marchstep=0.01,
+ device='cpu')
-plt.imshow(x_fista1.array)
-plt.title('FISTA1')
-plt.show()
-
-plt.semilogy(criter1)
-plt.show()
-#%%
-# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm
opt = {'tol': 1e-4, 'iter': 100}
-x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt)
-
-plt.imshow(x_fbpd1.array)
-plt.title('FBPD1')
-plt.show()
-
-plt.semilogy(criter_fbpd1)
-plt.show()
-#%%
-opt_FBPD = {'tol': 1e-4, 'iter': 10000}
-# Now FBPD for least squares plus TV
-lamtv = 10.0
-gtv = TV2D(lamtv)
-x_fbpdtv, it_fbpdtv, timing_fbpdtv, criter_fbpdtv = FBPD(x_init,None,f,gtv,opt=opt_FBPD)
-
-plt.imshow(x_fbpdtv.array)
-plt.show()
+x_fista_rof, it1, timing1, criter_rof = FISTA(x_init, f, g_rof,opt)
-plt.semilogy(criter_fbpdtv)
+plt.figure()
+plt.subplot(121)
+plt.imshow(x_fista_rof.array)
+plt.title('FISTA ROF TV')
+plt.subplot(122)
+plt.semilogy(criter_rof)
plt.show()
+# Repeat for FGP variant.
+g_fgp = _FGP_TV_(lambdaReg = lamtv,
+ iterationsTV=50,
+ tolerance=1e-5,
+ methodTV=0,
+ nonnegativity=0,
+ printing=0,
+ device='cpu')
-# Run CGLS, which should agree with the FISTA0
-x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt )
-
-plt.imshow(x_CGLS.array)
-plt.title('CGLS')
-#plt.title('CGLS recon, compare FISTA0')
-plt.show()
+x_fista_fgp, it1, timing1, criter_fgp = FISTA(x_init, f, g_fgp,opt)
-plt.semilogy(criter_CGLS)
-plt.title('CGLS criterion')
+plt.figure()
+plt.subplot(121)
+plt.imshow(x_fista_fgp.array)
+plt.title('FISTA FGP TV')
+plt.subplot(122)
+plt.semilogy(criter_fgp)
plt.show()
-#%%
+# Compare all reconstruction and criteria
clims = (0,1)
cols = 3
-rows = 2
+rows = 1
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.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.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.as_array(),vmin=clims[0],vmax=clims[1])
+a.set_title('FBPD TV')
+imgplot = plt.imshow(x_fbpdtv.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.as_array(),vmin=clims[0],vmax=clims[1])
+a.set_title('FISTA ROF TV')
+imgplot = plt.imshow(x_fista_rof.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.as_array(),vmin=clims[0],vmax=clims[1])
-
-#current = current + 1
-#a=fig.add_subplot(rows,cols,current)
-#a.set_title('FBPD TV')
-#imgplot = plt.imshow(x_fbpdtv.as_array(),vmin=clims[0],vmax=clims[1])
+a.set_title('FISTA FGP TV')
+imgplot = plt.imshow(x_fista_fgp.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')
+imgplot = plt.loglog(criter_rof , label='ROF TV')
+imgplot = plt.loglog(criter_fgp, label='FGP TV')
b.legend(loc='right')
plt.show()
-#%%