I've been playing around with tomographic reconstruction algorithms recently. I already have nice working implementations of FBP, ART, a SIRT/SART-like iterative scheme and even using straight linear algebra (slow!). This question is not about any of those techniques; answers of the form "why would anyone do it that way, here's some FBP code instead" are not what I'm looking for.
The next thing I wanted to do with this programme was "complete the set" and implement the so-called "Fourier reconstruction method". My understanding of this is basically that you apply a 1D FFT to the sinogram "exposures", arrange those as radial "spokes of a wheel" in 2D Fourier space (that this is a useful thing to do follows directly from the central slice theorem), interpolate from those points to a regular grid in that 2D space, and then it should be possible to inverse Fourier-transform to recover the original scan target.
Sounds simple, but I haven't had much luck getting any reconstructions which look anything like the original target.
The Python (numpy/SciPy/Matplotlib) code below is about the most concise expression I could come up with of what I'm trying to do. When run, it displays the following:
Figure 1: the target
Figure 2: a sinogram of the target
Figure 3: the FFT-ed sinogram rows
Figure 4: the top row is the 2D FFT space interpolated from the Fourier-domain sinogram rows; the bottom row is (for comparison purposes) the direct 2D FFT of the target. This is the point at which I'm starting to get suspicious; the plots interpolated from the sinogram FFTs look similar to the plots made by directly 2D-FFTing the target... and yet different.
Figure 5: the inverse-Fourier transform of Figure 4. I'd have hoped this would be a bit more recognizable as the target than it actually is.
Any ideas what I'm doing wrong ? Not sure if my understanding of Fourier method reconstruction is fundamentally flawed, or there's just some bug in my code.
import math
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate
import scipy.fftpack
import scipy.ndimage.interpolation
S=256 # Size of target, and resolution of Fourier space
A=359 # Number of sinogram exposures
# Construct a simple test target
target=np.zeros((S,S))
target[S/3:2*S/3,S/3:2*S/3]=0.5
target[120:136,100:116]=1.0
plt.figure()
plt.title("Target")
plt.imshow(target)
# Project the sinogram
sinogram=np.array([
np.sum(
scipy.ndimage.interpolation.rotate(
target,a,order=1,reshape=False,mode='constant',cval=0.0
)
,axis=1
) for a in xrange(A)
])
plt.figure()
plt.title("Sinogram")
plt.imshow(sinogram)
# Fourier transform the rows of the sinogram
sinogram_fft_rows=scipy.fftpack.fftshift(
scipy.fftpack.fft(sinogram),
axes=1
)
plt.figure()
plt.subplot(121)
plt.title("Sinogram rows FFT (real)")
plt.imshow(np.real(np.real(sinogram_fft_rows)),vmin=-50,vmax=50)
plt.subplot(122)
plt.title("Sinogram rows FFT (imag)")
plt.imshow(np.real(np.imag(sinogram_fft_rows)),vmin=-50,vmax=50)
# Coordinates of sinogram FFT-ed rows' samples in 2D FFT space
a=(2.0*math.pi/A)*np.arange(A)
r=np.arange(S)-S/2
r,a=np.meshgrid(r,a)
r=r.flatten()
a=a.flatten()
srcx=(S/2)+r*np.cos(a)
srcy=(S/2)+r*np.sin(a)
# Coordinates of regular grid in 2D FFT space
dstx,dsty=np.meshgrid(np.arange(S),np.arange(S))
dstx=dstx.flatten()
dsty=dsty.flatten()
# Let the central slice theorem work its magic!
# Interpolate the 2D Fourier space grid from the transformed sinogram rows
fft2_real=scipy.interpolate.griddata(
(srcy,srcx),
np.real(sinogram_fft_rows).flatten(),
(dsty,dstx),
method='cubic',
fill_value=0.0
).reshape((S,S))
fft2_imag=scipy.interpolate.griddata(
(srcy,srcx),
np.imag(sinogram_fft_rows).flatten(),
(dsty,dstx),
method='cubic',
fill_value=0.0
).reshape((S,S))
plt.figure()
plt.suptitle("FFT2 space")
plt.subplot(221)
plt.title("Recon (real)")
plt.imshow(fft2_real,vmin=-10,vmax=10)
plt.subplot(222)
plt.title("Recon (imag)")
plt.imshow(fft2_imag,vmin=-10,vmax=10)
# Show 2D FFT of target, just for comparison
expected_fft2=scipy.fftpack.fftshift(scipy.fftpack.fft2(target))
plt.subplot(223)
plt.title("Expected (real)")
plt.imshow(np.real(expected_fft2),vmin=-10,vmax=10)
plt.subplot(224)
plt.title("Expected (imag)")
plt.imshow(np.imag(expected_fft2),vmin=-10,vmax=10)
# Transform from 2D Fourier space back to a reconstruction of the target
fft2=scipy.fftpack.ifftshift(fft2_real+1.0j*fft2_imag)
recon=np.real(scipy.fftpack.ifft2(fft2))
plt.figure()
plt.title("Reconstruction")
plt.imshow(recon,vmin=0.0,vmax=1.0)
plt.show()