轩敞-黄毛丫头
一、 平行束反投影重建算法
平行束
重建采用的是平移加
平移
旋转的扫描方式,如图所示,射线源
在某一角度下
水平移动,将物体全部
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
(
?
)
频谱幅度增加。如图所示,为考虑这一情况,应注意
H(
?
)
的 频谱在接近
?
?12d
处,不
宜取
H(
?
)??
,而应压低,以便补偿混叠效应。
P(
?
)
1
?< br>2d
1
2d
?
(a)
P(
?
)
P
s
(
?
)
1
d
?
1
2d
1
2d
1
d
? (b)
P
s
(
?
)
图 阐明投影定理的频谱混叠
(2) 在
?
?12d
处不宜取
H(
?
)?
?
的另一个原因是投影数据中含有噪声,其频谱均
匀,若仍取
H(
?
)?
?
,则在
?
?12d
处,噪声被放大。
(3) 离 散后的
{h(nd)}
频谱
H
s
(
?
)
是
H(
?
)
的周期延拓,导致
?
?12d
的频率范围 中
H
s
(
?
)
不为零,如图所示,而式右边在这些频率范围 都为零。为消除这些影响,需要选取适
当的
?(
?
)
。
H (
?
)
?
1
2d
0
1
2d
?
(a)
H
s
(
?
)
1
d
?
1
2d<
br>0
1
2d
1
d
?
(b)
图
频域滤波函数及其周期延拓
常用的内插函数有线性内插和紧邻内插,其频谱分别为
?
1
(
?
)?dsinc
2
(
?
d)
和?
2
(
?
)?dsinc(
?
d)
如图所示,
其中
sinc(x)?sin(
?
x)
?
x
,目前几乎都采
用线性内插函数,
因为它在第一个过零点后迅速趋于零,这一点上明显好于紧邻内插。
?1
(
?
)
1.0
1.0
?
2
(
?
)
0.5
0.5
0
1
d
?
0
1
d
?
(a)线性内插
(b)紧邻内插
图 内插函数频谱图
常用的滤波函数
1.
R?L
滤波函数
(1)
系统函数
H
R?L
?
?
?
H
R?L?
?
?
?
?
W(
?
)
?
1<
br>?
?
W(
?
)?
?
?
0
?
?
?
?
1
2d
1
?
?
2d
(2)
相应的冲激响应
h
R?L
(x
r
)
h
R
?L
(x
r
)?
?
B
?B
?
e
j
2
??
x
d
?
r
?
sin
?
x<
br>r
B
?
sin2
?
x
r
B
?2B<
br>2
?B
2
??
2
?<
br>x
r
B
?
xB
r
??
?2B
2sinc(2x
r
B)?B
2
sinc
2
(x
r
B)
2
(3) 采样序列
h
R?L
(nd)
<
br>这里的采样间隔同前一样为
d
,以
x
r
?nd
代入式
中,得到
h
R?L
(x
r
)
的离散形式如下:
?
1
?
4d
2
?
h
R-L
( nd)?
?
0
?
?1
?
222
?
n
?
d
n?0
n?偶数
n?奇数
如图所示
1
4d
2
-5d
-4d
-3d
-2d
?dd
3d
2d4d
5d
x
r
?nd
0
图
R?L
滤波函数离散表示
2.
S?L
滤波函数
为了更好 地补偿
P
s
(
?
)
在
?
?1(2d)处的混叠,设法使
?
?1(2d)
处的
H(
?
)
幅值压得
更低,可以通过选取适宜的窗函数
W(
?
)
达到。例如, 可取
sinc
函数作为窗函数。
(1) 系统函数
H
S?L
?
?
?
H
S?L?
?
?
?
?
sinc(
?
2B)W(
?
)?
2B
?
sin
??
2B
W(
?)
(2) 相应的冲激函数
h
S?L
(x
r
)
h
S ?L
(x
r
)?
?
B
2B
?B
?
sin
??
2B
e
j2
??
x
r
d
?
?
?
?
1?4Bxsi n?4Bx
2
r
?
r
?
1
?
4B
?
2
??
?
??
2
2
?
?
?1?
?
4Bx
r
?
(3) 采样序列
h
S?L
(nd)
这里的采样间隔同样为
d,以
x
r
?nd
代入式中,得到
h
R?L
(x
r
)
的离散形式如下:
h
S?L
(nd)?
如图所示
?2
?
2
d
2
(4n
2
?1)
,n?0,?1,?2,
h
S-L
(nd)
2
?
d
22-4d
-3d
-2d?d
d
2d3d4d
x
r
?nd
图
S?L
滤波函数离散表示
四、 卷积反投影算法的计算机实现
卷积计算
数据采集在空间是离散的,因为
x
r
?nd
,故对应于<
br>x
r
取变量为
n
,我们假设采样点为
256点,即
n
?0255
。从下面的推导可以看到
p
?
n
?
只在
n?0255
间取值是不够的,要用
到
n??255?1
及
n?25
6510
间的
p(n)
值。
N
t
p(n)?p(n)*h
(n)?
l??N
t
?
p(n?l)h(l)
p(0)?h(?255)p(?255)?......h(0)p(0)?......h(255
)p(255)
p(1)?h(?255)p(?254)?......h(0)p(1)?....
..h(255)p(256)
p(255)?h(?255)p(0)?......h(0)p(2
55)?......h(255)p(510)
这些值不像理论情况补以零值,而是按下法扩充:
p(0)?p(1)
2
p(254)?p(255)p?
2
p'?
所以,若采样点为
N
,则扩充后的点数为
3N?2
。
为了便于计算机实现,可将式作如下变换:
n??255?1
n?256510
p(n)?
l??256
?
256
p(n
?l)h(l)
?
?
p(n?l'?256)h(l'?256)
l?1511
这里引用新坐标
l'?l?256
,即
l??256
正好是
l'
的原点。我们有
p(l)?p(l'?256)?p
s
(l')
h(l)?h(l'?256)?h
s
(l')
于是式变为:
511
l'?1
p(n)?
?
p
s
(n?l')h
s
(l')
相应地式变为:
p(0)?h
s
(1)p
s
(1)?......h
s
(256)p
s
(
256)?......h
s
(511)p
s
(511)
p(1)?
h
s
(1)p
s
(2)?......h
s
(256)p<
br>s
(257)?......h
s
(511)p
s
(512)
p(255)?h
s
(1)p
s
(256)
?......h
s
(256)p
s
(511)?......h
s
(511)p
s
(766)
射束计算与内插
目前CT中常用的内
插方式是以离散滤波函数
h(n)
与投影数据
p(n,m)
作离散卷积,得<
br>到滤波后的投影数据
p(n,m)
,其中
m
为对应于旋转角度
?
的变量,
?
?
为每次旋转的角度,
?
?m?
?<
br>。然后进行线性内插,求得
p(x
r
,m)?p(n,m)*?(x
r
)
。
对于空间某点
(x
?
,y
j
)
在某视角
?
?
?
m
下,必有一
x
r,m
随之而定,如图所示:
j
x
r,m
?x
?
cos
?
?y
j
sin
?
y
(x
?
,y
j
)
?
x
r
?
0
x
i
图 射束计算
由于
(x
?
,y
j
)
是空间中任意一点坐标,故由式得到的x
r,m
并不一定正好为
d
的整数倍,
它可能位于
n<
br>0
d
与
(n
0
?1)d
之间,即:
x
r,m
?(n
0
?
?
)d,
我们采用线
性内插,则可得到:
0<
?
<1
p[(n
0
?
?
)d]?p(n
0<
br>)?
p
?
(n
0
?1)d
?
?p(n
0
d)
d
?(1?
?
)p(n
0
)?
?
p
?
(n
0
?1)d
?
(x
r,m
?n
0
d)
在前面的推导过程中,
x
r
以
x?0,y?0
为坐标原点,而按照射束计算的实际配置,以
直径的一端
作为原点为宜,这样在计算中可避免负值(如图所示,红线代表零号射束,圆表
示扫描所需最小直径的圆
,设其直径为
D
)。这一新的坐标记以
x
r
,这样对于任一像素(x
?
,y
j
)
及视角
?
,有:
其中M、N为待重建物体长和宽的象素个数
所以有:
x
r
?x<
br>i
cos
?
?y
j
sin
?
?(i?(M?
1)2)cos
?
?(j?(N?1)2)sin
?
x
r
?x
r
?D2
五、 流程图及程序
程序流程
读入模拟待重建图像
计算投影值,外
层角度循环
内层图像行列循环,从第一
行第一列遍历图像每一点
计算图像象元在旋转坐标系下的横坐
标,
再参照图所示圆进行坐标变换,得
到以圆的直径一端为原点的坐标
判断图像象元落在探测器上的位置
如果该象元被射线穿过,则将该象元
的像素值加到对应的探测器象元上
判断图像上所有行是否计算完毕
行数
否
是
加一
判断图像上所有列是否计算完毕
列数
否
加一
是
判断旋转角度是否达到180度
否
是
旋转至下
一角度
将每一角度下的投影值与
R-L滤波器进行卷积滤波
反投影,外层角度循环
内层图像行列循环,从第一
行第一列遍历图像每一点
计算图像上的点在旋转坐标系下
的横
坐标,再参照图所示圆进行坐标变换,
得到以圆的直径一端为原点的坐标
判断图
像点落在探测器上的位置,
将投影值进行内插,得到该点的投
影值,将投影值进行反投影
判断图像上所有行是否计算完毕
行数
否
加一
是
判断图像上所有列是否计算完毕
列数
否
加一
是
判断旋转角度是否达到180度
否
是
旋转至下
一角度
反投影结束,将反投影值进行
调整,压缩至0-255之间
附录一(公式推导):
1.
f(t)?0
有
n<
br>个互不相同的实根
t
1
,t
2
,???t
n
,则有
[2]
?
?
f(t)
?
?
?i
n
1
?
(t?t
i
)
f'(t
i
)
证明:把
f(t)
在
t?t
i
附近足够小的邻域内用泰勒级数展开,注意到
f(t
i
)?0
并忽
略高次项
1
f(t)?f(t
i
)?f'(t
i
)(t?
t
i
)?ft
i
)(t?t
i
)
2
???
??
2
?f'(
t
i
)(t?t
i
)
由尺度变换可知在
t?t
i<
br>附近
?
?
f(t)
?
表示为:
?
?
f(t)
?
?
?
?
f'(t
i
)(t?t
i
)
?
?
若
f'(t)
在
n
个单根t?t
i
(i?1,2???n)
处都不为零,则
?
?
f(t)
?
?
?
i?1
n
1
?
(t?t
i
)
f'(t
i
)
1
?
(t?t
i
)
f'(t
i
)
对应到本文第四页中式有
?
?
rc
os(
?
?
?
)
?
?
?
(
??
?
0
)
?
(
?
?
?
0)
?
rcos'(
?
??
0
)rsin'(
?
?
?
0
)
其中
?
0
是
a(
?
)?0
的唯一的根,
rco
s(
?
?
?
0
)?0
所以有,
sin(
?
?
?
0
)?1
,这样式变为
以上证明过程可以参考《信号与系统》 郑君里 应启珩 杨为理 第二章 77页
2.
f(x,y)?
?
?
rcos(
?
?
?)
?
?
?
(
?
?
?
0
)
1
?
x
2?y
2
为
x?y
平面上的函数,其二维傅里叶变换为
[3,4]
F(
?
,
?
)?
函数二维傅里叶变换的定义为:
1
??
?
?
22
F(
?
,
?
)?
?
采用极坐标,即:
则式变为:
??
??
?
??
??
f(x,y)e
-j2
?
(
?
x?
?
y)
dxdy
x?rcos
?
?
?
?<
br>cos
?
y?rsin
?
?
?
?
sin?
F(
?
,
?
)?
?
2
?
0
?
??
0
rf(r,
?
)e
-j2
??
rcos(
?
?
?
)
drd
?
因为
f(r,
?
)?1r
具有圆对称性,即
f
与
?
无关,于是写出:
f
R
(r)?f(r,
?
)
则:
2
?
??
F(
?
,
?
)?
?
利用贝塞尔函数关系式:
0
?
0
rf
R
(r)e
-j2
??
rcos(
?
?
?
)
drd
?
其中
J
0
(???)
是第一类零阶贝塞尔函数,于是:
?
2
?
0
e
-jacos(
?
?
?
)
d
?
?2
?
J
0
(a)
F(
?
,
?
)?2
?
?
rf
R<
br>(r)J
0
(2
??
r)dr
0
??
上式右端与
?
无关,从而可以写成:
F(
?
,
?
)?F
P
(
?
)
上式中
F
P
(
?
)
称为
f
R(r)
的零阶汉克尔变换,查常用的汉克尔变换表对于
f
R
(r)?1r
,其汉
克尔变换为
F
P
(
?
)?1
?
所以有:
F(
?
,
?
)?
参考文献:
1
??
?
?
22
[1] 庄天戈.《计算机在生物医学中的应用(第二版)》.科学出版社 第四章124-157
[2] 郑君里,应启珩,杨为理.《信号与系统(第二版)》.高等教育出版社 第二章 77
[3](美)布拉斯维尔(Bracewell,.).《傅里叶变换及其应用》人民邮电出版社
第十二章 286,292
[4] 刘培森.《应用傅里叶变换》北京理工大学出版社
第一章 19-22
附录二(与程序对应所用图形):
边界用图
图1.程序中模拟所画图形
图2. 每一探测器象元内射束排列,计算每条射束
图3. 内插所用图形
玄米-messenger
邮件英语-莩怎么读
华氏温度-聒噪怎么读
钱英文-凡尔赛和约
dure-腹黑是什么意思
留学咨询机构哪家靠谱-糜的拼音
可爱的日语-增字开头的成语
十六英语怎么读-制裁什么意思
-
上一篇:Photoshop制作创意漂亮公告板
下一篇:WPS表格快捷键