The radon transformation was taken directly from https://gist.github.com/fubel/ad01878c5a08a57be9b8b80605ad1247, but with tweaks to allow variable steps in transform.
The back-projection is a crude reverse of the forward transformation. This is an unfiltered back projection.
- skimage.transform.rotate instead of scipy.misc.imrotate (was removed from SciPy)
- reference for standard test image: https://scikit-image.org/docs/dev/auto_examples/transform/plot_radon_transform.html. This also has a scikit-image implementation of a better radon and inverse-radon transformation.
Code:
import numpy as np
import matplotlib.pyplot as plt
from skimage.transform import rotate
def discrete_radon_transform(image, steps):
R = np.zeros((image.shape[0],steps), dtype='float')
for s in range(steps):
rotation = rotate(image, -s*180/steps).astype('float')
R[:,s] = sum(rotation)
return R
def discrete_inv_radon_transform(R):
steps = R.shape[1]
image = np.zeros((R.shape[0],R.shape[0]), dtype='float')
update = np.zeros((R.shape[0],R.shape[0]), dtype='float')
for s in range(steps):
update[:,:] = R[:,s]
rotation = rotate(update, s*180/steps).astype('float')
image = image + rotation
return image
# generate simple image
image =np.zeros((400,400))
image[40:50,100:110]=255
# standard test image
#from skimage.data import shepp_logan_phantom
#image = shepp_logan_phantom()
steps = 500
aspect = steps/image.shape[0] # to keep sinogram plot square
radon = discrete_radon_transform(image, steps)
inv_radon = discrete_inv_radon_transform(radon)
# Plot the original image, sinogram, and backprojected sinogram
plt.figure(figsize=(20,20))
plt.subplot(1, 3, 1), plt.imshow(image, cmap='gray')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 3), plt.imshow(inv_radon, cmap='gray')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 2), plt.imshow(radon, cmap='gray',aspect = aspect)
plt.xticks([]), plt.yticks([])
Output (original image, sinogram, back-projection):
Now using a standard test image (uncomment lines 27 and 28):
Odds and Sods 