assignments.ipynb

Kernel: Python 3 (Anaconda)

Exercise 1

import numpy as np
import matplotlib.pyplot as plt

from skimage.io import imread
from skimage import data_dir
from skimage.transform import radon, rescale, iradon

# read phantom image
image = imread(data_dir + "/phantom.png", as_grey=True)

# create sinogram
theta = np.linspace(0., 180., max(image.shape), endpoint=False)
sinogram = radon(image, theta=theta, circle=True)

# plot
plt.subplot(121)
plt.imshow(image)
plt.title('phantom')

plt.subplot(122)
plt.imshow(sinogram)
plt.title('sinogram')

plt.show()

# reconstruction using FBP with 'ramp' filter. Other filters are: shepp-logan, cosine, hamming, hann or none
reconstruction_fbp = iradon(sinogram, theta=theta, circle=True, filter = 'ramp')
error = reconstruction_fbp - image
print('FBP rms reconstruction error: %.3g' % np.sqrt(np.mean(error**2)))

# plot
plt.subplot(121)
plt.imshow(reconstruction_fbp)
plt.title('FBP reconstruction')

plt.subplot(122)
plt.imshow(error)
plt.title('error')

plt.show()

FBP rms reconstruction error: 0.0308

Exercise 2

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt

# define matrix
K = np.array([[1,1,1,0,0,0,0,0,0],[0,0,0,1,1,1,0,0,0],[0,0,0,0,0,0,1,1,1],[1,0,0,1,0,0,1,0,0],[0,1,0,0,1,0,0,1,0],[0,0,1,0,0,1,0,0,1],[1,0,0,0,1,0,0,0,1],[0,1,0,0,0,1,0,0,0],[0,0,0,1,0,0,0,1,0]])

# true image
utrue = np.array([1.0,1.0,1.0,3.0,3.0,3.0,1.0,1.0,3.0])

# data with noise
sigma = 0.1
noise = np.random.random(9)
f = K@xtrue + sigma*noise

## least squares solution
uls = np.linalg.lstsq(K,f)[0]

print('reconstruction error :', np.linalg.norm(uls - utrue)/np.linalg.norm(utrue))

## constrained least squares
ulsc = opt.lsq_linear(K,f,bounds=(0,3)).x

print('reconstruction error :', np.linalg.norm(ulsc - utrue)/np.linalg.norm(utrue))

reconstruction error : 0.129847250492
reconstruction error : 0.0165785799813

## ART iterative reconstruction

# define projection of vector u \in R^9 on [1,3]^9
Pc = lambda u : np.minimum(np.maximum(u,1),3)

uart = np.zeros(9)

for k in range(50):
    i = k % 9
    uart = Pc(uart + ((f[i] - K[i,:]@uart)/(K[i,:]@K[i,:]))*K[i,:])
    
print('reconstruction error :', np.linalg.norm(uart - utrue)/np.linalg.norm(utrue))

reconstruction error : 0.0186726720804

## column action method: access the jth column as K[:,j]