一种基于空间不变特性的快速体积测量方法

摘要

本发明提供了一种基于空间不变特性的快速体积测量方法，所述方法包括：对所获取的图像进行特征点提取与匹配、三维重建；在三维点云中加入特征点间实际距离与在相对坐标系内与真实世界的距离比例，作为尺度信息；利用特征点三维坐标，构建Delaunay三角网格，加入边界边检测元素，拟合最优基准面；根据离散积分的思想，将三角网格投影到基准面上，求出三角网所包围的不规则物体的体积。即可利用“一个相机，一把尺”完成物体所有特征点间实际距离的测量，省去了实地测量的麻烦；同时在构建Delaunay三角网格时添加了检测边界边的数据结构，可以有效的得到所构成边界曲面的外围基准面，进一步用于体积计算。

说明书

技术领域

本发明涉及图像处理领域，特别涉及土地调查中利用近景摄影技术实现不规则体体 积自动求取的方法。

背景技术

土地是人类赖以生存和发展的宝贵自然资源，土地的持续利用是人类社会持续发展 的基石。在人口日益膨胀的今天，尤其是在我国，土地承受着巨大的压力，人地关系日 趋紧张。经过“十五”、“十一五”的努力，“天上看、地上查、网上管”等现代高科技 手段在土地监管以及土地测量领域已经全方位应用，航空、航天遥感技术和GPS全球定 位系统等空间技术在土地管理中也已经广泛应用，大大提高了土地信息获取的效率和准 确性，基本实现土地基础数据的充分共享，以及目标的精确测量，为社会经济发展决策 提供服务。

在现有的土地调查技术中，“天上看”主要是以卫星、无人机遥感为主，尽管在大 范围测量上有不可替代的作用，但是也存在先天的不足。其中，卫星遥感测量，周期长， 时效性差，精度不能达到较高的标准。无人机遥感测绘，成本高，在高程变化复杂的区 域，精度难以得到保障。“地上查”主要以GNSS/全站仪或是皮尺测量为主，这两种测 量方法都是属于与物体实际坐标有关的绝对测量。基于GNSS的测量方法流程复杂，对 工作人员要求较高，布点受地形影响较大，过度较慢，非常不适用于自动化测量及监测； 而使用皮尺测量，有大量的外业工作，只适用于小区域测量，并且对于不规则体测量、 人类不能到达或较危险区域，它的适用性较差。因此利用近景摄影技术进行土地测量， 能够为土地测量提供一条新的快速采集土地信息的新途径，能够广泛应用于三维山体测 量、滑坡监测、沙坑容积测量等许多不规则体测量计算的场合。它能构实现目标的精确 测量，具有速度较快、精度较高而均匀、成本低、不损坏原地物、不受气候及季节的限 制等优点，在土地测量中的应用具有积极的现实意义。

客观世界在空间上是三维的，而现有的图像采集装置所获取的图像是二维的。尽管 图像中含有某些形式的三维空间信息，这些信息包括物体边与边之间的几何关系、两幅 图像的视差关系、两幅图像中特征点的对应关系以及物体轮廓信息等等，但要真正在计 算机中使用这些信息进行进一步的应用处理，就必须采用三维重建技术从二维图像中合 理地提取并表达这些三维信息。结合基于图像的三维重建技术在土地测量中的应用具有 积极的现实意义，它有着传统手段无法比拟的优势，有着诸多便利之处。它可将大量的 外业测量工作移到室内完成,它具有速度较快、精度较高而均匀、成本低、不损坏原地 物、不受气候及季节的限制等优点。新型体积测量技术在土地测量中的应用广泛，能实 现精确测量的目标。由于基于图像的三维重建技术所需要的成本低，具有很大的灵活性， 又能达到简单的三维建模功能，因此，在需要真实感建模的场合，基于图像的建模无疑 具有很高的实用价值。

发明内容

发明目的：针对上述现有技术，提出一种基于空间不变特性的快速体积测量方法， 实现快速获取不规则物体体积。

技术方案：一种基于空间不变特性的快速体积测量方法，包括如下步骤：

