CN103049930B - 三维医学图像动态高斯体模型绘制方法 - Google Patents

Abstract

本发明公开了一种三维医学图像动态高斯体模型绘制方法，包括以下步骤：（a）把图像序列保存为三维纹理以用于后面的纹理检索采样；（b）由视点到图像序列的最表面外层像素引射线穿过体数据，对每条射线上进行等间隔采样，得到采样点处的密度值；（c）使用传递函数将采样点的密度值映射为颜色和不透明度；（d）根据当前采样点与焦点的距离和角度信息，构建高斯模型，对采样点的不透明度进行修正；（e）沿光线方向对所有采样点的不透明度和颜色进行累加，最终得到投射到屏幕一点的颜色值。该方法可以清晰显示焦点区域处的内部结构同时保留了焦点周围的组织结构，而不像传统的统一剖切，达到了焦点加上下文的效果。

Description

三维医学图像动态高斯体模型绘制方法

技术领域

本发明属于图像处理技术领域，具体来说是一种三维医学图像动态高斯体模型绘制方法。

背景技术

体绘制技术是医学可视化研究的一个重要方面，特别是在医疗诊断和辅助外科等方面，已经成为一种必不可少的工具。随着医学成像设备的快速发展，由设备产生的数据量也在逐步增加，而临床感兴趣的组织或部位，相对于整个体数据来说是较小的，例如头部的血管、肾脏或肝脏内的病变、肺部的肿瘤等等。医生通过对数据的观察，进而得出结论，是一个非常复杂的任务，必须考虑到可疑病变组织的解剖结构，包括病灶的大小、形状，以及其空间位置和周边组织结构的关联。利用体绘制技术辅助，医生可以对可疑病灶部位进三维的观察，可以有效的观察不同部位或组织的空间位置信息，但在很多情况下，直接体绘制的显示效果并不理想，如头部血管造影数据中，头部的血管是医生观察的重要组织，但头骨的密度值和造影的血管密度值非常接近，利用直接体绘制中一维传递函数的调节很难调节出好的显示效果。

采用剖切方法去除遮挡感兴趣区域的组织是比较直观的方法，但剖切方法很难保留感兴趣组织的上下文信息。为了改进上面的不足，weiskopf等人提出构建任意的剖切曲面，用户通过对剖切曲面的选择来实现感兴趣区域上下文信息的保留，但在交互性上明显不够灵活，且对感兴趣区域的选择仍不够理想；

基于焦点和上下文的可视化方法，周等人提出采用距离来增强和减弱不同组织，其主要思想是引入焦点(或感兴趣点)，通过当前采样点到焦点的距离来作用于当前采样点的不透明度；焦点的引入可以确定用户目标区域，有利于快速锁定三维空间中得焦点，周的方法解决了感兴趣区域定位的问题，但对感兴趣区域的显示效果并不是但好；StefanBruckner等提出保留上下文环境体绘制方法，主要是运用梯度构造出不透明度衰减函数，根据采样点与视点的距离和光照强度来对不透明度衰减函数进行衰减，可以在保留组织结构的同时达到透视组织内部的效果。StefanBruckner的方法采用了全局处理，所以得到的透视效果也是全局透视，对非感兴趣区域也产生了不必要的影响。

在K·恩格尔的“利用虚拟切割的背景中体积可视化的***及方法”专利中，提出首先切割区域为内外两部分，再分别使用两个不同的传递函数对内外两个区域进行调整，以实现显示保留上下文的类似虚拟切割的效果（在专利原文中把保留上下文翻译为背景）。其目的同本发明类似，都是为了能够观察到被遮挡住的感兴趣区域信息。K·恩格尔的实现方法首先要对数据场进行切割，然而要完成精准的切割，既能够保证感兴趣区域和周围的组织信息被完整的保留下来，又要把遮挡区域划分出来，这显然是一个相当费时与交互性较难的工作，大大的降低了临床应用性的可能。

发明内容

本发明的目的在于提供一种三维医学图像动态高斯体模型绘制方法，在清晰显示感兴趣区域的同时，保留与该感兴趣区域周围相关的结构信息，以克服诸如剖切等方法带来的交互复杂性和结构信息的丢失以及现有方法中目标区域的模糊等不良现象。

