平行束反投影重建

一、 平行束反投影重建算法平行束 重建采用的是平移加旋转的扫描方式，如图所示，射线源在某一角度下水平移动，将物体全部照射后旋转一角度，如此重复，在这个过程中探测器相应地运动以接收X 射线。

1、反投影重建算法的物理概念： 断层平面中某一点的密度值可以看作是这一平面内所有经过该点的射线的投影值之和（的均值）。

整幅重建图像可以看作是所有方向下的投影累加而成。

射线标号示于图中，像素值（代表密度）分别1x ,2x ,3x ,4x ，赋值如下:15x =,20x =,32x =,418x =根据投影的定义(某条射线投影值为该条射线穿过的所有的像素值之和)，每条射线的投影i p (1,2i =)为：1215p x x =+=, 23420p x x =+=,3137p x x =+= 42418p x x =+=, 532p x ==, 61423p x x =+=720p x ==根据反投影重建算法的物理意义，重建图像中各像素，得到：113635x p p p =++=，214723x p p p =++=， 323529x p p p =++=，424661x p p p =++=(a) 原图像像素值 (b)反投影重建后图像 (c)求平均后图像图 反投影示例重建后的图像如图(b)所示，可以看出原图像中像素值不为零的点反投影重建后仍较突出，但原图中像素值为零的点，经反投影重建后不再为零，即有伪迹。

有时为了使重建后图像的像素值更接近于原图的像素值，在求反投影时，把数据除以投影的数目（即射线数），图 平行束平移加旋转图 断层像素值和射线如图(c)所示。

因此有：,11pn k k ii px pn ==∑该式可作为反投影重建算法的计算式。

其中k x 表示像素k 的值，i k p ,表示经过像素k的第i 条射线投影，p n 表示图像内的射线条数。

图(a)表示空间中一个孤立点源A ，密度为1。

经过A 点的三条射线也示于图中。

射线束理论上可以很多，取三条示意。

不经过A 点的射线投影为零，经过A 点的射线投影值均为1，1231p p p ===。

(a) 孤立点源A 及三条射线 (b)相应的反投影重建图，有星状伪迹图 孤立点源的反投影重建及星状伪迹经反投影重建后，得到A 点的像素值为123()/31A f p p p =++=。

A 点以外的像素值原来为零，经反投影重建后不等于零，而是等于1/3。

所以，经反投影重建后的图像除保留A 点的像外，还有像素值为1/n 的灰雾背景，后者称为星状伪迹。

产生星状伪迹的原因在于：反投影的本质是把取自有限物体空间的投影均匀地回抹（反投影）到射线所及的无限空间的各点上，包括原先像素值为零的点。

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

图中，(,)x y 为固定坐标，(,)r r x y 为旋转坐标，(,)r θ为极坐标，设位于坐标原点0,0x y ==的点源(,)x y δ为x y -断面中唯一的像点，扫描方式为平移加旋转。

即射线先平行移动，从物体的一侧移向另一侧；然后旋转一个角度φ∆，如此继续，直到累计的旋转角达到180φ-∆。

为止。

为了计及这一几何要求我们设置一旋转坐标系统r r x y -，它绕原点转动使投影总是沿着r y 方向，r rx y -的原点与x y -的原点相同。

二者的夹角为φ，不同的φ代表不同的投射方向。

投影线的位置可由(,)r x φ完全确定。

设φ为离散取值，如i φφ=，则相应的投影为：()(,)(,)(,)()()()|(cos())i i r r i r r r r r r rr r r r i p x p x f x y dy x y dy x y dy x r φφφφδδδδδθφ+∞+∞-∞-∞+∞=-∞======-⎰⎰⎰其中cos()r x r θφ=-，这可以从图的几何关系很容易得出。

若n φφ=，则相应的投影为：cos n n p r φδθφ=-[()]根据反投影重建的定义式,点(,)r r x y 的图像在所述坐标系统中表示为：[][]cos()11111(,)(,)()|cos()1cos()i r i iiN N r r r x r ii i N ii f r f x y p x p x r N N p r φφφφθφφφφφθθφθφφπ=-======-=-∆∑∑∑其中，N φ为投影数，/N φφπ∆=。

若在有限区间内将射线增至不相重的无限条，即连续取投影，则有：[]01(,)cos()f r p r d πφθθφφπ=-⎰在忽略射束硬化的情况下，φ在(,2)ππ区间内的投影值等于(0,)π区间内的投影值。

