MathJax.Hub.Config({tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]}}); function MyAutoRun() {    var topp=$(window).height()/2; if($(window).height()>450){ jQuery(".outline_switch_td").css({ position : "fixed", top:topp+"px" }); }  }    window.onload=MyAutoRun; $(window).resize(function(){ var bodyw=$win.width(); var _leftPaneInner_width = jQuery(".rich_html_content #leftPaneInner").width(); var _main_article_body = jQuery(".rich_html_content #main_article_body").width(); var rightw=bodyw-_leftPaneInner_width-_main_article_body-25;   var topp=$(window).height()/2; if(rightw<0||$(window).height()<455){ $("#nav-article-page").hide(); $(".outline_switch_td").hide(); }else{ $("#nav-article-page").show(); $(".outline_switch_td").show(); var topp=$(window).height()/2; jQuery(".outline_switch_td").css({ position : "fixed", top:topp+"px" }); } }); 基于二进制GA的B样条重构曲线节点优化
  软件学报  2016, Vol. 27 Issue (10): 2488-2498   PDF    
基于二进制GA的B样条重构曲线节点优化
胡良臣, 寿华好     
浙江工业 大学理学院, 浙江 杭州 310023
摘要: 带法向约束的自由曲线曲面重构在光学反射面设计中起着至关重要的作用.为解决法向约束下的曲线重构问题提出了一种优化方案,使得重构出的曲线在逼近数据点的同时,亦能满足相应法向约束.首先,利用惩罚函数的方法将带法向约束的优化问题转化为无约束的优化问题.然后,引入二进制编码的遗传算法(GA),建立合适的适应度函数,自适应产生优化节点向量,如此迭代进化,直到产生令人满意的重构曲线为止.考虑到节点向量非递减的特性,而遗传算法在寻找最优节点向量的过程中有可能打乱节点向量的顺序,所以在建立适应度函数的时候将变量调整为无序有界变量.通过与传统最小二乘方法和粒子群智能优化方法的比较,所提方案在解决带法向条件约束的曲线重构问题上优势明显,且对于任意形状的曲线重构都行之有效.
关键词: 离散数据     曲线重构     遗传算法     法向约束     B样条     节点向量    
Binary GA for Knot Optimization of B-Spline Curve Reconstruction
HU Liang-Chen, SHOU Hua-Hao     
College of Science, Zhejiang University of Technology, Hangzhou 310023, China
Foundation item: National Natural Science Foundation of China (61572430, 61272309, 61472366)
Abstract: Freedom curve/surface reconstruction with normal constraints is crucial in optical reflecting surface design. In this paper a binary code based genetic algorithm for knot optimization scheme is proposed to reconstruct a B-spline curve that not only approximates the data points but also meets the corresponding normal constraints. First, the constrained optimization problem is transformed into an unconstrained optimization problem by means of penalty function method. Then, the binary code based genetic algorithm (GA) is applied to find the best knot vector after establishing a suitable fitness function. Finally, adaptive generation of optimal knot vector and iterative evolution result in a satisfactory reconstructed curve. Since knot vector is non decreasing, and genetic algorithm may disrupt the order of knot vector in searching for the optimal knot vector, a process is also built to adjust variables into disordered bounded variables in the fitness function. Test results and a comparison with the traditional least square method as well as modern particle swarm optimization method show that the proposed scheme for reconstructing B-spline curve with normal constraints is superior and effective on arbitrary shape of discrete data set.
Key words: discrete data     curve reconstruction     genetic algorithm     normal constraints     B-spline     knot vector    

随着计算机辅助设计的流行, 逆向工程(reverse engineering)[1]蓬勃发展, 数据拟合的相关方法越来越多地应用于几何造型以及图像分析等诸多领域.通过激光扫描、结构光源转换仪或者X射线断层成像等技术获取测量数据, 利用这些离散数据进行数据拟合, 可以对原模型或产品进行模型重建及功能恢复.但有些情况下, 获取的数据点除了坐标信息以外还可能包含一些形状约束条件, 例如在光学工程领域, 通过对带有法向约束的数据点的重构来设计光学反射面.

由于B样条曲线的各种优良的属性[2], 本文采用B样条曲线拟合带法向约束的离散数据点方法.B样条理论是由Schoenberg在1946年提出来的, 1972年De Boor、Cox分别独立给出了B样条的标准算法[2].Gordon和Riesenfeld把B样条应用于形状的描述, 从而提出了B样条方法.B样条基取代Bernstein基, 构造出B样条曲线, 在继承Bézier方法的一切优点的同时克服了Bézier方法存在的缺点, 很好地解决了局部控制的问题, 又轻而易举地解决了连续性拼接问题, 从而使得自由曲线曲面的构造成为可能.B样条方法在表示和设计自由曲线曲面上具有极大的优越性, 它不仅是最流行的形状数学描述的主流方法之一, 而且已成为关于工业产品的几何定义国际标准[2]的有理B样条方法的基础.这些也是我们本文中将其作为带法向数据点逼近函数的重要原因.