本发明的目的可通过以下的技术措施来实现：三维医学图像动态高斯体模型绘制方法，该方法包括如下步骤：

（a）将三维图像序列存储为三维纹理数据，用于后面的三维图像序列的密度值索引，并绘制三维图像序列的包围盒；

（b）计算视点位置：调用开放式图形库中的glGetFloatv函数，获得一个4×4的模型视图矩阵，再计算出该矩阵的逆矩阵，把该逆矩阵最后一行的四个量作为视点位置；由视点到三维图像序列包围盒最外层的每一个像素都引一条射线穿过三维图像体数据，沿每条射线等间隔采样图像序列，得到采样点处的密度值；

（c）使用传递函数将采样点的密度值映射为颜色和不透明度；

（d）任选三维纹理数据中的一点标定为焦点，根据当前采样点与焦点的距离和角度信息，构建高斯模型，对当前采样点处的不透明度进行修正；

（e）沿射线方向对所有采样点重复步骤（c）、（d）的处理，得到所有采样点的颜色和修正后的不透明度，并且将所有采样点的颜色和修正后的不透明度进行累加，并把累加颜色值作为当前射线对应的屏幕像素点的颜色值；

（f）对所有的射线做同样的处理，最终得到投射到屏幕上的所有像素点的颜色值。

本发明中，所述步骤（a）中，将三维图像序列存储为三维纹理数据时首先对三维图像序列数据进行归一化处理，并且最大纹理坐标被确定为满足条件：T_c=2ⁿ，c∈{x，y，z}，其中x表示冠状面的垂线，y表示矢状面的垂线，z表示横断面的垂线，T_c是在c方向的纹理坐标，n是满足2ⁿ≥S_c的最小整数，S_c是在c方向的图像序列的维度；所述步骤（a）中，绘制三维图像序列的包围盒是以所述T_c，c∈{x,y,z}分别为长、宽、高绘制的立方体。

本发明中，所述步骤（b）中对每条射线上进行等间隔采样，采样间距在交互与静态时采用不同的值，所述交互指用户对图像进行的旋转、缩放、平移、裁剪操作，静态指用户对图像进行的静止操作。本发明中，所述步骤（b）中对每条射线上进行等间隔采样，采样间距在交互与静态时采用不同的值，即当用户进行图像旋转、缩放、平移、裁剪等交互操作时采用的采样间隔是1/128mm，在静态时采用的采样间隔是1/1024mm，所述交互指用户对图像进行的旋转、缩放、平移、裁剪操作，静态指用户对图像进行的静止操作。采样间距在交互与静态时采用不同的值。在交互时使用较大的采样间距，而在静态时使用更精确的采样间距。本文在交互与静态时分别采用的采样间隔是1/1024mm和1/128mm，这样既保证了渲染质量，又有很好的交互效果。

本发明中，所述步骤（c）使用传递函数将采样点的密度值映射为颜色和不透明度后，把不透明度乘以采样点的归一化的梯度值作为新的不透明度，然后再进行步骤（d）的处理。

本发明中，所述步骤（d）中构建动态高斯函数的过程是：所述焦点与当前采样点连线构成的基准线和视线的夹角为以及焦点与当前采样点的距离为参数即基准线的长度，同时引入参数k，来构建动态高斯模型，把为高斯函数的输入，高斯函数的均值为0，高斯函数方差为输出作为采样点处的新的不透明度的衰减系数，可得到如下高斯函数方程：

${\mathrm{\δ}}_{p_i}=k(1-{\cos \mathrm{\θ}}_{P_i})$

$x_{P_i}=d_{P_i}$

$m\left(P_i\right)=G\left(P_i\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\δ}}_{P_i}}{e}^{-\frac{{x_{P_i}}^2}{2{\mathrm{\δ}}_{P_i}^2}}$

其中，k为引入的调节参数，主要用于调整方差，控制感兴趣区域大小；为焦点与当前采样点连线和其所在视线的夹角，即为当前采样点处高斯函数的方差，当前采样点到焦点的距离为以0为均值，为方差，为输入，m(P_i)是对当前采样点构建的动态高斯函数修正值。

