What's wrong with this code for tomographic reconstruction by the Fourier method?

Question

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()

Answer 1

OK I've cracked this finally.

The trick basically came down to putting some fftshift/ifftshifts in the right place so the 2D Fourier space representation was not wildly oscillatory and doomed to be impossible to interpolate accurately. At least that's what I think fixed it. Most of what limited understanding I do have of Fourier theory is based on the continuous integral formulation, and I always find the discrete domain and FFTs a bit... quirky.

While I find matlab code rather cryptic, I have to credit this implementation for at least giving me the confidence this reconstruction algorithm can be expressed reasonably compactly in this sort of an environment.

First I'll show results, then code:

Figure 1: a new, more complex target.

Figure 2: the sinogram (OK OK, it's the Radon transform) of the target.

Figure 3: the FFT-ed rows of the sinogram (plotted with DC at center).

Figure 4: the FFT-ed sinogram transformed to 2D FFT space (DC at center). Color is a function of absolute value.

Figure 4a: Zoom in on center of 2D FFT space just to show the radial nature of the sinogram data better.

Figure 5: Top row: the 2D FFT space interpolated from the radially arranged FFT-ed sinogram rows. Bottom row: the expected appearance from simply 2D FFT-ing the target.

Figure 5a: Zoom in on central region of the subplots in Fig5 to show these look to be in pretty good agreement qualitatively.

Figure 6: The acid test: inverse 2D FFT of the interpolated FFT space recovers the target. Lena's still looking pretty good despite everything we've just put her through (probably because there are enough sinogram "spokes" to cover the 2D FFT plane fairly densely; things get interesting if you reduce the number of exposure angles so this is no longer true).

Here's the code; brings up the plots in less than 15s on Debian/Wheezy's 64bit SciPy on an i7.

import math
import matplotlib
import matplotlib.pyplot as plt
import numpy as np

import scipy.interpolate
import scipy.fftpack
import scipy.misc
import scipy.ndimage.interpolation

S=256 # Size of target, and resolution of Fourier space
N=259 # Number of sinogram exposures (odd number avoids redundant direct opposites)

V=100 # Range on fft plots

# Convenience function
def sqr(x): return x*x

# Return the angle of the i-th (of 0-to-N-1) sinogram exposure in radians.
def angle(i): return (math.pi*i)/N

# Prepare a target image
x,y=np.meshgrid(np.arange(S)-S/2,np.arange(S)-S/2)
mask=(sqr(x)+sqr(y)<=sqr(S/2-10))
target=np.where(
    mask,
    scipy.misc.imresize(
        scipy.misc.lena(),
        (S,S),
        interp='cubic'
        ),
    np.zeros((S,S))
    )/255.0

plt.figure()
plt.title("Target")
plt.imshow(target)
plt.gray()

# Project the sinogram (ie calculate Radon transform)
sinogram=np.array([
        np.sum(
            scipy.ndimage.interpolation.rotate(
                target,
                np.rad2deg(angle(i)), # NB rotate takes degrees argument
                order=3,
                reshape=False,
                mode='constant',
                cval=0.0
                )
            ,axis=0
            ) for i in xrange(N)
        ])

plt.figure()
plt.title("Sinogram")
plt.imshow(sinogram)
plt.jet()

# Fourier transform the rows of the sinogram, move the DC component to the row's centre
sinogram_fft_rows=scipy.fftpack.fftshift(
    scipy.fftpack.fft(
        scipy.fftpack.ifftshift(
            sinogram,
            axes=1
            )
        ),
    axes=1
    )

plt.figure()
plt.subplot(121)
plt.title("Sinogram rows FFT (real)")
plt.imshow(np.real(sinogram_fft_rows),vmin=-V,vmax=V)
plt.subplot(122)
plt.title("Sinogram rows FFT (imag)")
plt.imshow(np.imag(sinogram_fft_rows),vmin=-V,vmax=V)

# Coordinates of sinogram FFT-ed rows' samples in 2D FFT space
a=np.array([angle(i) for i in xrange(N)])
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()

plt.figure()
plt.title("Sinogram samples in 2D FFT (abs)")
plt.scatter(
    srcx,
    srcy,
    c=np.absolute(sinogram_fft_rows.flatten()),
    marker='.',
    edgecolor='none',
    vmin=-V,
    vmax=V
    )

# Let the central slice theorem work its magic!
# Interpolate the 2D Fourier space grid from the transformed sinogram rows
fft2=scipy.interpolate.griddata(
    (srcy,srcx),
    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(np.real(fft2),vmin=-V,vmax=V)
plt.subplot(222)
plt.title("Recon (imag)")
plt.imshow(np.imag(fft2),vmin=-V,vmax=V)

# Show 2D FFT of target, just for comparison
expected_fft2=scipy.fftpack.fftshift(
    scipy.fftpack.fft2(
        scipy.fftpack.ifftshift(
            target
            )
        )
    )

plt.subplot(223)
plt.title("Expected (real)")
plt.imshow(np.real(expected_fft2),vmin=-V,vmax=V)
plt.subplot(224)
plt.title("Expected (imag)")
plt.imshow(np.imag(expected_fft2),vmin=-V,vmax=V)

# Transform from 2D Fourier space back to a reconstruction of the target
recon=np.real(
    scipy.fftpack.fftshift(
        scipy.fftpack.ifft2(
            scipy.fftpack.ifftshift(fft2)
            )
        )
    )

plt.figure()
plt.title("Reconstruction")
plt.imshow(recon,vmin=0.0,vmax=1.0)
plt.gray()

plt.show()

Update 2013-02-17: If you've been interested enough to wade through that lot, some more output from the self-study programme of which it was a part can be found in the form of this poster. The body of code in this repository may also be of interest (although note the code isn't nearly as streamlined as that above). I may try and repackage it as an IPython "notebook" at some point.

Answer 2

I don't know exactly where the problem is, but the slice theorem means that these two special cases should be true:

fft2(target)[0] = fft(sinogram[270])
fft2(target)[:,0] = fft(sinogram[0])

So follow through your code and try to find the point where these stop being equivalent, working forward from the sinogram and backwards from the generated 2D FFT.

This doesn't look right:

In [47]: angle(expected_fft2[127:130,127:130])
Out[47]: 
array([[-0.07101021,  3.11754929,  0.02299738],
       [ 3.09818784,  0.        , -3.09818784],
       [-0.02299738, -3.11754929,  0.07101021]])

In [48]: fft2_ = fft2_real+1.0j*fft2_imag

In [49]: angle(fft2_[127:130,127:130])
Out[49]: 
array([[ 3.13164353, -3.11056554,  3.11906449],
       [ 3.11754929,  0.        , -3.11754929],
       [ 3.11519503,  3.11056604, -2.61816765]])

The 2D FFT you're generating is rotated 90 degrees from what it should be?

I'd suggest working with magnitude and phase rather than real and imaginary, so you can see more easily what's happening:

(The white corners are -inf from doing log(abs(0)), they're not a problem)

Answer 3

I believe that the actual theoretical reason why the first solution did not work comes from the fact that the rotations are done with respect to the centers of the images, inducing an offset of [S/2, S/2], which means that each of the rows of your sinogram is not from 0 to S, but rather from -S/2 to S/2. In your example, the offset is actually offset = np.floor(S/2.). Note that this works for S even or odd, and is equivalent to what you did in your code S/2 (although being more explicit avoids issues, when S is a float, for instance).

My guess is that the phase delays this shift introduces in the Fourier transform (FT) are at the origin of what you talk about in your second message: the phases are messed up, and one needs to compensate that shift in order to be able to apply the inversion of the Radon transform. One needs to dig more into that theory in order to be sure of what is exactly needed for the inverse to work as expected.

To compensate that offset, you can either use fftshift as you did (which puts the center of each row at the beginning, and since using DFT actually corresponds to computing the Fourier Transform of a S-periodic signal, you end up with the right stuff), or explicitly compensate this effect in the complex Fourier transform, when computing the sinogram FT. In practice, instead of:

sinogram_fft_rows=scipy.fftpack.fftshift(
    scipy.fftpack.fft(
        scipy.fftpack.ifftshift(
            sinogram,
            axes=1
            )
        ),
    axes=1
    )

you can remove the ifftshift and multiply each row by a corrective vector:

offset = np.floor(S/2.)
sinogram_fft_rows = scipy.fftpack.fftshift(
    scipy.fftpack.fft(sinogram, axis=1)
    * (np.exp(1j * 2.* np.pi * np.arange(S) * offset / S)),
    axes=1)

This comes from the Fourier transform properties, when considering a time-shift (check the FT wikipedia page for the "shift theorem", and apply for the shift equal to - offset - because we put the image back around the center).

Likewise, you can apply the same strategy to the reconstruction, and replace the fftshift by the correction of the phases, in both dimensions, but in the other direction (compensating back):

recon=np.real(
    scipy.fftpack.ifft2(
        scipy.fftpack.ifftshift(fft2)
        *  np.outer(np.exp(- 1j * 2.* np.pi * np.arange(S) * offset / S),
                    np.exp(- 1j * 2.* np.pi * np.arange(S) * offset / S))
        )
    )

Well, this does not improve your solution, but rather sheds another light on the theoretical aspects of your question. Hope that helps!

Additionally, I am not so fond of using fftshift because it tends to mess around with the way the fft is computed. In this case, however, you need to put the center of the FT in the center of the image before the interpolation to get fft2 (or at least be careful when setting r - so you could make it completely fftshift-free!), and fftshift really comes handy there. I however prefer to keep the use of that function for visualisation purposes, and not within the computation "core". :-)

Best regards,

Jean-Louis

P.S.: have you tried to reconstruct the image without cropping the circle around? that gives a pretty cool blurring effect on the corners, would be nice to have such a feature in programs like Instagram, wouldn't it?