RADON Compute Radon transform.
    The RADON function computes the Radon transform, which is the
    projection of the image intensity along a radial line
    oriented at a specific angle.
 
    R = RADON(I,THETA) returns the Radon transform of the
    intensity image I for the angle THETA degrees. If THETA is a
    scalar, the result R is a column vector containing the Radon
    transform for THETA degrees. If THETA is a vector, then R is
    a matrix in which each column is the Radon transform for one
    of the angles in THETA. If you omit THETA, it defaults to
    0:179.
 
    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.
 
    [R,Xp] = RADON(...) returns a vector Xp containing the radial
    coordinates corresponding to each row of R.
 
    Class Support
    -------------
    I can be of class double or of any integer class. All
    other inputs and outputs are of class double.
 
    Remarks
    -------
    The radial coordinates returned in Xp are the values along
    the x-prime 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).

    See also IRADON, PHANTOM.

 
    Example
    -------
%%%
% Adapted from standard example of the RADON function
% by Tomas Svoboda
%%%
I = zeros(100,100);
I(20:85, 25:75) = 1;
I = edge(I,'canny');
theta = 0:179;
[HA,xp] = radon(I,theta);
figure(1)
clf
imshow(theta,xp,HA,[],'n')
xlabel('\theta (degrees)')
ylabel('x''')
colormap(hot), colorbar
axis on 
figure(2)
clf 
imshow(I)
axis on;

[imr,imc] = size(I);
cr = floor( (imr+1)/2 );
cc = floor( (imc+1)/2 );

[r,c] = find(HA==max(HA(:)));
xpmax = xp(r);
thmax = (pi/180)*(theta(c));
if abs(thmax)>0
  x= [-(cc-1):cc];
  y= (xpmax - x*cos(thmax))/sin(thmax);
else
  y= [-(cr-1):cr];
  x= (xpmax - y*sin(thmax))/cos(thmax);
end

figure(2)
hold on
plot(x+cc,imr-(y+cr)+1,'y-','LineWidth',2)
plot(cc,cr,'g*','MarkerSize',8)
axis on