数字图像处理第6章图像重建

本章重点：图像重建的基本概念傅立叶反投影重建方法卷积法重建方法代数重建方法第6章图像重建

6.1图像重建概述6.2图像重建的方法6.3傅里叶反投影重建6.4卷积法重建6.5代数重建方法6.6重建的优化问题6.7重建图像的显示6.8小结第6章图像重建

6.1图像重建概述

图像重建是指根据对物体的探测获取的数据来重新建立图像。用于重建图像的数据一般是分时、分步取得的。图像重建是图像处理中一个重要研究分支，其重要意义在于获取被检测物体内部结构的图像而不对物体造成任何物理上的损伤。它在各个不同的应用领域中显示出独特的重要性。例如：医疗放射学、核医学、电子显微、无线电雷达天文学、光显微和全息成像学及理论视觉等领域都多有应用。6.1图像重建概述根据成像光源的获取方式和成像机理的不同，可以将图像重建分为三种不同的检测模型:透射模型、发射模型和反射模型。如何从这些模型获取的信息和数据出发，建立被检测物体的内部结构图像是本章中要探讨的主要问题。图6-1图像重建采用的透射、反射、放射三种模型

6.1图像重建概述投影重建一般是指从一个物体对多条射线的透射投影重建二维图像的过程。当射线穿过物体时在检测器上得到的值就叫做射线的投影。根据投影可以初步了解组织对射线的吸收强度，但是不可能准确地判断物体内密度分布情况。因此，投影重建时需要一系列投影才能重建二维图像。6.1图像重建概述图6-2表示等强度的射线透过不同密度分布时的情况。图6-2等强度射线穿透不同组织的情况

6.1图像重建概述根据被用于图像重建的数据获取方式不同，可以分为透射断层成像、发射断层成像和反射断层成像。发射断层成像系统，一般是将具有放射性的离子(放射元素)注人物体内部，从物体外检测其经过物体吸收之后的放射量。反射断层成像系统，是利用射线入射到物体上，检测经物体散射(反射)后的信号来重建的。

如果从研究的图像维数来分，图像重建可以是针对一系列沿直线投影图来重建二维图像，也可以是由一系列二维图像重建波长不同，可以分为X射线成像、超声成像、微波成像、激光共焦成像等。6.1图像重建概述图像重建技术在许多科学领域得到广泛应用，其中最为显著的是医学方面的应用。根据原始数据获取方法及重建原理的不同可分为如下几种：放射断层重建成像(EmissionComputedTomograghy,ECT)透射断层重建成像(TransmissionComputedTomography，TCT)反射断层重建成像(ReflectionComputedTomography，RCT)核磁共振重建成像(MagniticResonanceImaging，MRI)6.1图像重建概述6.1图像重建概述图像重建经多年研究已取得巨大进展，产生了许多有效的算法，如：傅立叶反投影法、卷积反投影法、代数法、迭代法等，其中以卷积反投影法运用最为广泛。近年来，由于与计算机图形学相结合，把多个二维图像合成三维图像，并加以光照模型和各种渲染技术，已能生成各种具有强烈真实感的高质量三维人工合成图像。6.2图像重建的方法在各种图像重建算法中，计算机断层扫描技术又称计算机层析（CT）占有重要的地位。计算机断层扫描技术的功能是将人体中某一薄层中的组织分布情况，通过射线对该薄层的扫描、检测器对透射信息的采集、计算机对数据的处理，并利用可视化技术在显示器或其他介质上显示出来。这项技术的重要基础是投影切片定理：即对于任何一个三维(二维)物体，它的二维(一维)投影的傅立叶变换恰好是该物体的傅立叶变换的主体部分。6.2图像重建的方法投影切片定理给出了图像在空间域上对X轴的投影与在频率域u轴的切片之间的关系。如果投影并非是对X轴进行，而是对与空间域的X轴成任意的角度θ的方向进行投影。6.2图像重建的方法根据二维傅立叶变换的相似性定理，有：6.2图像重建的方法于是傅立叶变换可以变为：可以看出，如果f(x,y)变换一个角度，则f(x,y)的频谱也将旋转同样的角度。6.2图像重建的方法因此，f(x,y)在X轴上投影的变换即为F(u,v)在u轴上的取值，结合旋转性可得f(x,y)在与X轴成θ角的直线上投影的傅立叶变换正好等于F(u,v)沿与u轴成θ角的直线上的取值。6.2图像重建的方法投影性质是利用线扩展函数进行系统辨别和计算机断层造影术的基础。6.3傅里叶反投影重建6.3.1反投影重建的基本原理6.3.2重建公式的推导6.3.3重建公式的实用化6.3傅里叶反投影重建对一个三维物体的内部信息，可以利用物体按某个方向分解成多个薄片的了解，掌握其内部信息。现实中许多物体不能被轻易切开（例如人体），为了无损伤地探测内部信息，投影切片定理给出了很好的启示。投影切片定理给出了射线沿y轴方向穿透物体薄片对X轴投影的傅立叶变换与物体薄片的频域函数F(u,v)沿u轴的切片相等。利用二维傅立叶变换的旋转性质可知，如果围绕物体薄片，改变θ角得到多个投影，就可以获得该物体薄片在频域上相应各个方向的频谱切片，从而了解到该薄片的整个频谱。通过傅立叶反变换就能得到物体薄片在空间域中的图像。

