CN113610726B - 一种区域限定的改进三次b样条插值修补方法 - Google Patents

Abstract

本发明公开了一种区域限定的改进三次B样条插值修补方法，包括步骤：获取至少一侧带有间隙开放孔洞的牲畜全局三维点云；寻找中轴平面，将牲畜全局三维点云左右侧对称划分左右两侧点云，对中轴平面进行拟合；将带有间隙开放孔洞的一侧点云投影到XOY平面中，得到投影后的点集，以块切片方式一侧点云上孔洞的上边界、下边界；计算出全局三维点云的孔洞的上边界U和下边界D；以块切片方式得到腹部点云的下边界，并计算出腹部最低点Down，在孔洞的上边界U的基础上加h得到Top，从而将待修补区域确定为[Down,Top]；将待修补区域点云沿X轴纵向切片，将切片投影到YOZ平面上，对投影点云进行有序化处理，设置插值点密度，采用三次B样条曲线插值方法完成修补。

Description

一种区域限定的改进三次B样条插值修补方法

技术领域

本发明涉及计算机图像修复技术领域，更具体的，涉及一种区域限定的改进三次B样条插值修补方法。

背景技术

表型值选择是猪育种过程中常用的手段，体尺是一组描述特定体重阶段猪体外在特征特性、反映猪品种特征特性及反应某一阶段生长发育情况的量化指标。长期以来猪体尺性状由人工测量完成，测量过程容易造成牲畜驱赶应激、人工测量准确性差、且难以度量。随着猪育种和养殖规模越来越大，其生产效率越来越低，因此迫切需要建立更为高效、准确且无接触的智能方案。

目前，活体猪三维点云采集主要有固定和移动两种方式，固定式采集设备将牲畜限定在栏杆通道中自由行走状态获取点云；移动式采用手持等方式将采集设备靠近牲畜调整角度自动完成点云采集。固定式采集通常有带栏杆的通道约束牲畜活动范围，防止破坏采集设备。而栏杆会形成阻挡给采集牲畜点云带来大片空洞，点云被分割成多个部分，破坏其完整性。有大片空洞点云将影响牲畜胸围、臀围、腹围的测量以及点云分割和骨架提取等。点云孔洞既影响模型可视化，也影响点云的后期处理。

由于点云孔洞类型，针对不同类型孔洞的修复，国内外学者进行了大量的研究。目前基于散乱点云孔洞修补的方法主要有：基于投影特征平面的孔洞修补算法，即将整个点云数据都参数化或投影到平面上进行处理。该方法适用于小尺度孔洞以及曲率特征变化不大、形状规律的曲面修补，在对细节丰富的表面的孔洞修补中存表现不佳。如中国专利公开号：CN107464223A，公开日：2017-12-12，公开了一种基于切片的点云孔洞修补方法，具体按照以下步骤实施：步骤1：对点云模型按照一定方向切割并投影，将此切割得到的点集投影在与其上下切割面平行的中间平面上；步骤2：采用聚类算法对步骤1切割投影得到的带状点集抽取其“核心骨架”；步骤3：对步骤2聚类后的质心排序；步骤4：确定孔洞边界点和邻接点；步骤5：获取切割层修补点，实现对单层孔洞的修补；步骤6：获取整体修补点，实现三维孔洞模型的完整修补。

目前还有三维空间直接修补孔洞方法，如采用经纬网格修补点云算法；移动最小二乘法插值法对孔洞区域进行曲面拟合完成填充和修复；双向波前法进行岛屿孔洞的修补。三角网格孔洞修补算法，改算法对孔洞区域曲面实现三角化，形成孔洞补丁，将补丁缝合在模型孔洞区域，实现孔洞修补。该算法适用于单连通区域的缺陷孔洞修复，对缺失区域复杂的拓扑结构，其算法时间复杂度较高。然而现有孔洞修补算法集中于对尺寸小、边界光滑以及曲率变化不大的孔洞修补。

发明内容

本发明为了解决以上现有技术中存在不足的问题，提供了一种区域限定的改进三次B样条插值修补方法。

为实现上述本发明目的，采用的技术方案如下：

一种区域限定的改进三次B样条插值修补方法，所述的方法包括步骤如下：

S1：获取至少一侧带有间隙开放孔洞的牲畜全局三维点云；

S2：寻找中轴平面，将牲畜全局三维点云左右侧对称划分为左侧点云和右侧点云，采用最小二乘法对中轴平面进行拟合；

