In [ ]:

    

%reload_ext autoreload
%autoreload 2
%matplotlib inline
# %xmode verbose


    

In [ ]:

    

from glob import glob
import os

import numpy as np
import matplotlib.pyplot as plt
import astra
import alg
import RegVarSIRT
import pickle
import skimage
import skimage.io
import scipy as sp
import ssim
import pandas as pd
import pickle


    

In [ ]:

    

import matplotlib
font = {'size'   : 18}

matplotlib.rc('font', **font)


    

In [ ]:

    

def log_progress(sequence, every=None, size=None, name='Items'):
    from ipywidgets import IntProgress, HTML, VBox
    from IPython.display import display

    is_iterator = False
    if size is None:
        try:
            size = len(sequence)
        except TypeError:
            is_iterator = True
    if size is not None:
        if every is None:
            if size <= 200:
                every = 1
            else:
                every = int(size / 200)     # every 0.5%
    else:
        assert every is not None, 'sequence is iterator, set every'

    if is_iterator:
        progress = IntProgress(min=0, max=1, value=1)
        progress.bar_style = 'info'
    else:
        progress = IntProgress(min=0, max=size, value=0)
    label = HTML()
    box = VBox(children=[label, progress])
    display(box)

    index = 0
    try:
        for index, record in enumerate(sequence, 1):
            if index == 1 or index % every == 0:
                if is_iterator:
                    label.value = '{name}: {index} / ?'.format(
                        name=name,
                        index=index
                    )
                else:
                    progress.value = index
                    label.value = u'{name}: {index} / {size}'.format(
                        name=name,
                        index=index,
                        size=size
                    )
            yield record
    except:
        progress.bar_style = 'danger'
        raise
    else:
        progress.bar_style = 'success'
        progress.value = index
        label.value = "{name}: {index}".format(
            name=name,
            index=str(index or '?')
        )


    

In [ ]:

    

# make phantom
size = 128
mu1 = 0.006
mu2 = 0.005
mu3 = 0.004
phantom = np.zeros((size, size))
half_s = size / 2

y, x = np.meshgrid(range(size), range(size))
xx = (x - half_s).astype('float32')
yy = (y - half_s).astype('float32')
  
mask_ell1 = pow(xx + 0.1*size, 2)/np.power(0.35*size, 2) + pow(yy, 2)/np.power(0.15*size, 2) <= 1
mask_ell2 = pow(xx - 0.15*size, 2)/np.power(0.3*size, 2) + pow(yy - 0.15*size, 2)/np.power(0.15*size, 2) <= 1 
phantom[mask_ell1] = mu1
phantom[mask_ell2] = mu2
phantom[np.logical_and(mask_ell1, mask_ell2)] = mu3
phantom[int(0.15*size):int(0.35*size), int(0.2*size):int(0.5*size)] = mu3 
phantom[int(0.20*size):int(0.25*size), int(0.25*size):int(0.3*size)] = 0
phantom[int(0.30*size):int(0.35*size), int(0.35*size):int(0.4*size)] = mu1*10
phantom = 1e+1 * phantom

# make sinogram
n_angles = size//2
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, size, angles)
vg = astra.create_vol_geom((size, size))
sino = alg.gpu_fp(pg, vg, phantom)
print(sino.min(), sino.max())


    

In [ ]:

    

def reshape_sinogram(sinogram, x_resize, angle_resize, order=1):
    sinogram = sinogram[::angle_resize,:]
    if not x_resize == 1:
        sinogram =sp.ndimage.interpolation.zoom(sinogram,[1, 1./x_resize], order=order)
#     print(sinogram.shape)
    sinogram = sinogram/x_resize
    return sinogram
    
def load_data(path_sino, resize_x=None, resize_angle=None, order=1):
    if resize_x is None:
        x_resize = 4
    else:
        x_resize = resize_x
    
    if resize_angle is None:
        angle_resize = 4
    else:
        angle_resize = resize_angle
    
    sinogram = plt.imread(path_sino)
    if len(sinogram.shape) == 3:
        sinogram = sinogram[...,0]
    sinogram = np.flipud(sinogram)
    old_n_angles = sinogram.shape[0] 
    sinogram = reshape_sinogram(sinogram, x_resize, angle_resize, order=order)
     
    
    detector_cell = sinogram.shape[1]
    n_angles = sinogram.shape[0] 

    pixel_size = 2.50e-3*x_resize
    os_distance = 49.430 / pixel_size
    ds_distance = 225.315 / pixel_size
    angle_step = 0.2 / 180.0 * np.pi
    
    angles = np.arange(n_angles) * angle_step * angle_resize  
    angles = angles - (angles.max() + angles.min()) / 2

    if old_n_angles%angle_resize == 0:
        angles = angles - angle_step * (angle_resize-1) / 2.
    else:
        angles = angles - angle_step * (old_n_angles%angle_resize-1) / 2.
#     angles = angles - angle_step * (angle_resize-1) / 2.
    angles = angles + np.pi / 2

    vol_geom = astra.create_vol_geom(detector_cell, detector_cell)
    proj_geom = astra.create_proj_geom('fanflat', ds_distance / os_distance, detector_cell, angles,
                                                                        os_distance, (ds_distance - os_distance))
    return proj_geom, vol_geom, sinogram


    

In [ ]:

    

sinograms = glob('/home/makov/diskmnt/big/yaivan/noisy_data/noised_sino/BHI3_2.49um_*__sino0980.tif')
sinograms


    

In [ ]:

    

plt.gray()


    

In [ ]:

    

from functools import lru_cache


    

In [ ]:

    

@lru_cache(8)
def get_reference_rec():
    name = os.path.split(sinograms[0])[-1]
    pg, vg , sino_noise = load_data(sinograms[0], 1, 1)
    k = sino_noise.shape[0]/pg['DetectorWidth']**2/(np.pi/2)
    proj_id = astra.create_projector('cuda', pg, vg)
    W = astra.OpTomo(proj_id)
    rec_reference = alg.gpu_fbp(pg, vg, sino_noise)*k
    return rec_reference


    

In [ ]:

    

plt.figure(figsize=(7,5))
plt.hist(get_reference_rec().ravel(), bins=1000);
plt.vlines(7,0,10000)
plt.grid()


    

In [ ]:

    

def diff(x0,x1,norm):
    mask = x1>7
    d0 = x0[mask>0].ravel()
    d1 = x1[mask>0].ravel()
    return (sp.linalg.norm(d0-d1,norm)+sp.linalg.norm(x0[mask==0].ravel(),norm))/x0.size


    

