Search code examples
matlabimage-processingtomography-reconstruction

Moire pattern shows up in the reconstructed image while using radon and iradon

I used this script to reconstruct an image of the Shepp-Logan phantom.

Basically, it just simply used radon to get sinogram and used iradon to transform it back.

However, I found that a very obvious moire pattern can be seen when I adjust contrast. This is even more obvious if I use my CT image dataset.

Can anyone help me to understand this? Thanks!

img = phantom(512)*1000;
views = 576;
angles = [0:180/576:180-180/576];
sino = radon(img,angles);
img_rec = iradon(sino,angles);
imshow(img_rec,[]);

Full image after being adjusted contrast:

full image after being adjusted contrast

Regions with obvious moire pattern:

regions with obvious moire pattern

Solution

  • This may be happening because of some factors:

    • From the MATLAB documentation, iradon uses 'Ram-Lak' (known as ramp filter) filtering as default and does not use any windowing to de-emphasizes noise in the high frequencies. You stated "This is even more obvious if I use my CT image dataset", that is because there you have real noise in the images. The documentation itself advises to use some windowing:

    "Because this filter is sensitive to noise in the projections, one of the filters listed below might be preferable. These filters multiply the Ram-Lak filter by a window that de-emphasizes higher frequencies."

    • Other inconvenient is related to the projector itself. The built-in functions radon and iradon from MATLAB does not take into account the detector size and the x-ray length which cross the pixel. These functions are just pixel driving methods, i.e., they basically project geometrically the pixels in the detector and interpolate them.

    Possible solutions:

    There are more sophisticated projectors today as [1] and [2]. As I stated here, I implemented the distance-driven projector for 2D Computed Tomography (CT) and 3D Digital Breast Tomosynthesis (DBT), so feel free to use it for your experiments.

    For example, I generated 3600 equally spaced projections of the phantom with the distance-driven method, and reconstructed it with the iradon function using this line code:

    slice = iradon(sinogram',rad2deg(geo.theta));

    Here