卷积反投影重建二维

1、卷积反投影图像重建1 1 反投影重建基本介绍A设待重建图像为 a(x,y)a(x,y), ,它的二维傅氏变换为A(SI,CO2)=A(PQ)。根据中心切片定A理，A(P*)可通过 a(x,y)a(x,y)在不到视角由下的投影p/x)的一维傅氏变换求得。即：A(,2)=A(：,)=p(x)=P(：)=P(：,)待建图像：1a(r)=a(x,y)=F2A(1,2)二11A(1,.2)ei(X.2y)d.，1d.24 二2-：_：、一1.A二A(P*)ei2*rc0s(e)|P|dPd1,0-00=PJP(p*)ei2JTPrcos(8_*)|PdPd 巾，0-OO=/d广|pp(P,*)ei2jT

2、pcos(0)dP03H因为xr=rcos(8-*),所以有:xf2y=2”.(xcos二，ysin)=2二Pxr=2：rcos(-)同时：do1d2=JdPd4先来看该式的第二个积分：口归匕性=口耳-3dp=h(*)P(xr,)|xr=rcosf一)=g(x,)cosh-)-grcos(-),1式中:g(x,户h(x)p(x,(1.4)(1.4)式(3.10)(3.10)的物理意义是投影p(xr仲)经过传递函数为|H=|H=F F1h(x)的滤波器后得到的修正后的投影g(x2)在满足xr=rcos(8-旬时的值。将(3.11)(3.11)代入(3.8)(3.8),得到：Aa(r,B)=(gr

3、cos(日一4),4d(1.5)(1.5)称为滤波反投影方程，其物理意义是经过给定点(r,的所有滤波后的投影在=0元范围内的累加一反投影重建，得出(r,B)点的像素值。(1.1)2二cos2二sin-2nPsinO2nPcos二4二2P(1.2)(1.3)可见，滤波(卷积)反投影算法的具体包含三大步：把在固定视角下测得的投影p(x,鹤经过滤波，得到滤波后的投影g(xr,*)；(2)(2)对每一个把g(xr*)反投影于满足xr=rcos(8-4)的射线上的所有各点(方);(3)(3)将步骤(2)(2)中的反投影值对所有 006gn6gn 进行累加(积分)，得到重建后的图像。2 重建流程2.1首先

4、我们利用 phantomphantom()()函数产生一个头部幻影图像，用以检测二维重建算法，代码如下：I=phantom(256);subplot(2,2,1)imshow(I,);Utle(256*256 原始图像);效果图如图 1 1 所示,为一个大椭圆和几个小椭圆。2.2初始参数设置重建采用的是平移加旋转的扫描方式，射线源在某一角度下水平移动，将物体全部照射后旋转一角度，如此重复，在这个过程中探测器相应地运动以接收 X X 射线。根据此原理，将重建程序的初始参数设置如下：N,N=size(I);z=2*ceil(norm(size(I)-floor(size(I)-1)/2)-1)+3

5、;%radon 变换默认平移点数/角度Nt=360;%角度采样点数Nd=N;%平移数x=pi/180;%角度增量d=N/Nd;%平移步长theta=1:Nt;a=zeros(N);2.3产生无噪声投影数据256256原始图图1R,xp=radon(I,theta);e=floor(z-Nd)+2;R=R(e:(Nd+e-1),:);R1=reshape(R,256,360);radon(I,theta)产生 I I 投影，默认 z z 点/角度，即使指定 N N 点也是 z z 点.所以为避免重建图像放大或缩小，下面计算取投影时需补偿，补偿量 e e 如对 256256 的图像，补偿为 55,

6、55,即 pmllpmll 勺第 5555 个点作为计算用的第一个投影. .2.4添加噪声并将有噪声平行投影进行显小mm,nn=size(R1);di=lognrnd(0,0.15,mm,nn);R1=10*(R1-min(R1(:)/(max(R1(:)-min(R1(:);I0=1.5e5;%incidentphotons;decreasethisforsimulatinglowdosescansrand(state,0),randn(state,0);yi=poissrnd(I0*di.*exp(-R1)+3*randn(size(R1);ifany(yi(:)=0)warn(%dof%

7、dvaluesare0insinogram!,.sum(yi(:)=0),length(yi(:);endR1=log(I0./max(yi,0.01);%noisysinogramR1=max(R1,0);%显示ff=2;uu=22000;v=ff*exp(R1/uu);subplot(2,2,2)imagesc(RI);title(256*360 有噪声平行投影,);colormap(gray)colorbarQ=reshape(R1,256,360);效果图如图-2:2.5滤波器的选择与设计Bf2h(s)=J|5exp(/2;rSs)dS=最基本的从投影重建图像的滤波器:示)。19711

8、971 年提出的 R-LR-L 重建滤波器(下图中实线表空域表达式为: :图-3R-L-3R-L 滤波器示意图(2.1(2.1) )其中 B 为截至频率。若2,d 为空间采样间隔，可解256#360 有噪声平行投影50100150200250图-2isnr100200300151050出离散的滤波器空间脉冲响应:fm=04d2hd)=fg(.v,w离散化就S得r*2.A/Afg式泄)二pm：d,心)=Zg/QH用-k)-?nLMJr=其中从滤波因尸方采样得到Aft(HI)=考虑g(JW)线性内插为欧必)gG心)=治(川)+(；-叫心(川+1)-京(加)a其中mdsin+l)rf反投影枳分NT/

9、(.Y.1)=Zg(MSQ5十sinn-0在角坐标离散位NT-Ag|_(/Arcos?A+/AvsinnA),nAPI=O/,j0,h2,1t-代码设计：form=1:Ntpm=Q(:,m);u=conv(hl,pm);pm=u(k);Cm=(N-1)/2)*(1-cos(m-1)*x)-sin(m-1)*x);fori=1:Nforj=1:NXrm=Cm+(j-1)*cos(m-1)*x)+(i-1)*sin(m-1)*x);ifXrm(Nd-1)n=Nd-1;endp=(1-t)*pm(n)+t*pm(n+1);a(N+1-i,j)=a(N+1-i,j)+p;endendend重建后的图像

10、如图-3-3 所示：程序中还包含对重建结果的评价/归一化均方距离判据/归一化平均绝对距离判据以及程序的运行时间等.不再进行详细介绍。程序的最终运行效果如下：重建图像图-31.510.50050100第230ff的像素图-4 程序最终运行效果图3 分析反投影重建方法包括卷积反投影重建的缺点是会产生星状伪迹，原因分析如下：断层平面中某一点的密度值可以看作是这一平面内所有经过该点的射线的投影值之和（的均值）。整幅重建图像可以看作是所有方向下的投影累加而成。射线标号示于图 5 5 中，像素值（代表密度）分别 x x1,x,x2,X,X3,X,X4,赋值如下：x1=5,x2=0,x3=2,x4=18根据

11、投影的定义（某条射线投影值为该条射线穿过的所有的像素值之和） ，每条射线的投影p（i=1,2- -）为：Pl=xi*x2=5, ,p2=x3*x4=20, ,p3=xi+x3=7p4=x2x4=18,p5=x3=2, ,p6=x1x4=23P7=x2=0256*256原始图像256*360有噪声平行投影inr上6iflV:w|重建图像图-5 断层像素值和射线根据反投影重建算法的物理意义， 重建图像中各像素， 得到:xi=pi+p3+p6=35,x2=pi+p4+p7=23,x3=p2p3p5=29x4=p2p4P6=61,共 9 页第 8 页重建后的图像如图 6(b)6(b)所示，可以看出原图

12、像中像素值不为零的点反投影重建后仍较突出，但原图中像素值为零的点，经反投影重建后不再为零，即有伪迹。有时为了使重建后图像的像素值更接近于原图的像素值，在求反投影时，把数据除以投影的数目(即射线数)，如图 6(c)6(c)所示。因此有：1XkPk,inpi=1该式可作为反投影重建算法的计算式。其中 X Xk表示像素 k k 的值，Pk,i表示经过像素k的第i条射线投影，np表示图像内的射线条数。图 7(a)7(a)表示空间中一个孤立点源 A,A,密度为 1 1。经过 A A 点的三条射线也示于图中。射线束理论上可以很多， 取三条示意。 不经过 A A 点的射线投影为零， 经过 A A 点的射线投影值均为 1,1,P1=P2=P3=1。5021835262961(3.1)(3.1)(a)原图像像素值(b)反投影重建后图像图 6 反投影示例(c)求平均后图像图 7 孤立点源的反投影重建及星状伪迹经反投影重建后，得到 A A 点的像素值为fA=(R+P2+P3)/3=1。