滤波反投影程序设计报告.doc

《滤波反投影程序设计报告》
课程名称： 生物医学图像处理2
院 系： 生物医学工程
姓 名：
学 号：
完成日期： 2017年4月23日
一、设计目的
用Matlab实现平行束滤波反投影算法，比较不同滤波函数的效果。
二、实验原理
（一）图像重建模型——shepp Logan头模型

是图像重建标准体模，由10个位置、大小、方向、密度各异的椭圆组成，代表一个脑部断层。
（二）重建理论推导
中心切片定理是从投影图像重建图像的理论基础，表述为：某断层图像f(x,y)在视角为θ时得到的平行投影的一维傅里叶变换等于f(x,y)二维傅里叶变换F（w1,w2）过原点的一个垂直切片，且切片与轴w1相交成θ角。如下图所示。

公式表述为：F（wcos,wsin）=P(w,) ①
将在-坐标系中表达的F（w1,w2）引入新的极坐标系中，新坐标系与原坐标系原点重合，有w1=wcos,w2=wsin.
面元换算为dw1dw2=wdwd.
有 f(x,y)=
=
= +
②
注意到在

其傅里叶变换存在如下关系：P(w,
将上式代入②式，有 f(x,y)= ③
令③式内积分等于g(xcos+ysin),则有
g(xcos+ysin)=t= xcosθ+ysinθdw
实际上，g(xcos+ysin)就是投射角度为时的滤波投影，在t-s坐标系中表达时，转化为g（t,）=h(t)*p(t,),h(t)是传递函数H（w）=|w|的傅里叶逆变换，p(t,)是P（w,）的傅里叶逆变换。所以③式可写成
f(x,y)=θ ④
在图3.5中注意到
Xr=rcos()=xcos
是从原点出发的通过点（r,）的射线方程，④式可写为：
f(x,y)=
④⑤两式表明：f(x,y)在（x,y）处的重建，等于通过该点的所有角度下滤波投影的累加所得到的像素值，而Xr=rcos()=xcos的变化代表了所有平行投影射线。
（三）Radon变换
一个无限薄的切片内相对线性衰减系数的分布是由它的所有线积分的集合唯一决定，揭示了函数和投影之间的关系，若函数为f(x,y),则不同角度下的投影可写为
P(t,)= ⑥
（四）滤波函数
由于直接反投影法把取自有限空间的投影均匀回抹到了射线所及的无限空间的各个像素上，使得原来像素值为0的点不为0，从而产生星状伪迹，滤波反投影算法用人为设计的一维滤波函数对所得投影数据进行卷积，而后进行反投影和累加时，由于正负抵消，可一定程度上消除星状伪迹。
现在常用的滤波函数有R-L、S-L、Hanning、Hamming、Parzen、Sigwin窗函数，本次设计比较了R-L、S-L和Parzen 滤波函数的效果，发现Parzen滤波效果最好。
1.R-L滤波函数
R-L滤波函数频率响应为：

R-L滤波函数滤波计算简单，避免了大量的正弦、余弦计算，所重建图像轮廓清楚，空间分辨率高，但有Gibbs现象，表现为明显的振荡响应，特别是工件的边缘和衰减系数变化剧烈（即密度变化）时，尤为明显。

S-L滤波器对投影中的高频成分具有抑制作用，进而使重建图像的振荡响应减小，对含噪声的投影数据，重建质量较用RL的好，但在低频段不及R-l滤波函数好，这是由于H（w）在低频段也偏离了|w|的缘故
3.Parzen滤波函数

三、程序实现
程序包含FBP_GUI.m和dps.m两个m文件，其中投影、滤波和反投影过程均为自己编写，没有使用Matlab自带函数，附录中对过程有详细注释，程序大致流程为:
运行FBP_GUI.m，跳出界面，选择图像，进行相应处理后存入P，传入主函数dps().

投影角度和步长定义为1，利用Radon变换原理进行投影存入R，截取正弦图R1

构造滤波函数R-L、S-L和Parzen的离散采样序列存入h

将滤波函数h与R1卷积存入G

线性内插，将滤波投影值回抹存入a

计算三个图像重建精度判据d,r,e，文本输出

重建结束、点击别的按钮可进行新图像的重建、点击exit则关闭整个GUI界面。

四、运行结果
1）sheep Logan256*256图像重建

（2）自定义灰度图像重建

（3）自定义RGB图像重建

（4）不同滤波函数重建效果比较（以sheep logan图像为例）
使用的滤波函数
归一化均方距离判据d
归一化平均绝对距离判据r
最坏情况距离判e
Parzen
0.50487
0.7073
0.48133
S-L
0.63614
0.99191
0.49659
R-L
0.75607
1.2157
0.50068
可以看出不同滤波函数重建精度：Parzen优于S-L优于R-L。所以后续图像均采用Parzen滤波函数进行滤波处理。
五总结
本次程序设计完成了设计目的，投影和反投影过程没有使用自带函数，但是在处理像素比较大的图像时程运行时间比较长，是一个缺点，图像也存在一定误差，效果不是特别完美，不过也算是锻炼了自己的能力，对图像重建中对滤波反投影算法有了深刻的了解，增加了对这门课程的兴趣，期待在以后的学习中能进一步完善程序，缩短运行时间，减小图像误差。
六、附录——程序代码
①dps.m文件代码如下：
function [a]=dps(P)
tic;
%P=phantom(256);
%P=imread('gray1.jpg');
%P=rgb2gray(P);
[N,N]=size(P);
subplot(2,3,2);
imshow(P);
title([int2str(N),'*',int2str(N),'原始图像']);
%先进行自定义radon变换------------------------------------------------------------
thm=45; %45度时会出现最大尺寸
pre = imrotate(P,thm);
[mmax,nmax] = size(pre);
s=1;
%创建一个180*nmax的空白图片，用以存储投影后的线状图片
Final = zeros(180/s,nmax);%这里180代表180角度，每个角度投影成为一条线
t = 1;
for theta = 1:s:179
Protate = imrotate(P,theta); %对原图旋转一个角度，求和（线积分）
Pf = sum(Protate,1);
[mreal,nreal]=size(Pf); %计算实际尺寸
%确定起始点
if (nmax - nreal)/2-floor((nmax - nreal)/2) == 0
From = floor((nmax - nreal)/2 + 1);%总点数为偶数时
else
From = floor((nmax - nreal)/2) + 1;%总点数为奇数时
end
%确定结束点
End = floor((nmax-nreal)/2) + nreal;
%将一个角度Radon变换后的线状图存入结果图像的某一行
Final(180/s-t,From:End) = Pf; %从最底下一行开始存起
%上移一行
t = t + 1;
end
%再逆时针旋转
R=imrotate(Final,90);
subplot(2,3,3);
imshow(R,[]);
title('自定义投影后图像');
z=2*ceil(norm(size(P)-floor((size(P)-1)/2)-1))+3;% radon变换默认平移点数/角度
e=floor((z-N)/2)+2;
R1=R(e:(N+e-1),:);
[mm,nn]=size(R1);
subplot(2,3,4);
imagesc(R1);
title([int2str(mm),'*',int2str(nn),'正弦图']);
colormap(gray);
colorbar;
%滤波函数构造------------------------------------------------------------
g=1-N:N-1;
% d1;
% R-L滤波函数
% for i=1:2*N-1
% if g(i)==0
% h(i)=1/(4*(d^2));
% else if mod(g(i),2)==0
% h(i)=0;
% else
% h(i)=(-1)/(pi^2*d^2*(g(i)^2));
% end
% end
% end
%S-L滤波函数
% d=1;
% for i=1:2*N-1
% h(i)=2/(pi^2*d^2*(14*g(i)^2));
% end
%Parzen滤波函数
for i=1:2*N-1
h(i)=(24*pi*g(i)*cos(pi*g(i))-96*sin(pi*g(i))-48*pi*g(i)*cos(pi*g(i)/2)+384*sin(pi*g(i)/2)-2*(pi^3)*(g(i)^3)-72*pi*g(i))/(4*(pi^5)*(g(i)^5));
end
h(N)=0.0435;
%将投影与滤波函数卷积----------------------------------------------------
G=zeros(N,180);
for m=1:180
for i=1:N
b=0;
for n=1:N
b=b+h(N+n-i)*R1(n,m);
G(i,m)=b;
end
end
end
%投影滤波后线性内插进行图像重建----------------------------------------------
a=zeros(N); %重建图像初始化，每个像素点像素值为0
ns=(N+1)/2;
for m=1:180 %遍历每个投影角度

r=pi/180; %将角度换为弧度
for i=1:N
for j=1:N %遍历原图的每一个像素点
Xrm=(i-ns)*cos(m*r)+(j-ns)*sin(m*r)+ns-1; %坐标转换
if Xrm<1
n=1; %内插运算整数值
t=abs(Xrm)-floor(abs(Xrm)); %内插运算小数值
else
n=floor(Xrm);
t=Xrm-floor(Xrm);
end
if n>N-1
n=N-1;
end
c=(1-t)*G(n,m)+t*G(n+1,m); %内插后的滤波投影值
a(N+1-j,i)=a(N+1-j,i)+c; %像素值的叠加
end
end
end
%将P、a的像素值变换到0-1之间
P=double(P);
P=P-min(min(P));
kk=max(max(P));
for i=1:N
for j=1:N
P(i,j)=P(i,j)/kk;
end
end
a=double(a);
a=a-min(min(a));
kkk=max(max(a));
for i=1:N
for j=1:N
a(i,j)=a(i,j)/kkk;
end
end
subplot(2,3,5);
imshow(a,[]); %重建图形的显示
title('滤波反投影重建图像');
%重建图像质量评价--------------------------------------------------------
%计算归一化均方距离判据d;
ave=0;
for x=1:N
for y=1:N
ave=ave+P(x,y);
end
end
ave=ave/(N*N);
x1=0;
x2=0;
for x=1:N
for y=1:N
x1=x1+(P(x,y)-a(x,y))^2;
x2=x2+(P(x,y)-ave)^2;
end
end
evaluate_d=sqrt(double(x1/x2));
%计算归一化平均绝对距离判据 r;
x3=0;
x4=0;
for x=1:N
for y=1:N
x3=x3+abs(P(x,y)-a(x,y));
x4=x4+abs(P(x,y));
end
end
evaluate_r=x3/x4;
%计算最坏情况距离判据e
ns=floor(N/2)-1;
g=zeros(ns);
for x=1:ns
for y=1:ns
T=(P(2*x,2*y)+P(2*x+1,2*y)+P(2*x,2*y+1)+P(2*x+1,2*y+1))/4;
R=(a(2*x,2*y)+a(2*x+1,2*y)+a(2*x,2*y+1)+a(2*x+1,2*y+1))/4;
g(x,y)=abs(T-R);
end
end
evaluate_e=max(max(g));
%结果文本显示------------------------------------------------------------
o=ones(500,1000);
subplot(2,36);
imshow(o,[]);
s_title=['图像重建精度判据如下:'];
text(0,0,s_title,'Fontsize',14);
s=num2str(toc);
s_one=['run time = ' s ' s;'];
text(0,100,s_one,'FontSize',10);
s=num2str(evaluate_d);
s_two=['归一化均方距离判据d=' s ';'];
text(0,200,s_two,'Fontsize',10);
s=num2str(evaluate_r);
s_three=['归一化平均绝对距离判据r' s ';'];
text(0,300,s_three,'Fontsize',10);
s=num2str(evaluate_e);
s_four=['最坏情况距离判据e=' s ';'];
text(0,400,s_four,'Fontsize',10);
toc
②FBP_GUI.m文件代码如下：
function varargout = FBP_GUI(varargin)
% FBP_GUI MATLAB code for FBP_GUI.fig
% FBP_GUI, by itself, creates a new FBP_GUI or raises the existing
% singleton*.
%
% H = FBP_GUI returns the handle to a new FBP_GUI or the handle to
% the existing singleton*.
%
% FBP_GUI('CALLBACK',hObject,eventData,handles,...) calls the local
% function named CALLBACK in FBP_GUI.M with the given input arguments.
%
% FBP_GUI('Property','Value',...) creates a new FBP_GUI or raises the
% existing singleton*. Starting from the left, property value pairs are
% applied to the GUI before FBP_GUI_OpeningFcn gets called. An
% unrecognized property name or invalid value makes property application
% stop. All inputs are passed to FBP_GUI_OpeningFcn via varargin.
%
% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one
% instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES
% Edit the above text to modify the response to help FBP_GUI
% Last Modified by GUIDE v2.5 15-Apr-2017 14:24:03
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @FBP_GUI_OpeningFcn, ...
'gui_OutputFcn', @FBP_GUI_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
% --- Executes just before FBP_GUI is made visible.
function FBP_GUI_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 FBP_GUI (see VARARGIN)
% Choose default command line output for FBP_GUI
handles.output = hObject;
% Update handles structure
guidata(hObject, handles);
% UIWAIT makes FBP_GUI wait for user response (see UIRESUME)
% uiwait(handles.figure1);
% --- Outputs from this function are returned to the command line.
function varargout = FBP_GUI_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;
% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject handle to pushbutton1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
I=phantom(256);
[A]=dps(I);
% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject handle to pushbutton2 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
[fnpn,fi]=uigetfile('*.*','选择图片');
lujing=[pn fn];%得到待重建图像存储路径
I=imread(lujing);
[A]=dps(I);
% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject handle to pushbutton3 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
[fn,pn,fi]=uigetfile('*.*','选择图片');
lujing=[pn fn];%得到待重建图像存储路径
I=imread(lujing);
I=rgb2gray(I);
[A]=dps(I);
% --- Executes on button press in pushbutton4.
function pushbutton4_Callback(hObject, eventdata, handles)
% hObject handle to pushbutton4 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
close all
. . . .