Learning Tomography: Fourier Rec. - Part II

In this series dedicated to Direct Fourier Reconstruction technique, here is the second (and final) part of the implementation corresponding to the filling of the 2D Fourier space.


Figure 1 summarizes the last two steps, we need to get the reconstruction image.

Fig.1: From the stack of 1D FTs, computation of the reconstruction image by filling a 2D Fourier space according to the central slice theorem.

According to Fig.1, the algorithm is composed of two parts:

1. Inserting each 1D-FT  into the 2D Fourier space.
2. Inverse transform of the 2D Fourier space to get the final 2D reconstruction image.

1- Step #1: Polar to Cartesian conversion

The formula to convert polar to cartesian coordinates [Wikipedia] is:
 
x = radius * cos(θ)
y = radius * sin(θ)with radius ∈ [-width/2;+width/2.0]

Moreover, as the origin of the 2D FT space is the center of the image, we must shift the x- and y-coordinates  ...

x = xorg + radius * cos(θ)
y = yorg + radius * sin(θ)

with xorg = 2DFT_space.width /2.0
and yorg = 2DFT_space.height /2.0

Now, the synopsis of this algorithm can be written in pseudo-code as follows:

for θ=0 to FT_sinogram.height
  for radius=-FT_sinogram.width/2 to FT_sinogram.width/2
    {x,y} = compute_cartesian_coordinates_from(radius,θ)
    value = FT_sinogram.getPixel(radius,θ)
    2DFT_space.setPixel(nearest(x),nearest(y),value)
  end for
end for

Note: Obviously, the x and y calculated values are floating point values and we need integer X- and Y-coordinates to correctly set the pixel in the Fourier image. Here, I choose the nearest integer value by using the Math.round(...) JavaScript function.

2- Step #2: Transforming the 2D Fourier Space into an image

Finally, the 2D Fourier space (Fig. 2A) must be transformed in real space by using the inverse 2D Fast Hartley Transform available in ImageJ (class FHT). To get a correct computation, it is required to swap the quadrants before running the inverse FHT (Fig.2B) and then, to swap again (Fig.2C) in real space to display the final 2D reconstruction (Fig. 2D).

Fig.2: Inverse Fourier Transform to get the 2D reconstruction.

3- Script

Here is the script...

+++ JavaScript IJ snippet +++ +++ End of  JavaScript IJ snippet +++

4- Results

From the sinogram [Link] transformed by FourierRecSinogram.js [Link], this script fills the Fourier space (Fig.3A), then swaps the quadrants before computing the inverse of the Hartley transform (Fig. 3BC).

Fig.3: Conversion of the 2D Fourier space into real space. Intermediate steps composed of 'swap quadrants' are required to get the final image (Fig.3).

Then, by swapping once more the image of Fig. 3C,  the image of Fig. 4. is finally correctly rendered...

Fig.4: Direct Fourier Reconstruction of Lena.

In conclusion, we were able to reconstruct in frequency space, the test image from the sinogram. Of course, the result isn't entirely satisfactory because even though Lena is recognizable without blur effect, the background appears 'foggy'.
Improving this first attempt of Direct Fourier Reconstruction is relatively easy and this will be the subject of a next post. To be continued ...

<< Previous: 1D-FT of sinogram

>> Next: Improvement: Zero padding

Other crazybiocomputing posts

Further readings are available in ...

• Learning Tomography Series  [Link]
• Image Processing TOC [Link]