步骤(1)，采用单台摄像机获取物体的多幅图像序列，并对所述摄像机进行标定；

步骤(2)，基于双目立体视觉技术，提取所述多幅图像上的特征点后，对获取的多 张图像序列进行立体匹配，找出任意两幅图像中同名像素点的对应坐标关系；

步骤(3)，任意选取所述图像序列中的两幅图像，以拍摄其中一幅图像的摄像机为 坐标原点建立相对坐标系，在所述相对坐标系中对所述步骤(2)提取的多幅图像序列 上的所有特征点进行三维重建，得到整个物体的三维点云结构；

步骤(4)，在所述整个物体的三维点云结构中，选择并测量两个特征点所对应物体 上的实际距离，计算得到所述实际距离与这两个特征点在所述相对坐标系中的距离的比 值，利用所述比值计算得到三维点云结构中所有特征点间的实际距离；

步骤(5)，依据重建的三维点云结构建立物体的Delaunay三角网模型，根据所述 Delaunay三角网模型获取计算物体体积的基准面，根据所述基准面以及所有特征点间的 实际距离计算物体的体积。

作为本发明的优选方案，依据三角网生长算法建立所述Delaunay三角网模型,并在 构建所述Delaunay三角网模型时加入检测三角网边界边的数据成员useCount，所述数据 成员useCount用于记录构建的每个三角形的每条边作为基础边来构建新三角形的次数； 完成所有三角形建立步骤后，选择数据成员useCount为1的所有边作为Delaunay三角网 模型的边界边，再将由所有边界边构成的封闭图形拟合成一个光滑平面，所述光滑平面 作为计算物体体积的基准面。

作为本发明的优选方案，所述步骤(2)具体为：

(21)，对获取的图像进行预处理，所述预处理包括图片平滑、锐化处理；

(22)，特征点提取：对任意两幅图像，先通过一个角点响应函数CRF去判断各个 待测点是否为FAST角点，所述CRF函数如下：

N＝Σf_CRF(I(p)，I(x))        (1)

$f_{\mathrm{CRF}}=\left(\begin{array}{cc}1& |I\left(x\right)-I\left(p\right)|>{ϵ}_d\\ 0& \mathrm{others}\end{array}\right)---\left(2\right)$

式(1)中，N为响应函数CRF值的和，即FAST角点的个数；I(x)为待测点周围 的任意一点的灰度值；I(p)表示当前待测点的灰度值；p表示当前待测点；式(2)中， ε_d为一个中心待测点与圆周任一点像素灰度差值的阈值；

设定阈值ε_d，提取出超过阈值N_f数目的FAST角点；然后利用Harris角点的评价 函数找到前N_f个较好的FAST角点；最后使用金字塔算法得到多尺度的待测图像，确定 出最终的FAST角点；

(23)，采用BRIEF描述子对提取的FAST角点进行描述；

(24)，利用ORB特征点匹配算法对选取的FAST角点进行立体匹配与解算，得到 任意两幅图片特征点匹配的基础矩阵F。

作为本发明的优选方案，所述步骤(3)具体为：

(31)，任意选取两幅图像，以拍摄其中一幅图像的摄像机为坐标原点，以它的摄 像机坐标系作为相对坐标系，计算得到拍摄这两幅图像的摄像机的位置关系，即摄像机 的外部参数；

(32)，结合摄像机标定时得到的摄像机内部参数，根据步骤(2)所得的这两幅图 像的对应基础矩阵F，计算得到其本质矩阵E；

(33)，对本质矩阵E进行奇异值分解，得到摄像机外部参数(R|t)的候选值，建 立拍摄这两幅图像的两个不同位置的相机的投影矩阵P₁和P₂：

P₁＝K[I O]         (3)

P₂＝K[R t]       (4)

其中I是3x3的单位矩阵；O是3x1的全零矩阵；K为摄像机内部参数矩阵；R为3x3 的旋转矩阵；t为三维平移向量；R、t均为相机外部参数；