In [ ]:

    

def l2_grad(data):
    g_data = np.gradient(data)
    return sp.linalg.norm((g_data[0]**2+g_data[1]**2),1)**0.5/data.size


    

In [ ]:

    

from skimage.measure import compare_ssim as ssim


    

In [ ]:

    

X,Y=np.meshgrid(np.arange(4000)-2000, np.arange(4000)-2000)
mask = (X**2+Y**2 <2000**2).astype('float32')


    

In [ ]:

    

sinograms[0]


    

In [ ]:

    

# import alg
# pg,vg, vsino = load_data(sinogram, 1,1)
# vsino = vsino.astype('float32')
# proj_id = astra.create_projector('cuda', pg, vg)
# rec = alg.gpu_cgls(pg,vg, vsino, 500)
# # print()
# # D = 1.0 * (0.01+ (vsino/np.percentile(vsino, 90))**2)
# # D /= D.mean()


    

In [ ]:

    

# plt.figure(figsize=(14,10))
# # plt.subplot(121)
# plt.imshow(rec[1200:1600,800:1200], vmin=0, vmax=75)
# plt.colorbar()
# # plt.axis('tight')
# # plt.subplot(122)
# # plt.imshow(vsino)
# # plt.colorbar()
# # plt.axis('tight')
# plt.show()


    

In [ ]:

    

rec_reference = get_reference_rec()
def show_pic(rec_fbp, rec_rec, rec_v, rec_reference, name = None):
    plt.figure(figsize=(14,14))
    plt.subplot(231)
    plt.title('FBP')
    plt.imshow(rec_fbp[1200:1600,800:1200], vmin=0, vmax=75)
    plt.subplot(232)
    plt.title('VCGLS+TV')
    plt.imshow(rec_rec[1200:1600,800:1200], vmin=0, vmax=75)
    plt.subplot(233)
    plt.title('Reference')
    plt.imshow(rec_reference[1200:1600,800:1200], vmin=0, vmax=75)
    plt.title('Reference image')

    plt.subplot(223)
    plt.plot(rec_fbp[1200:1600,1000], label='FBP')
    plt.plot(rec_rec[1200:1600,1000], label='VCGLS+TV')
    plt.plot(rec_reference[1200:1600,1000], label='Reference')
    plt.legend(loc=0)
    plt.grid()

    plt.subplot(224)
    plt.plot(rec_fbp[1400,800:1200], label='FBP')
    plt.plot(rec_rec[1400,800:1200], label='VCGLS+TV')
    plt.plot(rec_reference[1400,800:1200], label='Reference')
    plt.legend(loc=0)
    plt.grid()
    if name is not None:
        plt.savefig('{}.png'.format(name))
    plt.show()

#     plt.figure(figsize=(12,7))
#     plt.subplot(121)
#     plt.plot(rec_v['alpha'], label='Alpha')
#     plt.legend()
#     plt.grid()
#     plt.subplot(122)
#     plt.plot(rec_v['energy'], label='Energy')
#     plt.grid()
#     plt.legend()
#     plt.show()


    

In [ ]:

    

from skimage import filters
x = get_reference_rec().copy()
x[x>np.percentile(x,80)] = 0
val = filters.threshold_otsu(x)

print(val)


    

In [ ]:

    

from skimage import filters
import alg
import RegVarSIRT
def make_vcgls_TV(sinogram, x_resize, ang_resize, n_iter=10):
    name = os.path.split(sinogram)[-1]
    short_name=name.split('_')[2]
    pg,vg, vsino = load_data(sinogram, x_resize, ang_resize)
    vsino = vsino.astype('float32')
    proj_id = astra.create_projector('cuda', pg, vg)
    print()
    D = 1.0 * (0.1 + (vsino/np.percentile(vsino, 90))**2)
    D /= D.mean()
#     D = np.ones_like(vsino)
    Div = (1.0 / D).astype('float32')
    
    W = astra.OpTomo(proj_id)
    x0 = np.zeros((vsino.shape[1], vsino.shape[1]),dtype='float32')

    k = vsino.shape[0]/pg['DetectorWidth']**2/(np.pi/2)
    rec_fbp = alg.gpu_fbp(pg, vg, vsino)*k  # fixit
    
#     x0 = rec_fbp
#     x = x0.copy()
#     x[x>np.percentile(x,80)] = np.percentile(x,80)
#     val = filters.threshold_otsu(x)
    
    for i in range(50):
        rec_v = RegVarSIRT.run(W, vsino, x0, Div, Lambda=0, eps=1.e-10, n_it=n_iter, step='CG', normalize=True)
        rec_rec = rec_v['rec']
        x0 = rec_rec
#         x0[x0<val] = 0
        show_pic(rec_fbp, rec_rec, rec_v, get_reference_rec(), '{}_{}'.format(short_name, i))
    
#     if x_resize>=4 and ang_resize>=2:
#         rec_v = RegVarSIRT.run(W, vsino, x0, Lambda=100, eps=1.e-10, n_it=n_iter, step='CG', normalize=True)
#         rec_rec = rec_v['rec']


#         for l in log_progress([100, 100, 100, 50, 50, 10, 10, 5, 4, 3, 2, 1]):
#             rec_v = RegVarSIRT.run(W, vsino, rec_fbp, Div, Lambda=l, eps=1.e-10, n_it=n_iter, step='CG', normalize=True)
#             rec_rec = rec_v['rec']
#     else:
#         rec_rec = np.zeros_like(rec_fbp)
    return rec_rec, rec_fbp, rec_v


    

In [ ]:

    

rec_reference = get_reference_rec()
for sinogram in sinograms:
    print(sinogram)
#     rec_rec, rec_fbp, rec_v = 
    _, _, _ =make_vcgls_TV(sinogram, 1, 1, n_iter=10)

# plt.figure(figsize=(14,14))
# plt.subplot(231)
# plt.title('FBP')
# plt.imshow(rec_fbp[1200:1600,800:1200], vmin=0, vmax=75)
# plt.subplot(232)
# plt.title('VCGLS+TV')
# plt.imshow(rec_rec[1200:1600,800:1200], vmin=0, vmax=75)
# plt.subplot(233)
# plt.title('Reference')
# plt.imshow(rec_reference[1200:1600,800:1200], vmin=0, vmax=75)
# plt.title('Reference image')

