一种基于粒子群算法的B样条曲线拟合方法及系统与流程

本发明涉及计算机辅助设计领域，具体涉及一种基于粒子群算法的B样条曲线拟合方法及系统。

背景技术：

在工业领域中，物理原型通常建在带有塑料和其他材料的车间中，这通常会造成大量浪费，并且很难再次修改。

逆向工程(又称逆向技术)，是一种能够解决上述技术问题的方式，其能够将基于产品实物，反向得到产品对应的数字模型，数字模型比物理模型或产品实物更便于地进行修改与分，可见一定程度上能够满足上述需求。

B样条曲线拟合，因其灵活性与可控性，已成为逆向工程中常用的一种数据点拟合方式。而且，为使得到的B样条曲线是最优的，科研领域对B样条曲线进行了深入研究，比如最近，研究人员将他们的注意力转向了元启发式技术，这些技术已被广泛应用于解决传统优化算法难以解决的困难优化问题。关于节点矢量的问题可以被视为多变量多模态非线性优化问题，并且可以通过优化算法来解决。Yoshimoto等人提出了一种方法，通过将节点位置的候选者作为基因并使用遗传算法解决问题，将原始的连续非线性优化问题转换为离散的组合优化问题。Sarfraz将这项工作从开曲线扩展到闭曲线，并把字体和符号作为案例进行研究。然而，遗传算法有两个主要缺点：缺乏局部搜索能力和早熟收敛。因此，这种近似方法很难找到全局最优值。为了避免遗传算法的缺陷，Ulker提出了一种不同的方法，把遗传算法替换成人工免疫算法，通过将节点位置的候选者作为抗体以解决该优化问题。上述这类方法都是基于将原始的连续优化问题转换为离散的组合优化问题。然而，这种转换过程会引入很大的离散误差，这使得它们对于真实数据集的拟合结果不准确且不可靠。

为此，本发明提供一种基于粒子群算法的B样条曲线拟合方法及系统，用于解决上述技术问题。

技术实现要素：

针对现有技术的上述不足，本发明提供一种基于粒子群算法的B样条曲线拟合方法及系统，用于提高B样条曲线拟合的准确性与可靠性。

第一方面，本发明提供一种基于粒子群算法的B样条曲线拟合方法，包括步骤：

S1、采集目标物体的轮廓图像；

S2、依据所采集的轮廓图像，获取目标物体的轮廓数据点；

S3、基于所获取的轮廓数据点，计算并得到B样条曲线的初始控制点和节点矢量；

S4、构造初始控制点的坐标的初始种群，并利用粒子群算法对初始控制点进行优化，得到优化后的控制点；

S5、依据所得到的节点矢量和优化后的控制点，构建拟合的B样条曲线。

进一步地，在步骤S2中，依据所采集的轮廓图像，获取目标物体的轮廓数据点，包括：

步骤S21、按照自上到下、自左到右的顺序，对所述的轮廓图像进行搜索，获取搜索到的第一个边界点，并将该边界点记为轮廓数据点d1；

步骤S22、依序在轮廓数据点d1的右、右下、下、左下四个方向，对所述的轮廓图像进行搜索，获取搜索到的第一个边界点，并记为轮廓数据点d2；

步骤S23、沿顺时针方向在轮廓数据点d2的右、右下、下、左下、左、左上、上、右上各方向上，对所述的轮廓图像进行搜索，获取搜索到的第一个边界点，并记为轮廓数据点d3；

步骤S24、判断当前获取到的最新的轮廓数据点是否与轮廓数据点d1重合：

是，则搜索过程结束，并输出所搜索到的各轮廓数据点；

否，则继续执行步骤S25；

步骤S25、沿顺时针方向，继续在当前获取到的最新的轮廓数据点的右、右下、下、左下、左、左上、上、右上各方向上，对所述的轮廓图像进行搜索，获取搜索到的第一个边界点，并记为轮廓数据点di，i为当前搜索到的第i个轮廓数据点；之后转而继续执行步骤S24。

进一步地，在步骤S3中，基于所获取的轮廓数据点，计算并得到B样条曲线的初始控制点和节点矢量，包括：

步骤S31、对步骤S2中获取的轮廓数据点进行采样，得到一组采样数据点；

步骤S32、依据步骤S31中得到的采样数据点，计算并得到B样条曲线的初始控制点和节点矢量。

进一步地，在步骤S31中，对步骤S2中获取的轮廓数据点进行采样，得到一组采样数据点，包括：

步骤S311、计算步骤S2中所获取的每个轮廓数据点的曲率；

步骤S312、依据步骤S311中计算出的各曲率，对步骤S2中所获取的轮廓数据点进行筛选，筛选出曲率的绝对值大于预先设定的曲率阈值的轮廓数据点，记为候选优势点；

步骤S313、计算各候选优势点的弯曲度，并确定各候选优势点对应的支撑区间；

步骤S314、计算各候选优势点的评价曲率，对各候选优势点进行下述①和②两方面的过滤，得到过滤后的候选优势点，记为优势点：

①删除评价曲率小于ε的候选优势点；

②对于各相邻的两个候选优势点：若评价曲率均大于ε，则删除评价曲率最小的候选优势点；若评价曲率均等于ε，则删除支撑区间最小的候选优势点；

其中ε为预先设定的评价曲率阈值；

步骤S315、分别计算各相邻优势点之间的间距，对于间距超过预先给定的采样间距阈值的相邻优势点，增设其中点为优势点，直至任意两个相邻优势点之间的间距均小于所述的采样间距阈值；

步骤S316、采集步骤S314中得到的各优势点以及步骤S315中增设的各优势点，即得到所有的采样数据点。

进一步地，在步骤S32中，依据步骤S31中得到的采样数据点，计算并得到B样条曲线的初始控制点和节点矢量，包括步骤：