(34)，根据恢复出的投影矩阵P₁和P₂，在相对坐标系中利用SFM算法进行迭代 计算，重建出这两幅图像上特征点对应的三维点云结构，即得到特征点在相对坐标系内 的坐标信息；

(35)，在所述步骤(31)建立的相对坐标系下，按照步骤(31)至步骤(34)，对 其任意两幅图像做三维重建，得到整个物体在该相对坐标系下的三维点云结构。

作为本发明的优选方案，所述步骤(5)中，根据所述基准面以及所有特征点间的 实际距离计算物体的体积具体步骤为：

(1)，设Delaunay三角网模型中任意三角形的三个顶点坐标为：A(x₁，y₁，z₁)， B(x₂，y₂，z₂)，C(x₃，y₃，z₃)，使用向量法进行投影计算任意三角形的三个顶点投射到所述 基准面上的投影点坐标A₁(x₁',y₁',z₁')、B₁(x'₂,y'₂,z'₂)、C₁(x'₃,y'₃,z'₃)；

(2)，计算Delaunay三角网模型中任意三角形与所述基准面形成的单个三棱柱整 体体积：

21)，计算底面△A₁B₁C₁的面积：

$\overline{A_1{B}_1}=\sqrt{{(x_2^{\prime }-x_1^{\prime })}^2+{(y_2^{\prime }-y_1^{\prime })}^2+{(z_2^{\prime }-z_1^{\prime })}^2}=l_1,---\left(5\right)$

$\overline{A_1{C}_1}=\sqrt{{(x_3^{\prime }-x_1^{\prime })}^2+{(y_3^{\prime }-y_1^{\prime })}^2+{(z_3^{\prime }-z_1^{\prime })}^2}=l_2,---\left(6\right)$

$\overline{B_1{C}_1}=\sqrt{{(x_3^{\prime }-x_2^{\prime })}^2+{(y_3^{\prime }-y_2^{\prime })}^2+{(z_3^{\prime }-z_2^{\prime })}^2}=l_3;---\left(7\right)$

设$p=\frac{l_1+l_2+l_3}{2},$

底面△A₁B₁C₁的面积$S_{{\mathrm{\Delta A}}_1{B}_1{C}_1}=\mathrm{sqr}\left(\begin{array}{cccc}p& (p-l_1)& (p-l_2)& (p-l_3)\end{array}\right)---\left(8\right)$

22)，采用任意三角形的三个顶点的高程平均值来计算单个三棱柱的高h：

$\left(\begin{array}{c}h=\frac{1}{3}\times [\sqrt{{(x_1-x_1^{\prime })}^2+{(y_1-y_1^{\prime })}^2+{(z_1-z_1^{\prime })}^2}\\ +\sqrt{{(x_2-x_2^{\prime })}^2+{(y_2-y_2^{\prime })}^2+{(z_2-z_2^{\prime })}^2}\\ +\sqrt{{(x_3-x_3^{\prime })}^2+{(y_3-y_3^{\prime })}^2+{(z_3-z_3^{\prime })}^2}]\end{array}\right)---\left(9\right)$

23)，计算单个三棱柱的整体体积V为：

$V=S_{{\mathrm{\Delta A}}_1{B}_1{C}_1}·h---\left(10\right)$

(3)，将Delaunay三角网模型中所有的n个三角形对应生成的三棱柱的体积相加 合，即V_i为式(10)得到的单个三棱柱的整体体积；

(4)，将计算出体积V乘以步骤(4)得到的尺度信息的立方，得到物体的总体积。

