Radon Transform Documentation Matlab

Hello,Regards...

I Have a question about the Radon Transform...exactly the way which Matlab explain this graph:

"The algorithm first divides pixels in the image into four subpixels and projects each subpixel separately, as shown in the following figure.Each subpixel's contribution is proportionally split into the two nearest bins, according to the distance between the projected location and the bin centers. If the subpixel projection hits the center point of a bin, the bin on the axes gets the full value of the subpixel, or one-fourth the value of the pixel. If the subpixel projection hits the border between two bins, the subpixel value is split evenly between the bins."

I understand the Radon Transform get the product between the length(intersection) of each beam through the pixel and the pixel value (Proj=Area_intersection*pixel_value)to generate the projection to differentes angles and next the sum of every projection in every pixel (Sum along of one direction=Projection) but I don't understand this part of documentation of matlab....What refer to 4 subpixels? How can I divide 1 pixel into 4 ?? or Do It make reference to the arrounding area of every pixel??

Thank you!

Answer 1

[Image Analyst]

Think of the Radon transform like an x-ray. It makes a projection along a line. Now the digital image is spatially discrete -- quantized. So what if the projection is at an angle that doesn't line up with pixel rows and columns? What pixels are you going to project? The pixel that should be projected is mathematically at some floating point (x,y) location, and this is not usually exactly on an existing pixel center of your existing spatially discrete image. So to get that pixel you have to do an interpolation to find out what the pixel value should be, if there were a pixel there.

Answer 2

[Matt J]

I understand the Radon Transform get the product between the length(intersection) of each beam through the pixel and the pixel value (Proj=Area_intersection*pixel_value)

So, no, the radon transform is fundamentally a continuous-space operation. It is the integral through the continuous-space object along various lines. There are no "pixels" in a continuous-space operation.

What you have described is one way of discretizing the line integrals done by the radon transform when you have a pixelized object. It is the discretization used, for example, in Siddon's algorithm.

RL Siddon, "Fast calculation of the exact radiological path for a three-dimensional CT array". Med Phys. 2005

However, it is not the only way to discretize the Radon transform and not the method that Matlab's radon() command uses. Essentially, radon uses an approximate forward projection that is the transpose of the backprojection (minus the filtering) done by iradon.

What refer to 4 subpixels? How can I divide 1 pixel into 4 ??

It is not to be done by you. It is done internally by radon(). Basically, radon does a 2x2 upsampling of the input image before doing the forward projection. It is meant simply to get a less coarse forward projection result.

[Matt J]

However, it is not the only way to discretize the Radon transform and not the method that Matlab's radon() command uses.

And note that there are shortcomings to the radon() command's discretization as compared to Siddon's method, which is why many people do not use radon(), except for didactic purposes. For example, the square below is a smooth object, or at least it should have very smooth line integrals. However, we can see in the plot of its projection at 45 degrees (according to radon) that there are noticeable bumps and jags. You wouldn't see that with Siddon's algorithm.

>> P=radon(A,45);
>> plot(P,'0-');ylabel 'Projection';xlabel 'Detector Element';

[Julián Mauricio Ruiz Echeverri]

Thank you! I don't knew about the "Siddon"... so, in order to the discretization RTransf is better to use for example the Siddon algorithm to obtain a more accurate projection? and Does what Matlab use dor that discretization RTransf ?

[Image Analyst]

Matt, how did you get that? Because I can't seem to reproduce that shape:

imageWidth = 300;
squareWidth = 32;
A = zeros(imageWidth, imageWidth, 'double');
r1 = imageWidth/2 - squareWidth/2
r2 = r1 + squareWidth
A(r1:r2, r1:r2) = 4.5;
subplot(2, 1, 1);
imshow(A);
title('Original Image');
impixelinfo;
axis('on', 'image');

subplot(2, 1, 2);
P = radon(A, 45);
plot(P,'bo-');
ylabel 'Projection at 45';
xlabel 'Detector Element';
title('Radon Transform at 45 Degrees');
grid on;
xlim([180, 250]);

There's maybe a slight jiggle between x=208 and x=209 but not as severe as what you showed.

[Priyash Singh]

You might see more severe jiggles if you increase the size of the central square.

[Matt J]

Yes, for example with,

squareWidth = 197;

you will see something much closer to what I originally posted.