CN100589126C - 利用阻光度传递函数的ct图像体素成像方法 - Google Patents

Abstract

本发明提供了一种可以实现三维体数据中的多种物质分类的阻光度传递函数的构造方法。首先对试件进行CT扫描，获得二维CT图像序列，由二维CT图像序列中得到体数据，建立体数据直方图并对直方图进行分段，在每段中求取一个阈值作为阻光度传递函数的分段点，确定阻光度传递函数，确定物体中待显示的部分并依据阻光度传递函数对每个体素进行阻光度赋值，利用错切—变形算法进行体素成像，最后显示物体图像。

Description

利用阻光度传递函数的CT图像体素成像方法

(一)技术领域

本发明涉及一种CT图像体素成像方法，具体涉及一种利用阻光度传递函数的CT图像体素成像方法。

(二)背景技术

在CT图像体素成像过程中，为了在最终的可视化图像中以不同的颜色正确地表示出共存的多种物质的分布，就需要对体数据进行分类，找出数据与不同的物质之间的对应关系。物质的分类是通过阻光度传递函数实现的，具体体现在需要对待显示的数据集合中的每个体素指定一个阻光度值，它表示物体的透明程度，其值介于0.0到1.0之间。为体素指定阻光度是一个非常复杂的运算，是通过阻光度传递函数来实现的。根据体数据的分类结果将某个结构中的点赋予高的阻光度从而使这部分结构可见，相反，让不感兴趣的部分变成透明。

传统的KSW熵(Kapur J N，Sahoo P K，Wong AKC.A new method for gray-level picturethresholding using the entropy of the histogram.Computer Vision，Graphics，and Image Processing，1985，29(3)：273～285P)和全局熵算法(N.R.Pal.IEEE.SMC-21，No.5，1991.)是对二维图像进行单阈值分割，无法处理CT三维体数据，更不能用于多种物质的分类。

(三)发明内容

本发明的目的在于提供一种可以实现CT三维体数据中的多种物质分类的利用阻光度传递函数的CT图像体素成像方法。

本发明的目的是这样实现的：

首先对试件进行CT扫描，获得二维CT图像序列，由二维CT图像序列中得到体数据，建立体数据直方图并对直方图进行分段，在每段中求取一个阈值作为阻光度传递函数的分段点，确定阻光度传递函数，确定物体中待显示的部分并依据阻光度传递函数对每个体素进行阻光度赋值，利用错切-变形算法进行体素成像，最后显示物体图像；

所述的建立体数据直方图是建立CT三维体数据直方图，其方法为：设三维体数据的尺寸为C×N×G，f(x_l，j，k)表示体素x_l，j，k的灰度值，将三维体数据的直方图h(u)定义为：