S3：将带有间隙开放孔洞的一侧点云投影到XOY平面中，得到投影后的点集，以块切片方式得到带有间隙开放孔洞的一侧点云上孔洞的上边界L_upper、下边界L_lower；计算出全局三维点云的孔洞的上边界U和下边界D；

S4：以块切片方式得到腹部点云的下边界，并计算出腹部最低点Down，在孔洞的上边界U的基础上加h得到Top，从而将待修补区域确定为[Down,Top]，其中Down为下边界，Top为上边界；

S5：将待修补区域点云沿X轴纵向切片，将切片投影到YOZ平面上，对投影点云进行有序化处理，设置插值点密度，采用三次B样条曲线插值方法完成修补。

优选地，还包括步骤S6，为了便于后序测量工作且利于数据存储，采用体素滤波下采样算法对修补后的点云进行下采样处理。

进一步地，步骤S2，设全局三维点云数据为P＝{p₁，p₂,p₃，…，p_n}，其中点云p_i＝(x_i,y_i,z_i)，平面的一般方程：

z＝a₀x+a₁y+a₂ (2)

式中，a₀、a₁表示系数，a₂表示常数；

对于点云p_i(x_i，y_i,z_i)，则使得S表达式取得最小值：

若满足S最小，则满足式(4)

有

通过求解式(5)方程组，得到a₀，a₁,a₂的值，从而得到平面方程z＝a₀x_i+a₁y+a₂，即为所拟合的中轴平面。

再进一步地，步骤S3中，设带有间隙开放孔洞的一侧点云数据为L，点云个数为m，L＝{l₁,l₂,l₃,……,l_m},l_i＝(x_i,y_i,z_i)，作XOY平面投影，投影后点集为L_p＝{lp₁,lp₂,lp₃，…，lp_n}，lp_i＝(x_i,y_i，0)；

沿X轴方向切片，切片数量为B，切片大小slice；

切片j点集表示为S_j＝{s₁,s₂,s₃,……,s_n},s_k＝(x_k,y_k,0),x_k∈[min(x_i)+(j-1)*slice，min(x_i)+j*slice]，j∈(1，2，3，…，B)，k∈(1,2,3,…,n)；

将S_j中数据点根据y值进行排序，使得

y_k+1>y_k，求出max(y_k+1-y_k)，当y_t+1-y_t＝max(y_k+1-y_k)时，[y_t+1,y_t]即为当前切片中孔洞区域，记录为[b_j1,b_j2]；

从而计算得出一侧点云的孔洞上边界L_upper，下边界L_lower

再进一步地，根据孔洞的上边界L_upper、下边界L_lower计算出全局三维点云的孔洞的上边界U和下边界D，具体如下：

再进一步地，步骤S4，设一侧点云的腹部点云的下边界为L_a，设置切片数量B₂＝4，根据式(6)计算切片大小slice，选取第三片切片所处的X轴区间为[min(x_i)+2*slice,min(x_i)+3*slice]，求出第三片切片范围内的点的最小y值为L_a＝min(y_i)，根据式(9)计算出腹部的最低点Down如下：

式中，R_a表示另一侧点云的腹部点云的下边界。

再进一步地，步骤S5，对于修补区域点云H＝{h₁,h₂,h₃,…,h_a},h_i＝(x_i,y_i,z_i)，切片的数量为Blocks，沿着X轴方向切片，根据式(10)得出切片大小为slice

式中，i∈(1,2,3,…,a)；

将切片投影到YOZ平面上，其中第j切片为Q_j＝{q₁,q₂,q₃,……,q_b},q_i＝(x_i,y_i,z_i)。

再进一步地，对投影点云进行有序化处理，具体如下：因待修补点云切片投影从左到右呈现先递减后递增的趋势，所以以中轴平面为分界，对切片投影进行左右划分，将划分后的左侧点云降序排序，右侧点云升序排序，使得在左侧点云中，

y_i+1＜y_i，而右侧点云中，

y_i+1>y_i，由此完成有序化处理。

再进一步地，设定期望点云之间最小间距d，对于所拟合的三次曲线g，通过四个控制点P₀＝f₀(0,y₀,z₀),P₁＝f₁(0,y₁,z₁),P₂＝f₂(0,y₂,z₂),P₃＝f₃(0，y₃，z₃)以及t＝0,t＝1计算出起初其起点g₀，终点g₁：