6.3.1反投影重建的基本原理傅立叶反投影重建方法是最简单的一种变换重建方法。最早于1974年由Shepp和Logan提出，该方法是建立在“投影切片定理”这一理论基础之上的。傅立叶反投影重建方法应首先建立好以连续实函数形式给出的数学模型，然后利用反变换公式求解未知量，并适当调节反变换公式以适应在离散、有噪声干扰等条件下的应用需求。6.3.2重建公式的推导图6-5与图6-6给出了二维函数投影在空间域中和变换域中的旋转对应关系。6.3.2重建公式的推导由图可见，从（x,y）坐标系变换到（s,t）坐标系的旋转坐标变换公式为：若射线源的射线径向穿过被检测横截面。并向与X轴成θ角的S轴方向投影，则投影函数6.3.2重建公式的推导则投影函数的傅立叶变换可以写成在θ角不同的各个方向上获得空间域上的投影数据，根据投影切片定理在变换域上得到对应的切片数据。然后利用下式进行傅立叶反变换：

6.3.2重建公式的推导如果需要重建三维实体图像，很容易推广到三维：6.3.3重建公式的实用化

要使上述的反变换方法付诸实践，需要将连续函数形式转为离散形式。为方便后边的讨论，改写公式(6-7)为：其中ρ为在空间域的投影轴，θ为投影轴与原坐标x轴所成的夹角。由投影切片定理可知空间域投影轴从x变到ρ，转过θ角度，在变换域中切片也从u轴转到R轴，也旋转了θ角度。投影函数的傅立叶变换由下式给出：6.3.3重建公式的实用化

6.3.3重建公式的实用化

式(6-18)

可写成极坐标（R,θ）形式：6.3.3重建公式的实用化

用傅立叶变换方法重建图像的步骤主要分为三步：（1）首先根据(6-16)式计算出N个θ方向上投影几何的傅立叶变换，即求出多个不同θ值的F(R,θ)

；（2）利用(6-20)式，求出f’(x,y,θ);

（3）根据(6-21)式对于所有θ方向上的变换结果f’(x,y,θ)求和，便可得到重建图像的像素值。6.3.3重建公式的实用化

实现图像重建时首先需对对象进行离散采样和处理。图6-7实际对象的空间采样和投影情况6.3.3重建公式的实用化

将前面的连续方程(6-16)式改写为离散形式:6.3.3重建公式的实用化

6.3.3重建公式的实用化

6.3.3重建公式的实用化

