1、B样条基本概念1.1 B样条表示其中为控制顶点共n1个控制点控制点对应的参数序列即节点矢量。根据节点矢量的节点分布情况不同可分为均匀B样条曲线节点矢量中节点沿参数轴等距分布所有节点区间长度准均匀B样条曲线两端节点重复度k1内节点均匀分布一般非均匀B样条曲线节点矢量任意分布只需要在数学上成立节点序列非递减两端节点重复度≤k1内节点重复度≤k节点个数nk2。对于开曲线一般把两端节点取重复度k1即。剩下只需确定内节点值共n-k个。1.2 节点向量设置确定内节点值的Hartley-Judd算法其中为控制多边形的边长于是得到节点值2、B样条曲线求值和求导的德布尔算法理论给定控制顶点、次数k、节点矢量后就定义了一条k次B样条曲线。若给出曲线定义域内一参数值可以用德布尔递推方法计算出曲线上对应的点。德布尔算法的几何意义从控制多边形开始进行k-1层割角。如下图所示。求值公式为(这里的r是重复度)求r阶导矢先递推计算r阶导矢下的控制点再按照k-r阶曲线求值。经过r级递推计算出第r级的k-r1个中间顶点然后计算这些顶点与相关节点所定义的k-r次B样条曲线上参数为u的点。代码实现function [px, py] GetBPr(r, k, u, xVertex, yVertex, aNode, KnotMark) %------------------------------------------------------------ % 用德布尔算法求k次平面B样条曲线上参数u处的r阶导矢 % % 输入: % r - 导矢阶数 % k - B样条次数 % u - 参数 % xVertex - 控制点x坐标 % yVertex - 控制点y坐标 % aNode - 节点矢量 % KnotMark - 节点导矢选择 % 0: 右导矢 % 1: 左导矢 % % 输出: % px,py - u处r阶导矢 % %------------------------------------------------------------ %% 1. 判断阶数 if r k px 0; py 0; return; end %% 2. 查找u所在节点区间 nNode length(aNode); i nNode-k-1; %matlab最低索引1 while i 1 if u aNode(i) break; end i i-1; end if u aNode(k1) i k1; end %% 3. 一次B样条特殊情况 if (k1 r1) ... (u~aNode(i) || (uaNode(i) KnotMark0)) px (xVertex(i)-xVertex(i-1)) / ... (aNode(i1)-aNode(i)); py (yVertex(i)-yVertex(i-1)) / ... (aNode(i1)-aNode(i)); return; end if (k1 r1 KnotMark1) px (xVertex(i-1)-xVertex(i-2)) / ... (aNode(i)-aNode(i-1)); py (yVertex(i-1)-yVertex(i-2)) / ... (aNode(i)-aNode(i-1)); return; end %% 4. 计算节点重复度 if ik1 ilength(xVertex) Multiple 1; ii i; while u aNode(ii-1) Multiple Multiple 1; ii ii - 1; end % 对于节点重复的情况(例如u4u5u6)选择左边的导矢i4 if k1 KnotMark1 ... uaNode(i) i~k1 i i - Multiple; end end %% 5. 取影响当前点的k1个控制点 tempx zeros(k1,1); tempy zeros(k1,1); index 1; for ji-k:i tempx(index)xVertex(j); tempy(index)yVertex(j); indexindex1; end %% 6. 求r阶导控制点 if k1 for l1:r for ji-k:i-l jjj-(i-k)1; du aNode(jk1)-aNode(jl); beta (k-l1)/du; tempx(jj)beta*(tempx(jj1)-tempx(jj)); tempy(jj)beta*(tempy(jj1)-tempy(jj)); end end end %% 7. de Boor算法 for l1:(k-r) for ji-k:i-l-r jjj-ik1; duaNode(jk1)-aNode(jrl); if du0 alpha0; else alpha(u-aNode(jrl))/du; end tempx(jj)... (1-alpha)*tempx(jj)... alpha*tempx(jj1); tempy(jj)... (1-alpha)*tempy(jj)... alpha*tempy(jj1); end end %% 8. 输出 pxtempx(1); pytempy(1); end3、B样条最小二乘拟合理论预先计算数据点的参数值和节点矢量。设给定m1个有序数据点…逼近曲线的次数为k寻找一条k次B样条曲线满足端点约束即两端数据点其余数据点在最小二乘意义上被逼近即目标函数是关于n-1个控制顶点的最小值。注意不是在曲线上距离最近的点。设使目标函数最小则使其对于n-1个控制顶点的导数等于零。则于是是一个以为未知量的线性方程。令得到n-1个线性方程。写成矩阵形式代码实现没考虑首末端点的约束实现结果和matlab的spap2一致% B样条最小二乘拟合函数 function [D,U_2] getD_B_approx(NumCtrl,k,points) % 输入 % k曲线次数 % t节点矢量 % x数据点参数序列 % M数据点坐标每一列是一个数据点 %% 数据点弦长参数化 dusqrt(diff(points(1,:)).^2diff(points(2,:)).^2); points_u[0 cumsum(du)]; points_upoints_u/points_u(end); %% 获取节点矢量 U getU(NumCtrl,k,points_u,points); U_2 getU_2(NumCtrl,k,points_u); %% 最小二乘拟合 D k_bspline_approx(k,U_2,points_u,points); end % 计算基函数 % function [y,x] k_bspline_basis(j,k,t,x) % B-spline basis function value B(j,n) at x. % % Input arguments: % j: 第j个基函数 % interval index, 0 j numel(t)-n % k: 次数 % B-spline order (2 for linear, 3 for quadratic, etc.) % t: 节点矢量 % knot vector % x (optional): 所求值的位置 % value where the basis function is to be evaluated % % Output arguments: % y: % B-spline basis function value, nonzero for a knot span of n % Copyright 2010 Levente Hunyadi nk1; validateattributes(j, {numeric}, {nonnegative,integer,scalar}); validateattributes(n, {numeric}, {positive,integer,scalar}); validateattributes(t, {numeric}, {real,vector}); assert(all( t(2:end)-t(1:end-1) 0 ), ... Knot vector values should be nondecreasing.); if nargin 4 x linspace(t(n), t(end-n1), 100); % allocate points uniformly else validateattributes(x, {numeric}, {real,vector}); end assert(0 j j numel(t)-n, ... Invalid interval index j %d, expected 0 j %d (0 j numel(t)-n)., j, numel(t)-n); y k_bspline_basis_recurrence(j,k,t,x); end function y k_bspline_basis_recurrence(i,k,t,u) y zeros(size(u)); if k 0 % 求下一阶N_i,k-1(u) b k_bspline_basis(i,k-1,t,u); % (u-u_i)/(u_ik-u_i) dn u - t(i1); dd t(ik1) - t(i1); % 第一项 if dd ~ 0 % indeterminate forms 0/0 are deemed to be zero y y b.*(dn./dd); end % 求下一阶N_i1,k-1(u) b k_bspline_basis(i1,k-1,t,u); % (u_ik1-u)/(u_ik1-u_i1) dn t(ik11) - u; dd t(ik11) - t(i11); % 第二项 if dd ~ 0 y y b.*(dn./dd); end elseif t(i2) t(end) % treat last element of knot vector as a special case % k0最低级的时候 % N_i,0 y(t(i1) u u t(i2)) 1; else % 最后一个节点特殊处理没有右边界 y(t(i1) u) 1; end end function [B,x] k_bspline_basismatrix(k,t,x) % B-spline basis function value matrix B(n) for x. % % Input arguments: % n: % B-spline order (2 for linear, 3 for quadratic, etc.) % t: % knot vector % x (optional): % an m-dimensional vector of values where the basis function is to be % evaluated % % Output arguments: % B: % a matrix of m rows and numel(t)-n columns % Copyright 2010 Levente Hunyadi nk1; if nargin 2 B zeros(numel(x),numel(t)-n); for j 0 : numel(t)-n-1 B(:,j1) k_bspline_basis(j,k,t,x); %第j1列 end else [b,x] k_bspline_basis(0,k,t); B zeros(numel(x),numel(t)-n); B(:,1) b; for j 1 : numel(t)-n-1 B(:,j1) k_bspline_basis(j,k,t,x); end end end function D k_bspline_approx(k,t,x,M) % B-spline curve control point approximation with known knot vector. % % Input arguments: % k: % B-spline order (2 for linear, 3 for quadratic, etc.) % t: % knot vector % x: % B-spline values corresponding to which data points are observed % M: % d-by-m matrix of observed data points, possibly polluted with noise, % d is typically 2 for plane, 3 for space, or 3 or 4, respectively, if % weights are present % % Output arguments: % D: % d-by-n matrix of control points % Copyright 2010 Levente Hunyadi validateattributes(k, {numeric}, {positive,integer,scalar}); validateattributes(t, {numeric}, {real,vector}); validateattributes(x, {numeric}, {real,vector}); validateattributes(M, {numeric}, {real,2d}); B k_bspline_basismatrix(k,t,x); Q M * B; D Q / (B*B); end function U getU_2(NumCtrl,k,u) % 输入 % NumCtrl控制点数量 % k曲线次数 % u数据点参数序列 % points数据点坐标 % 数据点数量 num_dataPoint length(u); % 从数据点均匀挑选 index round(linspace(1,num_dataPoint,NumCtrl)); index unique(index); % 防止重复索引 tu(index); pk1; t t(:); n length(t)-p; if p2, error(message(SPLINES:AVEKNT:wrongk)) elseif n0, error(message(SPLINES:AVEKNT:toofewknots)) elseif p2, tstar reshape(t(1[1:n]),1,n); else temp repmat(t,1,p-1); temp sum(reshape([temp(:);zeros(p-1,1)],np1,p-1).)/(p-1); tstar temp(1[1:n]); end U[zeros(1,k1),tstar,ones(1,k1)]; end