采用平均间隔的方式在拟合曲线采集数据点进行曲线弧长的估算，设置t＝0.1求得进行采点曲线上10个点，计算相邻点之间线段长度并求和得到Length，以此来替代所拟合曲线的弧长

通过式(16)计算插值密度Density：

式中，d d表示期望的点云之间最小间距、Length表示曲线

的近似长度；

三次B样条曲线的基函数：

将三次B样条曲线的基函数代入到B样条曲线的表达式中，得到三次B样条曲线方程

P(t)＝P₀*F_0，3(t)+P₁*F_1,3(t)+P₂*F_2,3(t)+P₃*F_3,3(t)#(14)

式中，t表示修补的使用三次B样条曲线修补的参数；P₀,P₁,P₂,P₃分别代表四个控制点；

修补后点云数据集H表示如下：

式中，t_j表示根据密度Density计算出的值，P(t_j)表示根据三次B样条曲线计算的修补值。

再进一步地，步骤S6，采用体素滤波下采样算法对修补后的点云进行下采样处理，具体的，对于修补后的点云数据集H＝{h₁，h₂，h₃,…，h_m}，h_i＝(x_i，y_i，z_i)，设定体素包围盒的长、宽、高分别为length，width，high，并分别对应XYZ三个维度；

对于任意体素包围盒V_j，对任意h_i∈V_j，i≤m，包围盒V_j中有且仅有代表点v_j(x，y，z)

式中，k＝1、2…m。

本发明的有益效果如下：

本发明提出区域限定的改进三次B样条近似拟合修补方法，比传统的三次B样条近似拟合修补得到的点云数据均匀性和逼近性能更好。通过对不同四肢动物的点云孔洞模型数据和养殖场活猪点云孔洞数据进行腹部修补，通过分析修补结果可见本发明更具有稳健性和鲁棒性。

附图说明

图1是本实施例改进三次B样条曲线孔洞修补方法的实现流程。

图2是本实施例相机视角与世界坐标系示意图。

图3是本实施例待配准平面示意图。

图4是本实施例带孔洞的全局三维点云的猪体投影示意图。

图5是本实施例拟合中轴平面的示意图。

图6是本实施例切片投影的示意图。

图7是本实施例修补孔洞区域的示意图。

图8是本实施例待修补点云切片的投影示意图。

图9是本实施例待修补点云限定的修补范围的示意图。

图10是本实施例待修补点云有序化后修补范围的示意图。

图11是本实施例三次B样条曲线的示意图。

图12是本实施例修补后修补区域投影的示意图。

图13是本实施例修补后点云完整切片投影的示意图。

图14是本实施例模型点云数据示意图。

图15是本实施例测量位置及其切片示意图。

图16是本实施例各修补算法对猪体模型修补效果示意图。

图17是本实施例各修补算法估计值与人工测量值误差变化。

图18是本实施例极坐标变换修补算法腹围估计值与人工测量值的比较结果。

图19是本实施例曲面拟合算法腹围估计值与人工测量值的比较结果。

图20是本实施例三角化修补算法腹围估计值与人工测量值的比较结果。

图21是本实施例本文算法腹围估计值与人工测量值的比较结果。

图22是本实施例实验所用点云数据的三维效果。

图23是本实施例极坐标变换修补算法对实验所用点云数据处理的结果。

图24是本实施例曲面拟合修补算法对实验所用点云数据处理的结果。

图25是本实施例三角化修补算法对实验所用点云数据处理的结果。

图26是本实施例所述的方法对实验所用点云数据处理的结果。

具体实施方式

下面结合附图和具体实施方式对本发明做详细描述。

实施例1

如图1所示，一种区域限定的改进三次B样条插值修补方法，所述的方法包括步骤如下：

步骤S1：获取至少一侧带有间隙开放孔洞的牲畜全局三维点云。所述的牲畜包括猪、牛、羊、马、狗等四肢类动物，本实施例以活体猪为例进行详细说明。

本实施例的实验数据为养殖场采集自由行走状态下活体猪点云以及标准四肢类动物点云数据集。养殖场通过固定式采集设备将牲畜限定在栏杆通道中自由行走状态获取点云，具体如下：采集通道有3个KinectV2深度相机分别从左、右、上3个方向获取不同视角局部点云，通过点云去噪、配准融合等预处理。

