B-Spline Curves coefficients - division by zero (code in DELPHI)

Question

I was trying to implement the following recursive formula to my code

but to my surprise it turns out that after implementing this to DELPHI, I get an error due to division by zero. I am 98% sure that my knot vector is correctly calculated, which in a way means there shouldn't be any divisions by zero. I am 70% sure that the recursive formula is correctly implemented, for that reason I am posting my code here:

program project1;

uses
SysUtils;

Type
  TRealPoint = record
    x: single;
    y: single;
  end;

type
  TSample = Class(TObject)
    public
      KnotVector: array of single;
      FitPoints: array of TRealPoint;
      Degree: integer;
      constructor Create; overload;
      function Coefficient(i, p: integer; Knot: single): single;
      procedure GetKnots;
      destructor Destroy; overload;
  end;

constructor TSample.Create;
begin
  inherited;
end;

function TSample.Coefficient(i, p: integer; Knot: single): single;
var
  s1, s2: single;
begin
   If (p = 0) then
   begin
     If (KnotVector[i] <= Knot) And (Knot < KnotVector[i+1]) then Result := 1.0
     else Result := 0.0;
   end
   else
   begin
     s1 := (Knot - KnotVector[i])*Coefficient(i, p-1, Knot)/(KnotVector[i+p] - KnotVector[i]); //THIS LINE ERRORS due to division by zero ???
     s2 := (KnotVector[i+p+1]-Knot)*Coefficient(i+1,p-1,Knot)/(KnotVector[i+p+1]-KnotVector[i+1]);
     Result := s1 + s2;
   end;
end;

procedure TSample.GetKnots();
var
  KnotValue: single;
  i, MaxKnot: integer;
begin
  // KNOTS
  KnotValue:= 0.0;
  SetLength(KnotVector, Length(FitPoints) + 1 + Degree);
  MaxKnot:= Length(KnotVector) - (2*Degree + 1);
  for i := Low(KnotVector) to High(KnotVector) do
  begin
    if i <= (Degree) then KnotVector[i] := KnotValue / MaxKnot
    else if i > Length(FitPoints) then KnotVector[i] := KnotValue / MaxKnot
    else
    begin
      KnotValue := KnotValue + 1.0;
      KnotVector[i] := KnotValue / MaxKnot;
    end;
  end;
end;

destructor TSample.Destroy;
begin
   inherited;
end;

var
  i, j: integer;
  Test: TSample;
  N: array of array of single;
begin
  Test := TSample.Create;
  //define degree
  Test.Degree := 3;
  //random fit points
  j := 15;
  SetLength(Test.FitPoints, j + 1 + Test.Degree);
  For i := Low(Test.FitPoints) to High(Test.FitPoints) do
  begin
    Test.FitPoints[i].x := Random()*2000;
    Test.FitPoints[i].y := Random()*2000;
  end;
  //get knot vector
  Test.GetKnots;
  //get coefficients
  SetLength(N, j+1, j+1);
  For j := Low(N) to High(N) do
  begin
    For i := Low(N[j]) to High(N[j]) do
      begin
        N[j, i] := Test.Coefficient(i,3,Test.KnotVector[j]);
        write(floattostrf(N[j,i], ffFixed, 2, 2) + ', ');
      end;
    writeln();
  end;
  readln();
  Test.Free;
end.

Basically I'm not sure how to continue. I would need the values of matrix N (see this link) of basis coefficients but somehow using the formula from this link leads me to division by zero.

So... Is there a totally different way how to calculate those coefficients or what is the problem here?

UPDATE

Instead of using my own idea i tried to implement the algorithm from here as suggested by Dsm in the comments. As a result, there is no more divison by zero, but the result is totally unexpected anyways.

For n + 1 = 10 random fit points with spline degree 3 the basis matrix N (see link) is singular - as seen from the attached image.

Instead of that I would expect the matrix to be band matrix. Anyway, here is my updated code:

program project1;

uses
SysUtils;

Type
  TRealPoint = record
    x: single;
    y: single;
  end;

type
  TMatrix = array of array of double;

type
  TSample = Class(TObject)
    public
      KnotVector: array of double;
      FitPoints: array of TRealPoint;
      SplineDegree: integer;
      Temp: array of double;
      A: TMatrix;
      procedure GetKnots;
      function GetBasis(Parameter: double): boolean;
      procedure FormBasisMatrix;
  end;

procedure TSample.GetKnots();
var
  i, j: integer;
begin
  // KNOTS
  //https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/INT-APP/PARA-knot-generation.html
  SetLength(KnotVector, Length(FitPoints) + SplineDegree + 1);
  for i := Low(KnotVector) to High(KnotVector) do
  begin
    if i <= SplineDegree then KnotVector[i] := 0
    else if i <= (High(KnotVector) - SplineDegree - 1) then KnotVector[i] := (i - SplineDegree) / (Length(FitPoints) - SplineDegree)
    else KnotVector[i] := 1;
  end;
end;

function TSample.GetBasis(Parameter: double): boolean;
var
  m, d, k: integer;
  FirstTerm, SecondTerm: double;
