From 0178ae2bf4c84649a28e87ad20f33808ae6eebee Mon Sep 17 00:00:00 2001 From: Jakob Jorgensen Date: Thu, 8 Mar 2018 23:30:48 +0000 Subject: Import also FBPD from algs, move tv demo to wip --- Wrappers/Python/test/simple_demo_tv.py | 100 --------------------------------- Wrappers/Python/wip/simple_demo.py | 2 +- Wrappers/Python/wip/simple_demo_tv.py | 100 +++++++++++++++++++++++++++++++++ 3 files changed, 101 insertions(+), 101 deletions(-) delete mode 100644 Wrappers/Python/test/simple_demo_tv.py create mode 100644 Wrappers/Python/wip/simple_demo_tv.py (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/test/simple_demo_tv.py b/Wrappers/Python/test/simple_demo_tv.py deleted file mode 100644 index 0ae550a..0000000 --- a/Wrappers/Python/test/simple_demo_tv.py +++ /dev/null @@ -1,100 +0,0 @@ - - -import sys - -sys.path.append("..") - -from ccpi.framework import * -from ccpi.reconstruction.algs import * -from ccpi.reconstruction.funcs import * -from ccpi.reconstruction.ops import * -from ccpi.reconstruction.astra_ops import * -from ccpi.reconstruction.geoms import * - -import numpy as np -import matplotlib.pyplot as plt - -test_case = 2 # 1=parallel2D, 2=cone2D - -# Set up phantom -N = 128 - -x = np.zeros((N,N)) -x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 1.0 -x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 2.0 - -plt.imshow(x) -plt.show() - -vg = VolumeGeometry(N,N,None, 1,1,None) - -Phantom = VolumeData(x,geometry=vg) - -# 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 = SinogramGeometry('parallel', - '2D', - angles, - det_num,det_w) -elif test_case==2: - pg = SinogramGeometry('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, 'gpu') - -# 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 = VolumeData(np.zeros(x.shape),geometry=vg) - -# Now least squares plus 1-norm regularization -lam = 1 -g0 = TV2D(lam) - - -# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm -opt = {'tol': 1e-4, 'iter': 10000} -x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt) - -plt.imshow(x_fbpd1.array) -plt.show() - -plt.semilogy(criter_fbpd1) -plt.show() \ No newline at end of file diff --git a/Wrappers/Python/wip/simple_demo.py b/Wrappers/Python/wip/simple_demo.py index 21796b8..da617c1 100644 --- a/Wrappers/Python/wip/simple_demo.py +++ b/Wrappers/Python/wip/simple_demo.py @@ -2,7 +2,7 @@ #sys.path.append("..") from ccpi.framework import VolumeData -from ccpi.reconstruction.algs import FISTA +from ccpi.reconstruction.algs import * from ccpi.reconstruction.funcs import Norm2sq, Norm1 from ccpi.reconstruction.astra_ops import AstraProjectorSimple from ccpi.reconstruction.geoms import VolumeGeometry, SinogramGeometry diff --git a/Wrappers/Python/wip/simple_demo_tv.py b/Wrappers/Python/wip/simple_demo_tv.py new file mode 100644 index 0000000..0ae550a --- /dev/null +++ b/Wrappers/Python/wip/simple_demo_tv.py @@ -0,0 +1,100 @@ + + +import sys + +sys.path.append("..") + +from ccpi.framework import * +from ccpi.reconstruction.algs import * +from ccpi.reconstruction.funcs import * +from ccpi.reconstruction.ops import * +from ccpi.reconstruction.astra_ops import * +from ccpi.reconstruction.geoms import * + +import numpy as np +import matplotlib.pyplot as plt + +test_case = 2 # 1=parallel2D, 2=cone2D + +# Set up phantom +N = 128 + +x = np.zeros((N,N)) +x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 1.0 +x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 2.0 + +plt.imshow(x) +plt.show() + +vg = VolumeGeometry(N,N,None, 1,1,None) + +Phantom = VolumeData(x,geometry=vg) + +# 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 = SinogramGeometry('parallel', + '2D', + angles, + det_num,det_w) +elif test_case==2: + pg = SinogramGeometry('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, 'gpu') + +# 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 = VolumeData(np.zeros(x.shape),geometry=vg) + +# Now least squares plus 1-norm regularization +lam = 1 +g0 = TV2D(lam) + + +# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm +opt = {'tol': 1e-4, 'iter': 10000} +x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt) + +plt.imshow(x_fbpd1.array) +plt.show() + +plt.semilogy(criter_fbpd1) +plt.show() \ No newline at end of file -- cgit v1.2.3