# plt.subplot(223)
# plt.plot(rec_fbp[1200:1600,1000], label='FBP')
# plt.plot(rec_rec[1200:1600,1000], label='VCGLS+TV')
# plt.plot(rec_reference[1200:1600,1000], label='Reference')
# plt.legend(loc=0)
# plt.grid()

# plt.subplot(224)
# plt.plot(rec_fbp[1400,800:1200], label='FBP')
# plt.plot(rec_rec[1400,800:1200], label='VCGLS+TV')
# plt.plot(rec_reference[1400,800:1200], label='Reference')
# plt.legend(loc=0)
# plt.grid()
# plt.show()

# plt.figure(figsize=(12,7))
# plt.subplot(121)
# plt.plot(rec_v['alpha'], label='Alpha')
# plt.legend()
# plt.grid()
# plt.subplot(122)
# plt.plot(rec_v['energy'], label='Energy')
# plt.grid()
# plt.legend()
# plt.show()


# plt.figure(figsize=(14,10))
# plt.subplot(131)
# plt.imshow((rec_fbp-rec_rec)[1200:1600,800:1200], vmin=-5, vmax=5)
# plt.title('FPB-VCGLS')
# plt.subplot(132)
# plt.title('REF-VCGLS')
# plt.imshow((rec_reference-rec_rec)[1200:1600,800:1200], vmin=-5, vmax=5)
# plt.subplot(133)
# plt.imshow((rec_reference-rec_fbp)[1200:1600,800:1200], vmin=-5, vmax=5)
# plt.title('FBP-REF')
# plt.show()

# plt.figure(figsize=(10,10))
# plt.imshow((rec_fbp-rec_rec), vmin=-5, vmax=5)
# plt.title('FPB-VCGLS')
# plt.show()

# plt.figure(figsize=(10,10))
# plt.title('REF-VCGLS')
# plt.imshow((rec_reference-rec_rec), vmin=-5, vmax=5)
# plt.show()

# plt.figure(figsize=(10,10))
# plt.imshow((rec_reference-rec_fbp), vmin=-5, vmax=5)
# plt.title('FBP-REF')
# plt.show()

# plt.figure(figsize=(10,10))
# plt.imshow(rec_reference, vmin=0, vmax=75)
# plt.title('REF')
# plt.show()


    

In [ ]:

    

# plt.figure(figsize=(14,14))
# plt.subplot(231)
# plt.title('FBP')
# plt.imshow(rec_fbp[1200:1600,800:1200], vmin=0, vmax=75)
# plt.subplot(232)
# plt.title('VCGLS+TV')
# plt.imshow(rec_rec[1200:1600,800:1200], vmin=0, vmax=75)
# plt.subplot(233)
# plt.title('Reference')
# plt.imshow(rec_reference[1200:1600,800:1200], vmin=0, vmax=75)
# plt.title('Reference image')

# plt.subplot(223)
# plt.plot(rec_fbp[1200:1600,1000], label='FBP')
# plt.plot(rec_rec[1200:1600,1000], label='VCGLS+TV')
# plt.plot(rec_reference[1200:1600,1000], label='Reference')
# plt.legend(loc=0)
# plt.grid()

# plt.subplot(224)
# plt.plot(rec_fbp[1400,800:1200], label='FBP')
# plt.plot(rec_rec[1400,800:1200], label='VCGLS+TV')
# plt.plot(rec_reference[1400,800:1200], label='Reference')
# plt.legend(loc=0)
# plt.grid()
# plt.show()

# plt.figure(figsize=(12,7))
# plt.subplot(121)
# plt.plot(rec_v['alpha'], label='Alpha')
# plt.legend()
# plt.grid()
# plt.subplot(122)
# plt.plot(rec_v['energy'], label='Energy')
# plt.grid()
# plt.legend()
# plt.show()


# plt.figure(figsize=(14,10))
# plt.subplot(131)
# plt.imshow((rec_fbp-rec_rec)[1200:1600,800:1200], vmin=-5, vmax=5)
# plt.title('FPB-VCGLS')
# plt.subplot(132)
# plt.title('REF-VCGLS')
# plt.imshow((rec_reference-rec_rec)[1200:1600,800:1200], vmin=-5, vmax=5)
# plt.subplot(133)
# plt.imshow((rec_reference-rec_fbp)[1200:1600,800:1200], vmin=-5, vmax=5)
# plt.title('FBP-REF')
# plt.show()

# plt.figure(figsize=(10,10))
# plt.imshow((rec_fbp-rec_rec), vmin=-5, vmax=5)
# plt.title('FPB-VCGLS')
# plt.show()

# plt.figure(figsize=(10,10))
# plt.title('REF-VCGLS')
# plt.imshow((rec_reference-rec_rec), vmin=-5, vmax=5)
# plt.show()

# plt.figure(figsize=(10,10))
# plt.imshow((rec_reference-rec_fbp), vmin=-5, vmax=5)
# plt.title('FBP-REF')
# plt.show()

plt.figure(figsize=(10,10))
plt.imshow(rec[1200:1600,800:1200], vmin=0, vmax=75)
plt.title('REC')
plt.show()


    

In [ ]:

    




    

In [ ]:

    

res={}
x_range = [1,]
ang_range = [1,]

rec_reference = get_reference_rec()*mask
for sinogram in log_progress([sinograms[0],]):
    name = os.path.split(sinogram)[-1]
    short_name=name.split('_')[2]
    res[short_name] = {}
    for method in ['fbp', 'vcgls']:
        res[short_name][method] = {}
#         res[short_name][method]['shennon'] = np.zeros((len(x_range), len(ang_range)))
        res[short_name][method]['l2'] = np.zeros((len(x_range), len(ang_range)))
#         res[short_name][method]['l2_grad'] = np.zeros((len(x_range), len(ang_range)))
        res[short_name][method]['ssim'] = np.zeros((len(x_range), len(ang_range)))
        res[short_name][method]['corr'] = np.zeros((len(x_range), len(ang_range)))
        res[short_name][method]['rec'] = np.zeros((len(x_range), len(ang_range),4000,4000))
        
    for ix, x_resize  in enumerate(x_range):
        for iang, ang_resize in enumerate(ang_range):
            print(short_name, x_resize, ang_resize)
            
            n_iter=10
                
            rec_rec, rec_fbp = make_vcgls_TV(sinogram, x_resize, ang_resize, n_iter=n_iter)
            
            if not x_resize == 1:
                rec_fbp = sp.ndimage.interpolation.zoom(rec_fbp, x_resize, order=2)
            rec_fbp = rec_fbp.astype('float32')
            rec_fbp*=mask
            
            if not x_resize == 1:
                rec_rec = sp.ndimage.interpolation.zoom(rec_rec, x_resize, order=2)
            rec_rec = rec_rec.astype('float32')
            rec_rec*=mask
            
            for method in ['fbp', 'vcgls']:
                if method == 'fbp':
                    t_rec = rec_fbp
                elif method == 'vcgls':
                    t_rec = rec_rec
                    
