卷积反投影重建(二维)

卷积反投影图像重建

1 反投影重建基本介绍

设待重建图像为),(y x a ，它的二维傅氏变换为^

12(,)(,)A A ωωρφ=。根据中心切片定理，^

(,)A ρφ可通过),(y x a 在不同视角φ下的投影()r p x φ的一维傅氏变换求得。即：

待建图像：

12^

1212()1212

2

^2cos()0

2cos()

02cos()0

(,)(,)(,)1(,)4(,)(,)(,)i x y i r i r i r a r a x y F A A e d d A e d d P e

d d d P

e d ωωππρθφπ

πρθφπ

πρθφθωωωωωωπρφρρφ

ρφρρφ

φρρφρ

-∞∞

+-∞-∞

∞--∞

∞

--∞

∞--∞

======⎰⎰

⎰⎰⎰

⎰

⎰⎰

[]

(1.1)

因为cos()r x r θφ=-，所以有：

122(cos sin )22cos()r x y x y x r ωωπρφφπρπρθφ+=+==-

同时：

12d d J d d ωωρφ= 11222//2cos 2sin 4//2sin 2cos J ωρωφπφ

πρφ

πρωρ

ωφπφ

πρφ

∂∂∂∂-=

=

=∂∂∂∂ (1.2)

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

[]

22cos()

cos()cos()cos()(,)(,)|()(,)|(,)|cos(),r r r r i x i r x r r r x r r x r P e

d P

e d h x p x g x g r πρπρθφθφθφθφρρφρρρφρ

φφθφφ∞

∞

-=--∞

-∞

=-=-==*==-⎰

⎰ (1.3)

式中：

(,)()(,r r r

g x h x p x φφ=* (1.4) 式(3.10)的物理意义是投影(,)r p x φ经过传递函数为1[()]r F h x ρ=的滤波器后得到的修正后的投影(,)r g x φ在满足cos()r x r θφ=-时的值。将(3.11)代入(3.8)，得到：

^

(,)[cos(),]a r g r d π

θθφφφ=-⎰ (1.5)

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

^

121(,)(,)()()(,)r A A F p x P

P φφωωρφρρφ⎡⎤====⎣⎦

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

(1) 把在固定视角i φ下测得的投影(,)r p x φ经过滤波，得到滤波后的投影(,)r g x φ； (2) 对每一个i φ，把(,)r i g x φ反投影于满足cos()r i x r θφ=-的射线上的所有各点(,)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 射线。根据此原理,将重建程序的初始参数设置如下:

[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);