begin
  //http://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-curve-coef.html
  Result := False;
  //initialize to 0
  SetLength(Temp, Length(FitPoints));
  For m := Low(Temp) to High(Temp) do Temp[m] := 0.0;
  //special cases
  If Abs(Parameter - KnotVector[0]) < 1e-8 then
  begin
    Temp[0] := 1;
  end
  else if Abs(Parameter - KnotVector[High(KnotVector)]) < 1e-8 then
  begin
    Temp[High(Temp)] := 1;
  end
  else
  begin
    //find knot span [u_k, u_{k+1})
    for k := Low(KnotVector) to High(KnotVector) do if Abs(KnotVector[k] - Parameter) < 1e-8 then break;
    Temp[k] := 1.0;
    for d := 1 to SplineDegree do
    begin
      Temp[k - d] := (KnotVector[k + 1] - Parameter) * Temp[k - d + 1] / (KnotVector[k + 1] - KnotVector[k - d + 1]);
      for m := k - d + 1 to k - 1 do
      begin
        FirstTerm := (Parameter - KnotVector[m]) / (KnotVector[m + d] - KnotVector[m]);
        SecondTerm := (KnotVector[m + d + 1] - Parameter) / (KnotVector[m + d + 1] - KnotVector[m + 1]);
        Temp[m] := FirstTerm * Temp[m] + SecondTerm * Temp[m + 1];
      end;
      Temp[k] := (Parameter - KnotVector[k]) * Temp[k] / (KnotVector[k + d] - KnotVector[k]);
    end;
  end;
  Result := True;
end;

procedure TSample.FormBasisMatrix;
var
  i, j: integer;
begin
  SetLength(A, Length(FitPoints), Length(FitPoints));
  for j := Low(A) to High(A) do
  begin
    for i := low(A[j]) to High(A[j]) do //j - row, i - column
    begin
      If GetBasis(KnotVector[j + SplineDegree]) then A[j, i] := Temp[i];
    end;
  end;
end;


var
  i, j, iFitPoints: integer;
  Test: TSample;
  N: array of array of single;
begin
  Test := TSample.Create;
  //define degree
  Test.SplineDegree := 3;
  //random fit points
  iFitPoints := 10;
  SetLength(Test.FitPoints, iFitPoints);
  For i := Low(Test.FitPoints) to High(Test.FitPoints) do
  begin
    Test.FitPoints[i].x := Random()*200;
    Test.FitPoints[i].y := Random()*200;
  end;
  //get knot vector
  Test.GetKnots;
  //get B-Spline basis matrix
  Test.FormBasisMatrix;
  // print matrix
  for j := Low(Test.A) to High(Test.A) do
  begin
    for i := Low(Test.A) to High(Test.A) do write(FloatToStrF(Test.A[j, i], ffFixed, 2, 2) + ', ');
    writeln();
  end;
  readln();
  Test.Free;
end.

Answer 1

This does not appear to be the complete answer, but it may help you on your way, and the result is closer to what you expect, but as I say, not completely there.

First of all the knots do not look right to me. The knots appear to form a 'ramp' function (clamped line), and though I can't work out if 'm' has any specific value, I would expect the function to be continuous, which yours is not. Making it continuous gives better results, e.g.

procedure TSample.GetKnots();
var
  i, j: integer;
  iL : integer;
begin
  // KNOTS
  //https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/INT-APP/PARA-knot-generation.html
  iL := Length( FitPoints );
  SetLength(KnotVector, iL + SplineDegree + 1);
  // set outer knot values and sum used to geterate first internal value
  for i := 0 to SplineDegree - 1 do
  begin
    KnotVector[ i ] := 0;
    KnotVector[ High(KnotVector)-i] := 1;
  end;
  // and internal ones
  for i := 0 to High(KnotVector) - 2* SplineDegree + 1 do
  begin
    KnotVector[ SplineDegree + i - 1] := i / (iL - 1);
  end;
end;

I introduced iL = Length( Fitpoints ) for convenience - it is not important.

The second issue I spotted is more of a programming one. In the GetBasis routine, you evaluate k by breaking a for loop. The problem with that is that k is not guaranteed to persist outside the loop, so your use of it later is not guaranteed to succeed (although it may)

Finally, in the same place, your range determination is completely wrong in my opinion. You should be looking for parameter to lie in a half open line segment, but instead you are looking for it to lie close to an endpoint of that line.

Putting these two together

   for k := Low(KnotVector) to High(KnotVector) do if Abs(KnotVector[k] - Parameter) < 1e-8 then break;

should be replaced by

k1 := 0;
for k1 := High(KnotVector) downto Low(KnotVector)  do
begin
  if Parameter >= KnotVector[k1] then
  begin
   k := k1;
   break;
  end;
end;

where k1 is an integer.

I can't help feeling that there is a plus 1 error somewhere, but I can't spot it.

Anyway, I hope that this helps you get a bit further.

Answer 2

To build recursive pyramid for coefficient calculation at intervals, you have to start top level of recursion (inner loop of calculations) from the first real (not duplicate) knot index:

 For i := Test.Degree...

Also check the last loop index.

P.S. You can remove constructor and destructor from class description and implementation if they have nothing but inherited.