Ask an expert. Trust the answer.

Your academic and career questions answered by verified experts

Radon transform matrix representation

Date: 2023-01-04 14:18:04

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. 

Answers:

 

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

 


Why Matlabhelpers ?

Our Matlab assignment helpers for online MATLAB assignment help service take utmost care of your assignments by keeping the codes simple yet of high-quality. We offer the most reliable MATLAB solutions to students pursuing their Computer Science course from the Monash University, the University of Sydney, the University of New South Wales, the University of Melbourne; to name a few. Approach us today for best Matlab solutions online!

whatsApp order on matlabhelpers.com

telegram order on matlabsolutions.com