在输入图像为点源的情况下，由及可得：(,x y δ(,)h x y 重建后图像原像图 “反投影重建”成像系统图平移加旋转扫描方式所用坐标系[]01(,)cos()h r r d πθδθφφπ=-⎰因为，[]00()()'()a a δφφδφφ-=(此公式推导可参见附录一)，其中，0φ是()0a φ=的唯一的根。

令()cos()a r φθφ=-，则有：00()cos()0a r φθφ=-=因为0r ≠，所以0sin()1θφ-=，于是有：[]0000()11(,)cos()sin()1111sin()h r r d d r r r ππδφφθδθφφφππθφπθφπ-=-=-===-⎰⎰可见，相应于反投影重建算法的系统，它的点扩展函数不是δ函数，系统不是完美的。

虽然在0r =处能反映原图是点源的情况，但在0r ≠处，像素值0≠，上式定量地描述了反投影重建算法星状伪迹的本质。

若原像为(,)x y μ，则将原像取投影后再按反投影算法得到的重建图像为：(,)(,)(,)f x y x y h x y μ=**，其中**表示二维卷积。

要去除反投影算法的星状伪迹，可以在输出端加一滤波器，使加了滤波器后的反投影重建等效成像系统的点扩展函数(PSF)为(,)x y δ。

设滤波器的PSF 为(,)q x y ，相应的传递函数为[]2(,)(,)Q F q x y ξη=，要求(式的推导可参见附录一)：(,)1ξη=或：(,)Q ξηπρ=这是一个二维滤波器，实现起来较麻烦。

若ρ的变化范围扩至∞，则根本不能实现。

但它提供了去除星状伪迹的一个方向。

(,)(,)q x y x y δ=三、 滤波反投影重建算法反投影重建算法的缺点是引入星状伪迹，即原来图像中密度为零的点，重建后不一定为零，从而使图像失真。

去除伪迹的方法如图所示。

(a)(b)(a)方案之一，先“反投影”再二维滤波；(b) 方案之二，先对投影数据作一维滤波，再“反 投影”图 星状伪迹的去除投影定理投影定理（或中心切片定理）：在非衍射源情况下，某图像(,)f x y 在视角为φ时投影()r p x φ的一维傅氏变换给出(,)f x y 的二维傅氏变换^12(,)(,)A A ωωρφ=的一个切片。

切片与1ω轴相交成φ角，且通过坐标原点。

即：^1()(,)|r F p x A φφρφ⎡⎤=⎣⎦固定 图阐明投影定理 证明：()f x,y 的二维傅立叶变换为12()12,)()-j x y A(f x,y e dxdy ωωωω+∞+∞+-∞-∞=⎰⎰由图可以得到如下几何关系cos sin sin cos r r x x y y φφφφ⎡⎤⎡⎤⎡⎤=⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦输出图像输出图像1p φ以r y 为投影轴的投影为：()()()r r r r r r p x f x,y dy f x ,y dy φ+∞+∞-∞-∞==⎰⎰它的一维傅里叶变换为：1(cos sin )[()]()()()r r -j2x r r r-j2x r r r r r -j2x y F p x p x e dx f x ,y e dx dy f x,y e dxdyπρφφπρπρφφ+∞-∞+∞+∞-∞-∞+∞+∞+-∞-∞===⎰⎰⎰⎰⎰式中采用2πρ可以使后续反变换计算中的系数得以简化。

在推导上式的过程中，我们利用了下列关系：()()r r r f x,y f x ,y =以及cos sin sin cos r rr r r r x x xy dx dy dxdy dxdy y y xyφφφφ∂∂∂∂===∂∂-∂∂ 由式可知，当12,ωω的取值满足条件122cos 2sin ωπρφωπρφ=⎫⎬=⎭时，有121212()1cos sin 12cos sin [()]()|(,)|-j x y r 2222F p x f x,y e dxdy A ωωφωπρφωπρφωπρφωπρφωω∞∞+=-∞-∞=====⎰⎰上式说明：12,ωω不是独立的而是受到式的约束，其值必须局限于直线21(tg )ωφω=上。

根据投影定理，投影图像重建问题原则上可按如下流程求解： （1）采集不同视角下的投影；（2）求出各投影的一维傅氏变换，此即图像二维傅氏变换过原点的各切片，理论上是连续的无穷多片；（3）将上述各切片汇集成图像的二维傅氏变换； （4）求二维傅氏反变换的重建图像。

卷积反投影重建（即滤波反投影）待重建图像为),(y x a ，它的二维傅氏变换为^12(,)(,)A A ωωρφ=。

根据中心切片定理，^(,)A ρφ可通过),(y x a 在不同视角φ下的投影()r p x φ的一维傅氏变换求得。

即：待建图像：^121(,)(,)()()(,)r A A F p x PP φφωωρφρρφ⎡⎤====⎣⎦12^1212()12122^2cos()02cos()02cos()0(,)(,)(,)1(,)4(,)(,)(,)i x y i r i r i r a r a x y F A A e d d A e d d P e d d d P e d ωωππρθφππρθφππρθφθωωωωωωπρφρρφρφρρφφρρφρ-∞∞+-∞-∞∞--∞∞--∞∞--∞======⎰⎰⎰⎰⎰⎰⎰⎰[]因为cos()r x r θφ=-，所以有：122(cos sin )22cos()r x y x y x r ωωπρφφπρπρθφ+=+==-同时：12d d J d d ωωρφ=11222//2cos 2sin 4//2sin 2cos J ωρωφπφπρφπρωρωφπφπρφ∂∂∂∂-===∂∂∂∂先来看该式的第二个积分：[]22cos()cos()cos()cos()(,)(,)|()(,)|(,)|cos(),rrr r i x i r x r r r x r r x r P e d P e d h x p x g x g r πρπρθφθφθφθφρρφρρρφρφφθφφ∞∞-=--∞-∞=-=-==*==-⎰⎰式中：(,)()(,)r r r g x h x p x φφ=*式的物理意义是投影(,)r p x φ经过传递函数为1[()]r F h x ρ=的滤波器后得到的修正后的投影(,)r g x φ在满足cos()r x r θφ=-时的值。