步骤S321、采用积累弦长参数化方法，对步骤S31中得到的各采样数据点进行参数化，得到各采样数据点的参数值；

步骤S322、基于所得到的各采样数据点的参数值，采用代数平均值法得到所述的节点矢量；

步骤S323、基于所得到的节点矢量及各采样数据点的参数值，采用最小二乘法计算并得到B样条曲线的初始控制点。

进一步地，在步骤S4中，构造初始控制点的坐标的初始种群，并利用粒子群算法对初始控制点进行优化，得到优化后的控制点，包括：

步骤S41、记P＝{pi(xi,yi)|i＝0,1,…,n}为步骤S3中得到的所有的初始控制点形成的点集，构造初始控制点的坐标的初始种群POP＝{POP1,POP2,…,POPN},N为初始种群中个体的数量，每一个个体POPi均由在向量v＝(x0,x1,…,xn,y0,y1,…,yn)的每一项分别加上区间[-1,1]上的任意随机数值后得到，pi为步骤S3中得到的第i+1个初始控制点，(xi,yi)为初始控制点pi的坐标，向量v＝(x0,x1,…，xn,y0,y1,…,yn)由上述点集P中所有的初始控制点的坐标形成；

步骤S42、对于初始种群中的每一个个体POPi，设置其初始速度的每个分量为区间[-5,5]上的随机数；

步骤S43、初始化个体最优解：设置种群中每一个个体对应的个体最优解为初始种群中当前个体的当前值，计算每个个体的个体最优解的适应度；

步骤S44、初始化全局最优解：将步骤S43计算所得的个体的适应度最小的个体，设置为全局最优解gbest，并计算该全局最优解gbest的适应度；

步骤S45、计算种群中每个个体的适应度，分别将每个个体的适应度与对应个体的个体最优解的适应度进行比较，若当前个体的适应度小于当前个体的个体最优解的适应度，则更新当前个体的个体最优解、并更新该当前个体的个体最优解为当前个体的当前值，否则不更新该当前个体的个体最优解；之后继续执行步骤S45’；

步骤S45’、获取初始种群当前的种群最优解、以及计算该获取的种群最优解的适应度，并将当前获取的种群最优解的适应度与当前的全局最优解gbest的适应度进行比较，若当前的种群最优解的适应度小于当前的全局最优解gbest的适应度，则更新当前的全局最优解gbest、并更新该当前的全局最优解gbest为该当前的种群最优解，否则不更新当前的全局最优解gbest；

步骤S46、更新种群中每个个体的速度和位置；

步骤S47、判断是否达到预先设定的最大迭代次数：

是，则输出当前最新的全局最优解，并继续执行步骤S48；

否，则转而继续执行步骤S45；

步骤S48、将步骤S47输出的全局最优解调整为对应的二维坐标集合P′＝{p′i(x′i,y′i)|i＝0,1,…,n}，继而得到各优化后的控制点，其中，p′i为初始控制点pi对应的优化后的控制点，p′i(x′i,y′i)为优化后的控制点p′i的二维坐标。

另外，本发明还提供一种基于粒子群算法的B样条曲线拟合系统，包括：

轮廓图像采集单元，用于采集目标物体的轮廓图像；

数据点获取单元，用于依据轮廓图像采集单元采集的轮廓图像，获取目标物体的轮廓数据点；

初始控制点和节点矢量获取单元，用于基于数据点获取单元获取的轮廓数据点，计算并得到B样条曲线的初始控制点和节点矢量；

控制点优化单元，用于通过构造初始控制点的坐标的初始种群，利用粒子群算法对初始控制点和节点矢量获取单元获取到的初始控制点进行优化，得到优化后的控制点；

B样条曲线拟合单元，与所述的控制点优化单元以及初始控制点和节点矢量获取单元分别相连，用于依据所得到的节点矢量和优化后的控制点，构建拟合的B样条曲线。

进一步地，所述的数据点获取单元，包括：

第一数据点获取模块，与所述的轮廓图像采集单元相连，用于按照自上到下、自左到右的顺序，对轮廓图像采集单元采集的轮廓图像进行搜索，获取搜索到的第一个边界点，并将该边界点记为轮廓数据点d1；

第二数据点获取模块，与所述的第一数据点获取模块相连，用于依序在轮廓数据点d1的右、右下、下、左下四个方向，对所述的轮廓图像进行搜索，获取搜索到的第一个边界点，记为轮廓数据点d2；

第三数据点获取模块，与所述的第二数据点获取模块相连，用于沿顺时针方向在轮廓数据点d2的右、右下、下、左下、左、左上、上、右上各方向上，对所述的轮廓图像进行搜索，获取搜索到的第一个边界点，记为轮廓数据点d3；

第四数据点获取模块，用于在当前获取到的最新的轮廓数据点的右、右下、下、左下、左、左上、上、右上各方向上，对所述的轮廓图像进行搜索，获取搜索到的第一个边界点，记为轮廓数据点di，i为当前搜索到的第i个轮廓数据点；

判断模块，与所述的第三数据点获取模块、第四数据点获取模块分别相连，用于判断当前获取到的最新的轮廓数据点是否与轮廓数据点d1重合，并且：在判断结果为是时，用于控制数据点搜索过程结束，并控制输出所搜索到的各轮廓数据点；在判断结果为否时，用于调用第四数据点获取模块继续获取下一个轮廓数据点。

进一步地，初始控制点和节点矢量获取单元包括：

数据点采样模块，用于对数据点获取单元获取的轮廓数据点进行采样，得到一组采样数据点；

初始控制点和节点矢量计算模块，用于依据数据点采样模块得到的采样数据点，计算并得到B样条曲线的初始控制点和节点矢量；

其中，数据点采样模块对数据点获取单元获取的轮廓数据点进行采样，得到一组采样数据点的方法包括：

