I am looking for a MATLAB solution to generate the matrix representation of a discrete Radon transform (DRT). That is, given a vectorized version of an MxN image, X, I'd like to generate the matrix R such that R*X(:) is a DRT of the image. In MATLAB, I am expecting it to look something like the following: 

>> X = 2D_Image_Of_Size_MxN;
>> R = DRT_Matrix_Of_Size_LPxMN;
>> DRT = reshape( R * X(:), L, P );

I know there are several ways to define a DRT, so I'll just say that I am looking for a normal or standard or not-too-out-of-the-ordinary implmentation. 

function [ R rho theta ] = radonmatrix( drho, dtheta, M, N )
% radonmatrix - Discrete Radon Trasnform matrix
%
% SYNOPSIS
%   [ R rho theta ] = radonmatrix( drho, dtheta, M, N )
%
% DESCRIPTION
%   Returns a matrix representation of a Discrete Radon
%   Transform (DRT).
%
% INPUT
%   drho     Radial spacing the the DRT.
%   dtheta   Angular spacing of the DRT (rad).
%   M        Number of rows in the image.
%   N        Number of columns in the image.
%
% OUTPUT
%   R        LP x MN DRT matrix. The values of the L and
%            P will depend on the radial and angular spacings.
%   rho      Vector of radial sample locations.
%   theta    Vector of angular sample locations (rad).
%

% For each angle, we define a set of rays parameterized
% by rho. We then find the pixels on the MxN grid that
% are closest to each line. The elements in R corresponding
% to those pixels are given the value of 1.

% The maximum extent of the region of support. It's for
% rho = 0 and theta = pi/4, the line that runs caddy-corner.
W = sqrt( M^2 + N^2 );

rho = -W/2 : drho : W/2;
theta = 0 : dtheta : 180 - dtheta;

L = length( rho );
P = length( theta );

R = false( L*P, M*N );

% Define a meshgrid w/ (0,0) in the middle that
% we can use a standard coordinate system.
[ mimg nimg ] = imggrid( 1, 1, [ M N ] );

% We loop over each angle and define all of the lines.
% We then just figure out which indices each line goes
% through and put a 1 there.
for ii = 1 : P

  phi = theta(ii) * pi/180;

  % The equaiton is rho = m * sin(phi) + n * cos(phi).
  % We either define a vector for m and solve for n
  % or vice versa. We chose which one based on angle
  % so that we never g4et close to dividing by zero.
  if( phi >= pi/4 && phi <= 3*pi/4 )

    t =  -W : min( 1/sqrt(2), 1/abs(cot(phi)) ) : +W;
    T = length( t );

    rhom = repmat( rho(:), 1, T );
    tn = repmat( t(:)', L, 1 );
    mline = ( rhom - tn * cos(phi) ) ./ sin(phi);

    for jj = 1 : L
      p = round( tn(jj,:) - min( nimg ) ) + 1;
      q = round( mline(jj,:) - min( mimg ) ) + 1;  
      inds = p >= 1 & p <= N & q >= 1 & q <= M;
      R( (ii-1)*L + jj, unique( sub2ind( [ M N ], q(inds), p(inds) ) ) ) = 1;
    end

  else

    t =  -W : min( 1/sqrt(2), 1/abs(tan(phi)) ) : +W;
    T = length( t );

    rhon = repmat( rho(:)', T, 1 );    
    tm = repmat( t(:), 1, L );
    nline = ( rhon - tm * sin(phi) ) ./ cos(phi);

    for jj = 1 : L
      p = round( nline(:,jj) - min( nimg ) ) + 1;
      q = round( tm(:,jj) - min( mimg ) ) + 1;  
      inds = p >= 1 & p <= N & q >= 1 & q <= M;
      R( (ii-1)*L + jj, unique( sub2ind( [ M N ], q(inds), p(inds) ) ) ) = 1;
    end

  end

end

R = double( sparse( R ) );

return;

Here is the imggrid function used in the above.

function [ m n ] = imggrid( dm, dn, sz )
% imggrid -- Returns rectilinear coordinate vectors
%
% SYNOPSIS
%   [ m n ] = imggrid( dm, dn, sz )
%
% DESCRIPTION
%   Given the sample spacings and the image size, this
%   function returns the row and column coordinate vectors
%   for the image. Both vectors are centered about zero.
%
% INPUT
%   dm     Spacing between rows.
%   dn     Spacing between columns.
%   sz     2x1 vector of the image size: [ Nrows Ncols ].
%
% OUTPUT
%   m      sz(1) x 1 row coordinate vector.
%   n      1 x sz(2) column coordinate vector.

M = sz(1);
N = sz(2);

if( mod( M, 2 ) == 0 )
  m = dm * ( ceil( -M/2 ) : floor( M/2 ) - 1 )';
else
  m = dm * ( ceil( -M/2 ) : floor( M/2 ) )';
end

if( mod( N, 2 ) == 0 )
  n = dn * ( ceil( -N/2 ) : floor( N/2 ) - 1 );
else
  n = dn * ( ceil( -N/2 ) : floor( N/2 ) );
end