#                 res[short_name][method]['shennon'][ix, iang] = skimage.measure.shannon_entropy(
#                     (t_rec-t_rec.min())/(t_rec.max()-t_rec.min())) 
                res[short_name][method]['corr'][ix, iang] = np.correlate(t_rec.ravel(),rec_reference.ravel())
                res[short_name][method]['l2'][ix, iang] = diff(t_rec, rec_reference, 2)
#                 res[short_name][method]['l2_grad'][ix, iang] = l2_grad(t_rec)
                res[short_name][method]['ssim'][ix, iang] = ssim(t_rec, rec_reference)
                res[short_name][method]['rec'][ix, iang] = t_rec
            
            if True:
                plt.figure(figsize=(14,14))
                plt.subplot(231)
                plt.title('FBP')
                plt.imshow(rec_fbp[1200:1600,800:1200], vmin=0, vmax=75)
                plt.subplot(232)
                plt.title('VCGLS+TV')
                plt.imshow(rec_rec[1200:1600,800:1200], vmin=0, vmax=75)
                plt.subplot(233)
                plt.title('Reference')
                plt.imshow(rec_reference[1200:1600,800:1200], vmin=0, vmax=75)
                plt.title('Reference image')

                plt.subplot(223)
                plt.plot(rec_fbp[1200:1600,1000], label='FBP')
                plt.plot(rec_rec[1200:1600,1000], label='VCGLS+TV')
                plt.plot(rec_reference[1200:1600,1000], label='Reference')
                plt.legend(loc=0)
                plt.grid()

                plt.subplot(224)
                plt.plot(rec_fbp[1400,800:1200], label='FBP')
                plt.plot(rec_rec[1400,800:1200], label='VCGLS+TV')
                plt.plot(rec_reference[1400,800:1200], label='Reference')
                plt.legend(loc=0)
                plt.grid()
#                 plt.savefig('view_{}_{}_{}.png'.format(short_name, x_resize, ang_resize))
                #     plt.close()
                plt.show()
            else:
                plt.figure(figsize=(14,14))
                plt.subplot(221)
                plt.title('FBP')
                plt.imshow(rec_fbp[1200:1600,800:1200], vmin=0, vmax=75)
                plt.subplot(222)
                plt.title('Reference')
                plt.imshow(rec_reference[1200:1600,800:1200], vmin=0, vmax=75)
                plt.title('Reference image')

                plt.subplot(223)
                plt.plot(rec_fbp[1200:1600,1000], label='FBP')
                plt.plot(rec_reference[1200:1600,1000], label='Reference')
                plt.legend(loc=0)
                plt.grid()

                plt.subplot(224)
                plt.plot(rec_fbp[1400,800:1200], label='FBP')
                plt.plot(rec_reference[1400,800:1200], label='Reference')
                plt.legend(loc=0)
                plt.grid()
#                 plt.savefig('view_{}_{}_{}.png'.format(short_name, x_resize, ang_resize))
                #     plt.close()
                plt.show()


    

In [ ]:

    

np.correlate(t_rec.ravel(),rec_reference.ravel())


    

In [ ]:

    

df = pd.DataFrame(columns=['Name', 'Method', 'Angles reduction', 'Detector reduction',
                           'SSIM', 'L2', 'Corr',
                           'Time', 'Exposure', 'Data', 'Q'])
for k,v  in res.items():
    kt = np.prod([float(n) for n in k.split('x')])
    exposure = kt
    for kk, vv in v.items():
        for (x,y), r in np.ndenumerate(vv['ssim']):
            if np.mean(vv['rec'][x,y])==0:
                continue
                
            df.loc[len(df)] = ([k, kk, float(ang_range[y]), float(x_range[x]), 
                                r, vv['l2'][x,y], vv['corr'][x,y],
                                kt/ang_range[y], exposure, vv['rec'][x,y],
#                                 (r*vv['shennon'][x,y]+1)**0.25-1
#                                 -(vv['l2'][x,y])*1e6*vv['l2_grad'][x,y]
                                vv['corr'][x,y]*r
                               ]
                              )


    

In [ ]:

    

df.to_pickle('big_data.pkl')


    

In [ ]:

    

df.sort_values(by='Q', ascending=False, inplace=True);
numb = 0
for x in list(df.iterrows())[:]:
    title = 'Q:{:0.4f} {} {} D:{} A:{}'.format(
        x[1]['Q'],
        x[1]['Method'],
        x[1]['Name'],
        int(x[1]['Detector reduction']),
        int(x[1]['Angles reduction'])        
    )
    numb = numb+1
#     plt.figure(figsize=(14,14))
#     plt.subplot(221)
#     plt.title(title)
#     plt.imshow(x[1]['Data'][1200:1600,800:1200], vmin=0, vmax=75)
#     plt.subplot(222)
#     plt.imshow(rec_reference[1200:1600,800:1200], vmin=0, vmax=75)
#     plt.title('Reference image')
    
#     plt.subplot(223)
#     plt.plot(x[1]['Data'][1200:1600,1000], label='Reconstruction')
#     plt.plot(rec_reference[1200:1600,1000], label='Reference')
#     plt.legend(loc=0)
#     plt.grid()
    
#     plt.subplot(224)
#     plt.plot(x[1]['Data'][1400,800:1200], label='Reconstruction')
#     plt.plot(rec_reference[1400,800:1200], label='Reference')
#     plt.legend(loc=0)
#     plt.grid()
# #     plt.savefig(title+'.png')
# #     plt.close()
#     plt.show()
    skimage.io.imsave('Pic.A1.{}_{}.tif'.format(numb, title.replace(' ','_').replace(':','_')), x[1]['Data'].astype('float32'))


    

In [ ]:

    

skimage.io.imsave('Pic.5.1.1_refrerence_fbp_rec.tif',rec_reference)


    

In [ ]:

    

x[1]['Data'].dtype


    

In [ ]:

    

X,Y=np.meshgrid(np.arange(4000)-2000, np.arange(4000)-2000)
mask = X**2+Y**2 <2000**2


    

In [ ]:

    

