PDF (letter size)

## Matlab implementation to illustrate central slice theorem and back projection using inverse radon transform

Nov 10,2008   Compiled on September 7, 2023 at 11:58pm

### 1 Introduction

Central slice theorem says that if we make a projection of a 2D image on a projection line, and take the 1D Fourier transform (say A) of the projection itself, and then take a slice (say B) from the 2D Fourier transform of the image itself, then A=B. When taking the slice from the 2D Fourier transform it has to be done on a slice through the center of the 2D spectrum, and along a line parallel to the projection line mentioned above.

This is a simple Matlab implementation to illustrate the above. It allows the user to select a number of images, and the angle of the projection and then it ﬁnd the projection (using radon transform) and shows the 1D FT of this projection, and then it ﬁnds the 2D FT of the image itself.

The 2D FT is shown in a 3D view to allow the user to rotate it.

This program also illustrate back projection by allowing the user to select the number of degrees to project the image at, and then back project all the resulting projecting to reconstruct the original image (as is done in CT scans). We see that the more angles used, the closer the resulting image will be to the original image.

The 2D FT of the reconstructed image is also shown as number of degrees is increased allowing one to see the slices being added and the ﬁnal 2D FT convergence to the original image 2D FT.

Streaks in the image are noticed as it is being constructed. The streaks become less noticeable as more angles added. The inverse radon transform is used for back projection.

### 2 Animation

Here is a screen shot of the Matlab program.

### 3 source code

Current version is here zip ﬁle. This is a zip ﬁle which contains all the ﬁles needed. Download the zip ﬁle and extract it. It will created a folder. Simply add this folder to your MATLAB PATH, then from the Matlab console type the command

     nma_projection

This will bring up the GUI.

### 4 Soruce code listing

#### 4.1 nma_projection.m

function varargout = nma_projection(varargin)
%
%  illustrate central slice theorem and reconstruction
%  of 2D images from backprojections using inverse radon transform
%
%  work related to Math 597 project, summer 2008
%  by Nasser Abbasi 6/2/08
%