所述的点云去噪具体如下：猪体点云噪声点通常分为两类，一类是远离目标点云的远距离离群点，其分布往往较为稀疏，表现为三维空间中一小团聚在一起而又随机分布的孤立点，采用基于距离的区域增长方法能够滤除远距离离群点；另一类是靠近目标点云的近距离离群点，对于近距离离群点的消除采用非迭代双阈值去噪算法。

所述的配准融合具体如下：由于3个深度相机从上、左、右三个方向以不同视角抓取局部点云，因此采用点云配准融合实现将不同视角的点云通过旋转和平移转换到同一个坐标系。本实施例采用基于静物长方体的猪体点云配准算法，以长宽高分别为60cm、50cm、50cm的标准长方体作为标定物，放置于点云获取通道，各相机坐标系与对应位置关系如图2所示。以上方相机坐标系作为世界坐标系，分别与两侧坐标系进行两两配准。采用PCL随机采样一致算法提取出各视角相机所捕获点云中的标定物平面点云，从中确定对应方向上标定物平面的四角关键点，通过关键点定位标定物平面在世界坐标系中的位置如图3所示。

将需要配准的两个拟合平面先通过平移变换到同一对公共点，接着分别进行两次旋转对齐，通过第一次旋转变换使得两个拟合平面存在公共边，在第二次旋转变换中使得两个拟合平面完全垂直，最终两个拟合平面完全还原实际物理位置关系。通过两次齐次旋转获取到相应的齐次变换矩阵R与T，根据公式(4)将坐标变换到同一坐标系当中。

其中，R为旋转变换矩阵，T为平移向量，[x₀ y₀ z₀]^T为待配准的原始坐标，[x₁ y₁z₁]^T为相应变换后的坐标。

由此本实施例可以获取以两侧均带有间隙开放孔洞的牲畜全局三维点云，如图4所示。实验中对养殖场20头活体猪进行数据采集，每一活体猪均采集4个完整三维点云，共得到80头活体猪点云数据。标准四肢类动物点云数据为普林斯顿Model Net、佐治亚理工学院的大型几何模型档案馆的标准恐龙、马、牛、狗以及豹子点云。

步骤S2：猪体呈现身体左右两侧对称分布，寻找中轴平面，将猪体全局三维点云左右侧对称划分为左侧点云和右侧点云，采用最小二乘法对中轴平面进行拟合，具体如下：

设全局三维点云数据为P＝{p₁，p₂，p₃，…，p_n}，其中点云p_i＝(x_i，y_i，z_i)，平面的一般方程：

z＝a₀x+a₁y+a₂ (2)

式中，a₀、a₁表示系数，a₂表示常数；

对于点云p_i(x_i,y_i,z_i)，要使得平面更贴近于点云，则使得S表达式取得最小值：

若满足S最小，则满足式(4)

有

通过求解式(5)方程组，得到a₀,a₁,a₂的值，从而得到平面方程z＝a₀x_i+a₁y+a₂，即为所拟合的中轴平面。如图5所示，灰色部分表示拟合中轴平面。

步骤S3：拟合中轴平面将全局三维点云划分为左侧和右侧点云，简称为左点云以及右点云。左右划分降低两边点云的互相干扰，更易于边界确定。左点云、右点云分别投影到XOY平面中，栏杆遮挡而空缺的点云其相邻数据点间坐标Y值会出现突变，由此来确定孔洞边界位置。点云散乱且分布无规律，单点X坐标值对应的点云数量少且无序，进而以块切片方式获取数量稳定有分布规律性的点云数据。因此将带有间隙开放孔洞的一侧点云投影到XOY平面中，得到投影后的点集，以块切片方式得到带有间隙开放孔洞的一侧点云上孔洞的上边界L_upper、下边界L_lower。具体如下：

以左点云为例，设左点云数据为L，点云个数为m，L＝{l₁,l₂,l₃,……,l_m}，l_i＝(x_i,y_i,z_i)，作XOY平面投影，投影后点集L_p＝{lp₁,lp₂,lp₃,…,lp_n}，lp_i＝(x_i,y_i,0)。

沿X轴方向切片，切片数量为B，切片大小slice。

切片j点集表示为S_j＝{s₁,s₂,s₃,……,s_n},s_k＝(x_k,y_k,0),x_k∈[min(x_i)+(j-1)*slice,min(x_i)+j*slice],j∈(1,2,3,…,B)，将S_j中数据点根据y值排序，使得