这样便求得了数字计算的全部公式。实际执行时仍分三步：（1）据投影数据，对极轴变量m取一维傅立叶变换，即计算(6-25)式，并据(6-28)式求出（2）由(6-29)式或(6-30)式计算出的傅立叶反变换f’(pε,qε,n△θ)（3）由(6-31)式计算出重建图像f(pε,qε)6.4卷积法重建在(6-19)式中，当用FFT计算投影数据的傅立叶变换时，投影数据总被有限截断。当ρ的采样间隔为d时，在变换域R的变化范围为从-1/2d到1/2d，于是投影反变换重建公式可以近似写成：

6.4卷积法重建

由前面的(6-20)和(6-16)式可知6.4卷积法重建由式(6-35)

可以看出，要对已经得到的投影数据实现图像重建，则可以采取两个步骤：首先将投影数据g(ρ,θ)先和脉冲响应为式(6-34)

的滤波器进行卷积，然后由式(6-21)

对不同旋转角θ求和，就能实现图像重建。这就是采用卷积法进行图像重建的基本思路和方法。式(6-34)

表现出来的恰好是频率响应为｜R｜的滤波器，通常称之为ρ滤波器。6.5代数重建方法代数重建技术就是事先对未知图像的各像素给予一个初始估值，然后利用这些假设数据去计算各射束穿过对象时可能得到的投影值（射影和），再用它们和实测投影值进行比较，根据差异获得一个修正值，利用这些修正值，修正各对应射线穿过的诸像素值。如此反复迭代，直到计算值和实测值接近到要求的精确度为止。

6.5代数重建方法

具体实施步骤：（l）对于未知图像各像素均给予一个假定的初始值，从而得到一组初始计算图像；（2）根据假设图像，计算对应各射线穿过时，应得到的各个相应投影值z1*,z2*,……zn*；（3）将计算值z1*,z2*,……zn*和对应的实测值z1,z2,……zn进行比较，然后取对应差值zi

－zi*作为修正值；（4）用每条射线之修正值修正和该射线相交之诸像素值；（5）用修正后的象素值重复l～4各步，直到计算值和实测值之差，即修正值小到所期望的值为止．只要所测得的射线投影值z1，z2，……zn组成一个独立的集合，那么代数重建便将收敛于唯一解。6.6重建的优化问题考虑无噪声情况，把图像表达成按行堆积的向量：（6-36）式中，fi表示第i行向量，i＝1，2，…，N，所以f为N2×1维之列向量，它在转角为θk射线方向的投影为

（6-37）即假设也有N个分量，对于p个θ转角，可得到一个按每个角度投影堆积而成的向量（6-38）6.6重建的优化问题设投影值是图像像素之线性组合，即

区域Di（i=1,2,…,pN）表示对投影值gk,l有贡献的所有图像像素所决定的区域。于是式（6-39）又可表示成：

式中6.6重建的优化问题当p=N时，方程（6-40）有唯一解，当p<N时，方程（6-40）是超定的或者是过定的，它有不定解或无穷多种解。当P>N时，方程（6-40）是欠定的，此时方程组没有解，该方程组称为不相容的（或矛盾的）线性方程组。通常p<N，即方程组一般是超定的或是多解的，因此，涉及如何解相容线性方程组的问题。希望从这些解中找出一个满足一定优化条件的最优解，从而使重建后图像最接近原图像f。6.6重建的优化问题一种优化准则可以选为：重建图像像素和其邻域象素灰度相近。这是一种平滑准则，即准则函数可以选择

式中，C为8邻点平滑矩阵或基于二阶导数或四阶导数之平滑矩阵，它是一个半正定阵。因此有，当各点象素值相同时，它等于零。所以，准则函数（6-41）式未必能得到唯一解。如果再加上能量条件，即选择

作为约束条件，便可得到重建图像的唯一解。6.6重建的优化问题为取（6-40）式条件下的（6-42）式的最小值，使用拉格朗日乘数法，做辅助函数：6.6重建的优化问题两边乘B得：于是式中，#号表示伪逆（对应P＜n的情况，当P＝n时即为逆）。将其代入（6-45）式，最后得到：式中6.6重建的优化问题式（6-47）给出了无噪声情况下的最优化重建。在有噪声情况下，获得的投影向量中还会有噪声，即投影向量为：

