Radon-Transformation Matrix-Darstellung

Question

Ich bin auf der Suche nach einer MATLAB Lösung, um die Matrix-Darstellung eines diskreten Radon transform (DRT) zu generieren. Das heißt, bei einer vektorisierten Version eines MxN-Bildes, X, möchte ich die Matrix R so erzeugen, dass R*X(:) eine DRT des Bildes ist. In MATLAB, erwartete ich es bin so etwas wie folgt aussehen:

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

Ich weiß, es gibt mehrere Möglichkeiten, eine DRT zu definieren, also werde ich nur sagen, dass ich für einen normalen oder Standard Suche oder nicht-out-of-the-ordinary 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; 

Hier ist die imggrid Funktion, die oben verwendet wurde.

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,[]));