有益效果：本发明的一种基于空间不变特性的快速体积测量技术，在双目立体视觉 技术的基础上，仅利用“一个相机、一把尺”，借助空间三维重建，完成物体所有特征 点间实际距离的测量；利用一台相机获取在不同方位拍摄的多幅图像序列，算出拍摄每 两张相片的相机的相对位置，建立相对坐标系，在相对坐标系内，利用立体匹配算法和 三维重建算法对图像的特征点进行立体匹配与物体三维重建；对重建后的某两个特征点 加入相对约束(两个特征点之间实际的相对距离)，获得物体在相对坐标系内与物体实 际尺寸的比例，作为尺度信息，从而完成利用物体三维重建达到摄影测量的目的。再利 用三维重建之后特征点的三维坐标信息，构建一种自动基准面检测的Delaunay三角网 格算法，根据网格信息求取合适基准面，用于进行不规则体体积快速计算。相对于以往 摄影测量中的绝对位置测量修改为相对位置测量，本发明将通过对空间地物在相对坐标 系中进行三维重建，再利用尺度信息获取物体各点之间实际的距离，从而实现体积的快 速计算技术。技术人员只需要照相机、皮尺等测量工具，通过图像处理及三维重建技术， 便可以实现快速体积测量。

在本发明中，利用特征点的二维坐标信息，改进现有的Delaunay构网算法，在三 角网生长算法的基础上构建二维的Delaunay三角网格。在构网过程中，为了弥补 Delaunay三角网格算法在体积计算领域的局限性，本方法加入了检测三角网边界边的数 据成员useCount。若该数据成员值为0，说明对应的线段还没有在构建三角形的过程中 被使用过；若该数据成员值为1，则说明对应的线段仅仅在一个构建的三角形中被使用 过，并不涉及其他三角形，这样的三角边即为整个三角网格的边界边；若该数据成员值 为2，则说明这条线段在构建三角形的过程中被使用了两次，即不是三角网格的边界边。 构造了上述数据成员，则通过该数据成员的值即可将构成三角网边界的线段全部检测出 来，再将其位置记录下来之后，便可获得这些折线段所构成的一个面。但由于其所构成 的这个面并不是一个真正意义上的面，仅仅是一系列相连的折线段首尾接合所围成的一 个封闭图形，因此，本方法采用最小二乘法，将其拟合成一个光滑的平面，作为不规则 体体积计算的一个基准面，进一步用于Delaunay三角网所包围部分的体积计算。利用 三维重建后物体在相对坐标系内的三维坐标信息、尺度信息，以及通过三角网格拟合出 的基准面位置，将外围物体投影在基准面上，利用积分求体积的方法对三角网格所包围 的不规则体进行体积计算。

附图说明

图1是本发明的系统流程示意图。

具体实施方式

下面结合附图对本发明做更进一步的解释。

如图1所示，一种基于空间不变特性的快速体积测量方法，包括如下步骤：

步骤(1)，采用单台摄像机获取物体的多幅图像序列，并对所述摄像机进行标定：

首先，利用一台摄像机获取物体的图像序列，建立场景中物体在摄像机坐标系中与 物体在图像平面上对应像点之间的关系,具体为：利用一台摄像机在不同的位置，从不 同角度获取同一物体的多幅图像，所获取的每一幅图像均为物体的整体图像,每次小角 度的改变相机位置进行拍摄，直至所摄的图像序列包含物体的所有信息，即所获取的多 幅图像需要覆盖整个被测物体；然后进行摄像机标定与解算，获取摄像机的内部参数。

步骤(2)，基于双目立体视觉技术，提取特征点后，对获取的多张图像序列进行立 体匹配，以此找出任意两幅图像中同名像素点的对应坐标关系，即通过计算物体的空间 点在两幅图像中的视差来获得其三维坐标值，具体步骤如下：

(21)，对获取的图像进行预处理来提高图像质量，预处理包括图片平滑、锐化处 理；

(22)，特征点提取：对任意两幅图像，依据物体的物理特征来获取合适的特征点， 即将角点、边缘点等易于分辨的特征信息作为后续特征点匹配的基元。本实施例采用 FAST(Features From Accelerated Segment Test)特征点检测方法提取特征点：先通过一 个角点响应函数CRF去判断各个待测点是否为FAST角点，其待测点即易于分辨的特征 信息点，CRF函数如下：

N＝Σf_CRF(I(p)，I(x))     (1)

$f_{\mathrm{CRF}}=\left(\begin{array}{cc}1& |I\left(x\right)-I\left(p\right)|>{ϵ}_d\\ 0& \mathrm{others}\end{array}\right)---\left(2\right)$