重建时的约束条件就是使这一噪声尽可能小，此时，辅助函数可以改写成：

6.6重建的优化问题极小化上式，便可求得最优重建图像f，令：于是式中，为正定非奇异阵，所以逆存在。式（6-49）给出了有噪声情况下原图像的最优化重建。6.7重建图像的显示6.7.1三维图像重建的体绘制6.7.2三维图像重建的面绘制

6.7.1三维图像重建的体绘制三维图像重建类似于三维图像显示技术，属于“三维空间数据可视化”的问题。在三维空间的规则数据场中，只要给出三维网格某一角点的空间位置和某一数据点的序号，即可根据网格间距所对应的距离求出该点的空间位置。(a)规则三维空间网格与体元(b)头部切片图像图6-10规则三维空间网格、体元及系列头部切片图像6.7.1三维图像重建的体绘制两类不同的三维空间数据场可视化方法：第一类方法称为面绘制技术首先由三维空间数据场构造中间几何图元，如小三角形、小曲面等，然后再由传统的计算机图形学技术实现面绘制，加上光照模型，阴影处理，使得重建的三维图像极具真实感。第二类方法称为体绘制或直接体绘制方法是不需要构造中间几何图元，而是直接由三维数据场产生屏幕上的二维图像，由体元进行绘制。体绘制可以在空间域上进行，也可以在频率域上进行。6.7.1三维图像重建的体绘制直接体绘制的实质是将已采集到计算机上的三维离散数据场重新采样，然后按照一定的规则转换为图像显示缓存中的二维离散信号。三维离散数据场或者三维图像数据集合称为物体空间，显示图像的屏幕称为图像空间。在各种体绘制方法中最常见的有两种：图像空间到物体空间的方法物体空间到图像空间的方法6.7.1三维图像重建的体绘制

1、图像空间到物体空间的方法：确定新采样点：由显示图像屏幕上的每个像素点位置向物体空间发出射线，该射线与物体空间相交，交点即为物体空间上新的采样点。求新采样点的值：需选择适当的重构元素，对离散的三维数据场进行卷积运算，重构连续的原始信号。并根据重采样的奈奎斯特频率极限对三维连续数据场进行重新采样，得出重采样点的灰度值。在重采样的基础上，进行图像合成，即计算全部采样点对屏幕像素的贡献，得到屏幕上每一个像素点的光强度值。该算法又称为图像空间扫描的体绘制方法。6.7.1三维图像重建的体绘制6.7.1三维图像重建的体绘制2、从物体空间到图像空间的方法：对物体空间的三维空间网格点上的数据，逐层、逐行、逐个地计算它对屏幕像素的贡献，并加以合成，形成最后的图像。关键问题：如何计算物体空间三维网格的数据对二维屏幕像素的贡献？即确定物体空间每一个采样点对屏幕上多个点的贡献。6.7.1三维图像重建的体绘制6.7.1三维图像重建的体绘制实现光线投射体绘制方法是一个三维离散数据场的重采样和图像合成的过程，有以下几个步骤：

(1)对物体空间的三维离散数据场进行处理，包括去除各层图像的噪声、冗余数据。根据数据值的不同进行简单组织分类，并给不同属性的数据赋予不同的颜色和不透明度值。

(2)从图像空间显示屏幕上的每一个像素点根据设定的观察方向发出一条投射光线。这条射线穿过三维数据场，沿着这条射线选择k个等距离的采样点，当采样点刚好落在原来物体空间的网格点上，就用原来的灰度值来代替；当采样点不落在网格点时，就用该点邻近8个网格点的值进行插值得出新采样点的不透明度值和颜色值。6.7.1三维图像重建的体绘制(3)为了增强三维逼真效果，突出显示不同组织的边界面，可以采样表面并进行明暗计算。(4)计算每条射线对屏幕像素点的贡献，即沿每条射线方向从前向后或从后向前根据各采样点的颜色值和不透明值按照一定的规则加以合成，得出屏幕像素点的最终灰度值或颜色值。