plt.figure(figsize=(10,10))
plt.imshow(x[1]['Data']*mask, vmin=0, vmax=75)


    

In [ ]:

    

sinograms


    

In [ ]:

    

# def make_vcgls_TV(sinogram, x_resize, ang_resize):
#     pg,vg, vsino = load_data(sinogram, x_resize, ang_resize)
#     proj_id = astra.create_projector('cuda', pg, vg)
#     D = 3000.0 * (0.05 + (vsino/vsino.max())**2)
#     Div = 1.0 / D

#     W = astra.OpTomo(proj_id)
#     x0 = np.zeros((vsino.shape[1], vsino.shape[1]))

#     k = vsino.shape[0]/pg['DetectorWidth']**2/(np.pi/2)
#     rec_fbp = alg.gpu_fbp(pg, vg, vsino)*k  # fixit

#     rec_v = RegVarSIRT.run(W, vsino, x0, Lambda=100, eps=1.e-10, n_it=10, step='CG', normalize=True)
#     rec_rec = rec_v['rec']


#     for l in log_progress([100, 100, 100, 50, 50, 10, 10, 5, 4, 3, 2, 1]):
#         rec_v = RegVarSIRT.run(W, vsino, rec_fbp, Div, Lambda=l, eps=1.e-10, n_it=10, step='CG', normalize=True)
#         rec_rec = rec_v['rec']

#     return rec_rec, rec_fbp

# rec_rec, rec_fbp = make_vcgls_TV('/home/makov/diskmnt/big/yaivan/noisy_data/noised_sino/BHI3_2.49um_1x0.7__sino0980.tif',
#                                  4,4)

# rec_rec = sp.ndimage.interpolation.zoom(rec_rec, 4, order=2)
# rec_rec*=mask

# rec_fbp = sp.ndimage.interpolation.zoom(rec_fbp, 4, order=2)
# rec_fbp*=mask
            
# plt.figure(figsize=(14,14))
# plt.subplot(231)
# plt.title('FBP')
# plt.imshow(rec_fbp[1200:1600,800:1200], vmin=0, vmax=75)
# plt.subplot(232)
# plt.title('VCGLS+TV')
# plt.imshow(rec_rec[1200:1600,800:1200], vmin=0, vmax=75)
# plt.subplot(233)
# plt.title('Reference')
# plt.imshow(rec_reference[1200:1600,800:1200], vmin=0, vmax=75)
# plt.title('Reference image')

# plt.subplot(223)
# plt.plot(rec_fbp[1200:1600,1000], label='FBP')
# plt.plot(rec_rec[1200:1600,1000], label='VCGLS+TV')
# plt.plot(rec_reference[1200:1600,1000], label='Reference')
# plt.legend(loc=0)
# plt.grid()

# plt.subplot(224)
# plt.plot(rec_fbp[1400,800:1200], label='FBP')
# plt.plot(rec_rec[1400,800:1200], label='VCGLS+TV')
# plt.plot(rec_reference[1400,800:1200], label='Reference')
# plt.legend(loc=0)
# plt.grid()
# #     plt.savefig(title+'.png')
# #     plt.close()
# plt.show()


    

In [ ]:

    

plt.figure(figsize=(10,10))
plt.imshow(rec_rec, vmin=0, vmax=75)


    

In [ ]:

    

from skimage.filters import median


    

In [ ]:

    

plt.figure(figsize=(10,10))
plt.imshow(rec_fbp, vmin=0, vmax=75)


    

In [ ]:

    

plt.figure(figsize=(14,14))
plt.subplot(221)
plt.title('FBP')
plt.imshow(rec_fbp[1200:1600,800:1200], vmin=0, vmax=75)
plt.subplot(222)
plt.title('VCGLS+TV')
plt.imshow(rec_rec[1200:1600,800:1200], vmin=0, vmax=75)
plt.title('Reference image')

plt.subplot(223)
plt.plot(rec_fbp[1200:1600,1000], label='FBP')
plt.plot(rec_rec[1200:1600,1000], label='CGLS')
plt.legend(loc=0)
plt.grid()

plt.subplot(224)
plt.plot(rec_fbp[1400,800:1200], label='FBP')
plt.plot(rec_rec[1400,800:1200], label='CGLS')
plt.legend(loc=0)
plt.grid()
#     plt.savefig(title+'.png')
#     plt.close()
plt.show()


    

In [ ]:

    

df.sort_values(by='Q', ascending=False, inplace=True);


    

In [ ]:

    

def df_to_markdown(df, float_format='%.3g'):
    """
    Export a pandas.DataFrame to markdown-formatted text.
    DataFrame should not contain any `|` characters.
    """
    from os import linesep
    return linesep.join([
        '|'.join(df.columns),
        '|'.join(4 * '-' for i in df.columns),
        df.to_csv(sep='|', index=False, header=False, float_format=float_format)
    ]).replace('|', ' | ')


    

In [ ]:

    

# pdd = pd.DataFrame(df, columns=['Q', 'Name',  'Angles reduction', 'Detector reduction',
#                            'SSIM', 'L2', 'L2_grad', 'Shennon', 
#                            'Time', 'Exposure'])
# print(df_to_markdown(pdd))
# pdd


    

In [ ]:

    

df.sort_values(by='SSIM', ascending=False, inplace=True);

plt.figure(figsize=(10,10))
tdf = df.loc[df['Method'] == 'fbp']
plt.scatter(tdf['Exposure'], tdf['Q'],
        c=tdf['Angles reduction'], s = tdf['Detector reduction']*10,
        cmap=plt.cm.viridis)

tdf = df.loc[df['Method'] == 'vcgls']
plt.scatter(tdf['Exposure'], tdf['Q'],
        c=tdf['Angles reduction'], s = tdf['Detector reduction']*10,
        marker = '<',
        cmap=plt.cm.viridis)

plt.xlabel('Exposure')
plt.ylabel('Quality')
plt.grid()
plt.colorbar(label = 'Angles reduction')

plt.show()


    

In [ ]:

    

df.sort_values(by='SSIM', ascending=False, inplace=True);
plt.gray()
plt.figure(figsize=(10,10))
ttdf = df.loc[df['Time'] < 1]
tdf = ttdf.loc[df['Method'] == 'fbp']
plt.scatter(tdf['Time'], tdf['Q'],
        c=tdf['Angles reduction'], s = tdf['Detector reduction']*10,
        cmap=plt.cm.gray, label = 'FBP')

