CN1885349A

CN1885349A - Point cloud hole repairing method for three-dimensional scanning

Info

Publication number: CN1885349A
Application number: CN 200610088244
Authority: CN
Inventors: 达飞鹏; 朱春红
Original assignee: Southeast University
Current assignee: Southeast University
Priority date: 2006-07-05
Filing date: 2006-07-05
Publication date: 2006-12-27

Abstract

The invention discloses a three-dimensional scanning point cloud hole filling method capable of providing curved surface fitting accuracy. Select the coordinate point P _s around the point cloud hole as the preliminary fitting surface, and use the least square method to fit the selected coordinate point P _s to obtain the preliminary fitting surface S(u, v, w). According to the approximation method based on the Newton iteration method, the fitting surface is further optimized to improve the fitting accuracy. Take the line based on the curvature and other parameters on the fitted surface, and then take points on the line and other parameters to find the point to fill the hole, so as to achieve smooth fitting and filling of the hole. The invention adopts the surface fitting method to obtain a curved surface that accurately fits the scattered point sets around the hole, and can guarantee certain smoothness and shape retention.

Description

Hole filling method in point cloud of 3D scanning

技术领域technical field

本发明涉及一种对三维图形的修补方法，尤其涉及一种三维扫描的点云孔洞填补方法。The invention relates to a method for repairing three-dimensional graphics, in particular to a method for filling holes in point clouds of three-dimensional scanning.

背景技术Background technique

逆向工程(Reverse Engineering，RE)技术是20世纪80年代后期出现在先进制造领域里的新技术，其一般包括四个基本环节：三维形体检测与转换(物理数据的获得)、数据预处理(点云处理、识别、多视拼接)，CAD模型的建立(曲面重构)、CAM制件成型，其基本流程图如附图1所示。在三维形体检测与转换的过程中，通过三维数字化扫描仪对实物模型表面进行三维快速扫描测量，在满足离散采样速度和数据质量的前提下，获取产品的三维离散数据，点云孔洞的出现造成了这些数据的不完整，因此对孔洞进行补偿即对数据进行预处理是逆向工程中非常关键的承上启下的一环，直接影响重构成功与否和CAD模型的质量，对其后续环节起着非常关键的制约作用。本发明主要是涉及到在逆向工程的数字化过程中，用三维扫描系统(附图2)获取产品点云模型出现孔洞时的一种自动补偿方法。Reverse Engineering (Reverse Engineering, RE) technology is a new technology that appeared in the field of advanced manufacturing in the late 1980s. It generally includes four basic links: three-dimensional shape detection and conversion (obtaining physical data), data preprocessing (point Cloud processing, recognition, multi-view splicing), CAD model establishment (surface reconstruction), CAM part molding, the basic flow chart is shown in Figure 1. In the process of 3D shape detection and conversion, the 3D digital scanner is used to perform 3D rapid scanning measurement on the surface of the physical model, and the 3D discrete data of the product can be obtained under the premise of satisfying the discrete sampling speed and data quality. The incompleteness of these data is eliminated, so the compensation for holes, that is, the preprocessing of data, is a very critical link in reverse engineering, which directly affects the success of the reconstruction and the quality of the CAD model, and plays a very important role in the subsequent links. key constraints. The present invention mainly relates to an automatic compensation method when a hole appears in a product point cloud model acquired by a three-dimensional scanning system (accompanying drawing 2) during the digitalization process of reverse engineering.

近年来，点云孔洞的填补算法在国内外均取得了很大的进展，已发表了相当数量的文献，其中一些算法获得了较为广泛的应用，如基于能量优化和细分、三角网格模型、网格曲面模型中孔洞的填补算法。这些算法需先对三维扫描直接获得的点云模型做一定的前期处理或者对点云孔洞进行边界识别，实时性不强，复杂度也比较高。在实际的工程应用中，对出现的问题应能及时解决。In recent years, point cloud hole filling algorithms have made great progress at home and abroad, and a considerable number of documents have been published, some of which have been widely used, such as energy optimization and subdivision, triangular mesh model , Hole filling algorithm in mesh surface model. These algorithms need to do some pre-processing on the point cloud model directly obtained by 3D scanning or identify the boundary of the point cloud hole, the real-time performance is not strong, and the complexity is relatively high. In practical engineering applications, problems should be solved in time.

文献“Minimal energy surfaces using parametric splines.”(Gregory E.Fasshauer，LarryL.Schumarker.Computer Aided Geometric Design，1996，13：45～79)通过求解以“变形能量函数”为基础的优化目标函数，实现对孔洞的填补，保证了一定的光顺性，然而这类方法大多需要细分曲面，并且对多个子曲面片进行拼接，因此对曲面边界的连续性要求较高，通常很难达到曲面间的二阶连续。由于现有的三角网格模型的相关算法比较成熟，现有的算法大多都是基于此模型的，例如文献“A study of stereolithography file errors and repair”(Leong K F，ChuaC K，Ng Y M.International Journal of Advanced Manufacturing，1996，12：415～422)将孔洞填补归结为一个空间多边形的三角剖分问题。但大多在构造新三角片时仅仅采用原有的孔洞多边形顶点，没有增加新的三角片顶点，因而难以获得较好的用于填补孔洞的三角片形状，填补效果不够理想。并且这类算法对于复杂的点云数据编程建立三角网格模型运算量较大，在填补孔洞的过程中还需要改变三角网格的拓扑结构，因此具有对三角网格模型修改和再设计能力不足的缺点，复杂度高，这些缺点都限制了它在实践中的应用。对一类多边形孔洞的填充算法不适用于一般的点云模型。对于一般曲面网格的填充算法均需要准确的获得孔洞边界信息，难度大，不易实现。通过对点云孔洞周围散乱点集进行曲面拟合补偿孔洞的算法过程中，散乱点的参数化必不可少。通常有均匀参数化、向心参数化和累积弦长参数化方法，而这些方法主要针对呈拓扑矩形阵列的数据点，对无规则分布的点云数据，需要对其排序，难度也比较大。The literature "Minimal energy surfaces using parametric splines." (Gregory E. Fasshauer, Larry L. Schumarker. Computer Aided Geometric Design, 1996, 13: 45-79) realizes the optimization objective function based on the "deformation energy function". The filling of the holes ensures a certain level of smoothness. However, most of these methods require surface subdivision and splicing of multiple sub-surface pieces. Therefore, the continuity of the surface boundaries is high, and it is usually difficult to achieve two-dimensionality between surfaces. order continuous. Since the relevant algorithms of the existing triangular mesh model are relatively mature, most of the existing algorithms are based on this model, for example, the literature "A study of stereolithography file errors and repair" (Leong K F, ChuaC K, Ng Y M. International Journal of Advanced Manufacturing, 1996, 12: 415-422) attributed hole filling to a triangulation problem of a spatial polygon. However, most of them only use the original hole polygon vertices when constructing new triangles, without adding new triangle vertices, so it is difficult to obtain a better triangle shape for filling holes, and the filling effect is not ideal. Moreover, this type of algorithm has a large amount of computation for programming complex point cloud data to establish a triangular mesh model. In the process of filling holes, the topology of the triangular mesh needs to be changed, so it has insufficient ability to modify and redesign the triangular mesh model. The shortcomings of high complexity limit its application in practice. The filling algorithm for a class of polygonal holes is not suitable for general point cloud models. For the filling algorithm of general surface mesh, it is necessary to obtain the boundary information of holes accurately, which is difficult and difficult to realize. The parameterization of scattered points is essential in the algorithm process of compensating holes by surface fitting to the set of scattered points around holes in the point cloud. There are usually uniform parameterization, centripetal parameterization, and cumulative chord length parameterization methods, and these methods are mainly for data points in a topological rectangular array. For irregularly distributed point cloud data, they need to be sorted, which is relatively difficult.

传统的三角曲面插值算法也得到了较为广泛的应用，同时也获得了较好的效果。例如文献“An Adaptive Method for Smooth Surface Approximation to Scattered 3D Points.”(Park H，KimK.ComputerAided Design，1995，27(12)：929～939)中提出的构造光滑曲面拟合散乱点的自适应插值算法。这类插值算法虽然速度比较快，但有其致命的缺点，即：1只能处理单值的非封闭曲面；2要假定偏导矢在域边界曲线上呈线性分布；3重构的曲面依赖于估算出的偏导矢，所以曲面插值方法对离散点的逼近效果并不好，并且计算量大效率低。The traditional triangular surface interpolation algorithm has also been widely used and achieved good results. For example, the paper "An Adaptive Method for Smooth Surface Approximation to Scattered 3D Points." (Park H, KimK. Computer Aided Design, 1995, 27(12): 929-939) proposes an adaptive interpolation method for fitting scattered points on a smooth surface algorithm. Although this type of interpolation algorithm is relatively fast, it has its fatal shortcomings, namely: 1 can only deal with single-valued non-closed surfaces; 2 it must assume that the partial derivative vector is linearly distributed on the domain boundary curve; 3 the reconstructed surface depends on Because of the estimated partial derivative vector, the surface interpolation method does not have a good approximation effect on discrete points, and the calculation load is large and the efficiency is low.

在实际的应用过程中，用三维扫描系统得到的通常是具有海量散乱点的点云原始模型。此时就不可避免的存在孔洞现象，对这些孔洞进行及时的填补是关键。为此，本发明中使用一种新算法对孔洞进行补偿。如何对散乱点参数化以及提高曲面的拟合精度是此类算法中的难点。本发明中提出的对散乱点进行参数化的方法以及对曲面的迭代逼近优化方法能达到对孔洞的光滑拟合填充。In the actual application process, the original point cloud model with a large number of scattered points is usually obtained by the 3D scanning system. At this time, there will inevitably be holes, and timely filling of these holes is the key. For this reason, a new algorithm is used in the present invention to compensate for holes. How to parameterize the scattered points and improve the fitting accuracy of the surface is the difficulty in this kind of algorithm. The method for parameterizing scattered points and the iterative approximation optimization method for curved surfaces proposed in the invention can achieve smooth fitting and filling of holes.

发明内容Contents of the invention

本发明提供一种能够提供曲面拟合精确度的三维扫描的点云孔洞填补方法。The invention provides a point cloud hole filling method of three-dimensional scanning that can provide surface fitting accuracy.

本发明采用如下技术方案：The present invention adopts following technical scheme:

一种基于三角域贝塞尔曲面的三维扫描的点云孔洞填补方法：A point cloud hole filling method based on 3D scanning of triangular Bezier surfaces:

第一步：在点云孔洞周围且在屏幕坐标平面内设定一个三角形ABC，该三角形ABC的区域范围能使点云孔洞及其周围的点向屏幕坐标平面的投影落入三角形ABC内，并将投影落入三角形ABC内的点作为补孔时拟合曲面的点P_s(s＝0，1，…，m-1)，根据拟合曲面的点P_s在三角形ABC平面的投影P_s’位置计算其曲面参数化坐标(u_s，v_s，w_s)，u_s＝(ΔAP_s’B面积)/(ΔABC面积)、v_s＝(ΔAP_s’C面积)/(ΔABC面积)、w_s＝(ΔBP_s’C面积)/(ΔABC面积)，将拟合曲面的点P_s(s＝0，1，…，m-1)的坐标及其曲面参数化坐标(u_s，v_s，w_s)代入n次Bezier曲面方程并用最小二乘法得到n次Bezier曲面的控制点，从而得到初步拟合的曲面S(u，v，w)；Step 1: Set a triangle ABC around the point cloud hole and in the screen coordinate plane. The area of the triangle ABC can make the projection of the point cloud hole and its surrounding points to the screen coordinate plane fall into the triangle ABC, and The point that the projection falls into the triangle ABC is used as the point P _s (s=0, 1, ..., m-1) of the fitted surface when filling the hole, according to the projection P _s of the point P _s of the fitted surface on the triangle ABC plane 'Position is calculated by its surface parameterized coordinates (u _s , v _s , w _s ), u _s = (ΔAP _s 'B area)/(ΔABC area), v _s = (ΔAP _s 'C area)/(ΔABC area) , w _s =(ΔBP _s 'C area)/(ΔABC area), the coordinates of the point P _s (s=0,1,...,m-1) of the fitted surface and its surface parameterized coordinates (u _s , v _s , w _s ) are substituted into the nth degree Bezier surface equation and the control points of the nth degree Bezier surface are obtained by the least square method, so as to obtain the preliminary fitting surface S(u, v, w);

第二步：求出各个拟合曲面的点P_s(s＝0，1，…，m-1)到曲面的距离向量d(u，v，w)及曲面分别对对应点参数方向的偏微分S_u(u，v，w)，S_v(u，v，w)，S_w(u，v，w)，令：Step 2: Calculate the distance vector d(u, v, w) from the point P _s (s=0, 1, ..., m-1) of each fitted surface to the surface and the deviation of the surface to the parameter direction of the corresponding point Differentiate S _u (u, v, w), S _v (u, v, w), S _w (u, v, w), let:

$\{\begin{matrix} f f ((u u,, v v,, w w)) = = d d ((u u,, v v,, w w)) \cdot &Center Dot; {S S}_{u u} ((u u,, v v,, w w)) = = 00 \\ g g ((u u,, v v,, w w)) = = d d ((u u,, v v,, w w)) \cdot &Center Dot; {S S}_{v v} ((u u,, v v,, w w)) = = 00 \\ h h ((u u,, v v,, w w)) = = d d ((u u,, v v,, w w)) \cdot &Center Dot; {S S}_{w w} ((u u,, v v,, w w)) = = 00 \end{matrix}$

拟合曲面的点P_s的参数化坐标为(u_s，v_s，w_s)，以(u_s，v_s，w_s)为初始估计值，根据牛顿迭代法求解以上方程组，有：The parameterized coordinates of the point P _s of the fitted surface are (u _s , v _s , w _s ), taking (u _s , v _s , w _s ) as the initial estimated value, and solving the above equations according to the Newton iterative method, there are:

Hσ^T＝κ^T Hσ ^T =κ ^T

式中：In the formula:

σ＝(δu，δv，δw)，其中，δu，δv，δw为曲面u，v，w三个方向上的迭代步长。σ=(δu, δv, δw), where δu, δv, δw are the iteration step sizes in the three directions of surface u, v, w.

κ＝-(f(u_s，v_s，w_s)，g(u_s，v_s，w_s)，h(u_s，v_s，w_s))κ=-(f(u _s , v _s , w _s ), g(u _s , v _s , w _s ), h(u _s , v _s , w _s ))

$H h = = [\begin{matrix} {f f}_{u u} & {f f}_{v v} & {f f}_{w w} \\ {g g}_{u u} & {g g}_{v v} & {g g}_{w w} \\ {h h}_{u u} & {h h}_{v v} & {h h}_{w w} \end{matrix}] = = [\begin{matrix} {| | | | {S S}_{u u} | | | |}^{22} + + d d \cdot &Center Dot; {S S}_{uu u u} & {S S}_{v v} \cdot \cdot {S S}_{u u} + + d d \cdot &Center Dot; {S S}_{uv uv} & {S S}_{w w} \cdot &Center Dot; {S S}_{u u} + + d d \cdot &Center Dot; {S S}_{uw uw} \\ {S S}_{u u} \cdot &Center Dot; {S S}_{v v} + + d d \cdot &Center Dot; {S S}_{vu v u} & {| | | | {S S}_{v v} | | | |}^{22} + + d d \cdot \cdot {S S}_{vv vv} & {S S}_{w w} \cdot \cdot {S S}_{v v} + + d d \cdot &Center Dot; {S S}_{vw vw} \\ {S S}_{u u} \cdot &Center Dot; {S S}_{w w} + + d d \cdot &Center Dot; {S S}_{wu wu} & {S S}_{v v} \cdot &Center Dot; {S S}_{w w} + + d d \cdot &Center Dot; {S S}_{wv wv} & {| | | | {S S}_{w w} | | | |}^{22} + + d d \cdot \cdot {S S}_{ww ww} \end{matrix}]$

式中，d＝d(u_s，v_s，w_s)；f_u，f_v，f_w，g_u，g_v，g_w，h_u，h_v，h_w分别表示在点(u_s，v_s，w_s)处相应的向量对u，v，w的一阶偏导数；S_uu，S_uv，S_uw，S_vv，S_vu，S_vw，S_ww，S_wu，S_wv是曲面S(u，v，w)在点(u_s，v_s，w_s)处分别对u，v，w的二阶偏导数。In the formula, d=d(u _s , v _s , w _s ); f _u , f _v , f _w , g _u , g _v , g _w , h _u , h _v , _h _w respectively represent , v _s , w _s ) corresponding to the first-order partial derivative of the vector pair u, v, w; S _uu , S _uv , S _uw , S _vv , S _vu , S _vw , S _ww , S _wu , S _wv are The second order partial derivatives of surface S(u, v, w) at point (u _s , v _s , w _s ) with respect to u, v, w respectively.

则可得：Then you can get:

$[\begin{matrix} {| | | | {S S}_{u u} | | | |}^{22} + + d d \cdot \cdot {S S}_{uu u u} & {S S}_{v v} \cdot \cdot {S S}_{u u} + + d d \cdot \cdot {S S}_{uv uv} & {S S}_{w w} \cdot \cdot {S S}_{u u} + + d d \cdot \cdot {S S}_{uw uw} \\ {S S}_{u u} \cdot \cdot {S S}_{v v} + + d d \cdot &Center Dot; {S S}_{vu v u} & {| | | | {S S}_{v v} | | | |}^{22} + + d d \cdot &Center Dot; {S S}_{vv vv} & {S S}_{w w} \cdot &Center Dot; {S S}_{v v} + + d d \cdot &Center Dot; {S S}_{vw vw} \\ {S S}_{u u} \cdot \cdot {S S}_{w w} + + d d \cdot &Center Dot; {S S}_{wu wu} & {S S}_{v v} \cdot &Center Dot; {S S}_{w w} + + d d \cdot \cdot {S S}_{wv wv} & {| | | | {S S}_{w w} | | | |}^{22} + + d d \cdot &Center Dot; {S S}_{ww ww} \end{matrix}] [\begin{matrix} δu δ u \\ δv δ v \\ δw δw \end{matrix}] = = - - [\begin{matrix} f f (({u u}_{s the s},, {v v}_{s the s},, {w w}_{s the s})) \\ g g (({u u}_{s the s},, {v v}_{s the s},, {w w}_{s the s})) \\ h h (({u u}_{s the s},, {v v}_{s the s},, {w w}_{s the s})) \end{matrix}]$

根据三角Bezier曲面的定义可知δu+δv+δw＝0，即δw＝-δu-δv，因此得到如下迭代方程组：According to the definition of triangular Bezier surface, it can be known that δu+δv+δw=0, that is, δw=-δu-δv, so the following iterative equations are obtained:

$[\begin{matrix} {| | | | {S S}_{u u} | | | |}^{22} - - {S S}_{w w} \cdot \cdot {S S}_{u u} + + d d \cdot &Center Dot; (({S S}_{uu u u} - - {S S}_{uw uw})) & {S S}_{v v} \cdot \cdot {S S}_{u u} - - {S S}_{w w} \cdot &Center Dot; {S S}_{u u} + + d d (({S S}_{uv uv} - - {S S}_{uw uw})) \\ {S S}_{u u} \cdot &Center Dot; {S S}_{v v} - - {S S}_{w w} \cdot &Center Dot; {S S}_{v v} + + d d (({S S}_{vu v u} - - {S S}_{vw vw})) & {| | | | {S S}_{v v} | | | |}^{22} - - {S S}_{w w} \cdot &Center Dot; {S S}_{v v} + + d d \cdot &Center Dot; (({S S}_{vv vv} - - {S S}_{vw vw})) \\ {S S}_{u u} \cdot \cdot {S S}_{w w} - - {| | | | {S S}_{u u} | | | |}^{22} + + d d \cdot &Center Dot; (({S S}_{wu wu} - - {S S}_{ww ww})) & {S S}_{v v} \cdot \cdot {S S}_{w w} - - {| | | | {S S}_{w w} | | | |}^{22} + + d d \cdot \cdot (({S S}_{wv wv} - - {S S}_{ww ww})) \end{matrix}] [\begin{matrix} δu δu \\ δv δ v \end{matrix}] = = - - [\begin{matrix} f f (({u u}_{s the s},, {v v}_{s the s},, {w w}_{s the s})) \\ g g (({u u}_{s the s},, {v v}_{s the s},, {w w}_{s the s})) \\ h h (({u u}_{s the s},, {v v}_{s the s},, {w w}_{s the s})) \end{matrix}]$

其中，δu，δv为u，v两个方向上的迭代步长。S_uu，S_uv，S_uw，S_vv，S_vu，S_vw，S_ww，S_wu，S_wv是曲面S(u，v，w)在点(u_s，v_s，w_s)处分别对u，v，w的二阶偏导数，迭代求解直到 $1 / m Σ_{s = 0}^{m - 1} | | d (u_{s} + δu, v_{s} + δv, w_{s} + δw) - d (u_{s}, v_{s}, w_{s}) | | \leq ϵ,$ ε为预置曲面拟合精度，从而得到最终确定的三角Bezier曲面S’(u，v，w)；Among them, δu and δv are the iteration step sizes in the two directions of u and v. S _uu , S _uv , S _uw , S _vv , S _vu , S _vw , S _ww , S _wu , S _wv are surfaces S(u, v, w) at points (u _s , v _s , w _s ) respectively For the second-order partial derivatives of u, v, w, iteratively solve until $1 / m Σ_{the s = 0}^{m - 1} | | d (u_{the s} + δ u, v_{the s} + δv, w_{the s} + δw) - d (u_{the s}, v_{the s}, w_{the s}) | | \leq ϵ,$ ε is the preset surface fitting accuracy, so as to obtain the final triangular Bezier surface S'(u, v, w);

第三步：在所述的三角Bezier曲面S’(u，v，w)上基于曲率的等参数取线，再在线上等参数取点，用于填补点云孔洞，在求曲率的过程中首先需对三维点云进行删格划分，该删格划分为构造一个点云数据的最小外接正方体，其两两垂直的3条边分别与笛卡儿坐标系的3个坐标轴平行，沿三个坐标轴方向划分成边长为L空间六面体立方删格，其次在拟合曲面的点P_s的27个邻近子包围盒内求出其k个邻近点，然后设拟合曲面的点P_s和其k个邻近点组成集合K(P_s)，S(P_s)为拟合曲面的点P_s的k个邻近点最小二乘拟合平面，令 P为拟合曲面的点P_s的k个邻近点集合K(P_s)的形心，称为拟合曲面的点P_s的中心点，该中心点为：Step 3: On the triangular Bezier surface S'(u, v, w), take a line based on the parameters such as curvature, and then take points on the line and other parameters to fill the hole in the point cloud. In the process of calculating the curvature Firstly, the 3D point cloud needs to be divided into grids. The grid is divided into the smallest circumscribed cube for constructing a point cloud data. The three vertical sides are respectively parallel to the three coordinate axes of the Cartesian coordinate system. The direction of each coordinate axis is divided into a hexahedral cubic grid whose side length is L, and then its k adjacent points are found in the 27 adjacent sub-bounding boxes of the point P _s of the fitted surface, and then the point P _s of the fitted surface is set and its k neighboring points to form a set K(P _s ), S(P _s ) is the least squares fitting plane of k neighboring points P _s of the fitting surface, let P be the point P _s of the fitting surface The centroid of the set of k adjacent points K(P _s ), is called the center point of the point P _s of the fitted surface, and the center point is:

$\overset{&OverBar; &OverBar;}{P P} = = \frac{11}{((k k + + 11))} \underset{{P P}_{s the s} &Element; &Element; K K (({P P}_{s the s}))}{Σ Σ} {P P}_{s the s}$

设拟合曲面的点P_s的第j个邻近点到最小二乘平面S(P_s)的距离为d_j，到 P的距离为λ_j，那么对拟合曲面的点P_s的第j个点存在函数f_j(P_s)，该函数f_j(P_s)为：Suppose the distance from the jth adjacent point of the point P _s of the fitted surface to the least square plane S(P _s ) is d _j , and the distance to P is λ _j , then for the jth adjacent point of the point P _s of the fitted surface There is a function f _j (P _s ) at each point, and the function f _j (P _s ) is:

${f f}_{j j} (({P P}_{s the s})) = = \frac{{d d}_{j j}}{{λ λ}_{j j}}$

那么拟合曲面的点P_s的曲率函数可以表示为Then the curvature function of the point P _s of the fitted surface can be expressed as

$f f (({P P}_{s the s})) = = \frac{11}{k k} {Σ Σ}_{j j = = 11}^{k k} {f f}_{j j} (({P P}_{s the s}))$

根据曲率函数f(P_s)求出，拟合曲面的点P_s的平均曲率ρ_s，同理求出整个点云的平均曲率ρ₀。令d为整个点云的平均点距，则取点间隔Δω＝ρ₀×d/ρ_s，在Bezier曲面片上取点时，首先在曲面一个参数方向上以Δω等间隔取等参数曲线，再对每一条等参数曲线在另一个参数方向上以Δω等间隔取点，求得填补孔洞的点。According to the curvature function f(P _s ), the average curvature ρ _s of the point P _s of the fitted surface is calculated, and the average curvature ρ ₀ of the entire point cloud is calculated similarly. Let d be the average point distance of the entire point cloud, then take the point interval Δω=ρ ₀ ×d/ρ _s , when taking points on a Bezier surface patch, firstly take isoparametric curves at equal intervals of Δω in one parameter direction of the surface, and then For each isoparametric curve, take points at equal intervals of Δω in the direction of another parameter, and obtain the point to fill the hole.

本发明主要用于对三维扫描系统中各种具有复杂曲面形状的点云原始模型中孔洞进行填补的应用场合。利用本发明中的曲面拟合算法，即可得到一张精确逼近孔洞周围散乱点的曲面，随后在面上根据曲率的变化取点可以即可实现对孔洞的光滑填补。该方法主要有以下优点：The invention is mainly used in the application occasion of filling holes in various point cloud original models with complex curved surface shapes in a three-dimensional scanning system. Using the surface fitting algorithm in the present invention, a curved surface that accurately approximates the scattered points around the hole can be obtained, and then the holes can be filled smoothly by taking points on the surface according to the change of curvature. This method mainly has the following advantages:

(1)由于传统的曲面插值方法对离散点的逼近效果并不好，本发明采用曲面拟合的方法能得到一张精确拟合孔洞周围散乱点集的曲面。(1) Since the traditional surface interpolation method does not have a good approximation effect on discrete points, the present invention adopts a surface fitting method to obtain a surface that accurately fits the scattered point sets around the hole.

(2)采用基于三角域上的曲面拟合方法，比四边域的曲面更容易表现无规则的表面，能保证一定的光顺性和保形性。(2) Using the surface fitting method based on the triangular domain, it is easier to express the irregular surface than the surface of the quadrilateral domain, and can guarantee a certain degree of smoothness and shape retention.

(3)相比较多数算法中的网格化处理，本发明主要针对三维扫描系统得到的原始点云模型，不需要任何的网格化处理，适用性广，速度快。(3) Compared with the grid processing in most algorithms, the present invention is mainly aimed at the original point cloud model obtained by the 3D scanning system, does not need any grid processing, and has wide applicability and high speed.

(4)拟合曲面时，交互选取孔洞周围的散乱点，所以当孔洞边界形状比较复杂时，也同样适用。(4) When fitting the surface, the scattered points around the hole are selected interactively, so it is also applicable when the shape of the hole boundary is complex.

(5)对散乱点参数化时，无需拟合孔洞的边界曲线，由于避免了孔洞的边界识别，因此对于绝大多数具有任意形状的孔洞也同样适用。(5) When parameterizing scattered points, there is no need to fit the boundary curve of the hole. Since the recognition of the boundary of the hole is avoided, it is also applicable to most holes with arbitrary shapes.

(6)由于对散乱点进行曲面拟合采用了初步拟合和迭代逼近优化两步方法，能够使拟合的曲面获得更高的精度。(6) Since the two-step method of preliminary fitting and iterative approximation optimization is adopted for surface fitting of scattered points, the fitted surface can obtain higher precision.

(7)可以根据实际需要，在曲面拟合过程中人为控制曲面的次数及拟合的精度，避免了不必要的计算冗余。(7) According to the actual needs, the number of times and the fitting accuracy of the surface can be artificially controlled during the surface fitting process, avoiding unnecessary redundant calculation.

(8)在面上取点时，考虑曲率变化，能实现对孔洞的光滑补偿，一般在点云中，点云密度大的地方，曲率变化也比较大，因此，基于曲率的等参数取点能光滑的填补孔洞。(8) When taking points on the surface, consider the curvature change to achieve smooth compensation for the hole. Generally, in the point cloud, where the point cloud density is high, the curvature change is also relatively large. Therefore, based on curvature and other parameters to select points Can fill holes smoothly.

(9)操作过程比较简单，只要交互选定孔洞周围的离散点，接下来的的步骤都可自动完成，速度快。并且该方法具有很强的通用性。(9) The operation process is relatively simple, as long as the discrete points around the hole are interactively selected, the next steps can be automatically completed, and the speed is fast. And the method has strong versatility.

此外要注意的是，选取散乱点时，尽量在孔洞周围大范围内选取，可以避免自适应增加散乱点数量，以进一步提高算法的效率。In addition, it should be noted that when selecting scattered points, try to select them in a wide range around the hole, which can avoid adaptively increasing the number of scattered points and further improve the efficiency of the algorithm.

附图说明Description of drawings

图1是逆向工程流程图。Figure 1 is a flow chart of reverse engineering.

图2是光栅式三维扫描系统组成图。Figure 2 is a composition diagram of the raster type three-dimensional scanning system.

图3是算法整体流程图。Figure 3 is the overall flowchart of the algorithm.

图4是选取圆盖点云中的孔洞附近的点集实例示意图。Fig. 4 is a schematic diagram of an example of selecting a point set near a hole in a dome point cloud.

图5散乱点参数化流程图。Figure 5 Flow chart of parameterization of scattered points.

图6是曲面初步拟合流程图。。Figure 6 is a flow chart of the preliminary fitting of the curved surface. .

图7是对圆盖点云中的孔洞附近的点集拟合曲面实例示意图。Fig. 7 is a schematic diagram of an example of fitting a surface to a point set near a hole in a dome point cloud.

图8是是曲面迭代逼近优化流程图Figure 8 is a flow chart of surface iterative approximation optimization

图9是对圆盖点云中的孔洞填补实例示意图。Fig. 9 is a schematic diagram of an example of hole filling in a dome point cloud.

图10是删格化图。Figure 10 is a gridded graph.

图11是曲率函数求解图。Fig. 11 is a solution diagram of the curvature function.

图12是面上取点填补孔洞流程图。Fig. 12 is a flow chart of filling holes by taking points on the surface.

具体实施方式Detailed ways

$\{\begin{matrix} f f ((u u,, v v,, w w)) = = d d ((u u,, v v,, w w)) \cdot \cdot {S S}_{u u} ((u u,, v v,, w w)) = = 00 \\ g g ((u u,, v v,, w w)) = = d d ((u u,, v v,, w w)) \cdot \cdot {S S}_{v v} ((u u,, v v,, w w)) = = 00 \\ h h ((u u,, v v,, w w)) = = d d ((u u,, v v,, w w)) \cdot \cdot {S S}_{w w} ((u u,, v v,, w w)) = = 00 \end{matrix}$

Hσ^T＝κ^T Hσ ^T =κ ^T

式中：In the formula:

$H h = = [\begin{matrix} {f f}_{u u} & {f f}_{v v} & {f f}_{w w} \\ {g g}_{u u} & {g g}_{v v} & {g g}_{w w} \\ {h h}_{u u} & {h h}_{v v} & {h h}_{w w} \end{matrix}] = = [\begin{matrix} {| | | | {S S}_{u u} | | | |}^{22} + + d d \cdot \cdot {S S}_{uu u u} & {S S}_{v v} \cdot \cdot {S S}_{u u} + + d d \cdot &Center Dot; {S S}_{uv uv} & {S S}_{w w} \cdot &Center Dot; {S S}_{u u} + + d d \cdot &Center Dot; {S S}_{uw uw} \\ {S S}_{u u} \cdot &Center Dot; {S S}_{v v} + + d d \cdot &Center Dot; {S S}_{vu v u} & {| | | | {S S}_{v v} | | | |}^{22} + + d d \cdot &Center Dot; {S S}_{vv vv} & {S S}_{w w} \cdot &Center Dot; {S S}_{v v} + + d d \cdot &Center Dot; {S S}_{vw vw} \\ {S S}_{u u} \cdot \cdot {S S}_{w w} + + d d \cdot &Center Dot; {S S}_{wu wu} & {S S}_{v v} \cdot &Center Dot; {S S}_{w w} + + d d \cdot \cdot {S S}_{wv wv} & {| | | | {S S}_{w w} | | | |}^{22} + + d d \cdot &Center Dot; {S S}_{ww ww} \end{matrix}]$

则可得：Then you can get:

$[\begin{matrix} {| | | | {S S}_{u u} | | | |}^{22} + + d d \cdot &Center Dot; {S S}_{uu u u} & {S S}_{v v} \cdot \cdot {S S}_{u u} + + d d \cdot &Center Dot; {S S}_{uv uv} & {S S}_{w w} \cdot \cdot {S S}_{u u} + + d d \cdot &Center Dot; {S S}_{uw uw} \\ {S S}_{u u} \cdot &Center Dot; {S S}_{v v} + + d d \cdot &Center Dot; {S S}_{vu v u} & {| | | | {S S}_{v v} | | | |}^{22} + + d d \cdot &Center Dot; {S S}_{vv vv} & {S S}_{w w} \cdot &Center Dot; {S S}_{v v} + + d d \cdot \cdot {S S}_{vw vw} \\ {S S}_{u u} \cdot &Center Dot; {S S}_{w w} + + d d \cdot \cdot {S S}_{wu wu} & {S S}_{v v} \cdot &Center Dot; {S S}_{w w} + + d d \cdot &Center Dot; {S S}_{wv wv} & {| | | | {S S}_{w w} | | | |}^{22} + + d d \cdot &Center Dot; {S S}_{ww ww} \end{matrix}] [\begin{matrix} δu δu \\ δv δv \\ δw δw \end{matrix}] = = - - [\begin{matrix} f f (({u u}_{s the s},, {v v}_{s the s},, {w w}_{s the s})) \\ g g (({u u}_{s the s},, {v v}_{s the s},, {w w}_{s the s})) \\ h h (({u u}_{s the s},, {v v}_{s the s},, {w w}_{s the s})) \end{matrix}]$

$[\begin{matrix} {| | | | {S S}_{u u} | | | |}^{22} - - {S S}_{w w} \cdot &Center Dot; {S S}_{u u} + + d d \cdot &Center Dot; (({S S}_{uu u u} - - {S S}_{uw uw})) & {S S}_{v v} \cdot &Center Dot; {S S}_{u u} - - {S S}_{w w} \cdot &Center Dot; {S S}_{u u} + + d d (({S S}_{uv uv} - - {S S}_{uw uw})) \\ {S S}_{u u} \cdot \cdot {S S}_{v v} - - {S S}_{w w} \cdot \cdot {S S}_{v v} + + d d (({S S}_{vu v u} - - {S S}_{vw vw})) & {| | | | {S S}_{v v} | | | |}^{22} - - {S S}_{w w} \cdot &Center Dot; {S S}_{v v} + + d d \cdot \cdot (({S S}_{vv vv} - - {S S}_{vw vw})) \\ {S S}_{u u} \cdot &Center Dot; {S S}_{w w} - - {| | | | {S S}_{u u} | | | |}^{22} + + d d \cdot \cdot (({S S}_{wu wu} - - {S S}_{ww ww})) & {S S}_{v v} \cdot &Center Dot; {S S}_{w w} - - {| | | | {S S}_{w w} | | | |}^{22} + + d d \cdot &Center Dot; (({S S}_{wv wv} - - {S S}_{ww ww})) \end{matrix}] [\begin{matrix} δu δu \\ δv δv \end{matrix}] = = - - [\begin{matrix} f f (({u u}_{s the s},, {v v}_{s the s},, {w w}_{s the s})) \\ g g (({u u}_{s the s},, {v v}_{s the s},, {w w}_{s the s})) \\ h h (({u u}_{s the s},, {v v}_{s the s},, {w w}_{s the s})) \end{matrix}]$

其中，δu，δv为u，v两个方向上的迭代步长。S_uu，S_uv，S_uw，S_vv，S_vu，S_vw，S_ww，S_wu，S_wv是曲面S(u，v，w)在点(u_s，v_s，w_s)处分别对u，v，w的二阶偏导数，迭代求解直到 $1 / m Σ_{s = 0}^{m - 1} | | d (u_{s} + δu, v_{s} + δv, w_{s} + δw) - d (u_{s}, v_{s}, w_{s}) | | \leq ϵ,$ ε为预置曲面拟合精度，从而得到最终确定的三角Bezier曲面S’(u，v，w)；Among them, δu and δv are the iteration step sizes in the two directions of u and v. S _uu , S _uv , S _uw , S _vv , S _vu , S _vw , S _ww , S _wu , S _wv are surfaces S(u, v, w) at points (u _s , v _s , w _s ) respectively For the second-order partial derivatives of u, v, w, iteratively solve until $1 / m Σ_{the s = 0}^{m - 1} | | d (u_{the s} + δu, v_{the s} + δ v, w_{the s} + δw) - d (u_{the s}, v_{the s}, w_{the s}) | | \leq ϵ,$ ε is the preset surface fitting accuracy, so as to obtain the final triangular Bezier surface S'(u, v, w);

第三步：在所述的三角Bezier曲面S’(u，v，w)上基于曲率的等参数取线，再在线上等参数取点，用于填补点云孔洞，在求曲率的过程中首先需对三维点云进行删格划分，该删格划分为构造一个点云数据的最小外接正方体，其两两垂直的3条边分别与笛卡儿坐标系的3个坐标轴平行，沿三个坐标轴方向划分成边长为L空间六面体立方删格，其次在拟合曲面的点P_s的27个邻近子包围盒内求出其k个邻近点，然后设拟合曲面的点P_s和其k个邻近点组成集合K(P_s)，S(P_s)为拟合曲面的点P_s的k个邻近点最小二乘拟合平面，令 P为拟合曲面的点P_s的k个邻近点集合K(P_s)的形心，称为拟合曲面的点P_s的中心点，该中心点为：Step 3: On the triangular Bezier surface S'(u, v, w), take a line based on the parameters such as curvature, and then take points on the line and other parameters to fill the hole in the point cloud. In the process of calculating the curvature Firstly, the 3D point cloud needs to be divided into grids. The grid is divided into the smallest circumscribed cube for constructing a point cloud data. The three vertical sides are respectively parallel to the three coordinate axes of the Cartesian coordinate system. The direction of each coordinate axis is divided into a hexahedral cubic grid whose side length is L, and then its k adjacent points are found in the 27 adjacent sub-bounding boxes of the point P _s of the fitted surface, and then the point P _s of the fitted surface is set and its k neighboring points to form a set K(P _s ), S(P _s ) is the least squares fitting plane of k neighboring points P _s of the fitting surface, let P be the point P _s of the fitting surface The centroid of the set of k adjacent points K(P _s ), is called the center point of the point P _s of the fitting surface, and the center point is:

${f f}_{j j} (({P P}_{s the s})) = = \frac{{d d}_{j j}}{{λ λ}_{j j}}$

下面参照附图，对本发明加以详细描述：Below with reference to accompanying drawing, the present invention is described in detail:

在逆向工程中，面对的是密集的点云散乱数据。拟合曲面时，如果曲面的对象边界和形状极其复杂时，一般不便直接运用常规的曲面构造方法，而Bezier曲面由于其构造灵活，边界适应性好，具有构造复杂形状的潜力，而且其本身具有良好的性质：光滑性、局部性和保形性。基于三角域曲面的拟合方法最适合表现无规则型面的物体，特别是人面，地貌等自然物体以及玩具等产品，基于四边域参数曲面的拟合方法，通常要求数据点是有序的，这个条件比较苛刻，综合以上，因此本发明采用三角域上的Bezier曲面对孔洞部分进行拟合填充。总体算法流程见附图3。In reverse engineering, what is faced is dense point cloud scattered data. When fitting a surface, if the object boundary and shape of the surface are extremely complex, it is generally inconvenient to use the conventional surface construction method directly. However, the Bezier surface has the potential to construct complex shapes due to its flexible structure and good boundary adaptability. Good properties: smoothness, locality and shape retention. The fitting method based on triangular domain surfaces is most suitable for objects that represent irregular surfaces, especially natural objects such as human faces and landforms, and toys and other products. The fitting method based on quadrilateral domain parametric surfaces usually requires that the data points be in order , this condition is relatively harsh, and based on the above, the present invention uses the Bezier surface on the triangular domain to fit and fill the hole part. The overall algorithm flow is shown in Figure 3.

本发明主要涉及以下四方面的内容：The present invention mainly relates to the following four aspects:

1)散乱点的选定及参数化过程1) Selection and parameterization process of scattered points

在圆盖点云孔洞屏幕二维坐标系附近选择不共线的三点A、B、C，构成需要初步拟合的三角Bezier曲面片的三个角点，将孔洞周围的点向这三点所构成的三角形平面投影，如果点落在三角形平面内，则认为是拟合曲面时所需的点P_s(s＝0，1，…，m-1)，见附图4。Select three non-collinear points A, B, and C near the two-dimensional coordinate system of the hole screen of the dome point cloud to form the three corner points of the triangular Bezier surface patch that needs to be initially fitted, and move the points around the hole to these three points For the formed triangular plane projection, if the point falls within the triangular plane, it is considered as the point P _s (s=0, 1, . . . , m-1) required for fitting the curved surface.

根据选取的点P_s在三角形ABC平面的投影位置计算其曲面参数化坐标。即如果P_s在三角形平面ABC内相对应的投影点为P_s’，那么P_s的参数化坐标(u_s，v_s，w_s)分别为u_s＝area(ΔAP_s’B)/area(ΔABC)、v_s＝area(ΔAP_s’C)/area(ΔABC)、w_s＝area(ΔBP_s’C)/area(ΔABC)。算法流程图见附图5。According to the projected position of the selected point P _s on the triangle ABC plane, its surface parameterized coordinates are calculated. That is, if the corresponding projection point of P _s in the triangular plane ABC is P _s ', then the parameterized coordinates (u _s , v _s , w _s ) of P _s are respectively u _s = area(ΔAP _s 'B)/area (ΔABC), v _s =area(ΔAP _s 'C)/area(ΔABC), w _s =area( _ΔBP s'C)/area(ΔABC). The flowchart of the algorithm is shown in Figure 5.

2)最小二乘自适应初步拟合2) Least square adaptive preliminary fitting

三角域ABC上n次Bezier曲面的表达式为：The expression of n degree Bezier surface on triangular domain ABC is:

$S S ((u u,, v v,, w w)) = = {Σ Σ}_{i i = = 00}^{n no} {Σ Σ}_{j j = = 00}^{n no - - i i} {B B}_{ijk ijk}^{n no} ((u u,, v v,, w w)) {b b}_{ijk ijk}$

式中 $B_{ijk}^{n} (u, v, w) = \frac{n!}{i! j! k!} u^{i} v^{j} w^{k} &GreaterEqual; 0 \{\begin{matrix} 0 \leq u, v, w \leq 1, u + v + w = 1 \\ 0 \leq i, j, k \leq n, i + j + k = n \end{matrix}$ In the formula $B_{ijk}^{no} (u, v, w) = \frac{no!}{i! j! k!} u^{i} v^{j} w^{k} &Greater Equal; 0 \{\begin{matrix} 0 \leq u, v, w \leq 1, u + v + w = 1 \\ 0 \leq i, j, k \leq no, i + j + k = no \end{matrix}$

三角域上n次Bezier曲面的表达式中未知量为P_s的参数化坐标(u_s，v_s，w_s)，以及(n+1)(n+2)/2个控制顶点，(u_s，v_s，w_s)根据上步已知，所以只要求出(n+1)(n+2)/2个控制顶点坐标b_ijk(0≤i，j，k≤n，i+j+k＝n)，即可确定出一个三角Bezier曲面。The parameterized coordinates (u _s , v _s , w _s ) of the unknown quantity in the expression of the Bezier surface of _degree n on the triangular domain, and (n+1)(n+2)/2 control vertices, (u _s , v _s , w _s ) are known according to the previous step, so only (n+1)(n+2)/2 control vertex coordinates b _ijk (0≤i, j, k≤n, i+j) are required +k=n), a triangular Bezier surface can be determined.

根据三角Bezier曲面的性质可知三角Bezier曲面的三个角点正好是其控制多面体网格的三个控制顶点，那么所选三角形的三个顶点A、B、C即可做为三个已知控制顶点b_00n，b_0n0，b_n00，把P_s(s＝0，1，…，m-1)的坐标代入三角Bezier曲面方程，构成以(n+1)(n+2)/2-3个未知控制顶点为未知量的m个线性方程组。在所选三角形内部的三维坐标点数m一般大于(n+1)(n+2)/2-3，在不满足此条件即待求方程组系数矩阵奇异时，以三角形ABC的重心为基点，自适应扩大三角形ABC的面积来增加m的个数。最后采用最小二乘法求解此过约束方程组，确定出曲面其它的未知控制顶点，得到初步拟合的曲面S(u，v，w)。算法流程图见附图6。According to the properties of the triangular Bezier surface, the three corner points of the triangular Bezier surface are just the three control vertices of its control polyhedral mesh, then the three vertices A, B, and C of the selected triangle can be used as the three known control vertices Vertices b _00n , b _0n0 , b _n00 , put the coordinates of P _s (s=0, 1,..., m-1) into the triangular Bezier surface equation, and form the equation of (n+1)(n+2)/2-3 A system of m linear equations whose unknown control vertices are unknown quantities. The number of three-dimensional coordinate points m inside the selected triangle is generally greater than (n+1)(n+2)/2-3. When this condition is not satisfied, that is, when the coefficient matrix of the system of equations to be obtained is singular, the center of gravity of the triangle ABC is used as the base point. Adaptively expand the area of triangle ABC to increase the number of m. Finally, the least square method is used to solve the over-constrained equations, and other unknown control vertices of the surface are determined, and the preliminary fitting surface S(u, v, w) is obtained. The flowchart of the algorithm is shown in Figure 6.

3)曲面逼近优化3) Surface approximation optimization

经过最小二乘初步拟合后，散乱点集P_s与曲面S(u，v，w)的距离误差还比较大，因此需要进一步修正控制顶点，改变控制多面体的形状，从而改变曲面的形状，使曲面逼近点集。本发明采用了以牛顿迭代法为基础的迭代逼近算法，提高拟合精度，最终确定出曲面S’(u，v，w)，见附图7。After preliminary fitting by least squares, the distance error between the scattered point set P _s and the surface S(u, v, w) is relatively large, so it is necessary to further correct the control vertices and change the shape of the control polyhedron, thereby changing the shape of the surface. Make the surface approximate the point set. The present invention adopts an iterative approximation algorithm based on Newton's iterative method to improve the fitting accuracy and finally determine the curved surface S'(u, v, w), see Figure 7.

具体步骤如下：Specific steps are as follows:

(1)求出各个散乱点到曲面的距离向量d(u，v，w)及曲面分别对对应点参数方向的偏微分S_u(u，v，w)，S_v(u，v，w)，S_w(u，v，w)，组成待求方程组如下：(1) Calculate the distance vector d(u, v, w) from each scattered point to the surface and the partial differential of the surface to the parameter direction of the corresponding point S _u (u, v, w), S _v (u, v, w ), S _w (u, v, w), form the system of equations to be obtained as follows:

(2)采用基于牛顿迭代法的方法求解此方程组，最终可求得的迭代方程组为：(2) Using the method based on Newton's iterative method to solve this system of equations, the final iterative system of equations that can be obtained is:

$[\begin{matrix} {| | | | {S S}_{u u} | | | |}^{22} - - {S S}_{w w} \cdot \cdot {S S}_{u u} + + d d \cdot \cdot (({S S}_{uu u u} - - {S S}_{uw uw})) & {S S}_{v v} \cdot \cdot {S S}_{u u} - - {S S}_{w w} \cdot \cdot {S S}_{u u} + + d d \cdot \cdot (({S S}_{uv uv} - - {S S}_{uw uw})) \\ {S S}_{u u} \cdot \cdot {S S}_{v v} - - {S S}_{w w} \cdot \cdot {S S}_{v v} + + d d \cdot \cdot (({S S}_{vu v u} - - {S S}_{vw vw})) & {| | | | {S S}_{v v} | | | |}^{22} - - {S S}_{w w} \cdot \cdot {S S}_{v v} + + d d \cdot &Center Dot; (({S S}_{vv vv} - - {S S}_{vw vw})) \\ {S S}_{u u} \cdot \cdot {S S}_{w w} - - {| | | | {S S}_{w w} | | | |}^{22} + + d d \cdot \cdot (({S S}_{wu wu} - - {S S}_{ww ww})) & {S S}_{v v} \cdot \cdot {S S}_{w w} - - {| | | | {S S}_{w w} | | | |}^{22} + + d d \cdot \cdot (({S S}_{wv wv} - - {S S}_{ww ww})) \end{matrix}] [\begin{matrix} δu δu \\ δv δ v \end{matrix}] = = - - [\begin{matrix} f f (({u u}_{s the s},, {v v}_{s the s},, {w w}_{s the s})) \\ g g (({u u}_{s the s},, {v v}_{s the s},, {w w}_{s the s})) \\ h h (({u u}_{s the s},, {v v}_{s the s},, {w w}_{s the s})) \end{matrix}]$

(3)迭代求解以上方程组，当拟合曲面的点数少于迭代方程组的未知数个数时，及方程组系数奇异时，同样以三角形ABC的重心为基点，自适应扩大三角ABC的面积来增加m的个数，重复(2)。(3) Iteratively solve the above equations. When the number of points of the fitted surface is less than the number of unknowns of the iterative equations, and the coefficients of the equations are singular, the center of gravity of the triangle ABC is also used as the base point to adaptively expand the area of the triangle ABC. Increase the number of m, repeat (2).

算法流程图见附图8。The flowchart of the algorithm is shown in Figure 8.

4)面上取点4) Take points on the surface

在所述的三角Bezier曲面S’(u，v，w)上基于曲率的等参数取线，再在线上等参数取点，用于填补点云孔洞，见附图9。在求曲率的过程中首先需对三维点云进行删格划分，该删格划分为构造一个点云数据的最小外接正方体，其两两垂直的3条边分别与笛卡儿坐标系的3个坐标轴平行，沿三个坐标轴方向划分成边长为L空间六面体立方删格，见附图10，其次在拟合曲面的点P_s的27个邻近子包围盒内求出其k个邻近点，然后设拟合曲面的点P_s和其k个邻近点组成集合K(P_s)，S(P_s)为拟合曲面的点P_s的k个邻近点最小二乘拟合平面，令P为拟合曲面的点P_s的k个邻近点集合K(P_s)的形心，称为拟合曲面的点P_s的中心点，该中心点为：On the triangular Bezier surface S'(u, v, w), take a line based on the curvature isoparameters, and then take a point on the line and other parameters to fill the hole in the point cloud, see Figure 9. In the process of calculating the curvature, it is first necessary to divide the three-dimensional point cloud into grids. The grid is divided into the smallest circumscribed cube for constructing a point cloud data. The coordinate axes are parallel, and along the direction of the three coordinate axes, it is divided into a _hexahedron cubic grid whose side length is L. point, and then set the point P _s of the fitting surface and its k neighboring points to form a set K(P _s ), and S(P _s ) is the least squares fitting plane of the k neighboring points of the point P _s of the fitting surface, Let P be the centroid of the k adjacent point set K(P _s ) of the point P _s of the fitted surface, called the center point of the point P _s of the fitted surface, which is:

设拟合曲面的点P_s的第j个邻近点到最小二乘平面S(P_s)的距离为d_j，到 P的距离为λ_j，见附图11，那么对拟合曲面的点P_s的第j个点存在函数f_j(P_s)，该函数f_j(P_s)为：Suppose the distance from the jth adjacent point of the point P _s of the fitting surface to the least square plane S(P _s ) is d _j , and the distance to P is λ _j , see attached drawing 11, then for the point of the fitting surface There is a function f _j (P _s ) at the jth point of P _s , and the function f _j (P _s ) is:

${f f}_{j j} (({P P}_{s the s})) = = \frac{{d d}_{j j}}{{λ λ}_{j j}}$

算法流程图见附图12。The flowchart of the algorithm is shown in Figure 12.

Claims

1, a kind of point cloud hole filling method of three-dimensional scanning, it is characterized in that:

Step 1: Set a triangle ABC around the point cloud hole and in the screen coordinate plane. The area of the triangle ABC can make the projection of the point cloud hole and its surrounding points to the screen coordinate plane fall into the triangle ABC, and The point that the projection falls into the triangle ABC is used as the point P _s (s=0, 1, ..., m-1) of the fitted surface when filling the hole, according to the projection P _s of the point P _s of the fitted surface on the triangle ABC plane 'Position is calculated by its surface parameterized coordinates (u _s , v _s , w _s ), u _s = (ΔAP _s 'B area)/(ΔABC area), v _s = (ΔAP _s 'C area)/(ΔABC area) , w _s =(ΔBP _s 'C area)/(ΔABC area), the coordinates of the point P _s (s=0,1,...,m-1) of the fitted surface and its surface parameterized coordinates (u _s , v _s , w _s ) are substituted into the nth degree Bezier surface equation and the control points of the nth degree Bezier surface are obtained by the least square method, so as to obtain the preliminary fitting surface S(u, v, w);

Step 2: Calculate the distance vector d(u, v, w) from the point P _s (s=0, 1, ..., m-1) of each fitted surface to the surface and the deviation of the surface to the parameter direction of the corresponding point Differentiate S _u (u, v, w), S _v (u, v, w), S _w (u, v, w), let:

\{\begin{matrix} f f ((u u,, v v,, w w)) = = d d ((u u,, v v,, w w)) \cdot &Center Dot; {S S}_{u u} ((u u,, v v,, w w)) = = 00 \\ g g ((u u,, v v,, w w)) = = d d ((u u,, v v,, w w)) \cdot &Center Dot; {S S}_{v v} ((u u,, v v,, w w)) = = 00 \\ h h ((u u,, v v,, w w)) = = d d ((u u,, v v,, w w)) \cdot &Center Dot; {S S}_{w w} ((u u,, v v,, w w)) = = 00 \end{matrix}

According to Newton's iterative method, the following iterative equations are obtained:

[\begin{matrix} {| | | | {S S}_{u u} | | | |}^{22} - - {S S}_{w w} \cdot \cdot {S S}_{u u} + + d d \cdot &Center Dot; (({S S}_{uu u u} - - {S S}_{uw uw})) {S S}_{v v} \cdot &Center Dot; {S S}_{u u} - - {S S}_{w w} \cdot &Center Dot; {S S}_{u u} + + d d \cdot &Center Dot; (({S S}_{uv uv} - - {S S}_{uw uw})) \\ {S S}_{u u} \cdot &Center Dot; {S S}_{v v} - - {S S}_{w w} \cdot &Center Dot; {S S}_{v v} + + d d \cdot &Center Dot; (({S S}_{vu v u} - - {S S}_{uw uw})) {| | | | {S S}_{v v} | | | |}^{22} - - {S S}_{w w} \cdot \cdot {S S}_{v v} + + d d \cdot \cdot (({S S}_{vv vv} - - {S S}_{vw vw})) \\ {S S}_{u u} \cdot \cdot {S S}_{w w} - - {| | | | {S S}_{w w} | | | |}^{22} + + d d \cdot \cdot (({S S}_{wu wu} - - {S S}_{ww ww})) {S S}_{v v} \cdot \cdot {S S}_{w w} - - {| | | | {S S}_{w w} | | | |}^{22} + + d d \cdot \cdot (({S S}_{wv wv} - - {S S}_{ww ww})) \end{matrix}] [\begin{matrix} δu δ u \\ δv δ v \end{matrix}] = = - - [\begin{matrix} f f (({u u}_{s the s},, {v v}_{s the s},, {w w}_{s the s})) \\ g g (({u u}_{s the s},, {v v}_{s the s},, {w w}_{s the s})) \\ h h (({u u}_{s the s},, {v v}_{s the s},, {w w}_{s the s})) \end{matrix}]

Among them, δu and δv are the iterative step sizes in the directions of u and v, and S _uu , S _uv , S _uw , S _vv , S _vu , S _vw , S _ww , S _wu , and S _wv are surface S(u, v, w) at the point (u _s , v _s , w _s ) with respect to the second order partial derivatives of u, v, w respectively,

iteratively solve until

1 / m Σ_{the s = 0}^{m - 1} | | d (u_{the s} + δ u, v_{the s} + δv, w_{the s} + δw) - d (u_{the s}, v_{the s}, w_{the s}) | | \leq ϵ,

ε is the preset surface fitting accuracy, so as to obtain the final triangular Bezier surface S'(u, v, w);

Step 3: On the triangular Bezier surface S'(u, v, w), take a line based on the parameters such as curvature, and then take points on the line and other parameters to fill the hole in the point cloud. In the process of calculating the curvature Firstly, the 3D point cloud needs to be divided into grids. The grid is divided into the smallest circumscribed cube for constructing a point cloud data. The three vertical sides are respectively parallel to the three coordinate axes of the Cartesian coordinate system. The direction of each coordinate axis is divided into a hexahedral cubic grid whose side length is L, and then its k adjacent points are found in the 27 adjacent sub-bounding boxes of the point P _s of the fitted surface, and then the point P _s of the fitted surface is set and its k neighboring points to form a set K(P _s ), S(P _s ) is the least squares fitting plane of k neighboring points P _s of the fitting surface, let P be the point P _s of the fitting surface The centroid of the set of k adjacent points K(P _s ), is called the center point of the point P _s of the fitted surface, and the center point is:

\overset{&OverBar; &OverBar;}{P P} = = \frac{11}{((k k + + 11))} \underset{{P P}_{s the s} &Element; &Element; K K (({P P}_{s the s}))}{Σ Σ} {P P}_{s the s}

Suppose the distance from the jth adjacent point of the point P _s of the fitted surface to the least square plane S(P _s ) is d _j , and the distance to P is λ _j , then for the jth adjacent point of the point P _s of the fitted surface There is a function f _j (P _s ) at each point, and the function f _j (P _s ) is:

{f f}_{j j} (({P P}_{s the s})) = = \frac{{d d}_{j j}}{{λ λ}_{j j}}

Then the curvature function of the point P _s of the fitted surface can be expressed as

f f (({P P}_{s the s})) = = \frac{11}{k k} {Σ Σ}_{j j = = 11}^{k k} {f f}_{j j} (({P P}_{s the s}))

According to the curvature function f(P _s ), the average curvature of the point P _s of the fitted surface is obtained In the same way, find the average curvature of the entire point cloud

Let d be the average point distance of the entire point cloud, then take the point interval

Δω = \overset{&OverBar;}{ρ_{0}} \times \overset{&OverBar;}{d} / \overset{&OverBar;}{ρ_{the s}},

When taking points on a Bezier surface patch, first take isoparametric curves at equal intervals of Δω in one parameter direction of the surface, and then take points at equal intervals of Δω for each isoparametric curve in the direction of another parameter to obtain the point to fill the hole .

2. The point cloud hole filling method of three-dimensional scanning according to claim 1, wherein when the number of points of the fitting surface is less than the number of control vertices of the surface, the area of the set triangle ABC is enlarged so that the number of points of the fitting surface is greater than The surface controls the number of vertices.

3. The point cloud hole filling method of three-dimensional scanning according to claim 1, wherein when the number of points of the fitted surface is less than the number of unknowns of the iterative equation system, the area of the set triangle ABC is enlarged so that the fitted surface The number of points is greater than or equal to the number of unknowns in the iterative equation system.