式(1)中，N为响应函数CRF值的和，即FAST角点的个数；I(x)为待测点周围 的任意一点的灰度值；I(p)表示当前待测点的灰度值；p表示当前待测点；式(2)中， ε_d为一个中心待测点与圆周任一点像素灰度差值的阈值。

通过设定阈值ε_d，提取出超过阈值N_f数目的FAST角点；然后利用Harris角点 的评价函数找到前N_f个较好的FAST角点；最后使用金字塔算法得到多尺度的待测 图像，从而确定出最终要得到的FAST角点。

(23)，采用BRIEF(Binary Robust Independent Element Feature)描述子对提取的 FAST角点进行描述。

(24)，利用ORB(Oriented FAST and rotated BRIEF)特征点匹配算法对选取的FAST 角点进行立体匹配与解算，得到任意两幅图片特征点匹配的基础矩阵F；即多幅图像中， 任意每两幅图像组合均得到其对应的一个基础矩阵F。

步骤(3)，任意选取两幅图像，以拍摄其中一幅图像的摄像机为坐标原点建立相对 坐标系，在该相对坐标系中对步骤(22)最终确定的所有FAST角点进行三维重建，得 到整个物体的三维点云结构，具体为：

(31)，任意选取两幅图像，以拍摄其中一幅图像的摄像机为坐标原点，以它的摄 像机坐标系作为相对坐标系，计算得到拍摄这两幅图像的摄像机的位置关系，即摄像机 的外部参数；

(32)，结合摄像机标定时计算得到的摄像机内部参数，根据步骤(2)所得的这两 幅图像的对应基础矩阵F，计算得到其本质矩阵E；

(33)，对本质矩阵E进行奇异值分解，得到摄像机外部参数(R|t)的候选值，建 立拍摄这两幅图像的两个不同位置的相机的投影矩阵P₁和P₂：

P₁＝K[I O]      (3)

P₂＝K[R t]      (4)

其中I是3x3的单位矩阵；O是3x1的全零矩阵；K为摄像机内部参数矩阵；R为3x3 的旋转矩阵；t为三维平移向量；R、t均为相机外部参数；

(34)，根据恢复出的投影矩阵P₁和P₂，在相对坐标系中利用SFM(Structure From  Motion)算法进行迭代计算，重建出这两幅图像上特征点对应的三维点云结构，即得到 特征点在相对坐标系内的坐标信息。

由于本方案只使用了一台相机，故相机的内部参数是不变的；因此，在上述步骤(31) 建立的相对坐标系下，按照步骤(31)至步骤(34)的方法，对其任意两幅图像做三维 重建，得到整个物体在该相对坐标系下的三维点云结构。

步骤(4)，在已知了所有选取的特征点在相对坐标系内的坐标信息后，计算出各个 特征点之间的相对距离信息；然后在特征点中选取两个特征点，通过用皮尺丈量真实物 体上着两个特征点的实际距离，实际距离与这两个特征点在相对坐标系内的相对距离的 比值，即作为尺度信息。由此尺度信息可从三维点云结构中的任意两个特征点之间的相 对距离计算出其在真实世界中的实际距离。

步骤(5)，依据重建的三维点云结构建立物体的Delaunay三角网模型，基于该 Delaunay三角网模型和所有特征点间的实际距离计算物体的体积，具体步骤如下：

(51)，将重建后的三维点云结构投影到二维平面上，得到重建出来的特征点的二 维坐标信息，然后依据三角网生长算法构建Delaunay三角网模型。构建Delaunay三角 网格，即构建三种相关的数据结构：点、线、面，拓扑数据结构如下：

(1)点(T_point)类封装了离散采样点的三维坐标，构造出点的数据结构；

(2)线(T_Line)类封装了三角形的边的信息，其中每条边的两个端点的数据可 以存储采样点在容器中的序号。

(3)面(T_tri)类封装了生成的每个三角形的信息，三角形的边信息以序号类型 放入类中，可以随时调用，数据结构由边扩展到了面，用于最终三角网格的生成。

(4)三角网(T_TIN)类包含了构建Delaunay三角网格的核心函数，有了前面点、 线、面数据结构的构建基础，可以方便构建出Delaunay三角网格。

构建上述三种具有拓扑关系的数据结构，能够便于对重建点云的整合和管理，利用 三角网生长算法的思想进行三角网格的构建。

在构建Delaunay三角网模型时加入检测三角网边界边的数据成员useCount，数据 成员useCount用于记录构建的每个三角形的每条边作为基础边来构建新三角形的次数； 完成所有三角形建立步骤后，选择数据成员useCount为1的所有边作为Delaunay三角网 模型的边界边。具体步骤如下：

(a)，首先将重建出来的特征点的二维坐标信息放入初始数组容器中，并进行编号 保存；

(b)，以编号为0和1的点作为起始点，按照三角网判别法则构建一个寻找最大最 小角的函数来找到满足Delaunay构网条件的第三个点，作为2号点；

(c),以步骤(b)中确立的3点作为初始三角形，并判断该三角形的三条边是否 为可拓展边：若为可拓展边，则以可拓展边的两顶点作为基准再次运用寻找最大最小角 的函数找寻第三个点，并与第三个点构成新的三角形，从而完成新三角网的拓展；若不 是可拓展边，则停止拓展并给予该边标记，该不可拓展的边即是边界边；

(d)，每遍历一条边，该边的数据成员useCount的值就相应的增加1，直到Delaunay 三角网构网结束后，统计useCount的值，若该数据成员值为0，说明这条线段还没有在 构建三角形中被使用过；若该数据成员为1，则说明该线段仅仅在一个三角形中被使用， 即为整个三角网的边界边；若该数据成员为2，则说明这条线段在构建三角形的过程中 被使用了两次，即非网格的边界边。

(52)，完成所有三角形构建步骤后，根据数据成员useCount得到Delaunay三角网 模型的所有边界边，再将由所有边界边构成的封闭图形拟合成一个光滑平面，该光滑平 面作为计算物体体积的基准面，具体步骤如下：

在获取了所有的边界边信息后，找到所有边界边的两个端点的编号，并记录其端点 的三维坐标信息；设待拟合的最优基准面的空间平面方程的表达式为：

ax+by+cz-d＝0      (5)

其中，a、b、c、d均为方程系数，x、y、z为基准面坐标值；

要获得最佳拟合平面，应在条件a²+b²+c²＝1的约束下，获取式(6)中e的最小 值，即误差最小值：

$e={\Sigma }_{i=0}^n{({\mathrm{ax}}_i+{\mathrm{by}}_i+{\mathrm{cz}}_i-d)}^2---\left(6\right)$

其中，x_i、y_i、z_i分别为空间平面上任意一点的坐标值；n为空间平面上点的个数； 式(6)即等式约束的优化问题，利用拉格朗日乘数法求函数极值，令 构造优化模型f(a，b，c，d，λ)：

其中，λ为拉格朗日乘数；

令f分别对a，b，c，d求偏导，并令偏导数为零，整理如下：

由$\frac{{\partial }_f}{{\partial }_d}=0$得：

$d={\Sigma }_{i=0}^n({\mathrm{ax}}_i+{\mathrm{by}}_i+{\mathrm{cz}}_i)/n---\left(8\right)$

令$\bar{x}={\Sigma }_{i=0}^n{x}_i/n,\bar{y}={\Sigma }_{i=0}^n{y}_i/n,\bar{z}={\Sigma }_{i=0}^n{z}_i/n,$带入得：

令${\mathrm{\Delta x}}_i=x_i-\bar{x},{\mathrm{\Delta y}}_i=y_i-\bar{y},{\mathrm{\Delta z}}_i=z_i-\bar{z},$由$\frac{{\partial }_f}{{\partial }_a}=0$得：

$2{\Sigma }_{i=0}^n({\mathrm{a\Delta x}}_i+{\mathrm{b\Delta y}}_i+{\mathrm{c\Delta z}}_i){\mathrm{\Delta x}}_i-2\mathrm{\lambda a}=0---\left(10\right)$

同理，由$\frac{{\partial }_f}{{\partial }_b}=0$和$\frac{{\partial }_f}{{\partial }_c}=0$得：

$2{\Sigma }_{i=0}^n({\mathrm{a\Delta x}}_i+{\mathrm{b\Delta y}}_i+{\mathrm{c\Delta z}}_i){\mathrm{\Delta x}}_i-2\mathrm{\lambda b}=0---\left(11\right)$

$2{\Sigma }_{i=0}^n({\mathrm{a\Delta x}}_i+{\mathrm{b\Delta y}}_i+{\mathrm{c\Delta z}}_i){\mathrm{\Delta x}}_i-2\mathrm{\lambda c}=0---\left(12\right)$

将式(10)，式(11)和式(12)整理成矩阵形式，得：

$\left(\begin{array}{ccc}{\Sigma }_{i=0}^n{{\mathrm{\Delta x}}_i}^2& {\Sigma }_{i=0}^n{\mathrm{\Delta x}}_i{\mathrm{\Delta y}}_i& {\Sigma }_{i=0}^n{\mathrm{\Delta x}}_i{\mathrm{\Delta z}}_i\\ {\Sigma }_{i=0}^n{\mathrm{\Delta x}}_i{\mathrm{\Delta y}}_i& {\Sigma }_{i=0}^n{{\mathrm{\Delta y}}_i}^2& {\Sigma }_{i=0}^n{\mathrm{\Delta y}}_i{\mathrm{\Delta z}}_i\\ {\Sigma }_{i=0}^n{\mathrm{\Delta x}}_i{\mathrm{\Delta z}}_i& {\Sigma }_{i=0}^n{\mathrm{\Delta y}}_i{\mathrm{\Delta z}}_i& {\Sigma }_{i=0}^n{{\mathrm{\Delta z}}_i}^2\end{array}\right)\left(\begin{array}{c}a\\ b\\ c\end{array}\right)=\lambda \left(\begin{array}{c}a\\ b\\ c\end{array}\right)---\left(13\right)$

由上述矩阵形式看出，可以将等式约束的优化问题转化为求解矩阵特征值和特征向 量的问题，令$A=\left(\begin{array}{ccc}{\Sigma }_{i=0}^n{{\mathrm{\Delta x}}_i}^2& {\Sigma }_{i=0}^n{\mathrm{\Delta x}}_i{\mathrm{\Delta y}}_i& {\Sigma }_{i=0}^n{\mathrm{\Delta x}}_i{\mathrm{\Delta z}}_i\\ {\Sigma }_{i=0}^n{\mathrm{\Delta x}}_i{\mathrm{\Delta y}}_i& {\Sigma }_{i=0}^n{{\mathrm{\Delta y}}_i}^2& {\Sigma }_{i=0}^n{\mathrm{\Delta y}}_i{\mathrm{\Delta z}}_i\\ {\Sigma }_{i=0}^n{\mathrm{\Delta x}}_i{\mathrm{\Delta z}}_i& {\Sigma }_{i=0}^n{\mathrm{\Delta y}}_i{\mathrm{\Delta z}}_i& {\Sigma }_{i=0}^n{{\mathrm{\Delta z}}_i}^2\end{array}\right),$则上式化简为：

