基于距离变换分水岭算法的两维传递函数体数据识别方法

摘要

本发明公开一种基于距离变换分水岭算法的两维传递函数体数据识别方法，首先通过遍历属于用户已设置的传递函数的体素，并对其进行填充和取补理，然后依据其对应的距离变换的反值利用分水岭算法对这些体素进行分类，识别出重叠于传递函数空间的感兴趣结构。由于分水岭方法过分依赖计算所得的距离值，导致感兴趣结构的过分割造成误分的现象或者无法分割开在空间上相连的结构，因此依据体素在感兴趣结构内部的距离特点修正距离值来解决这个问题。本发明在没有增加传递函数空间新维数的前提下利用距离变换分水岭算法将每个结构识别出来，并可去除不感兴趣结构对其的遮挡。甚至当这些结构在空间上相连且灰度值相近时，也能得到令人满意的识别效果。

说明书

技术领域

本发明属于计算机图形学中的三维可视化领域，具体地，涉及一种识别体数据 中感兴趣结构的新方法，即利用距离变换的分水岭算法识别出重叠于两维传递函数 空间的感兴趣结构的方法。

背景技术

伴随着现代科学技术的不断进步，新的影像采集技术(例如：CT、MRI等)逐步在医 学上的应用，提高了医生检测疾病能力的同时，也产生了大量的需要分析的医学三维体 数据。因为直接体绘制(Direct Volume Rendering)的绘制结果是直接由整个体数据中的所 有体素所产生的，即产生的是包含在体数据中不同结构细节信息的整个体数据的图像。 所以不仅可以为医生提供感兴趣的器官的信息和其与周围组织的位置信息，还能提供更 加丰富的感兴趣器官或者病灶的内部信息，从而为医生理解疾病的发展和判断病情提供 了很好的依据。因此体绘制方法是分析和可视化医学体数据的一种非常实用的方法，并 且在临床医学上具有重要的应用价值。

由于传递函数(Transfer Function)的光学属性决定了体绘制中的最终绘制结果：颜色 值(RGB值)用于直观的区别显示不同感兴趣结构，阻光度值用于决定哪些体素显示和显 示的透明度值，而在绘制结果中每个像素的RGB值是由传递函数赋予每个体素的光学 属性值累加得到的。所以，传递函数是体绘制中的起决定性作用的步骤。为了可视化体 数据中用户感兴趣的不同组织，需要通过在传递函数空间中设置和调节不同的分类器组 件来识别出这些结构。由于传递函数决定了感兴趣结构在体绘制中的最终显示效果，所 以其一直都是体绘制研究中的热点也是一个难点，并且引起了世界各国的科研人员中越 来越多的关注和研究。

为了增加传递函数空间的识别能力，有些研究人员将更多的新属性值(例如距离、 均值、方差、尺度和形状描述子等)加入到传递函数空间中，从而构成新的高维空间。 但随着空间维数的不断增加，计算这些属性的所需时间也在不断增大，而且更不利于用 户直观的理解空间中的散点区域与感兴趣结构的对应关系。即使采用不同的降维方式 (例如平行坐标、自组织映射等)，也不利于没有相关数学背景的用户设计出合适的传递 函数。而且对于不同类型的体数据，如何选择不同的属性值从而构成合适的高维传递函 数空间更利于感兴趣结构的识别，这本身就是一个值得深入研究的问题。因此，虽然高 维传递函数空间可以增强空间的识别能力，但其如何更好的实用还需要经过深入研究。

由于两维传递函数空间不仅有较好的分类识别能力，还有直观的利于用户理解的交 互界面，她是现在流行也是实用的一类传递函数空间。现在常用的两维传递函数空间有： Kindlmann et.al提出的灰度-梯度幅值两维传递函数，et.al提出的LH空间， Haidach et.al提出的基于统计量的两维特征空间等等。由于灰度-梯度幅值两维传递函数 更加的直观和易于理解，其为现在已应用于实际的一种传递函数空间。因为在构成这些 两维传递函数空间的过程中，其中的每个像素对应着体数据中属于该属性对的体素的个 数，也就是说，即使是不同结构的体素也可能对应传递函数空间中的同一个位置(例如 图4、图6中所示)，所以空间自身的性质决定了其没有空域信息，而这个性质决定了重 叠于特征空间的结构是无法通过设计传递函数来消除的，尤其当这些感兴趣结构在空间 上连接到一起时，更是难于辨别。

发明内容