将屏幕上各像素点的颜色值都计算出来后，就形成了一幅图像。6.7.2三维图像重建的面绘制三维图像重建的面绘制技术与AutoCAD或动画制作软件3DS等三维应用软件的差别：在AutoCAD、3DS等三维制作系统中，用户只需要组合基本图元，调整光源位置、强度参数，物体表面的光滑度，材料吸收光系数等参数就可以把三维物体逼真地在二维屏幕上显示出来；面绘制技术在构造出中间图元之后，再由传统的计算机图形学技术实现面绘制，可以利用现有的图形软硬件实现绘制功能，速度比体绘制要快。面绘制算法中比较典型有：立方体步进法，简称为MC法(MarchingCubes)，四面体步进法又称为MT法(MarchingTetrahedra)等。6.7.2三维图像重建的面绘制立方体步进法的基本思想：将三维数据网格分成许多体元根据物体表面特征的信息或等值面的信息给出物体等值面的相关参数值逐个测试体元的8个角点是否在等值面上，通过线性插值得出体元中的哪些点在等值面上，连接这些点即可得到用三角形成多边形来代替立方体把所有这些三角形或多边形连接起来便得到该三维数据场的三维表面信息。6.7.2三维图像重建的面绘制二维情况下MC的基本思想为了简单起见，我们先用二维数据网格来解释。下图最左边图中的二维网格数据的正方形类似于三维数据场的体元，实心圆的轮廓就是准备要近似用直线来模拟的边界信息。图6-13二维网格中数据以直线代替曲线边界6.7.2三维图像重建的面绘制从左往右、从上往下逐行逐个扫描正方形，对每个正方形的4个顶点都计算是否在实心圆内，从而得出这些顶点的状态，对于小正方形的每一条边来说，其两端的端点状态总是如下几种可能：图6-14正方形任一边两个端点的状态6.7.2三维图像重建的面绘制只有出现上图(c)或(d)情况下，才可以确定此圆的边界经过这条边的某个位置，用星点表示。如此扫描所有的小正方形之后，便可得到与网格线相交的点。连接所有的点，即可得到该圆的轮廓信息。为了使得轮廓更加光滑，一方面减少网格尺寸的大小，另一方面在网格大小给定的情况下，不是直线连接这些交点，而是非线性插值得到更多的中间点，再进行连接得到轮廓信息。6.7.2三维图像重建的面绘制三维数据场中，确定的是物体表面轮廓信息。需以下几个步骤：

1.确定包含表面轮廓面的体元

在三维数据场中，等值面与体元相交的是三角形面或多边形面，体元中的8个角点状态要复杂多了。为了描述体元与等值面相交的状况，对每个角点赋予一个函数值，代表角点的某种状态，如0代表该角点位于等值面之外，该角点所对应的特征函数值小于等值面参数C0；1代表该角点位于等值面之内，该角点所对应的特征函数值大于或等于等值面参数C0。6.7.2三维图像重建的面绘制如果某体元中的一条边的一个端点在等值面之内，而另一个端点在等值面之外，那么，这条边必定与所求等值面相交。如果一个体元的8个角点位置用A，B，…．F表示，其8个角点对应的状态为1，0，0，0，0，0，0，0，则等值面必定与AB、AD、AE相交，如下图所示：图6-15角点的状态及其体元与等值面的相交6.7.2三维图像重建的面绘制由于每个体元有8个角点，每个角点有0，1两种可能状态。因此每个体元的8个角点共有28＝256个不同的状态。要从256种不同的状态中求出每个体元的等值面是一件非常困难的事。可以利用不同的对称性： （1）绕3条坐标轴的任一轴旋转。 （2）沿3条坐标轴的任一轴镜向反转。 考虑这些对称性之后，就可用15种基本组合反映一个体元中8个角点可能存在的全部256个状态。从8个角点的状态就可以判断出所求等值面将与哪些体元相交，穿过体元的等值面的形状等信息。6.7.2三维图像重建的面绘制2.等值面与体元相交的边界交点当三维离散数据场的密度很高、体元很小时，可以假定等值面与体元边界的交点用各角点函