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
construct a system matrix for projections of an image of size 50x50 with angles ranging from ð=0, 10, 20,..., 170°
Step1:
Step 2:
#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)
$\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)
Compute the forward projection of the Shepp-Logan image.
f = image.reshape( )
g =
#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)
intuition behind update formula:
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
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>
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
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)
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)
a) plot one of the ion trajectories (order = 'F')
b) plot the sum of all ion trajectories
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.