% NMA_PROJECTION M-file for nma_projection.fig
%      NMA_PROJECTION, by itself, creates a new NMA_PROJECTION or raises the existing
%      singleton*.
%
%      H = NMA_PROJECTION returns the handle to a new NMA_PROJECTION or the handle to
%      the existing singleton*.
%
%      NMA_PROJECTION('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in NMA_PROJECTION.M with the given input arguments.
%
%      NMA_PROJECTION('Property','Value',...) creates a new NMA_PROJECTION or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before nma_projection_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to nma_projection_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%

% Edit the above text to modify the response to help nma_projection

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
'gui_Singleton',  gui_Singleton, ...
'gui_OpeningFcn', @nma_projection_OpeningFcn, ...
'gui_OutputFcn',  @nma_projection_OutputFcn, ...
'gui_LayoutFcn',  [] , ...
'gui_Callback',   []);
if nargin && ischar(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT
end

% --- Executes just before nma_projection is made visible.
function nma_projection_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to nma_projection (see VARARGIN)

% Choose default command line output for nma_projection
handles.output = hObject;

set(handles.groupBtnTag,'SelectionChangeFcn',...
@groupBtnTag_SelectionChangeFcn);
set(handles.colormapTag,'SelectionChangeFcn',...
@colormapTag_SelectionChangeFcn);

% Update handles structure
set(handles.figure1,'Name',...
'Computed tomography simulation by Nasser Abbasi, CSUF. EE 518, Digital Signal Processing');
guidata(hObject, handles);

process(hObject,0);

end

% --- Outputs from this function are returned to the command line.
function varargout = nma_projection_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;
end

%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%
function groupBtnTag_SelectionChangeFcn(hObject, eventdata)
process(hObject,0);
end

%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%
function colormapTag_SelectionChangeFcn(hObject, eventdata)
handles   = guidata(hObject);
h         = get(handles.colormapTag);
hColorMap = h.SelectedObject;
switch get(hColorMap,'Tag')
case 'jetTag'
colormap(jet);
case 'grayTag'
colormap(gray);
end

end

%%%%%%%%%%%
%
%%%%%%%%%%%
function process(hObject,doSimulation)

%retrieve GUI data, i.e. the handles structure
handles = guidata(hObject);

nPixels  = 600;
c        = 22
clip     = 10;

h         = get(handles.groupBtnTag);
hselected = h.SelectedObject;

%Now obtain colormap
h          = get(handles.colormapTag);
hColorMap  = h.SelectedObject;
switch get(hColorMap,'Tag')
case 'jetTag'
colormap(jet);
case 'grayTag'
colormap(gray);
end

%Now obtain the angle of projection and display projection
h     = get(handles.angleTag);
angle = str2double(h.String);

switch get(hselected,'Tag')   % Get Tag of selected object
case 'blackDiskTag'

case 'verticalBarTag'
[A,p] = makeVerticalBar(nPixels);

case 'randomDiskTag'

case 'lenaTag'
fileName = 'lena.jpg';

case 'lungTag'
fileName = 'lung.jpg';
[nRow,nCol] = size(A);
d = min(nRow,nCol);
A = imresize(A, [d d]);

case 'blobsTag'
fileName     = 'blobs.gif';
[nRow,nCol]  = size(A);
d            = min(nRow,nCol);
A            = imresize(A, [d d]);

otherwise
% Code for when there is no match.
end

axes(handles.originalImageAxes);
imagesc(A);
axis(handles.originalImageAxes,'image');

%now display projection.
axes(handles.projectionAxes);
stairs(p);

%now make fft of projection
axes(handles.normalProjectionTransformAxes);
YProjection = fft(p);
YshiftedProjection = fftshift(YProjection);
plot(abs(YshiftedProjection));

%now make 2D fft of original image
Y2D = fft2(double(A));
Y2Dshifted = abs(fftshift(Y2D));
axes(handles.twoDTransformAxes);
imagesc(c*log(1+Y2Dshifted));

xlabel('v'); ylabel('u');
axis(handles.twoDTransformAxes,'image');

%
%display the 2D spectrum on 3D
make3dspectrum(handles.FFT2on3DAxes,abs(Y2Dshifted),clip,c);

if doSimulation
backProjection(hObject,A,clip,c);
end

end
%%%%%%%%%%%
%
%%%%%%%%%%%
function disk=makeDisk(isRandom,nPixels,r)
BLACK =0;
WHITE =255;

disk    = zeros(nPixels);
xCenter = ceil(nPixels/2);
yCenter = ceil(nPixels/2);

for i=1:nPixels
for j = 1:nPixels
xReal = i-xCenter;
yReal = j-yCenter;
distance = sqrt(xReal^2+yReal^2);
if distance>r
disk(i,j) = BLACK;
else
if isRandom
disk(i,j) = WHITE*rand;
else
disk(i,j) = WHITE;
end
end
end
end
end

%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%
function angleTag_Callback(hObject, eventdata, handles)
% hObject    handle to angleTag (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of angleTag as text
%        str2double(get(hObject,'String')) returns contents of angleTag as a double

process(hObject,0);
end

%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%

% --- Executes during object creation, after setting all properties.
function angleTag_CreateFcn(hObject, eventdata, handles)
% hObject    handle to angleTag (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
end

%%%%%%%%%%
%
%%%%%%%%%%
function [A,p]=makeVerticalBar(nPixels)
A=[0 0 1 0 0;
0 0 1 0 0;
0 0 1 0 0;
0 0 1 0 0;
0 0 1 0 0];

p=[1;
1;
1;
1;
1];
end

%%%%%%%%%%%%
%
%%%%%%%%%%%%%
function nAnglesTag_Callback(hObject, eventdata, handles)
% hObject    handle to nAnglesTag (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of nAnglesTag as text
%        str2double(get(hObject,'String')) returns contents of nAnglesTag as a double

end

% --- Executes during object creation, after setting all properties.
function nAnglesTag_CreateFcn(hObject, eventdata, handles)
% hObject    handle to nAnglesTag (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
end

% --- Executes on button press in BtnTag.
function BtnTag_Callback(hObject, eventdata, handles)
% hObject    handle to BtnTag (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%Now obtain how many angles to use for backprojection
process(hObject,1);
end
%%%%%%%%%%
%
%%%%%%%%%%%
function backProjection(hObject,A,clip,c)

%retrieve GUI data, i.e. the handles structure
handles = guidata(hObject);

%Now obtain how many angles to use for backprojection
h       = get(handles.nAnglesTag);
nAngles = str2double(h.String);

angles = zeros(nAngles,1);
delta=180/(nAngles+1)
for i=1:nAngles
angles(i)=delta*i;
end

for i = 1:nAngles
imagesc(I);

xlabel(sprintf('number of angles [%d]',i));

make3dspectrum(handles.FFT2on3DcurrentAxes,abs(fftshift(fft2(I))),clip,c)
end

end

% --- Executes on selection change in iradonInterpTag.
% hObject    handle to iradonInterpTag (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns iradonInterpTag contents as cell array
%        contents{get(hObject,'Value')} returns selected item from iradonInterpTag
end

% --- Executes during object creation, after setting all properties.
% hObject    handle to iradonInterpTag (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
end

% --- Executes on selection change in iradonFilterTag.
% hObject    handle to iradonFilterTag (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns iradonFilterTag contents as cell array
%        contents{get(hObject,'Value')} returns selected item from iradonFilterTag
end

% --- Executes during object creation, after setting all properties.
% hObject    handle to iradonFilterTag (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
end


#### 4.2 nma_build.m

function nma_build_HW5()

list = dir('*.m');

if isempty(list)
fprintf('no matlab files found\n');
return
end

for i=1:length(list)
name=list(i).name;
fprintf('processing %s\n',name)
p0 = fdep(list(i).name,'-q');
[pathstr, name_of_matlab_function, ext] = fileparts(name);

%make a zip file of the m file and any of its dependency
p1=dir([name_of_matlab_function '.fig']);
if length(p1)==1
files_to_zip =[p1(1).name;p0.fun];
else
files_to_zip =p0.fun;
end

zip([name_of_matlab_function '.zip'],files_to_zip)

end

end


#### 4.3 make3dspectrum.m

%****************************************************
%       Copyright (C) 2008 Nasser Abbasi
%  Free to use and modify for academic and research only
%***************************************************

function make3dspectrum(axesHandle,F,clip,c)

axes(axesHandle);
L      = size(F,1);
[i,j]  = find(F>(clip/100)*(max(max(F))));
d      = F;
d(sub2ind(size(F),i,j)) = 0;
r      = round(.15*L);
extent = r:L-r;
mesh(extent,extent,c*log(1+d(extent,extent)));
axis('tight');
xlabel('u');
ylabel('v');
view(-45,50);
end


%****************************************************
%       Copyright (C) 2008 Nasser Abbasi
%  Free to use and modify for academic and research only
%***************************************************

switch k
case 1
the_filter = 'Ram-Lak';
case 2
the_filter = 'Shepp-Logan';
case 3
the_filter =  'Cosine' ;
case 4
the_filter = 'Hamming' ;
case 5
the_filter = 'Hann' ;
case 6
the_filter = 'None';
end