步骤S311、计算数据点获取单元获取的每个轮廓数据点的曲率；

步骤S312、依据步骤S311中计算出的各曲率，对数据点获取单元获取的轮廓数据点进行筛选，筛选出曲率的绝对值大于预先设定的曲率阈值的轮廓数据点，记为候选优势点；

步骤S313、计算各候选优势点的弯曲度，并确定各候选优势点对应的支撑区间；

步骤S314、计算各候选优势点的评价曲率，对各候选优势点进行下述①和②两方面的过滤，获取过滤后的候选优势点，并记为优势点：

①删除评价曲率小于ε的候选优势点；

②对于各相邻的候选优势点：若评价曲率均大于ε，则删除评价曲率最小的候选优势点；若评价曲率均等于ε，则删除支撑区间最小的候选优势点；

其中ε为预先设定的评价曲率阈值；

步骤S315、分别计算各相邻优势点之间的间距，对于间距超过预先给定的采样间距阈值的相邻优势点，增设其中点为优势点，直至任意两个相邻优势点之间的间距小于所述的采样间距阈值；

步骤S316、采集步骤S314中得到的各优势点以及步骤S315中增设的各优势点，即采样得到的所有的采样数据点；

所述的初始控制点和节点矢量计算模块，依据数据点采样模块得到的采样数据点，计算并得到B样条曲线的初始控制点和节点矢量，实现方法包括步骤：

步骤S321、采用积累弦长参数化方法，对数据点采样模块得到的各采样数据点进行参数化，得到各采样数据点的参数值；

步骤S322、基于所得到的各采样数据点的参数值，采用代数平均值法得到所述的节点矢量；

步骤S323、基于所得到的节点矢量及各采样数据点的参数值，采用最小二乘法计算并得到B样条曲线的初始控制点。

进一步地，所述的控制点优化单元，其构造初始控制点的坐标的初始种群，并利用粒子群算法对初始控制点进行优化，得到优化后的控制点，具体实现方法包括：

步骤S41、记P＝{pi(xi,yi)|i＝0,1,…,n}为初始控制点和节点矢量获取单元得到的所有的初始控制点形成的点集，构造初始控制点的坐标的初始种群POP＝{POP1,POP2,…,POPN},N为初始种群中个体的数量，每一个个体POPi均由在向量v＝(x0,x1,…,xn,y0,y1,…,yn)的每一项分别加上区间[-1,1]上的任意随机数值后得到，pi为初始控制点和节点矢量计算模块中得到的第i+1个初始控制点，(xi,yi)为初始控制点pi的坐标，向量v＝(x0,x1,…,xn,y0,y1,…,yn)由上述点集P中所有的初始控制点的坐标形成；

步骤S42、对于初始种群中的每一个个体POPi，设置其初始速度的每个分量为区间[-5,5]上的随机数；；

步骤S43、初始化个体最优解：设置种群中每一个个体对应的个体最优解为初始种群中当前个体的当前值，计算每个个体的个体最优解的适应度；

步骤S44、初始化全局最优解：将步骤S43计算所得的个体的适应度最小的个体，设置为全局最优解gbest，并计算该全局最优解gbest的适应度；

步骤S45、计算种群中每个个体的适应度，分别将每个个体的适应度与对应个体的个体最优解的适应度进行比较，若当前个体的适应度小于当前个体的个体最优解的适应度，则更新当前个体的个体最优解、并更新该当前个体的个体最优解为当前个体的当前值，否则不更新该当前个体的个体最优解；之后继续执行步骤S45’；

步骤S45’、获取初始种群当前的种群最优解、以及计算该获取的种群最优解的适应度，并将当前获取的种群最优解的适应度与当前的全局最优解gbest的适应度进行比较，若当前的种群最优解的适应度小于当前的全局最优解gbest的适应度，则更新当前的全局最优解gbest、并更新该当前的全局最优解gbest为该当前的种群最优解，否则不更新当前的全局最优解gbest；

步骤S46、更新种群中每个个体的速度和位置；

步骤S47、判断是否达到预先设定的最大迭代次数：

是，则输出当前最新的全局最优解，并继续执行步骤S48；

否，则转而继续执行步骤S45；

步骤S48、将步骤S47输出的全局最优解调整为对应的二维坐标集合P′＝{p′i(x′i,y′i)|i＝0,1,…,n}，继而得到各优化后的控制点，其中，p′i为初始控制点pi对应的优化后的控制点，p′i(x′i,y′i)为优化后的控制点p′i的二维坐标。

本发明的有益效果在于：

本发明提供的基于粒子群算法的B样条曲线拟合方法及系统，均通过构造初始控制点的坐标的初始种群，并利用粒子群算法对初始控制点进行优化，从而得到优化后的控制点，继而结合节点矢量拟合出目标物体轮廓的B样条曲线，可见能够实现控制点位置的优化调整，避免了引入离散误差，使本发明能适应于真实数据集拟合，提高了B样条曲线拟合的准确度和可靠性。

此外，本发明设计原理可靠，结构简单，具有非常广泛的应用前景。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，对于本领域普通技术人员而言，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明一个实施例的方法的示意性流程图。

图2是本发明一个实施例的系统的示意性框图。

具体实施方式

为了使本技术领域的人员更好地理解本发明中的技术方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。

实施例1：

图1是本发明一个实施例的方法的示意性流程图。其中，图1执行主体可以为一种装置或系统。

如图1所示，该方法100包括：

步骤110，采集目标物体的轮廓图像；

步骤120，依据所采集的轮廓图像，获取目标物体的轮廓数据点；

步骤130，基于所获取的轮廓数据点，计算并得到B样条曲线的初始控制点和节点矢量；

步骤140，构造初始控制点的坐标的初始种群，并利用粒子群算法对初始控制点进行优化，得到优化后的控制点；

