利用三维欧氏距离变换实现ct图像中软组织显示的方法

专利名称：利用三维欧氏距离变换实现ct图像中软组织显示的方法
技术领域：
本发明涉及模式识别与信息处理技术，特别涉及利用三维欧氏距离变换实现医学图像中软组织显示方法。
本世纪七十年代，世界第一台X射线CT实验样机诞生于英国，它是由工程师Housefield设计制造的。近三十年来，CT技术日趋成熟，至今已发展出了五代产品，尤其是螺旋CT的诞生，不仅减少了扫描时间，还使图像质量大幅提高。其在医学领域的应用更是开创了无创伤诊断的历史新纪元，取得了巨大的成功。然而现有的影像技术通常是获得人体某一断层的影像数据，然后医生通过观察胶片或屏幕进行诊断。但是无论是胶片还是屏幕显示，医务人员所观察到的仍然是二维图像，而且由于成像设备的局限性、信号源的衰减和局部体效应多因素导致伪影和噪声生成，影响了图像质量，因此诊断结果的准确性在很大程度上取决于医生的临床经验和专业知识。据国内一家权威媒体报导，国内目前一些医院的误诊率竟高达70％，这种状况令人堪忧，解决方案之一就是借助影像医学图像处理与分析系统。影像医学图像处理与分析也是近几年兴起的一门新兴学科，如今方兴未艾。
随着计算机及其相关技术的迅速发展以及图形图像技术的日益成熟，借助图形、图像技术的有力手段，影像医学图像质量和显示方法可以得到极大的改善，从而使得诊断水平可以借助于图像处理与分析手段来得到极大的改善，这不仅可以利用现有的医学影像设备极大地提高医学临床诊断水平，而且能为医学培训、医学研究与教学、计算机辅助临床外科手术等提供电子实现手段，为医学的研究与发展提供坚实的基础，具有无可估量的价值。
医学图像的三维可视化是该学科的核心组成，它利用一系列的二维切片图像进行三维模型重建和显示，是进行定量分析的前提。三维可视化实现中有两种绘制技术表面绘制和直接体绘制。面绘制最大的特点是需要先对二维数据场进感兴趣对象进行分割和三维重建，生成对象边缘等值的曲面表示，再采用光照模型绘制图像。而体绘制是将三维体数据中的体素看成一个个半透明物质，并分类赋予其一定的颜色和阻光度，由光线穿过整个数据场，进行颜色合成，得到最终的绘制结果。
CT在骨骼的显示上取得了优于其他任何设备的效果；然而受到一些技术上、原理上的限制，目前的CT对软组织显示的能力远远不如MRI。众所周知，MRI设备昂贵，并有许多检查禁忌证。如果能增强软组织的显示就有可能替代MRI。患者可以不必既作CT检查，又作MRI检查。医生也可以从同一次检查的数据中，既看到骨骼信息又看到软组织信息。这将是件非常有临床实际意义的工作。
我们知道，皮肤在人体的表面，它由表皮、真皮以及皮下组织三部分构成。其中表皮是由形状不同、大小不一的鳞状上皮细胞由基底发育而成。真皮在表皮之下，含有丰富的血管和神经末梢。皮下组织位于真皮下部。由结缔组织纤维束和大量脂肪细胞所构成，故又称作皮下脂肪层，纤维束中包含有血管、淋巴、神经等。由于皮肤、脂肪组织的密度与肌肉的软组织的密度差别较小，在医学CT图像中各种软组织灰度差别不大，因此采用传统的分割方法如阈值法不易区分出肌肉、神经血管等不同的组织，不易重建出这些组织边缘的等值面。因此，CT图像中软组织的显示通常采用体绘制的方法。但是，在原有体数据中内部肌肉、神经血管受表皮及皮下脂肪的阻挡，也无法清晰显示。要进行软组织显示，就需要在原始数据基础上去除一定深度的表面皮肤和皮下脂肪组织对应的体素，再在新数据场中进行体绘制从而显示出皮下组织的血管和肌肉等软组织。
为实现上述目的，利用三维欧氏距离变换实现CT图像中软组织显示的方法包括(1)提取步骤，从医学图像中提取物体轮廓；(2)距离变化步骤，使用三维欧氏距离快速变换生成距离图像；(3)处理步骤(剥皮处理)，对三维数据处理生成满足一定深度和灰度的新的三维数据；(4)显示步骤，对新生成的三维数据进行分类和体绘制。
本发明的利用三维欧氏距离变换实现CT图像中软组织显示方法，能够达到对皮肤外轮廓的准确识别和软组织深度的精确定位，可清晰地再现皮下血管、肌肉与骨骼的空间解剖关系，在医学领域具有重要的应用价值。并且由于采用了自己设计的三维欧式变换算法和基于体绘制的三维数据场多表面显示方法，运算速度快，能满足医生实时交互的要求。因此，此方法具有高可信度、可应用性和可采纳性。
图3是三维欧氏距离快速变换法示意图；图4是基于体绘制的三维数据场多表面显示方法示意图；图5是灰度梯度加权函数；其中，f(v)表示体元的灰度值，w1(v)为灰度加权函数；g(v)表示体元中心处的梯度值，w2(v)为梯度加权函数。
图6是软组织显示效果；其中，实验数据为人体头部CT图像，共58个切片，间距1.5mm，分辨率为512×512×16位，实验结果如下图所示

