计算机断层重建(CT)是一个比较热门的领域,这篇文章简单介绍了反投影方法的重建过程。
参考资料:冈萨雷斯,《数字图像处理》,电子工业出版社。
直接反投影
该方法是沿着射线来的方向把一维信号反投影回去,可以想象成把投影穿过图像区域反“涂抹”回去。注意到相隔180°的投影互为镜像,因此,为了产生重建所要求的所有投影,只需要按照圆周的一半角度增量来考虑即可。下图是不同投影角度数量下直接反投影法的结果
从图中可以看到,直接反投影法有明显的伪影,也被称之为“晕环”。伪影的亮度比物体第,但比背景高。越多的投影使得伪影亮度越低,但会带来模糊。在CT重建中,去除伪影是一个重要的问题。
投影与Radon变换
Radon变换描述了投影产生时的数学过程。考虑直线
y
=
a
x
+
b
y=ax+b
y=ax+b,将其用极坐标来表示可写为:
x
cos
θ
+
y
sin
θ
=
ρ
x\cos\theta+y\sin\theta=\rho
xcosθ+ysinθ=ρ
对于图像
f
(
x
,
y
)
f(x,y)
f(x,y),沿该线积分为:
g
(
ρ
,
θ
)
=
∫
−
∞
∞
∫
−
∞
∞
f
(
x
,
y
)
δ
(
x
cos
θ
+
y
sin
θ
−
ρ
)
d
x
d
y
g(\rho,\theta)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)\delta(x\cos\theta+y\sin\theta-\rho)dxdy
g(ρ,θ)=∫−∞∞∫−∞∞f(x,y)δ(xcosθ+ysinθ−ρ)dxdy
其中
δ
\delta
δ为示性函数,或者是和示性函数差不多的冲激函数。离散情况下可写为:
g
(
ρ
,
θ
)
=
∑
0
M
−
1
∑
0
N
−
1
f
(
x
,
y
)
δ
(
x
cos
θ
+
y
sin
θ
−
ρ
)
g(\rho,\theta)=\sum_{0}^{M-1}\sum_{0}^{N-1}f(x,y)\delta(x\cos\theta+y\sin\theta-\rho)
g(ρ,θ)=0∑M−10∑N−1f(x,y)δ(xcosθ+ysinθ−ρ)
进一步可以得到直接反投影的数学表示:
f
(
x
,
y
)
=
∫
0
π
g
(
ρ
,
θ
)
d
θ
f(x,y)=\int_0^{\pi}g(\rho,\theta)d\theta
f(x,y)=∫0πg(ρ,θ)dθ
滤波反投影法(FBP)
傅里叶切片定理
关于
ρ
\rho
ρ投影的一维傅里叶变换为
G
(
ω
,
θ
)
=
∫
−
∞
∞
g
(
ρ
,
θ
)
e
−
j
2
π
ω
ρ
d
ρ
G(\omega,\theta)=\int_{-\infty}^{\infty}g(\rho,\theta)e^{-j2\pi\omega\rho}d\rho
G(ω,θ)=∫−∞∞g(ρ,θ)e−j2πωρdρ
傅里叶变换可以看成从空间域到频率域上的转换,这里
ω
\omega
ω为频率变量。将
g
g
g带入上式中可得:
G
(
ω
,
θ
)
=
∫
−
∞
∞
∫
−
∞
∞
∫
−
∞
∞
f
(
x
,
y
)
δ
(
x
cos
θ
+
y
sin
θ
−
ρ
)
e
−
j
2
π
ω
ρ
d
x
d
y
d
ρ
=
∫
−
∞
∞
∫
−
∞
∞
f
(
x
,
y
)
[
∫
−
∞
∞
δ
(
x
cos
θ
+
y
sin
θ
−
ρ
)
e
−
j
2
π
ω
ρ
d
ρ
]
d
x
d
y
=
∫
−
∞
∞
∫
−
∞
∞
f
(
x
,
y
)
e
−
j
2
π
ω
(
x
cos
θ
+
y
sin
θ
)
d
x
d
y
=
F
(
ω
cos
θ
,
ω
sin
θ
)
\begin{aligned} G(\omega,\theta)&=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)\delta(x\cos\theta+y\sin\theta-\rho)e^{-j2\pi\omega\rho}dxdyd\rho \\ &=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)\left[\int_{-\infty}^{\infty}\delta(x\cos\theta+y\sin\theta-\rho)e^{-j2\pi\omega\rho}d\rho\right]dxdy \\ &=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)e^{-j2\pi\omega(x\cos\theta+y\sin\theta)}dxdy \\ &=F(\omega\cos\theta,\omega\sin\theta) \end{aligned}
G(ω,θ)=∫−∞∞∫−∞∞∫−∞∞f(x,y)δ(xcosθ+ysinθ−ρ)e−j2πωρdxdydρ=∫−∞∞∫−∞∞f(x,y)[∫−∞∞δ(xcosθ+ysinθ−ρ)e−j2πωρdρ]dxdy=∫−∞∞∫−∞∞f(x,y)e−j2πω(xcosθ+ysinθ)dxdy=F(ωcosθ,ωsinθ)
其中
F
(
u
,
v
)
F(u,v)
F(u,v)为
f
(
x
,
y
)
f(x,y)
f(x,y)的二维傅里叶变换,即:
F
(
u
,
v
)
=
∫
−
∞
∞
∫
−
∞
∞
f
(
x
,
y
)
e
−
j
2
π
(
u
x
+
v
y
)
d
x
d
y
F(u,v)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)e^{-j2\pi(ux+vy)}dxdy
F(u,v)=∫−∞∞∫−∞∞f(x,y)e−j2π(ux+vy)dxdy
这说明了,对投影 g g g的一维傅里叶变换,等于对图像 f f f二维傅里叶变换在某一条线下的切片。这就是傅里叶切片定理。
平行射线下的滤波反投影重建
傅里叶逆变换:
f
(
x
,
y
)
=
∫
−
∞
∞
∫
−
∞
∞
F
(
u
,
v
)
e
j
2
π
(
u
x
+
v
y
)
d
u
d
v
f(x,y)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}F(u,v)e^{j2\pi(ux+vy)}dudv
f(x,y)=∫−∞∞∫−∞∞F(u,v)ej2π(ux+vy)dudv
令
u
=
ω
cos
θ
,
v
=
ω
sin
θ
u=\omega\cos\theta,v=\omega\sin\theta
u=ωcosθ,v=ωsinθ,则有
f
(
x
,
y
)
=
∫
0
2
π
∫
−
∞
∞
F
(
ω
cos
θ
,
ω
sin
θ
)
e
j
2
π
ω
(
x
cos
θ
+
y
sin
θ
)
ω
d
ω
d
θ
f(x,y)=\int_{0}^{2\pi}\int_{-\infty}^{\infty}F(\omega\cos\theta,\omega\sin\theta)e^{j2\pi\omega(x\cos\theta+y\sin\theta)}\omega d\omega d\theta
f(x,y)=∫02π∫−∞∞F(ωcosθ,ωsinθ)ej2πω(xcosθ+ysinθ)ωdωdθ
应用傅里叶切片定理,以及
G
(
ω
,
θ
+
π
)
=
G
(
−
ω
,
θ
)
G(\omega,\theta+\pi)=G(-\omega,\theta)
G(ω,θ+π)=G(−ω,θ)可得
f
(
x
,
y
)
=
∫
0
2
π
∫
−
∞
∞
G
(
ω
,
θ
)
e
j
2
π
ω
(
x
cos
θ
+
y
sin
θ
)
ω
d
ω
d
θ
=
∫
0
π
[
∫
−
∞
∞
∣
w
∣
G
(
ω
,
θ
)
e
j
2
π
ω
ρ
d
ω
]
ρ
=
x
cos
θ
+
y
sin
θ
d
θ
\begin{aligned} f(x,y)&=\int_{0}^{2\pi}\int_{-\infty}^{\infty}G(\omega,\theta)e^{j2\pi\omega(x\cos\theta+y\sin\theta)}\omega d\omega d\theta \\ &=\int_{0}^{\pi}\left[\int_{-\infty}^{\infty}|w|G(\omega,\theta)e^{j2\pi\omega\rho}d\omega\right]_{\rho=x\cos\theta+y\sin\theta}d\theta \end{aligned}
f(x,y)=∫02π∫−∞∞G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ=∫0π[∫−∞∞∣w∣G(ω,θ)ej2πωρdω]ρ=xcosθ+ysinθdθ
方括号里可以看成一维傅里叶逆变换,只是附加了 ∣ ω ∣ |\omega| ∣ω∣项,这一项可以看为一个斜坡滤波器。因为 ∣ ω ∣ |\omega| ∣ω∣在无穷处都到达 + ∞ +\infty +∞,所以傅里叶反变换无定义。实践中常常通过加窗的方式,使它在定义的范围之外为0。
简单的方法是通过方波限制,但方波有不希望的振铃特性,这会为图像重建带来伪影。因此可以考虑平滑窗,比如汉明窗或者韩窗:
h
(
w
)
=
{
c
+
(
c
−
1
)
cos
2
π
ω
M
−
1
0
⩽
ω
⩽
(
M
−
1
)
0
其他
h(w)=\left\{ \begin{aligned} &c+(c-1)\cos\frac{2\pi\omega}{M-1} \qquad &0\leqslant\omega\leqslant(M-1) \\ &0 \qquad &其他 \end{aligned}\right.
h(w)=⎩
⎨
⎧c+(c−1)cosM−12πω00⩽ω⩽(M−1)其他
总结: 获得完整反投影图像是由如下步骤得到的:
- 计算每个投影的一维傅里叶变换。
- 用滤波函数 ∣ ω ∣ |\omega| ∣ω∣乘以每个傅里叶变换,即乘以一个合适的窗。
- 得到每个滤波后的变换的一维傅里叶反变换。
- 对所有的一维反变换积分(求和)。
卷积与傅里叶反变换
考虑傅里叶反变换可以写成卷积的形式,可以有
f
(
x
,
y
)
=
∫
0
π
[
∣
ω
∣
G
(
ω
,
θ
)
e
j
2
π
ω
ρ
d
ω
]
ρ
=
x
cos
θ
+
y
sin
θ
d
θ
=
∫
0
π
[
s
(
ρ
)
⋆
g
(
ρ
,
θ
)
]
ρ
=
x
cos
θ
+
y
sin
θ
d
θ
=
∫
0
π
[
∫
−
∞
∞
g
(
ρ
,
θ
)
s
(
x
cos
θ
+
y
sin
θ
−
ρ
)
d
ρ
]
d
θ
\begin{aligned} f(x,y)&=\int_0^{\pi}\left[|\omega|G(\omega,\theta)e^{j2\pi\omega\rho}d\omega\right]_{\rho=x\cos\theta+y\sin\theta}d\theta \\ &=\int_0^{\pi}\left[s(\rho)\star g(\rho,\theta)\right]_{\rho=x\cos\theta+y\sin\theta}d\theta \\ &=\int_0^{\pi}\left[\int_{-\infty}^{\infty}g(\rho,\theta)s(x\cos\theta+y\sin\theta-\rho)d\rho\right]d\theta \end{aligned}
f(x,y)=∫0π[∣ω∣G(ω,θ)ej2πωρdω]ρ=xcosθ+ysinθdθ=∫0π[s(ρ)⋆g(ρ,θ)]ρ=xcosθ+ysinθdθ=∫0π[∫−∞∞g(ρ,θ)s(xcosθ+ysinθ−ρ)dρ]dθ
其中
⋆
\star
⋆表示卷积,
s
(
ρ
)
s(\rho)
s(ρ)表示
∣
ω
∣
|\omega|
∣ω∣的傅里叶逆变换。
注:
- 重建过程中不需要存储所有的反投影图像,单个求和运算仅被最后的反投影图像更新。
- 斜坡滤波器在频率域的直流项归零,故每幅反投影图像的均值将为零。这意味着每幅反投影图像都将有负像素和正像素。简单的解决办法是所有像素减去最小值再放缩。
f ′ = K [ f − min f max ( f − m i n f ) ] f'=K\left[\frac{f-\min{f}}{\max{(f-min{f})}}\right] f′=K[max(f−minf)f−minf]