CN103258342A - 基于灰度-3d susan算子两维直方图体可视化方法 - Google Patents

Abstract

本发明公开一种基于灰度-3D SUSAN算子两维直方图体可视化方法，步骤为：(1)读入三维体数据，首先计算出体数据中每个体素的3D SUSAN算子，然后计算出灰度-3DSUSAN算子两维传递函数空间；(2)根据对灰度-3D SUSAN算子特征空间直方图的理解，设置二维传递函数分类器组件并根据期望的绘制效果做出调整；(3)根据确定的传递函数所赋予每个体素的光学属性值绘制出最终结果。本发明首次将具有抗噪能力的SUSAN算子作为新属性引入到传递函数设计过程中，构造出两维灰度-3D SUSAN算子传递函数；原理简单，易于实现；更容易设计出合适的传递函数，得到用户满意的绘制图像。

Description

基于灰度-3D SUSAN算子两维直方图体可视化方法

技术领域

本发明属于计算机图形学中的三维可视化领域，具体地，涉及一种识别体数据中感兴趣组织的新方法，即基于灰度-3D SUSAN算子两维直方图体可视化方法。

背景技术

直接体绘制(Direct Volume Rendering)是可视化三维体数据的一种非常的重要方法。相对于传统的医生观察二维切片的方式，体绘制不止将二维图像序列以三维体绘制结果的形式直观展示出来，而且还可以为用户提供更加丰富的解剖信息，从而识别出器官或病灶的形状和大小以及与周围组织的空间位置关系，所以在医学的临床诊断上具有非常重大的实用价值。

在应用体绘制可视化三维数据的过程中，为了识别出感兴趣的组织，通常在显示时需要抑制其周围不感兴趣的组织结构和遮挡住感兴趣组织的结构，这就需要对不同的组织设置不同的阻光度传递函数；为了使用户感兴趣的器官或病灶具有更直观的表现力，从而更清楚地理解其细节信息，需要为每个组织指定合适的颜色值。这种为体数据中的每个体素设置不同光学属性(例如颜色、阻光度等)的过程就是设计传递函数(TransferFunction)的过程，由于其在直接体绘制中的重要性地位，在20世纪90年代被VolumeVisualization会议列为可视化领域里十大关键性技术之一，作为医学可视化研究中的一个热点和难点，越来越多的受到世界各国研究人员的关注。

传递函数作为直接体绘制的核心技术，按采用的属性个数的不同可以分为一维传递函数、二维传递函数和高维传递函数。在医学影像成像(例如CT、MRI、PET等)时由于噪声和部分容积效应等因素对最终测得体数据的影响，造成体素的灰度值有些扩散，甚至不同组织可能含有具有相同灰度值的体素，所以不易通过简单的灰度值这种一维传递函数识别。通过引入越来越多的特征属性虽然可以提高传递函数空间的分类能力，但由于不同属性值之间并非完全不相关的，所以对于不同的体数据选择几个特征属性以及哪几个属性构成传递函数空间才能达到最好的分类效果，现在并没有一个统一的标准；并且随着特征空间的维数增加，用户对其的理解能力和如何设计合适的传递函数，也是一个比较辣手而又无法避免的难题。

因为加入了新的分类属性值，使得二维传递函数具有比一维传递函数更强的分类能力，而且相对于高维特征空间，其更便于用户的理解和交互。SUSAN算子是Smith和Brady在《SUSAN-A New Approach to Low Level Image Processing》中提出的特征提取方法，利用统计模板中像素与核心像素灰度值的比较结果，得到USAN区域。有别于传统的基于微分的边缘检测算子，整个统计过程是一个积分过程，所以能对噪声数据的处理可能得到较好的效果。基于此，本发明设计出一种基于灰度-3D SUSAN算子的新两维传递函数，将传统的SUSAN算法扩展到三维，将其作为直接体绘制中传递函数的一个新特征属性，与灰度值构成一个新两维传递函数空间，并将其用于三维数据可视化中。

发明内容

针对现有技术中的缺陷，本发明的目的是为了使用户更好的提取体数据中感兴趣结构信息而提出一种基于灰度-3D SUSAN算子的新识别方法，尤其对受噪声影响比较严重的体数据。

