Backprojection from projection matrix using MATLAB

Question

I've a 256x256 projection matrix. each row is a projection taken with equal angles. i need to generate the original image with backprojection using matlab and I am not really familiar with matlab. Can you suggest me any code samples or alghorithms? I've found some similar codes i couldn't generate the original image using them.

Answer 1

This should be relatively simple with the iradon command, if you have the Image Processing Toolbox. If you don't, it will be a bit tougher because you need to roll your own version of that. Apparently you can't use this, but for what it's worth I get an image if I use:

I = iradon(Pteta',linspace(0,179,size(Pteta,1));

So, how can you do this yourself? I'll try to help you along the way without giving you the answer - this is homework after all!

First, think about your 0-degree projection. Imagine the axis you're projecting on has units 1,256. Now imagine back projection of these coordinates across your image, it would look something like this:

Similarly, think a 90-degree projection would look like this:

Cool, we can get these matrices by using [X, Y] = meshgrid(1:256);, but what about off-axis projections? Just think of the distance along some angled line, like converting polar/Cartesian coordinates:

theta = 45 % projection angle in degrees
t = X*cosd(theta) + Y*sind(theta);

And it works!

There is a problem, though! Notice the values go up over 350 now? Also it's sort of off-center. The coordinates now exceed the length of our projections because the diagonal of a square is longer than the side. I'll leave it to you to figure out how to resolve this, but figure the final image will be smaller than the initial projections, and you may need to use different units (-127 to 128 instead of 1 to 256).

Now you can just index your projections for those angles to backproject the actual values across the image. Here we have a second problem, though, because the values are not integers! We could just round them, this is called nearest-neighbor interpolation, but it doesn't give the best results.

proj = Pteta(angle,:);
% add projection filtering here
t = X*cosd(theta) + Y*sind(theta);
% do some rounding/interpolating to make t all integers
imagesc(proj(t));

For our off-center version, that gives us this image, or something similar:

Now you just need to do this for every angle, and add them all up.