基于投影的图像重建

基于投影的图像重建

1 雷登变换和中心投影定理

1.1 雷登变换

以电脑断层扫描技术(CT )为代表的各种层析成像方法对被测对象产生的并不是其直接三维图像，而是从一组投影中重建而成的该物体截面的二维图像。对于CT 而言，影像的性质就是线性衰减系数(,)x y μ。相应的X 射线沿某条路径的变换可由线积分来给出1

0ln (/)I I =path

?

()x μds (1)

其中ds 是路径长度微元，()x μ(,)μ=x y 是描述被成像物体的线性衰减系数二维分布的函数。数学上可将一个投影用雷登变换来描述。雷登变换是一个将一二维函数向极坐标空间',()x θ上的投影变换(见图1)，这一变换表达为2

2'

()x d P θ+∞

-∞=?x 2g ()x (δx n '

)-x (2)

其中n 即(cos ,sin )θθ，是一个与投影射束正交的单位向量；x 是一个在投影射束方向上

的向量。'

()x P θ称为一个投影；同时为了一般性考虑，影像的性质用二维函数()g x 来表

示。狄拉克δ函数(()x δ)“挑选出”与投影平行的直线路径。显然对于平行投影的数据，我们得到

'()x P θπ

+-= '()x P θ (3) 记住基本的投影几何中变换'=x '()x x 由以下两个式子给出： 'c o s s i n x x y θθ=+

'sin cos y x y θθ=-+

1

或者

0l n (/)I I =p a t h

?

()μx d x

2

2x dxdy d ≡

现在因为旋转是一个保面积变换3，所以我们在'',()x y 坐标系的投影为

''

2()x d P θ

+∞

-∞

=?'x g '()x '''()x x δ-

''''(,)dy g x y +∞

-∞

=?

(4)

图 1：图示为平行投影射束图

后一个等式可以等效写成为 ''

'()()x dy g x P θ+∞

-∞

=? (5)

1.2 傅里叶切片定理

投影的傅里叶变换()S θω为

3

或者该换的函数行列式可记作 ''(,)(,)

x y x y ??=1

'''

exp ()[]()dx x i x S P θθωω∞

-∞

-=?

=2

d ∞

-∞?x g '()x 'exp []i x ω-

=2

d ∞

-∞?x g ()x exp [i ω-x n ] (6)

其中我们可以利用等式(5)得到最后两个表达式。 现在有g ()x 的二维傅里叶变换，即,()u v G G =()u 记为 G ()u =

2

d

∞

-∞?x g ()x exp [i -u x ] (7)

变换到傅氏域中的极坐标系(,)ωθ

cos u ω=θ

sin v ω=θ

我们有

,()G ωθ=

2

d

∞

-∞?x g ()x e x p [i ω-x n ] (8)

比较等式(6)和等式(8)我们有以下结果

()S θω=,()ωθG (9) 这样，投影的傅里叶变换和原物体在与投影射束正交的射束方向上的频谱是完全相同的。如图2所示，这个重要的结果被称为傅里叶切片定理或者中心切片定理。换句话说，我们已经在与投影成正交的傅氏空间中确定了一组傅里叶系数。

图 2：图示为傅里叶切片定理图 2. 平行射束数据连续滤波反投影的等式

等式(8)的傅里叶逆变换为

g ()x 220

014d d π

πθω∞=??,()G ωωθexp [i ωx n ]

220014d d ππ

θω∞=??,()G ωωθe x p c o s s i n [{()()}]i x y ωθθ+ (10) 该积分式可被分解为两个积分的和

g ()x 2

001

4d d ππ

θω∞

=??,()G ωθexp cos sin [{()()}]i x y ωθθ+

+2

14π

θωπ

∞

??d d (,)ωωθπ+G exp cos sin [{()()}]i x y ωθπθπ+++

现有cos cos ()()θπθ+=-，sin sin ()()θπθ+=-所以可得到,,()()G G ωθπωθ+=-(见等式(8)),所以有

g ()x 2

01

4d d ππθω

∞

-∞=

??,()G ωωθexp [i ω-x n ] (11)

从傅里叶切片定理(等式(9))可以得到

g ()x 20

014d d π

πθω∞=

??()S θωωe x p [i ω-x n ]

20

14d Q π

θπ

θ=?()x n (反投影) (12) 其中

θQ ()x n 2

14()θωωωπ∞

-∞=

?d S exp [i ω-x n ]

1

12π-=

F

[()]S θωω

1

1

12{[πω--=*F F F

'

()]}x P θ

1

12[]ωπ

-=*F '

()x P θ

'

()()h x P θ=* (滤波) (13)

1

-F

表示傅里叶变换的逆变换，*表示卷积而h 1

([]/(2))ωπ-F

则是滤波脉冲响

应。这样，等式(13)就代表对斜坡滤波器ω产生的投影起到的滤波作用。这个经过滤波的投影θQ 将会对图像(,)x y 中所有沿投影的点起到相同的作用，进而可以被看作为是反投影(见等式(12))。这样(,)g x y 的值将由在180°范围内所有这样的投影的和构成(见图3)。因此可将这样一个重建的过程命名为滤波反投影。

3 有限频带的作用

3.1 对一个点的重建 考虑以下的图像

g ()x δ=()x (14) 显然有

?g ()x 2d x δ=?

()x 2

d x (15)

图 3：图示为平行射束反投影的图标(见等式(12))

从等式(5)可以我们得到投影

()P θω'dy δ∞

-∞=?

'

()x

'()x δ= (对于所有的θ而言) (16)

这样对投影的傅里叶变换为

()θωS '''e x p []()x x d x

