CN102749053A - 基于三维可视化和蒙特卡罗方法的体积测量方法 - Google Patents

Abstract

一种计算机信息处理技术领域的基于三维可视化和蒙特卡罗方法的体积测量方法，通过对被测物体的序列图像进行三维重建，并将重建的结果置于创建得到的虚拟盒内部，然后在虚拟盒内部产生若干个随机点，最后统计三维重建后的被测物体内部的随机点的个数，并根据虚拟盒的体积和随机点的总数按比例求取被测物体的体积大小。本发明实现了对不规则物体体积的测量，且测量对象不需要进行分割、分块等处理，结果精度相对较高。

Description

基于三维可视化和蒙特卡罗方法的体积测量方法

技术领域

本发明涉及的是一种计算机信息处理技术领域的方法，具体是一种虚拟环境下主要针对不规则、不封闭物体的基于三维可视化和蒙特卡罗方法的体积测量方法。 

背景技术

体积大小的测量在很多领域中都有着重要的意义，近年来，由于二维信息的局限性，科学计算可视化技术得到迅速发展，由此引起的可视化技术的应用越来越引起人们的重视。可视化工具包VTK是美国kitware公司开发的，主要针对的是二维和三维图形图像处理和可视化，它包含了可视化的多种方法，支持跨平台使用。 

目前对于体积的测量的方法主要是针对规则的物体，一般有现成的体积计算公式可用，对不规则物体体积的测量也有现成的方法。现有的测量不规则物体体积方法主要有：排水法、剖分法、经线法、二维超声法和注水法等。排水法是金标准，针对的对象是离体的组织或器官，不适宜人体内部器官或组织的测量；剖分法是最早提出的测量体积的方法，它将三角形网格先剖分成四面体然后求和，由于该方法复杂和低效，现在已经很少使用；经线法主要是针对离体的器官和病灶；二维超声是通过测出器官的长度、宽度和厚度，再按照标准的椭圆公式进行计算得出器官的体积，该方法简单方便，但是精度不高，不适用于相对精确的测量；注水法针对的对象是完全封闭的对象，虽然能相对精确的求得体积的大小，但在使用上存在一定的局限性。对求取不规则不封闭物体的方法却很少。 

蒙特卡罗方法是一种以概率统计理论为指导的数值计算方法，属于计算数学的一个分支，是利用随机数来解决问题的方法。具体的思路是：给定一个多边形网格，首先把它转化为三角形网格，然后在网格的包围盒外面放置一个简单的多边形容器，在容器内产生大量的随机点，然后判断有多少随机点在待测多边形的网格中，根据在多边形中的随机点的个数和总共的随机点数的比率，再结合多边形的容器的体积，就能得到待测多边形网格在容器内部的部分的体积。 

用蒙特卡罗方法求解时，最简单的情况是模拟一个事件A，它发生的概率是p。设置一个随机变量ξ，若事件A发生则ξ值为1，否则ξ值为0。那么随机变量ξ的数学期望E(ξ)＝1*p+0*q＝p，其中q＝1-p，即一次实验中事件A发生的概率，ξ的方差为E(ξ-E(ξ))²＝p-p²＝pq。假设在N次实验中事件A出现的次数是v，该频数v也是一个随机变量，其数学期望为E(v)＝Np，方差为σ²(v)＝Npq。令 

表示观察频率， 

按照 大数定理，当N值充分大时， 

$\stackrel{\‾}{p}=\frac{v}{N}\≈E\left(\mathrm{\ξ}\right)=p$ ........................................(式1) 

成立的概率等于1，由上述模型得到的频率 

和p近似相等，说明了频率是收敛的，并且收敛于此概率，可利用样本方差： 

${\mathrm{\σ}}^2\left(p\right)=\frac{\stackrel{\‾}{p}(1-\stackrel{\‾}{p})}{N-1}$ ........................................(式2) 

作为理论方差σ²(p)的相似估计。Monte Carlo方法的收敛性分析： 

蒙特卡罗通常采用某个随机变量X子样的算术平均值作为所求问题I的近似解： 

$\stackrel{\‾}{x_n}=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{x}_n$ ................................................(式3) 

其中：x₁，x₂，x₃，x₄，......，x_n是X子样函数，由柯尔莫哥罗夫加强大数定理得，当E(X)＝I时，算术平均值 

将以概率1收敛于I，按照中心极限定理，对于任何λ_α＞0时有 

