python image-processing fft interpolation complex-numbers

Interpolate matrix of complex numbers

I have want to rotate a row of complex numbers (which actually is a 1D FFT of a line of the Radon transform), I use imrotate in Matlab but I don't think the interpolation is doing what it should.

The goal is to reproduce the conversion from Radon to image space with the Projection-slice theorem.

^{(Image from Wikipedia)}

I need to take each row of the Radon transform and rotate it according to its angle and put it at the corresponding angle in a 2D matrix. Once this is done I can get the 2D ifft2 to recover the image (in theory). This is the goal. Anyone can help?

I thought using imrotate, but maybe that isn't the right thing? The goal is to map the FFT'd rows of the Radon transform to their correct position in a circle as shown in the figure above.

This is the actual result with rotate and nearest neighbour interpolation. The result on the right should be the usual SheppLogan phantom.

import numpy as np
import matplotlib.pyplot as plt
from skimage.io import imread
from skimage.data import shepp_logan_phantom
from skimage.transform import radon, rescale
from skimage.transform import iradon
from skimage.transform import rotate
import cv2 


x=shepp_logan_phantom()
x=cv2.resize(x, (128,128), interpolation = cv2.INTER_AREA)
theta=np.linspace(0,180,len(x))
R=radon(x,theta)

temp_=np.zeros((128,128)).astype(np.complex128)
fullFft2D=np.zeros((128,128)).astype(np.complex128)

for i in range(len(theta)):
    temp_[63,:]=np.fft.fftshift(np.fft.fft(R[:,i])).T
    fft_real=rotate(np.real(temp_),theta[i],order=0)
    fft_imag=rotate(np.imag(temp_),theta[i],order=0)
    fullFft2D += fft_real+1j*fft_real
    temp_=np.zeros((128,128)).astype(np.complex128)


plt.imshow(np.fft.fftshift(np.abs(np.fft.ifft2(np.fft.ifftshift(fullFft2D)))))

I have implemented what you (@Luengo) said as:

res=np.zeros((128,128))
tmp_=np.zeros((128,128)).astype(np.complex128)
for i in range(len(theta)):
    kspace_row = np.fft.fftshift(np.fft.fft(R[:,i])).T
    tmp_[63,:] = kspace_row
    res +=  rotate(np.abs(np.fft.ifft(np.fft.fftshift(tmp_))),-theta[i])

plt.imshow(res)

But it doesn't work (I m probably missing something?)

Solution

Rotating a single line in a 2D discrete image is really hard. You always end up with a rough approximation, interpolation doesn't help much.

The process you intend to follow is (I've added the filtering):

For each projection in the Radon transform:
- apply the FFT
- apply the wedge filter
- write this as a row though the origin of a 2D complex image
- rotate this image to match the orientation of the projection
- accumulate the result in an output frequency-domain image by summation
Apply the 2D IFFT to the frequency-domain image to obtain the reconstructed image

Because we know that the IFFT operation commutes with the summation, it is possible to move the IFFT operation into the loop:

For each projection in the Radon transform:
- apply the FFT
- apply the wedge filter
- write this as a row though the origin of a 2D complex image
- rotate this image to match the orientation of the projection
- Apply the 2D IFFT
- accumulate the result in an output spatial-domain image by summation

Also the rotation and the IFFT operations commute, so the above is identical to:

For each projection in the Radon transform:
- apply the FFT
- apply the wedge filter
- write this as a row though the origin of a 2D complex image
- Apply the 2D IFFT
- rotate this image to match the orientation of the projection
- accumulate the result in an output spatial-domain image by summation

In this latter case, we are rotating a spatial-domain image that is smooth; it is not a single line drawn in an otherwise empty image, it is a fully band-limited function that can be interpolated into properly. The rotation result is much better in this case.

This latter process is almost the same as what the back projection algorithm does. We can further realize that the 2D IFFT of an image with a single row of data through the origin (the rest of the image is all zero) is the same as taking a 1D IFFT, and replicating that across all rows of the image. This saves quite a bit of computation.

Here is some code. The first method would be (a few fixes from OP's code, but the output is still not recognizable!):

import numpy as np
import matplotlib.pyplot as plt
from skimage.io import imread
from skimage.data import shepp_logan_phantom
from skimage.transform import radon, rescale
from skimage.transform import iradon
from skimage.transform import rotate
import cv2 

x = shepp_logan_phantom()
x = cv2.resize(x, (128,128), interpolation = cv2.INTER_AREA)
theta = np.linspace(0, 180, len(x), endpoint=False)
R = radon(x, theta)

filter = np.abs(np.fft.fftfreq(128))

fullFft2D = np.zeros((128,128)).astype(np.complex128)
for i in range(len(theta)):
    temp_ = np.zeros((128,128)).astype(np.complex128)
    temp_[64,:] = np.fft.fftshift(filter * np.fft.fft(R[:,i]))
    fft_real = rotate(np.real(temp_), theta[i], order=0, center=(64,64))
    fft_imag = rotate(np.imag(temp_), theta[i], order=0, center=(64,64))
    fullFft2D += fft_real + 1j*fft_imag

y = np.fft.ifft2(np.fft.ifftshift(fullFft2D))
plt.imshow(np.fft.fftshift(y.real)); plt.show()

Fixes include: (1) The origin in the fftshifted frequency domain of size 128 is at 64, not 63. (2) The rotation is explicitly performed around the origin. (3) OP had a typo: fft_real + 1j* fft_real. (4) Added the wedge filtering. (5) Not including 180 degrees in the Radon transform (as it's identical to 0 degrees). (6) Using the real part of the IFFT, not the absolute value.

When computing through the frequency domain, if you expect a real-valued result but get non-trivial (trivial==almost zero) imaginary component, something is wrong. In the code above, the imaginary component is non-trivial. This is the result of the rotation of data that cannot be interpolated into properly. The rotation just destroys the changes of success.

The latter method would be:

y = np.zeros((128,128))
for i in range(len(theta)):
    tmp_ = np.zeros((128,128)).astype(np.complex128)
    tmp_[0,:] = filter * np.fft.fft(R[:,i])
    y += rotate(np.fft.ifft2(tmp_).real, -theta[i], center=(64,64))

plt.imshow(y); plt.show()

This code is simplified somewhat because we don't need to use fftshift, we can write the line directly in the origin as expected by the FFT (row 0). The result produced reproduces the phantom correctly.