对于型值点带有法向、切矢亦或者是曲率等条件约束的曲线插值或逼近问题, 已经有很多研究成果.诸如Piegl等人[2]提出的点-切向约束的B样条曲线插值算法, Gofuku等人[3]提出的点-切向/点-法向的B样条曲线插值几何算法, Abbas等人[4]提出的满足数据点及法向、曲率约束的3次B样条插值的构造算法, Shoichi等人[5]提出的切向曲率约束的3次均匀B样条插值的几何算法, 星蓉生等人[6]提出的3次均匀B样条曲线插值数据点及其切矢的PIA算法以及胡巧莉等人[7]提出的基于几何构造的法向无偏移的带法向约束的3次均匀B样条曲线插值算法.

除了传统的基于反求线性方程系统、几何构造等方法来求解曲线拟合问题外, 利用各种人工智能算法求解以及优化曲线拟合问题是当今一大热门, 比如Yoshimoto等人[8]利用实编码的遗传算法求解无约束的拟合数据的样条曲线, 该算法主要是在曲线拟合的基础上, 优化节点向量的个数与位置以减小拟合误差.在本文中, 我们把遗传算法用到法向约束下曲线拟合问题中, 优化内节点位置.鉴于[0, 1]的求解空间较小, 本文采用二进制编码方式相对实编码来说更加合适.Adi等人[9]用粒子群优化(PSO)方法进行NURBS曲线拟合, 取代解方程组求解控制点, 把控制点作为待求解变量, 但是误差会比较大.Gálvez等人[10]借鉴遗传算法在曲线拟合问题上的应用利用PSO自适应确定节点向量长度以及内节点个数优化拟合曲线.Kang等人[11]将稀疏优化的方法应用到寻找最优拟合曲线节点的问题上, 利用样条曲线p次导数为分段常数这一信息, 建立稀疏优化模型求解最佳节点数量与位置.相对遗传算法来说, 稀疏优化在计算上更加复杂, 但在求解最佳节点数量上优于遗传算法.

受单纯曲线拟合中节点优化的启发, 我们利用二进制编码的遗传算法(GA)解决带有法向约束的曲线逼近问题.众所周知, 一条好的拟合曲线关键在于节点向量的确定, 一个好的节点向量对于实现完美的曲线逼近起到至关重要的作用.本文将数据点带有法向约束的曲线逼近问题转化为基于二进制编码遗传算法的无约束最优化问题, 将节点向量自由化, 自适应地寻找最优节点向量, 最终产生最优逼近曲线.具体流程框架如图 1所示.

Fig. 1 Flowchart of the GA-based B-spline curve reconstruction algorithm with normal constraint 图 1 基于GA的法向约束下样条曲线重构算法流程图

1 问题描述及模型建立

由于B样条曲线所具有的几何不变性、凸包性、局部支撑性、变差减缩性等诸多的优良性质, 所以我们取B样条曲线作为数据点的重构曲线.

一条p次B样条曲线可以表示为

$ S(u) = \sum\limits_{j = 0}^n {{N_{j, p}}(u){P_j}} $ (1)

它的一阶导矢为

$ S'(u) = p\sum\limits_{j = 0}^n {\left( {\frac{{{N_{j, p-1}}(u)}}{{{u_{j + p}}-{u_j}}}-\frac{{{N_{j + 1, p - 1}}(u)}}{{{u_{j + p + 1}} - {u_{j + 1}}}}} \right){P_j}} $ (2)

其中, $ {\rm{\{ }}{P_j}{\rm{\} }}, j = 0, 1, ..., n$是样条逼近曲线的n+1个控制顶点, ${N_{j, p}}(u) $表示p次B样条基函数, 其定义由de Boor和Cox递推[2]给出:

