| |
ELG 7173 - Backprojection
The following function makeproj calculates
an interpolating linear projection matrix for an
angle a
function proj= makeproj( a, x, y)
if ~all(diff([size(x),size(y)])==0);
error('x and y must be square');
end
rmax= max([abs(x(:));abs(y(:))]);
spc = max(abs([mean(mean(diff(x'))), mean(mean(diff(y ))) ]));
plen= size(x,1);
%create x indices of projection ray
xx=x*sin(a) + y*cos(a);
xidx= (xx + rmax) / spc + 1;
xi_l= floor(xidx); xi_h = xi_l+1;
xi_il= (xi_h-xidx); xi_ih= (xidx-xi_l);
%create y indices of projection ray
yy=x*cos(a) - y*sin(a);
yidx= (yy + rmax) / spc + 1;
yi_l= floor(yidx); yi_h = yi_l+1;
yi_il= (yi_h-yidx); yi_ih= (yidx-yi_l);
% keep elements within bounds
kp = (xx.^2 + yy.^2) < (rmax - spc);
pnum = ones(plen,1) * (1:plen); pnum = pnum(kp);
% create sparse interpolation matrix
posn = [ yi_l(kp), yi_h(kp), yi_l(kp), yi_h(kp)] + ...
plen*([ xi_l(kp), xi_l(kp), xi_h(kp), xi_h(kp)]-1);
aprx = [xi_il(kp).*yi_il(kp), xi_il(kp).*yi_ih(kp), ...
xi_ih(kp).*yi_il(kp), xi_ih(kp).*yi_ih(kp)];
proj = sparse( [pnum,pnum,pnum,pnum], posn, aprx,plen,plen^2);
% Generate a sample image
spc=.025; rlim=1;
[x,y]= meshgrid(-rlim:spc:rlim,-rlim:spc:rlim);
plen= size(x,1);
img = (x.^2 + y.^2) > rlim;
img( x >.45 & x<.65 & y>-.05 & y<.45) =1;
img( x >-.55 & x<-.25 & y>.45 & y<.65) =1;
imagesc(img);
|
|
% shape of projection matrix for a 45° projection angle
pr45= makeproj(45*(pi/180),x,y);
imagesc(full(pr45));
|
|
% calculate projections proj sample projections
i=1; proj= zeros(plen,6);
for ang= [0:30:150];
prm= makeproj(ang*(pi/180),x,y);
proj(:,i) = prm* img(:);
i=i+1;
end
plproj=proj/max(proj(:)) + ones(plen,1)*[0:5];
plot(x(1,:),plproj);
|
|
% backprojection reconstruction imbp
imbp= zeros(size(img));
i=1;
for ang= [0:30:150];
prm= makeproj(ang*(pi/180),x,y);
pp= proj(:,i);
imbp= imbp + reshape( prm' * proj(:,i) , plen, plen);
i=i+1;
end
imbp = imbp + max(imbp(:))*( x.^2 + y.^2 > rlim );
imagesc(imbp);
|
Reconstruction with 6 (left) and 60 (right) projections
|
% filtered backprojection reconstruction imbp
imbp= zeros(size(img));
% Calculate convolution filter xc (assign #2)
i=1; for ang= [0:30:150];
prm= makeproj(ang*(pi/180),x,y);
pp= conv2( proj(:,i), xc', 'same');
imbp= imbp + reshape( prm' * pp , plen, plen);
i=i+1;
end
imbp = imbp + max(imbp(:))*( x.^2 + y.^2 > rlim );
imagesc(imbp);
|
Reconstruction with 6 (left) and 60 (right) projections
|
Last Updated:
$Date: 2005-03-21 11:57:07 -0500 (Mon, 21 Mar 2005) $
|
Engineering