本发明依据SUSAN算法在进行边缘检测时具有对局部噪声不敏感、较好的抗噪能力的特性，将传统的SUSAN算法扩展到三维，作为直接体绘制中传递函数的一个新属性值，并结合体数据本身的灰度值构成基于灰度-3D SUSAN算子的新两维特征空间，该传递函数空间可以为用户提供更好的理解体数据中结构信息的能力，利于用户依据自己的理解设计出合适的传递函数，绘制出满足其期望的体绘制结果。

本发明所述基于灰度-3D SUSAN算子两维直方图体可视化方法主要包括如下步骤：

(1)读入三维体数据，计算体数据中每个体素的3D SUSAN算子，其计算公式为：

$\mathrm{SUSAN}\left({\stackrel{\&RightArrow;}{r}}_i\right)=\left\{\begin{array}{ccc}g-\mathrm{USAN}\left({\stackrel{\&RightArrow;}{r}}_i\right)& \mathrm{if}& \mathrm{USAN}\left({\stackrel{\&RightArrow;}{r}}_i\right)<g\\ 0& & \mathrm{otherwise}\end{array}\right.,$

$\mathrm{USAN}\left({\stackrel{\&RightArrow;}{r}}_i\right)=\underset{|{\stackrel{\&RightArrow;}{r}}_j-{\stackrel{\&RightArrow;}{r}}_i|\&le;\mathrm{\&delta;}}{\mathrm{\&Sigma;}}{e}^{-{\left(\frac{I\left({\stackrel{\&RightArrow;}{r}}_j\right)-I\left({\stackrel{\&RightArrow;}{r}}_i\right)}{T}\right)}^6}$

其中

是位于点的灰度值，δ代表体素周围某一小邻域的尺度， $|{\stackrel{\&RightArrow;}{r}}_j-{\stackrel{\&RightArrow;}{r}}_i|=\sqrt{{(r_{\mathrm{jx}}-r_{\mathrm{ix}})}^2+{(r_{\mathrm{jy}}-r_{\mathrm{iy}})}^2+{(r_{\mathrm{jz}}-r_{\mathrm{iz}})}^2},$ T为灰度差阈值，g＝kNUM_max为几何阈值，k为调节参数，NUM_max为体素δ邻域内所能取得最大值，

是位于点

周围球形邻域δ内各体素的灰度值与

的相似度的和值；。

然后统计体数据中具有每个灰度-3D SUSAN算子数据对(I,SUSAN)的体素的个数IS_NUM(I,SUSAN)，并按照下面计算最终在两维直方图中的显示亮度值b_value(I,SUSAN)：

$b\_\mathrm{value}(I,\mathrm{SUSAN})=\frac{\log (\mathrm{IS}\_\mathrm{NUM}(I,\mathrm{SUSAN}))}{{\mathrm{ISLOG}}_{\max }}$

其中ISLOG_max为灰度-3D SUSAN算子数据对经对数处理之后的最大值，最后由灰度值作为新特征空间X轴的属性值，3D SUSAN算子作为Y轴属性值，构成一个新的两维传递函数空间。

(2)根据对灰度-3D SUSAN算子特征空间直方图的理解，设置二维传递函数分类器组件并根据期望的绘制效果做出调整；

为绘制出体数据中某个感兴趣的结构，用户通过对灰度-3D SUSAN算子两维直方图中区域与该感兴趣组织的对应关系，在特征空间上放置某个传递函数分类器组件，然后根据用户期望得到的绘制结果调节所设置的分类器组件，并观察由于调整组件而绘制的三维图像是否为期望的结果。分类器组件可以根据不同的需要设置不同形状，在本发明中采用的为四边形的组件，该组件有4个控制点和4个相互之间的连线组成，在调整组件时，既可以单独调节某个控制点，也可以调节两点之间的连线。

(3)在通过步骤(2)的交互操作之后，可以得到确认的传递函数，然后根据设置的传递函数赋于每个体素的光学属性值依据从前往后的累积过程绘制出三维效果图。

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

(1)本发明将具有较好的抗噪能力的3D SUSAN算子作为一个新的属性值引入到直接体绘制的传递函数设计中。

(2)本发明通过引入3D SUSAN算子这一新的属性值，构造出一个全新两维传递函数空间：灰度-3D SUSAN算子特征空间。

(3)本发明的原理比较简单，易于理解，便于在其他基于硬件和基于软件的体绘制算法中实现；

(4)本发明建立的新两维特征空间比传统的两维灰度-梯度幅值传递函数空间具有更好的表现力，提高了用户对体数据中结构信息的理解能力，利于用户设置合适的传递函数。

(5)由于SUSAN算子具有对局部噪声不敏感、较好的抗噪能力的特性，本发明对噪声数据(例如MR数据)，可以取得比传统两维特征空间更好的绘制效果。

附图说明

现在利用参照附图的方式对本发明的非限制性实施例进行详细说明，本发明的优势和特点将会变得更容易理解，其中：

图1是本发明计算两维灰度-3D SUSAN算子传递函数空间的流程图。

图2为通过本发明计算出的3D SUSAN算子在CT手数据实施例的某一片上显示结果。

图3是本发明计算每个体素的3D SUSAN算子的流程图。

图4a-图4b是利用本发明的两维灰度-3D SUSAN算子特征空间绘制出的最终结果以及作为对比实验的两维灰度-梯度幅值特征空间绘制出的最终结果的CT手数据应用实例。

图5a-图5b是利用本发明的特征空间绘制出的最终结果以及作为对比实验的两维灰度-梯度幅值特征空间绘制出的最终结果的MRA血管数据应用实例。

具体实施方式

下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进。这些都属于本发明的保护范围。

本实施例提供一种基于灰度-3D SUSAN算子两维直方图体可视化方法，具体步骤过程如下：

(1)整个计算两维灰度-3D SUSAN算子传递函数的流程图可以参照图1。读入三维体数据，计算体数据中每个体素的3D SUSAN算子，计算过程可以参考图3来理解，位于点

的计算公式为：

$\mathrm{SUSAN}\left({\stackrel{\&RightArrow;}{r}}_i\right)=\left\{\begin{array}{ccc}g-\mathrm{USAN}\left({\stackrel{\&RightArrow;}{r}}_i\right)& \mathrm{if}& \mathrm{USAN}\left({\stackrel{\&RightArrow;}{r}}_i\right)<g\\ 0& & \mathrm{otherwise}\end{array}\right.,$

$\mathrm{USAN}\left({\stackrel{\&RightArrow;}{r}}_i\right)=\underset{|{\stackrel{\&RightArrow;}{r}}_j-{\stackrel{\&RightArrow;}{r}}_i|\&le;\mathrm{\&delta;}}{\mathrm{\&Sigma;}}{e}^{-{\left(\frac{I\left({\stackrel{\&RightArrow;}{r}}_j\right)-I\left({\stackrel{\&RightArrow;}{r}}_i\right)}{T}\right)}^6}$

其中T在本实施例中设置为20，δ＝4，g＝200。然后统计每个灰度-3D SUSAN算子数据对(I,SUSAN)，在实施例的体数据中所含体素的个数IS_NUM(I,SUSAN)，并按照下面的仿射变换计算每个数据对在两维直方图中的显示亮度值b_value(I,SUSAN)：

$b\_\mathrm{value}(I,\mathrm{SUSAN})=\frac{\log (\mathrm{IS}\_\mathrm{NUM}(I,\mathrm{SUSAN}))}{{\mathrm{ISLOG}}_{\max }}$

最后构成一个以灰度值作为X轴属性、3D SUSAN算子作为Y轴属性的新两维传递函数空间。

所述参数δ是计算3D SUSAN算子采用的球形模板的半径，这样得到的是各向同性的相应，其大小依据侧重计算量还是侧重所求精度而选取不同的值。

所述参数T是灰度差阈值，既可帮助检测USAN区域的最小值，又可确定可去除的最大噪声。

所述参数g为几何阈值，针对噪声数据设置不同的大k值大小，给出最优消除噪声的能力。

所述的步骤(1)中统计体数据中具有每个灰度-3D SUSAN算子数据对(I,SUSAN)的体素的个数。

所述的步骤(1)中对其个数进行对数处理并除以对数的最大值，作为最终显示直方图的亮度值。

由图2某一切片数据可以看出，越接近两个物质的边界，3D SUSAN算子的值就越大，所以可用于判断体素是否位于边界区域。本发明的特征空间比传统的两维灰度-梯度幅值传递函数空间处理体数据，尤其是噪声数据可以获得更多有用信息，从而使用户根据自己的需要更容易设计出合适的传递函数。

(2)在灰度-3D SUSAN算子特征空间的交互界面上放置传递函数分类器组件，然后根据希望得到的绘制结果调节所设置的分类器组件，可以单独调节某个控制点，也可以调节两点之间的连线，然后观察由于调整组件而绘制的三维图像是否为期望的结果。

(3)在通过步骤(2)的操作过程之后，可以设计出合适的传递函数，然后根据设置的传递函数赋于每个体素的光学属性值绘制出三维效果图。为了增加最终绘制结果的立体感和真实性，本发明将Blinn-phong光照模型加入到体绘制过程中，其模型公式为：

$I_{\mathrm{BP}}=K_{\mathrm{am}}{I}_{\mathrm{am}}+K_{\mathrm{di}}{I}_{\mathrm{di}}(\stackrel{\&RightArrow;}{N}\&CenterDot;\stackrel{\&RightArrow;}{L})+K_{\mathrm{sp}}{I}_{\mathrm{sp}}{(\stackrel{\&RightArrow;}{N}\&CenterDot;\stackrel{\&RightArrow;}{H})}^n$

$\stackrel{\&RightArrow;}{H}=\frac{\stackrel{\&RightArrow;}{L}+\stackrel{\&RightArrow;}{V}}{|\stackrel{\&RightArrow;}{L}+\stackrel{\&RightArrow;}{V}|}$

其中，K_am、K_di和K_sp分别为环境光反射系数、漫反射系数和镜面反射系数，I_am、I_di和I_sp分别为环境光的颜色、漫反射光的颜色和镜面反射光的颜色，

是表面单位法向量，

为指向光源的单位向量，

是指向视点的单位向量，n是高光参数。

本发明在VS2008开发环境下利用OpenGL和Cg语言采用上述技术方案实现了整个算法，其中在CPU中计算属性值和特征空间，而直接体绘制的累积过程是在GPU中运算的。在运用时计算机的配置为NVIDIA GeForce GT630M显卡，2GBRAM，Intel Core i32.20GHz，Windows7操作系统。

以一个手臂的CT数据(492×240×155)和一个血管的MRA数据(256×256×256)作为应用实例来验证本发明所提出新方法的有效性和实用性。

图4a的下面部分和图5a的下面部分是本发现的两维灰度-3D SUSAN算子特征空间的直方图交互界面（灰度值作为X轴属性值，3D SUSAN算子作为Y轴），作为对比实验的是传统两维传递函数灰度-梯度幅值位于图4b的下半部分和图5b的下半部分（灰度值作为X轴，梯度幅值为Y轴）。在CT手臂数据中，通过对比图4a的下半部分图像和图4b的下半部分图像可以明显的得出本发明提出的传递函数边界信息更明显，更利于用户在特征空间识别出感兴趣区域，并且观察设置本发明图4a下半部分的传递函数得到的体绘制手臂图像(图4a上半部分)与设置图4b下半部分的传递函数而体绘制的图像(图4b上半部分)可以明显的得出，本发明的结果要比传统两维传递函数绘制的手臂骨骼图像真实，而且没有出现不需要的散乱点。在MRA血管数据中，当噪声对数据的影响非常严重时，在传统的灰度-梯度幅值直方图(图5b的下半部分)中的数据对点扩散比较严重，对于用户，尤其是没有相关体绘制数学背景和传递函数设计经验的人来说，无法确定在整个空间中如何设置合适的传递函数，而本发明的交互界面如图5a的下半部分所示，可非常明显的展示出体数据结构信息在灰度-3D SUSAN算子特征空间的聚类关系，即使是没有经验的用户也可以通过观察直方图信息很容易设计出合适的传递函数(图5a的下半部分)，绘制出期望的三维效果(图5a的上半部分)。而且通过对最终绘制结果的比较，本发明在MRA血管数据中的最终绘制结果要比传统传递函数的绘制结果(图5b的上半部分)更加真实和直观。

通过本发明在CT数据和MRA数据的应用实例，与传统两维灰度-梯度幅值的比较，可以得出从本发明的特征空间，用户可以得到更多有用的结构信息，更容易设计出使用户满意的合适传递函数，绘制更直观的图像，尤其是对于受噪声污染比较严重的MR数据。同时，本发明的原理比较简单，交互操作也非常直观，对医学数据的三维可视化是有非常重要的实用价值的。

由以上实施例可以看出，由于SUSAN算法在进行边缘检测时对局部噪声不敏感，具有较好的抗噪能力，本发明将传统的SUSAN算法扩展到三维，作为直接体绘制中传递函数的一个新属性值，结合体数据本身的灰度值构成一个新两维传递函数空间，然后根据用户在该特征空间设计的传递函数，绘制出使人满意的体绘制结果。本发明首次将具有抗噪能力的SUSAN算子作为新属性引入到传递函数设计过程中，构造出两维灰度-3D SUSAN算子传递函数；本发明的原理简单，易于在基于硬件和基于软件的体绘制方法中实现；对于噪声数据，尤其是MR数据本发明要比传统的两维灰度-梯度幅值传递函数空间获得更多有用的结构信息，从而使用户更容易设计出合适的传递函数，得到用户满意的绘制图像。

以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变形或修改，这并不影响本发明的实质内容。

Claims (7)

1.一种基于灰度-3D SUSAN算子两维直方图体可视化方法，其特征在于，包括以下步骤：

(1)读入三维体数据，首先计算出体数据中每个体素的3D SUSAN算子，然后计算出灰度-3D SUSAN算子两维传递函数空间；

所述计算出体数据中每个体素的3D SUSAN算子，具体为：

$\mathrm{SUSAN}\left({\stackrel{\&RightArrow;}{r}}_i\right)=\left\{\begin{array}{cc}g-\mathrm{USAN}\left({\stackrel{\&RightArrow;}{r}}_i\right)& \mathrm{ifUSAN}\left({\stackrel{\&RightArrow;}{r}}_i\right)<g\\ 0& \mathrm{otherwise},\end{array}\right.$

$\mathrm{USAN}\left({\stackrel{\&RightArrow;}{r}}_i\right)=\underset{|{\stackrel{\&RightArrow;}{r}}_j-{\stackrel{\&RightArrow;}{r}}_i|\&le;\mathrm{\&delta;}}{\mathrm{\&Sigma;}}{e}^{-{\left(\frac{I\left({\stackrel{\&RightArrow;}{r}}_j\right)-I\left({\stackrel{\&RightArrow;}{r}}_i\right)}{T}\right)}^6}$

其中

是位于点

的灰度值，δ代表体素周围某一小邻域的尺度， $|{\stackrel{\&RightArrow;}{r}}_j-{\stackrel{\&RightArrow;}{r}}_i|=\sqrt{{(r_{\mathrm{jx}}-r_{\mathrm{ix}})}^2+{(r_{\mathrm{jy}}-r_{\mathrm{iy}})}^2+{(r_{\mathrm{jz}}-r_{\mathrm{iz}})}^2},$ T为灰度差阈值，g＝kNUM_max为几何阈值，k为调节参数，NUM_max为体素δ邻域内所能取得最大值，

是位于点

周围球形邻域δ内各体素的灰度值与的相似度的和值；

然后统计体数据中具有每个灰度-3D SUSAN算子数据对(I,SUSAN)的体素的个数，对其个数进行对数处理并除以对数的最大值，作为最终显示直方图的亮度值；

最后由灰度值作为新特征空间X轴的属性值，3D SUSAN算子作为Y轴属性值，构成一个新的两维传递函数空间；

(2)根据对灰度-3D SUSAN算子特征空间直方图的理解，设置二维传递函数分类器组件并根据期望的绘制效果做出调整；

(3)根据确定的传递函数所赋予每个体素的光学属性值绘制出最终结果。

2.根据权利要求1所述的基于灰度-3D SUSAN算子两维直方图体可视化方法，其特征在于，所述参数δ是计算3D SUSAN算子采用的球形模板的半径，这样得到的是各向同性的相应，其大小依据侧重计算量还是侧重所求精度而选取不同的值。

3.根据权利要求1所述的基于灰度-3D SUSAN算子两维直方图体可视化方法，其特征在于：所述参数g为几何阈值，依据不同的数据对其选择不同的调节参数k，给出最优消除噪声的能力。

4.根据权利要求1所述的基于灰度-3D SUSAN算子两维直方图体可视化方法，其特征在于所述的步骤(2)中放置二维传递函数分类器组件，给出不同形状的分类组件。

5.根据权利要求4所述的基于灰度-3D SUSAN算子两维直方图体可视化方法，其特征在于所述二维传递函数分类器组件采用的为四边形的组件，该组件有4个控制点和4个相互之间的连线组成，在调整组件时，能单独调节某个控制点或者调节两点之间的连线。

6.根据权利要求1所述的基于灰度-3D SUSAN算子两维直方图体可视化方法，其特征在于所述的步骤(3)中绘制最终结果图时，加入Blinn-Phong光照模型。

7.根据权利要求6所述的基于灰度-3D SUSAN算子两维直方图体可视化方法，其特征在于所述的步骤(3)中绘制最终结果图时，在光线投射的累积过程中采用从前向后的累积方式。

Images