步骤150，依据所得到的节点矢量和优化后的控制点，构建拟合的B样条曲线。

可选地，作为本发明一个实施例，在步骤120中，依据所采集的轮廓图像，获取目标物体的轮廓数据点，包括：

步骤S21、按照自上到下、自左到右的顺序，对所采集的轮廓图像进行搜索，获取搜索到的第一个边界点，并将该边界点记为轮廓数据点d1；之后执行步骤S22；

步骤S22、依序在轮廓数据点d1的右、右下、下、左下四个方向，对所述的轮廓图像进行搜索，获取搜索到的第一个边界点，并记为轮廓数据点d2；之后执行步骤S23；

步骤S23、沿顺时针方向在轮廓数据点d2的右、右下、下、左下、左、左上、上、右上各方向上，对所述的轮廓图像进行搜索，获取搜索到的第一个边界点，并记为轮廓数据点d3；之后执行步骤S24；

步骤S24、判断当前获取到的最新的轮廓数据点是否与轮廓数据点d1重合：

是，则搜索过程结束，并输出所搜索到的各轮廓数据点；

否，则继续执行步骤S25；

步骤S25、沿顺时针方向，继续在当前获取到的最新的轮廓数据点的右、右下、下、左下、左、左上、上、右上各方向上，对所述的轮廓图像进行搜索，获取搜索到的第一个边界点，并记为轮廓数据点di，i为当前搜索到的第i个轮廓数据点；之后转而执行步骤S24。

综上，可顺序获取到目标物体的轮廓数据点d1、d2、…、di、…。

可选地，作为本发明一个实施例，在步骤S3中，基于所获取的轮廓数据点，计算并得到B样条曲线的初始控制点和节点矢量，包括：

步骤S31、对步骤S2中获取的轮廓数据点进行采样，得到一组采样数据点；

步骤S32、依据步骤S31中得到的采样数据点，计算并得到B样条曲线的初始控制点和节点矢量。

其中，在步骤S31中，对步骤S2中获取的轮廓数据点进行采样，得到一组采样数据点，包括：

步骤S311、计算步骤S2中所获取的每个轮廓数据点的曲率K；

步骤S312、依据步骤S311中计算出的各曲率K，对步骤S2中所获取的轮廓数据点依次进行筛选，筛选出曲率K的绝对值大于预先设定的曲率阈值的轮廓数据点，记为候选优势点；

步骤S313、计算各候选优势点的弯曲度，并确定各候选优势点对应的支撑区间；

步骤S314、计算各候选优势点的评价曲率，对各候选优势点进行下述①和②两方面的过滤，得到过滤后的候选优势点,记为优势点：

①删除评价曲率小于ε的候选优势点；

②对于各相邻的候选优势点：若评价曲率均大于ε，则删除评价曲率最小的候选优势点；若评价曲率均等于ε，则删除支撑区间最小的候选优势点；

其中ε为预先设定的评价曲率阈值；

步骤S315、分别计算各相邻优势点之间的间距，对于间距超过预先给定的采样间距阈值的相邻优势点，增设其中点为优势点，直至任意两个相邻优势点之间的间距小于所述的采样间距阈值；

步骤S314中得到的各优势点以及步骤S315中增设的各优势点，即为采样得到的所有的采样数据点。

需要说明的是，步骤S312对步骤S2中顺序获取的轮廓数据点依序进行筛选，所选出的候选优势点也是有顺序的；步骤S314对步骤S312中选出的候选优势点进行筛选，若删除步骤S2中获取的轮廓数据点中的任意一个候选优势点，则删除的该候选优势点则成为一个间隔点，删除的该候选优势点的上一个候选优势点和下一个候选优势点则不属于相邻的候选优势点。比如步骤S312中选出的候选优势点有d1、d2、d3、d4、…、di、…，候选优势点d3因其评价曲率小于ε被删除，则候选优势点d2和d4不属于步骤S314中所述的相邻的候选优势点。

另外需要说明的是，以步骤S314中得到的优势点有d1、d2、d4、…为例：若相邻优势点d1和d2之间的间距大于所述的采样间距阈值，且点d1-2为相邻优势点d1和d2之间的中点，则基于步骤S315，点d1-2也被增设为优势点，成为一个采样数据点，此时优势点d1-2与d1成为新增加的相邻优势点，优势点d1-2与d2也成为新增加的相邻优势点，若新增加的相邻优势点之间的间距大于所述的采样间距阈值，则继续增设其中点为优势点，直至各相邻优势点之间的间隔小于所述的采样间距阈值。

另外需要说明的是，本实施例中所涉及的各“相邻优势点之间的间距”，指的是相邻优势点对应候选优势点之间的候选优势点的数量；相对应地，本实施例中所涉及的“采样间距阈值”，表示预先设定的相邻优势点对应候选优势点之间的候选优势点的数量阈值。另外具体实现时，本领域技术技术人员还可依据实际情况，采用相应的距离长度替换上述候选优势点的数量。

其中，在步骤S32中，依据步骤S31中得到的采样数据点，计算并得到B样条曲线的初始控制点和节点矢量，包括步骤：

步骤S321、采用积累弦长参数化方法，对步骤S31中得到的各采样数据点进行参数化，得到各采样数据点的参数值，分别记为其中i表示第i+1个采样数据点，表示第i+1个采样数据点的参数值，步骤S31中得到的采样数据点的总数为m+1个；

步骤S322、基于所得到的各采样数据点的参数值，采用代数平均值法得到节点矢量，记U＝[u0,u1,…,ui,…,un+k+1],ui表示第i+1个采样数据点对应的节点值；

步骤S323、基于所得到的节点矢量及各采样数据点的参数值，采用最小二乘法计算并得到B样条曲线的初始控制点。

可选地，作为本发明一个实施例，在步骤S4中，构造初始控制点的坐标的初始种群，并利用粒子群算法对初始控制点进行优化，得到优化后的控制点，包括：

