ラドン変換行列の表現

Question

Radon transform（DRT）の行列表現を生成するために、MATLABソリューションを探しています。つまり、MxN画像のベクトル化されたバージョンXが与えられた場合、Rという行列を生成して、R*X(:)が画像のDRTであるようにしたいとします。 MATLABで、私はそれは、次のようなものを見て期待しています：

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

私はDRTを定義するには、いくつかの方法があります知っているので、私は、私は通常か標準を探していますということだけ言いますよまたは普及していません。埋め込みです。

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; 

ここでは、上記で使用されたimggrid関数です。

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