卷积反投影重建(二维)

1、计算机图像处理技术结课(报告)卷积反投影图像重建1反投影重建基本介绍A设待重建图像为a(x, y)，它的二维傅氏变换为 A( .1, .2)=A(几)。根据中心切片定A理，A( ：?,)可通过a(x,y)在不同视角 '下的投影p (Xr)的一维傅氏变换求得。即：A( 1，2)= A(，)二 F1 | P (xr )二 P ()= P()待建图像：_1a(rj)二 a(x,y)二 F?一A(匕)厶 i i A(2)ei( 1X '2y)d24 二=j" A(P)ei2用rcos(日P|d Pd©(1.1)0=r P(pf)ei2“PrcosZ)| p|d Pd

2、©0 -=fd 叮 °°|p P(Pf)ei2jlFrcosdd pJ0因为xrcosh -)，所以有：x r：'2y =2 J(xcos：；亠 ysin ) =2二x, =2：rcosf )同时：d创dco2 = J d Pd©2 cos2二 sin '-2 兀 Ps in©2 兀 Pcos©(1.2)先来看该式的第二个积分：门护(匕怩2叫叫戸=彷|时力严lecosg dP二皿片)以片，)lx严 cos-J二g(x,)" cosd-)=grcos( -),式中:g(x , ) h (x ) p (x ,(1

