平行束反投影重建

一、 平行束反投影重建算法


平行束 重建采用的是平移加
平移

旋转的扫描方式，如图所示，射线源
在某一角度下 水平移动，将物体全部
X射线管

照射后旋转一角度，如此重复，在这
个过程中探测器相应地运动以接收
X射线。
1、反投影重建算法的物理概念：
断层平面中某一点的密度值可
以看作是这一平面内 所有经过该点
的射线的投影值之和（的均值）。
整幅重建图像可以看作是所有
方向下的投影累加而成。

探测器

平移

图 平行束平移加旋转

(3) (4)
射线标号示于图中，像素值（代
表密度）分别
x
1
,
x
2
,
x
3
,
x
4
，
赋值如下:
x
1
?5
,
x
2
?0
,
x
3
?2
,
x
4
?18

)为：
根据投影的定义(某条射线投影值为该条射线穿过的所
有的像素值之和)，每 条射线的投影
p
i
(
i?1,2
(1)
5
20
18
(5)
(6)
(7)
p
1
?x
1
?x
2
?5
,
p
2
?x
3
?x
4
?20
,
p
3
?x
1
?x
3
?7

p
4
?x
2
?x
4
?18
,
p
5
?x
3
?2
,
p
6
?x
1
?x
4
?23

(2)
p
7
?x
2
?0

根据反投影重建算法的物理意义，重建图像中各像素，
得到：
图 断层像素值和射线
x
1
?p
1
?p
3
?p
6
?35
，
x
2
?p
1
?p
4
?p
7
?23
，
x
3
?p
2
?p
3
?p
5
?29
，
x
4
?p
2
?p
4
?p
6
?61


5
2
0
18

35
29
26
61

5
4.1
3.3
8.7

(a) 原图像像素值

(b)反投影重建后图像

(c)求平均后图像

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

如图(c)所示。
因此有：
1
x
k
?
n
p
的第
i
条射线投影，
n
p
表示图像内的射线条数。
?
p
i?1
n
p
k,i

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

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

A
经反投影重建后，得到A点的像素值为
f
A
?(p
1
?p
2
?p
3
)3?1
。A点以外的像素 值原
来为零，经反投影重建后不等于零，而是等于13。所以，经反投影重建后的图像除保留A
点的像外，还有像素值为1n的灰雾背景，后者称为星状伪迹。产生星状伪迹的原因在于：
反投影的本质 是把取自有限物体空间的投影均匀地回抹（反投影）到射线所及的无限空间的
各点上，包括原先像素值为 零的点。

二、 反投影重建算法的数学描述
我们把“取投影” 、“反投影重建” 、“重建后图像” 这些环节看作是一个以原像为
输入，重建后图像为输出的成像系统，如图所示，先来求该系统的点扩展函 数PSF。
原像
?
(x,y)

y
r
取投影 反投影重建
重建后图像
h(x,y)
图 “反投影重建”成像系统
y

图中，
(x
r
,y
r
)
为旋 转坐标，
(x,y)
为固定坐标，
(r,
?
)
为极坐标，设 位于坐标原点
x?0,y?0
的点源
?
(x,y)
为
f(r ,
?
)
x?y
断面中唯一的像点，扫描方式为平移加旋转。即射
线先 平行移动，从物体的一侧移向另一侧；然后旋转一
个角度
?
?
，如此继续，直 到累计的旋转角达到
180
。
??
?
x
r
r
?
?
0
为止。为了计及这一几何要求我们设置一旋转坐标系统
x
x
r
?y
r
，它绕原点转动使投影总是沿着
y
r
方向 ，
x
r
?y
r
的原点与
x?y
的原点相同。二者的 夹角为
?
，不同的
?
代表不同的投射方向。投影线的位置可由
(x< br>r
,
?
)
完全确
定。
图平移加旋转扫描方式所用坐标系
设
?
为离散取值，如
?
?
?
i
，则相应的投影为：
p
?
i
( x
r
)?p(x
r
,
?
i
)?
?
??
??
??
??
f
r
(x
r
,y
r
)dy
r
?
?
?
(x
r
,y
r
)dy
r
??
??

?
?
?
(x
r
)
?
(y
r
)dy
r< br>?
?
(x
r
)|
?
?
?
i
?
?
(rcos(
?
?
?
i
))
其中x
r
?rcos(
?
?
?
)
，这可以从图的几 何关系很容易得出。
若
?
?
?
n
，则相应的投影为： < br>p
?
n
?
?
[rcos(
?
?
?< br>n
)]

根据反投影重建的定义式,点< br>(x
r
,y
r
)
的图像在所述坐标系统中表示为：
1
f(r,
?
)?f(x
r
,y
r
)?
N
?
?
1
N
?
1
p(x)|?
?
?
i
rx
r
?rcos(
?
?
?
i
)
N
?
i?1
i
N
?
?
p
?x
?
rcos(
?
?
?
)
?
i?1< br>i
N
?
i

?
?
p
?< br>?
rcos(
?
?
?
)
?
?
?i?1
i
其中，
N
?
为投影数，
?
?
?
?
N
?
。若在有限区间内将射线增至不相重的无限条，即连续取
投 影，则有：
f(r,
?
)?
1
?
?
?
0
p
?
?
rcos(
?
?
?
)
?< br>d
?

在忽略射束硬化的情况下，
?
在
(
?
,2
?
)
区间内的投影值等于
(0,
?
)
区间内的投影值。
在输入图像为点源的情况下，由及可得： < /p>

h(r,
?
)?
因为，
?
?
a(?
)
?
?
1
?
?
?
?
rco s(
?
?
?
)
?
d
?

0
?
?
(
?
?
?
0
)
( 此公式推导可参见附录一)，其中，
?
0
是
a(
?
)?0< br>的唯一的根。
a'(
?
0
)
令
a(
?)?rcos(
?
?
?
)
，则有：
a(
?0
)?rcos(
?
?
?
0
)?0