$P\left(\right|\stackrel{\&OverBar;}{x_n}-I|<\frac{{\mathrm{\&lambda;}}_{\mathrm{\&alpha;}}\mathrm{\&sigma;}}{\sqrt{N}})\&ap;\frac{2}{\sqrt{2\mathrm{\&pi;}}}{\&Integral;}_0^{{\mathrm{\&lambda;}}_{\mathrm{\&alpha;}}}{e}^{-\frac{1}{2}{t}^2}\mathrm{dt}=1-\mathrm{\&alpha;}$ ..........(式4) 

表明，这个不等式 

$|\stackrel{\&OverBar;}{x_n}-I|<\frac{{\mathrm{\&lambda;}}_{\mathrm{\&alpha;}}\mathrm{\&sigma;}}{\sqrt{N}}$ ...........................................(式5) 

近似的和1-α的概率相等。在通常情况下，当α很小时，如α＝0.01时，α称为显著水平，1-α称之为置信水平，α称为随机变量X的标准差。当σ≠0时，蒙特卡罗方法的误差ε为： 

$\mathrm{\&epsiv;}=\frac{{\mathrm{\&lambda;}}_{\mathrm{\&alpha;}}\mathrm{\&sigma;}}{\sqrt{N}}$ ......................................................(式6) 

在上式中正态差λ_α与α是一一对应的关系，其对应关系可以用下面N(0，1)积分公式进行表示： 

$=\frac{1}{\sqrt{2\mathrm{\&pi;}}}{\&Integral;}_{-\&infin;}^{{\mathrm{\&lambda;}}_{\mathrm{\&alpha;}}}{e}^{-\frac{1}{2}{t}^2}\mathrm{dt}=1-\frac{\mathrm{\&alpha;}}{2}$ ..........................................(式7) 

分析公式 

可以得出该方法的误差精度ε是由σ和 

决定的。在σ不变的情况下，要提高精度一位数则要增加100倍的工作量；另外如果将σ减小十倍，则减少100倍的工作量。在实际中不能为了提高精度人为的随意改变参数，因为既要考虑精度也要考虑工作量和现实的条件。 

经过对现有技术的检索发现，浙江大学周海燕硕士论文《基于三维重建的肺部肿瘤的精确 测量》，该文献首先分析了肺部肿瘤CT图像的特点、三维重建各个步骤所涉及的方法以及肿瘤任意两点的直径测量所涉及的方法，以此为基础，对肺部肿瘤三维重建和测量进行了具体的方案研究与验证，实现了基于三维重建的肺部肿瘤的精确测量。包括：分析肺部肿瘤CT图像的特点，分析肺部肿瘤三维重建的各个步骤，包括：CT图像预处理、肺部肿瘤的分割、绘制与显示等。分析肿瘤直径测量所涉及的三维交互，包括：三维交互拾取和三维交互旋转。基于ITK与VTK等开发库，实现肺部肿瘤的三维重建，并在三维空间内进行任意两点的直径测量。通过临床医生验证，该方法只实现了在三维可视化的前提下的空间测量，比如空间距离的测量，没有涉及到体积大小的测量。 

中国专利文献号CN101756710A，公开日2010-06-30，记载了一种“颅内血肿体积的测量方法”，该方法1)CT图像的第一步分割中，颅骨及颅骨外非脑组织的去除利用了阈值法、区域生长及图像逻辑运算，由于颅骨在CT图像中与其它组织在灰度上差异很大，因此采用经验阈值很容易得到颅骨边缘，不需要复杂耗时的检测方法，区域生长采用四邻域生长，遍布整个脑组织与血肿，最后对图像进行“异或”和“乘”逻辑运算，得到只有脑组织与血肿的CT图像；2)CT图像的第二步分割中，对于脑组织与血肿的分割，考虑到CT图像可避免噪声干燥，采用二维熵阈值化的方法，用像素灰度及其邻域灰度为参数来分割图像；3)采用改进的遗传方法来搜寻最优分割阈值使二维熵准则函数最大；4)采用了“自适应的大变异算子”，通过判断适应度函数平均值与最大值比值，最终决定采用变异概率和交叉概率。但该技术的缺点在于：该方法过程复杂，需要进行两次分割处理，类似于文中提到的面积法求体积，工作量很大。文中提到自动分割的概念，也未对其实现过程详细介绍，最重要的一点就是该发明中未体现出分割后的图像数据及测量数据，缺乏一定的认可性。 

