Rappresentazione matrice di trasformazione radon

Question

Sto cercando una soluzione MATLAB per generare la rappresentazione matrice di un discreto Radon transform (DRT). Cioè, data una versione vettoriale di un'immagine MxN, X, mi piacerebbe generare la matrice R tale che R*X(:) sia un DRT dell'immagine. In MATLAB, mi aspetto che a guardare qualcosa di simile al seguente:

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

So che ci sono diversi modi per definire un DRT, quindi mi limiterò a dire che io sto cercando un normale o di serie o non-troppo-ordinario implmentation.

Answer 1

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; 

Questa è la funzione imggrid utilizzata in precedenza.

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 

Answer 2

image = phantom(); 
projection = radon(image); 
R = linsolve(reshape(image,1,[]), reshape(projection,1,[]));