Radon transform in matlab

Hi Everyone, Is there any numerically efficient package to compute the Radon transform and the back projection, probably using the non-uniform forier transform!

actually what I am looking for is a forward and inverse projector for tomography problem. I tried the Matlab Radon function but I noticed that if you apply it iteratively(as in the code below) the Image will be bleared! are there any solution for this problem?

x = phantom('Modified Shepp-Logan',256);
[Sinogram,xp] = radon(x,1:0.1:180);
for i=1:10
I = iradon(Sinogram, 1:0.1:180);
[Sinogram,xp] = radon(I,1:0.1:180);
end
imtool(I);

Answer 1

[Matt J]

but I noticed that if you apply it iteratively(as in the code below) the Image will be bleared!

Reconstruction is always lossy because all the projection and DSP operations involved are discretized, rather than occuring in ideal, continuous space. If you iteratively project and reconstruct, the loss accumulates.

[Alan AF]

Thank you Matt for the clarification, but how can we minimize this effect? and would you recommend an accurate and fast Radon transform?

I tried NUFFT by Fessler but unfortunately I could not know how to implement the Radon Transform using FFT?

[Matt J]

but how can we minimize this effect?

By not iterating. As Bjorn also said, the iterations do not have any obvious purpose except to degrade the result.

Jeff Fessler's library is what I would have recommended. He has a library of different forward projectors to choose from.

[Alan AF]

Hi Matt, what I am trying to do is to implement an iterative reconstruction algorithm for CT tomography by using approach as in the pseudo-code below:

Sinogram= under_sampled_radon_transform
I0 = iradon(Sinogram, undersampled_theta);% initial guess
for i=1:10
I=Enhance(I0);% apply operation in the pixel domain
[I_Sinogram,xp] = radon(I,1:0.1:180);
new_Sinogram=(I_Sinogram+Sinogram);% update solution the sinogram domain
I = iradon(new_Sinogram, 1:0.1:180);
end
imtool(I);

the problem is that I don't have an efficient forward/inverse projector! I tried the NUFFT by Fessler, it seemed fast but the code was complicated for me to understand. I did some attempts but still need an expert opinion, would you consider having a look on the code that I written?

[Matt J]

Does it matter that it's a black box? Many built-in MATLAB commands are black boxes, as well. Why not just run some tests to see if the output is what you expect from a given input?

[Alan AF]

yes true, I already did that, and here what happens:

after each iteration, the fine features of the images get higher pixel values till every thing disappear in the imageon the 10th iteration (each with, Radon and FBP being applied on.) it seamed that the filter that being used is causing such a problem. I don't know what to do, this is what I need some help with Fesler's NUFFT!

[Matt J]

What makes you think it is the NUFFT's fault? Did you try a different forward projector (e.g., MATLAB's radon command)? Is the result different?

[Alan AF]

Yes, I have tried the Matlab radon function as below, and the result was a blurred image, however the NUFFT code results was a blank image with vertical lines?! would you plz have a look on NUFFT code that I wrote ( download )

x = phantom('Modified Shepp-Logan',256);
[Sinogram,xp] = radon(x,1:0.1:180);
for i=1:30
I = iradon(Sinogram, 1:0.1:180,'linear','Ram-Lak');
[Sinogram,xp] = radon(I,1:0.1:180);
end
imtool(I);

[Matt J]

It won't mean anything to me. I've never worked with NUFFT projectors.

Answer 2

[Bjorn Gustavsson]

Why would you want to repeat the process in the first place? If you have data for every 0.1 degree (according to your calculation of the sinogram) then the iradon function calculates the proper inverse (you might select different filters for the transform and such), therefore you'll only enhance the problems with the inverse (blurring, pixeling, aliasing and possibly noise amplification if you have noise in your real data and whatnot), put an image display inside the loop together with a drawnow and see how things evolve...

If you want to look at some iterative algorithm you have to explain what you want to study.

HTH,

[Alan AF]

HI Bjorn, what I am tring to do is to implement an iterative reconstuction algorithem for CT tomography by using aproach as in the pusedo-code below:

Sinogram= under_sampled_radon_transform
I0 = iradon(Sinogram, undersampled_theta);% initial guess
for i=1:10
I=Enhance(I0);% apply operation in the pixel domain
[I_Sinogram,xp] = radon(I,1:0.1:180);
new_Sinogram=(I_Sinogram+Sinogram);% update solution the sinogram domain
I = iradon(new_Sinogram, 1:0.1:180);
end
imtool(I);

the problem is that I don't have an efficient forward/inverse projector! I tried the NUFFT by Fessler, it seemed fast but the code was a black box for me this is why I am asking for help! Thank you

[Bjorn Gustavsson]

Since I have no idea what your Enhance function does I cant see what you want to achieve. Perhaps you want to do some kind of high-pass filtering or the other? If that operation is isotropic or linear enough you might get away with transforming that operation into some kind of filter-function that you could pass into iradon (or edit the function, I dont remember what inputs you can give that function.) It looks a bit strange that you want to use a much denser/more complete angular radon-transform in your iteration since your initial data is apparently under-sampled (in what way? sparse angles or restricted angular coverage?)

So to be able to have further suggestions we need more information...

HTH

[Alan AF]

The Enhance function is a non-linear operation that is applied to the image (under reconstruction) in the image domain, basically it is a kind of operator that uses prior information (trained dictionary, Total variation minimization,...) about the objects in that image in-order to decrease the effect of missing wedge problem and the noise in a limited angle tomography problem. The reason I used the denser/more complete angular radon-transform (0.1) is to decrees the artefact that might arise because of lower angular steps.

I tried the NUFFT code of Jeff Fessler's, results was a blank image with vertical lines?! would you please consider having a look on that code ( download ) I am not sure why the high frequencies is being magnified in the back projector! probably it is a scaling issue?

[Matt J]

probably it is a scaling issue?

You should look at the forward projected views and see if/how the values are globally scaled relative to what radon produces.

Some things to think about: When using radon, what scaling are you doing to account for the physical sizes (in millimeters) of the image pixels. Similarly, what is NUFFT doing?

[Alan AF]

Hi Matt would you elaborate more about the scaling for the physical sizes, do you mean resolution?

[Matt J]

The radon transform is a line integral through the image. Roughly, it works by finding the intersection lengths of projection rays with pixels and multiplies by the length of intersection. But to compute the length of intersection, you have to know the pixel dimensions, e.g., in millimeters. Matlab's RADON command just assumes the pixels dimensions are 1x1 unitless. You would have to scale the result by the pixel lengths to get things in real units. I don't know whether you are doing that currently

But the Fessler NUFFT code might be drawing pixel dimension information from somewhere already, e.g., input you gave it.