y_k+1>y_k，求出max(y_k+1-y_k)，当y_t+1-y_t＝max(y_k+1-y_k)时，[y_t+1,y_t]即为当前切片中栏杆区域，记录为[b_j1,b_j2]。根据式(7)计算得出左点云的栏杆上边界L_upper，下边界L_lower

同理可以求出右点云的栏杆上边界R_upper、下边界R_lower，采用公式(8)计算全局点云的栏杆的上边界U和下边界D，如图6所示，两条平行线所确定的栏杆区域(孔洞区域)的位置，不同颜色深浅代表不同切片区域。

步骤S4：以块切片方式得到腹部点云的下边界，并计算出腹部最低点Down，在孔洞的上边界U的基础上加h得到Top，从而将待修补区域确定为[Down,Top]，其中Down为下边界，Top为上边界。

对于腹部点云范围的确定，需将左点云与右点云沿X轴方向分为平均切分四片，分别选取第三片的最低点分别作为左下边界和右下边界。

具体的，以左点云为例，设需要拟合的腹部下边界为L_a，设置切片数量B₂＝4，根据式(6)计算切片大小slice，第三片切片所处的X轴区间为[min(x_i)+2*slice,min(x_i)+3*slice]，求出该切片范围内的点的最小y值为L_a＝min(y_i)。因此也求得右点云的腹部下边界为R_a。通过式(9)计算出腹部最低点Down：

在栏杆位置(孔洞位置)上边界U的基础上加h＝0.1得到Top，并将修补区域确定为[Down,Top]，其中，Down为下边界，而Top为上边界，如图7所示，中两平行线间区域代表所确定的修补范围。

步骤S5：将待修补区域点云沿X轴纵向切片，将切片投影到YOZ平面上，对投影点云进行有序化处理，根据待修补效果设置插值点密度，采用三次B样条曲线插值修补。

对于修补区域点云H＝{h₁,h₂,h₃,…,h_a},h_i＝(x_i,y_i,z_i)，切片的数量为Blocks，沿着X轴方向切片，根据式(10)得出切片大小为slice。

将切片投影到YOZ平面上，其中第切片j为Q_j＝{q₁,q₂,q₃,……,q_b},q_i＝(x_i,y_i,z_i)。如图8所示，表示完整切片的投影，如图9所示，表示待修补区域投影，此时点云的分布杂乱无章，需要对点云进行有序化处理，因待修补点云切片投影从左到右呈现先递减后递增的趋势，所以以中轴平面为分界，对切片投影进行左右划分，将划分后的左侧点云降序排序，右侧点云升序排序，使得在左侧点云中，

y_i+1＜y_i，而右侧点云中，

y_i+1>y_i，如图10所示，表示有序化后的待修补点云投影。

对于排序后的点云，采用三次B样条曲线的方法进行点的修补，其中B样条曲线的表达式为：

其中，P_i(i＝0,1,2,…,n)表示空间中n+1个控制点集；t＝{t₀,t₁,……,t_m}是m+1个矢量；F_i,k(t)为k阶或k-1次B样条曲线基函数。

B样条曲线的基函数为：

当k＝3时，则有三次B样条曲线的基函数：

将三次B样条曲线的基函数代入到B样条曲线的表达式中，得到三次B样条曲线方程：

P(t)＝P₀*F_0,3(t)+P₁*F_1,3(t)+P₂*F_2,3(t)+P₃*F_3,3(t) (14)

在三次B样条曲线中，每四个点能够确定一条三次曲线，图11为三次B样条曲线作用示意图，拟合其中P₁,P₂,P₃,P₄为四个控制点，M₁,M₂分别代表直线P₀P₂和P₁P₃的中点，所拟合曲线的起终点S,E分别位于P₁M₁与P₂M₂的三分之一处。

在传统的三次B样条曲线中，插值过程中采用统一的t值。而本实施例中所处理的猪体点云切片数据中，数据点与数据点之间的距离并不一致，使用统一间隔的t值对点云进行插值，结果会呈现出原始点云密集处的点云加倍密集，而稀疏处依旧稀疏，且如此会徒增修补过程中运算量而达不到想要的修补效果。因此本实施例所述的方法提出一种自适应的t值算法，具体的，设定期望点云之间最小间距d，对于所拟合的三次曲线g，通过四个控制点P₀＝f₀(0,y₀,z₀),P₁＝f₁(0,y₁,z₁),P₂＝f₂(0,y₂,z₂),P₃＝f₃(0,y₃,z₃)以及t＝0,t＝1计算出起初其起点g₀，终点g₁：