因为
r?0
，所以
sin(
?
?
?
0
) ?1
，于是有：
h(r,
?
)?
11
?
1
?
?
0
?
?
rcos(
?
?
?
)
?
d
?
?
1
?
?
0
?
?
(
?
?
?
0
)
d
?
rsin(
?
?
?
0
)
111
???
?
rs in(
?
?
?
0
)
?
r
?
x2
?y
2

可见，相应于反投影重建算法的系统， 它的点扩展函数不是
?
函数，系统不是完美的。
虽然在
r?0
处能反 映原图是点源的情况，但在
r?0
处，像素值
?0
，上式定量地描述了反投影重建算法星状伪迹的本质。
若原像为
?
(x,y)
，则将原像取投 影后再按反投影算法得到的重建图像为：
f(x,y)?
?
(x,y)??h(x,y )
，其中
??
表示二维卷积。
要去除反投影算法的星状伪迹，可以在输出端 加一滤波器，使加了滤波器后的反投影重
建等效成像系统的点扩展函数(PSF)为
?
(x,y)
。设滤波器的PSF为
q(x,y)
，相应的传递函数
为
Q(
?
,
?
)?F
2
?
q(x,y)
?< br>，要求(式的推导可参见附录一)：
1
?
x?y
做二维傅氏变换，得 ：
22
??q(x,y)?
?
(x,y)


1
??
?
?
22
Q(
?
,
?
)? 1

或：
Q(
?
,
?
)?
??
2
?
?
2
?
??
< br>这是一个二维滤波器，实现起来较麻烦。若
?
的变化范围扩至
?
，则根 本不能实现。
但它提供了去除星状伪迹的一个方向。

三、 滤波反投影重建算法
反投影重建算法的缺点是引入星状伪迹，即原来图像中密度为零的点，重建后不一 定为
零，从而使图像失真。去除伪迹的方法如图所示。
输入图像
（原图像）
取投影 反投影重建
(a)
二维滤波器
??
输出图像
（同原图像）
输入图像
（原图像）
取投影 一维滤波器
(b)
反投影重建
输出图像
（同原图像）
(a)方案之一，先“反投影”再二维滤波；(b) 方案之二，先对投影数据作一维滤波，再“反
投影”

图 星状伪迹的去除
投影定理
投影定理（或中心切片定理）：
在非衍 射源情况下，某图像
f(x,y)
在视角为
?
时投影
p
?< br>(x
r
)
的一维傅氏变换给出
f(x,y)
的二维傅氏变换< br>A(
?
1
,
?
2
)?A(
?
,?
)
的一个切片。切片与
?
1
轴相交成
?
角， 且通过坐
标原点。即：
^
^
F
?
p(x
r
)
?
?
?A(
?
,
?
)|
?
固 定

1
?
?
一
p?
(x
r
)
x
r
维
傅
氏
变< br>换
A(
?
1
,
?
2
)
y
y
r
?
2
?
?
x
r
二维傅立
?
?
0
x
?
1
叶变换

图阐明投影定理
证明：
f(x,y)
的二维傅立叶变换为
A(
?
1
,
?
2
)?
?
由图可以得到 如下几何关系