步骤S41、记P＝{pi(xi,yi)|i＝0,1,…,n}为步骤S3中得到的所有的初始控制点形成的点集，构造初始控制点的坐标的初始种群POP＝{POP1,POP2,…,POPN},N为初始种群中个体的数量，每一个个体POPi均由在向量v＝(x0,x1,…,xn,y0,y1,…,yn)的每一项分别加上区间[-1,1]上的任意随机数值后得到，pi为步骤S3中得到的第i+1个初始控制点，(xi,yi)为初始控制点pi的坐标，向量v＝(x0,x1,…,xn,y0,y1,…,yn)由上述点集P中所有的初始控制点的坐标形成；

步骤S42、对于初始种群中的每一个个体POPi，设置其初始速度的每个分量为区间[-5,5]上的随机数；

步骤S43、初始化个体最优解：设置种群中每一个个体对应的个体最优解为初始种群中当前个体的当前值，计算每个个体的个体最优解的适应度；

步骤S44、初始化全局最优解：将步骤S43计算所得的个体的适应度最小的个体，设置为全局最优解gbest，并计算该全局最优解gbest的适应度；

步骤S45、计算种群中每个个体的适应度，分别将每个个体的适应度与对应个体的个体最优解的适应度进行比较，若当前个体的适应度小于当前个体的个体最优解的适应度，则更新当前个体的个体最优解、并更新该当前个体的个体最优解为当前个体的当前值，否则不更新该当前个体的个体最优解；之后继续执行步骤S45’；

步骤S45’、获取初始种群当前的种群最优解、以及计算该获取的种群最优解的适应度，并将当前获取的种群最优解的适应度与当前的全局最优解gbest的适应度进行比较，若当前的种群最优解的适应度小于当前的全局最优解gbest的适应度，则更新当前的全局最优解gbest、并更新该当前的全局最优解gbest为该当前的种群最优解，否则不更新当前的全局最优解gbest；

步骤S46、更新种群中每个个体的速度和位置；

步骤S47、判断是否达到预先设定的最大迭代次数：

是，则输出当前最新的全局最优解，并继续执行步骤S48；

否，则转而继续执行步骤S45；

步骤S48、将步骤S47输出的全局最优解调整为对应的二维坐标集合P′＝{p′i(x′i,y′i)|i＝0,1,…,n}，继而得到各优化后的控制点，其中，p′i为初始控制点pi对应的优化后的控制点，p′i(x′i,y′i)为优化后的控制点p′i的二维坐标。

为了便于对本发明的理解，下面以本发明的基于粒子群算法的B样条曲线拟合方法的原理，结合实施例中对目标物体轮廓的B样条曲线进行拟合的过程，对本发明提供的基于粒子群算法的B样条曲线拟合方法做进一步的描述。

具体的，以目标物体A为例，基于本发明所述的基于粒子群算法的B样条曲线拟合方法，拟合目标物体A的轮廓的k阶B样条曲线，具体拟合步骤包括：

步骤101，采集目标物体A的轮廓图像。

步骤102，依据所采集的轮廓图像，获取目标物体A的轮廓数据点，形成点集D＝{di(xi,yi)|i＝1,2,…,r}，其中di为所获取的第i个轮廓数据点，di(xi,yi)为轮廓数据点di的二维坐标，r为依据所采集的轮廓图像所获取的目标物体A的轮廓数据点的总数。

步骤103，基于步骤102中所获取的轮廓数据点，计算并得到B样条曲线的初始控制点和节点矢量。

首先，对步骤102中获取的轮廓数据点进行采样：

(1)按如下式①计算D＝{di(xi,yi)|i＝1,2,…,r}中每个轮廓数据点di(xi,yi)的曲率：

式①中，和可定义为：

(2)基于(1)中计算出的各曲率K，对步骤102中获取的轮廓数据点进行筛选，并筛选出曲率K的绝对值大于预先设定的曲率阈值a的轮廓数据点di(xi,yi)，记为候选优势点。

其中，在该步骤(2)中，曲率K小于或等于预先设定的曲率阈值a的轮廓数据点，则被删除。另外，所述曲率阈值a的取值区间范围一般为[2,5]。

(3)从上述(2)中得到的候选优势点中筛选优势点：

a)确定支撑区间。

首先依据经验确定支撑区间的最大值为L,然后计算每个候选优势点的弯曲度bij,计算公式如下：

bij＝max(|(xi-j-xi)+(xi+j-xi)|,|(yi-j-yi)+(yi+j-yi)|)；

接着计算每个候选优势点对应的支撑区间vi：随着使j从1逐步增大到L，在此过程中：如果bij＜bi(j+1),则j就是第i个候选优势点的支撑区间vi；否则如果bij＝bi(j+1)，那么最大的满足这个条件的j就是第i个候选优势点的支撑区间vi；否则依据经验确定的最大的支撑区间L就是第i个候选优势点的支撑区间vi；

其中，i表示筛选出的第i个候选优势点,bij表示第i个候选优势点的弯曲度，j是第i个候选优势点的支撑区间的大小，j＝1,2,…,L。

b)确定了支撑区间后,计算各候选优势点的评价曲率：

其中，式③中i表示第i个候选优势点，bvi表示第i个候选优势点的评价曲率，ki表示上述a)中确定的第i个候选优势点的支撑区间，bij表示第i个候选优势点的弯曲度。

c)从各候选优势点中筛选出优势点：

对各候选优势点进行下述①和②两方面的过滤，过滤后得到候选优势点均为所述的优势点：

①删除评价曲率＜ε的候选优势点；

②对于各相邻的候选优势点：若评价曲率均大于ε，则删除评价曲率最小的候选优势点；若评价曲率均等于ε，则删除支撑区间最小的候选优势点。

