Tutorial 2: System matrix

22/11/2022

email: ines.butz@physik.uni-muenchen.de

Getting started

for instructions on how to install python and jupyter please refer to slides from Tutorial 1

you can find the code for this tutorial here: https://gitlab.physik.uni-muenchen.de/Ines.Butz/tutorials-inverse-problems-and-machine-learning-in-medical-physics

if you already cloned the repository: in the command line, navigate to your repository folder, then "git pull"

alternatively, download the file Tutorial2- System matrix.ipynb directly and open it with jupyter.

additional packages to be installed in your environment:

• tqdm (this is just for making progress bars, you can also run the code without if you remove the tqdm calls.)

"conda install -c conda-forge tqdm"

Constructing a system matrix

import numpy as np
from skimage.data import shepp_logan_phantom
from skimage.transform import downscale_local_mean
import matplotlib.pyplot as plt
from scipy.ndimage import rotate
from tqdm.auto import tqdm

let's start with $\rho = 1$ and $\vartheta = 0$:

• create an image with first column = 1, rest = 0
• reshape the image: write all columns one after the other, in one single long column

im = np.zeros((4,4))
im[:,0] = 1

plt.imshow(im)

<matplotlib.image.AxesImage at 0x7f128dbfd990>

how to reshape a numpy array: https://numpy.org/doc/stable/reference/generated/numpy.reshape.html

im_flat = im.reshape((16,1), order = 'F')

plt.imshow(im_flat)

<matplotlib.image.AxesImage at 0x7f128dab6fd0>

let's do the same for $\rho = 1$ and $\vartheta = 45°$:

im = np.zeros((4,4))
im_rot = rotate(im, -45, reshape = False)
im_rot[:,0] = 1
im = rotate(im_rot, 45, reshape = False)

plt.imshow(im)

<matplotlib.image.AxesImage at 0x7f128dc67f10>

im_flat = im.reshape((16,1), order = 'F')

plt.imshow(im_flat)

<matplotlib.image.AxesImage at 0x7f13024f1b10>

Exercise 1:

construct the system matrix for $\rho = 1, ..., 4$ and $\vartheta = 0, 45 °$

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

#generate an empty system matrix, where you later store every computed column
n_pixels =
n_rays =
sys_matrix_T =

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

sys_matrix_column = 0
for theta in np.array([0,45]):
    for rho in np.arange(0,4):
        sys_matrix_column += 1
plt.figure()
plt.imshow(sys_matrix_T)

Exercise 2:

construct a system matrix for projections of an image of size 100x100 with angles ranging from $\vartheta = 0,1,...179°$

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

n_pixels =
n_rays =
sys_matrix_T =

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

sys_matrix_column = 0
for theta in tqdm(np.arange(0,180)):
    for rho in np.arange(0,100):
        sys_matrix_column += 1
plt.figure()
plt.imshow(sys_matrix_T)

Calculating a forward projection

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

(100, 100)

A = sys_matrix_T.T
f = image.reshape((A.shape[1], 1), order = 'F')
g = np.dot(A,f)
print(g.shape)

recap from last tutorial: compute projections by integrating

projection = np.sum(image, axis = 1)
fig, axs = plt.subplots(nrows = 1, ncols = 2, sharey = True, figsize = (20,10))
axs[0].imshow(image, cmap = 'gray')
axs[1].plot(projection, np.arange(np.shape(image)[0]))
#this just makes sure that the plots have the same height
axs[1].set_ylim(99,0)
axs[0].set_ylabel('image rows')
axs[0].set_xlabel('image columns')
axs[1].set_ylabel('image rows')
axs[1].set_xlabel('projection value')
asp = np.abs(np.diff(axs[1].get_xlim())[0]) / np.abs(np.diff(axs[1].get_ylim())[0])
axs[1].set_aspect(asp)

Exercise 3

visualize the same projection (computed using system matrix):

• we changed definition of $\vartheta = 0°$ $\rightarrow$ get projection from $\vartheta = 90°$ $
• which columns of the system matrix does that correspond to?
• caution: ordering of the columns

projection = g[90*100:91*100][::-1]

fig, axs = plt.subplots(nrows = 1, ncols = 2, figsize = (20,10))
axs[0].imshow(image, cmap = 'gray')
axs[1].plot(projection, np.arange(np.shape(image)[0]))
#this just makes sure that the plots have the same height
axs[1].set_ylim(99,0)
axs[0].set_ylabel('image rows')
axs[0].set_xlabel('image columns')
axs[1].set_ylabel('image rows')
axs[1].set_xlabel('projection value')
asp = np.abs(np.diff(axs[1].get_xlim())[0]) / np.abs(np.diff(axs[1].get_ylim())[0])
axs[1].set_aspect(asp)

Calculating a backprojection

Recap from analytical reconstruction:

• cannot directly invert projection operation
• in continous/ adequately sampled scenario:
  • can take 'detour' via Fourier Slice Theorem to reconstruct original image
  • alternatively: equivalent to filtered backprojection

Filtered Backprojection:
$f(x,y) = \int_{\vartheta=0}^{\pi}{g(\rho, \vartheta)\ \ast k_{ramp}(\rho)|_{\rho = xcos(\vartheta)+ysin(\vartheta)}}$ $\rightarrow$ discrete: $f(x,y) = \sum_{\vartheta=0}^{\pi}{a_{x,y}(\vartheta)g(\rho, \vartheta)\ \ast k_{ramp}(\rho)|_{\rho = xcos(\vartheta)+ysin(\vartheta)}}$

Backprojection:
$f(x,y) = \int_{\vartheta=0}^{\pi}{g(\rho, \vartheta)|_{\rho = xcos(\vartheta)+ysin(\vartheta)}}$ $\rightarrow$ discrete: $f(x,y) = \sum_{\vartheta=0}^{\pi}{a_{x,y}(\vartheta)g(\rho, \vartheta)|_{\rho = xcos(\vartheta)+ysin(\vartheta)}}$

we can see the same in system matrix formalism:

• $g = Af \rightarrow f = A^{-1}g$, but cannot invert system matrix ($A$ not invertible!)
• instead, we can compute 'backprojection' as $\tilde{f} = A^T g$

compare to discrete backprojection that we already knew:
$f(x,y) = \sum_{\vartheta=0}^{\pi}{a_{x,y}(\vartheta)g(\rho, \vartheta)|_{\rho = xcos(\vartheta)+ysin(\vartheta)}}$
$x,y\ \hat{=} 1,..., J$ (remember reshape operation)
$\vartheta \hat{=} 1,...,I$
$\rightarrow$ same as $f = A^T g$!

Exercise 4:

compute the backprojection of the projected Shepp-Logan data computed above $\rightarrow$ do you understand why it is blurred?

f_tilde =
f_tilde = f_tilde.reshape( , , order = 'F')

plt.imshow(f_tilde, cmap = 'gray')

Bonus: Algebraic Reconstruction Technique (ART)

Excercise 5:

Using the system matrix $A = \{a_{i,j}\}$, $i = 1,..,I; j = 1,...,J$, and the projections $g$ of the Shepp-Logan Phantom computed above, implement the ART algorithm for 100 iterations. Initilialize with a zero image.

f0 = np.zeros(np.dot(A.T, g).shape)
fn = f0.squeeze()
for k in tqdm(np.arange(   )):
    for i in np.arange(    ):
image_recon = fn.reshape(  ,  , order = 'F')
plt.figure()
plt.imshow(image_recon, cmap = 'gray')