中国专利文献号CN101601585，公开日2009-12-16，记载了一种“基于CT增强扫描技术的肝脏分段体积测量方法”，该技术包括以下步骤：1)图像分割：采用图像处理方法对CT扫描后获得的数据进行处理，获取肝脏器官的表面轮廓以及肝脏内部的管道区域；2)三维细化：获取肝脏管道中心线；3)管道树分级：将管道中心线处理成为一(棵)树形结构，即结点的数目等于边的数目加一；4)管道投影：将肝脏数据投影到一个投影平面，将三维的肝脏分段转换为一个二维的分类；5)曲线拟合：将血管中心线处理成为平滑的血管中心曲线；6)体积计算：根据血管中心曲线将肝脏分段，并计算每一段的体积。但该技术的缺点在于：该方法是将三维肝脏分段转换成二维，是在二维的基础上进行求取像素的点数来进行体积计算，不是针对整体求体积，而是通过部分求整体的体积，这样误差相对较大。 

发明内容

本发明针对现有技术存在的上述不足，提出一种基于三维可视化和蒙特卡罗方法的体积测 量方法，实现了对不规则物体体积的测量，且测量对象不需要进行分割、分块等处理，结果精度相对较高。 

本发明是通过以下技术方案实现的，本发明首先通过对被测物体的序列图像进行三维重建，并将重建的结果置于创建得到的虚拟盒内部，然后在虚拟盒内部产生若干个随机点，最后统计三维重建后的被测物体内部的随机点的个数，并根据虚拟盒的体积和随机点的总数按比例求取被测物体的体积大小。 

本发明具体包括以下步骤： 

第一步：对被测物体的序列图像进行三维重建，具体为：先对序列图像进行抽取等值面处理、拼接三角面片处理和数据映射处理，然后采用经典的面绘制方法-移动立方体的方法实现三维重建。 

所述的抽取等值面处理采用的是移动立方体(Marching Cube)面绘制方法，通过将一系列的二维切片数据视作三维数据场，并从中将具有特定阈值的物质抽取出来以拓扑的方式连接成三角面片，具体步骤包括： 

i)每次读两张切片图像，形成一个层单元(Layer)； 

ii)两切片图像上下相对应的四个点构成一个立方体单元(Cube)，如图1所示； 

iii)先按照从左至右、从前到后的顺序将该层单元中的立方体单元进行像素切片并形成新的立方体单元；再按从下到上的顺序将对应各个立方体单元进行像素切片并形成新的立方体单元，则对应n个层单元的图像共进行了n-1次像素切片。 

第二步：将第一步得到的三维重建模型置于一构建出的虚拟盒内部，具体为：通过SetProp3D()函数将重建后的结果添加到创建出的虚拟盒内部，并进行参数的设置。SetProp3D()函数是VTK中的常用函数类，主要是对重建后的物体进行操作，在这里作用就是将重建后的物体添加到虚拟盒中。 

第三步：采用伪随机数发生器函数来产生由随机数构成坐标的随机点，并判断随机点是否落在三维重建模型内部； 

所述的随机点是指随机函数产生的三维坐标所对应的点。 

所述的随机数通过以下方式得到：采用基于ANSI(美国国家标准学会)标准的随机数发生器函数生成随机序列，具体为： 

3.1)给无符号整型函数srand()提供一个种子(seed)； 

3.2)调用函数rand()，它会根据srand()提供的种子返回随机数； 

3.3)根据需要可以多次调用rand()函数并且随时可以为srand()函数提供新的种子，从而得到不同的新的随机数； 

所述的判断随机点采用夹角之和判断法、叉积判断法或交点计数判断法中任一得以实现。第四步：根据第三步得到的判断结果进行统计，得到三维重建模型内随机点的个数，并根据数值比例关系即可求得待测不规则物体的体积，具体为： 其中：V₁是被求体积，V是虚拟盒的体积，M是生成的随机点的总数，m是判断后落入被测体内部随机点的数目。 

本发明实现了对不规则物体体积的测量，且测量对象不需要进行分割、分块等处理，结果精度相对较高。 

附图说明

图1为立方体单元示意图。 

图2为实施例效果图。 

图中：(a)为正面示意图；(b)为背面示意图，测量得到对象的体积为：6076093.00立方毫米。 

图3为实施例中点在多边体内外的判断示意图； 

