Cone-Beam Geometry X
Z
Y
螺旋锥形束 CT
(From G,Wang)
心脏成像动态成像
4D CT
4D CT
4D CT 检测器
Multi-source cone-beam CT
Detector
Source
(From G,Wang)
小动物成像 Micro-CT
Completeness Conditions for Cone-
Beam CT
Completeness Condition,If on every plane
that intersects the object there exists at least
one cone-beam source point,then one can
reconstruct the object,(B.D.Smith1985)
The possible structures satisfying the conditions are:
Works on Cone-Beam CT
Kirillov A.A,(1961)—Completeness condition
Tuy H.K,(1983)-- Completeness condition
Feldkamp L.A.(1984)---Practical,Incomplete projection
Smith B.D.(1985)-- Completeness condition
Grangeat P.(1987)—Reconstruction Algorithm
Ge Wang,(1991)—Reconstruction Algorithm
Defrise M.Clark R.(1994)-- Reconstruction Algorithm
Kudo H,Saito T.(1991,1994)– Reconstruction Algorithm
CT重建算法
近似重建算法
迭代重建算法
精确重建算法
Tam-Danielsson Geometry
Tam,Three-dimensional computerized tomography
scanning method and system for large objects with
smaller area detector,US Patent 5,390,112,1995
Detector surface is limited to
a cylindrical section between
two consecutive helical turns
Every point is on one and
only one PI-line
Pi-Line
Perfect Mosaic
Tam,Samarasekera,Sauer,Exact cone-beam CT with a spiral scan,Phys,Med,Biol,43:847-855,1998
Kudo,Noo,Defrise,Cone-beam filtered-backprojection algorithm for truncated helical data,
Phys,Med,Biol,43:2885-2909,1998
(From G,Wang)
Katsevich Theorem (2002)
dsdxsqyDqsyxxf sqf
xI PI
s i n
1|)),,(),((
|)(|
1
2
1)( 2
0)(
),(),(),( xsuxsxse ),(s in),(c o s),,( xsexsxs
)(0sy
e
),( xsu
),(0xs?
)(1sy
)(2sy
)(xf
Detector
Plate
Source
Object
Pi-Line
(From G,Wang)
Noo’s Short-Scan (2002)
Integral on any line through a region
(From G,Wang)
我国 CT的发展
1) 85年第二代 CT,上海计算技术研究所等
2) 安科公司
3) Siemens
4) GE
5)中国人有才智、有能力造好 CT。
投影投影投影投影投影投影投影投影
月食、日食
日晷(计时)
几何投影,坐标分量(位臵)
魔术
电影
猴子捞月
杯弓蛇影
顾影自怜
希腊美少年?Narcissus( 水仙 )? narcissi( 自恋 )
投影投影投影经典断层成像投影谁先匹配准,奖励 20万通缉令投影投影反投影重建算法的一般步骤:
原像取投影 反投影重建重建后图像反投影重建算法投影重建算法的基本内容:
“断层平面中某一点的密度值可看作这一平面内所有经过该点的射线投影之和(的平均值)”。
式中 xk表示象素的值,pk,i为经过象素的第 i条射线投影。
可以这样理解:
其中 f(x,y,z)是身体组织密度。
,
1
1 pn
k k i
p i
xp
n?
(,,)lp f x y z d l
反投影重建算法的物理概念算法举例
1
1 3 7 1 1 1 5 2p x x x x
2
2 2 6 1 0 1 4 6p x x x x3
3 5 6 7 8 7p x x x x4
4 9 1 0 1 1 1 2 1p x x x x
5
5 1 6 1 1 1 6 5p x x x x
6
6 4 7 1 0 1 3 3p x x x x
算法举例
根据反投影算法
x1=p5 = 5
x6=p2+p3+p5=18
…
平均化处理,除以投影线数目
xi=xi/6
0 0 0 0
0 5 2 0
0 1 0 0
0 0 0 0
5 6 2 3
7 18 12 7
1 10 8 1
3 6 2 5
0.83 1 0.33 00.5
1.16 3 2 1.16
0.06 1.66 1.33 0.16
0.5 1 0.33 0.83
反投影重建后原像素值再除以投影线数,平均化断层平面中某一点的密度值可看作这一平面内所有经过该点的射线投影之和的平均值
12
3
4
5
6
星状伪迹
反投影重建后,原来为 0的点不再为 0,形成伪迹
0 0 0 0
0 5 2 0
0 1 0 0
0 0 0 0
0.83 1 0.33 0.5
1.16 3 2 1.16
0.06 1.66 1.33 0.16
0.5 1 0.33 0.83
原像素值 再除以投影线数,平均化星状伪迹
我们考虑孤立点源反投影重建,
中心点 A经 n条投影线投影后,投影值均为 1:
p1=p2=...=pn=1
因此重建后而其他点均为 1/n
这类伪迹成为星状伪迹
1/n 1/n 1/n
1/n 1 1/n
1/n 1/n 1/n
12
1 (,,,) 1
Anf p p pn
0 0 0
0 1 0
0 0 0
星状伪迹
产生星状伪迹的原因在于:反投影重建的本质是把取自有限物体空间的射线投影均匀地回抹
(反投影 )到射线所及的无限空间的各点之上,
包括原先像素值为零的点反投影重建算法的数学描述我们设臵一旋转坐标系统 Xr-Yr,它绕原点转动使投影线总是沿着 Yr方向。 Xr-Yr的原点与 X-Y的原点相重。两者的夹角为 φ,不同的 φ代表不同的投射方向。投影线的位臵可由
(Xr,φ)完全确定。 空间的任一点的位臵可用 (X,Y),(Xr,Yr)或极坐标 (r,θ)表出。
反投影重建算法的数学描述
先设 φ为离散取值,则投影为
根据反投影重建的定义式,点 的图像在所述坐标系统表示为式中,,为投影数
1 1( ) ( ) ( ) | c o sr r r rx y d y x r
1
1 1( ) (,) (,) (,)r r r r r r r r rp x p x f x y d y x y d y
,rrxy
1
11
1
(,),
11
c os c os
i
ii
N
r r r r
i
NN
ii
ii
f r f x y p x
N
p r p r
N
/ N N?
若在有限区间内射线增至不相重的无限条,即连续投影,
则式 (4.12)过渡到更一般的连续情况下的反投影表达式:
在输入图像为点源的情况下,由式 (4.10)及式 (4.13)可得
0
1,c o sf r p r d?
0
00
0
22
0
11
,c o s
s in
1 1 1 1 1
s in
h r r d d
r
rr xy
反投影重建算法的数学描述
可见,反投影重建算法相应的系统的扩展函数不是 δ函数,系统不是完美的.式定量地描述了反投影重建算法星状伟迹的本质.
要除去反投影算法的星状伪迹,我们可以在输出端加一滤波器,使加了滤波器后的反投影重建成像系统 PSF
= δ(x,y).使滤波器的 PSF为 q(x,y),相应的传递函数为
Q(ξ,η),这里
我们要求:
2,,Q F q x y
22
1 * *,(,)q x y x y
xy
反投影重建算法的数学描述
**表示二维卷积,对上式取二维傅立叶变换得
或
这是一只二维滤波器,实现起来比较麻烦.若
ρ的变换范围可扩至 ∞,根本不能实现,但不管怎样,它提供了去除星状伪迹的一个努力方向.
221,1Q
22,Q
反投影重建算法的数学描述
反投影重建算法要求对任一点,找出经过该点的所有射线,从而求得过该点的射线投影的均值。
从图中可以看出,过给定点的射线 都在以 r为直径,以 为圆心的圆周上,即正弦图
(,)r?
(,)rx?
c o s ( )rxr
1(,)
2r?
正弦图正弦图
在 空间中,在给定 情况下,
我们得到的是一条正弦曲线。图像中的每一点,对应于空间中的一条正弦曲线,
全部图像对应于一簇正弦曲线族,这就是正弦图。
(,)rx? (,)r?
正弦图反投影重建算法的实现
正弦图中,射线位臵用 描述,扫描操作由角增量
Δ,步距 d步进;若正弦图的 Δ,d不是足够小,图中的某些点 可能找不到一点( mΔ,nd)与其对应;
若令 Фi= m Δ,那么,
一般地,故不能直接得到,须内插。
1
0
(,) 1 / [ c o s( ) ]mM
m
f r N p m r m
c o s ( )r m n d [ c o s ( ) ]mp r m
(,)rx?
(,)r?
内插内插方法
1.紧邻内插,在 时,若则以 中的内容 直接代替
。 否则,用 代替
2.线形内插:如果,则
[ c o s ( ) ]mp r m
( ) { [ ( 1 ) ] ( ) } * [ c o s ( ) ] /m m mp n d p n d p n d r m n d d
m m
c o s ( ) c o s ( ) ( 1 )r m n d r m n d
(,)m nd?
[ c o s ( ) ]mp r m
()mp nd?
[ ( 1 ) ]mp n d
[ c o s ( ) ]mp r m
c os( ) ( 1 )nd r m n d
滤波反投影算法
去除星状伪迹的两种方法:
取投影取投影 一维滤波器 反投影重建二维滤波器反投影重建输入图像(原图像)
输入图像
(原图像)
输出图像
(同原图像 )
输出图像
(同原图像 )
图 4.18 星状尾迹的去除
( b)
( a)
第一种方法中的二维滤波器较难实现。
第二种方法交换了滤波与反投影的顺序,
此时投影数据是一维的,易于实现。此方法的理论基础就是我们即将提到的投影定理(中心切片定理)
滤波反投影算法投影定理
某一图像 在视角 时投影 的 1D傅立叶变换给出 的 2D傅立叶变换的一个切片,切片与 轴相交成 角,且通过坐标原点。即
投影定理的 证明,p138- 139
1?[ ( ) ] (,) |rF p x F 固 定
12?(,) (,)FF(,)f x y
(,)f x y
()rpx
1?
投影定理卷积反投影重建
待建图像为 a(x,y),它的 2D傅氏变换为,
由中心切片定理:
12?(,) (,)AA
12?(,) (,) [ ( ) ]rA A F p x
( ) (,)PP
卷积反投影重建
待建图像
其中
12 1 2? (,) (,) [ (,) ]a r a x y F A
2 c o s ( )
0
| | (,) jrd P e d d
2 c o s ( )| | (,) jrP e d
c o s ( )
c o s ( )
( ) (,) |
(,) |
[ c o s( ),]
r
xrr
r r x r
r
h x p x
gx
gr
卷积反投影重建
上式的物理意义是:
投影 经过传递函数为的滤波器滤波后得到的修正后的投影在 处的取值。后者恰是通过给定点的射线方程。
| | [ ( ) ]rF h x
(,)rgx?
c o s ( )rxr
(,)rpx?
卷积反投影重建
可得
物理意义:
经过给定点 的所有滤波后的(射线)
投影在 0- 范围内的累加 ——反投影重建,得出 点的像素值。
0
(,) [ c o s ( ),]a r g r d
(,)r?
(,)r?
卷积反投影重建
称“滤波反投影方程”,它体现了滤波(卷积)
反投影算法的三个步骤:
1 把在固定视角 下测得的投影 经过滤波,得到滤波后投影 ;
2 对每一个,把 反投射于满足的射线上的所有各点 ;
3 将步骤 2中的反投影值对所有 进行累加(积分),得重建后的图像。
i?
i?
(,)ripx?
(,)rigx?
(,)rigx?
(,)r?
c o s ( )rixr
0
0
(,) [ c o s ( ),]a r g r d
卷积反投影重建平行射束卷积反投影重建算法 的 计算机实现
滤波函数与内插函数的选取原则
常用滤波函数举例
卷积反投影算法计算机实现步骤成像的特点
1.投影数据的高频 (空间域频率 )分量幅度很小;
2.投影数据的采集是天然离散的;
投影数据滤波函数
3.存在噪声。
r
r
r x nd
r x nd r
n
p ( x,) | { p ( nd,) }
{ h( nd) } =h( x ) | h ( nd ) ( x nd )
| | B 1 / ( 2 d )
重建后的图象密度为:
滤波后的投影为:
c o s ( )0? (,) (,) (,) * ( ) rr r x ra r a x y p x h x d
c os sin
0
(,) ( )
r
N N
r x x m y m
m n N
d p nd m h x nd
(,) ( )(,) N r
nN
r p n d m h x n dp x m
(,) ( ) * ( ) ( ) * ( )N r r r
nN
p n d m x n d h x p n d h x?
方式一:根据给定,应用内插,得到在 轴上连续取值的,再求出 。则此时的 为怎么求?
(,)rp x m?
rx
()rhx ()rhx
()h nd
(,)rp x m?
方式二:先求采样点 n’d处滤波投影值,再内插求任意 处的滤波投影值:
rx
'',) (,) ( )(
N
nN
d m p n d m h n n dpn
',) (,) * ( )( rrm p n d m xpx
',0,1,2,3,nm
( ) [ ( ) ( )] * ( ) { ( ) } * ( )
N
r r r r
nN
h x h n d x n d x h n d x?
怎么求?
由上述各式得:
可见两方案是等同的。
,) (,) * ( ) * ( )( rrm p n d m h n d xpx
(,) * ( )rp n d m h x
(,)rp x m?
两种内插函数
紧邻内插:
线形内插:
1
1 / 2
( ) 0,5 / 2
0 / 2
r
rr
r
xd
x x d
xd
2
1
()
0
r
r
r
r
x
xd
dx
xd
滤波函数的选择
为了醒目起见,可将式 (4.39)中的 m 作为常数,于是有:
及
同理,把 看成常数,由及得:
) ( ) * ( ) * ( )( rrp n d h n d xpx
( ) * ( )rp n d h x?
(,) (,) * ( )r r rg x p x h x
(,) (,)rrg x p x
) ( ) * ( )( r r rp x h xpx?
对上式求傅立叶变换:
若采样间隔 d足够小,则:
所以 可表示为:
其中 是带宽为 1/(2d)的窗函数,即:
若不考虑离散和噪声,矩形窗显然是最好的选择:
1( ) 0,
2P d
()H()HW
W?
10,2W d
( ) ( ) ( )rp x P HF
1
1,
2( ) ( )
1
0,
2
r
dWW
d
若考虑离散和噪声,则:
若:
则:
可得如下等式:
( ) ( ) ( ) ( )rrp x p n d h n d xF F F F
( ) ( ) ( ) ( )s r rp n d P p x x n d
FF
( ) ( )sh n d HF
( ) ( )rxF
( ) ( ) ( ) ( )r s sp x P HF
( ) ( ) ( ) ( ) ( )ssP H P H
为使上式成立,或尽量减少误差,须注意以下几点:
1.选择 时,宜注意 的频谱在接近处,不宜取,而应压低,以便补偿混迭效应。
()H? ()H? 1 /( 2 )d
()H
2.在 处,不宜取 的另一理由是:
投影数据中含有噪声,其频谱均匀,高频不衰减,若仍取,则在 处,从而很大,起放大作用,显然不利。
3.离散后的 频谱 是 的周期延拓,
导致 的频率范围中 不为零,而上式左边在这些频率范围都为零。为消除这些影响,需选取适当的 。
1 /( 2 )d ()H
()H 1 /( 2 )d ()H? ()sH?
()h nd ()sH? ()H?
1 /( 2 )d ()sH?
()
总结
选择 及内插函数 很重要。
可调整 范围内的误差;
可解决 区间 的问题,
多用线性内插函数。
( ) | | ( )HW ()
()W? | | 1 /( 2 )d
() | | 1 /( 2 )d ( ) 0sH
紧邻内插函数的频谱 线性内插函数的频谱常用滤波函数举例
1,R—L滤波函数
(i)系统函数
式中,B=1/2d,且:
()RLH
( ) ( ) 2RLH W r e c t B
1,
0,2
B
r e c t
BB
(ii)相应的冲激响应 ()
R L rhx?
2() rB j a x
R L r Bh x e d
2 2 22 si n ( 2 ) si n ( )rrB c x B B c x B
(iii)采样序列
采样间隔为 d,对应的最高不失真空间频率为 1/(2d)
的离散形式如下:
偶数奇数如将 的离散表示进行线性内插,则得另一连续的空域函数,是 的一次近似。
()RLh nd?
()R L rhx?
2
2 2 2
1
,
04
( ) 0,
1
,
RL
nd
h n d n
n
nd?
()RLh nd?
()RLh nd? ()R L rhx?
2.S—L滤波函数
选取 sinc函数作为窗函数。于是有:
相应的冲激响应 为:
2( ) s i n ( / 2 ) s i n / 22SL BH c B r e c t BB
()S L rhx?
22( ) s i n
2
r
B jx
S L r B
Bh x e d
B
2
2
1 4 si n( 4 )14
2()
2 1 ( 4 )
rr
r
Bx BxB
Bx
如令 y=4BXr,则得:
2
2
1 sin ( )14
2( ) ( )
21SL
yyB
hy
y
采样序列
对 进行均匀采样,得:
将这一离散形式进行线性内插,则得连续的空域函数,
()SLh nd?
()SLh nd?
2 2 2
2( ),0,1,2,
( 4 1 )SLh n d ndn
()S L rhx?
2( ) ( ) ( ) * ( )S L r S L r r r
n
h x h x x n d x?
对经线性内插后的函数 进行傅立叶变换得:
用 S—L滤波函数重建的图像其振荡响应减小,对含噪声数据的重建质量也较 R—L滤波函数情况为好,
但在低频段不及 R—L滤波函数的重建质量高。
()S L rhx?
21 sin ( / 2 )
( ) sin
2 ( / 2 )S L r
ddhx
dd
F
计算机实现滤波反投影的主要步骤
卷积计算
射束计算与内插
反投影重建
CT值的窗口处理国产 XDN-1型 CT装臵主要技术参数
扫描方式:平移 /旋转,每次旋转 1度,共转 180
度
每度平行扫描采样点数,256点
图象矩阵,256× 256
CT值范围,-1024~ 1023
显示灰度,64级
断层数:每次扫描得到两个断层面
滤波后的投影函数定义式 (I):
离散形式 (II):
注意:
(1)滤波函数 h理论上是无限长的,但实际应用中不可能无限计算下去
(2)长卷积可以使用 FFT软件快速求解,但实际中使用硬件实现 (II)式,
速度更快。
(3) 投影 p的序列长度为 256(0到 255)。但根据 (II)式,p需要向两侧延拓到长度 766(-255到 510)。在 n=-255~ -1间补 p’=[p(0)+p(1)]/2;在
n=256~ 510间补 p’’=[p(254)+p(255)]/2;补平均值相比补零能减小偶然误差。
' ' '(,) (,) * ( ) ( ) ( )
r m r m r r r r rp x p x h x p x x h x d x
(,) (,) * ( ) (,) ( )
N l
llN
p n m p n m h n p n l m h l
卷积计算
考虑到滤波函数是对称的,m固定不变,滤波后投影的定义式可以进一步改写为 (III)式:
为了方便计算机实现,将 (III)式改写为 (IV)式:
继续变换可以避免负地址:
2 5 5 2 5 5
2 5 5 2 5 5
(,) ( ) ( ) ( ) ( ) ( )
ll
p n m p n p n l h l p n l h l
511
1
( ) ( 2 5 6 ) ( 2 5 6 )
l
p n p n l h l
卷积计算
511
'
''
1
( ) ( ) ( )ss
l
p n p n l h l
卷积计算
卷积计算的硬件如下图
RAM1存放未经滤波的投影数据 p(n),RAM2存放滤波函数的离散值 h(n),RAM4存放卷积后投影序列。
对于 RAM1,平移扫描采样点数为 256,加上 p’和 p’’,两层总共
516个数值,512不够,所以采用 1K单元,每单元长为 16位 (因为 p
的离散值必为非负值 )。
对于 RAM2,滤波函数序列 511点,所以采用 512× 32位。
每计算一点的卷积后投影总共需要 511次乘法和 510次加法。乘法用 16× 16乘法器实现,加法用运算累加器实现。
射束计算与内插
由于 是空间中任意一点的像素坐标,故按 算得的 值并不正好为 d的整数倍,不能对应 RAM4中的某一地址,
故需进行内插。
,c os s inr m i m j mx x y
(,)ijxy
,rmx
当 m确定 (也即 Φ确定 )后,其 xr可能在 n0d和 (n0+1)d之间,
取定 d=1,取定 m,有以下线性内插公式,
射束计算的任务就是求出 n0和 δ值:从而得到如下的公式:
0 0 0( ) ( 1 ) ( ) ( 1 )p n p n p n
01
11
c o s s i n ( ) c o s ( ) s i n
22
1
( 1 ) c o s ( 1 ) s i n ( c o s s i n )
2
r
NN
x x y i j
N
ij
射束计算与内插
上式表示相应于 的 介于第 n0号射束与 (n0+1)
号射束之间,与 n0号射束相距 δ。
在实际的计算中,先求出 (i,j)=(1,1)的值,以后当 i或 j增加时,只需要再迭加 cosΦ或 sin Φ即可。由于 Φ的离散值只有 180个,所以可以设臵 RAM3存储全部的 180组
cos Φ,sin Φ和 C值,通过以下框图实现内插:
0
1
2
1
( 1 ) c os ( 1 ) si n ( 1 c os si n )
2
( 1 ) c os ( 1 ) si n
rr
N
xx
N
ij
i j C n
rx
射束计算与内插
(,)ijxy
射束计算与内插实际应用的离散形式:
180
,
1
(,) [ (,),]rm
m
f i j p x i j m
(,) (,) * ( )r r rp x p x h x
0
(,) [ c o s ( ),]f r p r d
反投影重建反投影重建
累加求和部分由累加器和暂存图象的 WRAM完成。其中
WRAM为 64K× 24位。 64K是由于图象总共有
256× 256=64K个象素点; 24位则表明了重建图象的颜色质量。
总结整个卷积反投影重建计算机实现的流程:
(1)在当前 Φ下,扫描求出投影 p,对投影 p进行卷积滤波,结果存在 RAM4。
(2)对滤波后投影进行与图象坐标地址匹配的内插。
(3)内插结果累加到存储图象的 WRAM中。
(4) Φ步进,重复前三项步骤,直到 Φ=180。
(5)当 Φ=180,WRAM的图象累加完成,整个图象的全部象素的重建结束。
因为全部象素同时完成重建,所以全过程属于投影驱动法。
CT值的窗口处理
CT装臵的主要优点是引进了 CT数,并可对 CT
值开窗以选择一段动态范围,对应于监视器上的全部灰度,从而大大提高密度鉴别力。
常用的窗口处理包括:
单窗
SIGMA窗
双窗
自适应窗单窗
CT值的窗口处理
CT值的窗口处理
双窗:可以同时观察两种 CT值相差很大的不同组织的情况。双窗之间要在断层的平面区域之间划出界限,否则会导致混淆。
SIGMA窗,CT值到灰度值的映射曲线类似反转的 ∑。可以认为是双窗的一种变种。需要由有经验的医生诊断鉴别。
自适应窗:单窗的分段实现,可以同时观察两种 CT值差别很大的不同组织,而且不会出现区域的混淆。采用的方法是在单窗的基础上弃除没有诊断价值的一段 CT值范围。
扇束投影重建算法及其实现
笔束平移 /旋转扫描扫描时间长,只宜用于头颅扫描,
不能用来观察肺等动态脏器。
速度慢的根源:
1、采用单 X射线源单检测器;
2、既有平移运动又有旋转运动,运动方式改变频繁。
提高数据采集速度的途径:
1、扇束扫描:单 X射线源多检测器
2、机架连续旋转扇形扇形扇形各代 CT比较第一代 第二代 第三代 第四代 第五代扫描方式 笔束扫描 扇束扫描运动方式 平移 /旋转方式 连续扫描方式扫描时间 3min 10s~2min 2.8s~10s 1s~10s 更快主要用途 头颅扫描 全身扫描,观察除心脏外的脏器可用于血管造影和心脏造影第一代,笔束平移旋转第二代,扇束平移 /旋转第三代,扇束,连续旋转,单 X射线源多检测器同步旋转第四代,扇束,连续旋转方式第五代,C-100超高速 CT
扇束投影重建算法的分类
重排算法把一个视图中采得的扇形数据重新组合成平行的射线投影数据,然后用上一节的算法重建。
直接重建算法不必数据重排,只需适当加权即可运用与上一节类似的算法重建扇形射线的型式
等角射线型? 等距射线型等角射线扇束数据的直接重建算法
12DD
0S
,X射线源
:检测器阵所在弧线
:中心射线00SD
扇形位臵由该中心射线与 y轴交角? 确定,同一扇形中的射线 由?
规定,射线的绝对位臵由(?,?)唯一确定。
0SE
等角射线扇束数据的直接重建算法我们的任务是:给定重建物体断面图象 a(r,?)。
故射线 的位臵所对应的射线投影为:
现在要用到 旋转角的扫描,由扫描装臵的对称性,
可知:
(,)fp
(,) ( sin,) (,)frp x p D p
(,) (,)rrp x p x
(,) ( s i n,)rxD
,0 2,mm
0SE
360。
s i nrxD;
重建算法的推导
1.改写反投影重建表达式
2.改换变量以适应扇束投影数据
2
c o s ( )0
1(,) (,) * ( ) |
2 rr r x ra r p x h x d
0
2
0 c o s ( )0 a r c s in
1(,) ( si n '','') [ si n ( '') ] | | | ''
2
m
m r L
a r p D h L J d d
0
2
0 c o s ( )0 a r c s in
1 ( '',) [ si n ( '') ] | c o s '' ''
2
m
m
f
r
L
p h L D d d
重建算法的推导重建算法的推导
3,重臵滤波函数变量
4,用扇束投影数据与扇形变量表示重建公式
其中
2
0
00 22
0
( '')[ s in ( '') ] ( '')
s in ( '')
h L h
L
0
22
220
1(,) [ (,) c o s ] * ( ) |
2 s in
fa r p D h d
L
22 2 si n( )L D r D r
0
c o s ( )a r c s i n r
L
22
c o s ( )a r c s i n
2 s i n ( )
r
D r D r
扇束算法平行束算法结论
扇束情况下的重建算法较为复杂,但实质没有改变。可采用平行束情况下的算法实现,只需加以适当地修正即可。
扇束扫描程序
Z
Y
螺旋锥形束 CT
(From G,Wang)
心脏成像动态成像
4D CT
4D CT
4D CT 检测器
Multi-source cone-beam CT
Detector
Source
(From G,Wang)
小动物成像 Micro-CT
Completeness Conditions for Cone-
Beam CT
Completeness Condition,If on every plane
that intersects the object there exists at least
one cone-beam source point,then one can
reconstruct the object,(B.D.Smith1985)
The possible structures satisfying the conditions are:
Works on Cone-Beam CT
Kirillov A.A,(1961)—Completeness condition
Tuy H.K,(1983)-- Completeness condition
Feldkamp L.A.(1984)---Practical,Incomplete projection
Smith B.D.(1985)-- Completeness condition
Grangeat P.(1987)—Reconstruction Algorithm
Ge Wang,(1991)—Reconstruction Algorithm
Defrise M.Clark R.(1994)-- Reconstruction Algorithm
Kudo H,Saito T.(1991,1994)– Reconstruction Algorithm
CT重建算法
近似重建算法
迭代重建算法
精确重建算法
Tam-Danielsson Geometry
Tam,Three-dimensional computerized tomography
scanning method and system for large objects with
smaller area detector,US Patent 5,390,112,1995
Detector surface is limited to
a cylindrical section between
two consecutive helical turns
Every point is on one and
only one PI-line
Pi-Line
Perfect Mosaic
Tam,Samarasekera,Sauer,Exact cone-beam CT with a spiral scan,Phys,Med,Biol,43:847-855,1998
Kudo,Noo,Defrise,Cone-beam filtered-backprojection algorithm for truncated helical data,
Phys,Med,Biol,43:2885-2909,1998
(From G,Wang)
Katsevich Theorem (2002)
dsdxsqyDqsyxxf sqf
xI PI
s i n
1|)),,(),((
|)(|
1
2
1)( 2
0)(
),(),(),( xsuxsxse ),(s in),(c o s),,( xsexsxs
)(0sy
e
),( xsu
),(0xs?
)(1sy
)(2sy
)(xf
Detector
Plate
Source
Object
Pi-Line
(From G,Wang)
Noo’s Short-Scan (2002)
Integral on any line through a region
(From G,Wang)
我国 CT的发展
1) 85年第二代 CT,上海计算技术研究所等
2) 安科公司
3) Siemens
4) GE
5)中国人有才智、有能力造好 CT。
投影投影投影投影投影投影投影投影
月食、日食
日晷(计时)
几何投影,坐标分量(位臵)
魔术
电影
猴子捞月
杯弓蛇影
顾影自怜
希腊美少年?Narcissus( 水仙 )? narcissi( 自恋 )
投影投影投影经典断层成像投影谁先匹配准,奖励 20万通缉令投影投影反投影重建算法的一般步骤:
原像取投影 反投影重建重建后图像反投影重建算法投影重建算法的基本内容:
“断层平面中某一点的密度值可看作这一平面内所有经过该点的射线投影之和(的平均值)”。
式中 xk表示象素的值,pk,i为经过象素的第 i条射线投影。
可以这样理解:
其中 f(x,y,z)是身体组织密度。
,
1
1 pn
k k i
p i
xp
n?
(,,)lp f x y z d l
反投影重建算法的物理概念算法举例
1
1 3 7 1 1 1 5 2p x x x x
2
2 2 6 1 0 1 4 6p x x x x3
3 5 6 7 8 7p x x x x4
4 9 1 0 1 1 1 2 1p x x x x
5
5 1 6 1 1 1 6 5p x x x x
6
6 4 7 1 0 1 3 3p x x x x
算法举例
根据反投影算法
x1=p5 = 5
x6=p2+p3+p5=18
…
平均化处理,除以投影线数目
xi=xi/6
0 0 0 0
0 5 2 0
0 1 0 0
0 0 0 0
5 6 2 3
7 18 12 7
1 10 8 1
3 6 2 5
0.83 1 0.33 00.5
1.16 3 2 1.16
0.06 1.66 1.33 0.16
0.5 1 0.33 0.83
反投影重建后原像素值再除以投影线数,平均化断层平面中某一点的密度值可看作这一平面内所有经过该点的射线投影之和的平均值
12
3
4
5
6
星状伪迹
反投影重建后,原来为 0的点不再为 0,形成伪迹
0 0 0 0
0 5 2 0
0 1 0 0
0 0 0 0
0.83 1 0.33 0.5
1.16 3 2 1.16
0.06 1.66 1.33 0.16
0.5 1 0.33 0.83
原像素值 再除以投影线数,平均化星状伪迹
我们考虑孤立点源反投影重建,
中心点 A经 n条投影线投影后,投影值均为 1:
p1=p2=...=pn=1
因此重建后而其他点均为 1/n
这类伪迹成为星状伪迹
1/n 1/n 1/n
1/n 1 1/n
1/n 1/n 1/n
12
1 (,,,) 1
Anf p p pn
0 0 0
0 1 0
0 0 0
星状伪迹
产生星状伪迹的原因在于:反投影重建的本质是把取自有限物体空间的射线投影均匀地回抹
(反投影 )到射线所及的无限空间的各点之上,
包括原先像素值为零的点反投影重建算法的数学描述我们设臵一旋转坐标系统 Xr-Yr,它绕原点转动使投影线总是沿着 Yr方向。 Xr-Yr的原点与 X-Y的原点相重。两者的夹角为 φ,不同的 φ代表不同的投射方向。投影线的位臵可由
(Xr,φ)完全确定。 空间的任一点的位臵可用 (X,Y),(Xr,Yr)或极坐标 (r,θ)表出。
反投影重建算法的数学描述
先设 φ为离散取值,则投影为
根据反投影重建的定义式,点 的图像在所述坐标系统表示为式中,,为投影数
1 1( ) ( ) ( ) | c o sr r r rx y d y x r
1
1 1( ) (,) (,) (,)r r r r r r r r rp x p x f x y d y x y d y
,rrxy
1
11
1
(,),
11
c os c os
i
ii
N
r r r r
i
NN
ii
ii
f r f x y p x
N
p r p r
N
/ N N?
若在有限区间内射线增至不相重的无限条,即连续投影,
则式 (4.12)过渡到更一般的连续情况下的反投影表达式:
在输入图像为点源的情况下,由式 (4.10)及式 (4.13)可得
0
1,c o sf r p r d?
0
00
0
22
0
11
,c o s
s in
1 1 1 1 1
s in
h r r d d
r
rr xy
反投影重建算法的数学描述
可见,反投影重建算法相应的系统的扩展函数不是 δ函数,系统不是完美的.式定量地描述了反投影重建算法星状伟迹的本质.
要除去反投影算法的星状伪迹,我们可以在输出端加一滤波器,使加了滤波器后的反投影重建成像系统 PSF
= δ(x,y).使滤波器的 PSF为 q(x,y),相应的传递函数为
Q(ξ,η),这里
我们要求:
2,,Q F q x y
22
1 * *,(,)q x y x y
xy
反投影重建算法的数学描述
**表示二维卷积,对上式取二维傅立叶变换得
或
这是一只二维滤波器,实现起来比较麻烦.若
ρ的变换范围可扩至 ∞,根本不能实现,但不管怎样,它提供了去除星状伪迹的一个努力方向.
221,1Q
22,Q
反投影重建算法的数学描述
反投影重建算法要求对任一点,找出经过该点的所有射线,从而求得过该点的射线投影的均值。
从图中可以看出,过给定点的射线 都在以 r为直径,以 为圆心的圆周上,即正弦图
(,)r?
(,)rx?
c o s ( )rxr
1(,)
2r?
正弦图正弦图
在 空间中,在给定 情况下,
我们得到的是一条正弦曲线。图像中的每一点,对应于空间中的一条正弦曲线,
全部图像对应于一簇正弦曲线族,这就是正弦图。
(,)rx? (,)r?
正弦图反投影重建算法的实现
正弦图中,射线位臵用 描述,扫描操作由角增量
Δ,步距 d步进;若正弦图的 Δ,d不是足够小,图中的某些点 可能找不到一点( mΔ,nd)与其对应;
若令 Фi= m Δ,那么,
一般地,故不能直接得到,须内插。
1
0
(,) 1 / [ c o s( ) ]mM
m
f r N p m r m
c o s ( )r m n d [ c o s ( ) ]mp r m
(,)rx?
(,)r?
内插内插方法
1.紧邻内插,在 时,若则以 中的内容 直接代替
。 否则,用 代替
2.线形内插:如果,则
[ c o s ( ) ]mp r m
( ) { [ ( 1 ) ] ( ) } * [ c o s ( ) ] /m m mp n d p n d p n d r m n d d
m m
c o s ( ) c o s ( ) ( 1 )r m n d r m n d
(,)m nd?
[ c o s ( ) ]mp r m
()mp nd?
[ ( 1 ) ]mp n d
[ c o s ( ) ]mp r m
c os( ) ( 1 )nd r m n d
滤波反投影算法
去除星状伪迹的两种方法:
取投影取投影 一维滤波器 反投影重建二维滤波器反投影重建输入图像(原图像)
输入图像
(原图像)
输出图像
(同原图像 )
输出图像
(同原图像 )
图 4.18 星状尾迹的去除
( b)
( a)
第一种方法中的二维滤波器较难实现。
第二种方法交换了滤波与反投影的顺序,
此时投影数据是一维的,易于实现。此方法的理论基础就是我们即将提到的投影定理(中心切片定理)
滤波反投影算法投影定理
某一图像 在视角 时投影 的 1D傅立叶变换给出 的 2D傅立叶变换的一个切片,切片与 轴相交成 角,且通过坐标原点。即
投影定理的 证明,p138- 139
1?[ ( ) ] (,) |rF p x F 固 定
12?(,) (,)FF(,)f x y
(,)f x y
()rpx
1?
投影定理卷积反投影重建
待建图像为 a(x,y),它的 2D傅氏变换为,
由中心切片定理:
12?(,) (,)AA
12?(,) (,) [ ( ) ]rA A F p x
( ) (,)PP
卷积反投影重建
待建图像
其中
12 1 2? (,) (,) [ (,) ]a r a x y F A
2 c o s ( )
0
| | (,) jrd P e d d
2 c o s ( )| | (,) jrP e d
c o s ( )
c o s ( )
( ) (,) |
(,) |
[ c o s( ),]
r
xrr
r r x r
r
h x p x
gx
gr
卷积反投影重建
上式的物理意义是:
投影 经过传递函数为的滤波器滤波后得到的修正后的投影在 处的取值。后者恰是通过给定点的射线方程。
| | [ ( ) ]rF h x
(,)rgx?
c o s ( )rxr
(,)rpx?
卷积反投影重建
可得
物理意义:
经过给定点 的所有滤波后的(射线)
投影在 0- 范围内的累加 ——反投影重建,得出 点的像素值。
0
(,) [ c o s ( ),]a r g r d
(,)r?
(,)r?
卷积反投影重建
称“滤波反投影方程”,它体现了滤波(卷积)
反投影算法的三个步骤:
1 把在固定视角 下测得的投影 经过滤波,得到滤波后投影 ;
2 对每一个,把 反投射于满足的射线上的所有各点 ;
3 将步骤 2中的反投影值对所有 进行累加(积分),得重建后的图像。
i?
i?
(,)ripx?
(,)rigx?
(,)rigx?
(,)r?
c o s ( )rixr
0
0
(,) [ c o s ( ),]a r g r d
卷积反投影重建平行射束卷积反投影重建算法 的 计算机实现
滤波函数与内插函数的选取原则
常用滤波函数举例
卷积反投影算法计算机实现步骤成像的特点
1.投影数据的高频 (空间域频率 )分量幅度很小;
2.投影数据的采集是天然离散的;
投影数据滤波函数
3.存在噪声。
r
r
r x nd
r x nd r
n
p ( x,) | { p ( nd,) }
{ h( nd) } =h( x ) | h ( nd ) ( x nd )
| | B 1 / ( 2 d )
重建后的图象密度为:
滤波后的投影为:
c o s ( )0? (,) (,) (,) * ( ) rr r x ra r a x y p x h x d
c os sin
0
(,) ( )
r
N N
r x x m y m
m n N
d p nd m h x nd
(,) ( )(,) N r
nN
r p n d m h x n dp x m
(,) ( ) * ( ) ( ) * ( )N r r r
nN
p n d m x n d h x p n d h x?
方式一:根据给定,应用内插,得到在 轴上连续取值的,再求出 。则此时的 为怎么求?
(,)rp x m?
rx
()rhx ()rhx
()h nd
(,)rp x m?
方式二:先求采样点 n’d处滤波投影值,再内插求任意 处的滤波投影值:
rx
'',) (,) ( )(
N
nN
d m p n d m h n n dpn
',) (,) * ( )( rrm p n d m xpx
',0,1,2,3,nm
( ) [ ( ) ( )] * ( ) { ( ) } * ( )
N
r r r r
nN
h x h n d x n d x h n d x?
怎么求?
由上述各式得:
可见两方案是等同的。
,) (,) * ( ) * ( )( rrm p n d m h n d xpx
(,) * ( )rp n d m h x
(,)rp x m?
两种内插函数
紧邻内插:
线形内插:
1
1 / 2
( ) 0,5 / 2
0 / 2
r
rr
r
xd
x x d
xd
2
1
()
0
r
r
r
r
x
xd
dx
xd
滤波函数的选择
为了醒目起见,可将式 (4.39)中的 m 作为常数,于是有:
及
同理,把 看成常数,由及得:
) ( ) * ( ) * ( )( rrp n d h n d xpx
( ) * ( )rp n d h x?
(,) (,) * ( )r r rg x p x h x
(,) (,)rrg x p x
) ( ) * ( )( r r rp x h xpx?
对上式求傅立叶变换:
若采样间隔 d足够小,则:
所以 可表示为:
其中 是带宽为 1/(2d)的窗函数,即:
若不考虑离散和噪声,矩形窗显然是最好的选择:
1( ) 0,
2P d
()H()HW
W?
10,2W d
( ) ( ) ( )rp x P HF
1
1,
2( ) ( )
1
0,
2
r
dWW
d
若考虑离散和噪声,则:
若:
则:
可得如下等式:
( ) ( ) ( ) ( )rrp x p n d h n d xF F F F
( ) ( ) ( ) ( )s r rp n d P p x x n d
FF
( ) ( )sh n d HF
( ) ( )rxF
( ) ( ) ( ) ( )r s sp x P HF
( ) ( ) ( ) ( ) ( )ssP H P H
为使上式成立,或尽量减少误差,须注意以下几点:
1.选择 时,宜注意 的频谱在接近处,不宜取,而应压低,以便补偿混迭效应。
()H? ()H? 1 /( 2 )d
()H
2.在 处,不宜取 的另一理由是:
投影数据中含有噪声,其频谱均匀,高频不衰减,若仍取,则在 处,从而很大,起放大作用,显然不利。
3.离散后的 频谱 是 的周期延拓,
导致 的频率范围中 不为零,而上式左边在这些频率范围都为零。为消除这些影响,需选取适当的 。
1 /( 2 )d ()H
()H 1 /( 2 )d ()H? ()sH?
()h nd ()sH? ()H?
1 /( 2 )d ()sH?
()
总结
选择 及内插函数 很重要。
可调整 范围内的误差;
可解决 区间 的问题,
多用线性内插函数。
( ) | | ( )HW ()
()W? | | 1 /( 2 )d
() | | 1 /( 2 )d ( ) 0sH
紧邻内插函数的频谱 线性内插函数的频谱常用滤波函数举例
1,R—L滤波函数
(i)系统函数
式中,B=1/2d,且:
()RLH
( ) ( ) 2RLH W r e c t B
1,
0,2
B
r e c t
BB
(ii)相应的冲激响应 ()
R L rhx?
2() rB j a x
R L r Bh x e d
2 2 22 si n ( 2 ) si n ( )rrB c x B B c x B
(iii)采样序列
采样间隔为 d,对应的最高不失真空间频率为 1/(2d)
的离散形式如下:
偶数奇数如将 的离散表示进行线性内插,则得另一连续的空域函数,是 的一次近似。
()RLh nd?
()R L rhx?
2
2 2 2
1
,
04
( ) 0,
1
,
RL
nd
h n d n
n
nd?
()RLh nd?
()RLh nd? ()R L rhx?
2.S—L滤波函数
选取 sinc函数作为窗函数。于是有:
相应的冲激响应 为:
2( ) s i n ( / 2 ) s i n / 22SL BH c B r e c t BB
()S L rhx?
22( ) s i n
2
r
B jx
S L r B
Bh x e d
B
2
2
1 4 si n( 4 )14
2()
2 1 ( 4 )
rr
r
Bx BxB
Bx
如令 y=4BXr,则得:
2
2
1 sin ( )14
2( ) ( )
21SL
yyB
hy
y
采样序列
对 进行均匀采样,得:
将这一离散形式进行线性内插,则得连续的空域函数,
()SLh nd?
()SLh nd?
2 2 2
2( ),0,1,2,
( 4 1 )SLh n d ndn
()S L rhx?
2( ) ( ) ( ) * ( )S L r S L r r r
n
h x h x x n d x?
对经线性内插后的函数 进行傅立叶变换得:
用 S—L滤波函数重建的图像其振荡响应减小,对含噪声数据的重建质量也较 R—L滤波函数情况为好,
但在低频段不及 R—L滤波函数的重建质量高。
()S L rhx?
21 sin ( / 2 )
( ) sin
2 ( / 2 )S L r
ddhx
dd
F
计算机实现滤波反投影的主要步骤
卷积计算
射束计算与内插
反投影重建
CT值的窗口处理国产 XDN-1型 CT装臵主要技术参数
扫描方式:平移 /旋转,每次旋转 1度,共转 180
度
每度平行扫描采样点数,256点
图象矩阵,256× 256
CT值范围,-1024~ 1023
显示灰度,64级
断层数:每次扫描得到两个断层面
滤波后的投影函数定义式 (I):
离散形式 (II):
注意:
(1)滤波函数 h理论上是无限长的,但实际应用中不可能无限计算下去
(2)长卷积可以使用 FFT软件快速求解,但实际中使用硬件实现 (II)式,
速度更快。
(3) 投影 p的序列长度为 256(0到 255)。但根据 (II)式,p需要向两侧延拓到长度 766(-255到 510)。在 n=-255~ -1间补 p’=[p(0)+p(1)]/2;在
n=256~ 510间补 p’’=[p(254)+p(255)]/2;补平均值相比补零能减小偶然误差。
' ' '(,) (,) * ( ) ( ) ( )
r m r m r r r r rp x p x h x p x x h x d x
(,) (,) * ( ) (,) ( )
N l
llN
p n m p n m h n p n l m h l
卷积计算
考虑到滤波函数是对称的,m固定不变,滤波后投影的定义式可以进一步改写为 (III)式:
为了方便计算机实现,将 (III)式改写为 (IV)式:
继续变换可以避免负地址:
2 5 5 2 5 5
2 5 5 2 5 5
(,) ( ) ( ) ( ) ( ) ( )
ll
p n m p n p n l h l p n l h l
511
1
( ) ( 2 5 6 ) ( 2 5 6 )
l
p n p n l h l
卷积计算
511
'
''
1
( ) ( ) ( )ss
l
p n p n l h l
卷积计算
卷积计算的硬件如下图
RAM1存放未经滤波的投影数据 p(n),RAM2存放滤波函数的离散值 h(n),RAM4存放卷积后投影序列。
对于 RAM1,平移扫描采样点数为 256,加上 p’和 p’’,两层总共
516个数值,512不够,所以采用 1K单元,每单元长为 16位 (因为 p
的离散值必为非负值 )。
对于 RAM2,滤波函数序列 511点,所以采用 512× 32位。
每计算一点的卷积后投影总共需要 511次乘法和 510次加法。乘法用 16× 16乘法器实现,加法用运算累加器实现。
射束计算与内插
由于 是空间中任意一点的像素坐标,故按 算得的 值并不正好为 d的整数倍,不能对应 RAM4中的某一地址,
故需进行内插。
,c os s inr m i m j mx x y
(,)ijxy
,rmx
当 m确定 (也即 Φ确定 )后,其 xr可能在 n0d和 (n0+1)d之间,
取定 d=1,取定 m,有以下线性内插公式,
射束计算的任务就是求出 n0和 δ值:从而得到如下的公式:
0 0 0( ) ( 1 ) ( ) ( 1 )p n p n p n
01
11
c o s s i n ( ) c o s ( ) s i n
22
1
( 1 ) c o s ( 1 ) s i n ( c o s s i n )
2
r
NN
x x y i j
N
ij
射束计算与内插
上式表示相应于 的 介于第 n0号射束与 (n0+1)
号射束之间,与 n0号射束相距 δ。
在实际的计算中,先求出 (i,j)=(1,1)的值,以后当 i或 j增加时,只需要再迭加 cosΦ或 sin Φ即可。由于 Φ的离散值只有 180个,所以可以设臵 RAM3存储全部的 180组
cos Φ,sin Φ和 C值,通过以下框图实现内插:
0
1
2
1
( 1 ) c os ( 1 ) si n ( 1 c os si n )
2
( 1 ) c os ( 1 ) si n
rr
N
xx
N
ij
i j C n
rx
射束计算与内插
(,)ijxy
射束计算与内插实际应用的离散形式:
180
,
1
(,) [ (,),]rm
m
f i j p x i j m
(,) (,) * ( )r r rp x p x h x
0
(,) [ c o s ( ),]f r p r d
反投影重建反投影重建
累加求和部分由累加器和暂存图象的 WRAM完成。其中
WRAM为 64K× 24位。 64K是由于图象总共有
256× 256=64K个象素点; 24位则表明了重建图象的颜色质量。
总结整个卷积反投影重建计算机实现的流程:
(1)在当前 Φ下,扫描求出投影 p,对投影 p进行卷积滤波,结果存在 RAM4。
(2)对滤波后投影进行与图象坐标地址匹配的内插。
(3)内插结果累加到存储图象的 WRAM中。
(4) Φ步进,重复前三项步骤,直到 Φ=180。
(5)当 Φ=180,WRAM的图象累加完成,整个图象的全部象素的重建结束。
因为全部象素同时完成重建,所以全过程属于投影驱动法。
CT值的窗口处理
CT装臵的主要优点是引进了 CT数,并可对 CT
值开窗以选择一段动态范围,对应于监视器上的全部灰度,从而大大提高密度鉴别力。
常用的窗口处理包括:
单窗
SIGMA窗
双窗
自适应窗单窗
CT值的窗口处理
CT值的窗口处理
双窗:可以同时观察两种 CT值相差很大的不同组织的情况。双窗之间要在断层的平面区域之间划出界限,否则会导致混淆。
SIGMA窗,CT值到灰度值的映射曲线类似反转的 ∑。可以认为是双窗的一种变种。需要由有经验的医生诊断鉴别。
自适应窗:单窗的分段实现,可以同时观察两种 CT值差别很大的不同组织,而且不会出现区域的混淆。采用的方法是在单窗的基础上弃除没有诊断价值的一段 CT值范围。
扇束投影重建算法及其实现
笔束平移 /旋转扫描扫描时间长,只宜用于头颅扫描,
不能用来观察肺等动态脏器。
速度慢的根源:
1、采用单 X射线源单检测器;
2、既有平移运动又有旋转运动,运动方式改变频繁。
提高数据采集速度的途径:
1、扇束扫描:单 X射线源多检测器
2、机架连续旋转扇形扇形扇形各代 CT比较第一代 第二代 第三代 第四代 第五代扫描方式 笔束扫描 扇束扫描运动方式 平移 /旋转方式 连续扫描方式扫描时间 3min 10s~2min 2.8s~10s 1s~10s 更快主要用途 头颅扫描 全身扫描,观察除心脏外的脏器可用于血管造影和心脏造影第一代,笔束平移旋转第二代,扇束平移 /旋转第三代,扇束,连续旋转,单 X射线源多检测器同步旋转第四代,扇束,连续旋转方式第五代,C-100超高速 CT
扇束投影重建算法的分类
重排算法把一个视图中采得的扇形数据重新组合成平行的射线投影数据,然后用上一节的算法重建。
直接重建算法不必数据重排,只需适当加权即可运用与上一节类似的算法重建扇形射线的型式
等角射线型? 等距射线型等角射线扇束数据的直接重建算法
12DD
0S
,X射线源
:检测器阵所在弧线
:中心射线00SD
扇形位臵由该中心射线与 y轴交角? 确定,同一扇形中的射线 由?
规定,射线的绝对位臵由(?,?)唯一确定。
0SE
等角射线扇束数据的直接重建算法我们的任务是:给定重建物体断面图象 a(r,?)。
故射线 的位臵所对应的射线投影为:
现在要用到 旋转角的扫描,由扫描装臵的对称性,
可知:
(,)fp
(,) ( sin,) (,)frp x p D p
(,) (,)rrp x p x
(,) ( s i n,)rxD
,0 2,mm
0SE
360。
s i nrxD;
重建算法的推导
1.改写反投影重建表达式
2.改换变量以适应扇束投影数据
2
c o s ( )0
1(,) (,) * ( ) |
2 rr r x ra r p x h x d
0
2
0 c o s ( )0 a r c s in
1(,) ( si n '','') [ si n ( '') ] | | | ''
2
m
m r L
a r p D h L J d d
0
2
0 c o s ( )0 a r c s in
1 ( '',) [ si n ( '') ] | c o s '' ''
2
m
m
f
r
L
p h L D d d
重建算法的推导重建算法的推导
3,重臵滤波函数变量
4,用扇束投影数据与扇形变量表示重建公式
其中
2
0
00 22
0
( '')[ s in ( '') ] ( '')
s in ( '')
h L h
L
0
22
220
1(,) [ (,) c o s ] * ( ) |
2 s in
fa r p D h d
L
22 2 si n( )L D r D r
0
c o s ( )a r c s i n r
L
22
c o s ( )a r c s i n
2 s i n ( )
r
D r D r
扇束算法平行束算法结论
扇束情况下的重建算法较为复杂,但实质没有改变。可采用平行束情况下的算法实现,只需加以适当地修正即可。
扇束扫描程序