本发明步骤（d）中还可以对得到的高斯函数的输出作用于采样点处的新的不透明度，从而实现对采样点处的新的不透明度的修正： $\mathrm{\α}\left(P_i\right)={\mathrm{\α}}_{\mathrm{tf}}\left(f_{P_i}\right)\·{\left|\right|g_{p_i}\left|\right|}_{[0..1]}m\left(P_i\right)$

其中，为当前采样点经过传递函数映射得到的不透明度，为当前采样点的归一化的梯度值，m(P_i)是对当前采样点构建的动态高斯函数修正值，最终得到该采样点处用于累加计算的不透明度α(P_i)。

基准线表示当前采样点到焦点的空间位置关系，基准线的长度越短，当前采样点越靠近感兴趣区域的中心，应该大幅度衰减感兴趣区域前的遮挡部分。如图3，在两条光线上的两个采样点S₁和S₂，这两个采样点的视场角相同，于是它们就具有相同的高斯曲线形式，由图中可看到，采样点1的基准线长度较小，距离焦点较近，就会得到较大的衰减值，反之，采样点2因为距离焦点较远而衰减变弱。认为当衰减值低于虚线标注区域以下时（一般取峰值的10%），则不再进行衰减，采样点2刚好位于虚线上，则采样点2为感兴趣区域的边界点。

本发明把基准线与视线的夹角称之为视场角。视场角越小，所有位于感兴趣区域内的采样点被限制在一个很小的范围内，也就是说，感兴趣区域变小了。如图4，在两条光线上的两个采样点S₁和S₂，这两个采样点的基准线长度相同，采样点1的视场角较小，构建的高斯函数均值较小，曲线很陡峭。比较两个采样点会发现，对于同一基准线长度，采样点1的衰减值较小，刚好位于感兴趣区域的边界上，即感兴趣区域范围为x≤d；而采样点2仍具有较大的衰减值，要衰减到感兴趣区域的边界处，可以对x＞d的更大区域进行衰减。采样点的感兴趣区域范围大于采样点1的。

为了控制感兴趣区域的大小，本发明引入参数k，即调整高斯函数的均值为k·(1-cosθ)。通过调整参数k，控制透视窗口的大小。当k增大时，高斯函数的方差也随之增大，高斯曲线变的就越平缓，感兴趣区域变大。实验统计，参数k的范围一般在[0,2]之间。按照一条视线上对数据场采样得到采样点的先后顺序排列，每个采样点得到的参数和也具有一定的连续性，不会出现过大的跳跃现象，这样就保证了，被衰减部位或组织的边缘比较平缓，不会出现剖切显示中边界粗糙的现象，可以提高成像效果。

所述感兴趣区域的选择是通过把三维纹理数据中的一点标定为焦点，再以该点为中心来确定一个区域，该区域大小应该是可调节的。

本发明中，所述步骤（e）中在着色器程序中对颜色和不透明度进行累加计算，使用两层while嵌套循环进行迭代计算。

本发明中，所述步骤（f）在计算最终屏幕显示的颜色时，只是位于图像序列的包围盒内的像素点才会进行投射光线，并计算最终像素点颜色，位于包围盒之外的像素点颜色直接设置为背景色。

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

（1）本发明所述方法与周等人提出的方法相比较，引入了一个全新的高斯函数模型，该模型同时考虑到以焦点与当前采样点连线和视线间的角度信息和当前采样点与焦点的距离信息。通过引入该复杂模型，在感兴趣区域清晰显示的同时，保留与该感兴趣区域相关结构信息，克服剖切带来的复杂交互和结构信息的丢失以及现有方法中目标区域的模糊现象，大大提升了成像质量。在感兴趣区域之外的其它地方保留了原有的渲染效果，而把衰减效果较好的应用在感兴趣区域，同时在感兴趣区域边界处会有一个逐渐衰减减弱的过渡区域，使得整体渲染效果变得很自然。

（2）本发明与weiskopf的方法相比较，简化了用户交互操作，用户只简单确定焦点区域，通过简单参数调节，便可以达到很好的对焦点区域的透视效果。

（3）本发明与K恩格尔的发明相比较，通过引入动态高斯模型，只要确定一个焦点，即可以自动的去除遮挡区域，显示出感兴趣区域，同时保留下周围的组织信息，明显比K恩格尔的发明更加的方便，易操作。且K恩格尔引入了两个传递函数分别用来调整分割区域内部和外部，而本发明保留了仅仅使用一个传递函数，可以大大的节省程序在运行时的内存空间，提升处理速度。

附图说明

图1为本发明中光线投射算法三维示意图：对屏幕上的像素点投射一条光线，穿过体纹理，沿着光线对体纹理进行采样；

图2为本发明中光线投射算法二维示意图；

图3为本发明中采样点的视场角(θ)相同，基准线长度(d)不同的情况下构建的高斯函数模型；

图4为本发明中采样点的视场角(θ)不同，基准线长度(d)相同的情况下构建的高斯函数模型。

具体实施方式

以下是采用本发明三维医学图像动态高斯体模型绘制方法对一套头部的CT数据进行成像的具体过程：

实验采用的数据来源为OsiriX网站；GPU为QuadroFX1700显卡。

步骤1，读入CT图像，图像的大小为512×512×460，首先对图像序列数据进行归一化处理，计算最大纹理坐标为满足条件：T_c=2ⁿ，c∈{x,y，z}。T_c是在c方向的纹理坐标，n是满足2ⁿ≥S_c的最小整数，S_c是在c方向的图像序列的维度。把图像序列保存为三维纹理以用于三维体数据的密度值索引。

步骤2，由开放式图形库（OpenGL，OpenGraphicsLibrary）的glGetFloatv函数可以获得一个4X4的模型视图矩阵，再计算出该矩阵的逆矩阵invMVM[4X4]，把该逆矩阵的最后一行的四个量作为视点位置，如下：

(invMVM[12],invMVM[13],invMVM[14],invMVM[15])。

如图1和图2的光线投射算法示意图，由视点发射出光线进入三维数据场，记录下第一个进入数据场的采样点，称为入射点。使用传递函数把采样点处得到的灰阶转换成颜色和不透明度。根据中心差值算法计算出当前采样点的梯度向量和梯度大小，使用梯度大小对不透明度进行线性调整，以增强边界效果；同时采用Blinn-Phong光照模型，以增强显示效果。

三维连续数据场用f(P_n)表示，P_n表示三维数据场中采样点的坐标，在采样点P_i处的密度值和梯度为 $f_{P_i},g_{P_i}=\▿f\left(p_i\right),{\left|\right|g_{P_i}\left|\right|}_{[0..1]}$ 为归一化后梯度的模。使用传递函数将当前采样点的密度值映射为颜色和不透明度分别为 $c\left(P_i\right)=c_{\mathrm{tf}}\left(f_{P_i}\right),$ $\mathrm{\α}\left(P_i\right)={\mathrm{\α}}_{\mathrm{tf}}\left(f_{P_i}\right).$

为了增强边界效果，做如下优化处理：同时为了增强显示效果，把不透明度再乘以如下光照模型，光照模型计算公式为： $L\left(P_i\right)=k_a+k_d\·\left|\right|\hat{L}\·{\hat{g}}_{P_i}\left|\right|+k_e\·{\left(\right||\hat{H}\·\hat{g}\left(P_i\right)|\left|\right)}^m,$ 其中k_a、k_d、k_e分别为环境光系数、漫反射光系数、镜面反射光系数，m是镜面反射光汇聚指数，为视线入射方向向量，为半角向量。

步骤3，采用高斯函数来构建不透明度衰减函数m(P_i)。如图3和图4左图所示，首先由用户选定感兴趣点(即焦点)，可以通过三维定位标定于数据场中的一点，在图3和图4中感兴趣点为红色的点；黑色的点是当前视线上的采样点，采用光线投射算法进行渲染的过程，是沿每条视线，对采样点的光学属性积分的过程，为了降低非感兴趣组织采样点对积分过程的影响，以焦点与当前采样点连线和视线的夹角，图3和图4中为S₁与红色焦点的连线，与视线的夹角θ₁，以及当前采样点S₁与焦点的距离d₁为参数，来构建高斯函数，图3和图4右图，即为构建的高斯函数。具体的构建高斯函数的方法是，以距离d₁为高斯函数的输入，高斯函数的均值为0，高斯函数方差为(1-cosθ₁)。

因为该高斯函数有两个自变量，以下使用控制变量法对该高斯函数进行说明，即保证一个变量不变，讨论另一个变量变化所引起的作用。

对于具有相同的视场角，而基准线的长度不同的两个采样点，如图3所示，采样点1同采样点2的视场角相同，则两者对应的高斯曲线相同，但基准线越长，距离焦点越远，对该采样点的衰减就越小，见S₂；对于基准线的长度相同，而视场角不同的两个采样点，如图4所示，在两条光线上的两个采样点S₁和S₂，这两个采样点的基准线长度相同，采样点1的视场角较小，构建的高斯函数均值较小，曲线很陡峭。比较两个采样点会发现，对于具有相同的基准线长度的采样点，采样点1的衰减值较小，刚好位于感兴趣区域的边界上，即感兴趣区域范围为x≤d；而采样点2仍具有较大的衰减值，要衰减到感兴趣区域的边界处，可以对x＞d的更大区域进行衰减。采样点的感兴趣区域范围大于采样点1的。具体公式表示如下：

${\mathrm{\δ}}_{p_i}=k(1-{\cos \mathrm{\θ}}_{P_i})$

$x_{P_i}=d_{P_i}$

$m\left(P_i\right)=G\left(P_i\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\δ}}_{P_i}}{e}^{-\frac{{x_{P_i}}^2}{2{\mathrm{\δ}}_{P_i}^2}}---\left(2\right)$

k为引入的调节参数，主要用于调整方差，控制感兴趣区域大小；为焦点与当前采样点连线和其所在视线的夹角，即为当前采样点处高斯函数的方差，当前采样点到焦点的距离为以0为均值，为方差，构建高斯函数，为输入，然后将该高斯函数的输出作用于不透明度，从而实现对不透明度的修正：

$\mathrm{\α}\left(P_i\right)={\mathrm{\α}}_{\mathrm{tf}}\left(f_{P_i}\right)\·{\left|\right|g_{p_i}\left|\right|}_{[0..1]}m\left(P_i\right)$

步骤4，沿着光线穿过体数据场，对数据场不断采样，直到判断光线投射出数据场为止，如图1和图2所示。判断方法是通过比较当前采样点处的x，y，z三个方向纹理坐标是否大于1或小于0，只要有一个成立就判定为光线投射出了数据场。沿光线方向采用离散近似积分由前向后对采样点的光学特性进行合成,最后得到该视线在成像平面上一点的颜色值和不透明度值。离散化后具体的合成公式如下：

α_i=α_i-1+α(P_i)·(1-α_i-1)

c_i=c_i-1+c(P_i)·α(P_i)·(1-α_i-1)

其中α(P_i)和c(P_i)为当前采样点的不透明度和颜色值，是经过以上步骤2和步骤3优化后的颜色和不透明度，α_i-1、c_i-1为视线进入当前采样点之前累积的不透明度和颜色值，α_i、c_i为视线穿过当前采样点后的累积不透明度和颜色值。

步骤5，沿射线方向对所有采样点的不透明度和颜色进行累加，并把累加颜色值作为当前射线对应的屏幕像素点的颜色值；

步骤6，对所有的射线做同样的处理，最终得到投射到屏幕上的所有像素点的颜色值。

本发明的实施方式不限于此，根据本发明的上述内容，按照本领域的普通技术知识和惯用手段，本发明中不透明度衰减函数也可以用其它类似高斯函数的形式来构造，其中的参数取值也可以选择像素灰度、局部熵、光照信息等。因此在不脱离本发明上述基本技术思想前提下，本发明还可以做出其它多种形式的修改、替换或变更，均落在本发明权利保护范围之内。

Claims (6)

1.三维医学图像动态高斯体模型绘制方法，该方法包括如下步骤：

(a)将三维图像序列存储为三维纹理数据，用于后面的三维图像序列的密度值索引，并绘制三维图像序列的包围盒；

(b)计算视点位置：调用开放式图形库中的glGetFloatv函数，获得一个4×4的模型视图矩阵，再计算出该矩阵的逆矩阵，把该逆矩阵最后一行的四个量作为视点位置；由视点到三维图像序列包围盒最外层的每一个像素都引一条射线穿过三维图像体数据，沿每条射线等间隔采样图像序列，得到采样点处的密度值；

(c)使用传递函数将采样点的密度值映射为颜色值和不透明度；

(d)任选三维纹理数据中的一点标定为焦点，根据当前采样点与焦点的距离和角度信息，构建高斯模型，对当前采样点处的不透明度进行修正；

所述步骤(d)中构建动态高斯函数的过程是：所述焦点与当前采样点连线构成的基准线和视线的夹角为以及焦点与当前采样点的距离为参数即基准线的长度，同时引入参数k，来构建动态高斯模型，把为高斯函数的输入，高斯函数的均值为0，高斯函数方差为输出作为采样点处的新的不透明度的衰减系数，得到如下高斯函数方程：

${\mathrm{\δ}}_{p_i}=k(1-{\mathrm{cos\θ}}_{P_i})$

$x_{P_i}=d_{P_i}$

$m\left(P_i\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\δ}}_{P_i}}{e}^{-\frac{{x_{P_i}}^2}{2{\mathrm{\δ}}_{P_i}^2}}$

其中，k为引入的调节参数，用于调整方差，控制感兴趣区域大小；P_i为当前采样点，即为当前采样点处高斯函数的方差，m(P_i)是对当前采样点构建的动态高斯函数修正值；

(e)沿射线方向对所有采样点重复步骤(c)、(d)的处理，得到所有采样点的颜色值和修正后的不透明度，并且将所有采样点的颜色值和修正后的不透明度进行累加，并把累加颜色值作为当前射线对应的屏幕像素点的颜色值；

(f)对所有的射线做同样的处理，最终得到投射到屏幕上的所有像素点的颜色值。

2.根据权利要求1所述的三维医学图像动态高斯体模型绘制方法，其特征在于：所述步骤(a)中，将三维图像序列存储为三维纹理数据时首先对三维图像序列数据进行归一化处理，并且最大纹理坐标被确定为满足条件：T_c＝2ⁿ，c∈{x,y,z}，其中x表示冠状面的垂线，y表示矢状面的垂线，z表示横断面的垂线，T_c是在c方向的纹理坐标，n是满足2ⁿ≥S_c的最小整数，S_c是在c方向的图像序列的维度；所述步骤(a)中，绘制的三维图像序列的包围盒是以所述T_c，c∈{x,y,z}分别为长、宽、高绘制的立方体。

3.根据权利要求1所述的三维医学图像动态高斯体模型绘制方法，其特征在于：所述步骤(b)中对每条射线上进行等间隔采样，采样间距在交互与静态时采用不同的值，所述交互指用户对图像进行的旋转、缩放、平移、裁剪操作，静态指用户对图像进行的静止操作。

4.根据权利要求1所述的三维医学图像动态高斯体模型绘制方法，其特征在于：所述步骤(c)使用传递函数将采样点的密度值映射为颜色值和不透明度后，把不透明度乘以采样点的归一化的梯度值作为新的不透明度，然后再进行步骤(d)的处理。

5.根据权利要求1所述的三维医学图像动态高斯体模型绘制方法，其特征在于：所述步骤(e)中在着色器程序中对颜色值和不透明度进行累加计算，使用两层while嵌套循环进行迭代计算。

6.根据权利要求1所述的三维医学图像动态高斯体模型绘制方法，其特征在于：所述步骤(f)在计算最终屏幕显示的颜色值时，位于图像序列的包围盒内的像素点才会进行投射光线，并计算最终像素点颜色值，位于包围盒之外的像素点颜色值直接设置为背景色。