卷积反投影重建(二维).doc

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

卷积反投影图像重建

1 反投影重建基本介绍[1]
设待重建图像为，它的二维傅氏变换为。根据中心切片定理，可通过在不同视角下的投影的一维傅氏变换求得。即：

待建图像：
(1.1)

因为，所以有：

同时：

(1.2)
先来看该式的第二个积分：
(1.3)
式中：
(1.4)
式(3.10)的物理意义是投影经过传递函数为的滤波器后得到的修正后的投影在满足时的值。将(3.11)代入(3.8)，得到：
(1.5)
称为滤波反投影方程，其物理意义是经过给定点的所有滤波后的投影在~范围内的累加—反投影重建，得出点的像素值。

可见，滤波（卷积）反投影算法的具体包含三大步：

(1) 把在固定视角下测得的投影经过滤波，得到滤波后的投影；
(2) 对每一个，把反投影于满足的射线上的所有各点；
(3) 将步骤(2)中的反投影值对所有<进行累加（积分），得到重建后的图像。
2 重建流程
2.1 首先我们利用phantom（）函数产生一个头部幻影图像，用以检测二维重建算法，代码如下：
I=phantom(256);
subplot(2,2,1)
imshow(I,[]);
title('256*256原始图像');
效果图如图1所示， 图1
为一个大椭圆和几个小椭圆。
2.2 初始参数设置
重建采用的是平移加旋转的扫描方式，射线源在某一角度下水平移动，将物体全部照射后旋转一角度，如此重复，在这个过程中探测器相应地运动以接收X射线。根据此原理,将重建程序的初始参数设置如下:
[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(:)));
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 滤波器的选择与设计
最基本的从投影重建图像的滤波器：1971年提出的R-L重建滤波器(下图中实线表示)。
图-3 R-L 滤波器示意图
空域表达式为:
(2.1)

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

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 重建图像

代码设计:
for m=1:Nt
pm=Q(:,m);
u=conv(hl,pm);
pm=u(k);
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
程序中还包含对重建结果的评价/归一化均方距离判据/归一化平均绝对距离判据以及程序的运行时间等.不再进行详细介绍。程序的最终运行效果如下：

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

整幅重建图像可以看作是所有方向下的投影累加而成。平移
射线标号示于图5中，像素值（代表密度）分别,,,，
赋值如下:,,,

图5 图像重建像素值

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

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

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

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

该式可作为反投影重建算法的计算式。其中表示像素的值，表示经过像素的第条射线投影，表示图像内的射线条数。
图7(a)表示空间中一个孤立点源A，密度为1。经过A点的三条射线也示于图中。射线束理论上可以很多，取三条示意。不经过A点的射线投影为零，经过A点的射线投影值均为1，
。

(a) 孤立点源A及三条射线 (b)相应的反投影重建图，有星状伪迹
图 7 孤立点源的反投影重建及星状伪迹

经反投影重建后，得到A点的像素值为。A点以外的像素值原来为零，经反投影重建后不等于零，而是等于1/3。所以，经反投影重建后的图像除保留A点的像外，还有像素值为1/n的灰雾背景，后者称为星状伪迹。产生星状伪迹的原因在于：反投影的本质是把取自有限物体空间的投影均匀地回抹（反投影）到射线所及的无限空间的各点上，包括原先像素值为零的点。
参考文献
[1] 滤波反投影,Zhengdong Zhou.2012
共9页第11页