1

I'm doing an exercise to reconstruct an CT-scan image from its sinogram (given).

From the sinogram, implementing the central slice theorem, I did the 1D FFT (after zero padding)

sino_fft = fftshift(fft(fftshift(sinogram)));

and recover the back-projection through the 2D inverse FFT

backprojection = fftshift(ifft2(fftshift(sino_fft)));

Now I should convert the polar back projection to cartesian coordinate and reconstruct the grayscale image. I'm not sure how do it. backprojection is a double complex, so I converted into x-y with

theta = 0:180/(size(backprojection,1)-1):180;
[x,y] = pol2cart(theta,abs(backprojection));

and at this point I have a scattered data set. The hint from the exercise is to sum iteratively from all the angles to get the recovered image. But I don't understand. Just do image = x+y? In this case I get an empty gray image.

enter image description here

Edit after @Christoph Rackwitz comment:

If I remove the inner fftshift from sino_fft and backprojection I get a splitted fourier transform.

enter image description here

while with my code I get a centered fourier transform.

enter image description here

Shika93
  • 635
  • 4
  • 24
  • This seems to be some sort of tomography from your words, but you have not specify which one. That would help. Also it would help describe which algorithm are you trying to implement in particular. i.e. a backprojection of CT (which your words kinda fit) is not a 2D IFFT. – Ander Biguri Jun 18 '23 at 10:13
  • Yes sorry. It's a CT scan tomography. I added few lines. I've a sort of guide for the exercise saying: compute 1D FFT of the sinogram, filter the result to improve the quality (I skipped it in my question), compute the 2D IFFT to get the back projection and convert it from polar to cartesian to get the final image – Shika93 Jun 18 '23 at 10:48
  • I think the innermost fftshift is wrong. does not require application at all, and only makes sense in the frequency domain, not the space domain. the code is incomplete. `sino_filtered` is undefined. – Christoph Rackwitz Jun 18 '23 at 10:49
  • Removed `sino_filtered`. Is useless to the question, sorry. I copied-pasted from matlab. If i remove the inner ffshift I get a splitted image. I add a plot with your tip. – Shika93 Jun 18 '23 at 10:55
  • 1
    well that's progress. remember, fftshift may only be applied to spectra, not spatial domain, and it merely rolls _the spectrum_ so DC is in the center. this merely makes application of masks/filters easier, DC being in the center instead of the corners. – Christoph Rackwitz Jun 18 '23 at 11:30
  • I also saw this [topic](https://stackoverflow.com/questions/61260808/implementing-a-filtered-backprojection-algorithm-using-the-central-slice-theorem) doing something very similar. I skipped the filtering part in my question because it is not neccessary, but clearly it is needed before doing the ifft. I need to understand how from polar go to cartesian to check the result (because plotting the back-projected image is not helping) – Shika93 Jun 18 '23 at 11:40

0 Answers0