0

Due to the nature of my problem, I want to evaluate the numerical implementations of the Radon transform in Matlab (i.e. different interpolation methods give different numerical values).

while trying to code my own Radon, and compare it to Matlab's output, I found out that my radon projection sizes are different than Matlab's.

So a bit of intuition of how I compute the amount if radon samples needed. Let's do the 2D case.

The idea is that the maximum size would be when the diagonal (in a rectangular shape at least) part is proyected in the radon transform, so diago=sqrt(size(I,1),size(I,2)). As we dont wan nothing out, n_r=ceil(diago). n_r should be the amount of discrete samples of the radon transform should be to ensure no data is left out.

I noticed that Matlab's radon output is always even, which makes sense as you would want a "ray" through the rotation center always. And I noticed that there are 2 zeros in the endpoints of the array in all cases.

So in that case, n_r=ceil(diago)+mod(ceil(diago)+1,2)+2;

However, it seems that I get small discrepancies with Matlab.

A MWE:

 % Try: 255,256
 pixels=256;
 I=phantom('Modified Shepp-Logan',pixels);

 rd=radon(I,pi/4);
 size(rd,1)

 s=size(I);
 diagsize=sqrt(sum(s.^2));
 n_r=ceil(diagsize)+mod(ceil(diagsize)+1,2)+2

rd=

   367


n_r =

   365

As Matlab's Radon transform is a function I can not look into, I wonder why could it be this discrepancy.

Ander Biguri
  • 35,140
  • 11
  • 74
  • 120
  • Would it make sense to post your version of the radon transform so we can compare? – NoDataDumpNoContribution Feb 05 '15 at 09:29
  • @Trilarion I don't fell like is needed, as my point is to code Matlab's Radon transform in order to know how it does it, so I can modify that code if I want. Any further code in the radon transform would have not influence in the preallocated size. – Ander Biguri Feb 05 '15 at 13:14

2 Answers2

3

I took another look at the problem and I believe this is actually the right answer. From the "hidden documentation" of radon.m (type in edit radon.m and scroll to the bottom)

Grandfathered syntax

R = RADON(I,THETA,N) returns a Radon transform with the projection computed at N points. R has N rows. If you do not specify N, the number of points the projection is computed at is:

   2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))+3

This number is sufficient to compute the projection at unit intervals, even along the diagonal.

I did not try to rederive this formula, but I think this is what you're looking for.

Ander Biguri
  • 35,140
  • 11
  • 74
  • 120
purple arrow
  • 141
  • 2
1

This is a fairly specialized question, so I'll offer up an idea without being completely sure it is the answer to your specific question (normally I would pass and let someone else answer, but I'm not sure how many readers of stackoverflow have studied radon). I think what you might be overlooking is the floor function in the documentation for the radon function call. From the doc:

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)

For example, in a 20-by-30 image, the center pixel is (10,15).

This gives different behavior for odd- or even-sized problems that you pass in. Hence, in your example ("Try: 255, 256"), you would need a different case for odd versus even, and this might involve (in effect) padding with a row and column of zeros.

purple arrow
  • 141
  • 2
  • I though about that, indeed. And you are rigth, maybe nobody can answer my question hehe, thanks for the try! However, that mistmathc does not happen in all even numbers, so I dont understand why it happens some times. – Ander Biguri Feb 04 '15 at 21:42