Upload New File

fbp.m

+function [f,xf] = fbp(R,theta,omega,rho,b)
+%
+%   Inverse Radon transform (filtered backprojection).
+%
+%   on input:   R      Radon transform, 
+%                      each column corresponds to one angle theta
+%               theta  *all* angles at which data are taken
+%               omega  samples of the "s"-variable in Rf(theta,s)
+%               rho    region of interest
+%               b      banwidth of ramp filter
+%
+%   on output:  f      reconstructed function
+%               xf     x (and y) arguments of the corresponding columns/rows
+
+m = 128;   % pixel per row/column
+
+n = size(R,1);
+N = max(64,2^nextpow2(2*n));
+
+% ramp-Filter der Laenge N mit Bandbreite b
+% b = 3;
+k = 2*pi/(N*(omega(2)-omega(1)));
+nb = ceil(b/k);
+
+filter = 2*(0:N/2)/N;
+filter(nb:end) = 0;
+filter = [filter';filter(end-1:-1:2)'];
+filter = filter*ones(1,length(theta)); 
+
+
+
+R(N,1) = 0;  % fuellt array R durch Nullen auf zur Dimension N x length(theta)
+
+% Faltung:
+R = fft(R);
+R = R.*filter;
+R = ifft(R);
+
+R(n+1:end,:) = [];   % Rest wieder abschneiden
+
+clear filter
+
+% nun square of interest fuer backprojection:
+
+if omega(end)/sqrt(2) < rho 
+  error('region of interest is too big')
+end
+
+xf = (2*(0:m-1)/(m-1)-1)*rho;
+x = ones(m,1)*xf;
+y = x';
+
+% nun backprojection
+
+s = omega;
+
+cs = cos(theta);
+sn = sin(theta);
+
+f = zeros(m,m);
+
+for j = 1:length(theta)
+
+  p = R(:,j);
+  r = x*cs(j) + y*sn(j);
+  flocal = interp1(s,p,r(:),'linear');
+  flocal = reshape(flocal,m,m);
+
+  f = f + flocal;
+
+end
+
+f = f*pi/2/length(theta); 
+
+