## The Back Projection Code

```% backproject3.m

%% This is a MATLAB function that takes filtered back projections without
%% using the 'imrotate' command.  Here, we try to cut out all the loops.
%% PR is a matrix whose columns are the projections at each angle.
%% THETA is a row vector of the angles of the respective projections.
%%
%% Written by : Justin K. Romberg

function [BPI,M] = backproject3(PR, THETA)

% figure out how big our picture is going to be.
n = size(PR,1);
sideSize = n;

% filter the projections
filtPR = projfilter(PR);
%filtPR = filterplus(PR);
%filtPR = PR;

% convert THETA to radians
th = (pi/180)*THETA;

% set up the image
m = length(THETA);
BPI = zeros(sideSize,sideSize);

% find the middle index of the projections
midindex = (n+1)/2;

% set up x and y matrices
x = 1:sideSize;
y = 1:sideSize;
[X,Y] = meshgrid(x,y);
xpr = X - (sideSize+1)/2;
ypr = Y - (sideSize+1)/2;

% loop over each projection
%figure
%colormap(jet)
%M = moviein(m);
for i = 1:m
tic
disp(['On angle ', num2str(THETA(i))]);

% figure out which projections to add to which spots
filtIndex = round(midindex + xpr*sin(th(i)) - ypr*cos(th(i)));

% if we are "in bounds" then add the point
BPIa = zeros(sideSize,sideSize);
spota = find((filtIndex > 0) & (filtIndex <= n));
newfiltIndex = filtIndex(spota);
BPIa(spota) = filtPR(newfiltIndex(:),i);
%keyboard
BPI = BPI + BPIa;

toc

%imagesc(BPI)
%M(:,i) = getframe;
%figure(2)
%plot(filtPR(:,i));
%keyboard
end

BPI = BPI./m;
```
Back to the previous page.