卷积反投影重建(二维)

计算机图像处理技术结课（报告） 共9页第1页 卷积反投影图像重建

1 反投影重建基本介绍 设待重建图像为),(yxa，它的二维傅氏变换为^12(,)(,)AA。根据中心切片定理，^(,)A可通过),(yxa在不同视角下的投影()rpx的一维傅氏变换求得。即：

待建图像：

12^1212()12122^2cos()02cos()02cos()0(,)(,)(,)1(,)4(,)(,)(,)ixyiririraraxyFAAeddAeddPedddPed



















[] (1.1) 因为cos()rxr，所以有： 122(cossin)22cos()rxyxyxr

同时：

12ddJdd

112

22

//2cos2sin4//2sin2cosJ

 (1.2)

先来看该式的第二个积分：

22cos()cos()cos()cos()(,)(,)|()(,)|(,)|cos(),rrrrixirxrrrxrrxrPedPedhxpxgxgr









(1.3) 式中：(,)()(,)rrrgxhxpx (1.4)

^121(,)(,)()()(,)rAAFpxPP计算机图像处理技术结课（报告）

共9页第2页 式(3.10)的物理意义是投影(,)rpx经过传递函数为1[()]rFhx的滤波器后得到的修正后的投影(,)rgx在满足cos()rxr时的值。将(3.11)代入(3.8)，得到：

^

0(,)[cos(),]argrd (1.5)

称为滤波反投影方程，其物理意义是经过给定点(,)r的所有滤波后的投影在0~范围内的累加—反投影重建，得出(,)r点的像素值。

可见，滤波（卷积）反投影算法的具体包含三大步： (1) 把在固定视角i下测得的投影(,)rpx经过滤波，得到滤波后的投影(,)rgx； (2) 对每一个i，把(,)rigx反投影于满足cos()rixr的射线上的所有各点(,)r；

(3) 将步骤(2)中的反投影值对所有0<进行累加（积分），

得到重建后的图像。 2 重建流程 2.1 首先我们利用phantom（）函数产生一个头部幻影图像，用以检测二维重建算法，代码如下： I=phantom(256); subplot(2,2,1) imshow(I,[]); title('256*256原始图像'); 效果图如图1所示， 图1 为一个大椭圆和几个小椭圆。 2.2 初始参数设置 重建采用的是平移加旋转的扫描方式，射线源在某一角度下水平移动，将物体全部照射后旋转一角度，如此重复，在这个过程中探测器相应地运动以接收X射线。根据此原理,将重建程序的初始参数设置如下: 计算机图像处理技术结课（报告） 共9页第3页 [N,N]=size(I); 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 产生无噪声投影数据 [R,xp] = radon(I,theta); e=floor((z-Nd)/2)+2; R=R(e:(Nd+e-1),:); R1=reshape(R,256,360); radon(I,theta)产生I投影,默认z点/角度,即使指定N点也是z点.所以为避免重建图像放大或缩小,下面计算取投影时需补偿,补偿量e,如对256的图像,补偿为55,即pm的第55个点作为计算用的第一个投影。 2.4 添加噪声并将所有噪声平行投影进行显示 [mm,nn]=size(R1); di=lognrnd(0,0.15,mm,nn); R1= 10*(R1-min(R1(:)))/( max(R1(:))-min(R1(:))); 计算机图像处理技术结课（报告） 共9页第4页 I0 = 1.5e5; % incident photons; decrease this for simulating "low dose" scans rand('state', 0), randn('state', 0); yi= poissrnd(I0 * di.*exp(-R1))+3*randn(size(R1)); if any(yi(:) == 0) warn('%d of %d values are 0 in sinogram!', ... sum(yi(:)==0), length(yi(:))); end R1 = log(I0 ./ max(yi,0.01)); % noisy sinogram R1=max(R1,0); % 显示 ff=2; uu=22000; v=ff*exp(R1/uu); subplot(2,2,2) imagesc(R1); title('256*360有噪声平行投影'); colormap(gray) colorbar Q=reshape(R1,256,360); 效果图如图-2: 图-2 2.5 滤波器的选择与设计 计算机图像处理技术结课（报告） 共9页第5页 最基本的从投影重建图像的滤波器：1971年提出的R-L重建滤波器(下图中实线表示)。

图-3 R-L 滤波器示意图 空域表达式为:

(2.1) 其中B为截至频率。若，d为空间采样间隔，可解出离散的滤波器空间脉冲响应：

(2.2) 计算机图像处理技术结课（报告）

共9页第6页 MATLAB 代码如下: g=-(Nd/2-1):(Nd/2); for i=1:256 if g(i)==0 hl(i)=1/(4*d^2); else if mod(g(i),2)==0 hl(i)=0; else hl(i)=(-1)/(pi^2*d^2*(g(i)^2)); end end end k=Nd/2:(3*Nd/2-1);

2.6 重建图像 计算机图像处理技术结课（报告）

共9页第7页

代码设计: for m=1:Nt pm=Q(:,m); u=conv(hl,pm); pm=u(k); 计算机图像处理技术结课（报告） 共9页第8页 Cm=((N-1)/2)*(1-cos((m-1)*x)-sin((m-1)*x)); for i=1:N for j=1:N Xrm=Cm+(j-1)*cos((m-1)*x)+(i-1)*sin((m-1)*x); if Xrm<1 n=1; t=abs(Xrm)-floor(abs(Xrm)); else n=floor(Xrm); t=Xrm-floor(Xrm); end if n>(Nd-1) n=Nd-1; end p=(1-t)*pm(n)+t*pm(n+1); a(N+1-i,j)=a(N+1-i,j)+p; end end end 重建后的图像如图-3 所示: 图-3 程序中还包含对重建结果的评价/归一化均方距离判据/归一化平均绝对距离判据以及程序的运行时间等.不再进行详细介绍。程计算机图像处理技术结课（报告） 共9页第9页 序的最终运行效果如下：

图-4 程序最终运行效果图 3、分析 反投影重建方法包括卷积反投影重建的缺点是会产生星状伪迹，原因分析如下： 断层平面中某一点的密度值可以看作是这一平面内所有经过该点的射线的投影值之和（的均值）。

整幅重建图像可以看作是所有方向下的投影累加而成。射线标号示于图5中，像素值（代表密度）分别1x,2x,3x,4x， 赋值如下:15x,20x,32x,418x 计算机图像处理技术结课（报告） 共9页第10页

图5 图像重建像素值

根据投影的定义(某条射线投影值为该条射线穿过的所有的像素值之和)，每条射线的投影ip(1,2i)为：

1215pxx, 23420pxx,3137pxx 42418pxx, 532px,

61423pxx

720px

根据反投影重建算法的物理意义，重建图像中各像素，得到： 113635xppp，214723xppp，

323529xppp424661xppp，

5

218

0

35296126 54.18.73.3 (a) 原图像像素值 (b)反投影重建后图像 (c)求平均后图像 图 6 反投影示例

重建后的图像如图6(b)所示，可以看出原图像中像素值不为零的点反投影重建后仍较突出，但原图中像素值为零的点，经反投影重建后不再为零，即有伪迹。有时为了使重建后图像的像素值更接近于原图的像素值，在求反投影时，把数据除以投影的数目（即射线数），如图6(c)所示。

因此有：

,11pnkkiipxpn



(3.1)

50 218

(1)(2)

(3)(4)

(5)(6)(7)