$ {N_{j, 0}}(u) = \left\{ \begin{array}{l} 1{\rm{ }}, {\rm{ }}{u_j} \le u < {u_{j + 1}}\\ 0{\rm{ }}, {\rm{ }}其他 \end{array} \right. $ (3)
$ {N_{j,p}}(u) = \frac{{u - {u_j}}}{{{u_{j + p}} - {u_j}}}{N_{j,p - 1}}(u) + \frac{{{u_{j + p + 1}} - u}}{{{u_{j + p + 1}} - {u_{j + 1}}}}{N_{j + 1,p - 1}}(u),p \ge 1 $ (4)

规定$\frac{0}{0} = 0 $.式中p表示B样条次数, ${u_j} $为节点向量 ${\bf{U}}{\rm{\{ }}{u_j}{\rm{\}, }}j = 0, 1, ..., n + p + 1$中的元素.

给定M+1个数据点di, i=0, ..., M及其法向约束条件li, i=0, ..., M, 找一条曲线S(t)去逼近给定的数据di, i=0, ..., M且对于数据点在曲线上的对应参数处满足法向约束条件.假定与数据点di, i=0, ..., M对应的参数值是ti, 则带法向约束的B样条曲线逼近问题可以用两个方程组成的方程组表示.

$\left\{ {\begin{array}{*{20}{c}} {{d_i} = S({t_i}) + {\varepsilon _i}}\\ {S'({t_i}) \cdot {l_i} = 0{\rm{ }}} \end{array}} \right., {\rm{ }}i = 0, ..., M$ (5)

其中, εi为数据点di的拟合误差, S′(ti), i=0, 1, …, M为逼近曲线在参数值ti处的一阶导矢.

由B样条曲线的定义式(1)~式(4), 假设控制顶点集{Pi}, i=0, 1, …, n和节点向量U{uj}, j=0, 1, …, n+p+1为未知变量, 则式(5)中的带法向约束的B样条曲线逼近问题可转化为带有等式约束的高度非线性且多维多变量的最小值优化问题.

$ \left\{ \begin{array}{l} Q({P_0}, {P_1}, ..., {P_n}, U) = \min \left( {\sum\nolimits_{i = 0}^M {{{\left\| {{d_i}-S({t_i})} \right\|}^2}} } \right) = \min \left( {\sum\nolimits_{i = 0}^M {{{\left\| {{d_i}-\sum\nolimits_{j = 0}^n {{N_{j, p}}({t_i}){P_j}} } \right\|}^2}} } \right)\\ S'({t_i}) \cdot {l_i} = \left( {p\sum\nolimits_{j = 0}^n {\left( {\frac{{{N_{j, p-1}}({t_i})}}{{{u_{j + p}} - {u_j}}} - \frac{{{N_{j + 1, p - 1}}({t_i})}}{{{u_{j + p + 1}} - {u_{j + 1}}}}} \right){P_j}} } \right) \cdot {l_i} = 0, i = 0{\rm{, 1}}, ..., M \end{array} \right.{\rm{ }} $ (6)

其中, ${\left\| \cdot \right\|^2} $为平方范数.

2 法向约束下B样条曲线重构的实现

由于使用传统的计算方法很难实现对如式(6)所示的带约束的曲线逼近问题求解, 所以本文提出了一种解决此类问题的方案:首先使用一种在优化模型上加惩罚函数的方法将带约束的优化模型式(6)转化为无约束的问题, 继而建立一个合理而有效的适应度函数; 其次利用遗传算法自由化节点向量, 利用遗传算法所产生的节点向量为先决条件, 使用最小二乘法产生最优控制顶点, 如此循环迭代, 最终产生最优的节点向量以及控制顶点.

2.1 无约束问题的转化

为了实现式(1)中的等式约束S′(tili=0, 我们将允许严格的等式约束出现适当的误差, 即:S′(tili≈0.这样, 在误差允许的范围内, 可保证给定数据点的法向与逼近曲线的法向一个大致的贴合.这样相比于不带法向约束的曲线重构来说, 会使得重构出来的曲线更加接近原曲线.接下来我们将上式的平方和作为一个惩罚函数, 整理式(5)和式(6), 最终可得一个关于最小值优化的数学模型如下:

$ \begin{array}{l} Q({P_0}, {P_1}, ..., {P_n}, U, {C_0}) = F({P_0}, {P_1}, ..., {P_n}, U) + {C_0}G({P_0}, {P_1}, ..., {P_n}, U)\\ {\rm{ }} = \sum\limits_{i = 0}^M {{{\left\| {{d_i}-\sum\limits_{j = 0}^n {{N_{j, P}}({t_i}){P_j}} } \right\|}^2}} + {C_0}\sum\limits_{i = 0}^M {{{({l_i} \cdot S'({t_i}))}^2}} \end{array} $ (7)

C0作为惩罚系数, 在这里设为1.我们将利用遗传算法对式(7)变换后的适应度函数寻求最优解.若在数据点和法向处均能精确插值, 那么式(7)将有最小值0.

2.2 传统最小二乘法[2]及其解法改进

假设节点向量已知的情况下, 对式(6)中关于数据点误差式子中的变量{Pl}, l=1, …, n-1求偏导数:

$ \frac{{\partial f}}{{\partial {P_l}}} = \sum\limits_{i = 1}^{M-1} {\left( {-2{N_{l, p}}({t_i}){R_i} + 2{N_{l, p}}({t_i})\sum\limits_{j = 1}^{n-1} {{N_{j, p}}({t_i}){P_j}} } \right)}, $

其中, $ {R_i} = {d_i} - {N_{0, p}}({t_i}){d_0} - {N_{n, p}}({t_i}){d_M}.$

令偏导数等于0, 可得:

$ \sum\limits_{i = 1}^{M-1} {\sum\limits_{j = 1}^{n-1} {{N_{l, p}}({t_i}){N_{j, p}}({t_i}){P_j} = \sum\limits_{i = 1}^{M-1} {{N_{l, p}}({t_i}){R_i}} } } . $

l=1, 2, ..., n-1时, 上式成立, 则可得线性方程组:

$ ({{\bf{N}}^{\rm{T}}}{\bf{N}}){\bf{P}} = {{\bf{N}}^{\rm{T}}}{\bf{R}}. $

N是(M-1)×(n-1)的基函数矩阵:

$ {\bf{N}} = \left[ {\begin{array}{*{20}{c}} {{N_{1.p}}({t_1})}& \cdots &{{N_{n - 1.p}}({t_1})}\\ \vdots & \ddots & \vdots \\ {{N_{1.p}}({t_{M - 1}})}& \cdots &{{N_{n - 1.p}}({t_{M - 1}})} \end{array}} \right], $
$ {\bf{R}} = \left[{\begin{array}{*{20}{c}} {{R_1}}\\ \vdots \\ {{R_{M-1}}} \end{array}} \right], {\bf{P}} = \left[{\begin{array}{*{20}{c}} {{P_1}}\\ \vdots \\ {{P_{n-1}}} \end{array}} \right]. $

传统方法求解${\bf{P}} = {({{\bf{N}}^{\rm{T}}}{\bf{N}})^{-1}}{{\bf{N}}^{\rm{T}}}{\bf{R}} $, 需要刻意构造节点向量以至使得${{\bf{N}}^{\rm{T}}}{\bf{N}} $正定且条件数(判断矩阵病态与否的一种度量, 条件数越大越病态)是好的.但是本文中由于节点向量是不断变化的, 随着遗传算法的迭代并不能保证${{\bf{N}}^{\rm{T}}}{\bf{N}} $正定, 所以传统意义上, 采用形如$ {\bf{P}} = {({{\bf{N}}^{\rm{T}}}{\bf{N}})^{-1}}{{\bf{N}}^{\rm{T}}}{\bf{R}}$的方法并不能求解出P.为此, 本文将采取另一种思路, 即对矩阵N求Moore-Penrose广义逆N+, 这样无论节点向量如何变化都能保证方程有且有唯一的极小最小二乘解:P=N+R.

N+的具体求法分情况讨论如下:

(1)如果N为满秩矩阵, N+=N-1;

(2)如果N为行满秩矩阵, N+=NT(NNT)-1;

(3)如果N为列满秩矩阵, N+=(NTN)-1NT;

(4)如果N为降秩矩阵, 将N满秩分解为${\bf{N}} = {\bf{CD}}, {\bf{C}} \in {{\bf{C}}^{(M-1) \times r}}, {\bf{D}} \in {{\bf{D}}^{r \times (n-1)}} $, 然后分别求分解后矩阵的左逆和右逆:${\bf{C}}_{\rm{L}}^{-1} = {({{\bf{C}}^{\rm{T}}}{\bf{C}})^{-1}}{{\bf{C}}^{\rm{T}}}, {\bf{D}}_{\rm{R}}^{-1} = {{\bf{D}}^{\rm{T}}}{({\bf{D}}{{\bf{D}}^{\rm{T}}})^{ - 1}}, $最后可得极小最小二乘解${{\bf{N}}^ + } = {\bf{D}}_{\rm{R}}^{-1}{\bf{C}}_{\rm{L}}^{-1}. $

在数值实验中, 通过与传统最小二乘法对比, 将会发现本文提出的法向约束下的曲线重构方法更加有效.

2.3 节点向量的设置

在最小二乘拟合下, 影响一条完整的B样条重构曲线形状的因素有两点:(a)节点向量的选取; (b)数据点参数化的分布情况.不同的节点向量和数据点参数化的选择对曲线的形状有很大的影响, 不合理的节点和参数化将导致曲线逼近程度变坏.一般来说, 数据点参数化可以采用向心参数化方法[12], 这种方法也得到了广泛的认可.但是对于节点向量的选取, 时至今日, 解决这一问题的方法虽然很多, 但最佳的方案仍尚待研究.鉴于这种原因, 本文中针对节点向量的最优化选取做出了自适应的调整.

对于节点向量的选取, 具体方案如下:首先将控制点端点进行插值, 也就是说对于一个非递减的节点向量$ {\bf{U}}{\rm{\{ }}{u_k}{\rm{\} }}_{k = 0}^{n + p + 1}, {u_0} \le {u_1} \le ... \le {u_{n + p + 1}}$来说, 两端节点的重复度设置为B样条曲线的阶数p+1:${u_0} = {u_1} = ... = {u_p} = 0, $; 然后将首尾两端的数据点参数化设为0和1.如此首尾两端的数据点将成为重构B样条曲线的首尾控制顶点: ${d_0} = {P_0}, {d_M} = {P_n} $; 最后, 对于节点向量来说, 未知量就只剩下n+p个内节点, 对于这些未知的节点向量, 我们将建立与数据误差和法向误差相关的适应度函数, 然后利用遗传算法对其进行自适应的调整, 直到产生最佳的节点向量为止, 使得待重构曲线在逼近数据点的同时, 亦能满足数据点处的法向条件约束.

2.4 数据点参数化设置

利用参数曲线进行数据点重构, 首先必须决定数据点在参数曲线上的分布情况, 也就是数据点的参数化.数据点参数化方法有很多, 采取不同的参数化方法将会对重构曲线的形状产生较大的影响.这里我们采用向心弦长参数化[12]方法.相对于其他参数化方法, 这种方法在很好地模拟数据点分布情况的同时, 当数据点急剧变化时, 能够使得重构曲线效果更佳.参数ti通过数据点之间的位置关系进行估计, 估计公式为

$ \left\{ {\begin{array}{*{20}{c}} {{t_0} = 0, {\rm{ }}}\\ {{t_{i + 1}} = {t_i} + \frac{{\sqrt {\left\| {{d_{i + 1}}-{d_i}} \right\|} }}{{\sum\limits_{j = 0}^M {\sqrt {\left\| {{d_{j + 1}}-{d_j}} \right\|} } }}} \end{array}} \right., {\rm{ }}i = 0, 1, ..., M $ (8)

其中, 参数ti就是数据点di在其B样条逼近曲线上对应的参数值.

3 遗传算法原理及相关设置

遗传算法(genetic algorithm)是一种借鉴生物界物种在繁衍生存过程中的遗传进化规律(适者生存, 优胜劣汰的机制)演化而来的随机全局优化搜索方法.遗传算法的基本理论与基本方法起初是由美国的Holland教授1975年提出来的, 其最为突出的特点是直接对可行解本身进行仿生操作, 不存在函数的求导以及连续性等条件的限定; 具有与生俱来的并行性结构和较好的全局搜索寻优的能力; 自动确定和进行优化的搜索空间, 自适应地调整搜索的方向, 不需要人为干预以及确定规则.遗传算法的这些优良的性质, 已被广泛地应用于自适应控制、信号处理、人工生命、机器学习和组合优化等领域.它是现代有关智能计算中的关键技术.

现行的编码方式有很多, 例如:二进制编码、浮点编码(即实编码)、符号编码等.二进制编码方法是遗传算法使用最广泛的编码方式之一, 它使用二值符号集{0, 1}, 它所构成的个体基因型是一个二进制编码符号串.然后由多个个体基因型组成染色体(可行解), 二进制编码符号串的长度间接代表问题所要求的求解精度.二进制编码具有编码、解码简单易行, 交叉、变异等算子的操作易于实现等优点.符号编码适用于特定的问题, 通用性不够.浮点编码适用于待搜索的解空间较大的情况, 对于较小的搜索空间, 采用二进制编码效果比较好.在本文中由于与节点相关的优化问题中变量的取值范围为[0, 1], 求解空间比较小, 所以本文中采纳二进制编码的方式是合理的.

3.1 相关算子以及参数设置

下面我们给出本文中所采用的编码、个体评价方式以及各算子的操作.

下面我们给出本文中所采用的编码、个体评价方式以及各算子的操作.

a.编码:对于一个完整的n维可行解$\Im = {\rm{\{ }}\Im _1^{}, \Im _2^{}, ..., \Im _n^{}{\rm{\}, }} $其中每一维度的分量${\Im _i}, i = 1, ..., n $的取值范围为$[{a_i}, {b_i}] $.如果我们要求的精度是小数点后m位, 那么每个分量$ {\Im _i}, i = 1, ..., n$都要至少被分为$({b_i}-{a_i}) \times {10^m} $个部分.对于一个分量的二进制编码串位数(用${w_i} $来表示), 用如下公式来计算:

$ {2^{{w_i}}} < (b-a) \times {10^m} \le {2^{{w_i}}}-1. $

鉴于可行解的每一维分量的取值范围均为[0, 1], 搜索空间比较狭窄, 可行解需要较高的精度, 故编码位数取得比较大.在所有的实验中分量编码的位数均取为25, 那么对于一个完整的可行解$\Im $编码后的形式如下:

$ \begin{align} & \text{ }\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{25}\times n位 \\ & \overleftrightarrow{0100101011100101110100001,...,1000110101001110000110101} \\ & \overleftrightarrow{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{25}位\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\ \ \ \ \ \ \ \ \ \overleftrightarrow{\text{ }\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{25}位\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ } \\ \end{align} $

b.解码:从二进制转换成一个可用的十进制数值的方式如下:

$ {\Im _i} = {a_i} + (二制代表的十制值进实际进数) \times \frac{{{b_i}-{a_i}}}{{{2^{{w_i}}}-1}}. $

c.个体评价:适应度函数的建立是为了评价群体每个个体的好坏程度.由于在遗传算法进化过程中不能保证节点向量的有序性, 而一个无序的节点向量对于B样条曲线来说是不成立的, 所以结合式(7), 本文中特意设计适应度函数来规避这个问题, 具体如下:

$ \begin{array}{l} f(\Im )\\ \{ \\ {\rm{ }}U = sort(\Im );{\rm{ //}}对\Im各分量进行排序 \\ {\rm{ }}P = solve(U);{\rm{ //}}最小二乘法求控制点\\ {\rm{ }}val = - Q({P_0},{P_1},...,{P_n},U,{C_0});{\rm{ //}}式(7)取负值转化为最大值问题\\ {\rm{ return }}val;\\ \} \end{array} $

如上操作的目的是, 将变量排序刻画在函数体内与优化问题形成一个函数整体, 那么在遗传算法进化过程中就不必再考虑变量的顺序问题了.

与其他编码方式一样, 二进制编码的遗传算法同样包括了3个基本操作, 选择、交叉以及变异.

(1)选择:本文中设置种群代沟为0.9, 通过轮盘赌选择上一代种群中的90%的个体复制到下一代中.

(2)交叉:这里我们选择交叉概率为0.7的单点交叉方式, 并对本代个体数量交叉直到恢复原始种群规模.

(3)变异:我们取变异概率为pm=0.7/length, length为染色体结构的长度.

3.2 本文具体算法描述

总结以上关于遗传算法以及相关参数的概述, 给出以下关于带有法向约束的曲线拟合算法的描述.

Step 1:输入待拟合数据di, i=0, ..., M及其对应法向li, i=0, ..., M.

Step 2:数据点参数化, 以及有效的适应度函数(如第3.1.节c中描述)的建立.

Step 3:遗传算法相关参数的设置, 包括种群规模PSIZE、代沟Pg、交叉概率Pc、变异概率Pm、最大迭代次数MITER

Step 4:初始化二进制编码的种群$\zeta (t) = \left\{ {{\zeta _i}} \right\}_{i = 0}^{PSIZE-1}. $

Step 5:解码并通过最小二乘法计算控制顶点, 计算初始种群$\left\{ {{\zeta _i}} \right\}_{i = 0}^{PSIZE-1} $中每一个个体的适应度值.

Step 6:判断迭代终止的条件, 如果满足条件, 则输出最优节点向量和控制顶点并计算重构曲线, 否则转到Step 7继续进行.

Step 7:对ζ(t)按照轮盘赌方式执行选择算子产生ζ1(t).

Step 8:执行交叉算子单点交叉对ζ1(t)进行操作产生ζ2(t).

Step 9:变异算子对ζ2(t)进行操作产生ζ(t+1), 继续Step 5.

4 数值实验

例1:以参数曲线$\left\{ {\begin{array}{*{20}{c}} {x = r\cos (3t)}\\ {y = r\sin (3t)} \end{array}} \right. $作为采样函数, 其中, $ r = 10(1 + t), $在区间[0.2π]上以0.05为间隔均匀取点和相应法向, 取点个数为126, 种群的规模取40, 最大迭代次数为300, 样条曲线次数为3.

例2:以参数曲线$ \left\{ {\begin{array}{*{20}{c}} {x = {\rm{sin}}(0.75t)}\\ {y = \sin (t){\rm{ }}} \end{array}} \right.$作为采样函数, 在区间[-4π, 4π]上以0.2为间隔均匀取点和相应法向, 取点个数为126, 种群规模取40, 最大迭代次数为300, 样条曲线次数为3.

例3:以参数曲线$ \left\{ \begin{array}{l} x = 2(5\cos (t)-\cos (6t))\\ y = 2(3\sin (t)-\sin (4t)) \end{array} \right.$作为采样函数, 在区间$ [0, 2{\rm{\pi }}]$上以0.05为间隔均匀取点和相应法向, 取点个数为126, 种群规模取40, 最大迭代次数为300, 样条曲线次数为3.例4:以参数曲线$ \left\{ {\begin{array}{*{20}{c}} {x = (a + b)\cos (t)-h\cos \left( {\frac{{a + b}}{b}t} \right)}\\ {y = (a + b)\sin (t)-h\sin \left( {\frac{{a + b}}{b}t} \right)} \end{array}} \right.$作为采样函数, 其中, a=3, b=1/3, h=1, 在区间$[-{\rm{\pi }}, {\rm{\pi }}] $上以0.05为间隔均匀取点和相应法向, 取点个数为126, 种群规模取40, 最大迭代次数为300, 样条曲线次数为3.

本文选取了4个较为复杂的实例(如图 2所示), 第1个例子中数据点间隔变化比较均匀, 第2个例子和第4个例子是出于其多交叉点的考虑, 第3个例子局部变化较为尖锐, 这4个例子均具有一般的代表性.通过比较具有代表性的例子对所提出的算法的可行性进行了验证, 每个实例均做10次相同实验, 然后取均值并在表 1~表 4中分别列出具体数值.图 3图 4图 5图 6图 7图 8图 9图 10分别为4个实例的遗传算法优化节点向量的收敛图和曲线重构效果图, 据此可以比较直观地看出4个例子的实际收敛效果和拟合效果, 实际验证了本文算法对于解决此类问题是可行的.图 3图 5图 7图 9中种群均值的变化和解的变化表明不同的节点向量对曲线重构的影响还是相当明显的, 算法迭代的前期均能比较快速地收敛, 但是对于后期, 随着种群越来越靠近最优值, 收敛速度逐渐变慢, 所以想获取更高的曲线重构精度就必须尽量延长算法的迭代次数或者增加节点向量的变量数目.第2.2节中所描述的传统的最小二乘法为了保证NTN的正定性, 所以它的节点向量需要特意去构造以保证正定性, 但是一旦将节点向量固定也就意味着B样条曲线即固定, 丧失了自由性, 那么曲线的自由性就有所限制, 这样重构出来的曲线误差就会比较大.本文中将节点向量自由化之后, 通过遗传算法寻求节点数目(节点数目可人为自由选取)确定之后的最优节点向量, 加之以法向条件的约束, 使得与原扫描模型最匹配的重构曲线被构造出来.通过表 1~表 4中的数据表明, 本文所提出的节点优化的方案确实显示出了比较明显的优越性.

Fig. 2 Test examples 图 2 实验实例

Fig. 3 GA convergence effect for example 1 图 3 例1的GA收敛效果图

Fig. 4 GA-Based B-spline reconstruction with normal constraint for example 1 图 4 例1基于GA的法向约束下B样条重构效果图

Fig. 5 GA convergence effect for example 2 图 5 例2的GA收敛效果图

Fig. 6 GA-Based B-spline reconstruction with normal constraint for example 2 图 6 例2基于GA的法向约束下B样条重构效果图

Fig. 7 GA convergence effect for example 3 图 7 例3的GA收敛效果图

Fig. 8 GA-Based B-spline reconstruction with normal constraint for example 3 图 8 例3基于GA的法向约束下B样条重构效果图

Fig. 9 GA convergence effect for example 4 图 9 例4的GA收敛效果图

Fig. 10 GA-Based B-spline reconstruction with normal constraint for example 4 图 10 例4基于GA的法向约束下B样条重构效果图

Table 1 The comparison of average error for example 1 表 1 例1的平均误差比较

Table 2 The comparison of average error example 2 表 2 例2的平均误差比较

Table 3 The comparison of average error for example 3 表 3 例3的平均误差比较

Table 4 The comparison of average error for example4 表 4 例4的平均误差比较

另外, 为了凸显遗传算法的优越性, 我们还利用粒子群算法(PSO)对4个例子做了同样的实验.关于本实验PSO的相关参数的设置:与遗传算法一致, 种群规模取40, 迭代次数取300, 对于PSO特有的学习因子、惯性权重, 这里分别取2.0和0.7.由于在迭代过程中种群的个体可能会发生越界的情况(变量要求在[0, 1]之间), 我们将越界个体重新随机初始化[0, 1]之间的某个位置.通过表 1~表 4中的数据可以看出, 遗传算法相对粒子群算法来说在求解此类问题上较为优越.

在实验过程中, 上述所有实例的数据点参数化方法均选用向心参数化.为了显示不同的参数化方法对以上各实例重构效果的影响情况, 利用本文的方法在相同的实验条件下, 我们分别选取累加弦长参数法、均匀参数化方法以及本文所采纳的向心参数化方法进行误差比较.如表 5所示, 我们可以看出, 不同的参数化方法确实对实验结果有着不同的影响, 从数据来看均匀参数化方法效果最好, 主要原因可能是因为这4个例子的数据点分布取得比较均匀.另外可以看出, 本文所采纳的向心参数化方法比累加弦长参数化方法的重构效果要更好一些.当然, 这仅是对本文中所采纳的4个例子而言, 不同的参数化选择对于其他的重构情况的效果可能又会不一样.参数化的选择很大程度上取决于数据分布情况, 那么对于不同的参数化方法就有不同的适用场合, 这给我们一定的启发, 如果有合适的方法对参数化作自适应选择, 那么重构的效果可能会更加鲜明.

Table 5 The comparison of average error for different parameterization method 表 5 不同参数化方法平均误差比较

法向约束下二维曲线的重构可以快速应用到三维光学反射面的设计上.例如探照灯的聚光镜面设计问题.探照灯聚光镜面是一张旋转曲面, 它是由xOy坐标面上一条带法向约束的曲线绕x轴旋转而成.按照聚光镜性能要求要求, 在其旋转轴上一点O处发射光线, 经它反射后都与旋转轴平行.这里我们给出一个设计例子, 如图 11所示.

Fig. 11 The design of optical reflecting surface 图 11 光学反射面设计示意图

5 结论

带法向约束的自由曲线曲面重构在光学反射面设计中起到十分重要的作用.对于没有约束的曲线重构方法所构造出来的曲线, 虽然在扫描数据点的逼近程度上会比较令人满意, 但是由于没有约束条件, 可能会导致重构出来的曲线丢失原扫描模型的某些特征或性质.本文针对带有法向约束条件的离散数据集所提出来的这种法向约束下的曲线重构算法, 会在逼近数据点的同时满足在数据点处的法向约束, 因而相对于传统的曲线重构会使得重构出来的曲线效果更好.这种方案通过建立数据点和法向处的误差平方和函数, 继而以此作为遗传算法适应度函数的基础, 不断通过遗传算法的进化原理调整节点向量的位置, 使得B样条重构曲线数据误差和法向误差越来越小, 形状与原扫描模型越来越贴合.其中, 自适应的节点向量是本算法的核心步骤.而且本文所提出的方案对于解决法向约束下的任意扫描数据点集的重构都有效.通过实验数据与实验效果来看, 相比于传统的最小二乘方法以及粒子群算法来说, 本文所提方案确实有着比较明显的优越性.

进一步地, 本文所提出的方案可以很方便地推广到带约束的空间曲线或者曲面重构问题上去.带法向约束的二维曲线重构作为三维曲面重构的基础, 可以很自然地将这种思想延伸到带法向约束的三维曲面重构的问题上.而对于三维数据点集的曲线重构问题来说, 鉴于它不具有法向, 应用本文方案稍微加以修改就可以解决带有切向约束条件的三维空间曲线重构问题.这些都将成为未来我们继续探索的方向.

参考文献
[1] Várady T, Martin RR, Cox J. Reverse engineering of geometric models-An introduction. Computer-Aided Design, 1997, 29 (4) :255–268. [doi:10.1016/S0010-4485(96)00054-1]
[2] Piegl L, Tiller W. The NURBS Book. New York: Springer-Verlag, 1997 .
[3] Gofuku S, Tamura S, Maekawa T. Point-Tangent/Point-Normal B-spline curve interpolation by geometric algorithms. Computer-Aided Design, 2009, 41 (6) :412–422. [doi:10.1016/j.cad.2009.02.005]
[4] Abbas A, Nasri A, Maekawa T. Generating B-spline curves with points, normals and curvature constraints:A constructive approach. The Visual Computer, 2010, 26 (6) :823–829. [doi:10.1007/s00371-010-0441-2]
[5] Shoichi O, Ahmad N, Lin HW, Abdulwahed A, Yuki K, Takashi M. Uniform B-spline curve interpolation with prescribed tangent and curvature vectors. IEEE Trans.on Visualization and Computer Graphics, 2012, 18 (9) :1474–1487. [doi:10.1109/TPDS.2011.262]
[6] Xing RS, Pan RJ. PIA for uniform cubic B-spline curve interpolation with prescribed tangent vector. Journal of Fujian Normal University:Nature Science Edition, 2014, 30 (1) :25–32(in Chinese with English abstract). http://www.cnki.com.cn/Article/CJFDTOTAL-FJSZ201401006.htm
[7] Hu QL, Shou HH. Uniform cubic B-spline curve interpolation with normal constraint. Journal of Zhejiang University (SCIENCE EDITION), 2014, 41 (6) :619–623(in Chinese with English abstract). [doi:10.3785/j.issn.1008-9497.2014.06.002]
[8] Yoshimoto F, Harada T, Yoshimoto Y. Data fitting with a spline using a real-coded genetic algorithm. Computer-Aided Design, 2003, 35 (8) :751–760. [doi:10.1016/S0010-4485(03)00006-X]
[9] Adi DI, Shamsuddin SM, Ali A. Particle swarm optimization for NURBS curve fitting. In: Proc. of the 6th Int'l Conf. on Computer Graphics, Imaging and Visualization (CGIV 2009). Washington: IEEE Computer Society Press, 2009. 259–263. [doi: 10.1109/ CGIV.2009.64]
[10] Gálvez A, Iglesias A. Efficient particle swarm optimization approach for data fitting with free knot B-splines. Computer-Aided Design, 2011, 43 (12) :1683–1692. [doi:10.1016/j.cad.2011.07.010]
[11] Kang HM, Chen FL, Li YS, Deng JS, Yang ZW. Knot calculation for spline fitting via sparse optimization. Computer-Aided Design, 2014, 58 (8) :179–188. [doi:10.1016/j.cad.2014.08.022]
[12] Lee ETY. Choosing nodes in parametric curve interpolation. Computer-Aided Design, 1989, 21 (6) :363–370. [doi:10.1016/0010-4485(89)90003-1]
[6] 星蓉生, 潘日晶. 三次均匀B样条曲线插值数据点及其切矢的PIA算法. 福建师范大学学报:自然科学版, 2014,30 (1) :25–32. http://www.cnki.com.cn/Article/CJFDTOTAL-FJSZ201401006.htm
[7] 胡巧莉, 寿华好. 带法向约束的3次均匀B样条曲线插值. 浙江大学学报(理学版), 2014,41 (6) :619–623. [doi:10.3785/j.issn.1008-9497.2014.06.002]