图中：(a)为叉积判断法示意图，(b)夹角之和为2π判断法示意图，(c)夹角之和为0判断法示意图，(d)交点计数判断法示意图。 

图4为实施例流程示意图。 

具体实施方式

下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。 

实施例 

基于本发明使用具有普遍性，本发明中采用如图2所示的头颅作为体积测量的对象进行效果演示。 

如图4所示，本实施例包括以下步骤： 

第一步：对被测物体的序列图像进行三维重建，具体为：先对序列图像进行抽取等值面处理、拼接三角面片处理和数据映射处理，然后在本实施例中采用经典的面绘制方法-移动立方体的方法实现三维重建。 

所述的抽取等值面处理采用的是移动立方体(Marching Cube)面绘制方法，通过将一系列的二维切片数据视作三维数据场，并从中将具有特定阈值的物质抽取出来以拓扑的方式连接成三角面片，具体步骤包括： 

i)每次读两张切片图像，形成一个层单元(Layer)； 

ii)两切片图像上下相对应的四个点构成一个立方体单元(Cube)，如下图所示；(图1) 

iii)先按照从左至右、从前到后的顺序处理层中的立方体单元；再按从下到上顺序处理，如对应n个层单元的图像则处理n-1次，所述的处理是指进行切片中像素的处理，使之形成新的立方体。 

本实施例中，对三角面片的拼接是利用vtkStripper类将生成的三角面片拼接生成三角带形成等值面。 

本实施例中，采用VTK中的SetMapper()函数进行映射处理。 

第二步：通过SetProp3D()函数将重建后的结果添加到创建出的虚拟盒内部，并进行参数的设置。SetProp3D()函数是VTK中的常用函数类，主要是对重建后的物体进行操作，在这里作用就是将重建后的物体添加到虚拟盒中。此处虚拟盒是根据物体的大小来生成的，不需要刻意设置其高低，它的大小可以通过顶点坐标求出。 

第三步：采用伪随机数发生器函数来产生随机数，随机数的数目可以人为的确定，为接下来的计算提供方便，并判断点在物体的内外；此处点是指随机函数产生的三维坐标所对应的像素点；物体是指虚拟盒内部的三维构造出的对象模型。 

所述的随机数通过以下方式得到：采用基于ANSI(美国国家标准学会)标准的随机数发生器函数生成随机序列，具体为： 

3.1)给无符号整型函数srand()提供一个种子(seed)； 

3.2)调用函数rand()，它会根据srand()提供的种子返回随机数； 

3.3)根据需要可以多次调用rand()函数并且随时可以为srand()函数提供新的种子，从而得到不同的新的随机数； 

所述的判断采用以下任意一种方式实现： 

-夹角之和判断法：将任意一点和除该点以外的其他多个点组成的多边形的任意两个顶点进行连接，得到若干组两条相连线段的组合，当所有线段组合的夹角累加后为360°即2pi，则该点位于多边形内部。 

如图3所示，图3(b)中多边形Q₁Q₂...Q₅，将Q₀与各顶点相连接，其中Q_iQ₀Q_i+1＝β_i，这是∑β＝2π的情况；图3(c)中将R₀与R₁R₂...R₅分别相连，其中R_iR₀R_i+1＝γ_i，∑γ＝0的情况；图中Q_iQ₀Q_i+1＝β_i，R_iR₀R_i+1＝γ_i，当∑β＝2π，则判断点落在多边形内；当∑γ＝0，则判断点落在多边形的外部； 

-叉积判断法：该方法适用于凸多边形，点在凸多边形内的充要条件是叉积V_i×V_i+1(i＝1，2，......n)的符号相同，在三维空间内，此法仅适用于凸多边形，；如图3(a)中D₁D₂...D₆表示为多边形按逆时针排序的顶点，D₀为中心点，其中：多边形顶点按逆时针排序D₁D₂...D₆，中心点为D₀，V_i＝D_i-D₀，V_n+1＝V₁，δ_i为D_i到D_i+1的夹角，例如δ₁是D₁D₀D₂的夹角； 

交点计数判断法：当多边形为凹多边形时，可通过将任一待判断点沿同一方向作一射线至无穷远，然后求射线与多边形的交点，若交点个数为奇数，则点在多边形内，否则点在多边形外； 

如图3(d)中A、B、C、D为判断点，a、b、c、d分别对应表示射线的方向，根据交点的个数判断点在多边形内外。 

第四步：根据第三步得到的判断结果进行统计，并根据数值比例关系得到待测不规则物体的体积，具体为： 

