Fitting closed b-spline to 2D data

Question

I am trying to implement my own code in Matlab to fit a closed b-spline to a set of 2D data.

I have generated the 2D data in form of an ellipse and used De boor algorithm. The code does work fine when set up for open b-spline, but it got stuck with closed (periodic) b-spline.

Here are the steps:

1. Generate the ordered 2D data C [2xn] in form of a 2D ellipse
2. Generate the parameter vector t [nx1] as per the uniform spaced or the chord-length method.
3. Generate the knot vector u [m+dx1] (:# control points, d: b-spline order): a) Open b-spline: uniform-clamped method b) closed b-spline: uniform method

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

4. Estimate the b-spline blending function values for every generated parameter value from t. Store the result in matrix B [nxm].

Function to construct the b-spline basis matrix 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

Function to estimate the b-spline blending function values

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

5. Estimate and store the control point vectors p (in least square sense) in P

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

6. In case of closed b-spline wrap the the control point vectors

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

The result for open b-spline set up:

https://i.stack.imgur.com/lgvIl.jpg

The result for the closed (periodic) b-spline set up

https://i.stack.imgur.com/TBodq.jpg