??
??
?
??
??
f(x,y )e
-j(
?
1
x?
?
2
y)
dxdy< br>
?
x
r
?
?
cos
?

??
?
?
?
y
r
?
?
?sin< br>?
sin
?
??
x
?
?
y
?

cos
?
?
???

以
y
r
为投影轴的投影为：
p
?
(x
r< br>)?
?
它的一维傅里叶变换为：
??
??
f(x,y)dy
r
?
?
??
??
f
r
(x
r,y
r
)dy
r

F
1
[p
?
(x
r
)]?
?
?
?
??
??
??
p
?
(x
r
)e
-j2??
x
r
dx
r
??
??
??
??
??
??
??
??
f
r
(x
r,y
r
)e
-j2
??
x
r
dx
r< br>dy
r
f(x,y)e
-j2
??
(xcos
??ysin
?
)
dxdy

??
式中采用
2
??
可以使后续反变换计算中的系数得以简化。
在推导上式的过程中，我们利用了下列关系：
f(x,y)?f
r
(x
r
,y
r
)

以及
?x
r
?x
dx
r
dy
r
?
?y
r
?x
由式可知，当
?
1
,
?2
的取值满足条件
?x
r
cos
?
?y
?< br>?y
r
?sin
?
?y
sin
?
dxdy? dxdy

cos
?
?
1
?2
? ?
cos
?
?
?

?
2
?2
??
sin
?
?
时，有
F
1
[p
?
(x
r
)]?
?

?
????
?
?
f(x,y)e
-j(
?
1
x?
?
2
y)
dxdy|
?
1
?2??
cos
?
?
2
?2
??
sin
?
?A(
?
1
,
?
2
)|
?
1?2
??
cos
?
?
2
?2
??
si n
?

上式说明：
?
1
,
?< br>2
不是独立的而是受到式的约束，其值必须局限于直线
?
2
?(tg< br>?
)
?
1
上。
根据投影定理，投影图像重建问题原则上可按如下流程求解：
（1）采集不同视角下的投影；
（2）求出各投影的一维傅氏变换，此即图像二维傅氏变换过原点的各切片，理论上是连
续的无 穷多片；
（3）将上述各切片汇集成图像的二维傅氏变换；
（4）求二维傅氏反变换的重建图像。

卷积反投影重建（即滤波反投影） 待重建图像为
a(x,y)
，它的二维傅氏变换为
A(
?
1,
?
2
)?A(
?
,
?
)
。根据中心 切片定理，
^
A(
?
,
?
)
可通过
a(x ,y)
在不同视角
?
下的投影
p
?
(x
r
)
的一维傅氏变换求得。即：

待建图像：
^
A(
?
1
,
?
2
)?A(
?
,
?
)?F
1
?
?
(
?
)?P(
?
,
?)
?
p
?
(x
r
)
?
?
?P
^

a(r,
?
)?a(x,y)?F
2
?1
[
A(
?
1
,
?
2
)
]
?
1
4
?
?
0
2
^
??
?
^
?
??
????
A(
?
1
,
?
2
)e
i(
?
1
x?
?
2
y)
d
?
1
d
?
2

?
??
?
?
??
0
??
A(
?
,
?
)e
i2
??
rcos(
?
?
?
)?
d
?
d
?
??
P(
?
,
?
)e
i2
??
rcos(
?
?
?
)
?
d
?
d
?
?
??
?
?
d?
?
0
?
?
P(
?
,
?
)e
i2
??
rcos(
?
?
?
)
d
?
因为
x
r
?rcos(
?
?
?
)
，所以有：
?
1
x?
?
2
y?2
??
(xcos
?
?ysin
?
)?2
??
x
r
?2
??
rcos(
?
?
?
)

同时：
d
?
1
d
?
2
?Jd
?
d
?

J?
先来看该式的第二个积分：
?
?
1
?
?
?
?
2
?
?
?
?
1
?
?
2
?
cos
?
?
?
?
2
?
?
2
?
sin
?
?2
??
sin?
?4
?
2
?

2
??cos
?
?
?
??
?
P(
?
,
?
)e
i2
??
rcos(
?
?
?
)< br>d
?
?
?
?
P(
?
,
?
) e
i2
??
x
|
x?rcos(
?
?
?< br>)
d
?
r
?
??
r
?h(x
r)?p(x
r
,
?
)|
x
r
?rcos(?
?
?
)
?g(x
r
,
?
)|
x
r
?rcos(
?
?
?
)
?g
?rcos(
?
?
?
),
?
?
式中：
式的物理意义是投影
p(x
r
,
?
)
经过传递函数为
?
?F
1
[h(x
r
)]
的滤波器后得到的修正后
的投影
g(x
r
,
?
)
在满足
x
r?rcos(
?
?
?
)
时的值。将代入，得到：
a( r,
?
)?
?
g[rcos(
?
?
?
),
?
]d
?