tdf = ttdf.loc[df['Method'] == 'vcgls']
plt.scatter(tdf['Time'], tdf['Q'],
        c=tdf['Angles reduction'], s = tdf['Detector reduction']*10,
        marker = 's',
        cmap=plt.cm.gray, label = 'VCGLS+TV')

plt.xlabel('Time')
plt.ylabel('Quality')
plt.grid()
plt.colorbar(label = 'Angles reduction')
plt.legend(loc=0)
plt.savefig('Pic_2.png')
plt.show()


    

In [ ]:

    

df.sort_values(by='SSIM', ascending=False, inplace=True);
plt.gray()
plt.figure(figsize=(10,10))
ttdf = df
tdf = ttdf.loc[df['Method'] == 'fbp']
plt.scatter(tdf['Time'], tdf['Q'],
        c=tdf['Angles reduction'], s = tdf['Detector reduction']*10,
        cmap=plt.cm.gray, label = 'FBP')

tdf = ttdf.loc[df['Method'] == 'vcgls']
plt.scatter(tdf['Time'], tdf['Q'],
        c=tdf['Angles reduction'], s = tdf['Detector reduction']*10,
        marker = 's',
        cmap=plt.cm.gray, label = 'VCGLS+TV')

plt.xlabel('Time')
plt.ylabel('Quality')
plt.grid()
plt.colorbar(label = 'Angles reduction')
plt.legend(loc=0)
plt.savefig('Pic_1.png')
plt.show()


    

In [ ]:

    

df.sort_values(by='SSIM', ascending=False, inplace=True);
plt.gray()
plt.figure(figsize=(10,10))
ttdf = df
tdf = ttdf.loc[df['Method'] == 'fbp']
plt.scatter(tdf['Exposure'], tdf['Q'],
        c=tdf['Angles reduction'], s = tdf['Detector reduction']*10,
        cmap=plt.cm.gray, label = 'FBP')

tdf = ttdf.loc[df['Method'] == 'vcgls']
plt.scatter(tdf['Exposure'], tdf['Q'],
        c=tdf['Angles reduction'], s = tdf['Detector reduction']*10,
        marker = 's',
        cmap=plt.cm.gray, label = 'VCGLS+TV')

plt.xlabel('Exposure')
plt.ylabel('Quality')
plt.grid()
plt.colorbar(label = 'Angles reduction')
plt.legend(loc=0)
plt.show()


    

In [ ]:

    

df.sort_values(by='SSIM', ascending=False, inplace=True);

plt.figure(figsize=(10,10))
tdf = df.loc[df['Method'] == 'fbp']
for ix,x in tdf.iterrows():
    t = df.loc[df['Method'] == 'vcgls']
    t = t.loc[t['Name'] == x['Name']]
    t = t.loc[t['Angles reduction'] == x['Angles reduction']]
    t = t.loc[t['Detector reduction'] == x['Detector reduction']]
    if len(t)>0:
        plt.scatter(x['Q'], np.asarray(t['Q']), s=50, c=t['Exposure'])
        
plt.plot(np.arange(15,50),np.arange(15,50))
plt.xlabel('FBP quality')
plt.ylabel('VCGLS+TV quality')
plt.grid()
plt.colorbar(label = 'Time')
plt.savefig('fbp_vcgls.png')
plt.show()
# plt.scatter(df['Time'], df['Q'],
#         c=df['Angles reduction'], s = df['Detector reduction']*10,
#         cmap=plt.cm.viridis)

# plt.xlabel('Time')
# plt.ylabel('Quality')
# plt.grid()
# plt.colorbar(label = 'Angles reduction')

# plt.show()


    

In [ ]:

    

np.asarray(t['Q'])


    

In [ ]:

    

t = df.loc[df['Method'] == 'vcgls']
t = t.loc[t['Name'] == x['Name']]
t = t.loc[t['Angles reduction'] == x['Angles reduction']]
x['Angles reduction']
t


    

In [ ]:

    

tdf.


    

In [ ]:

    

df.sort_values(by='SSIM', ascending=False, inplace=True);

plt.figure(figsize=(10,10))
plt.scatter(df['SSIM'], df['Q'],
            cmap=plt.cm.viridis)

plt.xlabel('SSIM quality')
plt.ylabel('Q')
plt.grid()
# plt.colorbar(label = 'Angles reduction')

plt.show()


    

In [ ]:

    

df.sort_values(by='SSIM', ascending=False, inplace=True);

plt.figure(figsize=(10,10))
plt.scatter(df['SSIM'], df['L2'],
            cmap=plt.cm.viridis)

plt.xlabel('SSIM quality')
plt.ylabel('L2')
plt.grid()
# plt.colorbar(label = 'Angles reduction')

plt.show()


    

In [ ]:

    

df.sort_values(by='SSIM', ascending=False, inplace=True);

plt.figure(figsize=(10,10))
plt.scatter(df['SSIM'], df['L2_grad']*1000,
            cmap=plt.cm.viridis)

plt.xlabel('SSIM quality')
plt.ylabel('L2_grad * 1000')
plt.grid()
# plt.colorbar(label = 'Angles reduction')

plt.show()


    

In [ ]:

    

df.sort_values(by='SSIM', ascending=False, inplace=True);

plt.figure(figsize=(10,10))
plt.scatter(df['SSIM'], df['Time'],
            s=20*df['Detector reduction'], c=df['Angles reduction'],
            cmap=plt.cm.viridis)

plt.xlabel('SSIM quality')
plt.ylabel('Time, a.u.')
plt.grid()
plt.colorbar(label = 'Angles reduction')

plt.show()


    

In [ ]:

    

import matplotlib
font = {'family' : 'normal',
        'weight' : 'normal',
        'size'   : 18}

matplotlib.rc('font', **font)



plt.figure(figsize=(10,10))
plt.scatter(df['L2'], df['Time'],
            s=20*df['Detector reduction'], c=df['Angles reduction'],
            cmap=plt.cm.viridis)

plt.xlabel('L2 quality')
plt.ylabel('Time, a.u.')
plt.grid()
plt.colorbar(label = 'Angles reduction')

plt.show()


    

In [ ]:

    

df.sort_values(by='Q', ascending=False, inplace=True);


    

In [ ]:

    

files = glob('Q:*.png')
files = sorted(files, reverse=True)
# print(files)