$h\left(u\right)=\frac{\underset{i=1}{\overset{C}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\underset{k=1}{\overset{G}{\mathrm{\Σ}}}g\left(u\right)}{C\×N\×G};$

所述的直方图进行分段是对CT三维体数据直方图h(u)分段，其分段方法为：

首先计算h(u)的非零端点f_min和f_max，再将f_min～f_max平均分为n段，即

然后分别计算这n段的累积直方图

若则将该段进行扩展，即使得

具体步骤为：

步骤(1)：

按灰度将直方图h(u)，其中u＝0，1，...，255，平均分成n段，端点为f₁，f₂，...，f_n-1；

步骤(2)：

计算每段的累计直方图h_m，

$h_m=\underset{u=f_m}{\overset{f_{m+1}}{\mathrm{\Σ}}}h\left(u\right)$ m＝1，2，...，n-2

步骤(3)：

比较h_m与1/n+2的大小，如果满足下式，则转向步骤(5)，否则转向步骤(4)；

h_m≥1/n+2  m＝1，2，...，n-2

步骤(4)：

进而按下式进行分段扩展，首先判断是否是在首段或左段的累计直方图是否小于右段的累计直方图，若是，则向右扩展；否则，判断是否是在末段或右段的累计直方图是否小于等于左段的累计直方图，若是，则向左扩展；

f_m+1＝f_m+1+1  m＝1或h_m-1＜h_m+1

f_m-1＝f_m-1-1  m＝n-2或h_m-1≥h1_m+1

然后转向步骤(2)；

步骤(5)：

确定端点为f₁，f₂，...，f_n-1；

用上述分段方法将CT三维体数据直方图分成n段后，在每段中求取一个阈值作为阻光度传递函数的分段点，分别为t₁，t₂，...t_n；

设f(x_i，j，k)、

α(x_i，j，k)分别为体素x_i，j，k的灰度值、梯度值和阻光度值；假设在体数据中有n种待显示的物质，第v种待显示的物质的灰度范围是t_v～t_v+1，其中t_v＜t_v+1，v＝1，2，…，n；与t_v对应的阻光度值为α_v；

所述的确定阻光度传递函数选择分段函数，模型为：

$\▿f\left(x_{i,j,k}\right)=\left[\begin{array}{c}\frac{1}{2}(f\left(x_{i+1,j,k}\right)-f\left(x_{i-1,j,k}\right)),\frac{1}{2}(f\left(x_{i,j+1,k}\right)-f\left(x_{i,j-1,k}\right))\\ ,\frac{1}{2}(f\left(x_{i,j,k+1}\right)-f\left(x_{i,j,k-1}\right))\end{array}\right]$

α(x_i，j，k)为f(x_i，j，k)和

的函数，其中n＝4。

本发明还有这样一些技术特征：

1、所述的在每段中求取一个阈值作为阻光度传递函数的分段点采用基于CT三维体数据直方图熵的多阈值分类方法。

2、所述的基于CT三维体数据直方图熵的多阈值分类方法中测量CT三维体数据直方图熵的方法为：

设体数据的各个分段的直方图为h_m(i)，

(1)KSW熵方法

假设试验结果由S个基本事件组成，第r个事件出现的概率为p_r，则整个试验结果获取的信息量H为：

$H=E\left(\mathrm{\ΔI}\right)=-\underset{r=1}{\overset{S}{\mathrm{\Σ}}}{p}_r\ln p_r$

三维体数据的总熵如下式所示，求使得H(t)最大的t值即为三维体数据的阈值；

$H\left(t\right)=\ln P_t(1-P_t)+\frac{H_t}{P_t}+\frac{H_{\max }-H_t}{1-P_t}$

式中：

$P_t=\underset{i=f_{\min }}{\overset{t}{\mathrm{\Σ}}}{h}_m\left(i\right)$

$H_{\max }=-\underset{i=f_{\min }}{\overset{f_{\max }}{\mathrm{\Σ}}}{h}_m\left(i\right)\ln h\left(i\right)$

$H_t=-\underset{i=f_{\min }}{\overset{t}{\mathrm{\Σ}}}{h}_m\left(i\right)\ln h\left(i\right)$

(2)全局熵方法

全局熵的定义为：

$H=\underset{i}{\mathrm{\Σ}}{h}_m\left(i\right)e^{1-h\left(i\right)}$

其中的三维体数据的总熵如下式所示，求使得H(t)最小的t值即为体数据的阈值；

$H\left(t\right)=\underset{i=t+1}{\overset{f_{\max }}{\mathrm{\Σ}}}[h_m\left(i\right)/(P-P_A)]\×e^{1-[h_m\left(i\right)/(P-P_A)]}+\underset{i=f_{\min }}{\overset{t}{\mathrm{\Σ}}}(h_m\left(i\right)/P_A)\×e^{1-h_m\left(i\right)/P_A}$

式中：

$P_A=\underset{i=f_{\min }}{\overset{t}{\mathrm{\Σ}}}{h}_m\left(i\right)$

$P=\underset{i=f_{\min }}{\overset{f_{\max }}{\mathrm{\Σ}}}{h}_m\left(i\right).$

本发明的工作原理是：首先从二维CT图像序列中获得体数据，然后绘制体数据的直方图，对体数据直方图进行分段，在每一段中依据CT三维体数据直方图的KSW熵和全局熵构造阻光度传递函数，对体素进行阻光度赋值，利用错切-变形体绘制算法，显示物体中任意感兴趣部分的三维构型。

在目前阻光度传递函数的设计仍然较为困难。本发明所采用的阻光度传递函数为分段函数，其中分段点的确定是非常关键的，并对成像结果有很大的影响。本发明主要解决阻光度传递函数分段点的确定问题。本发明对传统的KSW熵和全局熵算法进行了改进，提出了两种阻光度传递函数的构造算法，实现了CT三维体数据中的多种物质分类。

(四)附图说明

图1为突出表面强度和物体实体感的阻光度函数示意图；

图2为三维体数据直方图分段方法流程；

图3为基于KSW熵的阻光度传递函数构造的三维体数据的直方图；

图4为基于KSW熵的阻光度传递函数构造的阻光度函数示意图；

图5为基于全局熵的阻光度传递函数构造的体数据的直方图；

图6为基于全局熵的阻光度传递函数构造的阻光度函数示意图。

(五)具体实施方式

下面结合附图和具体实施例对本发明作进一步的详细说明：

本发明的方法为：对试件进行CT扫描，获得二维CT图像序列，从二维CT图像序列中得到体数据。确定阻光度传递函数，建立体数据直方图并对直方图进行分段，在每段中求取一个阈值作为阻光度传递函数的分段点。确定物体中待显示的部分并依据阻光度传递函数对每个体素进行阻光度赋值，利用错切-变形算法显示物体中感兴趣的部分。

1、阻光度传递函数的构造

为了灵活地显示物体中任意感兴趣的部分，阻光度传递函数通常为分段函数。本发明采用的阻光度传递函数模型如下，该函数能突出表面强度和物体实体感。

设f(x_i，j，k)、

α(x_i，j，k)分别为体素x_i，j，k的灰度值、梯度值和阻光度值；假设在体数据中有n种待显示的物质，第v种待显示的物质的灰度范围是t_v～t_v+1，其中t_v＜t_v+1，v＝1，2，…，n；与t_v对应的阻光度值为α_v。

$\▿f\left(x_{i,j,k}\right)=\left[\begin{array}{c}\frac{1}{2}(f\left(x_{i+1,j,k}\right)-f\left(x_{i-1,j,k}\right)),\frac{1}{2}(f\left(x_{i,j+1,k}\right)-f\left(x_{i,j-1,k}\right))\\ ,\frac{1}{2}(f\left(x_{i,j,k+1}\right)-f\left(x_{i,j,k-1}\right))\end{array}\right]---\left(2\right)$

α(x_i，j，k)为f(z_i，j，k)和

的函数，三者的关系可以用图1表示，其中n＝4。

2、传递函数分段点的确定

在本发明中提出基于三维体数据直方图熵的多阈值分类方法，将三维体数据的直方图按其灰度大小分为n段，每一段中确定一个阈值，即可得到n个阈值，即阻光度传递函数中的n个分段点，从而将体数据分成n+1个类。

2.1、体数据直方图的计算

设体数据的尺寸为C×N×G，f(x_l，j，k)表示体素x_l，j，k的灰度值，本发明将三维体数据的直方图h(u)定义为：

$h\left(u\right)=\frac{\underset{l=1}{\overset{C}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\underset{k=1}{\overset{G}{\mathrm{\Σ}}}g\left(u\right)}{C\×N\×G}---\left(4\right)$

2.2、三维体数据直方图的分段

本发明中提出的h(u)分段方法如下。

首先计算h(u)的非零端点f_min和f_max，再将f_min～f_max平均分为n段，即

然后分别计算这n段的累积直方图

若

则为了保证每一个分段中的信息量，应该将该段进行扩展，即使得且最终保证每段的灰度范围不重叠。该方法的具体流程如图2所示。

步骤1：

按灰度将直方图h(u)，(u＝0，1，...，255)平均分成n段，端点为f₁，f₂，...，f_n-1。

步骤2：

计算每段的累计直方图h_m，

$h_m=\underset{u=f_m}{\overset{f_{m+1}}{\mathrm{\Σ}}}h\left(u\right)$ m＝1，2，...，n-2    (5)

步骤3：

比较h_m与1/n+2的大小，如果满足式(6)，则转向步骤5，否则转向步骤4。

h_m≥1/n+2  m＝1，2，...，n-2    (6)

步骤4：

进而按(7)式进行分段扩展。即判断是否是在首段或左段的累计直方图是否小于右段的累计直方图，若是，则向右扩展；否则，是否是在末段或右段的累计直方图是否小于等于左段的累计直方图，若是，则向左扩展。

f_m+1＝f_m+1+1  m＝1或h_m-1＜h_m+1

                                     (7)

f_m-1＝f_m-1-1  m＝n-2或h_m-1≥h_m+1

然后转向步骤2。

步骤5：

确定端点为f₁，f₂，...，f_n-1。

2.3、各分段中阈值的求取

信息论中用熵作不确定性的度量，或者说对随机试验结果的无知度量。如果一个试验结果由几个独立的基本事件组成，各基本事件出现的概率为已知，那么试验后所得到的信息量与人们已知的各基本事件的出现概率成反比。对于确定要发生的事件，没有获得任何信息；而对于最不确定是否发生的事件，则获得最多的信息。

目前还没有人将信息论中熵的概念应用于体数据的分类。在我们的算法中，需要测量体数据灰度直方图的熵，并由此自动找出最佳阈值对体数据进行分类。设体数据的各个分段的直方图为h_m(i)，最大灰度值为f_max，最小灰度值为f_min，则h_m(i)表明三维体数据中灰度值为i(f_min＜i＜f_max)的体素出现的概率，本发明中提出了两种测量三维体数据直方图熵的方法。

(1)KSW熵方法

假设试验结果由S个基本事件组成，第r个事件出现的概率为p_r，则整个试验结果获取的信息量H为：

$H=E\left(\mathrm{\ΔI}\right)=-\underset{r=1}{\overset{n}{\mathrm{\Σ}}}{p}_r\ln p_r---\left(8\right)$

这就是信息论中的熵。

本发明中的三维体数据的总熵如式(9)所示，求使得H(t)最大的t值即为三维体数据的阈值。

$H\left(t\right)=\ln P_t(1-P_t)+\frac{H_t}{P_t}+\frac{H_{\max }-H_t}{1-P_t}---\left(9\right)$

式中：

$P_t=\underset{i=f_{\min }}{\overset{t}{\mathrm{\Σ}}}{h}_m\left(i\right)---\left(10\right)$

$H_{\max }=-\underset{i=f_{\min }}{\overset{f_{\max }}{\mathrm{\Σ}}}{h}_m\left(i\right)\ln h\left(i\right)---\left(11\right)$

$H_t=-\underset{i=f_{\min }}{\overset{t}{\mathrm{\Σ}}}{h}_m\left(i\right)\ln h\left(i\right)---\left(12\right)$

图3-4为三维体数据的直方图及其基于KSW熵的阻光度传递函数。以大小为512×512×15的体数据为例，其中图3为体数据的直方图，由于不同灰度值所对应的象素的比例相差较大，为了使直方图更加清晰，对不同灰度段采用不同的显示比例，并将比例尺图中标出，图4为基于KSW熵的阻光度传递函数。

(2)全局熵方法

全局熵的定义为

$H=\underset{i}{\mathrm{\Σ}}{h}_m\left(i\right)e^{1-h\left(i\right)}---\left(13\right)$

本发明中的三维体数据的总熵如式(14)所示，求使得H(t)最小的t值即为体数据的阈值。

$H\left(t\right)=\underset{i=t+1}{\overset{f_{\max }}{\mathrm{\Σ}}}[h_m\left(i\right)/(P-P_A)]\×e^{1-[h_m\left(i\right)/(P-P_A)]}+\underset{i=f_{\min }}{\overset{t}{\mathrm{\Σ}}}(h_m\left(i\right)/P_A)\×e^{1-h_m\left(i\right)/P_A}---\left(14\right)$

式中：

$P_A=\underset{i=f_{\min }}{\overset{t}{\mathrm{\Σ}}}{h}_m\left(i\right)---\left(15\right)$

$P=\underset{i=f_{\min }}{\overset{f_{\max }}{\mathrm{\Σ}}}{h}_m\left(i\right)---\left(16\right)$

图5-6为三维体数据的直方图及其基于全局熵的阻光度传递函数。以大小为512×512×253的体数据为例，其中图5为体数据的直方图，由于不同灰度值所对应的象素的比例相差较大，为了使直方图更加清晰，对不同灰度段采用不同的显示比例，并将比例尺图中标出，图6为基于全局熵的阻光度传递函数。

Claims (3)

1、一种利用阻光度传递函数的CT图像体素成像方法，首先对试件进行CT扫描，获得二维CT图像序列，由二维CT图像序列中得到CT体数据，建立CT体数据直方图并对直方图进行分段，在每段中求取一个阈值作为阻光度传递函数的分段点，确定阻光度传递函数，确定物体中待显示的部分并依据阻光度传递函数对每个体素进行阻光度赋值，利用错切一变形算法进行体素成像，最后显示物体图像；其特征在于：

所述的建立CT体数据直方图是建立三维体数据直方图，其方法为：设CT三维体数据的尺寸为C×N×G，f(x_i，j，k)表示体素x_i，j，k的灰度值，将CT三维体数据的直方图h(u)定义为：

$h\left(u\right)=\frac{\underset{i=1}{\overset{C}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\underset{k=1}{\overset{G}{\mathrm{\Σ}}}g\left(u\right)}{C\×N\×G};$

所述的直方图进行分段是对CT三维体数据直方图h(u)分段，其分段方法为：

首先计算h(u)的非零端点f_min和f_max，再将f_min～f_max平均分为n段，即 $f_{\min }~\frac{f_{\max }+(n-1)f_{\min }}{n},\frac{f_{\max }+(n-1)f_{\min }}{n}~\frac{2\×f_{\max }+(n-2)f_{\min }}{n}...\frac{(n-1)f_{\max }+f_{\min }}{n}~f_{\max };$ 然后分别计算这n段的累积直方图

若 $\underset{u}{\mathrm{\&Sigma;}}h\left(u\right)<\frac{1}{n+2},$ 则将该段进行扩展，即使得 $\underset{u}{\mathrm{\&Sigma;}}h\left(u\right)\&GreaterEqual;\frac{1}{n+2};$ 具体步骤为：

步骤(1)：

按灰度将直方图h(u)，其中u＝0，1，..，255，平均分成n段，端点为f₁，f₂，..，f_n-1；

步骤(2)：

计算每段的累计直方图h_m，

$h_m=\underset{u=f_m}{\overset{f_{m+1}}{\mathrm{\&Sigma;}}}h\left(u\right)$ m＝1，2，..，n-2

步骤(3)：

比较h_m与1/n+2的大小，如果满足下式，则转向步骤(5)，否则转向步骤(4)；

$h_m\&GreaterEqual;\frac{1}{n+2}$ m＝1，2，..，n-2

步骤(4)：

进而按下式进行分段扩展，首先判断是否是在首段或左段的累计直方图是否小于右段的累计直方图，若是，则向右扩展；否则，判断是否是在末段或右段的累计直方图是否小于等于左段的累计直方图，若是，则向左扩展；

f_m+1＝f_m+1+1  m＝1或f_m-1＜h_m+1

f_m-1＝f_m-1-1  m＝n-2或h_m-1≥h_m+1

然后转向步骤(2)；

步骤(5)：

确定端点为f₁，f₂，...，f_n-1；

用上述分段方法将CT三维体数据直方图分成n段后，在每段中求取一个阈值作为阻光度传递函数的分段点，分别为t₁，t₂，...t_n；

设f(x_i，j，k)、

α(x_i，j，k)分别为体素x_i，j，k的灰度值、梯度值和阻光度值；假设在体数据中有n种待显示的物质，第v种待显示的物质的灰度范围是t_v～t_v+1，其中t_v＜t_v+1，v＝1，2，…，n；与t_v对应的阻光度值为α_v；

所述的确定阻光度传递函数选择分段函数，模型为：

$\&dtri;f\left(x_{i,j,k}\right)=\left[\begin{array}{c}\frac{1}{2}(f\left(x_{i+1,j,k}\right)-f\left(x_{i-1,j,k}\right)),\frac{1}{2}(f\left(x_{i,j+1,k}\right)-f\left(x_{i,j-1,k}\right))\\ ,\frac{1}{2}(f\left(x_{i,j,k+1}\right)-f\left(x_{i,j,k-1}\right))\end{array}\right]$

α(x_i，j，k)为f(x_i，j，k)和

的函数，其中n＝4。

2、根据权利要求1所述的利用阻光度传递函数的CT图像体素成像方法，其特征在于所述的在每段中求取一个阈值作为阻光度传递函数的分段点采用基于CT三维体数据直方图熵的多阈值分类方法。

3、根据权利要求2所述的利用阻光度传递函数的CT图像体素成像方法，其特征在于所述的基于CT三维体数据直方图熵的多阈值分类方法中测量CT三维体数据直方图熵的方法为：

设体数据的各个分段的直方图为h_m(i)，

(1)KSW熵方法

假设试验结果由S个基本事件组成，第r个事件出现的概率为p_r，则整个试验结果获取的信息量H为：

$H=E\left(\mathrm{\&Delta;I}\right)=-\underset{r=1}{\overset{S}{\mathrm{\&Sigma;}}}{p}_r\ln p_r$

三维体数据的总熵如下式所示，求使得H(t)最大的t值即为三维体数据的阈值；

$H\left(t\right)=\ln P_t(1-P_t)+\frac{H_t}{P_t}+\frac{H_{\max }-H_t}{1-P_t}$

式中：

$P_t=\underset{i=f_{\min }}{\overset{t}{\mathrm{\&Sigma;}}}{h}_m\left(i\right)$

$H_{\max }=-\underset{i=f_{\min }}{\overset{f_{\max }}{\mathrm{\&Sigma;}}}{h}_m\left(i\right)\ln h\left(i\right)$

$H_t=-\underset{i=f_{\min }}{\overset{t}{\mathrm{\&Sigma;}}}{h}_m\left(i\right)\ln h\left(i\right)$

(2)全局熵方法

全局熵的定义为：

$H=\underset{i}{\mathrm{\&Sigma;}}{h}_m\left(i\right)e^{1-h\left(i\right)}$

其中的三维体数据的总熵如下式所示，求使得H(t)最小的t值即为体数据的阈值；

$H\left(t\right)=\underset{i=t+1}{\overset{f_{\max }}{\mathrm{\&Sigma;}}}[h_m\left(i\right)/(P-P_A)]\&times;e^{1-[h_m\left(i\right)/(P-P_A)]}+\underset{i=f_{\min }}{\overset{t}{\mathrm{\&Sigma;}}}(h_m\left(i\right)/P_A)\&times;e^{1-h_m\left(i\right)/P_A}$

式中：

$P_A=\underset{i=f_{\min }}{\overset{t}{\mathrm{\&Sigma;}}}{h}_m\left(i\right)$

$P=\underset{i=f_{\min }}{\overset{f_{\max }}{\mathrm{\&Sigma;}}}{h}_m\left(i\right).$

Images