基于特殊窗函数的局部CT图像重建算法及其分析
合集下载
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第 18 卷 第 3 期 2009 年 9 月（1-9）
CT 理论与应用研究 CT Theory and Applications
文章编号：1004-4140（2009）03-0001-09
Vol.18, No.3 Sep., 2009
基于特殊窗函数的局部 CT图像重建算法及其分析
宋沛然，渠刚荣
（北京交通大学 理学院，北京 100044）
C( At) t3
dt
令，
R1
1 2 π2
π 0
d
t
l
Rf
(t
t0
,
)
Rflocal
(t
t0
,
)
r cos(
1 ) t0
t2
1 t2
dt
R2
1 4 π3A
π 0
d
t l
Rf
(t
t0 , )
Rflocal (t
t0 , )
C
Ar cos( ) t0 t r cos( ) t0 t 3
其中， k d0 ，d0 为径向采样步长，k 为边界外所取像素个数。我们用投影数据 Rflocal (t, ) 重建 x B(x0 , l) 的局部图像。于是有重建公式：
π
flocal (r, ) d Rflocal (t, ) qA r cos( ) t dt
0
（2）
3 误差分析
我们使用特殊窗函数[1]
FA
(l
)
0 A2 A2 l2
A2
e , A2 l2
0
,
l A l A
其中： 0 为系数归一化常数， l 为 ROI 半径， A 为窗函数带宽。特别地，当 A 1 时，
F1
(l
)
1
0 l
2
1
e
1 l 2
,
0
,
l 1 l 1
令 fglobal (r, ) 为全局重建函数； ferr (r, ) 为截断误差，则有：
ferr (r, ) ferr (r0 , )
π
d Rf (t, ) Rflocal (t, ) qA r cos( ) t qA (t0 t) dt
0
t t0 l
π
d Rf (t t0 , ) Rflocal (t t0 , ) qA r cos( ) t0 t qA (t) dt
1 2 π2
π 0
d
t l
r
cos(
1 )
t0
t 2
1 t2
2
dt
2
由文献[10]知，C1 (l, l ) 是个很小的数。若 x B(x0 , l) ，则 l r cos( ) t0 l 。于是，
c
R2
Rf Rflocal L 4π3 A
π
0 t l
1 (r cos( ) t0 t)3
Rf (t, ) f (l, ) dl L
（1）
重建一个二维函数 f (r, ) 。由式（1）把函数 f 映射到 Rf 的变换称为二维 Radon 变换。由 Rf 表示 f 的逆变换
1 2π2
π 0
t r cos( ) t
dt d
由此式可得到著名的 FBP 重建公式。为方便叙述，我们引入记号： D ，关于函数第一个变 量求偏导； H ，关于函数第一个变量做 Hilbert 变换； B ，反投影算子，则 Radon 逆变换 公式可写为：
0
t r0 cos( ) l
（3）
制作者(版权所有):《
》编辑部,
4
CT 理论与应用研究
当 A 1 时，由分部积分，
q1 (s)
1 2π2s2
C(s) 4π3s3
，
且有，
qA (s) A2 q1( As)
1 2 π2s2
C( As) 4 π3 A s3
A (s) 2
A
FA (l) sin(2 l s) dl 为所取正则化函数族； arc tan
0
r cos( ) h r sin( )
,
π 2
；
1
D r cos( )2 h r sin( )2 2 。
5 扇形束算法计算机实现
假设只在 的 M 个等间隔点处取投影值，其角度间隔为 ，且有 M 2π 成立；对于
4 扇形束重建
扇形束几何示意如图 2 所示。
图 2 扇形束几何 Fig.2 Fan-beam geometry
图 2 中， h 为光源到中心 O 距离；记 OP p 。有扇形束重建公式[12]：
f (r, )
1 4 π2
2π 0
1
1 D
g ( , )
1 D
g
(
,
)
sin( )
0
t l
因此，由（4）式，比较误差可写为：
1
2 π2
π
0
d
t
l
Rf
(t
t0
,
)
Rflocal
(t
t0
,
)
1
r cos( ) t0 t 2
1 t2
dt
1 4 π3 A
π 0
d
t l
Rf
(t
t0 , )
Rflocal (t
t0 , )
C
Ar cos( ) t0 t r cos( ) t0 t 3
1 t3
dt d
c
Rf Rflocal 4 π3 A
L
1
2
(
1
l)2
c
Rf Rflocal 2 π3 A 2
L
由上式可知 R2 依赖于 c ，因此依赖于选取的窗函数。此处 R2 是一个放大的估计，又因为式 （5）是震荡积分，由黎曼-勒贝格引理知，当 A 足够大时， R2 震荡趋于零。而在实际实现 时，取 A 1 d0 ，即窗的带宽是采样间隔的倒数。故比较误差 ferr (r, ) ferr (r0 , ) 的主要部 分由 R1 决定。因此，截断误差的主要表现是常数偏移。
关键词：局部重建；误差分析；扇形束；特殊窗函数
中图分类号：TP 391.41
文献标识码：A
在医学应用中，为了减少 X 射线对人体的损害，人们考虑局部感兴趣区域（Region of Interest，ROI）的图像重建。局部重建是使用 ROI 内的投影数据或增加少量 ROI 外的投影 数据，重建 ROI 的图像。2004 年，范惠荣等构造出了一种用于整体和局部图像重建的特殊 窗函数[1－2]。基于 Radon 变换的小波反演公式[3]，我们得到了图像重建的卷积反投影算法， 并验证了特殊选取的小波[3]可用于整体和局部图像重建。上述两种重建方法和本文中使用的 方法是近似的。
f
(r, )
1 2π
BHD Rf
(r, )
假设 f (r, ) 在 E2 上有紧支集（ E 为待重建区域半径），则 FBP 算法重建公式为：
π
E
f (r, ) d Rf (t, ) qA r cos( ) t dt
0
-E
A
其中， qA (s) 2 l FA (l) cos(2l s) dl ， FA (l) 为窗函数[12]。 0
因为，
ferr (r, ) fglobal (r, ) flocal (r, )
根据式（2），
π
fglobal (r, ) d Rf (t, ) qA r cos( ) t dt
0
π
ferr (r, ) d
Rf (t, ) Rflocal (t, ) qA r cos( ) t dt
C( At) t3
dt
制作者(版权所有):《
》编辑部,
3期
宋沛然等：基于特殊窗函数的局部 CT 图像重建算法及其分析
5
由 Schwartz 不等式， R1 C1 (l, l ) Rf Rflocal L2 ，其中：
1
C1(l, l )
2 局部重建
设 f (r, ) 在 l2 上 有 紧 支 集 ， 令 B(x0 , l) { x | x x0 | l } 是 待 重 建 局 部 区 域。一般，局部重建需用到局部区域外的部
分数据。我们根据 B(x0 , l ) 内的投影数据 Rf (t, ) 重建区域 B(x0 , l) 的图像（图 1）。
6
CT 理论与应用研究
18 卷
则有重建公式：
f (r, )
h 4 π2
2π 0
1 D2
q1( ) cos
q2 ( ) cos
g( ,
)
d
d
（7）
其中，q1(s) sA (s) sin2 s ，q2 (s) A (s) sA (s) sin s ，当 s 0 时，取 s 0 的极限值；
摘要：为了减少 X 射线对人体的伤害，局部图像重建成为人们研究的重点之一。本文研究了基 于一种特殊窗函数的局部图像重建算法。通常，直接用局部区域的投影数据重建局部图像，会 使重建图像产生常数偏移，我们将局部区域边缘的投影数据延拓为相应沿径向未知投影数据的 方法，改进了局部重建图像，并给出了基于特殊窗函数重建时，该方法的误差分析。另外，本 文将该算法应用于扇形束图像重建，并通过数值实验验证了此特殊窗函数应用于扇形束重建时 同样有效。
d d
其中，
g( , ) Rf ( p, ) Rf (h sin , )
（6）
为扇形束投影数据。假设对于 0 2π 和 （ 0 arcsin E h π 2 ）所有 g( , )
均已知，且 g( , ) 关于 和 的偏导数存在。
制作者(版权所有):《
》编辑部,
其中，
1
C(s) 3F1(l) l F1(l) sin(2 l s) dl 0 1 C(s) c 3F1(l) l F1(l) dl 0 1
C( As) 3F1(l) l F1(l) sin(2 l A s) dl 0
18 卷
（4） （5）
由式（3），比较误差为（为表述简便，令 t0 r0 cos( ) ）：
最近人们发展了精确局部重建算法。2004 年，Zou 和 Pan 建立了扇形束和锥形束反投 影滤波公式[4]，该公式对特殊的局部区域有精确重建；Noo 等基于滤波反投影算法的思想， 用经典的有限 Hilbert 反演公式导出了相应的局部重建算法[5]；张慧滔等将这些算法应用于 工业 CT[6]。这两种算法都依赖于特殊的局部重建区域。Ye[7－8]和 Yu 等[9]，指出已知待重建区 域中部分子区域的投影信息时，在锥形束扫描情况下可得到精确的重建图像，Yu 等给出了 上述算法的稳定格式重建算法。
每个投影，只在角度间隔为 的 2N 1 个等间隔角度处采样。所以，投影值 g n, m 是已
知的，其中 N n N, 0 m M 1 ， n 。由黎曼和逼近公式，计算式（7）的内层 积分离散为：
收稿日期：2009-04-03。 基金项目：国家自然科学基金（60772041）。
制作者(版权所有):《
》编辑部,
2
CT 理论与应用研究
18 卷
1 滤波反投影（Filtered Back Projection，FBP）算法重建公式
经典的计算机断层成像问题就是由投影值：
通常，只用局部区域内（或局部区域加部分局部区域外）的投影数据进行相应区域图 像的重建，会使得到的重建图像产生常数偏移[10]。根据已知的局部区域投影沿径向作常数 延拓[11]，可以部分地减少局部重建图像的常数偏移，我们将这个方法用于文献[1]中的算法， 给出了局部重建的误差分析，并用于相应的扇形束图像重建。数值实验表明，这个特殊窗 函数在扇形束重建中也同样能够得到很好的重建效果。
这 里 ， 对 于 任 意 Rf (t, ) ， 有 t r0 cos( ) l 。令 flocal (r, ) 表示局 部重建函数，则：
图 1 局部重建投影数据范围 Fig.1 The region of ROI projection
π
l r0 cos( )
flocal (r,) d
为未知投影区域的投影值[11]。
Rf (t, ), Rflocal (t, ) Rf (r0 cos( ) l , ),
r0 cos( ) (l ) t t r0 cos( ) (l )
r0 cos( ) l
Rf (r0 cos( ) l , ), t r0 cos( ) l
Rf (t, ) qA r cos( ) t dt
0
l r0 cos( )
制作者(版权所有):《
》编辑部,
3期
宋沛然等：基于特殊窗函数的局部 CT 图像重建算法及其分析
3
通常，直接用 B(x0 , l ) 内的投影数据重建局部区域图像，会使得到的重建图像产生 常数偏移[10]，影响重建效果。为减轻常数偏移，我们将已知投影区域边界上的投影值延拓
CT 理论与应用研究 CT Theory and Applications
文章编号：1004-4140（2009）03-0001-09
Vol.18, No.3 Sep., 2009
基于特殊窗函数的局部 CT图像重建算法及其分析
宋沛然，渠刚荣
（北京交通大学 理学院，北京 100044）
C( At) t3
dt
令，
R1
1 2 π2
π 0
d
t
l
Rf
(t
t0
,
)
Rflocal
(t
t0
,
)
r cos(
1 ) t0
t2
1 t2
dt
R2
1 4 π3A
π 0
d
t l
Rf
(t
t0 , )
Rflocal (t
t0 , )
C
Ar cos( ) t0 t r cos( ) t0 t 3
其中， k d0 ，d0 为径向采样步长，k 为边界外所取像素个数。我们用投影数据 Rflocal (t, ) 重建 x B(x0 , l) 的局部图像。于是有重建公式：
π
flocal (r, ) d Rflocal (t, ) qA r cos( ) t dt
0
（2）
3 误差分析
我们使用特殊窗函数[1]
FA
(l
)
0 A2 A2 l2
A2
e , A2 l2
0
,
l A l A
其中： 0 为系数归一化常数， l 为 ROI 半径， A 为窗函数带宽。特别地，当 A 1 时，
F1
(l
)
1
0 l
2
1
e
1 l 2
,
0
,
l 1 l 1
令 fglobal (r, ) 为全局重建函数； ferr (r, ) 为截断误差，则有：
ferr (r, ) ferr (r0 , )
π
d Rf (t, ) Rflocal (t, ) qA r cos( ) t qA (t0 t) dt
0
t t0 l
π
d Rf (t t0 , ) Rflocal (t t0 , ) qA r cos( ) t0 t qA (t) dt
1 2 π2
π 0
d
t l
r
cos(
1 )
t0
t 2
1 t2
2
dt
2
由文献[10]知，C1 (l, l ) 是个很小的数。若 x B(x0 , l) ，则 l r cos( ) t0 l 。于是，
c
R2
Rf Rflocal L 4π3 A
π
0 t l
1 (r cos( ) t0 t)3
Rf (t, ) f (l, ) dl L
（1）
重建一个二维函数 f (r, ) 。由式（1）把函数 f 映射到 Rf 的变换称为二维 Radon 变换。由 Rf 表示 f 的逆变换
1 2π2
π 0
t r cos( ) t
dt d
由此式可得到著名的 FBP 重建公式。为方便叙述，我们引入记号： D ，关于函数第一个变 量求偏导； H ，关于函数第一个变量做 Hilbert 变换； B ，反投影算子，则 Radon 逆变换 公式可写为：
0
t r0 cos( ) l
（3）
制作者(版权所有):《
》编辑部,
4
CT 理论与应用研究
当 A 1 时，由分部积分，
q1 (s)
1 2π2s2
C(s) 4π3s3
，
且有，
qA (s) A2 q1( As)
1 2 π2s2
C( As) 4 π3 A s3
A (s) 2
A
FA (l) sin(2 l s) dl 为所取正则化函数族； arc tan
0
r cos( ) h r sin( )
,
π 2
；
1
D r cos( )2 h r sin( )2 2 。
5 扇形束算法计算机实现
假设只在 的 M 个等间隔点处取投影值，其角度间隔为 ，且有 M 2π 成立；对于
4 扇形束重建
扇形束几何示意如图 2 所示。
图 2 扇形束几何 Fig.2 Fan-beam geometry
图 2 中， h 为光源到中心 O 距离；记 OP p 。有扇形束重建公式[12]：
f (r, )
1 4 π2
2π 0
1
1 D
g ( , )
1 D
g
(
,
)
sin( )
0
t l
因此，由（4）式，比较误差可写为：
1
2 π2
π
0
d
t
l
Rf
(t
t0
,
)
Rflocal
(t
t0
,
)
1
r cos( ) t0 t 2
1 t2
dt
1 4 π3 A
π 0
d
t l
Rf
(t
t0 , )
Rflocal (t
t0 , )
C
Ar cos( ) t0 t r cos( ) t0 t 3
1 t3
dt d
c
Rf Rflocal 4 π3 A
L
1
2
(
1
l)2
c
Rf Rflocal 2 π3 A 2
L
由上式可知 R2 依赖于 c ，因此依赖于选取的窗函数。此处 R2 是一个放大的估计，又因为式 （5）是震荡积分，由黎曼-勒贝格引理知，当 A 足够大时， R2 震荡趋于零。而在实际实现 时，取 A 1 d0 ，即窗的带宽是采样间隔的倒数。故比较误差 ferr (r, ) ferr (r0 , ) 的主要部 分由 R1 决定。因此，截断误差的主要表现是常数偏移。
关键词：局部重建；误差分析；扇形束；特殊窗函数
中图分类号：TP 391.41
文献标识码：A
在医学应用中，为了减少 X 射线对人体的损害，人们考虑局部感兴趣区域（Region of Interest，ROI）的图像重建。局部重建是使用 ROI 内的投影数据或增加少量 ROI 外的投影 数据，重建 ROI 的图像。2004 年，范惠荣等构造出了一种用于整体和局部图像重建的特殊 窗函数[1－2]。基于 Radon 变换的小波反演公式[3]，我们得到了图像重建的卷积反投影算法， 并验证了特殊选取的小波[3]可用于整体和局部图像重建。上述两种重建方法和本文中使用的 方法是近似的。
f
(r, )
1 2π
BHD Rf
(r, )
假设 f (r, ) 在 E2 上有紧支集（ E 为待重建区域半径），则 FBP 算法重建公式为：
π
E
f (r, ) d Rf (t, ) qA r cos( ) t dt
0
-E
A
其中， qA (s) 2 l FA (l) cos(2l s) dl ， FA (l) 为窗函数[12]。 0
因为，
ferr (r, ) fglobal (r, ) flocal (r, )
根据式（2），
π
fglobal (r, ) d Rf (t, ) qA r cos( ) t dt
0
π
ferr (r, ) d
Rf (t, ) Rflocal (t, ) qA r cos( ) t dt
C( At) t3
dt
制作者(版权所有):《
》编辑部,
3期
宋沛然等：基于特殊窗函数的局部 CT 图像重建算法及其分析
5
由 Schwartz 不等式， R1 C1 (l, l ) Rf Rflocal L2 ，其中：
1
C1(l, l )
2 局部重建
设 f (r, ) 在 l2 上 有 紧 支 集 ， 令 B(x0 , l) { x | x x0 | l } 是 待 重 建 局 部 区 域。一般，局部重建需用到局部区域外的部
分数据。我们根据 B(x0 , l ) 内的投影数据 Rf (t, ) 重建区域 B(x0 , l) 的图像（图 1）。
6
CT 理论与应用研究
18 卷
则有重建公式：
f (r, )
h 4 π2
2π 0
1 D2
q1( ) cos
q2 ( ) cos
g( ,
)
d
d
（7）
其中，q1(s) sA (s) sin2 s ，q2 (s) A (s) sA (s) sin s ，当 s 0 时，取 s 0 的极限值；
摘要：为了减少 X 射线对人体的伤害，局部图像重建成为人们研究的重点之一。本文研究了基 于一种特殊窗函数的局部图像重建算法。通常，直接用局部区域的投影数据重建局部图像，会 使重建图像产生常数偏移，我们将局部区域边缘的投影数据延拓为相应沿径向未知投影数据的 方法，改进了局部重建图像，并给出了基于特殊窗函数重建时，该方法的误差分析。另外，本 文将该算法应用于扇形束图像重建，并通过数值实验验证了此特殊窗函数应用于扇形束重建时 同样有效。
d d
其中，
g( , ) Rf ( p, ) Rf (h sin , )
（6）
为扇形束投影数据。假设对于 0 2π 和 （ 0 arcsin E h π 2 ）所有 g( , )
均已知，且 g( , ) 关于 和 的偏导数存在。
制作者(版权所有):《
》编辑部,
其中，
1
C(s) 3F1(l) l F1(l) sin(2 l s) dl 0 1 C(s) c 3F1(l) l F1(l) dl 0 1
C( As) 3F1(l) l F1(l) sin(2 l A s) dl 0
18 卷
（4） （5）
由式（3），比较误差为（为表述简便，令 t0 r0 cos( ) ）：
最近人们发展了精确局部重建算法。2004 年，Zou 和 Pan 建立了扇形束和锥形束反投 影滤波公式[4]，该公式对特殊的局部区域有精确重建；Noo 等基于滤波反投影算法的思想， 用经典的有限 Hilbert 反演公式导出了相应的局部重建算法[5]；张慧滔等将这些算法应用于 工业 CT[6]。这两种算法都依赖于特殊的局部重建区域。Ye[7－8]和 Yu 等[9]，指出已知待重建区 域中部分子区域的投影信息时，在锥形束扫描情况下可得到精确的重建图像，Yu 等给出了 上述算法的稳定格式重建算法。
每个投影，只在角度间隔为 的 2N 1 个等间隔角度处采样。所以，投影值 g n, m 是已
知的，其中 N n N, 0 m M 1 ， n 。由黎曼和逼近公式，计算式（7）的内层 积分离散为：
收稿日期：2009-04-03。 基金项目：国家自然科学基金（60772041）。
制作者(版权所有):《
》编辑部,
2
CT 理论与应用研究
18 卷
1 滤波反投影（Filtered Back Projection，FBP）算法重建公式
经典的计算机断层成像问题就是由投影值：
通常，只用局部区域内（或局部区域加部分局部区域外）的投影数据进行相应区域图 像的重建，会使得到的重建图像产生常数偏移[10]。根据已知的局部区域投影沿径向作常数 延拓[11]，可以部分地减少局部重建图像的常数偏移，我们将这个方法用于文献[1]中的算法， 给出了局部重建的误差分析，并用于相应的扇形束图像重建。数值实验表明，这个特殊窗 函数在扇形束重建中也同样能够得到很好的重建效果。
这 里 ， 对 于 任 意 Rf (t, ) ， 有 t r0 cos( ) l 。令 flocal (r, ) 表示局 部重建函数，则：
图 1 局部重建投影数据范围 Fig.1 The region of ROI projection
π
l r0 cos( )
flocal (r,) d
为未知投影区域的投影值[11]。
Rf (t, ), Rflocal (t, ) Rf (r0 cos( ) l , ),
r0 cos( ) (l ) t t r0 cos( ) (l )
r0 cos( ) l
Rf (r0 cos( ) l , ), t r0 cos( ) l
Rf (t, ) qA r cos( ) t dt
0
l r0 cos( )
制作者(版权所有):《
》编辑部,
3期
宋沛然等：基于特殊窗函数的局部 CT 图像重建算法及其分析
3
通常，直接用 B(x0 , l ) 内的投影数据重建局部区域图像，会使得到的重建图像产生 常数偏移[10]，影响重建效果。为减轻常数偏移，我们将已知投影区域边界上的投影值延拓