with open('Appendix.md', 'w') as f:
    f.write('#Приложение 1.#Приложение 1. Изображения участка реконструированной области отсортированные по качеству реконструкции.\n\n')
    for ifile, file in enumerate(files):
        Q, method, name, D,A = file.split(' ')
        Q = Q[2:]
        D = D[2:]
        A = A[2:-4]
        method = 'FBP' if 'fbp' in method else 'VCGLS+TV'
        legend = '![]({})\n\n'+\
                'Рисунок A1.{}. Сравнение участка реконструированной области методом {} с опорным восстановлением.'+ \
                '(Верх лево) Название эксперимента: {}. Достигнутое качество: {}. '+ \
                'Прореживание по углам: {}. Усреднение по пикселям детектора: {}. '+ \
                '(Верх право) Опорное восстановление.'+ \
                '(Низ) Вертикальное и горизонтальное сечение приведённых реконструкций.\n\n'
        f.write(legend.format(file, ifile+1, method,  name, Q, A,D))


    

In [ ]:

    

for k,v  in res.items():
    for kk, vv in v.items():
        name = '{}_{}'.format(k,kk)
        plt.figure(figsize=(10,8))
        X,Y=np.meshgrid(ang_range, x_range)
        plt.pcolormesh(X, Y, vv['ssim']**2*vv['shennon'], cmap=plt.cm.gray_r)
        plt.ylabel('Detector reduction')
        plt.xlabel('Angles reduction')
        plt.title(name)
        plt.colorbar(orientation='vertical')


    

In [ ]:

    

for sinogram in log_progress(sinograms):
    x_resize = 2
    ang_resize = 2
    name = os.path.split(sinogram)[-1]
    print(name)
    pg, vg , sino_noise = load_data(sinogram, x_resize, ang_resize)
    
    # estimate noise
    D = 1.0 * (0.05 + (sino_noise/sino_noise.max())**2)
    Div = 1.0 / D
    
    plt.figure(figsize=(15,15))
    plt.imshow(sino_noise, cmap='gray')
    plt.colorbar(orientation='horizontal')
    plt.show()
    
    proj_id = astra.create_projector('cuda', pg, vg)
    W = astra.OpTomo(proj_id)
    x0 = np.zeros((sino_noise.shape[1], sino_noise.shape[1]))
    
    k = sino_noise.shape[0]/pg['DetectorWidth']**2/(np.pi/2)
    rec_fbp = alg.gpu_fbp(pg, vg, sino_noise)*k*x_resize * 2  # fixit
    
    eps = 1e-10
    n_iter = 200

    Lambda=4.

    #SIRT
    rec = RegVarSIRT.run(W, sino_noise, np.zeros_like(x0), eps=eps, n_it=n_iter, step='steepest')
    rec['fbp'] = rec_fbp
    
    plt.figure(figsize=(15,10))
    plt.subplot(121)
    plt.imshow(rec['rec'], interpolation='nearest')
    plt.title('SIRT')
    plt.subplot(122)
    plt.semilogy(rec['energy'])
    plt.grid()
    plt.show()
    
    with open('{}_{}_{}_s.pkl'.format(name, x_resize, ang_resize),'wb') as f:
        pickle.dump(rec, f)
    del rec['x_k']
    %xdel rec

    #SIRT+TV
    rec = RegVarSIRT.run(W, sino_noise, np.zeros_like(x0), Lambda=Lambda, eps=eps, n_it=n_iter, step='steepest')
    rec['fbp'] = rec_fbp
    
    plt.figure(figsize=(15,10))
    plt.subplot(121)
    plt.imshow(rec['rec'], interpolation='nearest')
    plt.title('SIRT+TV')
    plt.subplot(122)
    plt.semilogy(rec['energy'])
    plt.grid()
    plt.show()
    
    with open('{}_{}_{}_s_tv.pkl'.format(name, x_resize, ang_resize),'wb') as f:
        pickle.dump(rec, f)
    del rec['x_k']
    %xdel rec

    #CGLS
    rec = RegVarSIRT.run(W, sino_noise, np.zeros_like(x0), eps=eps, n_it=n_iter, step='CG')
    rec['fbp'] = rec_fbp
    
    plt.figure(figsize=(15,10))
    plt.subplot(121)
    plt.imshow(rec['rec'], interpolation='nearest')
    plt.title('CGLS')
    plt.subplot(122)
    plt.semilogy(rec['energy'])
    plt.grid()
    plt.show()
    
    with open('{}_{}_{}_cg.pkl'.format(name, x_resize, ang_resize),'wb') as f:
        pickle.dump(rec, f)
    del rec['x_k']
    %xdel rec

    #CGLS+TV
    rec = RegVarSIRT.run(W, sino_noise, np.zeros_like(x0), Lambda=Lambda, eps=eps, n_it=n_iter, step='CG')
    rec['fbp'] = rec_fbp
    
    plt.figure(figsize=(15,10))
    plt.subplot(121)
    plt.imshow(rec['rec'], interpolation='nearest')
    plt.title('CGLS+TV')
    plt.subplot(122)
    plt.semilogy(rec['energy'])
    plt.grid()
    plt.show()
    
    with open('{}_{}_{}_cg_tv.pkl'.format(name, x_resize, ang_resize),'wb') as f:
        pickle.dump(rec, f)
    del rec['x_k']
    %xdel rec 
    
    rec = RegVarSIRT.run(W, sino_noise, np.zeros_like(x0), eps=eps, n_it=3, step='CG')
    x0 = rec['rec']
    
    #VSIRT
    rec = RegVarSIRT.run(W, sino_noise, x0, Div, eps=eps, n_it=n_iter, step='steepest')
    rec['fbp'] = rec_fbp
    
    plt.figure(figsize=(15,10))
    plt.subplot(121)
    plt.imshow(rec['rec'], interpolation='nearest')
    plt.title('VSIRT')
    plt.subplot(122)
    plt.semilogy(rec['energy'])
    plt.grid()
    plt.show()
    
    with open('{}_{}_{}_vs.pkl'.format(name, x_resize, ang_resize),'wb') as f:
        pickle.dump(rec, f)
        
    del rec['x_k']
    %xdel rec

    #VSIRT+TV
    rec = RegVarSIRT.run(W, sino_noise, x0, Div, Lambda=Lambda, eps=eps, n_it=n_iter, step='steepest')
    rec['fbp'] = rec_fbp
    
    plt.figure(figsize=(15,10))
    plt.subplot(121)
    plt.imshow(rec['rec'], interpolation='nearest')
    plt.title('VSIRT+TV')
    plt.subplot(122)
    plt.semilogy(rec['energy'])
    plt.grid()
    plt.show()
    
    with open('{}_{}_{}_vs_tv.pkl'.format(name, x_resize, ang_resize),'wb') as f:
        pickle.dump(rec, f)
    
    del rec['x_k']
    %xdel rec

    #VCGLS
    rec = RegVarSIRT.run(W, sino_noise, x0, Div, eps=eps, n_it=n_iter, step='CG')
    rec['fbp'] = rec_fbp
    
    plt.figure(figsize=(15,10))
    plt.subplot(121)
    plt.imshow(rec['rec'], interpolation='nearest')
    plt.title('VCGLS')
    plt.subplot(122)
    plt.semilogy(rec['energy'])
    plt.grid()
    plt.show()
    
    with open('{}_{}_{}_vcg.pkl'.format(name, x_resize, ang_resize),'wb') as f:
        pickle.dump(rec, f)
    
    del rec['x_k']
    %xdel rec

    #VCGLS+TV
    rec = RegVarSIRT.run(W, sino_noise, x0, Div, Lambda=Lambda, eps=eps, n_it=n_iter, step='CG')
    rec['fbp'] = rec_fbp
    
    plt.figure(figsize=(15,10))
    plt.subplot(121)
    plt.imshow(rec['rec'], interpolation='nearest')
    plt.title('VCGLS+TV')
    plt.subplot(122)
    plt.semilogy(rec['energy'])
    plt.grid()
    plt.show()
    
    with open('{}_{}_{}_vcg_tv.pkl'.format(name, x_resize, ang_resize),'wb') as f:
        pickle.dump(rec, f)
    
    del rec['x_k']
    %xdel rec


    

