diff options
| -rw-r--r-- | Wrappers/Python/ccpi/framework.py | 10 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/reconstruction/algs.py | 37 | ||||
| -rw-r--r-- | Wrappers/Python/wip/simple_demo.py | 61 |
3 files changed, 106 insertions, 2 deletions
diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py index 5a507d9..b2f8a7e 100644 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework.py @@ -425,6 +425,15 @@ class DataSet(object): repres += "Axis labels: {0}\n".format(self.dimension_labels) repres += "Representation: \n{0}\n".format(self.array) return repres + + def clone(self): + '''returns a copy of itself''' + + return type(self)(self.array, + dimension_labels=self.dimension_labels, + deep_copy=True, + geometry=self.geometry ) + @@ -895,4 +904,5 @@ if __name__ == '__main__': geom_type='parallel', pixel_num_v=3, pixel_num_h=5 , channels=2) sino = SinogramData(geometry=sgeometry) + sino2 = sino.clone()
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/reconstruction/algs.py b/Wrappers/Python/ccpi/reconstruction/algs.py index be0abb5..088b36e 100644 --- a/Wrappers/Python/ccpi/reconstruction/algs.py +++ b/Wrappers/Python/ccpi/reconstruction/algs.py @@ -21,6 +21,7 @@ import numpy import time from ccpi.reconstruction.funcs import BaseFunction +from ccpi.framework import SinogramData, VolumeData def FISTA(x_init, f=None, g=None, opt=None): @@ -126,3 +127,39 @@ def FBPD(x_init, f=None, g=None, h=None, opt=None): timing = numpy.cumsum(timing[0:it+1]); return x, it, timing, criter + +def CGLS(A,b,max_iter,x_init): + '''Conjugate Gradient Least Squares algorithm''' + + r = b.clone() + x = x_init.clone() + + d = A.adjoint(r) + + normr2 = (d**2).sum() + + timing = numpy.zeros(max_iter) + criter = numpy.zeros(max_iter) + + # algorithm loop + for it in range(0, max_iter): + + t = time.time() + + Ad = A.direct(d) + alpha = normr2/( (Ad**2).sum() ) + x = x + alpha*d + r = r - alpha*Ad + s = A.adjoint(r) + + normr2_new = (s**2).sum() + beta = normr2_new/normr2 + normr2 = normr2_new + d = s + beta*d + + # time and criterion + timing[it] = time.time() - t + criter[it] = (r**2).sum() + + return x, it, timing, criter + diff --git a/Wrappers/Python/wip/simple_demo.py b/Wrappers/Python/wip/simple_demo.py index da617c1..766e448 100644 --- a/Wrappers/Python/wip/simple_demo.py +++ b/Wrappers/Python/wip/simple_demo.py @@ -56,7 +56,7 @@ elif test_case==2: dist_center_detector=OrigDetec) # ASTRA operator using volume and sinogram geometries -Aop = AstraProjectorSimple(vg, pg, 'gpu') +Aop = AstraProjectorSimple(vg, pg, 'cpu') # Unused old astra projector without geometry # Aop_old = AstraProjector(det_w, det_num, SourceOrig, @@ -83,6 +83,7 @@ x_init = VolumeData(np.zeros(x.shape),geometry=vg) x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None) plt.imshow(x_fista0.array) +plt.title('FISTA0') plt.show() # Now least squares plus 1-norm regularization @@ -93,6 +94,7 @@ g0 = Norm1(lam) x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0) plt.imshow(x_fista1.array) +plt.title('FISTA') plt.show() plt.semilogy(criter1) @@ -103,7 +105,62 @@ 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.title('FBPD') plt.show() plt.semilogy(criter_fbpd1) -plt.show()
\ No newline at end of file +plt.show() + +# Run CGLS, which should agree with the FISTA0 +x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(Aop, b, 1000, x_init) + +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() + + +#%% +cols = 3 +rows = 2 +current = 1 +fig = plt.figure() +# projections row +a=fig.add_subplot(rows,cols,current) +a.set_title('phantom {0}'.format(numpy.shape(Phantom.as_array()))) +imgplot = plt.imshow(Phantom.as_array()) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('FISTA0') +imgplot = plt.imshow(x_fista0.as_array()) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('FISTA1') +imgplot = plt.imshow(x_fista1.as_array()) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('FBPD') +imgplot = plt.imshow(x_fbpd1.as_array()) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('CGLS') +imgplot = plt.imshow(x_CGLS.as_array()) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('criteria') +imgplot = plt.loglog(criter0 , label='FISTA0') +imgplot = plt.loglog(criter1 , label='FISTA1') +imgplot = plt.loglog(criter_fbpd1, label='FBPD') +imgplot = plt.loglog(criter_CGLS, label='CGLS') +a.legend(loc='right') +plt.show() +#%%
\ No newline at end of file |