为解决不同感兴趣结构在现有两维传递函数空间中存在重叠的问题(这种重叠 问题无法通过设置传递函数来识别)，本发明提出一种基于距离变换的分水岭算法的 两维传递函数的体数据识别方法，利用基于距离变换的分水岭算法具有对空间上相连 接的结构的分割能力，本发明将其引入到传递函数的设计过程中，又因为其对计算所得 的距离值比较敏感，从而导致感兴趣结构的过分割造成误分现象或者在空间上相连的结 构无法分割开，因此本发明根据体素在感兴趣结构内部的距离特点修正距离值来解决这 个问题。

为了实现上述目的，根据通过设置的传递函数空间得到是一个完整的边界的性质， 本发明首先对新构成的二值体数据进行填充操作，用以防止结构内的小的零值区域对距 离变换造成影响并最终影响到分水岭的分割效果。然后通过修正已计算的距离值，从而 消除分水岭算法对其的过分依赖性，从而可以到期望的分割效果，还可以去除其他不感 兴趣的组织对感兴趣结构的遮挡。

本发明所述基于距离变换分水岭算法的两维传递函数体数据识别方法，步骤如下：

(1)读入体数据，根据用户需要的两维传递函数空间计算其两个属性值，并由此构 建该两维传递函数空间；

在本步骤中，根据用户需要构建的两维传递函数空间对新读入的体数据的每个 体素计算所需的两个属性值x和y。然后依据这两个属性值构成用于用户交互的两维 传递函数空间。在该空间中的每个点亮度值是由体数据中同时具有这两个属性值的体素 的个数取对数得到的。如果体数据中所有具有属性值点对(x,y)的体素的个数为 I_num(x,y)，则其在两维传递函数空间中对应点的最终亮度值为I(x,y)：

$I(x,y)=\frac{\log \left(I_{\mathrm{num}}(x,y)\right)}{I_{\max }}$

其中，I_max为对所有属性值点对的个数经对数计算之后的最大值。

(2)为得到感兴趣的结构，用户通过在该传递函数空间设置传递函数的分类器组件， 然后在体数据中对属于该分类器组件的所有体素进行查找并标记，并对该体数据进行填 充操作和对其进行取补操作，得到新构成的二值体数据，对该新构成的二值体数据计算 其距离变换(Distance Transform)，然后依据该距离值的反值行分水岭算法得到分类结果。

在本步骤中，根据步骤(1)所构建的传递函数空间，用户依据需要在上面设置传递 函数的不同形状的分类组件，在调节已设置的传递函数时，观察所绘制结果的变化情况， 在不断的进行交互调整过程中，确定得到满意的可视化结果。但由于两维传递函数空间 在构成时，只是统计了具有相同属性对的体素的个数，所以其本身就不具有空间信息。 可是在很多情况下，不同结构可能具有相同的属性，即不同的感兴趣结构可能在传递函 数空间中占有相似的区域，这样无论在传递函数空间怎么设置传递函数都是无法识别出 这些结构的(如图4中的通过在传递函数空间设置分类器后三个不同感兴趣结构同时显 示出来，图6中脑部的不同结构也同时显示出来)。所以可以依据这些结构在空间上的 连通性，采用分水岭分割的方法将其识别出来。

根据用户已设置的分类器对体数据中的体素进行标记(属于该分类器的为1，其他为 0)。同时为了防止标记之后的体数据中结构内的小的零值区域对距离变换造成影响从而 影响到分水岭的最终分割效果，对新构成的二值体数据采用6邻域(图2所示)的方式进 行填充，并对填充之后的体数据取补，防止出现不同结构位于同一个集水盆区域。在对 新的二值体数据进行距离变换时，是计算体数据中的零值体素所在的位置到非零值体素 所在位置的最短距离，此处的距离是采用26邻域(图3所示)的欧式距离。例如两个体素 所在的位置分别为(x₁,y₁,z₁)和(x₂,y₂,z₂)，则它们的欧式距离D₁₂为：

$D_{12}=\sqrt{{(x_2-x_1)}^2+{(y_2-y_1)}^2+{(z_2-z_1)}^2}$

