2

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

Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
Chenyang Zhao
  • 21
  • 2
  • 1
    Help you to do what? – Adam May 18 '19 at 17:13
  • I need help on understanding how this pattern happen. I don't think this is normal reconstruction artifact. Because if it is, the artifact should decrease as the increasing of number of projections. However, no matter how many projections that I used, this pattern is still there. – Chenyang Zhao May 19 '19 at 03:25
  • I found there are saw-tooth artifact in the 45 and 135 degree projection data that was got from radon function. And this saw-tooth artifact only exist in these two projection data. This might be the source of the pattern in the image. Please try this code: ```matlab I = zeros(1000,1000); I(250:750, 250:750) = 1; theta = [0 45 90 135]; [R,xp] = radon(I,theta); figure;plot(R);legend('0°','45°','90°','135°'); ``` So, why there are saw-tooth artifact existing in 45 and 135 degree projection data? – Chenyang Zhao May 19 '19 at 03:34
  • Interesting. You're right. The diagonal projections from `radon` have some bad sampling artifacts. It looks like it uses a pretty naive sampling scheme. – Cris Luengo May 19 '19 at 05:04
  • "This is even more obvious if I use my CT image dataset." Is this a sonogram obtained in a CT scanner? Or are you using `radon` with that data set too? Also, I hope you are not trying to use `iradon` to reconstruct a real-world CT scan. `radon` and `iradon` are meant as a demonstration at best. Don't use these functions for anything serious. – Cris Luengo May 19 '19 at 05:17
  • No, I did simulation on my CT image data. So I used radon to get the sinogram. Then, I used my own reconstruction method to reconstruct image. I tried FBP, my own method, and iradon. All of them have the same issue, which makes me suspect the problem is in the projection obtaining part. – Chenyang Zhao May 20 '19 at 06:58
  • To have full fourier sampling of the object you need `pi*max(size(img))` angles, otherwise you are sampling less and you will have artifacts – Ander Biguri May 21 '19 at 10:24

1 Answers1

3

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

rvimieiro
  • 1,043
  • 1
  • 10
  • 17
  • I guess this saw-tooth patter that you noticed is due to the sub-pixel sampling, as stated [here](https://www.mathworks.com/help/images/ref/radon.html#f6-466342). So when the projection is at, e.g., 45 degrees, the attenuation contribution of the two pixel that are in the diagonal are counted for a single detector, while the other two split for the other two detectors. I drew this [picture](https://i.imgur.com/tSruRsA.png) as an example. – rvimieiro May 20 '19 at 12:54
  • That makes sense. Though I think the projection pixel distance is 1, in your figure it is sqrt(2)/2. – Cris Luengo May 20 '19 at 13:21
  • Thanks a lot. It makes sense to me. Can I understand what you said in this way that there is no interpolation at 45 and 135 degree along the diagonal directions of pixels? Do you think the distance driven projection method can improve? – Chenyang Zhao May 20 '19 at 15:56
  • There is an interpolation, however, the path on the diagonal is greater in the middle compared with the extremities, as I stated in the figure. That is why the saw-tooth is visible on the angular views. Let me know if it is clear to you. These methods that I mentioned are overall better, so you can test them and see if it is better for you (which I think it will be). – rvimieiro May 20 '19 at 19:50
  • Thanks. [link](https://github.com/rodrigovimieiro/OpenCodes/blob/master/Tomography/Computed%20Tomography%20(CT)/DistanceDrive/2D/NaiveDistanceDriven.m) Is this code for fanbeam projection or parallel beam? It looks like fanbeam, but you used iradon to reconstruct image. – Chenyang Zhao May 21 '19 at 05:49
  • I'm still confused. If it is simply caused by longer path of 45, why other angles that are very close to 45, for example 46, don't have the same issue? – Chenyang Zhao May 21 '19 at 08:47
  • 1
    @ChenyangZhao this is not the only reason. The filtered backprojection algorithm is only good in with infinite sampling which is not possible. This is why iterative methods are being pushed by researchers, because they can deal with this stuff better. – Ander Biguri May 21 '19 at 10:09
  • 1
    I used `iradon` because I set the geometry parameters to simulate a parallel acquisition. Refer to the [article](https://iopscience.iop.org/article/10.1088/0031-9155/49/11/024/meta) of DD, so you can understand better how it works. The saw-tooth is caused by this long path, however, there are other inconvenient as stated by @AnderBiguri. Let us know if it is clear to you. – rvimieiro May 21 '19 at 13:35