B样条曲线计算

1、B样条基本概念

1.1 B样条表示

其中，为控制顶点（共n+1个控制点），控制点对应的参数序列即节点矢量。

根据节点矢量的节点分布情况不同，可分为均匀B样条曲线（节点矢量中节点沿参数轴等距分布，所有节点区间长度），准均匀B样条曲线（两端节点重复度k+1，内节点均匀分布），一般非均匀B样条曲线（节点矢量任意分布，只需要在数学上成立，节点序列非递减，两端节点重复度≤k+1，内节点重复度≤k，节点个数n+k+2）。

对于开曲线，一般把两端节点取重复度k+1，即。剩下只需确定内节点值，共n-k个。

1.2 节点向量设置

确定内节点值的Hartley-Judd算法

其中，为控制多边形的边长，，于是得到节点值：

2、B样条曲线求值和求导的德布尔算法

理论：

给定控制顶点、次数k、节点矢量后，就定义了一条k次B样条曲线。若给出曲线定义域内一参数值，可以用德布尔递推方法计算出曲线上对应的点。

德布尔算法的几何意义：从控制多边形开始进行k-1层割角。如下图所示。

求值公式为：(这里的r是重复度)

求r阶导矢：

先递推计算r阶导矢下的控制点，再按照k-r阶曲线求值。经过r级递推，计算出第r级的k-r+1个中间顶点，然后计算这些顶点与相关节点所定义的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(k+1) i = k+1; end %% 3. 一次B样条特殊情况 if (k==1 && r==1) && ... (u~=aNode(i) || (u==aNode(i) && KnotMark==0)) px = (xVertex(i)-xVertex(i-1)) / ... (aNode(i+1)-aNode(i)); py = (yVertex(i)-yVertex(i-1)) / ... (aNode(i+1)-aNode(i)); return; end if (k==1 && r==1 && KnotMark==1) 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 i>k+1 && i<=length(xVertex) Multiple = 1; ii = i; while u == aNode(ii-1) Multiple = Multiple + 1; ii = ii - 1; end % 对于节点重复的情况(例如u4=u5=u6)，选择左边的导矢（i=4） if k>1 && KnotMark==1 && ... u==aNode(i) && i~=k+1 i = i - Multiple; end end %% 5. 取影响当前点的k+1个控制点 tempx = zeros(k+1,1); tempy = zeros(k+1,1); index = 1; for j=i-k:i tempx(index)=xVertex(j); tempy(index)=yVertex(j); index=index+1; end %% 6. 求r阶导控制点 if k>1 for l=1:r for j=i-k:i-l jj=j-(i-k)+1; du = aNode(j+k+1)-aNode(j+l); beta = (k-l+1)/du; tempx(jj)=beta*(tempx(jj+1)-tempx(jj)); tempy(jj)=beta*(tempy(jj+1)-tempy(jj)); end end end %% 7. de Boor算法 for l=1:(k-r) for j=i-k:i-l-r jj=j-i+k+1; du=aNode(j+k+1)-aNode(j+r+l); if du==0 alpha=0; else alpha=(u-aNode(j+r+l))/du; end tempx(jj)=... (1-alpha)*tempx(jj)+... alpha*tempx(jj+1); tempy(jj)=... (1-alpha)*tempy(jj)+... alpha*tempy(jj+1); end end %% 8. 输出 px=tempx(1); py=tempy(1); end

3、B样条最小二乘拟合

理论：

预先计算数据点的参数值和节点矢量。设给定m+1个有序数据点，，…，逼近曲线的次数为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：数据点坐标（每一列是一个数据点） %% 数据点弦长参数化 du=sqrt(diff(points(1,:)).^2+diff(points(2,:)).^2); points_u=[0 cumsum(du)]; points_u=points_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 n=k+1; 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-n+1), 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_i+k-u_i) dn = u - t(i+1); dd = t(i+k+1) - t(i+1); % 第一项 if dd ~= 0 % indeterminate forms 0/0 are deemed to be zero y = y + b.*(dn./dd); end % 求下一阶N_i+1,k-1(u) b = k_bspline_basis(i+1,k-1,t,u); % (u_i+k+1-u)/(u_i+k+1-u_i+1) dn = t(i+k+1+1) - u; dd = t(i+k+1+1) - t(i+1+1); % 第二项 if dd ~= 0 y = y + b.*(dn./dd); end elseif t(i+2) < t(end) % treat last element of knot vector as a special case % k=0最低级的时候 % N_i,0 y(t(i+1) <= u & u < t(i+2)) = 1; else % 最后一个节点特殊处理，没有右边界 y(t(i+1) <= 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 n=k+1; if nargin > 2 B = zeros(numel(x),numel(t)-n); for j = 0 : numel(t)-n-1 B(:,j+1) = k_bspline_basis(j,k,t,x); %第j+1列 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(:,j+1) = 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); % 防止重复索引 t=u(index); p=k+1; t = t(:); n = length(t)-p; if p<2, error(message('SPLINES:AVEKNT:wrongk')) elseif n<0, error(message('SPLINES:AVEKNT:toofewknots')) elseif p==2, tstar = reshape(t(1+[1:n]),1,n); else temp = repmat(t,1,p-1); temp = sum(reshape([temp(:);zeros(p-1,1)],n+p+1,p-1).')/(p-1); tstar = temp(1+[1:n]); end U=[zeros(1,k+1),tstar,ones(1,k+1)]; end