采用平均间隔的方式在拟合曲线采集数据点进行曲线弧长的估算，设置t＝0.1求得进行采点曲线上10个点，计算相邻点之间线段长度并求和得到Length，以此来替代所拟合曲线的弧长

计算插值密度。修补后点云数据集H的表达式如下：

式中，t_j表示根据密度Density计算出的值，P(t_j)表示根据三次B样条曲线的修补值；

通过式(16)计算插值密度Density：

式中，d表示期望的点云之间最小间距，Length表示曲线

的近似长度；

图12表示待修补区域经三次B样条插值修补后点云投影，图13表示修补完整切片投影。

步骤S6，修补后的点云数量的激增不便于后序测量工作且不利于数据存储，因此有必要进行点云下采样。本实施例采用体素滤波下采样算法对修补后的点云数据进行下采样处理，具体的，对于修补后的点云数据集H＝{h₁,h₂,h₃,…,h_m},h_i＝(x_i,y_i,z_i)，设定体素包围盒的长、宽、高分别为length,width,high，并分别对应XYZ三个维度。对于任意体素包围盒V_j，对任意h_i∈V_j,i≤m，包围盒V_j中有且仅有代表点v_j(x,y,z)。

为了验证本实施例所述的区域限定的改进三次B样条插值修补方法的技术效果，与现有技术进行对比分析，具体如下：

猪体点云孔洞修补，获取完整点云既有利于点云各部分分割，也有利于猪体尺测量，尤其是三围体尺测量。本实施例以修补后猪体点云数据计算腹围以验证算法的优劣，并与极坐标变换修补算法、曲面拟合修补算法以及三角化修补算法进行比较分析及腹围估算以检验本实施例所述的方法的稳健性。

本实施例首先对获取到的1:1猪体模型点云数据进行实验测量，图14表示采集到的带栏杆孔洞的三维猪体点云数据，共有数据点数量为73217个。

依据猪腹围、胸围、臀围的测量方法，确定其在三维点云数据中的测量位置，如图15所示，为三维数据中测量位置及其切片情况示意图。

分别采用极坐标变换修补、曲面拟合修补、三角化修补以及本实施例的方法分别进行修补实验。图16为四种修补算法修补结果展示，其中(a)(b)(c)(d)分别表示上述四种算法修补结果。

对修补后点云指定位置的切片计算腹围长度，其中各种修补算法计算腹围结果与人工测量的结果的对比情况如表1所示。表1结果显示本实施例的方法在对胸围、腹围以及臀围三个围度的估算上与人工测量值相比偏差最小。

表1猪体三维计算值

模型猪其姿态固定，无论是修补孔洞或是测量腹围都相对简单。而活体猪数据采集过程时易出现扭头、低头、抬头、弓背等各种姿态，使得猪躯体部分形态变化，给精确测量胸围、腹围以及臀围带来困难。

实验数据为对20头活猪，每一头猪在自由状态下通过采集通道获取不同姿态下的4组三维数据，共80组含孔洞猪体点云。每一头活猪均进行人工腹围测量。这80组点云数据分别采用极坐标变换修补、三角化修补、曲面拟合修补以及本文算法实现孔洞修补以及腹围估算。图17表示四种算法的腹围估算值与人工测量值之间误差变化情况，图18为极坐标变换修补算法腹围估计值与人工测量值的比较结果，图19为曲面拟合算法腹围估计值与人工测量值的比较结果，其中极坐标变换修补算法，图20为三角化修补算法腹围估计值与人工测量值的比较结果，其中极坐标变换修补算法，图21为本文算法腹围估计值与人工测量值的比较结果。其中极坐标变换修补算法、曲面拟合算法、三角化修补算法以及本文算法对腹围测算与人工测量结果绝对平均误差分别为3.32、2.93、5.03、1.91。相较于其它三种算法，本文算法修补的三维猪体点云所测算的腹围值最接近人工测量值，且在对同一头猪体的不同三维数据的腹围估算中表现更加平稳，说明相比于其它三种修补算法，本文算法具有更好的稳定性。

鲁棒性检验

