3

我正在尝试在 Matlab 中实现我自己的代码,以将封闭的 b 样条拟合到一组 2D 数据。

我以椭圆的形式生成了二维数据并使用了 De boor 算法。为开放 b 样条设置代码时,该代码确实可以正常工作,但它被封闭(周期性)b 样条卡住了。

以下是步骤:

  1. 以二维椭圆的形式生成有序二维数据C [2xn]
  2. 根据均匀间隔或弦长方法生成参数向量t [nx1]。
  3. 生成节点向量u [m+dx1] (:# control points, d: b-spline order): a) Open b-spline: uniform-clamped method b) closed b-spline: uniform
function u = npa_contfit_genknot(t,d,nmp,opt,varargin)
   switch opt
      case 'uniform-clamped'
         a = varargin{1}(1);
         b = varargin{1}(2);
         u = [a*ones(1,d),...
                  a:(b-a)/(nmp-d):b,...
                     b*ones(1,d)]';
      case 'uniform'
         a = varargin{1}(1);
         b = varargin{1}(2);
         u = linspace(a,b,nmp+d)';
      case 'average'
         u = conv(t, [zeros(1,d) ones(1,d)]/d, 'same');
         u = [zeros(d,1);u(d+1:end-1);ones(d,1)];
   end
end
  1. 从t估计每个生成的参数值的 b 样条混合函数值。将结果存储在矩阵 B [nxm] 中。

构造 b 样条基矩阵B的函数:

% de-boor algorithm
B   = zeros(nmt_i,nmp); 
for j=0:nmp-1
   [b,px] = npa_contfit_bsplinebasis(j,d,u,t_i);
    B(:,j+1) = b;
end

估计 b 样条混合函数值的函数

function [y,t] = npa_contfit_bsplinebasis(i,d,u,t)

   y = bspline_basis_recur(i,d,u,t);

   % nested function
   function y = bspline_basis_recur(i,d,u,t)
      y = zeros(size(t));
      if d > 1
          b = npa_contfit_bsplinebasis(i,d-1,u,t);
          dn = t - u(i+1);
          dd = u(i+d) - u(i+1);
          if dd ~= 0  % indeterminate forms 0/0 are deemed to be zero
              y = y + b.*(dn./dd);
          end
          b = npa_contfit_bsplinebasis(i+1,d-1,u,t);
          dn = u(i+d+1) - t;
          dd = u(i+d+1) - u(i+1+1);
          if dd ~= 0
              y = y + b.*(dn./dd);
          end
      elseif u(i+2) < u(end)  % treat last element of knot vector as a special case
         y(u(i+1) <= t & t < u(i+2)) = 1;
      else
         y(u(i+1) <= t) = 1;
      end
   end % end
end
  1. 在P中估计并存储控制点向量p(在最小二乘意义上)

P = (pinv(B)*C')';

  1. 在封闭 b 样条的情况下,将控制点向量包裹起来

P = [P(:,end-d+1:end),P(:,d+1:end-d)]

开放 b 样条设置的结果:

http://imgur.com/7IFJdZ0

闭合(周期性)b 样条设置的结果

http://imgur.com/qfFFGjF

4

0 回答 0