Find control points from B-Spline curve through set of data points in matlab

Question

I'm using data set with 200 data points that is used to draw B-Spline curve and I want to extract the 100 original control points from this curve to use it in one algorithm to solve one problem. The result of control points it's too small compared with the value of the data points of B-Spline curve so I don't know if I make something wrong in the following code or not I need help to know that because I must used these control points to complete my study in one algorithm

link of set of data points: https://drive.google.com/open?id=0B_2BUqaJptbqUkRWLWdmbmpQakk

Code :

% read data set

dataset = importdata("path of data set here"); 
x = dataset(:,1);  
y = dataset(:,2);

for i=1:200
controlpoints(i,1) = x(i);
controlpoints(i,2) = y(i);
controlpoints(i,3) = 0;


end

% Create Q with some points from originla matrix controlpoints ( I take only 103 points)
counter =1;
for i=1:200
   if (i==11) || (i==20) || (i==198)
      Q(counter,:) = F(i,:);
      counter = counter +1;
   end
   if ne(mod(i,2),0)
    Q(counter,:) = F(i,:);
    counter = counter+1;
   end

  end

I used Centripetal method to find control points from curve like the following picture

Complete my code:

% 2- Create Centripetal Nodes array from Q
   CP(1) = 0;
   CP(103) =1;

     for i=2:102
       sum = 0;
       for j=2:102
         sum = sum + sqrt(sqrt((Q(j,1)-Q(j-1,1))^2+(Q(j,2)-Q(j-1,2))^2));
       end
     CP(i) = CP(i-1) + (sqrt(sqrt((Q(i,1)-Q(i-1,1))^2+(Q(i,2)-Q(i-1,2))^2))/sum);
   end

   p=3;  % degree 

   % 3- Create U_K  array from CP array    

     for i=1:103
       U_K(i) = CP(i);
     end

To calculate control points we must follow this equation P=Qx(R') --> R' is inverse of R matrix so we must find R matrix then fins P(control points matrix) by the above equation. The following scenario is used to find R matrix

and to calculate N in B-Spline we must use these recursive function

Complete my code :

% 5- Calculate R_i_p matrix
  for a=1:100
    for b=1:100

       R_i_p(a,b) = NCalculate(b,p,U_K(a),U_K);
    end
  end

% 6- Find inverse of R_i_p matrix
  R_i_p_invers = inv(R_i_p);

% 7- Find Control points ( 100 points because we have curve with 3 degree  )


 for i=1:100
     for k=1:100
       PX(i) = R_inv(i,k) * Q(k,1);
       PY(i) = R_inv(i,k) * Q(k,2); 
     end
 end

 PX2 = transpose(PX);
 PY2 = transpose(PY);
 P = horzcat(PX2,PY2);  % The final control points you can see the values is very small compared with the original data points vlaues

My Recursive Function to find the previous R matrix:

function z = NCalculate(j,k,u,U)
    if (k == 1 )
        if ( (u > U(j)) && (u <= U(j+1)) )
           z = 1;
        else
           z = 0;
        end
    else
           z = (u-U(j)/U(j+k-1)-U(j)* NCalculate(j,k-1,u,U) ) + (U(j+k)-u/U(j+k)-U(j+1) *  NCalculate(j+1,k-1,u,U));
    end
  end

Really I need to this help so much , I tried in this problem from one week :(

Updated: Figure 1 for the main B-spline Curve , Figure 2 for the result control points after applied reverse engineering on this curve so the value is so far and so small compared with the original data points value

Updated(2): I made some updates on my code but the problem now in the inverse of R matrix it gives me infinite value at all time

 % 2- Create Centripetal Nodes array from Q
  CP(1) = 0;
  CP(100) =1;

  sum = 0;

    for i=2:100
      sum = sum + sqrt(sqrt((Q(i,1)-Q(i-1,1))^2+(Q(i,2)-Q(i-1,2))^2));
    end

    for i=2:99
        CP(i) = CP(i-1) + (sqrt(sqrt((Q(i,1)-Q(i-1,1))^2+(Q(i,2)-Q(i-1,2))^2))/sum);
    end

   % 3- Create U_K  array from CP array  

       for i=1:100
         U_K(i) = CP(i);
       end

  p=3;

  % create Knot vector
  % The first elements
     for i=1:p+1
         U(i) = 0;
     end

  % The last elements
     for i=100:99+p+1
         U(i) = 1;
     end

 % The remain elements
      for j=2:96
         sum = 0;
         for i=j:(j+p-1)
           sum = sum + U_K(i);
         end
         U(j+p) = (1/p)* sum;
      end

      % 5- Calculate R_i_p matrix
         for a=1:100
             for b=1:100

                  R_i_p(a,b) = NCalculate(b,p,U_K(a),U);
             end
         end

       R_i_p_invers = inv(R_i_p);

       % 7- Find Control points ( 100 points )


         for i=1:100
            for k=1:100
               if isinf(R_inv(i,k))
                  R_inv(i,k) = 0;
               end

              PX(i) = R_inv(i,k) * Q(k,1);
              PY(i) = R_inv(i,k) * Q(k,2); 
            end

           end


      PX2 = transpose(PX);
      PY2 = transpose(PY);
      P = horzcat(PX2,PY2);

Note: I updated my NCalculate recursive function to give me 0 if the result is NaN (not a number )

function z = NCalculate(j,k,u,U)
if (k == 1 )
    if ( (u >= U(j)) && (u < U(j+1)) )
        z = 1;
    else
        z = 0;
    end
else
    z = (u-U(j)/U(j+k-1)-U(j)* NCalculate(j,k-1,u,U) ) + (U(j+k)-u/U(j+k)-U(j+1) *  NCalculate(j+1,k-1,u,U));
end

if isnan(z)
    z =0;
end

end

Answer 1

I think there are a few dubious issues in your approach:

• First of all, if you try to create a b-spline curve interpolating 103 input points (and no other boundary conditions are imposed), the b-spline curve will have 103 control points regardless what degree the b-spline curve is.
• The U_K array is the parameter assigned to each input point. They are not the same as the knot sequence ti used by the Cox DeBoor recursive formula. If the b-spline curve is of degree 3, you shall have (103+3+1) knot values in the knot sequence. You can create the knot values in the following way:

0) Denote the parameters as p[i], where i = 0 to (n-1), p[0]=0.0 and n is number of points.
1) Create the knot values as

knot[0] = (p[1]+p[2]+p[3])/D  (where D is degree)
knot[1] = (p[2]+p[3]+p[4])/D
knot[2] = (p[3]+p[4]+p[5])/D 
...... 

These are the interior knot values. You should notice that p[0] and p[n-1] will not be used in this step. You will have (n-D-1) interior knots.

2) Now, add p[0] to the front of the knot values (D+1) times and add p[n-1] to the end of the knot values (D+1) times and you are done. At the end, you will have (N+D+1) knots in total.