1## Copyright (C) 2007 Alexander Barth <barth.alexander@gmail.com> 2## 3## This program is free software; you can redistribute it and/or modify it under 4## the terms of the GNU General Public License as published by the Free Software 5## Foundation; either version 3 of the License, or (at your option) any later 6## version. 7## 8## This program is distributed in the hope that it will be useful, but WITHOUT 9## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 10## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 11## details. 12## 13## You should have received a copy of the GNU General Public License along with 14## this program; if not, see <http://www.gnu.org/licenses/>. 15 16## -*- texinfo -*- 17## @deftypefn {Function File} {[@var{RT},@var{xp}] =} radon(@var{I}, @var{theta}) 18## @deftypefnx {Function File} {[@var{RT},@var{xp}] =} radon(@var{I}) 19## 20## Calculates the 2D-Radon transform of the matrix @var{I} at angles given in 21## @var{theta}. To each element of @var{theta} corresponds a column in @var{RT}. 22## The variable @var{xp} represents the x-axis of the rotated coordinate. 23## If @var{theta} is not defined, then 0:179 is assumed. 24## @end deftypefn 25 26function [RT,xp] = radon (I,theta) 27 28 ## Input checking 29 if (nargin == 0 || nargin > 2) 30 print_usage (); 31 elseif (nargin == 1) 32 theta = 0:179; 33 endif 34 35 if (! isnumeric (I) || ndims (I) != 2) 36 error("radon: I must be a MxN numeric matrix"); 37 endif 38 39 if (!isvector(theta)) 40 error("radon: second input must be a vector"); 41 endif 42 43 [m, n] = size (I); 44 45 # center of image 46 xc = floor ((m+1)/2); 47 yc = floor ((n+1)/2); 48 49 # divide each pixel into 2x2 subpixels 50 51 d = reshape (I,[1 m 1 n]); 52 d = d([1 1],:,[1 1],:); 53 d = reshape (d,[2*m 2*n])/4; 54 55 b = ceil (sqrt (sum (size (I).^2))/2 + 1); 56 xp = [-b:b]'; 57 sz = size(xp); 58 59 [X,Y] = ndgrid (0.75 - xc + [0:2*m-1]/2,0.75 - yc + [0:2*n-1]/2); 60 61 X = X(:)'; 62 Y = Y(:)'; 63 d = d(:)'; 64 65 th = theta*pi/180; 66 67 for l=1:length (theta) 68 # project each pixel to vector (-sin(th),cos(th)) 69 Xp = -sin (th(l)) * X + cos (th(l)) * Y; 70 71 ip = Xp + b + 1; 72 73 k = floor (ip); 74 frac = ip-k; 75 76 RT(:,l) = accumarray (k',d .* (1-frac),sz) + accumarray (k'+1,d .* frac,sz); 77 endfor 78 79endfunction 80 81%!test 82%! A = radon(ones(2,2),30); 83%! assert (A,[0 0 0.608253175473055 2.103325780167649 1.236538105676658 0.051882938682637 0]',1e-10) 84