滤波反投影程序设计报告

0《滤波反投影程序设计报告》课程名称：生物医学图像处理2院系：生物医学工程姓名：学号：完成日期：2017年4月23日1一、设计目的用Matlab实现平行束滤波反投影算法，比较不同滤波函数的效果。二、实验原理（一）图像重建模型——sheppLogan头模型是图像重建标准体模，由10个位置、大小、方向、密度各异的椭圆组成，代表一个脑部断层。（二）重建理论推导中心切片定理是从投影图像重建图像的理论基础，表述为：某断层图像f(x,y)在视角为θ时得到的平行投影的一维傅里叶变换等于f(x,y)二维傅里叶变换F（w1,w2）过原点的一个垂直切片，且切片与轴w1相交成θ角。如下图所示。2公式表述为：F（wcosθ,wsinθ）=P(w,θ)①将在𝜔1-𝜔2坐标系中表达的F（w1,w2）引入新的极坐标系ω−θ中，新坐标系与原坐标系原点重合，有w1=wcosθ,w2=wsinθ.面元换算为dw1dw2=wdwdθ.有f(x,y)=∫∫𝐹（𝑤𝑐𝑜𝑠𝜃，𝑤𝑠𝑖𝑛𝜃）𝑒𝑗2𝜋𝑤(𝑥𝑐𝑜𝑠𝜃+𝑦𝑠𝑖𝑛𝜃)𝑤𝑑𝑤𝑑𝜃∞02𝜋0=∫∫𝑃(𝑤,𝜃)𝑒𝑗2𝜋𝑤(𝑥𝑐𝑜𝑠𝜃+𝑦𝑠𝑖𝑛𝜃)𝑤𝑑𝑤𝑑𝜃∞02𝜋0=∫∫𝑃(𝑤,𝜃)𝑒𝑗2𝜋𝑤(𝑥𝑐𝑜𝑠𝜃+𝑦𝑠𝑖𝑛𝜃)𝑤𝑑𝑤𝑑𝜃∞0𝜋0+∫∫𝑃(𝑤,𝜃+𝜋)𝑒−𝑗2𝜋𝑤(𝑥𝑐𝑜𝑠𝜃+𝑦𝑠𝑖𝑛𝜃)𝑤𝑑𝑤𝑑𝜃∞0𝜋0②注意到在θ与θ+π两个角度下的平行束射线投影的投影存在如下关系：p(t,θ)=p(−t,θ+π)其傅里叶变换存在如下关系：P(w,θ+π)=P(−w,θ)将上式代入②式，有f(x,y)=∫∫𝑃(𝑤,𝜃)|𝑤|𝑒𝑗2𝜋𝑤(𝑥𝑐𝑜𝑠𝜃+𝑦𝑠𝑖𝑛𝜃)𝑑𝑤𝑑𝜃∞−∞𝜋0③令③式内积分等于g(xcosθ+ysinθ),则有g(xcosθ+ysinθ)=∫|𝑤|𝑃(𝑤,𝜃)𝑒𝑗2𝜋𝑤𝑡|+∞−∞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)=∫𝑔(𝑥𝜋0𝑐𝑜𝑠𝜃+𝑦𝑠𝑖𝑛𝜃)𝑑θ④在图3.5中注意到Xr=rcos(θ−φ)=xcosθ+ysinθ是从原点出发的通过点（r,θ）的射线方程，④式可写为：f(x,y)=∫𝑔[𝑟𝑐𝑜𝑠(𝜃−𝜑),𝜑]𝜋0𝑑φ⑤④⑤两式表明：f(x,y)在（x,y）处的重建，等于通过该点的所有角度下滤波投影的累加所得到的像素值，而Xr=rcos(θ−φ)=xcosθ+ysinθ的变化代表了所有平行投影射线。（三）Radon变换一个无限薄的切片内相对线性衰减系数的分布是由它的所有线积分的集合唯一决定，揭示了函数和投影之间的关系，若函数为f(x,y),则不同角度下的投影可写为P(t,θ)=∬𝑓(𝑥,𝑦)𝛿(𝑥𝑐𝑜𝑠𝜃+𝑦𝑠𝑖𝑛𝜃−𝑡)𝑑𝑥𝑑𝑦+∞−∞⑥（四）滤波函数由于直接反投影法把取自有限空间的投影均匀回抹到了射线所及的无限空间的各个像素上，使得原来像素值为0的点不为0，从而产生星状伪迹，滤波反投影算法用人为设计的一维滤波函数对所得投影数据进行卷积，而后进行反投影和累加时，由于正负抵消，可一定程度上消除星状伪迹。3现在常用的滤波函数有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滤波函数4三、程序实现程序包含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界面。5四、运行结果（1）sheepLogan256*256图像重建（2）自定义灰度图像重建6（3）自定义RGB图像重建（4）不同滤波函数重建效果比较（以sheeplogan图像为例）使用的滤波函数归一化均方距离判据d归一化平均绝对距离判据r最坏情况距离判eParzen0.504870.70730.48133S-L0.636140.991910.49659R-L0.756071.21570.50068可以看出不同滤波函数重建精度：Parzen优于S-L优于R-L。所以后续图像均采用Parzen滤波函数进行滤波处理。五、总结本次程序设计完成了设计目的，投影和反投影过程没有使用自带函数，但是在处理像素比较大的图像时程运行时间比较长，是一个缺点，图像也存在一定误差，效果不是特别完美，不过也算是锻炼了自己的能力，对图像重建中对滤波反投影算法有了深刻的了解，增加了对这门课程的兴趣，期待在以后的学习中能进一步完善程序，缩短运行时间，减小图像误差。7六、附录——程序代码①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:179Protate=imrotate(P,theta);%对原图旋转一个角度，求和（线积分）Pf=sum(Protate,1);[mreal,nreal]=size(Pf);%计算实际尺寸%确定起始点if(nmax-nreal)/2-floor((nmax-nreal)/2)==0From=floor((nmax-nreal)/2+1);%总点数为偶数时elseFrom=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('自定义投影后图像');8z=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%end%end%S-L滤波函数%d=1;%fori=1:2*N-1%h(i)=2/(pi^2*d^2*(1-4*g(i)^2));%end%Parzen滤波函数fori=1:2*N-1h(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));endh(N)=0.0435;%将投影与滤波函数卷积----------------------------------------------------G=zeros(N,180);form=1:180fori=1:Nb=0;9forn=1:Nb=b+h(N+n-i)*R1(n,m);G(i,m)=b;endendend%投影滤波后线性内插进行图像重建----------------------------------------------a=zeros(N);%重建图像初始化，每个像素点像素值为0ns=(N+1)/2;form=1:180%遍历每个投影角度r=pi/180;%将角度换为弧度fori=1:Nforj=1:N%遍历原图的每一个像素点Xrm=(i-ns)*cos(m*r)+(j-ns)*sin(m*r)+ns-1;%坐标转换ifXrm1n=1;%内插运算整数值t=abs(Xrm)-floor(abs(Xrm));%内插运算小数值elsen=floor(Xrm);t=Xrm-floor(Xrm);endifnN-1n=N-1;endc=(1-t)*G(n,m)+t*G(n+1,m);%内插后的滤波投影值a(N+1-j,i)=a(N+1-j,i)+c;%像素值的叠加endendend%将P、a的像素值变换到0-1之间P=double(P);P=P-min(min(P));kk=max(max(P));fori=1:Nforj=1:NP(i,j)=P(i,j)/kk;endenda=double(a);a=a-min(min(a));kkk=max(max(a));fori=1:N10forj=1:Na(i,j)=a(i,j)/kkk;endendsubplot(2,3,5);imshow(a,[]);%重建图形的显示title('滤波反投影重建图像');%重建图像质量评价--------------------------------------------------------%计算归一化均方距离判据d;ave=0;forx=1:Nfory=1:Nave=ave+P(x,y);endendave=ave/(N*N);x1=0;x2=0;forx=1