平行束反投影重建

一、平行束反投影重建算法平行束重建采用的是平移加旋转的扫描方式，如图1.1所示，射线源在某一角度下水平移动，将物体全部照射后旋转一角度，如此重复，在这个过程中探测器相应地运动以接收X射线。1、反投影重建算法的物理概念：断层平面中某一点的密度值可以看作是这一平面内所有经过该点的射线的投影值之和（的均值）。整幅重建图像可以看作是所有方向下的投影累加而成。射线标号示于图1.2中，像素值（代表密度）分别1x,2x,3x,4x，赋值如下:15x,20x,32x,418x根据投影的定义(某条射线投影值为该条射线穿过的所有的像素值之和)，每条射线的投影ip(1,2i)为：1215pxx,23420pxx,3137pxx42418pxx,532px,61423pxx720px根据反投影重建算法的物理意义，重建图像中各像素，得到：113635xppp，214723xppp，323529xppp，424661xppp521803529612654.18.73.3(a)原图像像素值(b)反投影重建后图像(c)求平均后图像图1.3反投影示例重建后的图像如图1.3(b)所示，可以看出原图像中像素值不为零的点反投影重建后仍较突出，但原图中像素值为零的点，经反投影重建后不再为零，即有伪迹。有时为了使重建后图像的像素值更接近于原图的像素值，在求反投影时，把数据除以投影的数目（即射线数），50218(1)(2)(3)(4)(5)(6)(7)X射线管平移平移探测器图1.1平行束平移加旋转图1.2断层像素值和射线如图1.3(c)所示。因此有：,11pnkkiipxpn(1.1)该式可作为反投影重建算法的计算式。其中kx表示像素k的值，ikp,表示经过像素k的第i条射线投影，pn表示图像内的射线条数。图1.4(a)表示空间中一个孤立点源A，密度为1。经过A点的三条射线也示于图中。射线束理论上可以很多，取三条示意。不经过A点的射线投影为零，经过A点的射线投影值均为1，1231ppp。(a)孤立点源A及三条射线(b)相应的反投影重建图，有星状伪迹图1.4孤立点源的反投影重建及星状伪迹经反投影重建后，得到A点的像素值为123()/31Afppp。A点以外的像素值原来为零，经反投影重建后不等于零，而是等于1/3。所以，经反投影重建后的图像除保留A点的像外，还有像素值为1/n的灰雾背景，后者称为星状伪迹。产生星状伪迹的原因在于：反投影的本质是把取自有限物体空间的投影均匀地回抹（反投影）到射线所及的无限空间的各点上，包括原先像素值为零的点。A(1)(2)(3)A二、反投影重建算法的数学描述我们把“取投影”、“反投影重建”、“重建后图像”这些环节看作是一个以原像为输入，重建后图像为输出的成像系统，如图2.1所示，先来求该系统的点扩展函数PSF。图2.2中，(,)xy为固定坐标，(,)rrxy为旋转坐标，(,)r为极坐标，设位于坐标原点0,0xy的点源(,)xy为xy断面中唯一的像点，扫描方式为平移加旋转。即射线先平行移动，从物体的一侧移向另一侧；然后旋转一个角度，如此继续，直到累计的旋转角达到180。为止。为了计及这一几何要求我们设置一旋转坐标系统rrxy，它绕原点转动使投影总是沿着ry方向，rrxy的原点与xy的原点相同。二者的夹角为，不同的代表不同的投射方向。投影线的位置可由(,)rx完全确定。设为离散取值，如i，则相应的投影为：()(,)(,)(,)()()()|(cos())iirrirrrrrrrrrrripxpxfxydyxydyxydyxr(2.1)其中cos()rxr，这可以从图2.2的几何关系很容易得出。若n，则相应的投影为：cosnnpr[()](2.2)根据反投影重建的定义式(1.1),点(,)rrxy的图像在所述坐标系统中表示为：cos()11111(,)(,)()|cos()1cos()iriiiNNrrrxriiiNiifrfxypxpxrNNpr(2.3)其中，N为投影数，/N。若在有限区间内将射线增至不相重的无限条，即连续取投影，则有：01(,)cos()frprd(2.4)在忽略射束硬化的情况下，在(,2)区间内的投影值等于(0,)区间内的投影值。在输入图像为点源的情况下，由(2.1)及(2.4)可得：取投影反投影重建(,)xy(,)hxy重建后图像原像图2.1“反投影重建”成像系统rx0r(,)frxyry图2.2平移加旋转扫描方式所用坐标系01(,)cos()hrrd(2.5)因为，00()()'()aa(此公式推导可参见附录一)，其中，0是()0a的唯一的根。令()cos()ar，则有：00()cos()0ar因为0r，所以0sin()1，于是有：0000220()11(,)cos()sin()11111sin()hrrddrrrxy(2.6)可见，相应于反投影重建算法的系统，它的点扩展函数不是函数，系统不是完美的。虽然在0r处能反映原图是点源的情况，但在0r处，像素值0，上式定量地描述了反投影重建算法星状伪迹的本质。若原像为(,)xy，则将原像取投影后再按反投影算法得到的重建图像为：(,)(,)(,)fxyxyhxy，其中表示二维卷积。要去除反投影算法的星状伪迹，可以在输出端加一滤波器，使加了滤波器后的反投影重建等效成像系统的点扩展函数(PSF)为(,)xy。设滤波器的PSF为(,)qxy，相应的传递函数为2(,)(,)QFqxy，要求(式2.7的推导可参见附录一)：做二维傅氏变换，得：221(,)1Q或：22(,)Q(2.8)这是一个二维滤波器，实现起来较麻烦。若的变化范围扩至，则根本不能实现。但它提供了去除星状伪迹的一个方向。221(,)(,)qxyxyxy(2.7)三、滤波反投影重建算法反投影重建算法的缺点是引入星状伪迹，即原来图像中密度为零的点，重建后不一定为零，从而使图像失真。去除伪迹的方法如图3.1所示。(a)(b)(a)方案之一，先“反投影”再二维滤波；(b)方案之二，先对投影数据作一维滤波，再“反投影”图3.1星状伪迹的去除3.1投影定理投影定理（或中心切片定理）：在非衍射源情况下，某图像(,)fxy在视角为时投影()rpx的一维傅氏变换给出(,)fxy的二维傅氏变换^12(,)(,)AA的一个切片。切片与1轴相交成角，且通过坐标原点。即：^1()(,)|rFpxA固定(3.1)图3.2阐明投影定理证明：()fx,y的二维傅立叶变换为12()12,)()-jxyA(fx,yedxdy(3.2)由图2.2可以得到如下几何关系cossinsincosrrxxyy(3.3)输入图像（原图像）输出图像（同原图像）取投影反投影重建二维滤波器输入图像（原图像）取投影输出图像（同原图像）反投影重建一维滤波器1()rpx12(,)A2rxyryrxx0二维傅立叶变换一维傅氏变换以ry为投影轴的投影为：()()()rrrrrrpxfx,ydyfx,ydy它的一维傅里叶变换为：1(cossin)[()]()()()rr-j2xrrr-j2xrrrrr-j2xyFpxpxedxfx,yedxdyfx,yedxdy(3.4)式中采用2可以使后续反变换计算中的系数得以简化。在推导上式的过程中，我们利用了下列关系：()()rrrfx,yfx,y以及cossinsincosrrrrrrxxxydxdydxdydxdyyyxy(3.5)由式(3.2)可知，当12,的取值满足条件122cos2sin(3.6)时，有121212()1cossin12cossin[()]()|(,)|-jxyr2222Fpxfx,yedxdyA(3.7)上式说明：12,不是独立的而是受到式(3.6)的约束，其值必须局限于直线21(tg)上。根据投影定理，投影图像重建问题原则上可按如下流程求解：（1）采集不同视角下的投影；（2）求出各投影的一维傅氏变换，此即图像二维傅氏变换过原点的各切片，理论上是连续的无穷多片；（3）将上述各切片汇集成图像的二维傅氏变换；（4）求二维傅氏反变换的重建图像。3.2卷积反投影重建（即滤波反投影）待重建图像为),(yxa，它的二维傅氏变换为^12(,)(,)AA。根据中心切片定理，^(,)A可通过),(yxa在不同视角下的投影()rpx的一维傅氏变换求得。即：待建图像：^121(,)(,)()()(,)rAAFpxPP12^1212()12122^2cos()02cos()02cos()0(,)(,)(,)1(,)4(,)(,)(,)ixyiririraraxyFAAeddAeddPedddPed[](3.8)因为cos()rxr，所以有：122(cossin)22cos()rxyxyxr同时：12ddJdd11222//2cos2sin4//2sin2cosJ(3.9)先来看该式的第二个积分：22cos()cos()cos()cos()(,)(,)|()(,)|(,)|cos(),rrrrixirxrrrxrrxrPedPedhxpxgxgr(3.10)式中：(,)()(,)rrrgxhxpx(3.11)式(3.10)的物理意义是投影(,)rpx经过传递函数为1[()]rFhx的滤波器后得到的修正后的投影(,)rgx在满足cos()rxr时的值。将(3.11)代入(3.8)，得到：^0(,)[cos(),]argrd(3.12)称为滤波反投影方程，其物理意义是经过给定点(,)r的所有滤波后的投影在0~范围内的累加—反投影重建，得出(,)r点的像素值。可见，滤波（卷积）反投影算法的具体包含三大步：(1)把在固定视角i下测得的投影(,)rpx经过滤波，得到滤波后的投影(,)rgx；(2)对每一个i，把(,)rigx反投影于满足cos()rixr的射线上的所有各点(,)r；(3)将步骤(2)中的反投影值对所有0进行累加（积分），得到重建后的图像。3.3滤波函数与内插函数的选取原则由式(3.10),(3.11)可知，先要把投影值(,)rpx做一维滤波。滤波器的传递函数为()H。可以证明这一理想的滤波函数是不可实现的。我们的目的是：选取这样的滤波函数，使它既可以实现，又有足够的精度。为此，结合成像的具体情况考虑。(1)投影数据的高频（空间分辨率）分量幅度很小；(2)投影数据的采集是