In [ ]:

    

%matplotlib inline


    

In [ ]:

    

import numpy as np 
import pylab as plt 

import scipy

import astra 
import tomopy


import alg 
import sirt 
import sirt_noise


    

In [ ]:

    

plt.gray()


    

In [ ]:

    

# make phantom

size = 128
ref_phantom = tomopy.misc.phantom.shepp2d(size=size).astype('float32')
ref_phantom = np.squeeze(ref_phantom)/np.max(ref_phantom)

pad = 20

phantom = np.zeros(np.asarray(ref_phantom.shape)+2*pad, dtype='float32')
phantom[pad:-pad, pad:-pad] = ref_phantom

# make sinogram
n_angles = 180
angles = np.arange(0.0, 180.0,  180.0 / n_angles)
angles = angles.astype('float32') / 180 * np.pi

pg = astra.create_proj_geom('parallel', 1.0, phantom.shape[0], angles)
vg = astra.create_vol_geom(phantom.shape)
sino = alg.gpu_fp(pg, vg, phantom)
sino = sino.astype('float32')


    

In [ ]:

    

D = np.ones_like(sino)


    

In [ ]:

    

plt.figure(figsize=(16,7))
plt.subplot(131)
plt.imshow(phantom)
plt.colorbar(orientation='horizontal')

plt.subplot(132)
plt.imshow(sino)
plt.axis('auto')
plt.colorbar(orientation='horizontal')

plt.subplot(133)
plt.imshow(D)
plt.axis('auto')
plt.colorbar(orientation='horizontal')

plt.show()


    

In [ ]:

    

mask = np.zeros_like(sino)
mask[:,32+pad:-32-pad] = 1
sino[mask == 0] = 0
D[mask==0] = 0


    

In [ ]:

    

plt.figure(figsize=(16,7))
plt.subplot(131)
plt.imshow(phantom)
plt.colorbar(orientation='horizontal')

plt.subplot(132)
plt.imshow(sino)
plt.axis('auto')
plt.colorbar(orientation='horizontal')

plt.subplot(133)
plt.imshow(D)
plt.axis('auto')
plt.colorbar(orientation='horizontal')

plt.show()


    

In [ ]:

    

# %%timeit
proj_id = astra.create_projector('cuda', pg, vg)
W = astra.OpTomo(proj_id)
eps = 1e-30

x0 = np.zeros_like(phantom)
#x0 = rec_1.copy()
rec = sirt_noise.run(W, sino, D, x0, eps, 100, 'steepest')
en_0 = rec['energy']
alpha_0 = rec['alpha']
rec_0 = rec['rec']

astra.projector.delete(proj_id)


    

In [ ]:

    

plt.figure(figsize=(16,16))
plt.subplot(221)
plt.imshow(rec_0)
plt.colorbar(orientation='horizontal')

plt.subplot(222)
plt.semilogy(en_0)
plt.grid()

plt.subplot(223)
plt.imshow((W*rec_0.ravel()).reshape((len(angles),-1)))
plt.colorbar(orientation='horizontal')
plt.show()


    

In [ ]:

    

pg = astra.create_proj_geom('parallel', 1.0, phantom.shape[0], angles[0])
vg = astra.create_vol_geom(phantom.shape)

x0 = np.zeros_like(phantom)

en_0 =[]
alpha_0 = []

ang_index=np.arange(len(angles))
for iter_numb in range(10):
    np.random.shuffle(ang_index)
    for ang_id in np.array_split(ang_index,len(angles)//10):
        pg['ProjectionAngles'] = angles[ang_id]
        proj_id = astra.create_projector('cuda', pg, vg)
        W = astra.OpTomo(proj_id)
        eps = 1e-30
        rec = sirt_noise.run(W, sino[ang_id], D[ang_id], x0, eps, 1, 'steepest')
        x0 = rec['rec']
        astra.projector.delete(proj_id)
    en_0.extend(rec['energy'])
    alpha_0.extend(rec['alpha'])
rec_0=x0

astra.projector.delete(proj_id)


    

In [ ]:

    

plt.figure(figsize=(16,16))
plt.subplot(221)
plt.imshow(rec_0)
plt.colorbar(orientation='horizontal')

plt.subplot(222)
plt.semilogy(en_0)
plt.grid()

pg['ProjectionAngles'] = angles
proj_id = astra.create_projector('cuda', pg, vg)
W = astra.OpTomo(proj_id)
        
plt.subplot(223)
plt.imshow((W*rec_0.ravel()).reshape((len(angles),-1)))
plt.colorbar(orientation='horizontal')
plt.show()


    

In [ ]: