CN104167013A - 一种用于突出显示体数据中目标区域的体绘制方法 - Google Patents

Abstract

本发明涉及一种用于突出显示体数据中目标区域的体绘制方法，属于计算机图形学中科学可视化直接体绘制技术领域。首先计算原始数据的多个特征，对特征数据进行评估选择合适的特征组成传输函数。用归一化割方法分割二维直方图图像对目标数据进行分离。如果目标分离不够精细，就继续选择合适的特征对当前结果进一步分离，从而逐步得到精细的分离结果。之后利用分离结果合成最终的传输函数。最后在体绘制时利用传输函数突出显示目标区域，得到理想的可视化效果。本发明方法强化了传输函数的特征区分能力，使用了基于二维直方图图像分割和层次聚类的交互方式，操作便捷。而且适用于多种数据。

Description

一种用于突出显示体数据中目标区域的体绘制方法

技术领域

本发明涉及一种用于突出显示体数据中目标区域的体绘制方法，属于计算机图形学中科学可视化直接体绘制技术领域。

背景技术

直接体绘制方法是探索三维图像整体信息和内部结构的一种有效手段。直接体绘制方法通过传输函数为不同的体素设置不同的颜色和不透明度，然后采用光学积分模型对体数据采样，得到绘制结果。其中传输函数的设计是决定体绘制结果好坏的关键。通过传输函数，用户可以从原始体数据中分离出用户感兴趣的区域，并对其突出显示。按使用的特征维度分类，传输函数的几种类型及相应存在的问题如下：

(1)以灰度作为特征的一维传输函数是目前应用最广的传输函数，它具有良好的交互性和较好的可视化效果。但对于具有相同灰度的不同特征，一维灰度传输函数并不能很好地区分。

(2)二维传输函数在科研领域十分流行。它凭借特征维度的优势相比一维传输函数拥有更强的特征区分能力。但是将三维以上的体素信息投射到二维传输函数上也会引起维度丢失，而且二维传输函数的交互不直观，操作不便捷。

(3)三维以上的传输函数具有更好的特征维度优势，特征区分能力更突出。但是其传输函数的交互非常复杂，无论在科研领域还是工业界都没有得到广泛的应用。

除此之外，还有很多应用于某种特定体数据的传输函数设计方法，如医学中用于突出血管的传输函数设计。也有很多强化交互模式的传输函数设计方法，如直接与体绘制结果交互的所见即所得的传输函数设计。

以上这些方法都有自己独特的优势，但不能同时具有简单的交互模式，强大的特征区分能力和多种数据应用范围这些特点。如何让传输函数对多种数据同时拥有强大的特征区分能力和交互的便捷性仍然是一个亟待解决的问题。

发明内容

本发明的目的是提出一种用于突出显示体数据中目标区域的体绘制方法，在二维传输函数的基础上，依次使用不同的特征对原始体数据进行处理，以求利用特征维度的拓展，逐步精确目标区域的分离。同时利用面向数据分类的特征空间度量差异性的评估方法，在分离目标区时下一二维传输函数特征维度的选取，使得选用特征之间的差异性较大，从而逐步构造合适的特征序列。

本发明提出的用于突出显示体数据中目标区域的体绘制方法，包括以下步骤：

(1)分别计算体数据的梯度特征、曲率特征、分型特征、光滑度、线型特征和三维不变矩，具体过程如下：

(1-1)设体数据中任意一个体素的梯度特征为G，G＝|g|，G为g的模，g为体数据中体素的梯度向量，f为体数据灰度特征在x、y、z三个方向上的函数映射，灰度特征为体素的数值特征，为的转置，则体数据的梯度特征为体数据中所有体素的梯度特征按照体素在体数据中的位置排列得到的一个三维矩阵；

(1-2)设体数据中任意一个体素的曲率特征为C_m，κ₁为体素所在位置的主曲率，κ₂为体素所在位置的副曲率，T为几何张量Z的迹，F为几何张量Z的弗罗贝尼乌斯范数，几何张量其中H_e是体素的海森矩阵，P_m是体素所在位置的等值面切平面矩阵，P_m＝I-n_vn_v ^t，I为单位矩阵，n_v为体素所在位置的表面向量，n_v ^t为n_v的转置，则体数据的曲率特征为体数据中所有体素的曲率特征按照体素在体数据中的位置排列得到的三维矩阵；

(1-3)设体数据中任意一个体素的分型特征为D，对体数据中任意一个体素v_i，分型特征 $D=\frac{{\mathrm{\Σ}}_{r=1}^{k_d}(X_r-\stackrel{\‾}{X})(Y_r-\stackrel{\‾}{Y})}{{\mathrm{\Σ}}_{r=1}^{k_d}{(X_r-\stackrel{\‾}{X})}^2},$ 其中X_r＝logr， $\stackrel{\‾}{X}={\mathrm{\Σ}}_{r=1}^{k_d}\frac{X_r}{k_d},Y_r=\log u\left(\stackrel{\‾}{B}(v_i,r)\right),$ $\stackrel{\‾}{Y}={\mathrm{\Σ}}_{r=1}^{k_d}\frac{Y_r}{k_d},u\left(\stackrel{\‾}{B}(v_i,r)\right)={\mathrm{\Σ}}_{v_j\∈U(v_i,r)}f\left(v_j\right),$ U(v_i，r)为体素v_i的邻域，r为邻域半径，f(v_j)为体素v_j的灰度特征，v_j为体素v_i的邻域U(v_i，r)中任意体素，表示体素v_i的以r为半径的邻域中所有体素的灰度特征的和，k_d的取值为2，则体数据的分型特征为体数据中所有体素的分型特征按照体素在体数据中的位置排列得到的三维矩阵；

