I have implemented the procedure to reconstruct an image from its Radon transform involving the Fourier Slice Theorem (FST). After application of the FST, the data is transformed from polar coordinates to cartesian, by regridding.

It works fine, except the reconstructed image is "zoomed out" (see picture). What causes this difference, and how can I resolve it, so that the phantom in the original, an reconstruction are exactly the same size?

enter image description here


N = 128;
f = phantom(N);

% Radon Transform
theta = 0:1:180;
[g, xp] = radon(f, theta);

% Fourier Slice Theorem
ghatpol = fftshift(fft(ifftshift(g)));

% Transform to cartesian coordinates
phi = -theta*2*pi/360;
[pX, pY] = pol2cart(phi, xp);
x = -N/2:1:N/2-1;
[X, Y] = meshgrid(x, x);
fhatcart = griddata(pX, pY, ghatpol, X, Y);

frec = ifftshift(ifft2(fhatcart));

subplot(121), imagesc(f)
subplot(122), imagesc(abs(frec)) 

MathWorks' image processing toolbox describes the Radon transform as the projection of the image intensity along a radial line oriented at a specific angle. Now, let the set of parallel lines spaced 1 px apart is drawn in a fixed direction (let it be vertical lines) and you rotate the image bounded by the square 128x128 units (in which the phantom slice is inscribed) and for each 1 degree of rotation take the projections of the image (let it be horizontal projection of rotated image), integrating the image intensity (optical density, luminosity, etc) along the parallel lines of your fixed line set. Your object under test, in all its rotation orientations, is bounded not by the 128x128 square, but by a circle of 128·sqrt(2) diameter. To make sure that no region of the source image is missing, the Radon transform is performed for the circle of 128·sqrt(2) diameter.

Now, you perform the reconstruction of all the Radon transform data. The meaningful data are located in the circle of 128·sqrt(2) diameter, not in the 128x128 square, and you inscribe the reconstructed image into in the 128x128 square as is seen in the lines

x = -N/2:1:N/2-1; [X, Y] = meshgrid(x, x); 
fhatcart = griddata(pX, pY, ghatpol, X, Y);  
frec = ifftshift(ifft2(fhatcart));

of your code.

As the result, the structures of the source image are zoomed out by ~1.414.

You can adjust the meshgrid

x = -N/2:1.41421:N/2-1;

and the other parameters accordingly, or zoom in the reconstructed image in other ways, if you need to keep a 100% zoom.

  • $\begingroup$ Thanks for the explanation! But now I am wondering, after adapting the x variable to sqrt(2) step size, my reconstruction now has a different size (90x90). How can I make it to be 128x128 and 100% zoom? In later steps, I want to calculate differences and errors between those images, so I need same grid size as well. $\endgroup$
    – Sim
    May 21 at 12:03
  • $\begingroup$ use x = -round(N/2*1.41421):1.41421:round((N/2-1)*1.41421) resulting in the 181x181 px image. I may be wrong in a N/2*1.41421 part, adjust this value to have the 181x181 image and cut the central part of 128x128 pixels. Verify that it is 100% zoom and use it for your purpose. $\endgroup$
    – V.V.T
    May 21 at 12:13
  • $\begingroup$ or maybe the line must be x = -round(N/2*1.41421):1:round((N/2-1)*1.41421) and you would cut the central part from it. Debug it on your own. Now that you know the reason of you trouble and the course to remedy, you can find the solution, ok? $\endgroup$
    – V.V.T
    May 21 at 12:21
  • $\begingroup$ Yes perfect! I think I got it! Thank you very much! $\endgroup$
    – Sim
    May 21 at 12:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.