Tutorial 3: Algebraic Reconstruction

import numpy as np
from tqdm.auto import tqdm
from scipy.ndimage import rotate
import matplotlib.pyplot as plt
from skimage.data import shepp_logan_phantom
from skimage.transform import downscale_local_mean
import time
from IPython.display import display, clear_output

Part 1: Algebraic Reconstruction

Recap: Creating a system matrix

Exercise (Recap):

construct a system matrix for projections of an image of size 50x50 with angles ranging from 𝜗=0, 10, 20,..., 170°
Step1:

• generate an empty system matrix, where you later store every computed column
• what is the shape of the system matrix? $\rightarrow$ $n_{rows} =$ #pixels, $n_{columns} =$ #rays

Step 2:

• loop over all angles $\vartheta$
  • loop over all trajectories $\rho$
    • calculate the pixel weights for each trajectory
    • flatten the image to a single column
    • write the column into the system matrix
• transpose the matrix

#Create a system matrix
n_pixels =
n_rays =
sys_matrix_T = np.zeros((n_pixels, n_rays))
sys_matrix_column = 0
for theta in tqdm(  ):
    for rho in np.arange(  ):
        sys_matrix_column += 1
A =

  0%|          | 0/18 [00:00<?, ?it/s]

plt.figure(figsize = (15,10))
plt.imshow(sys_matrix_T)

Recap: Forward- and Back-Projection

$\rightarrow$ number of projections (rays/ ions/ ...): $I$
$\rightarrow$ number of pixels: $J$

forward projection:

back projection:

image = shepp_logan_phantom()
image = downscale_local_mean(image, (8, 8))
print(image.shape)

(50, 50)

Exercise:

Compute the forward projection of the Shepp-Logan image.

f = image.reshape(  )
g =

Excursion: Updating plots within a loop

#Images 
nrows = 1
ncols = 1
fig, axs = plt.subplots(nrows, ncols, squeeze = False)
for k in np.arange(10):
    axs[0,0].clear()
    axs[0,0].imshow(np.random.rand(50,50))
    display(fig)
    clear_output(wait = True)
    plt.pause(0.1)

#Graphs
nrows = 1
ncols = 1
fig, axs = plt.subplots(nrows, ncols, squeeze = False)
#create empty lists
x = []
y = []
for k in np.arange(20):
    #write value into list
    x += [k]
    y += [4*k**2+3]
    axs[0,0].clear()
    axs[0,0].plot(x, y)
    display(fig)
    clear_output(wait = True)
    plt.pause(0.1)

Algebraic Reconstruction

intuition behind update formula:

• loop over number of iterations
  • update is calculated for each projection $i$:
    • compute forward projection $\tilde{g_i} =\vec{f^{(n)}}\cdot \vec{a_i}$
    • compare to measured projection $\epsilon_i = \tilde{g_i}-g_i$
    • compute correction: backproject error term $\epsilon_i$ ($\epsilon_i \cdot \vec{a_i}$), normalize by trajectory length (divide by $\vec{a_i} \cdot \vec{a_i}$)

f0 = np.zeros(np.dot(A.T, g).shape)
fn = f0.squeeze()
a_0 = A[0,:]
g0 = g[0]
print(np.shape(fn), 'image')
print(np.dot(fn, a_0), 'forward projection')
print((np.dot(fn, a_0)-g0)/(np.dot(a_0, a_0)), 'error term')
print(np.shape((np.dot(fn, a_0)-g0)/(np.dot(a_0, a_0))*a_0), 'back-projection of error term')

(2500,) image
0.0 forward projection
[6.53760315e-22] error term
(2500,) back-projection of error term

Exercise (Recap):

Implement the ART algortihm to reconstruct the 50x50 pixel Shepp-Logan image using the system matrix $A$ and projection data $g$ computed above. Initialize with a zero image and run for 80 iterations.

f0 =
fn = f0.squeeze()
for k in tqdm(  ):
    for i in np.arange(  ):

  0%|          | 0/80 [00:00<?, ?it/s]

plt.figure()
plt.imshow(fn.reshape(50,50, order = 'F'))

<matplotlib.image.AxesImage at 0x2277cf312a0>

Exercise:

Let's take a closer look...
For one iteration:
Plot the current reconstructed image, update the plot with every image update.

f0 = np.zeros(np.dot(A.T, g).shape)
fn = f0.squeeze()
fig, axs = plt.subplots(1,1, squeeze = False)
for i in np.arange(A.shape[0]):
    #here: compute fn
    axs[0,0].clear()
    axs[0,0].set_title('update '+ str(i)+' /'+str(A.shape[0]))
    #here: plot image, live -update

Exercise:

For each iteration:
Compute the error of the current reconstructed image with respect to the ground truth image (for example np.linalg.norm(current image - ground truth)). Plot the error as a function of iteration number (live plot).

f0 = np.zeros(np.dot(A.T, g).shape)
fn = f0.squeeze()
fig, axs = plt.subplots(1,1,squeeze = False)
axs[0,0].set_xlabel('iteration')
axs[0,0].set_ylabel('error wrt groud truth')
iterations = []
error = []
for k in   :
    for i in    :
        #here: compute fn
    iterations += [  ]
    error += [  ]
    axs[0,0].clear()
    axs[0,0].set_title('iteration '+ str(k+1)+' / 80')
    # here: plot error (live-plot)

Part 2: System matrix in ion imaging

from scipy.io import loadmat
import scipy

sl_rsp = loadmat('ict_Ines.mat')['ict']

plt.figure()
plt.imshow(sl_rsp, cmap = 'gray')
print(np.shape(sl_rsp))

(128, 128)

sys_matrix_T = scipy.sparse.load_npz('matrice2d_Angle1.npz').toarray()
print(sys_matrix_T.shape)

(16384, 9000)

Exercise:

a) plot one of the ion trajectories (order = 'F')
b) plot the sum of all ion trajectories

Bonus Exercise

Read the system matrix for each angle and construct the complete system matrix (for all angles).
Compute the sinogram for the RSP Shepp-Logan phantom.
Using the system matrix, and the computed sinogram, reconstruct an Ion CT using the Algebraic Reconstruction Technique.