通过对距离的局部极值计算可以很明显的知道：相对于感兴趣的结构，散乱点和其他不 感兴趣的组织所具有的局部极大值点都是很小的(图6左边的散乱结构、侧脑室结构 和眼睛)，这是比较明显的一个特性，事实上，也可以从这个特性把有些不感兴趣的 结构去掉。对经过距离变换之后的值采用取反操作，从而使位于结构内部的体素变为极 小值点，这样可以识别出这些结构。在对其进行分水岭算法时是依据其26邻域的连通 性。最终得到的新体数据中标记为0的为分水岭，其他标记数字对应不同的目标。对最 终得到的体数据，可以通过查找不同数字标记的体素来显示识别出的不同的结构。因为 最终获得的感兴趣结构在整个标记的体数据中所占的比例比较大，所以只显示含有体素 数量比较大的几个类别从而查找出感兴趣的结构(图5的三个不同感兴趣球状结构，图 7的侧脑室结构和两个眼球)。

(3)若出现过分割现象或无法分割开的在空间相连的结构，则通过修正距离值的方 法得到期望的分割效果，并识别出感兴趣的结构，去除其他不感兴趣结构对其的遮挡。

在本步骤中，由于基于距离变换的分水岭方法对所计算的距离值过分依赖，如 果直接对距离变换之后的体数据采用分水岭算法，很容易造成过分割现象(图6右边 将脑部结构分割成许多小的结构)。因此可以通过修正新体数据中的距离值来消除这些 现象。由步骤(2)中可以得知，位于感兴趣结构内部的距离值要比其他结构的值大很多， 所以可以依据这个特性来修正距离值。对于位于点的距离为：

$D_{\overrightarrow{r}}=\left(\begin{array}{cc}C& D_{\overrightarrow{r}}\ge C\\ D_{\overrightarrow{R}}& \mathrm{others}\end{array}\right)$

其中C为某一距离值，在本发明中C∈[2,5]。

在修正距离值时，通过对大于等于某一数值的距离值同时赋予该值，从而消除直接 应用分水岭算法造成的过分割或分割不完全的情况。

在去除其他不感兴趣的结构的遮挡时，可以将这些不感兴趣的结构设置为0，而只 显示感兴趣的结构。

采用本发明的方法，可以识别出重叠于传递函数空间的感兴趣的结构，即使是相互 连接且灰度值变化不大的感兴趣结构也可以很好的识别出来(图7左边侧脑室无法从其 他脑部结构中识别出来，而通过修正距离值的方法可以得到完整的侧脑室结构(右边)， 并且可以将两个眼球分别显示出来)。同时还可以去除其他不感兴趣的结构对其的遮挡： 通过查找不同数字标记的体素来只显示出的感兴趣的结构，从而消除其他组织和散点 对其的遮挡(图5每个球状结构的单独显示，图7中的侧脑室的单独显示，可以去除 包含其外面的那些遮挡结构)。

本发明利用基于距离变换的分水岭算法对重叠于两维传递函数的感兴趣结构进行 识别，即使是相互连接且灰度值变化不大的感兴趣结构也可以很好的识别出来，同时可 以解决不同的感兴趣结构之间的遮挡问题，去除不感兴趣的结构和散乱点。同时本发明 还可以消除分水岭算法所造成的过分割现象，从而可以得到完整结构的可视化结果。

与现有技术相比，本发明具有如下的有益效果：

(1)本发明在没有增加传递函数空间的维数的前提下，将空域位置信息引入到两维 传递函数中，因此没有增加计算属性值所需的时间和用户交互的复杂度；

(2)本发明根据感兴趣结构在空间相连通的性质，利用基于距离变换的分水岭算 法识别出这些重叠于传递函数空间的结构，并可以分别显示不同的感兴趣结构达到消 除不同结构之间的遮挡，从而有利于用户直观的观察体数据中的感兴趣结构；

(3)本发明可以识别出在空间上相互连接且灰度值变化不大的不同感兴趣结构， 这是现在所有传递函数所无法达到的效果；

(4)本发明可以消除直接对距离值采用分水岭算法所造成过分割现象和无法分 割开的在空间相连的结构，从而提高了分水岭算法的识别能力，进而得到更好地感兴趣 结构的显示效果；

(5)本发明的原理简单且易于理解，可以很容易的在现有的两维传递函数空间中实 现并取得满意的可视化效果。

附图说明

下面利用参照附图的方式对本发明的非限制性实施例进行说明，本发明的原理 及其优点将会更直观的体现，其中：

图1是本发明方法的流程图。

图2为体数据中6邻域示意图。

图3为体数据中26邻域示意图。

图4是模拟数据在灰度-梯度幅值传递函数空间(右边)显示，以及其中间一切片 的显示(左边)的应用实例。