$A=\left(\begin{array}{c}a\\ b\\ c\end{array}\right)=\lambda \left(\begin{array}{c}a\\ b\\ c\end{array}\right),$λ为矩阵A的特征值，而α＝[a b c]^T即为矩阵A的特征向量；由矩阵 特征值的求解方法可知：λ＝(Aα·α)/(α·α），由于a²+b²+c²＝1，故α·α＝1， 所以所以λ的最小值就是的最小 值。

所以求出矩阵A的最小特征值λ_min对应的特征向量α，就解算出了系数a，b，c 的值，带入式(8)即可求出d的值。

由此，利用带约束的最小二乘算法将边界边拟合成了一个最优基准面，将用于后面 的不规则体体积计算。

(53)，利用拟合得到的最优基准面以及步骤(4)得到的尺度信息，利用积分求体 积的方法求取物体的体积，即采用基于三棱柱体积计算的离散积分算法计算物体的体 积：计算由Delaunay三角网格剖分后的每一个三角形与最优基准面投影构成的小三棱 柱的体积，最后将所有三棱柱的体积进行加合，就得到了整个三角网所包围的物体的体 积，具体的步骤如下：

设Delaunay三角网模型中任意三角形的三个顶点坐标为：A(x₁，y₁，z₁)，B(x₂，y₂，z₂)， C(x₃，y₃，z₃)；三个顶点投射到最优基准面的空间平面上的投影点为A₁、B₁、C₁；已知最 优基准面的空间平面的方程为：ax+by+cz＝d；

(1)计算任意三角形的三个顶点投射到最优基准面上的投影点坐标A₁(x₁',y₁',z₁')、 B₁(x'₂,y'₂,z'₂)、C₁(x'₃,y'₃,z'₃)：

本实施例中使用向量法进行投影计算：

$G^{\prime }=G-\overrightarrow{n}*(\overrightarrow{\mathrm{QG}}*\overrightarrow{n})---\left(14\right)$

其中，G点为重建出三维坐标的每个特征点；向量为单位化的空间平面法向量；Q 点为空间平面上的任意一点；G′点为投影结果。

(2)计算Delaunay三角网模型中任意三角形与最优基准面形成的单个三棱柱整体 体积：

1)计算底面△A₁B₁C₁的面积：

$\overline{A_1{B}_1}=\sqrt{{(x_2^{\prime }-x_1^{\prime })}^2+{(y_2^{\prime }-y_1^{\prime })}^2+{(z_2^{\prime }-z_1^{\prime })}^2}=l_1,---\left(15\right)$

$\overline{A_1{C}_1}=\sqrt{{(x_3^{\prime }-x_1^{\prime })}^2+{(y_3^{\prime }-y_1^{\prime })}^2+{(z_3^{\prime }-z_1^{\prime })}^2}=l_2,---\left(16\right)$

$\overline{B_1{C}_1}=\sqrt{{(x_3^{\prime }-x_2^{\prime })}^2+{(y_3^{\prime }-y_2^{\prime })}^2+{(z_3^{\prime }-z_2^{\prime })}^2}=l_3,---\left(17\right)$

设$p=\frac{l_1+l_2+l_3}{2},$

底面△A₁B₁C₁的面积$S_{{\mathrm{\Delta A}}_1{B}_1{C}_1}=\mathrm{sqr}\left(\begin{array}{cccc}p& (p-l_1)& (p-l_2)& (p-l_3)\end{array}\right);---\left(18\right)$

2)计算单个三棱柱的高：

本方法中三棱柱的高采用三个点的高程平均值来估计，能够更加接近于真实体积， 计算三棱柱的高步骤如下：

$\left(\begin{array}{c}h=\frac{1}{3}\times [\sqrt{{(x_1-x_1^{\prime })}^2+{(y_1-y_1^{\prime })}^2+{(z_1-z_1^{\prime })}^2}\\ +\sqrt{{(x_2-x_2^{\prime })}^2+{(y_2-y_2^{\prime })}^2+{(z_2-z_2^{\prime })}^2}\\ +\sqrt{{(x_3-x_3^{\prime })}^2+{(y_3-y_3^{\prime })}^2+{(z_3-z_3^{\prime })}^2}]\end{array}\right)---\left(19\right)$

3)计算单个三棱柱的整体体积V为：

$V=S_{{\mathrm{\Delta A}}_1{B}_1{C}_1}·h---\left(20\right)$

以上过程即为每个小三角网格与基准面构成的小三棱柱的体积，将所有三角形生 成的三棱柱的体积相加合，即V_i即式(10)得到的单个三棱柱的整体体积； 最后，将计算出体积V乘以步骤(4)得到的尺度信息的立方，即可以得到物体的总体 积。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员 来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也 应视为本发明的保护范围。