本实施例采用三次B样条曲线拟合算法对多种标准动物点云进行孔洞修补。这些标准动物点云源自普林斯顿Model Net、佐治亚理工学院的大型几何模型档案馆的标准点云数据集，选取其中四只类动物数据dinosaur、horse、cow、dog、puma。通过随机函数为每一种动物点云制造栏杆孔洞以及腹部孔洞。图22表示实验所用点云数据的三维效果，其中(a)～(e)分别对应为dinosaur、horse、cow、dog、puma五种动物的标准点云展示效果，(f)～(j)分别为含有随机孔洞的dinosaur、horse、cow、dog以及puma点云数据。

分别用上述四种修补算法对含有随机孔洞的标准动物点云数据进行孔洞修补，图23、图24、图25、图26分别对应表示极坐标变换修补算法、曲面拟合修补算法、三角化修补算法以及本文修补算法修补效果展示。可以看出对于不同类型的四肢动物，各种修补算法的修补效果参差不齐，其中极坐标变换修补算法虽然对动物的腹部的修补效果良好，但容易受到动物的四肢以及头部的影响，导致修补效果大打折扣，曲面拟合算法能够较好的修补小范围的孔洞，但是对于大尺度的裂缝的修补效果不佳，三角化的修补算法在对侧面点云的平滑上表现良好，但是对于散乱的腹部点云的修补不尽人意。本文的算法无论是对于动物的侧面还是动物的腹部都有着较好的平滑和修补，极大地提升了对四肢动物点云的修补效果，具有较好的鲁棒性。

根据以上分析，本文提出区域限定的改进三次B样条近似拟合修补算法，比传统的三次B样条近似拟合修补得到的点云数据均匀性和逼近性能更好。通过对不同四肢动物的点云孔洞模型数据和养殖场活猪点云孔洞数据进行腹部修补，通过分析修补结果可见本文算法的稳健性和鲁棒性。

从修补视觉效果上，曲面拟合和本文算法较为接近原始点云。然对比腹围测量精度方面，本文算法均优于其他算法，具有较高的测量精度。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (9)

1.一种区域限定的改进三次B样条插值修补方法，其特征在于：所述的方法包括步骤如下：

S1：获取至少一侧带有间隙开放孔洞的牲畜全局三维点云；

S2：寻找中轴平面，将牲畜全局三维点云左右侧对称划分为左侧点云和右侧点云，采用最小二乘法对中轴平面进行拟合；

S3：将带有间隙开放孔洞的一侧点云投影到XOY平面中，得到投影后的点集，以块切片方式得到带有间隙开放孔洞的一侧点云上孔洞的上边界L_upper、下边界L_lower；计算出全局三维点云的孔洞的上边界U和下边界D；

S4：以块切片方式得到腹部点云的下边界，并计算出腹部最低点Down，在孔洞的上边界U的基础上加h得到Top，从而将待修补区域确定为[Down,Top]，其中Down为下边界，Top为上边界；

S5：将待修补区域点云沿X轴纵向切片，将切片投影到YOZ平面上，对投影点云进行有序化处理，设置插值点密度，采用三次B样条曲线插值方法完成修补；

步骤S3中，设带有间隙开放孔洞的一侧点云数据为L，点云个数为m，L＝{l₁,l₂,l₃,……,l_m},l_i＝(x_i,y_i,z_i)，作XOY平面投影，投影后点集为L_p＝{lp₁,lp₂,lp₃,…,lp_n},lp_i＝(x_i,y_i,0)；

沿X轴方向切片，切片数量为B，切片大小slice；

切片j点集表示为S_j＝{s₁,s₂,s₃,……,s_n},s_k＝(x_k,y_k,0),x_k∈[min(x_i)+(j-1)*slice,min(x_i)+j*slice],j∈(1,2,3,…,B)，k∈(1,2,3,…,n)；

将S_j中数据点根据y值进行排序，使得

求出max(y_k+1-y_k)，当y_t+1-y_t＝max(y_k+1-y_k)时，[y_t+1,y_t]即为当前切片中孔洞区域，记录为[b_j1,b_j2]；

从而计算得出一侧点云的孔洞上边界L_upper，下边界L_lower

2.根据权利要求1所述的区域限定的改进三次B样条插值修补方法，其特征在于：还包括步骤S6，为了便于后序测量工作且利于数据存储，采用体素滤波下采样算法对修补后的点云进行下采样处理。