其中：V₁是被求体积，V是虚拟盒的体积，M是生成的随机点的总数，m是判断后落入被测体内部随机点的数目。 

实验结果分析：如图1、2所示，测量对象为头颅，数据来源是可视化人项目。 

目前尚未找到成熟软件进行对不规则体积进行测量，本发明中采用对规则物体的测量，进行说明本发明数据的可靠性。用此方法进行求取规则球体的体积，对球体半径不断的改变，得出不同的结果值，如下表1所示： 

表1： 

其中：R是球体半径，经过半径的变化得出不同的测量结果值；测量值是通过该方法得出的结果；真实值是该球体真正的体积大小；误差值是真实值与测量值之差；误差率是通过下面公式得出的：误差率＝(真实值-测量值)/测量值。 

结果分析，测量值与真实值之间有一定的误差，但是都在允许的范围内。 

Claims (10)

1.一种基于三维可视化和蒙特卡罗方法的体积测量方法，其特征在于，首先通过对被测物体的序列图像进行三维重建，并将重建的结果置于创建得到的虚拟盒内部，然后在虚拟盒内部产生若干个随机点，最后统计三维重建后的被测物体内部的随机点的个数，并根据虚拟盒的体积和随机点的总数按比例求取被测物体的体积大小。

2.根据权利要求1所述的方法，其特征是，具体包括以下步骤：

第一步：对被测物体的序列图像进行三维重建；

第二步：将第一步得到的三维重建模型置于一构建出的虚拟盒内部；

第三步：采用伪随机数发生器函数来产生由随机数构成坐标的随机点，并判断随机点是否落在三维重建模型内部；

第四步：根据第三步得到的判断结果进行统计，得到三维重建模型内随机点的个数，并根据数值比例关系即可求得待测不规则物体的体积。

3.根据权利要求1或2所述的方法，其特征是，所述的重建是指：先对文件进行抽取等值面处理、拼接三角面片处理和数据映射处理，然后采用经典的面绘制方法-移动立方体的方法实现三维重建。

4.根据权利要求3所述的方法，其特征是，所述的抽取等值面处理采用的是移动立方体面绘制方法，通过将一系列的二维切片数据视作三维数据场，并从中将具有特定阈值的物质抽取出来以拓扑的方式连接成三角面片。

5.根据权利要求4所述的方法，其特征是，所述的抽取等值面处理的具体步骤包括：

i)每次读两张切片图像，形成一个层单元；

ii)两切片图像上下相对应的四个点构成一个立方体单元，如图1所示；

iii)先按照从左至右、从前到后的顺序将该层单元中的立方体单元进行像素切片并形成新的立方体单元；再按从下到上的顺序将对应各个立方体单元进行像素切片并形成新的立方体单元，则对应n个层单元的图像共进行了n-1次像素切片。

6.根据权利要求2所述的方法，其特征是，所述的随机数通过以下方式得到：采用基于ANSI标准的随机数发生器函数生成随机序列。

7.根据权利要求2所述的方法，其特征是，所述的随机数通过以下方式得到：

3.1)给无符号整型函数srand()提供一个种子(seed)；

3.2)调用函数rand()，它会根据srand()提供的种子返回随机数；

3.3)根据需要可以多次调用rand()函数并且随时可以为srand()函数提供新的种子，从而得到不同的新的随机数。

8.根据权利要求1或2所述的方法，其特征是，所述的判断随机点采用夹角之和判断法、叉积判断法或交点计数判断法中任一得以实现。

9.根据权利要求8所述的方法，其特征是，所述的夹角之和判断法是指：将任意一点和除该点以外的其他多个点组成的多边形的任意两个顶点进行连接，得到若干组两条相连线段的组合，当所有线段组合的夹角累加后为360°即2pi，则该点位于多边形内部；

所述的叉积判断法是指：点在凸多边形内的充要条件是叉积V_i×Vi₊₁(i＝1，2，……n)的符号相同；

所述的交点计数判断法是指：当多边形为凹多边形时，可通过将任一待判断点沿同一方向作一射线至无穷远，然后求射线与多边形的交点，若交点个数为奇数，则点在多边形内，否则点在多边形外。

10.根据权利要求1或2所述的方法，其特征是，所述的第四步具体为：

其中：V₁是被求体积，V是虚拟盒的体积，M是生成的随机点的总数，m是判断后落入被测体内部随机点的数目。

Images