0
^


g(x
r
,
?)?h(x
r
)?p(x
r
,
?
)

?
称为滤波反投影方程，其物理意义是经过给定点
(r,
?
)
的所有滤波后的投影在
?
?0
~
?
范围
内的累加—反投影 重建，得出
(r,
?
)
点的像素值。
可见，滤波（卷积）反投影算法的具体包含三大步：
(1) 把在固定视角
?
i
下测得的投影
p(x
r
,
?
)
经过滤波，得到 滤波后的投影
g(x
r
,
?
)
；
(2) 对每一 个
?
i
，把
g(x
r
,
?
i
)< br>反投影于满足
x
r
?rcos(
?
?
?
i< br>)
的射线上的所有各点
(r,
?
)
；
(3) 将步 骤(2)中的反投影值对所有
0
<
?
?
?
进行累加（积分） ，得到重建后的图像。

滤波函数与内插函数的选取原则
由式,可知， 先要把投影值
p(x
r
,
?
)
做一维滤波。滤波器的传递函 数为
H(
?
)?
?
。可以
证明这一理想的滤波函数是不可实 现的。我们的目的是：选取这样的滤波函数，使它既可以
实现，又有足够的精度。为此，结合成像的具体 情况考虑。
(1) 投影数据的高频（空间分辨率）分量幅度很小；
(2) 投影数据的采集是天然离散的；

(3) 存在噪声。
在物体尺寸有限的情况下，投影数据的空间分布是有限的，因而严格来说，其频带无限。
但若物体密度在 空间变化是平稳的，则高频分量幅度不大。若采样间隔
d
足够小,则在折叠
频率
12d
以上的频率分量可以不计。即：
P(
?
)?0
由于投影数据在空间上的天然离散性，我们有：
?
?
1

2d

p(x
r
,
?
)|
x
r
?rcos(?
?
?
)
?
?
p
?
nd,
?
?
?

其中，
d
为射 束平移的步距；
n
为正整数；
?
为某一固定的视角；
?
?< br>?
表示序列，与此相应，
滤波函数也取离散形式：
r
?
h( nd)
?
?h(x
r
)|
x?nd

滤波后的投影值为：
?
p(nd)
?
?
?
p?
nd
?
?
*
?
h(nd)
?

由于实际中物体是连续的，即
x
r
连续，但是由于采样是离散的， 不一定能 够正好落在
采样点位置上，为了得到任意一点的值，需要对式得到的值进行内插。内插函数为
? (x
r
)
，
则经内插后的任意一点滤波后投影值为：

p(x
r
)?{p(nd)}*{h(nd)}*?(x
r
)

如前所述，投影值的高频分量小，在折叠频率以上可以忽略不计，所以可以把滤波器强
制为带限 的而不影响结果，即：
H(
?
)?
?
W(
?
)< br>?
1
?
?
W(
?
)?
?
?
0
?
?
?
?
1
2d

1
?
?
2d
结合式与式，可得：
p(x
r
)?{p(nd)}*{h(nd)}*?(x
r
)?h(x
r
)?p(x
r
)

对式两边进行傅里叶变换得：
?
F
?
?
p(x
r
)
?
?F
?
{p( nd)}
?
F
?
{h(nd)}
?
?(
?
)
?H
s
(
?
)P
s
(
?
)?(
?
)=H(
?
)P(
?
)
其中：

??
k
??
1
??
P
s
(
?)?F
?
p(x
r
)
?
?
(x
r?nd)
?
?
?
P(
?
?)
d
??< br>??
d
??

????
k
? ?
1
H
s
(
?
)?F
?
h(x
r
)
?
?
(x
r
?nd)
?
?
?< br>H(
?
?)
d
??
??
d
??
即时 域中采样对应频域中平移再叠加。
为了使式成立，需要注意以下问题：
(1)
P
s
(
?
)
为各个相隔
1d
的
P(
?
)
的叠加（幅度乘一系数）。由于
P(
?
)
在
?
?12d
时
不是完全消失，只是幅值很小。叠加的结果，必有混叠效应，使在
?
接近
12d
处的
P
s
(
?
)