图像处理(第九章).doc

9. 图像的重建9.1 投影重建概述9.1.1投影重建方式利用投影重建方式工作的系统有很多种类，下面是几个例子。1.透射断层成像在透射断层成像（transmission computed tomography，TCT，简称CT）系统中，从发射源射出的射线穿透物体到达接受器。射线在通过物体时被物体吸收一部分，余下部分被接受器接受。由于物体各部分对射线的吸收不同，所以接受器获得的射线强度实际上反映了物体各部分对射线的吸收情况。TCT是常见和简单的CT。因常使用X射线（不会产生衍射），所以也称XCT。如果用代表射线源的强度，代表在沿射线方向物体点的线性衰减系数/因子，代表辐射的射线，代表穿透物体的射线强度，那么有：（9.1.1）如果物体时均匀的，则：（9.1.2）其中，代表穿透物体后的射线强度，代表没有物体时射线强度，是射线在物体内部的长度，代表物体的线性衰减系数。近年来CT所需的扫描时间迅速减少，而与此同时其空间分辨率不断增加。1972年CT的空间分辨率为，扫描时间需5分钟，而现在获取的空间分辨率的图像只需用不到1秒钟的时间。使用超声的CT也根据TCT的原理工作，不过发射源发射的是超声。由于超声并不是严格沿直线传播的，且物体对其的影响与内部结构有关，所以常需要使用非线性模型。2.发射断层成像在发射断层成像（emission computed tomography，ECT）系统中，发射源在物体内部（可将被检测的物体本身看作信号源）。一般是将具有放射性的离子注入物体内部，从物体外检测其放射出来的量。通过这种方法可以了解离子在物体内的运动情况和分布，从而可以检测到与生理有关的状况。现在常用的ECT主要有两种：PET（positron emission tomography）和SEPCT（single photon emission CT）。PET的历史可追溯到上世纪50年代放射性物质的成像开始得到使用。与其他截面成像技术（如CT和下面的MRI）相比，PET具有很大的敏感度，具有检测和显示纳摩尔（nanomolar）级的能力。图9.1.1给出PET成像系统的构成示意图。PET采用在衰减时放出正电子的放射性离子，放出的正电子很快与负电子相撞湮灭而产生一对光子并以相反的方向射出。所以相对放置的两个检测器接受到这两个光子就可以确定一条射线。检测器围绕物体呈环形分布，相对的两个检测器构成一组检测器对以检测有一对正负电子产生的两个光子。图9.1.1PET成像系统构成示意图如果两个光子被一对检测器同时记录下来，那么产生这两个光子的湮灭现象肯定发生在连接两个检测器的直线上。为避免射线的影响，一般要在检测到10万或更多的湮灭事件后才用断层重建方法计算正电子发射的轨迹。对事件的投影记录数据可表示为：（9.1.3）其中，是投影数据，是对射线的线性衰减系数，是同位素的分布函数。PET利用不同角度的1-D投影来重建2-D图像，如果利用不同角度的2-D投影则可重建3-D图像。目前，PET的空间分辨率已达2mm3mm（接近极限），有些成像系统有多达20000个检测器。SPECT是一种结合核医学成像技术和断层重建技术的成像技术。所成的像与PET类似都反映了成像物体的信息。两者给出的信息都基于放射性离子的空间分布。在SPECT中，任何在衰减中能产生射线的放射性离子都可使用。SPECT可利用寿命（半衰期）较长因而比较容易获得的试剂（而PET因为所用离子寿命短需要有加速器）。为确定射线方向，要用能阻止射线偏移的准直器来定向采集中子以确定射线方向。SPECT使用围绕物体从不同角度得到的2-D核医学图像，并利用从多个投影重建的方法提供3-D放射性分布信息。与TCT不同，SPECT使用的放射源是在物体的内部。从TCT获得的是不同材料的不同衰减系数，而从SPECT获得的是放射性物质在物体内的分布。图9.1.2给出SPECT成像系统的构成示意图。将放射性物质注入物体内，不同的材料（如组织或器官）吸收后会发射射线光子。一定方向的光子可穿过准直器到达晶体，在那里射线光子转化为能量较低的光子并由光电倍增器转化为电信号。这些电信号提供了光子与晶体作用的位置，放射性物质的3-D分布就转化为2-D的投影图像。图9.1.2SPECT成像系统构成示意图SPECT的敏感度比PET要低很多，原因是为确定射线方向需要使用铅制造的准直器（collimator），这样就限制了断层成像装置的立体角。一般发射断层成像的敏感度可表示为：（9.1.4）其中是检测器的面积，是检测器的效率（对SPECT是1而对PET是2），是衰减系数（一般在0.2到0.6之间），是断层的半径。一般PET的敏感度与SPECT的敏感度之比约为150除以分辨率，当分辨率为7.5mm时这个比值为20。3.反射断层成像反射断层成像（reflection CT，RCT）也是利用投影重建的原理工作的。常见的一个例子是雷达系统，其中的雷达图示物体反射的回波所产生的。例如在前视雷达（FLR）中，雷达发射器从空中向地面发射无线电波。雷达接受器在特定角度所接受到的回波强度是地面反射量在一个扫描阶段中的积分。图9.1.3给出非聚焦合成孔径雷达（synthetic aperture radar，SAR）成像的示意图。在合成孔径雷达成像中，雷达是运动的而目标是不动的（利用它们之间的相对运动产生大的相对孔径以提高横向分辨率）。设为雷达（载体）沿轴的运动速度，为有效积累时间，为电波波长。两个点目标沿雷达运动方向分布，目标位于雷达孔径正侧视（雷达波束指向与雷达运动方向互相垂直）的中心线上（轴），目标与目标间的位移量为。雷达与目标间的最近距离为，此时定义为时间零点，。设在前后距离的变化量为。当时。在目标处回波信号的双程（电波在天线和目标间来回传播）超前相位为：（9.1.5）在目标处回波信号的双程超前相位为：（9.1.6）如果发射信号频率足够高，回波信号可认为是连续的，此时可对时间段到进行积分来处理回波信号。进一步设在积分时间内为均匀发射，则在目标处的回波响应为：（9.1.7）4.磁共振成像磁共振成像（magnetic resonance imaging，MRI）在早期称为核磁共振（nuclear magnetic resonance，NMR）。它的工作原理简介如下。氢核以及其他具有奇数个质子或中子的核，包含具有一定磁动量或旋量的质子。如果把它们放在磁场中，它们就会像陀螺在地球重力场中一样在磁场中进动。一般情况下质子在磁场中是随意排列着的，当一个适当强度和频率的共振场信号作用于物体时，质子吸收能量并转向与磁场相交的朝向。如果此时将共振场信号除去，质子吸收的能量将释放并可被接受器检测到。根据检测到的信号就可以确定质子的密度。通过控制所用的共振场信号和磁场强度，可每次检测到沿通过目标中一条线的信号。换句话说，检测到的信号时MRI信号沿直线的积分。每个MRI成像系统都包含有磁场子系统、发射/接收子系统、计算机图像重建和显示子系统。在磁场子系统中，有纵向的主磁场、非均匀磁场和横向的射频磁场，在它们的共同作用下，成像物体会产生磁共振信号，该信号能被设备中的检测线圈所测得。假设射频场的作用时间远小于自旋原子核的横向和纵向弛豫时间常数，则可以将磁共振信号表示为：（9.1.8）其中是被横向和纵向弛豫时间常数等物理参数加权的原子核的自旋密度分布函数，是磁共振信号的射频接收线圈的灵敏度分布函数，是原子核自旋进动的拉莫尔（Larmor）频率的空间分布随时间变化的函数，这里，是原子核的旋磁比，是主磁场强度，是时变的非均匀磁场；是成像物体所处的空间区域。磁共振成像就是根据时变非均匀磁场和射频磁场及其激励而产生的磁共振信号来重建物体的自旋密度分布函数。数学上将磁共振成像看作是一个逆问题，即在已知，的条件下求解的积分方程（9.1.8）。最早提出的磁共振成像方法将非均匀磁场设计为线性梯度磁场，从而把积分方程（9.1.8）简化为Radon变换，然后通过求解Radon逆变换重建物体的自旋密度图像。磁共振成像的分辨率近年有了很大提高。在1980年时，分辨率约为若干个毫米，而现在已接近十分之一毫米。顺便指出，还有一种方法将线圈接收到的信号表示为物体中的旋量密度（density of spins）的函数。对具有旋量密度为的物体，其对应的空间频率信号为：（9.1.9）积分区域为整个物体。和构成一对傅里叶变换对。5.电阻抗断层成像电阻抗断层成像（electrical impedance tomography，EIT）采用交流电场对物体进行激励检测。这种方法对电导或电抗的改变比较敏感。通过使用低频率的电流注入物体内部并测量在物体外表处的电势场（根据电导率分布利用有限元计算场域边界电压），再采用图像重建算法就可以重建出物体内部区域的电导和电抗的分布或变化的图像（基于边界测量值估计场域电导率分布）。目前，EIT是仅有的能对电导成像的方法。由于不同的生物组织或器官在不同的生理、病理条件下其电阻抗特性不同，所以EIT图像能反映组织或器官携带的病理和生理信息。因为EIT不使用核素或射线，只需要较少的电流，无毒无害，可作为一种无损伤检测技术为病人进行长期、连续的图像监护。常用EIT系统的工作模式分两类：注入电流式和感应电流式。注入电流式EIT采用注入激励和测量技术对成像区域的阻抗分布信息进行测量：在物体外表驱动电极上施加恒定交流激励，成像区域的等效阻抗反映在不同测量电极上测到的电压信号的幅值和相位上。采用解调技术可解调出被测信号中反映成像区域阻抗分布的部分信息。感应电流式EIT采用与物体外表非接触的激励线圈进行交流激励，从而在成像目标内部产生感应电流（涡流），从外表检测到感应电流场并从而计算出内部相应的阻抗分布或变化。注入电流式EIT需要在体表安放多个激励/测量电极（复合电极），受电极尺寸的限制，激励电极数目不可能无限增加，使独立测量数受到影响。感应电流式因采用非接触的线圈激励模式，可通过增加线圈的数目或移动线圈的位置等方法增加独立测量数目。相对于注入电流式EIT，感应电流式技术需要额外的功率放大器进行驱动，而且在成像区域产生感应电流的同时，在信号检测的输入回路也产生较大的感应电动势，对有用信号的检测有些影响。例9.1.1注入电流式EIT图像示例图9.1.4给出两幅利用注入电流式EIT重建出的图像，其中灰度取决于各处电阻抗的数值。图9.1.4电阻抗断层成像结果从数学角度来讲，EIT和各种CT有类似处，因为两种情况下都需要处理数据以产生反映物体内部结构的图像，而且成像常针对穿过物体的一个2-D截面进行。区别是EIT借助电流的扩散来获得电导的分布，这点与CT，PET，SPECT，MRI都不同。EIT有一些很吸引人的特征，包括用于EIT成像的技术很安全，简便。相对于CT，PET，SPECT，MRI等技术，EIT的分辨率要差。EIT的分辨率依赖于电极的数量，但可以同时接触到物体的电极数量常受到许多限制。EIT成像方式的主要缺点源于其非线性的病态问题，如果测量有很小的误差就有可能导致对电导的计算产生很大的影响。综上所述，如果测量到的数据具有物体某种感兴趣物理特性在空间分布的积分形式，那么就需要用投影重建的方法来获得物体内部的，反映不同物理特性的图像。9.1.2投影重建原理1.基本模型下面先介绍一下简单的从投影重建图像的模型（见图9.1.5）。这里用图像代表某种物理量在2-D平面上的分布。将需要投影重建的物质材料限制在一个无限薄的平面上，使得重建图像在任意点的灰度值正比于射线投影到的那个点所固有的相对线性衰减系数。图9.1.5投影重建示意图为讨论方便（且也符合实际情况），设在一个以原点为圆心的单位圆外为0。现考虑有一条由发射源到接受器的直线在平面上与在内相交。这条直线可用两个参数来确定：它与原点的距离；它与轴的夹角。如果用表示沿直线对的积分，借助坐标变换可得：（9.1.10）这个积分就是沿方向的投影，其中积分限取决于、和。当是单位圆时，设积分上下限分别为和，则：（9.1.11）如果直线落在外（与不相交），则：（9.1.12）可见式（9.1.10）表示的积分方程是有定义并可计算的。在实际的投影重建环境中，用表示需要被重建的目标，由确定的积分路线对应一条从发射源到接受器的射线。接受器所获得的积分测量值是。在这些定义下，从投影进行重建可描述为：对给定的，要确定。从数学上讲就是要解积分方程（9.1.10）。需要注意，一个定义在有限区域的函数可被无穷个投影所唯一确定，但不一定能被任何有限个投影所唯一确定。2.3-D重建利用断层重建可以实现真正的3-D成像。目标用体素的3-D数组来表示，而独立的投影给出2-D的投影图（数组）。每个投影来源于一个点，满足锥束几何（cone-beam geometry），与用于单片的扇束方法（fan-beam method）对应。实用的3-D锥束扫描CT系统使用面阵检测器，可看作扇束系统的一种推广。视线方向要将离开平面和进入3-D的朝向（用两个极角表示）包括进去。这并不需要根据两个极角旋转目标，因为使用锥束几何可提供不同角度的3-D投影线，正如扇束方法在2-D时做的那样（见9.3.1小节）。最好的重建是用一系列尽可能均匀覆盖3-D朝向的视角实现的。具体成像时可以使用不同的几何结构形式，最简单的一种是绕单个轴旋转目标，如图9.1.6。这种方法的优点是旋转准确，因为此时重建质量依赖于旋转中心的一致性。但另一方面，重建体素的失真可能会比较明显，特别在平行于旋转轴的方向和接近目标上下端点的地方。这种单轴旋转的方法最常用于光、中子或射线断层扫描中，因为此时目标可能相当大且射线沿各个方向穿过目标所通过的距离比较接近。3-D锥束扫描CT系统使用了面阵检测器，可看作扇束系统的一种推广。沿物体轴线方向的图像分辨率可利用螺旋状的扫描来提高，此时让目标在旋转的同时也沿轴线向上移动，如图9.1.7。对小尺寸的工业目标，这是最好的几何结构形式。目前典型的射线CT截面上的分辨率可达小于1mm，但截面的厚度常达好几个毫米。利用螺旋CT，则不需要一层层成像而是在3个方向上都可获得均匀（各向同性）的分辨率。图9.1.6通过绕单个轴旋转目标来实现立体成像图9.1.7螺旋状扫描实现立体成像 在连续体积数据采集中，需要将物体缓慢地穿过检测场，借助插值可获得3-D的数据。目前螺旋CT图像质量主要由射线管所能发射的能量所限制。9.2 傅里叶反变换重建1.基本步骤和定义傅里叶反变换重建法是一种变换重建方法。变换重建方法主要包括以下三个步骤：（1） 建立数学模型，其中已知量和未知量都是连续实数的函数；（2） 利用反变换公式解未知量；（3） 调节反变换公式 以适应离散、有噪声应用的需求。注意上面第（2）步中在解未知量时理论上可以有多个等价的公式来解。而第（3）步中，由于离散化时可采用不同的挖，所以理论上等价的方法对实际数据应用的结果会不同。在具体应用中，所测量到的数据对应于许多 个离散点（s,）上的估计值，而所重建的图像也是一个离散的2-D数组。在以下的讨论中，先考虑在S和上都均匀采样的情况。设在N个相差的角度上进行投影，在每个方向（对应某个投影角度）上用M个间距为的射线测量，定义整数M+和M-为：（9.2.1）为了保证一系列射线覆盖单位圆，需要选和。此时为平行投影的射线数据。另一种与此对应的是扇形扫描所产生的扇束散开线数据，将在后面介绍。设图像区被一个直角网格覆盖，其中和用类似于式（9.2.1）的方法定义（K为X方向上的点数），和也用类似于式（9.2.1）的方法定义（L为Y方向上的点数）。根据这些定义，一个重建算法就是要通过MN个测量值估计出在KL个采样点的。2.博里叶变换投影定理变换方法的基础是博里叶变换投影定理。设是对应第一个变换量s的1-D博里叶变换，即：（9.2.2）F(X,Y)是的2-D傅里叶变换：（9.2.3）那么可以证明如下投影定理：（9.2.4）即以角进行投影的博里叶变换等于的博里叶变换在博里叶空间处的值。换句话说, 在与X轴成角的直线上投影的博立叶变换是的博里叶变换在朝上角上的一个截面。上述结果在第5章讨论Randon变换时也涉及到（参见5.5.1小节）。对投影(Randon变换)的1-D博里叶变换可得到定义在博里叶空间的极坐标网络。这样就需要进行插值以获得直角坐标系统中的，然后才能通过2-D博里叶反变换得到。整个过程可见图9.2.1，其中R代表Randon变换，F代表博里叶变换。尽管这个方法看起来最简单，由于它需要进行插值和2-D博里叶反变换，所以计算量会相当大。3博里叶反变换重建公式根据博里叶变换投影定理很容易就可以得到博里叶反变换重建的公式。对公式(9.2.4)两边在直角坐标系中取博里叶反变换： (9.2.5)注意到是的博里叶变换，所以上式给出，计算的一个重建公式。这也是第一个在投影重建中得到应用的技术。实用中需要加一个窗，以把积分区限制在博里叶XY平面上的一个有限区域W，这样就得到的一个带线逼近：(9.2.6)现在考虑计算。实际中只在一系列q = qn角取值（qn代表nDq）可用在一系列采样点(mDs, qn)对g()的求和得到：(9.2.7)如果令R = kDR（k为整数，DR为采样间距），取DR = 1/(MDs) ，则：(9.2.8)根据可对任意的插值出，然后由式（9.2.6）可知：(9.2.9)于是可由下式确定：(9.2.10)进一步，如果令Dx = 1/(UDX)，Dy = 1/(VDY) ，则：(9.2.11)图9.2.2（a）右半部分给出在（博里叶）XY平面上的位置，用“”表示。在这些采样位置上可用式（9.2.8）计算G。这些点在平面上形成一个极坐标模式，各点在极坐标系中分布在等距的圆周上，在每个圆周上也是等角度分布的。图9.2.2（a）左半部中“+”点表示位置，而需要知道的是这些点的博里叶变换以根据式（9.2.11）计算。这些点在平面上形成一个直角坐标网络。可根据在极坐标系中的已知值进行插值（参见3.4.2小节）以得到这些网络点的估算值。上式只给出的限个(U个k , V个l)估计值。对这些参数的双重求和可用FFT快速计算。图9.2.2（a）右半部的极坐标模式对博里叶平面的采样不是很有效，一种可能的改进是采用图9.2.2（b）的模式。在图9.2.2（b）的右半部，点在平面上仍形成一个极坐标模式，但这里的点分布的更为均匀。原来在三个圆周上的分布成为在六个圆周角上的分布，相邻两个圆周上的点是交错分布的，所以点的分布更为均匀，更接近如图9.2.2（b）左半部直角坐标网络上的分布。综上所述，基于博里叶反变换的重建技术主要有三个步骤：图9.2.2 博里叶空间的直角和极坐标网络（1） 对以角的投影进行1-D博里叶变换，见式（9.2.8）；（2） 在博里叶空间从极坐标向直角坐标插值；（3） 进行2-D博里叶反变换以得到重建图像，见式（9.2.11）。这里由于第三步需要用到2-D变换，所以不能根据所获得的部分投影数据重新建图像，必须在获得全部投影数据后再重建图像。4.模型重建为了检验重建公式的正确性和把握重建算法中各个参数对重建效果的影响，人们常设计并结合处各种模型（phantom,也称幻影）图像进行实验。一副常用的实验图像是Shepp-Logan头部模型图Shepp 1974。图9.2.3给出它的一副改进结果图（尺寸，256级灰度，其中各部分的参数见表9.2.1Toft 1996。原Shepp-Logan图的对比度比较小，图中小的椭圆看不大清楚，改进图调整了各个椭圆的密度，使得到各部分的对比度有所增强，视觉效果更好了一些。表9.2.1 改进的Shepp-Logan头部模型图的参数例9.2.1 博里叶反变换重建示例在实际重建中，需要在许多方向上获得足够多的投影以恢复空间图像。图9.2.4和图9.2.5给出借助图9.2.3的模型图像和博里叶反变换重建方法得到的一组例子。从图9.2.4（a）到图9.2.4(f)依次为对图9.2.3沿圆周等角度分别进行4次投影，8次投影，16次投影，32次投影，64次投影和90次投影得到的2-D频率空间图像。从图9.2.5（a）到图9.2.5（f）依次为图9.2.4（a）到图9.2.4(f)对应的重建结果图像。由图9.2.4可见频率空间分布随投影次数增加逐渐向中心聚集的情况。由图9.2.5可看出重建图像的质量（包括受到不均匀和模糊两方面的影响）随投影次数增加而改变的情况。图9.2.3 博里叶反变换重建得到的结果9.3 逆投影重建逆投影的原理是将从各个方向得到的投影逆向返回到该方向的各个位置，如果对多个投影方向中的每个方向都进行这样的逆投影，就有可能建立平面上的一个分布。逆投影重建中最典型的是卷积逆投影重建。9.3.1 卷积逆投影重建1. 连续公式推导卷积投影重建法也是一种变换重建方法，也可根据傅里叶变换投影定理推出。在极坐标系中取式（9.2.4）的反变换：(9.3.1)如将用式（9.2.2）代入就得到由重建的公式。实际应用时如同在傅里叶反变换重建法中一样也需要在傅里叶空间引入一个窗。根据采样定理，只能在有限带宽| R | 1/(2s)的情况下对进行估计。如果定义：(9.3.2)则代入式(9.3.1)并交换对s和R的积分次序，可得到：(9.3.3)也可将上式分解成以下两个顺序的操作来完成：(9.3.4)(9.3.5)其中是在角方向的投影与h（s）的卷积，可称为在角方向上卷积了的投影，而h（s）称为卷积函数。式(9.3.4)代表的过程是个卷积过程，而式(9.3.5)代表的过程称为逆投影。因为中的参数是一条以角通过（x,y）点的射线的参数，所以是所有与过的射线所对应的卷积后投影的积分。2. 离散计算实际中式（9.3.5）代表的逆投影过程可用下式近似计算：(9.3.6)对每个，需要对KL个计算。因为K和L一般都很大，所以如果直接计算其工作量是很大的。有一种实用的方法是对计算，然后根据M个值以插值方法获得KL个值。这样式(9.3.4)的卷积的离散域中由两步操作完成：先是一个离散卷积，其结果用表示；然后是一次插值，其结果用表示。它们分别由下面两式给出：(9.3.7)(9.3.8)其中是插值函数。例9.3.1 卷积逆投影重建示例卷积可看作一种滤波手段。卷积逆投影相当于对数据先滤波再将结果逆投影回来，这样可使模糊得到校正。图9.3.1给出借助图9.2.3的模型图像用卷积逆投影重建方法得到的一组例子。从图9.3.1（a）到图9.3.1（f）依次为对图9.2.3沿圆周等角度分别进行4次投影，16次投影，32次投影，64次投影和90次投影得到的重建结果。由这些图可看出当投影次数比较少时，重建图像中沿投影方向有很明显的亮线，这是逆投影造成的结果。不过这里对每个投影得到的数据借助1-D的卷积来滤波，再将这些数据像采集那样扩散回图像，并不需要像傅里叶方法那样存储复频率空间图。图9.3.1 卷积逆投影重建中投影数量的影响3. 扇束投影重建在实际应用中常需要尽量缩短投影的时间，这样可以减少由于物体在投影期间的运动而造成的图像失真以及对患者的伤害。一种有效的方法是使用一个发射器和一组接受器，这样就可同时获得多条投影线。这种投影方式称为扇束投影方式，常用的主要有两种测量几何类型，分别见图9.3.2（a）和（b）。图9.3.2（a）对应发射器和接受器一起绕被检测物体旋转的扇束系统，其中检测器和X射线管装在同一个运动架上绕物体旋转。图9.3.2（b）对应发射器绕被检测物体旋转而接受器固定的扇束系统，其中检测器为一个完整的固定环，运动的只有X射线管（在电子射线系统中，运动的是聚焦点）。要在扇束投影的情况下进行重建，可通过将中心投影转化为平行投影再用平行投影重建技术进行重建。下面讨论如何调整前面推导出的平行投影重建公式以用于这里扇束投影的情况。考虑接受器在一段弧上等角度间隔排列的情况（见图9.3.2（c）。前面讨论中用所指定的一条射线可看做是一组用指定的射线中一条，其中是该源射线与中心射线的离散角，是源与原点连线和Y轴的夹角，它确定了源的方向。线积分现在记为（对，D是源到原点的距离）。这里假设源处在物体的外部，所以对所有的，有0,这里为与目标相切的射线与中心射线的夹角。由图9.3.2（c）可得下面关系：（9.3.9）另外设U是从源到需重建点P的距离（P的位置可用r和表示），是从源到P的连线与中心发射线的夹角，则有：（9.3.10）（9.3.11）这样可以得到：（9.3.12）利用式（9.3.9）到式（9.3.12），可将式（9.3.3）写为：（9.3.13）将换成：（9.3.14）将式（9.3.2）代入，并把对积分的上下限换为和0,则上式可仿照从式（9.3.3）到式（9.3.4）和式（9.3.5）的分解转化成两步：（9.3.15）（9.3.16）4. 傅里叶反变换重建法和卷积逆投影法的比较傅里叶反变换重建法和卷积逆投影重建法都基于傅里叶变换投影定理。不同的是，在推导傅里叶反变换重建公式时，2-D傅里叶反变换是用极坐标表示的。尽管看起来它们源出一辙，但傅里叶反变换重建法实际中较少应用，而卷积逆投影法则大量实用。其两个主要原因如下。（1）卷积逆投影的基本算法很容易用软件和硬件实现，而且在数据质量高的情况下可重建出清晰的图像。而傅里叶反变换重建法由于需要2-D插值，所以实现不易，且重建的图像质量很差。但是傅里叶反变换重建法所需要的计算量比较小，所以当数据量和图像尺寸大时比较有吸引力。在射电天文学研究中，傅里叶反变换重建法得到广泛的应用。这是因为测量到的数据对应于目标空间分布的傅里叶变换采样点。另外在MRI中，投影的傅里中变换可直接测量到，所以可直接用傅里叶反变换重建法。（2）在平行投影时导出的卷积逆投影公式可以用不同的方法修改以适用于扇形扫描投影的情况（式（9.3.14）给出一个例子）。但在平行投影时导出的傅里叶反变换重建公式还不能在保持原有效率的条件下进行修改以适用于扇形扫描投影的情况。这时需要在投影空间利用2-D插值重新组织扇形投影数据以利用平行投影算法重建图像。9.3.2 其他逆投影重建1.逆投影滤波在逆投影滤波器中，先进行逆投影，然后进行滤波或卷积。对投影进行逆投影操作得到的结果是模糊了的图像，它是真实图像与进行2-D卷积的结果。令对投影进行逆投影操作得到的模糊图像为（参见第5章）：（9.3.17）真实图像与模糊图像由下式联系（代表2-D卷积）：（9.3.18）为推导上式，可从式（5.5.3）出发：（9.3.19）再应用逆投影得到：（9.3.20）注意在上式中1-D傅里叶反变换是在傅里叶反变换是在傅里叶空间对径向变量的操作，这表明必须要先转换为极坐标。变量q就是傅里叶空间的径向变量，。如果令的1-D傅里叶反变换是，那么：（9.3.21）如果将s映射为，得到： （9.3.22）其中，用和进行了替换，径向积分对q的正值部分进行。注意到等式右边的表达是2-D的傅里叶反变换：（9.3.23）根据卷积定理： （9.3.24）因为等号右边第1项等于等于，所以式（9.3.18）得到验证。根据式（9.3.19）,可通过取2-D傅里叶变换得到重建算法：（9.3.25）或可写成：（9.3.26）将b用替换，并取上式的2-D傅里叶反变换，得到：（9.3.27）这是用于对投影进行逆投影滤波以实现重建的基本公式。引进2-D窗口函数：（9.3.28）可把式（9.3.27）写成：（9.3.29）一旦确定了窗函数，通过2-D傅里叶反变换就可得到，再经过与逆投影进行2-D卷积就可实现重建。图9.3.3给出实现上述算法的一种方案流程，其中R代表Radon变换，F代表傅里叶变换，B代表逆投影。例9.3.2 逆投影滤波重建示例图9.3.4给出一组逆投影滤波重建的示例，其中所用的原始图像仍如图9.2.3所示。在图9.3.4中，沿圆周等角度分别进行4次投影，8次投影，16次投影，64次投影和90次投影得到的重建结果分别显示在图9.3.4（a）到图9.3.4（f）中。该结果比图9.2.5要好一些。 图9.3.4 逆投影滤波重建结果2.滤波逆投影滤波逆投影也称滤波投影的逆投影重建。该方法的基本思路是：每个投影中的衰减（由目标吸收造成）与沿每个投影线的目标结构有关，如果仅从一个投影是无法获得沿该方向各个位置的吸收的，但可把对吸收值的测量平均分布到该方向上，要是能对多个方向的每个方向都进行这样的分布，则将吸收值叠加起来就可得到反映目标结构的特征值。这种重建方法与收集足够多的投影并进行傅里叶空间重建是等价的但计算量要少得多。例9.3.3 滤波逆投影重建流程滤波逆投影重建方法可看作对Radon反变换的一种近似的计算机实现方法（将希尔伯特变换消去了）。它利用傅里叶反变换将对q函数的计算转化成对径向量s的计算。图9.3.5给出滤波投影的逆投影重建流程图，其中R代表Radon变换，F代表傅里叶变换，B代表逆投影。图9.3.6给出借助图9.2.3的模型图像用滤波逆投影重建方法得到的一组例子。从图9.3.6（a）到图9.3.6（f）依次为对图9.2.3沿圆周等角度分别进行4次投影，8次投影，16次投影，64次投影和90次投影得到的重建结果。由这些图可看出滤波逆投影与卷积逆投影的效果比较接近，但重建图像的对比度较强，尤其当投影次数比较少时更为明显。图9.3.6 图像空间的滤波逆投影重建示例需要指出，随着投影数量的增加，图像的主要结构越来越明显，但原来密度均匀的区域现在密度不再均匀而有些从中心向周围逐渐减弱。另外，原来比较清晰的区域边缘也变得有些模糊了。这里的原因主要是在滤波逆投影中从各个方向来的投影在中心区域比在周围区域叠加得比较密集。这种效果与失焦光学系统产生的效果有些类似，那里点扩散函数与频率的倒数（或与频率变换中心的距离）是成正比的。9.4 级数展开重建1. 模型重建一般是要获物体每个位置的密度（线性衰减系数）值，或者说是得到重建图像中每个像素的灰度值。重建的输入数据是沿每条射线的投影积分，在沿射线方向上每个像素位置对线性衰减系数贡献的求和值（可用实际路线的长度加权）就等于测量到的吸收数值。如图9.4.1所示，对水平射线，总吸收数值是在像素处的吸收值的总和。由图9.4.1可见，对每条射线的积分都能提供一个方程（投影已知），合起来构成一组齐次方程。方程组中未知数的数量就是图像平面中像素的数量，方程的数量就是线积分的数量，它一般也是沿每个投影方向上的检测器数量乘以投影的数量。上述问题可以看作是解一组齐次方程的问题。一般情况下，这组方程的数量很大，但由于有许多值为0（由于许多像素并不包含在特定的线积分中），所以这是一个稀疏的方程组。有许多计算方法可以解决这样一个问题。下面只介绍一个简单的，利用级数展开的方法以帮助理解这类方法的原理。这个方法一般称为代数重建方法，也有称迭代算法，优化技术等的。代数重建方法与前面变换或投影的方法不同，它从一开始就是在离散域中进行的。这可借助图9.4.2的示意图来介绍。图9.4.2 代数重建法示意图在图9.4.2中，要重建的目标放在一个直角坐标网格中，发射源和接收器都是点状的，它们之间的连线对应一条射线（设共有M条射线）。将每个像素按扫描次序排为1到N（N为网格总数）。在第j个像素中，射线吸收系数可认为是常数，第i条射线与第j个像素相交的长度，代表第j个像素沿第i条射线所做贡献的权值。如果用表示沿射线方向的总吸收的测量值，则：（9.4.1）写成矩阵形式：（9.4.2）其中y是测量矢量，x是图像矢量，MN矩阵A是投影矩阵。为了获得高质量的图像，M和N都需在105量级，所以A是个非常大的矩阵。但由于对每条射线来说，它只与很少的像素相交，因此A中常只有少于1%的元素不为0。2. 无松弛的代数重建技术下面先介绍一个简单的，一般称之为无松弛的代数重建技术。首先初始化一个图像矢量x(0)，然后用下式进行迭代： （9.4.3）上式中是一个矢量，“”表示内积。这个方法的思路是这样的：每次取1条射线，改变图象中与该线相交的象素的值，从而把当前的图象矢量x(k)更新为x(k+1)。具体运算中就是将测量值与由当前算得的投影数据的差正比于aij重新分配到各个象素上去。 这个方法可以有一个几何解释。式（9.4.3）中与第i个方程对应的所有解的集合构成一个超平面。它的法线矢量是。将x(k)投影到第i超平面上就得到。如果整个方程组有惟一解，则所有超平面有一个共同交点，迭代投影的最终结果就会到达这个点。3松弛的代数重建技术有一种松弛的代数重建技术可看作是上述方法的直接推广，它的迭代表达式如下：（9.4.4）其中，r是松弛系数，它控制了收敛的速度。R的取值一般在0到2之间，当r1时，式（9.4.4）变成式（9.4.3）；当r很小时，式（9.4.4）变成与传统的最小平方解等价。图9.4.3给出利用上述方法得到的一组结果示例。这里考虑采用扇形投影，共有25个检测器均匀分布在1/4圆周上。设图像区域为1616,每次投影均被完全覆盖。设沿投影角度为00,900和1800进行投影和重建，这样一共有75个方程，256个未知量。图9.4.3（a）给出原始像素数组，其中密度值在0，20之间。投影重建从初始条件全图为均匀密度10开始，图9.4.3（b）到图9.4.3（d）分别给出经过1次，5次和50次迭代的结果。由图9.4.3可见，图像中空的区域和内部正方形很快就重建出来了，目标的边缘随迭代次数增加而得到改善。重建的误差，特别在图像的四个角处（因为只有较少的射线通过）和内部正方形的四个角处（吸收值发生急剧变化，与平均值差距较大），是比较明显的。图9.4.3 代数重建方法中的迭代效果示例4.级数法的一些特点变换法比级数法所需的计算量要小得多，因而大多数实用系统都采用变换法。但与变换法相比，级数法有一些独特的优点：（1）由于在空域中进行重建比较容易调整以适应新的应用环境，所以级数法比较灵活，常在利用新物理原理和新数据采集方法的重建使用；（2）级数法能重建出相对较高对比度的图像（特别是对密度突变的材料）；（3）借助多次迭代可用于从较少投影（10）重建图像的工作；（4）比变换法更适合于ECT，不过吸收和放射分布都需考虑；（5）比变换法更适合于3-D重建问题（也是由于在空域中比较灵活，所以易推广）；（6）比变换法更适合于不完整投影情况，这是因为