其中，ε为预先设定的评价曲率阈值，取值范围一般为[2,5]。

(4)新增优势点：

通过上述步骤(3)得到一组优势点，分别计算该组优势点中各相邻的优势点之间的间距，对于间距超过预先给定的采样间距阈值t的相邻优势点，对应增设相邻优势点之间的中点为优势点，直至任意两个相邻优势点之间的间距小于所述的采样间距阈值t。

其中，采样间距阈值t的取值范围一般为[5,10]。

(5)采样数据点：采集上述(3)中得到的各优势点以及上述(4)中增设的各优势点，即得到所有的采样数据点。记集合Q＝{qi(xi,yi)|i＝0,1,…,m}为采样所得到的所有的采样数据点形成的集合。

接着，依据上述采样得到的采样数据点Q＝{qi(xi,yi)|i＝0,1,…,m}，获取B样条曲线的初始控制点和节点矢量：

(1)基于采样所得到的所有的采样数据点Q＝{qi(xi,yi)|i＝0,1,…,m}，采用规范积累弦长参数化法，计算采样数据点的参数值计算公式如下所示：

式中，Δqi＝qi+1-qi，为采样数据点之间的弦线矢量。

(2)节点配置以确定节点矢量U＝[u0,u1,…,uk,uk+1,…,un,un+1,…,un+k+1]：