i ωδ∞

-∞=-? 1= (对于所有的θ而言) (17)

利用等式(10)我们可以得到

?g

()x 2

20

1

4d d π

πθω∞

=??exp [i ωω]x n (18)

其中?g

()x 是对g ()x 的重建。考虑到重建的有限带宽的可积情况，我们得到 ?g

()x 201

4s

d ωπ

ω=?20

e x p [d i π

ωθω?]x n (19)

其中ωs 是空域中的波数，否则系统是不敏感的。所有在这个平面上的任意角度的像点都可以用极坐标来表示为

c o s ()x r θ=

s i n

()y r θ= 这样重建图像可由下式给出

?g

20

1,4()s

r d ωπφω=?20exp cos cos sin sin [{()()()()}]r d i πωθωφθφθ+?

2014s

d ωπ

ω=?20exp cos [()]d i r πωθωθφ-? 令ψθφ=-，可以得到

?g

201,4()s

r d ωπ

φω=?2exp cos [()]d i r πφ

φ

ωψωψ--?

(20)

等式(20)中的第二重积分是在一个完整的圆上，这样就有

?g

()x 20

1

4s

d ωπω=?20

e x p [d i π

ωψω?cos ()]ψx (21)

更进一步来说，因为cos()ψ是偶函数，则上式可以写成

?g

()x 2

21s

d ωπω=?0

e x p [d i π

ωψω?c o s ()]ψx (22)

幸运的是，对ψ积分的求解有一个标准的表解法(见Abromovitz and Stegun 中9.1.21，第360页)，由下式给出

001

e x p c o s ()[(

)]z i z d J π

θθπ

=

? (23)

其中J ν是零阶第一类贝塞尔函数(例见Kreyszig 中第208页)，这样得到

?g ()x 012s

d ωωπ

=?0(J ωω)x (24) 利用替换式t =ωx 上式可以由下式给出

?g

()x 2

2πω=s

x

s x

ω?d t 0()t t J (25)

幸运的是，这个积分的求解同样有一个标准的表解法(A&S 中11.3.20，第484页)，即

10

()z

t dt z J z J νννν-=? (26)

其中J 是ν阶第一类贝塞尔函数。这样，对于1ν=，s z ω=x 条件而言，我们可以得到

?g

()x 2πω=s

x

1J (s ω)x (27) 图4显示了对于不同的s ω值的结果。可以注意到一个很有趣的结果，这与由圆孔夫

琅禾费衍射得到的结果是一个形式。这样，有限带宽的作用即是对图像像素的模糊作用。

可以利用对等式(27)在2R 上的积分来验证这个结论是正确的。于是就有

20

d π

θ∞

??

x ?g

()x d x 10

()s s d J x x ωω∞

=? 1=

其中我们已经用过了标准结果，即

01()t dt J

ν

∞

=?

图 4：有限带宽对平面中心的图像象点重建的作用(见等式(27))。(a)2ω=s ,(b)4ω=s , (c)

6ω=s (如

果距离单位是mm ，则ωs 的单

位是rad mm -

)

3.2 一般情况

我们考取更加一般的情况。令g ()x 为图形性质的实际分布。例如在X 射线CT 的

g ()x =()x μ其中()x μ是线性吸收系数的二维分布函数。任何图像都可以描述成如下式

g ()x g ∞

-∞=?'()x '()δ-x x 2d 'x (28)

利用这个重建的图像，?g

()x (对于连续的情况)将和实际的图像g ()x 通过下式相联系起来

?g

()x [h g =()]x =g ∞

-∞?'()x '()δ-x x 2d 'x

g ∞

-∞=?'()x [h '()δ-x x ]2d 'x

从前面部分我们已经得到

[h '()δ-x x ]'

2πω=-s

x x 1J (s ω')-x x (29)

于是，我们得到

?g

()x 2ωπ

=s ''

'1()

()s x x x x g x J ω--?

2d 'x (30)