Simple v filtered back projection 

Matlab code

% Program to plot development of a sinogram, simple back projection, and 

% spatially-filtered back projection of a cartoon fly 

 

clear

close all

 

vid = VideoWriter('SimpleVfilteredBP.mp4','MPEG-4');

vid.Quality = 100;

vid.FrameRate = 60;

open(vid);

 

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');

set(0,'defaultfigurecolor',[1 1 1]); 

colormap bone

 

phantomData = imread('cartoonFly.png');

phantomData = double(phantomData - min(phantomData(:))); % set min of image to zero

% pad the image with zeros so we don't lose anything when we rotate.

padDims = ceil(norm(size(phantomData)) - size(phantomData)) + 2;

P       = padarray(phantomData,padDims);

 

pos1 = [0.23 0.5 0.27 0.45]; % Original

pos2 = [0.5 0.5 0.27 0.45];  % Sinogram

pos3 = [0.23 0.0 0.27 0.45]; % Simple back projection

pos4 = [0.5 0.0 0.27 0.45];  % Filtered back projection

 

%__________________________________________________________________________

 

% Plot original image

subplot('position',pos1)

imagesc(P);

title('Original','FontSize',20)

axis equal

axis off

% Set angular step size

freq = 2; % [1/degree]

thetas = 0:1/freq:180-1/freq;

 

%__________________________________________________________________________

 

subplot('position',pos2)

% Compute sinogram

numOfAngularProjections  = length(thetas);

numOfParallelProjections = size(P,1);

sinogram = zeros(numOfParallelProjections,numOfAngularProjections);

% Loop over the number of angles

for i = 1:length(thetas)

    % Rotate image

    tmpImage = imrotate(P,-thetas(i),'bilinear','crop');

    

    % Fill sinogram

    sinogram(:,i) = sum(tmpImage,2);

    imagesc(sinogram);

    title('Sinogram','FontSize',20)

    axis square

    axis off

    drawnow

    

    frame = getframe(gcf);

    writeVideo(vid,frame);

end

 

%__________________________________________________________________________

 

subplot('position',pos3)

% Simple backprojection

 

numOfParallelProjections = size(sinogram,1);

numOfAngularProjections  = length(thetas);

% Convert thetas to radians

thetas = (pi/180)*thetas;

% Set up the backprojected image

BPI = zeros(numOfParallelProjections,numOfParallelProjections);

% Find the middle index of the projections

midindex = floor(numOfParallelProjections/2) + 1;

% Set up the coords of the image

[xCoords,yCoords] = meshgrid(ceil(-numOfParallelProjections/2):ceil(numOfParallelProjections/2-1));

% Loop over each projection angle

for i = 1:numOfAngularProjections

    % Figure out which projections to add to which spots

    rotCoords = round(midindex + xCoords*sin(thetas(i)) + yCoords*cos(thetas(i)));

    % Check which coords are in bounds

    indices   = find((rotCoords > 0) & (rotCoords <= numOfParallelProjections));

    newCoords = rotCoords(indices);

    

    % Summation

    BPI(indices) = BPI(indices) + sinogram(newCoords,i)./numOfAngularProjections;

    imagesc(BPI)

    title('Simple backprojection','FontSize',20)

    axis equal

    axis off

    drawnow

    

    frame = getframe(gcf);

    writeVideo(vid,frame);    

end

 

%__________________________________________________________________________

 

subplot('position',pos4)

% Reconstruction with filtered back projection in spatial domain

 

BPI = zeros(numOfParallelProjections,numOfParallelProjections);

% Find the middle index of the projections

midindex = floor(numOfParallelProjections/2) + 1;

% Set up the coords of the image

[xCoords,yCoords] = meshgrid(ceil(-numOfParallelProjections/2):ceil(numOfParallelProjections/2-1));

% Set up filter: now for the spatial domain!!!

if mod(numOfParallelProjections,2) == 0

    halfFilterSize = floor(1 + numOfParallelProjections);

else

    halfFilterSize = floor(numOfParallelProjections);

end

filter = zeros(1,halfFilterSize);

filter(1:2:halfFilterSize) = -1./((1:2:halfFilterSize).^2 * pi^2);

filter = [fliplr(filter) 1/4 filter]; 

 

% Loop over each projection

for i = 1:numOfAngularProjections

    % Figure out which projections to add to which spots

    rotCoords = round(midindex + xCoords*sin(thetas(i)) + yCoords*cos(thetas(i)));

    % Check which coords are in bounds

    indices   = find((rotCoords > 0) & (rotCoords <= numOfParallelProjections));

    newCoords = rotCoords(indices);

    

    % Filter

    filteredProfile = conv(sinogram(:,i),filter,'same');

    % Summation

    BPI(indices) = BPI(indices) + filteredProfile(newCoords)./numOfAngularProjections;

    imagesc(BPI)

    title('Ram-Lak filtered backprojection','FontSize',20)

    axis square

    axis off

    drawnow

    

    frame = getframe(gcf);

    writeVideo(vid,frame);

end

close(vid);