图1)是头部的一张原始CT图像，图2)是图1)经距离变换得到的距离图像，图3)是深度为9.76mm灰度为0～4000的处理后头部CT图像，图4)是原始图像的三维重建图，图5)是用图3)的新图像三维重建头部图像，它是剥除了表面皮肤后的，其中颈部的血管和面颊的肌肉组织清晰可见，图6)是对图1)数据采用腐蚀膨胀方法处理后的重建图像。通过对比可以明显看出，采用距离变换处理后的图像对血管和软组织有更清晰的重建再现。发明的详细描述下面结合附图详细描述本发明的软组织显示方法。作为一种具体实现方案，由四个步骤的图像处理和显示过程组成，该方案的结构参见图1。这四个步骤分别是提取目标外轮廓、三维欧氏距离快速变换、剥皮处理和基于体绘制的多表面显示，下面逐一进行介绍。
步骤一提取目标外轮廓(分割)这一步的目的是为距离变换算法做预处理，将目标物体从背景中分割出来，也称为二值化的过程。由于我们最终是要将三维医学CT图像中的物体从外向内进行剥皮处理，要计算物体内各点到物体外轮廓的最近距离，必须先检测出目标的边缘。对于软组织显示，也就是检测出皮肤的外轮廓。
为了缩短运算时间，针对CT图像中背景的灰度值与目标物体有较大的差值，我们可以采用传统的阈值法和区域增长法，当然也可以采用活动轮廓、边缘检测等其他分割方法。
对于区域增长法，需要用户选择皮肤轮廓上的一个点作为种子点。阈值法的关键是阈值的选择，可以由用户选择区分背景与非背景的灰度阈值，也可用自动阈值法确定阈值。常见的自动阈值法有P-参数法，状态法，微分直方图法，判别分析法和可变阈值法。针对医学图像噪声多的特点，可以采用判别分析法。即在图像灰度值的直方图中，求得阈值t把灰度值的集合分成两组，使得两组得到最佳分离。最佳分离的标准是两组的平均值的方差和各组方差的比为最大。该方法在直方图中有两个波峰时，可作为状态法起作用；即使不存在波峰时也可求出阈值。设给定图像具有L级灰度值，阈值为k，k将图像的像素分成两组1，2。组1的像素数设为ω1(k)，平均灰度值为M1(k)，方差为σ1(k)；组2的像素数设为ω2(k)，平均灰度值为M2(k)，方差为σ2(k)。设全体像素的平均灰度值定为Mτ。则组内的方差σw2=ω1σ12+ω2σ22]]>组间的方差σB2=ω1(M1-Mτ)2+ω2(M2-Mτ)2=ω1ω2(M1-M2)2]]>最佳标准 值为最大，即 取最大值。
步骤二三维欧氏距离快速变换因为需要对医学CT图像中的软组织自外向内进行处理，即要计算皮肤内软组织各点到皮肤外轮廓的最近距离，所以我们以皮肤外轮廓之外的背景点为特征点(二值化值为0)，其余部分为非特征点(二值化值为1)。距离变换就成了求解图形中各像素到0-像素的最短距离的处理，离轮廓越远(靠中心越近)，变换后的像素特征值就取得越高。通过距离变换计算，我们可以得到距离图像，其中相对于轮廓内各点的特征值为该点到目标外轮廓的最近距离，而物体外轮廓之外的背景点对应距离图像中的值均为0。
距离是图像处理中十分重要的概念。距离的定义依据应用环境的特点，有不同的定义。距离变换是将二值文件转换成距离图像的变换。它是一个映射，此映射中的每个像素的值是在物体或背景中该点到其最近的特征像素的距离。在距离图中，像素的值是不与任何特征像素相交的最大圆半径。距离变换被广泛应用于图像处理和模式识别领域如目标细化，骨架抽取，形状的插值和匹配。
在实际计算中，常采用两种距离测度非欧氏距离和欧氏距离，前者常用的是城市街区(City Block)、切削(Chamfer)等距离，算法采用串行扫描实现距离变换，在扫描过程中传递最短距离信息，计算简单，但得到的仅是欧氏距离的一种近似值。后者采用行列交叉处理的方法，缩小搜索最近黑点的范围，算法复杂度达到O(n2)。考虑到医学图像对距离的精确测量有较高的要求，我们使用欧式距离变换。常规的三维欧式距离变换的算法复杂度为O(n6)，运算速度慢，实用性不高。为改善运算速度，减少运算量，我们在陈岐(陈岐.完全欧几里德距离变换的最优算法.计算机学报，1995，(18)8611-616)提出的二维欧氏距离变换算法基础上，设计实现了一种新的三维欧氏距离变换算法。
三维欧氏距离变换算法的基本思想是将一个n×n×n的三维二值图像分解成n个n×n的二维二值图像，我们将特征象素称为“黑点”，背景象素称为“白点”。此基本思想源于三维医学图像通常是由一组二维断层图像组成，从严格意义上讲不算是真正的三维图像，而且通常是规则数据场。我们首先采用二维欧氏距离变换算法对n个二值图像进行二维欧氏距离变换，得到每个象素在二维图像中的最近黑点。然后将每个二维图像中的象素与其它各层二维图像中的黑点进行距离比较，并通过优化方法来减少需要参与距离比较的二维图像层数和其中黑点的个数，从而找到每个象素在三维空间的最近黑点。将三维图像中每个象素的特征值赋为该象素到其最近黑点的距离，也就得到了三维欧氏距离图。该算法时间复杂度为O(n3logn)。
该三维欧氏距离快速变换算法的根本特点是采用了优化方法，减少需要参与距离比较的二维图像层数和其中黑点的个数，优化设计基于以下三个命题，下面将给出证明(参见图2)。
命题1设(i，j，k1)与(i，j，k2)为两个二维图像上具有相同行列位置的两个象素，且已知(i，j，k2)在其二维图像中的最近黑点为(a，b，k2)，则(i，j，k1)在(i，j，k2)所在二维图像中的最近黑点亦为(a，b，k2)。
证明对(i，j，k2)所在二维图像中的任意一个象素(z，w，k2)，根据已知条件知(z-i)2+(w-j)2≥(a-i)2+(b-j)2则(z-i)2+(w-j)2+(k2-k1)2≥(a-i)2+(b-j)2+(k2-k1)2即(i，j，k1)在(i，j，k2)所在二维图像中的最近黑点为(a，b，k2)。证毕。
由命题1可知，为找到一个象素(i，j，k)在另一个二维图像中的最近黑点，不必与该图像中所有黑点进行距离计算与比较，一次便可命中。为找到一个象素在三维空间中的最近黑点，本来需要将其与n个二维图像中的所有黑点做距离计算与比较，应用命题1，则只需对每个二维图像中的一个象素做一次距离计算，并作比较即可求得最近黑点。因此大大缩小了在三维空间中对最近黑点的搜索范围。
为减少参与距离计算、比较的二维图像的个数，我们给出命题2。
命题2设(i，j，k1)与(i，j，k2)为两个二维图像上具有相同行列位置的两个象素且k1＜k2，若V(i，j，k1)＝(a，b，c)，V(i，j，k2)＝(o，p，q)，其中V(x，y，z)表示求取象素(x，y，z)二维图像中的最近黑点的函数，则c≤q。
证明根据已知条件v(i，j，k1)＝(a，b，c)，V(i，j，k2)＝(o，p，q)，则(o-i)2+(p-j)2+(q-k1)2≥(a-i)2+(b-j)2+(c-k1)2(a-i)2+(b-j)2+(c-k2)2≥(o-i)2+(p-j)2+(q-k2)2，相加两个不等式，得(q-k1)2+(c-k2)2≥(c-k1)2+(q-k2)2，由已知条件k1＜k2可得c≤q。证毕。
如果已知V(i，j，k)＝(o，p，q)，由命题2可知，对其它二维图像上与(i，j，k)具有相同行列位置的象素(i，j，l)，若l＜k，则v(i，j，l)只能在第1至第q个二维图像上；若l＞k，则v(i，j，l)只能在第q至第n个二维图像上。推而广之，若n个具有相同行列位置的象素被d 等分，则这些象素的最近黑点V也被分成d个部分，求不同d等分象素的最近黑点只需到该点所处的相应部分去寻找即可。
由此，我们可以采用等分法缩小对最近黑点的搜索范围。我们首先对第n个二维图像上的象素在三维空间中寻找它们的最近黑点，根据命题1，这一步中每个象素都要做n次距离计算与比较；然后寻找第 个二维图像中的象素的最近黑点，对它的处理实际上是对n个二维图像进行二等分，应用命题2，为寻找它们的最近黑点需参与距离计算与比较的二维图像个数可大大减少；接着寻找两个等分中心图像即第 和第 个二维图像中象素的最近黑点，处理结束后n个二维图像被四等分；然后寻找四个等分中心图像即第 和 个二维图像中象素的最近黑点；依此类推，到最后寻找 个等分中心图像中的象素的最近黑点。对n个二维图像进行1，2，4，8，...， 等分，共进行log2n次等分后，每一个二维图像全部被处理过一次，则所有象素的最近黑点全被找到。
在上述过程中，等分越细，寻找等分中心二维图像上象素的最近黑点的搜索范围越小即参与距离计算与比较的二维图像个数越少。
命题3设n个二维图像已被d(d＝1，2，4，8，...， )等分，为寻找d个等分中心二维图像中相同行列位置的一组象素的最近黑点，需要参与距离计算的次数小于1.5n。
证明n个具有相同行列位置的象素被d等分，根据命题2，这些象素的最近黑点V也被分成d个部分，设其中第k部分元素个数为lk，则Σk=1dlk=n]]>，由于V在相邻两部分的边界上的元素同时属于两个部分，为第k个等分的中心象素求最近黑点需要使用V的相应部分的lk+1个元素进行距离计算并比较。对d个等分的所有中心象素求它们的最近黑点需要Σk=1d(lk+1)-1=n+d-1<1.5n]]>次距离计算。证毕。
应用命题3，求d个等分中心二维图像中所有象素的最近黑点，需要参与距离计算的次数小于1.5n3。
接下来我们给出三维完全欧氏距离变换的算法描述，参见图3步骤1对n个二维二值图像进行二维距离变换，求得每个象素(i，j，k)在其所在二维图像中的最近黑点，记为M(i，j，k)；步骤2第n个二维图像中的每个象素(i，j，n)依次与n-1个点M(i，j，k)(k＝1，2，...，n-1)进行距离计算与比较，距离最近的点即为其在三维空间中的最近黑点；步骤3设置等分数变量d初值为2；步骤4对n个二维图像进行d等分，等分平面的序号q依次为1dn,3dn,...,2d-1dn]]>，且q＜n；步骤5依次计算第q个二维图像中象素的最近黑点，q=1dn,3dn,...,2d-1dn:]]>第q个二维图像中的每个象素(i，j，q)与其最近黑点可能范围中的点M(i，j，k)依次进行距离计算与比较，距离最近的点即为其在三维空间中的最近黑点；步骤6若 ，则已扫描完n个二维图像，结束；否则在上一次等分基础上再二等分，即d＝2d，返回步骤4继续运算。
因为最近黑点是通过距离比较得到的，因此每个象素的在三维空间中最近黑点找到后，它们之间的距离也同时被求得。
根据陈岐的文献，一个n×n的二维图像其欧氏距离变换算法的时间复杂度为O(n2)，则第1步的时间复杂度为O(n3)。算法的主要计算在第5步，根据命题3，第5步需参与距离计算的次数小于1.5n3，结合第3步可知整个算法的时间复杂度为O(n3logn)。
步骤三剥皮处理我们对原始的医学CT图像结合第二步得到的距离图像进行处理。首先，用户可以通过界面交互式指定距离深度，不同的距离深度对应最后的显示效果也不一样，另外，用户也通过界面指定灰度和梯度的范围。所谓的“剥皮”处理就是将原始图像中距离特征值在指定深度范围内以及灰度和梯度在指定范围内的点去掉(即将该点的灰度值置为0)，这样就得到了一个新的三维医学图像(一组二维切片数据)。
步骤四基于体绘制的三维数据场多表面显示为了清晰地显示肌肉、神经血管等皮下软组织，我们在上一步得到的新三维图像场基础上，采用一种基于体绘制的三维数据场多表面显示方法进行三维显示，便于用户从多个角度进行立体观察，以及进行定量分析和深入处理。
传统的直接体绘制法将数据场中的体元看成一个半透明物质，并赋予其一定的颜色和阻光度，由光线穿过整个数据场，进行颜色合成。目前有三类直接体绘制方法光线投射法、投影成象法和频域变换法。光线投射法从屏幕上象素出发，发出光线穿过数据场，每条光线穿过数据场时进行采样和颜色累积，得到相关象素的颜色，直至形成最后的可视图。该方法成象质量较好，但速度慢；投影成象法沿某个投影方向将数据场中的体元逐个投影到屏幕上，屏幕上的每个象素将所获得的体元影响进行累积计算以得到它的颜色，直至形成最后的可视图。该方法的成象速度快，但难以进行光照计算，成象质量较差；频域变换法利用Fourier变换将三维的数据场空间转化为三维的频域空间，并由频域空间的二维切片获得三维数据场空间的可视图。该方法所生成的可视图不能反映空间域中颜色合成时的遮挡关系，使观察者难以判断物质分布的前后关系(王文成，吴恩华.用于体绘制的可变模板法.计算机学报，1997年7月，20(7)592-599.)。
直接体绘制法将数据场中的多种物质在一个可视图中显示，揭示它们的相互关系。但图像难免比较模糊，且由于遮挡关系，离视点较远的部分不易被观察和分析。
1988年Levoy(Levoy M.Display of surfaces from volume data.IEEEComputer Graphics and Applications，1988，8(3)29-37.)利用分类函数对三维数据场中的体元进行阻光度分配，并对物质边界面处的体元赋予高阻光度以突出显示物质表面，每一个阻光度不为0的体元对最终图像的显示都有所贡献。唐泽圣等(唐泽圣，袁骏.用图像空间为序的体绘制技术显示三维数据场.计算机学报，1994年11月，17(11)801-808.)对此方法进行了改进，以重构三维数据场代替原有的重构光亮度场，提高了可视图的质量并节约了存储空间，并用改进的方法来显示多等值面，但需要通过不断的采样来来判断等值面的位置，计算量很大。Udupa等(Udupa J K，Odhner D. Shell rendering.IEEE Computer Graphics &Application，1993，13(6)58-67.)采用投影成象法来绘制三维数据场并显示多等值面，但是对投影成象中光照效应差的问题仍没有解决，仅用体元中心的梯度值来代替整个体元的法向量。
在许多应用中，人们有时只关心不同物质间的边界部分，并由其反映整个数据场的全局状况。针对此特点，我们提出了一种基于体绘制的三维数据场多表面显示方法。为加快显示速度，我们不用分类函数对数据场中的体元进行阻光度分配(唐果，赵晓东，汪元美.模糊表面构造、阻光度映射及快速体绘制研究.电子学报，1999年1月，27(1)17-21.)，我们用灰度梯度加权提取出三维数据场中不同物质间的边界，只对这些边界上的体元进行颜色合成计算，因而可大大减少计算量。此外，还可根据显示需要对不同的边界赋予相应的阻光度，可灵活的反映不同物质边界在三维空间中的相互关系；在提高显示图像质量方面，传统的投影成象方法用体元中心的梯度值来代替整个体元的法向量，因而无法反映体元所影响的象素之间的差别。我们将边界上的体素作为不同物质的混合体，采用方向有关的三线性插值来计算视线方向与体素内等值面的交点，对体元所影响的多个象素存在着不同的交点，根据交点的法向量进行光照效应计算，从而可提高显示图像的质量；最后用投影成象法显示图像。
首先，我们采用灰度梯度加权提取出三维数据场中不同物质间的边界，根据显示的需要我们只对这些边界上的体元赋予相应的阻光度并进行光亮度合成计算，因而可大大减少计算量，提高显示的速度；我们将边界上的体元作为不同物质的混合体，采用与方向有关的三线性插值来计算视线方向与体素内等值面的交点，根据交点的法向量进行光照效应计算以提高显示图像的质量；最后用投影成象法显示最终的图像。
该方法由三个步骤组成物质边界提取，不同物质阻光度赋值，边界体元内等值面与视线交点的法向量计算，和光亮度合成计算，算法示意图参见图4。
步骤1物质边界提取在我们的方法中，我们只关心不同物质的边界面，而不考虑同一物质内部的体元对最后显示图像的贡献，因此可大大减少三维数据场的绘制时间。对物质边界的提取主要依据三维数据场的物理属性而定。对物质边界的提取属于图像分割的问题，本文不予重点讨论。我们对三维数据场进行灰度梯度加权处理，并根据要提取的物质边界数目设定相应的阈值，以判断提取不同的物质边界。加权处理的粗略函数如图5所示。
图5中f(v)表示体元的灰度值，g(v)表示体元中心处的梯度值，加权处理函数为w(v)＝w1(v)×w2(v)。如果我们将原始三维数据场表示为<V，f>，其中V为数据场中所有体元的集合，f为体元的灰度值。则经过灰度梯度加权后的三维数据场可表示为<V，w>。
步骤2不同物质阻光度赋值根据提取的物质边界数目，为不同类型的物质分配不同阻光度。一种物质的阻光度越大，透过该物质体元的光越少，位于其后的物质越不可见。相反，阻光度越小，该物质的透明度就越高，我们就可以同时看到前后不同的物质。根据希望突出显示哪种类型物质的需要，允许对阻光度的赋值进行调整。
对三维数据场，我们将其分解成二维切片数据的集合，并按切片顺序的方式来保存物质边界上的体元。对任一个切片上的边界体元，我们定义一种数据结构来描述它。
voxel[i]＝{x，y，z，tt，op}其中i表示切片的序号；(x，y，z)表示边界体素在三维数据场中的位置；tt表示物质边界的类型，如皮肤、肌肉或骨骼；op表示赋予该类型物质边界上的体元的阻光度。
需要指出的是上述物质边界具有一定的厚度。
步骤3边界体元内等值面与视线交点的法向量计算提取出不同物质边界后，我们可以采用投影成象法对边界上的体元进行快速显示。投影成象法的一个主要缺点是对每个体元难以获得有效的较准确的法向量进行光照效应计算，以突现物质分布的空间感。传统的方法多以中心差分获得体元中心点的梯度值来代替体元的法向量，这不能反映一个体元投射到屏幕上所影响的象素的变化。Webbor(Webber R.E.Ray tracing voxel based data via biquadratic local surface interpolation.TheVisual Computer，1990，6(1)8-15.)根据体元及其26个相邻体元的灰度值，用一个双二次函数在每个体元中插值出一个曲面，再依此曲面去计算法向量，这种方法能反映体元投影到屏幕上对不同象素影响的细微变化，但它的计算量太大，严重影响成象的速度。
边界上的体元是不同物质的混合体，在体元内部存在一个等值面。本发明采用一种与方向有关的三线性插值来求得视线方向与等值面的交点，并计算交点的法向量，视线穿过体元，对一个象素的影响只要计算交点处的光照效应就可以了。我们以体元的8个顶点处灰度值的平均值m作为阈值，以视线方向作为体元内插值的方向。设视线与体元的两个交点分别为p1和p2，其值为v1和v2，则这条射线与以m定义的等值面的交点为p＝(1-k)p1+kp2其中 k必须满足0≤k≤1，才能保证交点在体元内，否则我们认为该体元对射线到达的象素没有影响。计算出交点在体元内的位置后，我们根据体元的8个顶点处的梯度值做三线性插值计算出交点处的梯度值，代替该点的法向量进行光照效应计算，作为体元对其影响的象素的光亮度贡献。
步骤4交点的光照效应计算视线方向与边界体元内等值面的交点处的颜色值根据Phong光照模型求得，即C＝Ia+KdIl·(N·L)+KsIl·(N·V′)其中，Kd，Ks表示漫反射系数和镜面反射系数；Ia，Il表示环境光强度和主光源强度；N，L，V表示归一化的表面法向量、光源方向、和视线反方向；V′=V+L|V+L|]]>交点的表面法向量N可由步骤3求得。
步骤5投影平面上象素的颜色合成我们用投影成象法来快速显示边界上的体元，对任何视线方向，边界上的体元以逐个切片、逐行、逐列的方式投影到投影平面上。图像合成的方式采用从后往前法，其公式为Cout·αout＝Cnow·αnow+(1-αnow)·Cin·αinαout＝αnow+(1-αnow)·αin其中Cnow，αnow表示视线与体元内等值面交点处的颜色和阻光度，Cin，αin代表视线进入交点前的颜色和阻光度，Cout，αout表示叠加上交点处的颜色和阻光度后的相应值。
我们分配一个缓冲区IB，大小为投影平面上的象素数，初始值设为背景颜色，边界上的体元开始投影后，对每一个体元所影响的多个象素p，IB(P)更新为Cout·αout，直至所有边界体元全被处理过，则缓冲区中包含的即为最终显示图像的各个象素颜色。
经过实验表明，基于体绘制的三维数据场多表面显示方法显示图像的效果好，能够清楚的反映三维数据场中的全局状况，并且绘制的速度非常快，基本达到了实时显示的要求。
下面说明利用三维欧氏距离变换实现CT图像中软组织显示的具体实施过程。实验数据为人体头部CT图像，共58个切片，间距1.5mm，分辨率为512×512×16位。
1)首先通过数据接口读取数据。
2)点击分层显示按钮，进入软组织显示界面。
3)进入软组织分层显示界面后，如果已经有从前生成的二值文件可以选择“显示已有二值文件”；否则选择“新建二值文件”按钮。后者将打开生成二值文件界面，用户可以选择某种分割算法提取轮廓边界，分割结束后返回软组织分层显示界面。新建或打开已有的二值文件显示在原始图像下方。
4)在灰度范围和深度范围框中填入数据；按“分层显示”，系统将在距离变换后再进行剥皮处理，在右图会显示出处理后的二维图像。
5)保存新生成的一组二维切片数据。
6)打开新生成的CT图像，进入体绘制参数定义界面，可以通过分类函数将表面、骨骼、肌肉分出层次，然后采用体绘制的方法进行三维显示。软组织显示效果图参见图6。
上述实验结果与发明人对利用三维欧氏距离变换实现CT图像中软组织显示方法的理论分析结论一致，具有高可信度、可应用性和可采纳性。
权利要求
1.一种利用三维欧氏距离变换实现CT图像中软组织显示的方法包括步骤(1)提取步骤，从医学图像中提取物体轮廓；(2)距离变化步骤，使用三维欧氏距离快速变换生成距离图像；(3)处理步骤(剥皮处理)，对三维数据处理生成满足一定深度和灰度的新的三维数据；(4)显示步骤，对新生成的三维数据进行分类和体绘制。
2.按权利要求1所述的方法，其特征在于所述的提取步骤是距离变换步骤的的预处理步骤。
3.按权利要求1所述的方法，其特征在于所述的三维欧氏距离快速变换采用优化方法，所述的优化方法是减少需要参与距离比较的二维图像层数和其中“黑点”的个数。
4.按权利要求3所述的方法，其特征在于只需对每个二维图像中的一个像素作一次距离计算，通过比较后可获得最近黑点。
5.按权利要求3所述的方法，其特征在于对n个二维图像进行二等分。
6.按权利要求5所述的方法，其特征在于n个二维图像进行1，2，4，8...，n/2等分，共进行log2n次等分。
7.按权利要求6所述的方法，其特征在于为求得等分中心二维图像中所有像素的最近黑点，需要参与距离计算的次数小于1.5n3。
8.按权利要求1所述的方法，其特征在于所述的三维欧氏距离快速变换包括步骤1对n个二维二值图像进行二维距离变换，求得每个象素(i，j，k)在其所在二维图像中的最近黑点，记为M(i，j，k)；步骤2第n个二维图像中的每个象素(i，j，n)依次与n-1个点M(i，j，k)(k＝1，2，...，n-1)进行距离计算与比较，距离最近的点即为其在三维空间中的最近黑点；步骤3设置等分数变量d初值为2；步骤4对n个二维图像进行d等分，等分平面的序号q依次为1dn,3dn,...n2d-1dn]]>，且q＜n；步骤5依次计算第q个二维图像中象素的最近黑点，q=1dn,3dn,...,2d-1dn:]]>第q个二维图像中的每个象素(i，j，q)与其最近黑点可能范围中的点M(i，j，k)依次进行距离计算与比较，距离最近的点即为其在三维空间中的最近黑点；步骤6若 ，则已扫描完n个二维图像，结束；否则在上一次等分基础上再二等分，即d＝2d，返回步骤4继续运算。
9.按权利要求1所述的方法，其特征在于所述的对三维数据进行分类和体绘制包括步骤(1)物质边界提取；(2)不同物质阻光度赋值；(3)边界体元内等值面与视线焦点的法向量计算；(4)焦点的光照效应计算；(5)投影平面上像素的颜色合成。
全文摘要
一种利用三维欧氏距离变换实现CT图像中软组织显示的方法，包括(1)提取步骤，从医学图像中提取物体轮廓；(2)距离变化步骤，使用三维欧氏距离快速变换生成距离图像；(3)处理步骤(剥皮处理)，对三维数据处理生成满足一定深度和灰度的新的三维数据；(4)显示步骤，对新生成的三维数据进行分类和体绘制。本发明能够达到对皮肤外轮廓的准确识别和软组织深度的精确定位，可清晰地再现皮下血管、肌肉与骨骼的空间解剖关系，在医学领域具有重要的应用价值。并且由于采用了自己设计的三维欧式变换算法和基于体绘制的三维数据场多表面显示方法，运算速度快，能满足医生实时交互的要求。因此，此方法具有高可信度、可应用性和可采纳性。
文档编号G06T15/00GK1403057SQ0114213
公开日2003年3月19日 申请日期2001年9月13日 优先权日2001年9月13日
发明者田捷 申请人:田捷