3、.4)式(3.10)的物理意义是投影p(xr,)经过传递函数为Fdh(Xr)的滤波器后得到的修正后的投影g(x,)在满足xr二rcos()时的值。将(3.11)代入(3.8)，得到:Aa(r,日)=J0 gr cos(日一©),©d©(1.5)称为滤波反投影方程，其物理意义是经过给定点(r,R的所有滤波后的投影在 -0-：范围内的累加一反投影重建，得出(rj)点的像素值。计算机图像处理技术结课(报告)可见，滤波(卷积)反投影算法的具体包含三大步：(1)把在固定视角】下测得的投影P(Xr,)经过滤波，得到滤波后的投影g(Xr,)； 对每一个1，把g(xr, i)反投

4、影于满足xr rcos(J的射线上的所有各点(r,v)；(3)将步骤 中的反投影值对所有Ov乜二进行累加(积分)，得到重建后的图像。2重建流程256-2 图像图12.1首先我们利用phantom ()函数产生一个头部幻影图像，用以检测二维重建算 法，代码如下：I=pha ntom(256);subplot(2,2,1)imshow(l,);title('256*256 原始图像');效果图如图1所示，为一个大椭圆和几个小椭圆。2.2初始参数设置重建采用的是平移加旋转的扫描方式，射线源在某一角度下水平移动，将物体全部 照射后旋转一角度，如此重复，在这个过程中探测器相应地运动以接收

5、X射线。根据此原理,将重建程序的初始参数设置如下：N,N=size(l);z=2*ceil(norm(size(I)-floor(size(I)-1)/2)-1)+3;% radon 变换默认平移点数/角度Nt=360;%角度采样点数Nd=N; %平移数x二pi/180; %角度增量 d=N/Nd; %平移步长theta = 1:Nt;a=zeros(N);2.3产生无噪声投影数据共9页第4页计算机图像处理技术结课（报告）R,xp = rad on (l,theta);e二floor(z-Nd)/2)+2;R=R(e:(Nd+e-1),:);R仁 reshape(R,256,360);rado

6、n(I,theta)产生I投影，默认z点/角度，即使指定N点也是z点.所以为避免重建 图像放大或缩小,下面计算取投影时需补偿,补偿量e如对256的图像,补偿为55,即pm勺 第55个点作为计算用的第一个投影.2.4添加噪声并将有噪声平行投影进行显示mm, nn =size(R1);di=log nrn d(0,0.15,mm, nn);R1= 10*(R1-mi n(R1(:)/( max(R1(:)-mi n(R1(:);I0 = 1.5e5; % in cide nt phot ons; decrease this for simulati ng "low dose"

7、sca ns ran d('state', 0), randn ('state', 0);yi= poissrnd(I0 * di.*exp(-R1)+3*ra ndn (size(R1);if any (yi(:) = 0)warn ('%d of %d values are 0 in sinogram!', .sum(yi(:)=0), le ngth(yi(:);endR1 = log(I0 ./ max(yi,0.01); % n oisy si nogramR1=max(R1,0);%显示ff=2;uu=22000;v=ff*exp(R

8、1/uu);共9页第5页计算机图像处理技术结课（报告）15subplot(2,2,2)imagesc(RI);title('256*360有噪声平行投影');colormap(gray)colorbarQ=reshape(R1,256,360);效果图如图-2:256*360有噪声平行投影50100150200250100 200 300图-2-102.5滤波器的选择与设计最基本的从投影重建图像的滤波器：1971年提出的R-L重建滤波器（下图中实线表 示）。图-3 R-L滤波器示意图 空域表达式为：Bf2h(s) = J |ijexp(y2S5)</5 =*严、sm()a

9、*订s、 i snr()jrs=JIS-()2* J_ d1(2.1)其中B为截至频率。若, d为空间采样间隔，可解共9页第7页计算机图像处理技术结课（报告）出离散的滤波器空间脉冲响应:7m = 0=«0iv is eivn(2.2)丁丿=m is oddMATLAB弋码如下：g=-(Nd/2-1):(Nd/2);for i=1:256if g(i)=0hl(i)=1/(4*dA2);else if mod(g(i),2)=0hl(i)=0;elsehl(i)=(-1)/(piA2*dA2*(g(i)A2);endendendk=Nd/2:(3*Nd/2-1);2.6重建图像卷积 g

10、(斗 0) = gGE)"G)反投影=0离散化3)得-11<tn< gff(»?) = g(d,?a)= Y gK)h(m-k)買屮从滤波禺MG）采样得到h(in) = dh(itnl)考虑g (加)线性内插为£(/")g(兵 /?A) =乳(m) + (J - 加)衣(川 +1) - 或(山) a其中 wd <s< (m + l)rf反投影积分JV-1/ (.Y. 1 ) = A(ACQSf7 十 win M )，/? n-0il角坐标r(X. V)离散值i9J 0»±L ±2< -代码设计：

11、for m=1:Ntpm=Q(:,m);u=c on v(hl,pm);pm=u(k);Cm=(N-1)/2)*(1-cos(m-1)*x)-si n(m-1)*x);for i=1:Nfor j=1:NXrm二Cm+(j-1)*cos(m-1)*x)+(i-1)*si n(m-1)*x);if Xrm<1n=1;t=abs(Xrm)-floor(abs(Xrm);elsen=floor(Xrm); t=Xrm-floor(Xrm); end if n>(Nd-1)n=Nd-1;endp=(1-t)*pm( n)+t*pm( n+1); a(N+1-i,j)=a(N+1-i,j)+

12、p;end重建图像图-3endend重建后的图像如图-3所示:程序中还包含对重建结果的评价/归一化均方距离判据/归一化平均绝对距离判据 以及程序的运行时间等不再进行详细介绍。程序的最终运行效果如下：共9页第10页计算机图像处理技术结课（报告）2於360有噪声平行投影256*256原始图像50100150200250100 200 300151050重建图像1.50.550100第2301亍的像素图-5断层像素值和射线图-4程序最终运行效果图3分析反投影重建方法包括卷积反投影重建的缺点是会产生星状伪迹，原因分析如下： 断层平面中某一点的密度值可以看作是这一平面内所有经过该点的射线的投影值之和 （

13、的均值）。整幅重建图像可以看作是所有方向下的投影累加而成。射线标号示于图5中，像素值（代表密度）分别 , x2, x3, x4，赋值如下：x1 =5, x2 =0 , X3 =2 , x4 =18根据投影的定义（某条射线投影值为该条射线穿过的所 有的像素值之和），每条射线的投影Pi （ i 1,2）为：Pl = Xi x2 = 5 ,P = x3 x4 = 20 , p3 = xi x3 = 7p = x2 x4 = 18,P5 = X3 = 2 ,p6 =为 X4 = 23P7 = X2 二 0根据反投影重建算法的物理意义，重建图像中各像素，得到:X P1 P3 P6 =35，X2 =P1

14、P4 P7 =23，X3 =P2 P3 P5 =29 X4 =P2 P4 P6 =61，共9页第8页计算机图像处理技术结课(报告)50218(a)原图像像素值35262961(b)反投影重建后图像图6反投影示例8.7(c)求平均后图像共9页第13页重建后的图像如图6(b)所示，可以看出原图像中像素值不为零的点反投影重建后 仍较突出，但原图中像素值为零的点，经反投影重建后不再为零，即有伪迹。有时为了 使重建后图像的像素值更接近于原图的像素值，在求反投影时，把数据除以投影的数目 (即射线数)，如图6(c)所示。因此有：nPXknpi 二1Pk,i(3.1)共9页第#页共9页第#页该式可作为反投影重建算法的计算式。其中Xk表示像素k的值，Pk,i表示经过像素 k的第i条射线投影，np表示图像内的射线条数。图7(a)表示空间中一个孤立点源 A，密度为1。经过A点的三条射线也示于图中。 射线束理论上可以很多，取三条示意。不经过 A点的射线投影为零，经过 A点的射线 投影值均为1， p1 p2 p3 =1。(a)孤立点源A及三条射线(b