22/11/2022
email: ines.butz@physik.uni-muenchen.de
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:
"conda install -c conda-forge tqdm"
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$:
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>
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
n_pixels =
n_rays =
sys_matrix_T =
Step 2:
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)
construct a system matrix for projections of an image of size 100x100 with angles ranging from $\vartheta = 0,1,...179°$
Step1:
n_pixels =
n_rays =
sys_matrix_T =
Step 2:
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)
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)
visualize the same projection (computed using system matrix):
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)
Recap from analytical reconstruction:
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:
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$!
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')
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')