设定初始控制点的数量为n+1，对应的坐标参数形成集合P＝{pi(xi,yi)|i＝0,1,…,n，其中m>n，该k为逼近曲线的次数且该k≥1。

根据端点插值与曲线定义域要求，采用定义域两端节点为k+1重的重节点端点条件，也即固支条件，定义：u0＝u1＝…＝uk＝0,un+1＝un+1＝…＝un+k+1＝1。定义域共包含n-k+1个节点区间，其节点值的选取应反映数据点参数值的分布情况，可按以下决定：

设c是一个正实数，令

则定义域的内节点为：

按上述方法,即可得到上述节点矢量：

U＝[u0,u1,…,uk,uk+1,…,un,un+1,…,un+k+1]。

(3)基于上述得到的采样数据点的参数值和节点矢量U＝[u0,u1,…,uk,uk+1,…,un,un+1,…,un+k+1]，采用最小二乘法计算并得到B样条曲线的初始控制点P＝{pi(xi,yi)|i＝0,1,…,n}：

B样条曲线逼近，是试图寻找一条k次B样条曲线：

其中，pj是控制点，Nj,k(u)是基函数。

基于上述计算得到的节点矢量U＝[u0,u1,…，uk，uk+1，…，un，un+1，…，un+k+1]，定义上述基函数Nj，k(u)如下：

并且，由如上所述的固支条件，有p0＝q0，pn＝qm。

根据最小二乘算法，可以得到以n-1个控制点p1，p2，…，pn-1为未知量的n-1个方程的方程组：

(N^TN)D＝R

其中，N是(m-1)×(n-1)阶标量矩阵，

N^T是N的转置矩阵。R和D都是含n-1个矢量元素的列阵，

其中，

用高斯消元法，可解出D，从而得到所有的控制点P＝{pi(xi，yi)|i＝0，1，…，n}。

步骤104，构造初始控制点的坐标的初始种群，并利用粒子群算法对初始控制点进行优化，得到优化后的控制点:

①节点矢量U＝[u0，u1，…，uk，uk+1，…，un，un+1，…，un+k+1]保持不变，将所得到的n+1个初始控制点P＝{pi(xi，yi)|i＝0，1，…，n}展开成一维向量，即可得到由2(n+1)个坐标值组成的向量v＝(x0，x1，…，xn，y0，y1，…，yn)，采用在一维向量v的每一项分别加上[-1,1]之间的一个随机值的方法，构造初始控制点的坐标的初始种群POP＝{POP1，POP2，…，POPN},N为初始种群中个体的数量，每一个个体POPi均由在向量v＝(x0，x1，…，xn，y0，y1，…，yn)的每一项分别加上区间[-1,1]上的任意随机数值后得到；

②对于每一个粒子POPi，设置其初始速度的每个分量为[-5,5]上的随机值，如：第i个粒子POPi＝(xi1,xi2,…，xiD)，设置其初始速度为vi＝(vi1,vi2,…,viD)，其中vij∈[-5,5]，此处j＝1,2,…,D，D＝2(n+1)；

③初始化个体最优解：设置种群中每个个体对应的个体最优解为初始种群中当前个体的当前值，分别计算初始种群POP＝{POP1,POP2,…,POPN}中每个个体POPi的适应度，根据曲线的距离平方和误差,其适应度函数定义如下：

其中，qj是第j个采样数据点的坐标，是采样数据点qj参数化的值；

④初始化全局最优解：找到③中适应度值最小的个体POPt＝(xt1,xt2,…,xtD)，设置为全局最优解gbest＝POPt＝(xt1,xt2,…,xtD)，并计算该全局最优解gbest的适应度；

⑤设定最大迭代次数为100，在每一次迭代中：

分别计算种群中各个粒子(即初始种群POP的每个个体POPi)的适应度，表示第k次迭代粒子POPi的位置矢量，将每个粒子的适应度与其个体最优解pbesti的适应度进行比较，若当前粒子的适应度小于其个体最优解pbesti的适应度，则更新其个体最优解为当前的粒子，即否则不更新；

之后获取初始种群当前的种群最优解、以及计算该获取的种群最优解的适应度，并将当前获取的种群最优解的适应度与当前的全局最优解gbest的适应度进行比较，若当前的种群最优解的适应度小于当前的全局最优解gbest的适应度，则更新当前的全局最优解gbest、并更新该当前的全局最优解gbest为该当前的种群最优解，否则不更新当前的全局最优解gbest。

对于每一代粒子，它的第d维(1≤d≤D)速度和位置根据如下公式进行更新：

其中上述公式(3)和(4)中，W称为惯性权重，c1和c2是加速度常数，一般取c1＝c2∈[0,4]。r1和r2是两个随机函数，取值范围为[0,1]，以增加搜索随机性。表示第k次迭代粒子POPi飞行速度矢量的第d维分量，表示第k次迭代粒子POPi位置矢量的第d维分量。pbestid表示粒子POPi个体最优解的第d维分量，gbestd表示群体的全局最优解的第d维分量。

根据上述方法，经过多次迭代直至达到最大迭代次数，从而得到优化后控制点的最优值gbest＝xt＝(xt1,xt2,…,xtD)＝(x′0,x′1,…,x′n,y′0,y′1,…,y′n)，将其调整成二维坐标P′＝{p′i(x′i,y′i)|i＝0,1，…，n}。

步骤105，依据所得到的节点矢量和优化后的控制点，构建拟合的B样条曲线:

基于以上得到的n+k+2个节点矢量和n+1个优化后的控制点,由他们定义的k次B样条曲线，如下：

其中，pj是控制点，Nj，k(u)是定义在节点矢量U上的k次基函数，该基函数可由以下递推关系定义，即

将[0,1]区间均匀划分1000等份，得到T＝[t0，t1，…，t1000]，t0＝0，t1000＝1，将每个ti作为变量u输入(1)式，从而构造出一条优化后的B样条曲线。

综上可见，本发明基于目标物体的轮廓图像，反算出B样条曲线的节点矢量和初始控制点，再把初始控制点位置候选者作为种群的个体并使用粒子群算法进行优化，从而得到优化后的B样条曲线的控制点，之后重建出一个B样条曲线。可见本发明采用优化调整控制点的位置的方式，并没有像遗传算法那样将原问题转换成离散组合优化问题，避免了引入离散误差，使得本方法能适应于真实数据集拟合，一定程度上提高了曲线拟合的准确度和可靠性。

实施例2：

图2为本发明所述基于粒子群算法的B样条曲线拟合系统的一个实施例。该实施例中的基于粒子群算法的B样条曲线拟合系统与实施例1中的基于粒子群算法的B样条曲线拟合方法相对应。

参见图2，本发明提供了一种基于粒子群算法的B样条曲线拟合系统200，该系统200包括：

轮廓图像采集单元201，用于采集目标物体的轮廓图像；

数据点获取单元202，用于依据轮廓图像采集单元201采集的轮廓图像，获取目标物体的轮廓数据点；

初始控制点和节点矢量获取单元203，用于基于数据点获取单元202获取的轮廓数据点，计算并得到B样条曲线的初始控制点和节点矢量；

控制点优化单元204，用于通过构造初始控制点的坐标的初始种群，利用粒子群算法对初始控制点和节点矢量获取单元203获取到的初始控制点进行优化，得到优化后的控制点；

B样条曲线拟合单元205，与所述的控制点优化单元204以及初始控制点和节点矢量获取单元203分别相连，用于依据所得到的节点矢量和优化后的控制点，构建拟合的B样条曲线。

使用时，先通过轮廓图像采集单元201采集目标物体的轮廓图像，之后通过数据点获取单元202在轮廓图像采集单元201采集的轮廓图像上获取目标物体的轮廓数据点；之后由初始控制点和节点矢量获取单元203，基于数据点获取单元获取的轮廓数据点，计算并得到B样条曲线的初始控制点和节点矢量；之后由控制点优化单元204，通过构造初始控制点的坐标的初始种群，利用粒子群算法对初始控制点和节点矢量获取单元203获取到的初始控制点进行优化，得到优化后的控制点；之后通过B样条曲线拟合单元205，依据所得到的节点矢量和优化后的控制点，构建拟合的B样条曲线。

其中，所述的数据点获取单元202，包括：

第一数据点获取模块，与所述的轮廓图像采集单元201相连，用于按照自上到下、自左到右的顺序，对轮廓图像采集单元201采集的轮廓图像进行搜索，获取搜索到的第一个边界点，并将该边界点记为轮廓数据点d1；

第二数据点获取模块，与所述的第一数据点获取模块相连，用于依序在轮廓数据点d1的右、右下、下、左下四个方向，对上述轮廓图像进行搜索，获取搜索到的第一个边界点，记为轮廓数据点d2；

第三数据点获取模块，与所述的第二数据点获取模块相连，用于沿顺时针方向在轮廓数据点d2的右、右下、下、左下、左、左上、上、右上各方向上，对上述轮廓图像进行搜索，获取搜索到的第一个边界点，记为轮廓数据点d3；

第四数据点获取模块，用于在当前获取到的最新的轮廓数据点的右、右下、下、左下、左、左上、上、右上各方向上，对上述轮廓图像进行搜索，获取搜索到的第一个边界点，记为轮廓数据点di，i为当前搜索到的第i个轮廓数据点；

判断模块，与所述的第三数据点获取模块、第四数据点获取模块分别相连，用于判断当前获取到的最新的轮廓数据点是否与轮廓数据点d1重合，并且：在判断结果为是时，用于控制数据点搜索过程结束，并控制输出所搜索到的各轮廓数据点；在判断结果为否时，用于调用第四数据点获取模块继续获取下一个轮廓数据点。

可优选地，初始控制点和节点矢量获取单元203包括：

数据点采样模块，用于对数据点获取单元202获取的轮廓数据点进行采样，得到一组采样数据点；

初始控制点和节点矢量计算模块，用于依据数据点采样模块得到的采样数据点，计算并得到B样条曲线的初始控制点和节点矢量；

其中，数据点采样模块对数据点获取单元202获取的轮廓数据点进行采样，得到一组采样数据点，包括：

步骤S311、计算数据点获取单元202获取的每个轮廓数据点的曲率；

步骤S312、依据步骤S311中计算出的各曲率，对数据点获取单元202所获取的轮廓数据点进行筛选，筛选出曲率的绝对值大于预先设定的曲率阈值的轮廓数据点，记为候选优势点；

步骤S313、计算各候选优势点的弯曲度，并确定各候选优势点对应的支撑区间；

步骤S314、计算各候选优势点的评价曲率，对各候选优势点进行下述①和②两方面的过滤，得到过滤后的候选优势点，记为优势点：

①删除评价曲率＜ε的候选优势点；

②对于各相邻的候选优势点：若评价曲率均大于ε，则删除评价曲率最小的候选优势点；若评价曲率均等于ε，则删除支撑区间最小的候选优势点；

其中ε为预先设定的评价曲率阈值；

步骤S315、分别计算各相邻优势点之间的间距，对于间距超过预先给定的采样间距阈值的相邻优势点，增设其中点为优势点，直至任意两个相邻优势点之间的间距小于所述的采样间距阈值；

步骤S316、采集步骤S314中得到的各优势点以及步骤S315中增设的各优势点，即得到所有的采样数据点。

所述的初始控制点和节点矢量计算模块，依据数据点采样模块得到的采样数据点，计算并得到B样条曲线的初始控制点和节点矢量，实现方法包括步骤：

步骤S321、采用积累弦长参数化方法，对数据点采样模块得到的各采样数据点进行参数化，得到各采样数据点的参数值，分别记为其中i表示第i+1个采样数据点，表示第i+1个采样数据点的参数值，数据点采样模块得到的采样数据点的总数为m+1个；

步骤S322、基于所得到的各采样数据点的参数值，采用代数平均值法得到节点矢量，记U＝[u0，u1，…,ui,…,un+k+1],ui表示第i+1个采样数据点对应的节点值；

步骤S323、基于所得到的节点矢量及各采样数据点的参数值，采用最小二乘法计算并得到B样条曲线的初始控制点。

其中，所述的控制点优化单元204，其构造初始控制点的坐标的初始种群，并利用粒子群算法对初始控制点进行优化，得到优化后的控制点，具体实现方法包括：

步骤S41、记P＝{pi(xi,yi)|i＝0,1,…,n}为初始控制点和节点矢量获取单元203得到的所有的初始控制点形成的点集，构造初始控制点的坐标的初始种群POP＝{POP1,POP2，…，POPN},N为初始种群中个体的数量，每一个个体POPi均由在向量v＝(x0，x1,…,xn,y0,y1,…,yn)的每一项分别加上区间[-1,1]上的任意随机数值后得到，pi为初始控制点和节点矢量计算模块中得到的第i+1个初始控制点，(xi，yi)为初始控制点pi的坐标，向量v＝(x0，x1，…，xn，y0，y1，…，yn)由上述点集P中所有的初始控制点的坐标形成；

步骤S42、对于每一个个体POPi，设置其初始速度的每个分量为区间[-5,5]上的随机数；

步骤S43、初始化个体最优解：设置种群中每一个个体对应的个体最优解为初始种群中当前个体的当前值，计算每个个体的个体最优解的适应度；

步骤S44、初始化全局最优解：将步骤S43中计算所得的适应度最小的个体，设置为全局最优解gbest，并计算该全局最优解gbest的适应度；

步骤S45、计算种群中每个个体的适应度，分别将每个个体的适应度与对应个体的个体最优解的适应度进行比较，若当前个体的适应度小于当前个体的个体最优解的适应度，则更新当前个体的个体最优解、并更新该当前个体的个体最优解为当前个体的当前值，否则不更新该当前个体的个体最优解；之后继续执行步骤S45’；

步骤S45’、获取初始种群当前的种群最优解、以及计算该获取的种群最优解的适应度，并将当前获取的种群最优解的适应度与当前的全局最优解gbest的适应度进行比较，若当前的种群最优解的适应度小于当前的全局最优解gbest的适应度，则更新当前的全局最优解gbest、并更新该当前的全局最优解gbest为该当前的种群最优解，否则不更新当前的全局最优解gbest；

步骤S46、更新种群中每个个体的速度和位置；

步骤S47、判断是否达到预先设定的最大迭代次数：

是，则输出当前最新的全局最优解，并继续执行步骤S48；

否，则转而继续执行步骤S45；

步骤S48、将步骤S47输出的全局最优解调整为对应的二维坐标集合P′＝{p′i(x′i,y′i)|i＝0,1,…,n}，继而得到各优化后的控制点，其中，p′i为初始控制点pi对应的优化后的控制点，p′i(x′i,y′i)为优化后的控制点p′i的二维坐标。

因此，本实施例中所述的基于粒子群算法的B样条曲线拟合系统，所能达到的技术效果可以参见上文中的描述，此处不再赘述。

需要说明的是，本说明书中各个实施例之间相同相似的部分互相参见即可。尤其，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例中的说明即可。

需要说明的是，本发明中所涉及的“右下”方向为“右”方向顺时针旋转45°后得到的方向，“右上”方向为“右”方向逆时针旋转45°后得到的方向，“左下”方向为“左”方向逆时针旋转45°后得到的方向，“左上”方向为“左”方向顺时针旋转45°后得到的方向。

另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独存在，也可以两个或两个以上单元集成在一个单元中。

综上可见，本发明采用粒子群算法优化调整控制点的位置，避免了先对节点矢量进行优化再计算出控制点的方法的使用，一定程度上简化了计算过程，简单易行且收敛速度快，可见一定程度上提高了曲线拟合的准确度、可靠性及效率。

尽管通过参考附图并结合优选实施例的方式对本发明进行了详细描述，但本发明并不限于此。在不脱离本发明的精神和实质的前提下，本领域普通技术人员可以对本发明的实施例进行各种等效的修改或替换，而这些修改或替换都应在本发明的涵盖范围内/任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。