In [ ]:

    

from cycler import cycler
for sinogram in log_progress(sinograms):
    plt.figure(figsize=(8,8))
    name = os.path.split(sinogram)[-1]
    short_name=name.split('_')[2]
    prefixes = ['s','cg']
    prefixes.extend(['v'+p for p in prefixes])
    prefixes.extend([p+'_tv' for p in prefixes])
#     prefixes.extend([p+'_n' for p in prefixes])
    
    plt.rc('axes', prop_cycle=(cycler('color', ['r', 'g', 'b', 'y']*4)+
                           cycler('linestyle', ['-']*4+ ['--']*4+['-']*4+ ['--']*4)+
                           cycler('lw', [1]*8+ [2]*8)))
    
    for prefix in log_progress(prefixes):
        if not os.path.exists('{}_2_2_{}.pkl'.format(name, prefix)):
            print('NOT found {}'.format(prefix))
            continue

            
        with open('{}_2_2_{}.pkl'.format(name, prefix),'rb') as f:
            res = pickle.load(f)
            
        plt.semilogy(res['energy'], label='{}_{}'.format(short_name, prefix))
        %xdel res
        
    plt.title('Energy')
    plt.xlabel('Interation number')
    plt.ylabel('a.u.')
    plt.grid()
    plt.legend(loc=0)
    plt.savefig('{}_energy.png'.format(short_name))
    plt.show()


    

In [ ]:

    

quality = {}
for sinogram in log_progress(sinograms):    
    mask = rec_reference>15
    name = os.path.split(sinogram)[-1]
    short_name=name.split('_')[2]
    quality[short_name] = {}
    
    prefixes = ['s','cg']
    prefixes.extend(['v'+p for p in prefixes])
    prefixes.extend([p+'_tv' for p in prefixes])
    prefixes.extend([p+'_n' for p in prefixes])
    
    for prefix in log_progress(prefixes):
        quality[short_name][prefix] = {}
        if not os.path.exists('{}_2_2_{}.pkl'.format(name, prefix)):
            print('NOT found {}'.format(prefix))
            continue
            
        with open('{}_2_2_{}.pkl'.format(name, prefix),'rb') as f:
            res = pickle.load(f)
            res['fbp'] = res['fbp']/4
    
        l1 = []
        l2 = []
        for i in range(len(res['x_k'])):
            l1.append(diff(res['x_k'][i], rec_reference, 1))
            l2.append(diff(res['x_k'][i], rec_reference, 2))
        
        quality[short_name][prefix]['l1'] = l1
        quality[short_name][prefix]['l2'] = l2
        
        with open('quality_all.pkl','wb') as f:
            pickle.dump(quality, f)
            
        pos_l1 = np.argmin(l1)
        pos_l2 = np.argmin(l2)
        
        plt.figure(figsize=(10,5))
        plt.plot(l1/np.max(l1), label='l1')
        # plt.hlines(sp.linalg.norm(res['fbp'].ravel()-rec_reference.ravel(), 1)/np.max(l1), 
        #            0, len(l1), 'r', label='fbp L1')
        plt.plot(l2/np.max(l2), label='l2')
        # plt.hlines(sp.linalg.norm(res['fbp'].ravel()-rec_reference.ravel(), 2)/np.max(l2), 
        #            0, len(l2), label='fbp L2')
        plt.legend(loc=0)
        plt.ylabel('a.u.')
        plt.xlabel('Iteration')
        plt.title('{}  {}'.format(short_name, prefix))
        plt.grid()
        plt.savefig('quality_{}_{}.png'.format(short_name, prefix))
        plt.show()
        
        plt.figure(figsize=(15,15))
        plt.subplot(221)
        plt.imshow(res['x_k'][pos_l2][600:800,400:600], interpolation='nearest', vmin=0, vmax=150)
        plt.title('{}_{}_{}'.format(short_name, prefix,pos_l2))
        plt.subplot(222)
        plt.imshow(res['x_k'][-1][600:800,400:600], interpolation='nearest', vmin=0, vmax=150)
        plt.title('{}_{}_{}'.format(short_name, prefix,len(l2)))
        plt.subplot(223)
        plt.imshow(res['fbp'][600:800,400:600], interpolation='nearest', vmin=0, vmax=150)
        plt.title('{}_{}'.format(short_name, 'fbp'))
        plt.subplot(224)
        # plt.imshow(res['x_k'][50][600:800,400:600], interpolation='nearest', vmin=0, vmax=300)
        plt.imshow(rec_reference[600:800,400:600], interpolation='nearest', vmin=0, vmax=150)
        plt.title('{}_{}'.format(short_name, 'reference'))
        
        plt.savefig(plt.savefig('rec_{}_{}.png'.format(short_name, prefix)))
        plt.show()
        
        
        
#     plt.title('Energy')
#     plt.xtitle('Interation number')
#     plt.ytitle('a.u.')
#     plt.grid()
#     plt.legend(loc=0)
#     plt.savefig('{}_energy.png'.format(sinogram))
#     plt.show()