卷积反投影重建(二维)

卷积反投影图像重建

1 反投影重建基本介绍[1]

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

12(,)(,)A A ωωρφ=。根

据中心切片定理，^

(,)A ρφ可通过),(y x a 在不同视角φ下的投影()r

p x φ的

一维傅氏变换求得。即：

待建图像：

12^

1212()1212

2

^2cos()0

2

cos()

2

cos()

(,)(,)(,)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)

^

121(,)(,)()()(,)r

A A F p x P

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

式(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 θ点的像素值。

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

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

subplot(2,2,2)

imagesc(R1);

title('256*360有噪声平行投影');

colormap(gray)

colorbar

Q=reshape(R1,256,360);

效果图如图-2: 图-2 2.5 滤波器的选择与设计