(1-4)设体数据中任意一个体素的光滑度为S，对体数据中任意一个体素v_i，光滑度其中σ²为体素v_i的邻域U(v_i，r)中所有体素的灰度特征方差，体素v_i的邻域半径r取2，则体数据的光滑度为体数据中所有体素的光滑度按照体素在体数据中的位置排列得到的三维矩阵；

(1-5)设体数据中任意一个体素的线型特征为F_flit，对体数据中任意一个体素v_i，线型特征其中为权值，取值为任意非负实数，r为体素v_i的邻域半径， $B_{k_f}=\left\{\begin{array}{c}\stackrel{\‾}{G_{k_f}},G_{k_f}>\mathrm{Mean}\left(G_{k_f}\right)\\ {0,G}_{k_f}\≤\mathrm{Mean}\left(G_{k_f}\right)\end{array}\right.,G_{k_f}={\mathrm{\Σ}}_{v_j\∈L_{k_f}}f\left(v_j\right),$ 其中为均一化后的结果， 为体素v_i的邻域内以邻域边界为两端端点且通过v_i的线段，线段总数为12r²+1，v_j为线段L_k上的一个体素，体素v_i的邻域半径r取1，则体数据的线型特征为体数据中所有体素的线型特征按照体素在体数据中的位置排列得到的三维矩阵；

(1-6)设体数据中任意一个体素的三维不变矩为T，三维不变矩T＝μ₂₀₀+μ₀₂₀+μ₀₀₂，其中μ₂₀₀、μ₀₂₀和μ₀₀₂分别为中心距，对任意一个体素v_i，中心距计算通式为：

${\mathrm{\μ}}_{\mathrm{pqt}}=\underset{v_j\∈U(v_i,r)}{\mathrm{\Σ}}{(x_1-\stackrel{\‾}{x_1})}^p{(x_2-\stackrel{\‾}{x_2})}^q{(x_3-\stackrel{\‾}{x_3})}^t\·f(x_1,x_2,x_3){\mathrm{dx}}_1{\mathrm{dx}}_2{\mathrm{dx}}_3$

其中v_j是体素v_i邻域U(v_i，r)内的任意一个体素，x₁，x₂，x₃分别为体素v_j在体数据中三个维度的坐标值，分别为体素v_i的邻域U(v_i，r)的灰度特征重心在体数据中三个维度的坐标值，f(x₁，x₂，x₃)为v_j的灰度特征，体素v_i的邻域半径r取4，则体数据的三维不变矩为体数据中所有体素的三维不变矩按照体素在体数据中的位置排列得到的三维矩阵；

(1-7)将原始体数据进行备份，备份的原始体数据记为V_copy；

(2)分别计算上述步骤(1)中得到的六个特征的特征评估量E，E＝ω₁L_n-ω₂R_n+ω₃C_n，其中ω₁，ω₂，ω₃为权值，取值为任意非负实数，L_n、R_n、C_n分别为对特征距离评估量L、特征相关性评估量R、特征价值评估量C进行均一化处理的结果，L、R、C分别为：

其中为用k均值聚类算法对特征P₁聚类后的第i_l个类的重心，P₁为上一次评估时选取的特征评估量E为最大的特征，若当前评估为第一次评估，则特征P₁为灰度特征，k₁为用k均值聚类算法对特征P₁聚类后得到的类的个数，为用k均值聚类算法对特征P₂聚类后的第j_l个类的重心，P₂为待评估的特征，k₂为用k均值聚类算法对待评估的特征P₂聚类后得到的类的个数；

R＝H(P₁)+H(P₂)-H(P₁，P₂)，其中H(P₁)为特征P₁的信息熵，H(P₂)为待评估特征P₂的信息熵，H(P₁，P₂)为特征P₁和待评估的特征P₂的联合熵；

C＝H(P₁，P₂)-H(P₁)；

(3)从步骤(2)的六个特征中选取特征评估量E为最大的特征，记为P，用特征P与原始体数据的灰度特征组成一个二维传输函数，用该二维传输函数对原始体数据中的目标区域进行分离，得到原始体数据中位置在目标区域中的体素，分离过程如下：

(3-1)用特征P与体数据的灰度特征组成一个二维统计直方图，用归一化割方法对该直方图图像进行分割，得到多个二维统计直方图图像子区域，记该多个子区域组成的集合为Y0；

(3-2)计算原始体数据中所有体素在二维统计直方图中的坐标，根据该坐标，从原始体数据中分别找到与每个子区域相对应的体素，利用光线投射方法，分别对与每个子区域相对应的体素进行体绘制，得到多个体绘制图像；

(3-3)从子区域集合Y₀中得到最终子区域集合O_final，过程如下：

(3-3-1)初始化时，最终子区域集合为空集，对Y₀中所有子区域进行从1到w排序，其中w是子区域集合Y₀中子区域的数量；

(3-3-2)从子区域集合Y₀中取出标号最小的子区域y₀，若与子区域y₀相对应的体绘制图像中只包含目标区域，则将y₀放入最终子区域集合中，并将y₀从Y₀中删除，将删除y₀后的Y₀记为Y₁；若与子区域y₀相对应的体绘制图像中只包含非目标区域，则将y₀从Y₀中删除，将删除y₀后的Y₀记为Y₁，若与子区域y₀相对应的体绘制图像中既包含目标区域又包含非目标区域，则用归一化割方法对子区域y₀进行分割，得到两个新的子区域y₁和y₂，分别将y₁和y₂加入到子区域集合Y₀中，将y₀从Y₀中删除，得到新的子区域集合Y₁，Y₁中y₁标号为N_max+1，y₂标号为N_max+2，其中N_max为子区域集合Y₀中子区域标号的最大值；

重复上述步骤，若子区域集合Y变为空集，则记最终子区域集合为O_final，进行步骤(3-4)；若上述步骤重复次数超过30次，则将Y₃₀中剩余的所有子区域加入到最终子区域集合中，并记最终子区域集合为O_final，进行步骤(3-4)；

(3-4)根据步骤(3-2)得到原始体数据中所有体素的坐标，分别分离出坐标在最终子区域集合O_final中的每个子区域的所有体素，所有体素组成的集合记为V_r，按照体素在原始体数据中的坐标，将V_r中的所有体素加入到原始体数据中，将加入V_r中的所有体素后的原始体数据记为V_n，根据下式对V_n的灰度特征进行处理，得到V_n的新的灰度特征如下：

$f_{\mathrm{new}}\left(v_h\right)=\left\{\begin{array}{cc}{f}_{\mathrm{old}}\left(v_h\right),& v_h\∈V_r\\ 0,& v_h\∉V_r\end{array}\right.$

其中v_h为V_n中的任意一个体素，f_old(v_h)为V_n中v_h的原来的灰度特征，f_new(v_h)为V_n中v_h的新的灰度特征；

(3-5)对上述V_n中的体素进行判断，若V_n中的任意体素v_h满足v_h∈V_r，则进行步骤(4)，若V_n中存在一个体素v_h，满足则根据下式对步骤(1)的6种特征数据进行处理，得到与每种特征相对应的新的特征值：

$P_{\mathrm{feature}\_\mathrm{new}}\left(v_h\right)=\left\{\begin{array}{cc}{P}_{\mathrm{feature}\_\mathrm{old}}\left(v_h\right),& v_h\∈V_r\\ 0,& v_h\∉V_r\end{array}\right.$

其中v_h为V_n中任意一个体素，P_{feature_old}(v_h)表示v_h的原来特征值，P_{feature_new}(v_h)表示v_h的新特征值，返回步骤(2)；

(4)利用步骤(3)得到的体数据V_n和步骤(1)中的备份原始体数据V_copy，合成一个传输函数，利用该传输函数得到体绘制图像，具体过程如下：

(4-1)将步骤(3)得到的体数据V_n中的体素归为一个类，记为m₀；

(4-2)根据备份的原始体数据V_copy，将步骤(3)得到的体数据V_n中灰度特征变为0的体素的灰度特征值恢复为原来的灰度特征值，将恢复灰度特征值的体素组成一个集合，记为V_o，用V_o中所有体素的梯度特征和灰度特征组成一个二维统计直方图，采用归一化割方法将此直方图分割成e个子区域，分别记为a₁、a₂…a_e；

(4-3)计算上述V_o中所有体素在梯度特征和灰度特征组成的二维统计直方图中的坐标，根据V_o中所有体素在该二维统计直方图中的坐标，将V_o划分为e个类，分别记为m₁、m₂...m_b...m_e，b＝1，2，...，e，分类方法为：对V_o中任意一个体素v_i，若v_i在该二维统计直方图中的坐标在子区域a_b内部，则将v_i归到类m_b中；

(4-4)根据显示需要，为体数据V_n中的体素的类m₀设置不透明度和颜色，为类m₁、m₂…m_e设置不透明度和互相不同的颜色，不透明度取值范围为[0，1]，颜色采用红绿蓝颜色模型，红、绿、蓝颜色分量的强度值范围为[0，255]，根据体数据中每个类的颜色和不透明度，采用光线投射算法得到最终的突出显示目标区域的体绘制图像。

本发明提出的用于突出显示体数据中目标区域的体绘制方法，具有以下优点：

(1)本方法通过依次使用不同的特征对原始体数据进行处理，拓展了传输函数的特征维度，增强了传输函数的特征区分能力，从而更好地分离体数据中的目标区域，利用分离结果对目标区域突出显示，可以得到突出显示目标区域的体绘制图像。

(2)本方法在分离体数据的目标区域时，利用面向数据分类的特征空间度量差异性的评估方法，智能指导在分离目标区时二维传输函数特征的选取，避免了用户选择特征时的不确定性，使方法的流程智能化。

(3)本方法使用基于二维直方图图像分割和层次聚类的交互方式，计算和操作十分便捷。

(4)本方法具有数据适用范围广的潜质。

附图说明

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

具体实施方式

本发明提出的用于突出显示体数据中目标区域的体绘制方法，其流程框图如图1所示，包括以下步骤：

(1)分别计算体数据的梯度特征、曲率特征、分型特征、光滑度、线型特征和三维不变矩，具体过程如下：

(1-1)设体数据中任意一个体素的梯度特征为G，G＝|g|，G为g的模，g为体数据中体素的梯度向量，f为体数据灰度特征在x、y、z三个方向上的函数映射，灰度特征为体素的数值特征，为的转置，则体数据的梯度特征为体数据中所有体素的梯度特征按照体素在体数据中的位置排列得到的一个三维矩阵；

(1-2)设体数据中任意一个体素的曲率特征为C_m，κ₁为体素所在位置的主曲率，κ₂为体素所在位置的副曲率，T为几何张量Z的迹，F为几何张量Z的弗罗贝尼乌斯(Froberlius)范数，几何张量其中H_e是体素的海森(Hessian)矩阵，P_m是体素所在位置的等值面切平面矩阵，P_m＝I-n_vn_v ^t，I为单位矩阵，n_v为体素所在位置的表面向量，n_v ^t为n_v的转置，则体数据的曲率特征为体数据中所有体素的曲率特征按照体素在体数据中的位置排列得到的三维矩阵；

(1-3)设体数据中任意一个体素的分型特征为D，对体数据中任意一个体素v_i，分型特征 $D=\frac{{\mathrm{\Σ}}_{r=1}^{k_d}(X_r-\stackrel{\‾}{X})(Y_r-\stackrel{\‾}{Y})}{{\mathrm{\Σ}}_{r=1}^{k_d}{(X_r-\stackrel{\‾}{X})}^2},$ 其中X_r＝logr， $\stackrel{\‾}{X}={\mathrm{\Σ}}_{r=1}^{k_d}\frac{X_r}{k_d},Y_r=\log u\left(\stackrel{\‾}{B}(v_i,r)\right),$ $\stackrel{\‾}{Y}={\mathrm{\Σ}}_{r=1}^{k_d}\frac{Y_r}{k_d},u\left(\stackrel{\‾}{B}(v_i,r)\right)={\mathrm{\Σ}}_{v_j\∈U(v_i,r)}f\left(v_j\right),$ U(v_i，r)为体素v_i的邻域，r为邻域半径，f(v_j)为体素v_j的灰度特征，v_j为体素v_i的邻域U(v_i，r)中任意体素，表示体素v_i的以r为半径的邻域中所有体素的灰度特征的和，k_d的取值为2，则体数据的分型特征为体数据中所有体素的分型特征按照体素在体数据中的位置排列得到的三维矩阵；

(1-4)设体数据中任意一个体素的光滑度为S，对体数据中任意一个体素v_i，光滑度其中σ²为体素v_i的邻域U(v_i，r)中所有体素的灰度特征方差，体素v_i的邻域半径r取2，则体数据的光滑度为体数据中所有体素的光滑度按照体素在体数据中的位置排列得到的三维矩阵；

(1-5)设体数据中任意一个体素的线型特征为F_flit，对体数据中任意一个体素v_i，线型特征其中为权值，取值为任意非负实数，r为体素v_i的邻域半径， $B_{k_f}=\left\{\begin{array}{c}\stackrel{\‾}{G_{k_f}},G_{k_f}>\mathrm{Mean}\left(G_{k_f}\right)\\ {0,G}_{k_f}\≤\mathrm{Mean}\left(G_{k_f}\right)\end{array}\right.,G_{k_f}={\mathrm{\Σ}}_{v_j\∈L_{k_f}}f\left(v_j\right),$ 其中为均一化后的结果， 为体素v_i的邻域内以邻域边界为两端端点且通过v_i的线段，线段总数为12r²+1，v_j为线段L_k上的一个体素，体素v_i的邻域半径r取1，则体数据的线型特征为体数据中所有体素的线型特征按照体素在体数据中的位置排列得到的三维矩阵；

(1-6)设体数据中任意一个体素的三维不变矩为T，三维不变矩T＝μ₂₀₀+μ₀₂₀+μ₀₀₂，其中μ₂₀₀、μ₀₂₀和μ₀₀₂分别为中心距，对任意一个体素v_i，中心距计算通式为：

${\mathrm{\μ}}_{\mathrm{pqt}}=\underset{v_j\∈U(v_i,r)}{\mathrm{\Σ}}{(x_1-\stackrel{\‾}{x_1})}^p{(x_2-\stackrel{\‾}{x_2})}^q{(x_3-\stackrel{\‾}{x_3})}^t\·f(x_1,x_2,x_3){\mathrm{dx}}_1{\mathrm{dx}}_2{\mathrm{dx}}_3$

其中v_j是体素v_i邻域U(v_i，r)内的任意一个体素，x₁，x₂，x₃分别为体素v_j在体数据中三个维度的坐标值，分别为体素v_i的邻域U(v_i，r)的灰度特征重心在体数据中三个维度的坐标值，f(x₁，x₂，x₃)为v_j的灰度特征，体素v_i的邻域半径r取4，则体数据的三维不变矩为体数据中所有体素的三维不变矩按照体素在体数据中的位置排列得到的三维矩阵；

(1-7)将原始体数据进行备份，备份的原始体数据记为V_copy；

(2)分别计算上述步骤(1)中得到的六个特征的特征评估量E，E＝ω₁L_n-ω₂R_n+ω₃C_n，其中ω₁，ω₂，ω₃为权值，取值为任意非负实数，L_n、R_n、C_n分别为对特征距离评估量L、特征相关性评估量R、特征价值评估量C进行均一化处理的结果，L、R、C分别为：

其中为用k均值聚类算法对特征P₁聚类后的第i_l个类的重心，P₁为上一次评估时选取的特征评估量E为最大的特征，若当前评估为第一次评估，则特征P₁为灰度特征，k₁为用k均值聚类算法对特征P₁聚类后得到的类的个数，为用k均值聚类算法对特征P₂聚类后的第j_l个类的重心，P₂为待评估的特征，k₂为用k均值聚类算法对待评估的特征P₂聚类后得到的类的个数；

R＝H(P₁)+H(P₂)-H(P₁，P₂)，其中H(P₁)为特征P₁的信息熵，H(P₂)为待评估特征P₂的信息熵，H(P₁，P₂)为特征P₁和待评估的特征P₂的联合熵；

C＝H(P₁，P₂)-H(P₁)；

(3)从步骤(2)的六个特征中选取特征评估量E为最大的特征，记为P，用特征P与原始体数据的灰度特征组成一个二维传输函数，用该二维传输函数对原始体数据中的目标区域进行分离，得到原始体数据中位置在目标区域中的体素，分离过程如下：

(3-1)用特征P与体数据的灰度特征组成一个二维统计直方图，用归一化割(normalizedcut)方法对该直方图图像进行分割，得到多个二维统计直方图图像子区域，记该多个子区域组成的集合为Y₀；

(3-2)计算原始体数据中所有体素在二维统计直方图中的坐标，根据该坐标，从原始体数据中分别找到与每个子区域相对应的体素，利用光线投射方法，分别对与每个子区域相对应的体素进行体绘制，得到多个体绘制图像；

(3-3)从子区域集合Y₀中得到最终子区域集合O_final，过程如下：

(3-3-1)初始化时，最终子区域集合为空集，对Y₀中所有子区域进行从1到w排序，其中w是子区域集合Y₀中子区域的数量；

(3-3-2)从子区域集合Y₀中取出标号最小的子区域y₀，若与子区域y₀相对应的体绘制图像中只包含目标区域，则将y₀放入最终子区域集合中，并将y₀从Y₀中删除，将删除y₀后的Y₀记为Y₁；若与子区域y₀相对应的体绘制图像中只包含非目标区域，则将y₀从Y₀中删除，将删除y₀后的Y₀记为Y₁，若与子区域y₀相对应的体绘制图像中既包含目标区域又包含非目标区域，则用归一化割(normalized cut)方法对子区域y₀进行分割，得到两个新的子区域y₁和y₂，分别将y₁和y₂加入到子区域集合Y₀中，将y₀从Y₀中删除，得到新的子区域集合Y₁，Y₁中y₁标号为N_max+1，y₂标号为N_max+2，其中N_max为子区域集合Y₀中子区域标号的最大值；

重复上述步骤，若子区域集合Y变为空集，则记最终子区域集合为O_final，进行步骤(3-4)；若上述步骤重复次数超过30次，则将Y₃₀中剩余的所有子区域加入到最终子区域集合中，并记最终子区域集合为O_final，进行步骤(3-4)；

(3-4)根据步骤(3-2)得到原始体数据中所有体素的坐标，分别分离出坐标在最终子区域集合O_final中的每个子区域的所有体素，所有体素组成的集合记为V_r，按照体素在原始体数据中的坐标，将V_r中的所有体素加入到原始体数据中，将加入V_r中的所有体素后的原始体数据记为V_n，根据下式对V_n的灰度特征进行处理，得到V_n的新的灰度特征如下：

$f_{\mathrm{new}}\left(v_h\right)=\left\{\begin{array}{cc}{f}_{\mathrm{old}}\left(v_h\right),& v_h\∈V_r\\ 0,& v_h\∉V_r\end{array}\right.$

其中v_h为V_n中的任意一个体素，f_old(v_h)为V_n中v_h的原来的灰度特征，f_new(v_h)为V_n中v_h的新的灰度特征；

(3-5)对上述V_n中的体素进行判断，若V_n中的任意体素v_h满足v_h∈V_r，则进行步骤(4)，若V_n中存在一个体素v_h，满足则根据下式对步骤(1)的6种特征数据进行处理，得到与每种特征相对应的新的特征值：

$P_{\mathrm{feature}\_\mathrm{new}}\left(v_h\right)=\left\{\begin{array}{cc}{P}_{\mathrm{feature}\_\mathrm{old}}\left(v_h\right),& v_h\∈V_r\\ 0,& v_h\∉V_r\end{array}\right.$

其中v_h为V_n中任意一个体素，P_{feature_old}(v_h)表示v_h的原来特征值，P_{feature_new}(v_h)表示v_h的新特征值，返回步骤(2)；

(4)利用步骤(3)得到的体数据V_n和步骤(1)中的备份原始体数据V_copy，合成一个传输函数，利用该传输函数得到体绘制图像，具体过程如下：

(4-1)将步骤(3)得到的体数据V_n中的体素归为一个类，记为m₀；

(4-2)根据备份的原始体数据V_copy，将步骤(3)得到的体数据V_n中灰度特征变为0的体素的灰度特征值恢复为原来的灰度特征值，将恢复灰度特征值的体素组成一个集合，记为V_o，用V_o中所有体素的梯度特征和灰度特征组成一个二维统计直方图，采用归一化割(normalized cut)方法将此直方图分割成e个子区域，分别记为a₁、a₂…a₃；

(4-3)计算上述V_o中所有体素在梯度特征和灰度特征组成的二维统计直方图中的坐标，根据V_o中所有体素在该二维统计直方图中的坐标，将V_o划分为e个类，分别记为m₁、m₂...m_b...m₃，b＝1，2，...，e，分类方法为：对V_o中任意一个体素v_i，若v_i在该二维统计直方图中的坐标在子区域a_b内部，则将v_i归到类m_b中；

(4-4)根据显示需要，为体数据V_n中的体素的类m₀设置不透明度和颜色，为类m₁、m₂…m₃设置不透明度和互相不同的颜色，不透明度取值范围为[0，1]，颜色采用红绿蓝颜色模型，红、绿、蓝颜色分量的强度值范围为[0，255]，根据体数据中每个类的颜色和不透明度，采用光线投射算法得到最终的突出显示目标区域的体绘制图像。

本发明方法总体流程如下：

(1)计算原始体数据的灰度，梯度，分型特征，曲率，光滑度，线型特征和三维不变矩这7种特征；

(2)用距离，相关性和信息量三个特征评估方法对可变特征进行评估，并融合三个方法的评估结果选取合适的可变特征。选取出的特征与灰度特征组合构成传输函数；

(3)用normalized-cut方法分割二维直方图图像，逐步找到二维直方图中包含目标体素的区域，并分离出坐标在目标区域内部的所有体素。如果当前分离结果不够精细，就继续选择合适的特征对当前结果做进一步的分离，即用迭代的方式逐步得到精细的分离结果；

(4)把从原始体数据中分离出的每一个目标区域作为原始数据的一个类。原始体数据中经过分离后剩余的区域通过用normalized-cut方法分割直方图图像实现分类，即分割得到的每个图像子区域归为一类。在完成数据的分类后，对不同的类设置不同的颜色和不透明度，从而完成传输函数设计。

Claims (1)

1.一种用于突出显示体数据中目标区域的体绘制方法，其特征在于该方法包括以下步骤：

(1)分别计算体数据的梯度特征、曲率特征、分型特征、光滑度、线型特征和三维不变矩，具体过程如下：

(1-1)设体数据中任意一个体素的梯度特征为G，G＝|g|，G为g的模，g为体数据中体素的梯度向量，f为体数据灰度特征在x、y、z三个方向上的函数映射，灰度特征为体素的数值特征，为的转置，则体数据的梯度特征为体数据中所有体素的梯度特征按照体素在体数据中的位置排列得到的一个三维矩阵；

(1-2)设体数据中任意一个体素的曲率特征为C_m，κ₁为体素所在位置的主曲率，κ₂为体素所在位置的副曲率，T为几何张量Z的迹，F为几何张量Z的弗罗贝尼乌斯范数，几何张量其中H_e是体素的海森矩阵，P_m是体素所在位置的等值面切平面矩阵，P_m＝I-n_vn_v ^t，I为单位矩阵，n_v为体素所在位置的表面向量，n_v ^t为n_v的转置，则体数据的曲率特征为体数据中所有体素的曲率特征按照体素在体数据中的位置排列得到的三维矩阵；

(1-3)设体数据中任意一个体素的分型特征为D，对体数据中任意一个体素v_i，分型特征 $D=\frac{{\mathrm{\Σ}}_{r=1}^{k_d}(X_r-\stackrel{\‾}{X})(Y_r-\stackrel{\‾}{Y})}{{\mathrm{\Σ}}_{r=1}^{k_d}{(X_r-\stackrel{\‾}{X})}^2},$ 其中X_r＝logr， $\stackrel{\‾}{X}={\mathrm{\Σ}}_{r=1}^{k_d}\frac{X_r}{k_d},Y_r=\log u\left(\stackrel{\‾}{B}(v_i,r)\right),$ $\stackrel{\‾}{Y}={\mathrm{\Σ}}_{r=1}^{k_d}\frac{Y_r}{k_d},u\left(\stackrel{\‾}{B}(v_i,r)\right)={\mathrm{\Σ}}_{v_j\∈U(v_i,r)}f\left(v_j\right),$ U(v_i，r)为体素v_i的邻域，r为邻域半径，f(v_j)为体素v_j的灰度特征，v_j为体素v_i的邻域U(v_i，r)中任意体素，表示体素v_i的以r为半径的邻域中所有体素的灰度特征的和，k_d的取值为2，则体数据的分型特征为体数据中所有体素的分型特征按照体素在体数据中的位置排列得到的三维矩阵；

(1-4)设体数据中任意一个体素的光滑度为S，对体数据中任意一个体素v_i，光滑度其中σ²为体素v_i的邻域U(v_i，r)中所有体素的灰度特征方差，体素v_i的邻域半径r取2，则体数据的光滑度为体数据中所有体素的光滑度按照体素在体数据中的位置排列得到的三维矩阵；

(1-5)设体数据中任意一个体素的线型特征为F_flit，对体数据中任意一个体素v_i，线型特征其中为权值，取值为任意非负实数，r为体素v_i的邻域半径， $B_{k_f}=\left\{\begin{array}{c}\stackrel{\‾}{G_{k_f}},G_{k_f}>\mathrm{Mean}\left(G_{k_f}\right)\\ {0,G}_{k_f}\≤\mathrm{Mean}\left(G_{k_f}\right)\end{array}\right.,G_{k_f}={\mathrm{\Σ}}_{v_j\∈L_{k_f}}f\left(v_j\right),$ 其中为均一化后的结果， 为体素v_i的邻域内以邻域边界为两端端点且通过v_i的线段，线段总数为12r²+1，v_j为线段L_k上的一个体素，体素v_i的邻域半径r取1，则体数据的线型特征为体数据中所有体素的线型特征按照体素在体数据中的位置排列得到的三维矩阵；

(1-6)设体数据中任意一个体素的三维不变矩为T，三维不变矩T＝μ₂₀₀+μ₀₂₀+μ₀₀₂，其中μ₂₀₀、μ₀₂₀和μ₀₀₂分别为中心距，对任意一个体素v_i，中心距计算通式为：

${\mathrm{\μ}}_{\mathrm{pqt}}=\underset{v_j\∈U(v_i,r)}{\mathrm{\Σ}}{(x_1-\stackrel{\‾}{x_1})}^p{(x_2-\stackrel{\‾}{x_2})}^q{(x_3-\stackrel{\‾}{x_3})}^t\·f(x_1,x_2,x_3){\mathrm{dx}}_1{\mathrm{dx}}_2{\mathrm{dx}}_3$

其中v_j是体素v_i邻域U(v_i，r)内的任意一个体素，x₁，x₂，x₃分别为体素v_j在体数据中三个维度的坐标值，分别为体素v_i的邻域U(v_i，r)的灰度特征重心在体数据中三个维度的坐标值，f(x₁，x₂，x₃)为v_j的灰度特征，体素v_i的邻域半径r取4，则体数据的三维不变矩为体数据中所有体素的三维不变矩按照体素在体数据中的位置排列得到的三维矩阵；

(1-7)将原始体数据进行备份，备份的原始体数据记为V_copy；

(2)分别计算上述步骤(1)中得到的六个特征的特征评估量E，E＝ω₁L_n-ω₂R_n+ω₃C_n，其中ω₁，ω₂，ω₃为权值，取值为任意非负实数，L_n、R_n、C_n分别为对特征距离评估量L、特征相关性评估量R、特征价值评估量C进行均一化处理的结果，L、R、C分别为：

其中为用k均值聚类算法对特征P₁聚类后的第i_l个类的重心，P₁为上一次评估时选取的特征评估量E为最大的特征，若当前评估为第一次评估，则特征P₁为灰度特征，k₁为用k均值聚类算法对特征P₁聚类后得到的类的个数，为用k均值聚类算法对特征P₂聚类后的第j_l个类的重心，P₂为待评估的特征，k₂为用k均值聚类算法对待评估的特征P₂聚类后得到的类的个数；

R＝H(P₁)+H(P₂)-H(P₁，P₂)，其中H(P₁)为特征P₁的信息熵，H(P₂)为待评估特征P₂的信息熵，H(P₁，P₂)为特征P₁和待评估的特征P₂的联合熵；

C＝H(P₁，P₂)-H(P₁)；

(3)从步骤(2)的六个特征中选取特征评估量E为最大的特征，记为P，用特征P与原始体数据的灰度特征组成一个二维传输函数，用该二维传输函数对原始体数据中的目标区域进行分离，得到原始体数据中位置在目标区域中的体素，分离过程如下：

(3-1)用特征P与体数据的灰度特征组成一个二维统计直方图，用归一化割方法对该直方图图像进行分割，得到多个二维统计直方图图像子区域，记该多个子区域组成的集合为Y₀；

(3-2)计算原始体数据中所有体素在二维统计直方图中的坐标，根据该坐标，从原始体数据中分别找到与每个子区域相对应的体素，利用光线投射方法，分别对与每个子区域相对应的体素进行体绘制，得到多个体绘制图像；

(3-3)从子区域集合Y₀中得到最终子区域集合O_final，过程如下：

(3-3-1)初始化时，最终子区域集合为空集，对Y₀中所有子区域进行从1到w排序，其中w是子区域集合Y₀中子区域的数量；

(3-3-2)从子区域集合Y₀中取出标号最小的子区域y₀，若与子区域y₀相对应的体绘制图像中只包含目标区域，则将y₀放入最终子区域集合中，并将y₀从Y₀中删除，将删除y₀后的Y₀记为Y₁；若与子区域y₀相对应的体绘制图像中只包含非目标区域，则将y₀从Y₀中删除，将删除y₀后的Y₀记为Y₁，若与子区域y₀相对应的体绘制图像中既包含目标区域又包含非目标区域，则用归一化割方法对子区域y₀进行分割，得到两个新的子区域y₁和y₂，分别将y₁和y₂加入到子区域集合Y₀中，将y₀从Y₀中删除，得到新的子区域集合Y₁，Y₁中y₁标号为N_max+1，y₂标号为N_max+2，其中N_max为子区域集合Y₀中子区域标号的最大值；

重复上述步骤，若子区域集合Y变为空集，则记最终子区域集合为O_final，进行步骤(3-4)；若上述步骤重复次数超过30次，则将Y₃₀中剩余的所有子区域加入到最终子区域集合中，并记最终子区域集合为O_final，进行步骤(3-4)；

(3-4)根据步骤(3-2)得到原始体数据中所有体素的坐标，分别分离出坐标在最终子区域集合O_final中的每个子区域的所有体素，所有体素组成的集合记为V_r，按照体素在原始体数据中的坐标，将V_r中的所有体素加入到原始体数据中，将加入V_r中的所有体素后的原始体数据记为V_n，根据下式对V_n的灰度特征进行处理，得到V_n的新的灰度特征如下：

$f_{\mathrm{new}}\left(v_h\right)=\left\{\begin{array}{cc}{f}_{\mathrm{old}}\left(v_h\right),& v_h\∈V_r\\ 0,& v_h\∉V_r\end{array}\right.$

其中v_h为V_n中的任意一个体素，f_old(v_h)为V_n中v_h的原来的灰度特征，f_new(v_h)为V_n中v_h的新的灰度特征；

(3-5)对上述V_n中的体素进行判断，若V_n中的任意体素v_h满足v_h∈V_r，则进行步骤(4)，若V_n中存在一个体素v_h，满足则根据下式对步骤(1)的6种特征数据进行处理，得到与每种特征相对应的新的特征值：

$P_{\mathrm{feature}\_\mathrm{new}}\left(v_h\right)=\left\{\begin{array}{cc}{P}_{\mathrm{feature}\_\mathrm{old}}\left(v_h\right),& v_h\∈V_r\\ 0,& v_h\∉V_r\end{array}\right.$

其中v_h为V_n中任意一个体素，P_{feature_old}(v_h)表示v_h的原来特征值，P_{featutr_new}(v_h)表示v_h的新特征值，返回步骤(2)；

(4)利用步骤(3)得到的体数据V_n和步骤(1)中的备份原始体数据V_copy，合成一个传输函数，利用该传输函数得到体绘制图像，具体过程如下：

(4-1)将步骤(3)得到的体数据V_n中的体素归为一个类，记为m₀；

(4-2)根据备份的原始体数据V_copy，将步骤(3)得到的体数据V_n中灰度特征变为0的体素的灰度特征值恢复为原来的灰度特征值，将恢复灰度特征值的体素组成一个集合，记为V_o，用V_o中所有体素的梯度特征和灰度特征组成一个二维统计直方图，采用归一化割方法将此直方图分割成e个子区域，分别记为a₁、a₂…a_e；

(4-3)计算上述V_o中所有体素在梯度特征和灰度特征组成的二维统计直方图中的坐标，根据V_o中所有体素在该二维统计直方图中的坐标，将V_o划分为e个类，分别记为m₁、m₂...m_b...m_e，b＝1，2，...，e，分类方法为：对V_o中任意一个体素v_i，若v_i在该二维统计直方图中的坐标在子区域a_b内部，则将v_i归到类m_b中；

(4-4)根据显示需要，为体数据V_n中的体素的类m₀设置不透明度和颜色，为类m₁、m₂…m_e设置不透明度和互相不同的颜色，不透明度取值范围为[0，1]，颜色采用红绿蓝颜色模型，红、绿、蓝颜色分量的强度值范围为[0，255]，根据体数据中每个类的颜色和不透明度，采用光线投射算法得到最终的突出显示目标区域的体绘制图像。