From 00b6766211f6fa77a599d6925b6674b0339170fd Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Fri, 11 May 2018 11:06:29 +0100 Subject: Change incorrect file ext to py --- Wrappers/Python/wip/demo_simple_RGLTK.md | 214 ------------------------------- Wrappers/Python/wip/demo_simple_RGLTK.py | 214 +++++++++++++++++++++++++++++++ 2 files changed, 214 insertions(+), 214 deletions(-) delete mode 100644 Wrappers/Python/wip/demo_simple_RGLTK.md create mode 100644 Wrappers/Python/wip/demo_simple_RGLTK.py (limited to 'Wrappers') diff --git a/Wrappers/Python/wip/demo_simple_RGLTK.md b/Wrappers/Python/wip/demo_simple_RGLTK.md deleted file mode 100644 index 9f0a4c3..0000000 --- a/Wrappers/Python/wip/demo_simple_RGLTK.md +++ /dev/null @@ -1,214 +0,0 @@ - -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_ - -import numpy as np -import matplotlib.pyplot as plt - -test_case = 1 # 1=parallel2D, 2=cone2D - -# Set up phantom -N = 128 - - -vg = ImageGeometry(voxel_num_x=N,voxel_num_y=N) -Phantom = ImageData(geometry=vg) - -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() - -# Set up measurement geometry -angles_num = 20; # angles number - -if test_case==1: - angles = np.linspace(0,np.pi,angles_num,endpoint=False) -elif test_case==2: - angles = np.linspace(0,2*np.pi,angles_num,endpoint=False) -else: - NotImplemented - -det_w = 1.0 -det_num = N -SourceOrig = 200 -OrigDetec = 0 - -# 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 -b = Aop.direct(Phantom) -out2 = Aop.adjoint(b) - -#plt.imshow(b.array) -#plt.show() - -#plt.imshow(out2.array) -#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') - -opt = {'tol': 1e-4, 'iter': 100} - -x_fista_rof, it1, timing1, criter_rof = FISTA(x_init, f, g_rof,opt) - -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') - -x_fista_fgp, it1, timing1, criter_fgp = FISTA(x_init, f, g_fgp,opt) - -plt.figure() -plt.subplot(121) -plt.imshow(x_fista_fgp.array,cmap="BuPu") -plt.title('FISTA-FGP-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.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) - -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() - -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, Aop, b, opt ) - -plt.imshow(x_CGLS.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.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]) - -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]) - -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]) - -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') -b.legend(loc='right') -plt.show() -#%% diff --git a/Wrappers/Python/wip/demo_simple_RGLTK.py b/Wrappers/Python/wip/demo_simple_RGLTK.py new file mode 100644 index 0000000..9f0a4c3 --- /dev/null +++ b/Wrappers/Python/wip/demo_simple_RGLTK.py @@ -0,0 +1,214 @@ + +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_ + +import numpy as np +import matplotlib.pyplot as plt + +test_case = 1 # 1=parallel2D, 2=cone2D + +# Set up phantom +N = 128 + + +vg = ImageGeometry(voxel_num_x=N,voxel_num_y=N) +Phantom = ImageData(geometry=vg) + +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() + +# Set up measurement geometry +angles_num = 20; # angles number + +if test_case==1: + angles = np.linspace(0,np.pi,angles_num,endpoint=False) +elif test_case==2: + angles = np.linspace(0,2*np.pi,angles_num,endpoint=False) +else: + NotImplemented + +det_w = 1.0 +det_num = N +SourceOrig = 200 +OrigDetec = 0 + +# 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 +b = Aop.direct(Phantom) +out2 = Aop.adjoint(b) + +#plt.imshow(b.array) +#plt.show() + +#plt.imshow(out2.array) +#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') + +opt = {'tol': 1e-4, 'iter': 100} + +x_fista_rof, it1, timing1, criter_rof = FISTA(x_init, f, g_rof,opt) + +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') + +x_fista_fgp, it1, timing1, criter_fgp = FISTA(x_init, f, g_fgp,opt) + +plt.figure() +plt.subplot(121) +plt.imshow(x_fista_fgp.array,cmap="BuPu") +plt.title('FISTA-FGP-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.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) + +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() + +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, Aop, b, opt ) + +plt.imshow(x_CGLS.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.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]) + +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]) + +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]) + +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') +b.legend(loc='right') +plt.show() +#%% -- cgit v1.2.3