

In [1]:

    
%matplotlib inline

"""
Created on Wed Nov 19 00:30:14 2014

@author: talbot
"""

from __future__ import print_function

# projection operator
#

import numpy as np
from projections import build_projection_operator
from Imgdeblur import makespdiag, deblur_tikh_sparse, solve_tykhonov_sparse
from time import time
import matplotlib.pyplot as plt

import scipy as sp
import scipy.sparse
from scipy.misc import lena

import skimage
import skimage.transform

def project_operator(img):
    l = img.shape[0]
    H = build_projection_operator(l)
    return H
        
## Reconstruction with well-suited algorithm
def gen_tikhonov_reconstruct(sino, mylambda, H):
    '''
        This algorithm performs a general-purpose reconstruction
        with a damped least-square algorithm that does not require
        the computation of H'H.
        It is not suitable to the case where H is symmetric. In this
        case it is better to construct H'H and solve with a congugate
        gradient.
    '''
    t1 = time()
    tikh_sol = scipy.sparse.linalg.lsmr(H,sino,mylambda)
    t2 = time()
    print("Time Tikhonov = %f\n" % (t2-t1))
    return tikh_sol
   
def sym_tikhonov_reconstruct(sino, mylambda, H, prior='Lap'):
    ''' 
        Solves Tikhonov problems with a sparse symetric H.
        if H is non-symmetric but sparse, it is likely that
        H'H will not be sparse and this algorithm is inefficient.
        Acceptable priors are 'Lap', 'Grad', 'Id'
    '''
    t1 = time()
    ## we use an identity prior for comparison with gen_
    ## tikhonov_reconstruct because for this one it is the 
    ## only available prior.
    sol_tikh = deblur_tikh_sparse(sino, H, mylambda, 'Id')
    t2 = time()
    print("Time = %f\n" % (t2-t1))
    return sol_tikh
    
if __name__ == '__main__':
    try:
        import mkl
        have_mkl = True
        backend = 'anaconda+mkl'
        print("Running with MKL Acceleration with %d threads" % mkl.get_max_threads())
    except ImportError:
        have_mkl = False
        backend = 'anaconda'
        print("Running with normal backends")
    
    
    fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 5))    
    imglena = np.uint8(lena()) # 512x512 too big
    
    lenatiny = skimage.transform.rescale(imglena,0.5,3) # small image

    ax1.imshow(lenatiny, interpolation='nearest', cmap=plt.cm.gray)
    ax1.axis('off')

    Htiny = project_operator(lenatiny)
    sinotiny = Htiny.dot(lenatiny.ravel())[:,None]
    ## add noise
    sinotnoise = sinotiny + 0.2 * np.random.randn(*sinotiny.shape)

    ax2.imshow(sinotnoise.reshape(lenatiny.shape),interpolation='nearest',cmap="gray")
    ax2.axis('off')
    ## Solve with symmetric solver, here a bad idea
    ## sol_tikh = sym_tikhonov_reconstruct(sinotnoise.reshape(lenatiny.shape), \
    ##         0.2, Htiny)

    ## Solve with a more efficient solver
    ## we do not reshape here because we only use the identity prior
    sol_tikh = gen_tikhonov_reconstruct(sinotnoise, 10, Htiny)


    img_tikh=sol_tikh[0].reshape(lenatiny.shape)

    ax3.imshow(img_tikh, interpolation='nearest', cmap=plt.cm.gray)
    ax3.axis('off')









    



Running with MKL Acceleration with 2 threads
Time Tikhonov = 33.601303

Computed Tomography with Tikhonov Regularization