How the "radon" transform works in matlab????

I need to find the thickness of a segmented object that is in binary image. So i am trying to use 'radon' transform to find the intensity values with the value as "1" (white color) along the specified orientation. In this way i can find thickness at specified locations as i want. Before using this function 'radon' i want to understand how it works. But i am not able to understand how radon works. As an example i took one simple image and coded as follows.

 I = zeros(100,100);
I(1:25,1:25) = 1;
figure;imshow(I)
[R,xp] = radon(I, 90);
figure;plot(R)

I plotted R and checked for theta = 90,0. The plot shows impulse wave for both and there is a shift for theta = 0. The size® = 145. How does R takes 145, when the size of image is only 100? How could even for theta = 0, R gives some intensity values?? Can anyone give me a clear understanding on this radon??

Answer 1

[Matt J]

How does R takes 145, when the size of image is only 100?

RADON is choosing a grid size for R that will capture the projections of the entire 100x100 grid at all angles theta. At theta=45, the length of the grid's projection shadow is sqrt(2)*100 and so something like R=145 is required to cover it.

How could even for theta = 0, R gives some intensity values??

You have to explain why you don't think R should have intensity values. Your 25x25 square is viewable at all angles, theta.

[Divya]

Thank You Matt J. I have few more doubts for your answers

1) Why the length of the grid's projection is sqrt(2)*100? Is it by applying trigonometric ratios for the right angle triangle formed by the projection? If it so, i found it to be 100/sqrt(2). I applied cosine ratio to find the length of the projection. Because the help in matlab says that

" The radial coordinates returned in xp are the values along the x'-axis, which is oriented at theta degrees counterclockwise from the x-axis. The origin of both axes is the center pixel of the image, which is defined as floor((size(I)+1)/2) "

So i took the center of the image, I(50,50) as the origin of the axes and i assumed the angle is rotating in counterclockwise starting from the right side. In that way, for theta = 45, i get a right angled triangle on the right side.

I think i went wrong in taking the origin and orientation. Please correct me in this context

2) For your second answer, i think if i get the proper orientation and origin then i will have intensity values at all angles.

[Matt J]

Why the length of the grid's projection is sqrt(2)*100?

The length of the grid's projection will be the length of the grid's diagonal which is sqrt(2)*100. The radial R-axis coincides with the grid's diagonal at theta=45 and all pixels along it can contribute to the projection.

For your second answer, i think if i get the proper orientation and origin then i will have intensity values at all angles.

You already do have intensity values values at all angles. Whenever there is non-zero intensity in the grid, there will be non-zero intensity in all the radon projections (assuming the default grid size).

[Divya]

Hey Matt, i am very sorry. I din get you on this. I understood that all the pixels on the radial R-axis contribute to projection and the number of pixels is what 145 in the above example for theta=45. But what i need is from where the radial R-axis starts???(I mean where is the origin??) and in what direction does it rotates??? How the value for the length of the grid's projection or grid's diagonal has sqrt(2)*100??

[Matt J]

But what i need is from where the radial R-axis starts???(I mean where is the origin??) and in what direction does it rotates???

The origin of the R-axis is at the center, i.e., at R=73, and that point is aligned with the center of the image grid. It rotates counter-clockwise with theta. You can verify this by calling RADON for several theta and watching how the projections of the 25x25 square move with theta

How the value for the length of the grid's projection or grid's diagonal has sqrt(2)*100??

That's simple Pythagoras. The sides of the grid are of length 100. By Pythagoras, the diagonal length is therefore

 sqrt(100^2 +100^2) = sqrt(2)*100

[Divya]

OK... I got you.. Thank you Matt J

[Divya]

Matt, I have few more doubts By reading radon, (from page http://www.mathworks.in/help/images/radon-transform.html), i understood the working of radon and why the formula of R(theta) came. Lets take much simpler example

I=[1 1 1;
   1 1 1;
   1 1 1];

according to what you said, if theta = 45, the length is sqrt(2)*3. Then if theta = 0, the length of the projection should be 3 right!(since the projection is parallel to y' axis and the length of y' is 3) but when i am checking in matlab, it is giving me 7 values(first and last as 0). Why??

I want to know how the values of R are calculated. So i applied the formula as R = integration(f(x'cos(t)-y'sin(t),x'sin(t)+y'cos(t))dy) where t = theta(0). Is it the correct way to find the values of R??? I am not getting the right values in this way? Could you please clear me in this?

[Matt J]

according to what you said, if theta = 45, the length is sqrt(2)*3. Then if theta = 0, the length of the projection should be 3 right!

No, by default the length of the output will be approximately sqrt(2)*N independently of theta, because that is the worst case length. RADON uses theta=45 to determine the worst case length, because that is the projection view angle where the image grid casts the widest shadow. You can test this using the following code and you will see that the plot of L increases approximately as N.

 N=5:5:1000; 
  for i=1:length(N); 
    L(i)=length(radon(ones(N(i)),45))/sqrt(2); 
  end, 
 plot(N,L)

Note, by the way, that RADON accepts an 'output_size' argument if you don't like this default.

I want to know how the values of R are calculated.

You should read "doc radon" for that. RADON does not compute an exact line integral. It uses a rather crude approximation, in fact, which is why people who work in tomography largely do not use RADON. It should be pretty good for projecting large, piece-wise uniform object, though.

[Divya]

Thank you Matt. this is quite useful info. Problem is that i need to find thickness(of an object in binary format) using radon. Without knowing what the values and where it is coming from, it is quite hard for me to analyse the data(thickness). I looked into "doc radon" i couldn't get the answer there. Can you suggest me whether radon is helpful in finding out the thickness of an object?

[Matt J]

Why don't you try it on an object of a known thickness and see if it is accurate enough for your needs?

[Divya]

I am getting some weird numbers Matt. I am sure if i can know what the values in R signifies, i will be able to get the required output from R(if it is possible). Can you please help me out in understanding about the values of R?

[Matt J]

Well, I don't know what you're seeing and how it compares with what you expect to see. It sounds like you want R to give you the cross-sectional lengths of an object at different angles and the code you posted confirms that that is the case. When I run the code, R shows a peak value of 25 at theta=0, which is fine because lines parallel to the sides of the square cut through the square with length 25. At theta=45, I see a peak of about 35 which is the diagonal length of a square whose sides are length 25. So, your examples all seem to support that RADON is doing what you want it to do.

[Divya]

Yeah Matt, I got that R(ceil(size(R,1)/2),:) can be taken as the thickness as it gives the line integral covering complete pixel but not sub pixel values. Thank you for clearing my doubts