滤波反投影程序设计报告

《滤波反投影程序设计报告》

课程名称：生物医学图像处理2

院系：生物医学工程

姓名：

学号：

完成日期：2017年4月23日

一、设计目的

用Matlab实现平行束滤波反投影算法，比较不同滤波函数的效果。

二、实验原理

（一）图像重建模型——shepp Logan头模型

是图像重建标准体模，由10个位置、大小、方向、密度各异的椭圆组成，代表一个脑部断层。

（二）重建理论推导

中心切片定理是从投影图像重建图像的理论基础，表述为：某断层图像f(x,y)在视角为θ时得到的平行投影的一维傅里叶变换等于f(x,y)二维傅里叶变换F（w1,w2）过原点的一个垂直切片，且切片与轴w1相交成θ角。如下图所示。

公式表述为：F （wcos θ,wsin θ）=P(w,θ) ① 将在ω1-ω2坐标系中表达的F （w 1,w 2）引入新的极坐标系ω−θ中，新坐标系与原坐标系原点重合，有w1=wcos θ,w2=wsin θ. 面元换算为dw1dw2=wdwd θ.

有 f(x,y)= ∫∫F （wcosθ，wsinθ）e j2πw(xcosθ+ysinθ)wdwdθ∞

02π

0 =∫∫P(w,θ)e j2πw(xcosθ+ysinθ)wdwdθ∞02π

0 =∫∫P(w,θ)e j2πw(xcosθ+ysinθ)wdwdθ∞0π

0 +

∫∫P(w,θ+π)e −j2πw(xcosθ+ysinθ)wdwdθ∞0π

0 ② 注意到在θ与θ+π两个角度下的平行束射线投影的投影存在如下关系：

p (t,θ)=p (−t,θ+π)

其傅里叶变换存在如下关系：P(w,θ+π)=P(−w,θ)

将上式代入②式，有 f(x,y)= ∫∫P (w,θ)|w|e j2πw(xcosθ+ysinθ)dwdθ∞

−∞π

0 ③ 令③式内积分等于g(xcos θ+ysin θ),则有

g(xcos θ+ysin θ)=∫|w |P (w,θ)e j2πwt |+∞−∞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)= ∫g(x π

0cosθ+ysinθ)d θ ④ 在图3.5中注意到

Xr=rcos(θ−φ)=xcos θ+ysinθ

是从原点出发的通过点（r,θ）的射线方程，④式可写为：

f(x,y)= ∫g[rcos (θ−φ),φ]π0dφ ⑤ ④⑤两式表明：f(x,y)在（x,y ）处的重建，等于通过该点的所有角度下滤波投影的累加所得到的像素值，而Xr=rcos(θ−φ)=xcos θ+ysinθ的变化代表了所有平行投影射线。

（三）Radon 变换

一个无限薄的切片内相对线性衰减系数的分布是由它的所有线积分的集合唯一决定，揭示了函数和投影之间的关系，若函数为f(x,y),则不同角度下的投影可写为

P(t,θ)=∬f (x,y )δ(xcosθ+ysinθ−t )dxdy +∞

−∞ ⑥

（四）滤波函数

由于直接反投影法把取自有限空间的投影均匀回抹到了射线所及的无限空间的各个像素上，使得原来像素值为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自带函数，附录中对过程有详细注释，程序大致流程为:

四、运行结果

（1）sheep Logan256*256图像重建

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

（3）自定义RGB图像重建

（

滤波函数进行滤波处理。

五、总结

本次程序设计完成了设计目的，投影和反投影过程没有使用自带函数，但是在处理像素比较大的图像时程运行时间比较长，是一个缺点，图像也存在一定误差，效果不是特别完美，不过也算是锻炼了自己的能力，对图像重建中对滤波反投影算法有了深刻的了解，增加了对这门课程的兴趣，期待在以后的学习中能进一步完善程序，缩短运行时间，减小图像误差。

六、附录——程序代码

①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('自定义投影后图像');