滤波反投影程序设计报告Word格式.docx
《滤波反投影程序设计报告Word格式.docx》由会员分享，可在线阅读，更多相关《滤波反投影程序设计报告Word格式.docx（17页珍藏版）》请在冰豆网上搜索。
有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，文本输出
四、运行结果
（1）sheepLogan256*256图像重建
（2）自定义灰度图像重建
（3）自定义RGB图像重建
（4）不同滤波函数重建效果比较（以sheeplogan图像为例）
使用的滤波函数
归一化均方距离判据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;
fortheta=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;
%d=1;
%R-L滤波函数
%fori=1:
2*N-1
%ifg（i）==0
%h（i）=1/（4*（d^2））;
%elseifmod（g（i）,2）==0
%h（i）=0;
%else
%h（i）=（-1）/（pi^2*d^2*（g（i）^2））;
%end
%S-L滤波函数
%h（i）=2/（pi^2*d^2*（1-4*g（i）^2））;
%Parzen滤波函数
fori=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））;
h（N）=0.0435;
%将投影与滤波函数卷积----------------------------------------------------
G=zeros（N,180）;
form=1:
180
N
b=0;
forn=1:
b=b+h（N+n-i）*R1（n,m）;
G（i,m）=b;
%投影滤波后线性内插进行图像重建----------------------------------------------
a=zeros（N）;
%重建图像初始化，每个像素点像素值为0
ns=（N+1）/2;
180%遍历每个投影角度
r=pi/180;
%将角度换为弧度
N
forj=1:
N%遍历原图的每一个像素点
Xrm=（i-ns）*cos（m*r）+（j-ns）*sin（m*r）+ns-1;
%坐标转换
ifXrm<
1
n=1;
%内插运算整数值
t=abs（Xrm）-floor（abs（Xrm））;
%内插运算小数值
n=floor（Xrm）;
t=Xrm-floor（Xrm）;
ifn>
N-1
n=N-1;
c=（1-t）*G（n,m）+t*G（n+1,m）;
%内插后的滤波投影值
a（N+1-j,i）=a（N+1-j,i）+c;
%像素值的叠加
%将P、a的像素值变换到0-1之间
P=double（P）;
P=P-min（min（P））;
kk=max（max（P））;
fori=1:
P（i,j）=P（i,j）/kk;
a=double（a）;
a=a-min（min（a））;
kkk=max（max（a））;
a（i,j）=a（i,j）/kkk;
subplot（2,3,5）;
imshow（a,[]）;
%重建图形的显示
滤波反投影重建图像'
%重建图像质量评价--------------------------------------------------------
%计算归一化均方距离判据d;
ave=0;
forx=1:
for
升级会员 