图5利用两维灰度-梯度幅值传递函数空间(图4右边)绘制出的模拟体数据结构 (左边)的应用实例，将其作为本发明的一个对比实验，以及通过本发明方法得到的 主要类别的直方图(右下)和分别对应的感兴趣结构(右边)。

图6是利用两维灰度-梯度幅值传递函数空间(中间)绘制出的头部结构(左边)的 应用实例，将其作为本发明的一个对比实验，以及直接应用分水岭算法所得到的分 割结果，以其中一矢状面作为显示。

图7是利用本发明识别重叠于传递函数空间的感兴趣结构的分水岭方法绘制出 的感兴趣结构的应用实例，其中左边为C＝2的时候识别出结构，右边是C＝3识别 出来的完整侧脑室和眼睛的显示。

具体实施方式

下面结合具体实施例对本发明进行详细说明，该实施例将有助于本领域的技术人员 进一步理解本发明的原理和优点，但它不以任何形式限制本发明。在此需要指出，本领 域的技术人员在不脱离本发明原理的前提下做出的若干改进和变形，这些都属于本发明 的保护范围。

如图1所示，本实施例提供一种基于距离变换的分水岭算法的两维传递函数的体数 据识别方法。为了直观的说明本发明的实用性，本实施例将本发明的原理思路应用到现 在实用的两维灰度-梯度幅值传递函数空间中，具体实施步骤为：

(1)读入体数据，在计算每个体素的梯度幅值采用中心差值方法，位于点 的体素的梯度幅值为：

$\left|\overrightarrow{G}\left(\overrightarrow{r}\right)\right|=\sqrt{{\left(\frac{\partial \overrightarrow{G}\left(\overrightarrow{r}\right)}{\partial x}\right)}^2+{\left(\frac{\partial \overrightarrow{G}\left(\overrightarrow{r}\right)}{\partial y}\right)}^2+{\left(\frac{\partial \overrightarrow{G}\left(\overrightarrow{r}\right)}{\partial z}\right)}^2}$

其中偏导数的计算公式是：

$\frac{\partial \overrightarrow{G}\left(\overrightarrow{r}\right)}{\partial x}=\frac{I(r_{x+1},r_y,r_z)-I(r_{x-1},r_y,r_z)}{2}$

$\frac{\partial \overrightarrow{G}\left(\overrightarrow{r}\right)}{\partial y}=\frac{I(r_x,r_{y+1},r_z)-I(r_x,r_{y-1},r_z)}{2}$

$\frac{\partial \overrightarrow{G}\left(\overrightarrow{r}\right)}{\partial z}=\frac{I(r_x,r_y,r_{z+1})-I(r_x,r_y,r_{z-1})}{2}$

I(r_x,r_y,r_z)表示位于点的体素的灰度值。

然后以灰度值作为X轴属性，梯度幅值作为Y轴属性构建两维灰度-梯度幅值传 递函数空间，其每个点(x,y)的亮度值I(x,y)采用如下公式计算：

$I(x,y)=\frac{\log \left(I_{\mathrm{num}}(x,y)\right)}{I_{\max }}$

I_num(x,y)是体数据中包含具有属性值点对(x,y)的体素总个数，I_max是所有I_num经过对 数处理之后的最大值。

(2)为了得到感兴趣的结构，用户通过在步骤(1)中所构建的传递函数空间设置何时 的传递函数。但由于两维传递函数空间在构成时，只是统计了具有相同属性对的体素的 个数，即所有属于不同感兴趣结构的体素只要具有相同的属性值就会在特征空间中对应 同一个点，因此当不同感兴趣结构在传递函数空间重叠时，这些结构就无法通过设置传 递函数将其分开(图4中模拟数据的三个感兴趣结构都直接对应其传递函数空间的同一 个’拱形’结构)。

依据已经设置的分类器对体数据中的体素进行分类标记：属于的设为1，其他为0。 并对此标记的体数据利用6邻域(图2所示，中间为关注的体素，其周围有6个体素) 的方式进行填充，将填充之后的体数据进行取补操作。然后根据公式计算距离值：

$D_{12}=\sqrt{{(x_2-x_1)}^2+{(y_2-y_1)}^2+{(z_2-z_1)}^2}$

其中D₁₂表示体数据中的零值体素所在的位置(x₁,y₁,z₁)到所有非零值体素所在位置 (x₂,y₂,z₂)的最短距离，此处的距离是采用26邻域的欧式距离。对所求距离值进行取 反操作并根据体素的26邻域(图3所示，中间为关注的体素，其周围有26个体素)的 连通性对其进行分水岭算法。最终得到的体数据中标记为0的为分水岭，其他标记数字 为不同的目标，可以显示出不同感兴趣结构。