3.根据权利要求2所述的区域限定的改进三次B样条插值修补方法，其特征在于：步骤S2，设全局三维点云数据为P＝{p₁,p₂,p₃,…,p_n}，其中点云p_i＝(x_i,y_i,z_i)，平面的一般方程：

z＝a₀x+a₁y+a₂ (2)

式中，a₀、a₁表示系数，a₂表示常数；

对于点云p_i(x_i,y_i,z_i)，则使得S表达式取得最小值：

若满足S最小，则满足式(4)

有

通过求解式(5)方程组，得到a₀,a₁,a₂的值，从而得到平面方程z＝a₀x_i+a₁y+a₂，即为所拟合的中轴平面。

4.根据权利要求3所述的区域限定的改进三次B样条插值修补方法，其特征在于：根据孔洞的上边界L_upper、下边界L_lower计算出全局三维点云的孔洞的上边界U和下边界D，具体如下：

5.根据权利要求4所述的区域限定的改进三次B样条插值修补方法，其特征在于：步骤S4，设一侧点云的腹部点云的下边界为L_a，设置切片数量B₂＝4，根据式(6)计算切片大小slice，选取第三片切片所处的X轴区间为[min(x_i)+2*slice,min(x_i)+3*slice]，求出第三片切片范围内的点的最小y值为L_a＝min(y_i)，根据式(9)计算出腹部的最低点Down如下：

式中，R_a表示另一侧点云的腹部点云的下边界。

6.根据权利要求5所述的区域限定的改进三次B样条插值修补方法，其特征在于：步骤S5，对于修补区域点云H＝{h₁,h₂,h₃,…,h_a},h_i＝(x_i,y_i,z_i)，切片的数量为Blocks，沿着X轴方向切片，根据式(10)得出切片大小为slice

式中，i∈(1,2,3,…,a)；

将切片投影到YOZ平面上，其中第j切片为Q_j＝{q₁,q₂,q₃,……,q_b},q_i＝(x_i,y_i,z_i)。

7.根据权利要求6所述的区域限定的改进三次B样条插值修补方法，其特征在于：对投影点云进行有序化处理，具体如下：因待修补点云切片投影从左到右呈现先递减后递增的趋势，所以以中轴平面为分界，对切片投影进行左右划分，将划分后的左侧点云降序排序，右侧点云升序排序，使得在左侧点云中，

而右侧点云中，

由此完成有序化处理。

8.根据权利要求7所述的区域限定的改进三次B样条插值修补方法，其特征在于：设定期望点云之间最小间距d，对于所拟合的三次曲线g，通过四个控制点P₀＝f₀(0,y₀,z₀),P₁＝f₁(0,y₁,z₁),P₂＝f₂(0,y₂,z₂),P₃＝f₃(0,y₃,z₃)以及t＝0,t＝1计算出起初其起点g₀，终点g₁：

采用平均间隔的方式在拟合曲线采集数据点进行曲线弧长的估算，设置t＝0.1求得进行采点曲线上10个点，计算相邻点之间线段长度并求和得到Length，以此来替代所拟合曲线的弧长

通过式(12)计算插值密度Density：

式中，d表示期望的点云之间最小间距，Length表示曲线

的近似长度；

三次B样条曲线的基函数：

将三次B样条曲线的基函数代入到B样条曲线的表达式中，得到三次B样条曲线方程

P(t)＝P₀*F_0,3(t)+P₁*F_1,3(t)+P₂*F_2,3(t)+P₃*F_3,3(t)#(14)

式中，t表示修补的使用三次B样条曲线修补的参数，P₀,P₁,P₂,P₃分别代表四个控制点；

修补后点云数据集H表示如下：

式中，t_j表示根据密度Density计算出的值，P(t_j)表示根据三次B样条曲线的修补值。

9.根据权利要求8所述的区域限定的改进三次B样条插值修补方法，其特征在于：步骤S6，采用体素滤波下采样算法对修补后的点云进行下采样处理，具体的，对于修补后的点云数据集H＝{h₁,h₂,h₃,…,h_m},h_i＝(x_i,y_i,z_i)，设定体素包围盒的长、宽、高分别为length,width,high，并分别对应XYZ三个维度；

对于任意体素包围盒V_j，对任意h_i∈V_j,i≤m，包围盒V_j中有且仅有代表点v_j(x,y,z)

式中，k＝1、2…m。

Images