(3)通过步骤(2)的处理，用户会发现对于有些数据如果直接对所求的距离值运用分 水岭算法，会造成严重的过分割现象(图6所示，同一个结构被分割成了许多小的结构)， 而如果直接对距离值进行简单的处理可能会出现有些在空间上相连接的感兴趣结构无 法分开的现象（图7左边图所示）。因此，对于位于点的距离为：

$D_{\overrightarrow{r}}=\left(\begin{array}{cc}C& D_{\overrightarrow{r}}\ge C\\ D_{\overrightarrow{r}}& \mathrm{others}\end{array}\right)$

其中C为某一距离值，在本发明中C∈[2,5]，因为大部分感兴趣结构的内部体素距离 结构边缘的距离并非很大，这样可以取得较好的分割效果(图7右边所示)。

本实施例中的步骤(2)、(3)显示的感兴趣结构是所含体素数目较多的类别的结构。

本发明具体实施例是在VS2008平台下利用C++语言、OpenGL语言、Cg语言 和Matlab语言采用上述技术思路实现的，其中在CPU中计算属性值、构建传递函 数空间以及运行分水岭算法，而体绘制的光线投射过程是利用Cg语言在GPU中实 现的。所有实验的运行环境为：Windows 7操作系统，NVIDIA GeForce GT 630M显 卡，2GB RAM，Intel Core i3 2.20GHz。

利用一个模拟数据(256×256×256)和一个头部的MRI数据(256×256×124)作 为具体实施例来例证本发明的有效性和实用性。

图5(左边)是利用两维灰度-梯度幅值传递函数空间(图4右边)绘制出的三个球 状结构，可以发现当在传递函数空间设置完分类器之后，最终的绘制结果中同时出 现了这三个结构，并且其中有两者在空间相连而且灰度值比较接近(图4左边)。通 过本发明的基于距离变换的分水岭算法可以将其分开(右边)，并通过单独显示每个 感兴趣结构可以消除彼此之间的遮挡问题。但对于大部分真实数据，由于复杂的空 间位置关系很难简单的直接通过对距离值进行分水岭变换就可以取得满意的效果。

图6(左边)是在两维灰度-梯度幅值传递函数空间设置合适的传递函数(中间)绘 制出的脑部结构图，可以发现侧脑室结构完全被其他不感兴趣的结构所包围，几乎 无法被发现，而且两个眼球也与其他结构相混迭到一起。这些结构无法通过设置传 递函数去区分的。而如果直接通过对距离值进行分水岭变换也是无法得到感兴趣的 结构的(图6右边)，因为真实体数据结构的复杂性，会将同一个结构分割为许多小 的碎片状结构，无法给用户直观的显示和理解。而通过修正距离值的方法可以很好 的解决这个问题。当C＝2时，侧脑室结构和脑部的其他结构连接到一起，并且它们 的灰度值很接近，这样就无法识别出侧脑室结构(如图7左边所示)，但当C＝3就可以识 别出来的完整侧脑室(如图7右上所示)，通过单独将该结构单独显示，可以去除其他非 感兴趣结构的遮挡(图6左边)，从而利于用户直观的观察整个侧脑室结构，还可以识别 出与其他组织混迭到一起的两个眼球结构。

上述实施例的基于距离变换的分水岭算法的两维传递函数的体数据识别方法，可以 识别出在重叠于传递函数空间的不同感兴趣结构，无论是其是否在空间上相接，即使具 有相似灰度值的相连结构，本发明的方法也可以很好的将其识别出来。同时消除不同结 构之间的遮挡问题，并丢弃不感兴趣的散乱点和结构，本发明对分水岭算法的过分割现 象采用修正距离值的方法取得较好的效果。

综上所述，本发明对于在传递函数空间重叠的不同感兴趣结构，在没有增加空间维 数的前提下利用基于距离变换的分水岭算法将每个结构识别出来，即使是具有相似灰度 值且在空间上也相互连接的结构，也可以得到令人满意的效果；本发明利用修正距离值 的方法对由于分水岭算法的过分割问题造成的误分现象进行了修正，从而得到完整的感 兴趣结构。本发明不止原理简单，利于理解，而且更容易在现有的两维传递函数空间实 现，得到满意的可视化效果。

以上是对本发明的具体实施例进行了详细解释说明。在此需要指出的是，本发 明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各 种变形或修改都属于本发明的保护范围。