Home B-Spline Curves coefficients - division by zero (code in DELPHI)
Reply: 0

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

user6973
1#
user6973 Published in May 25, 2018, 9:07 am

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.
You need to login account before you can post.

About| Privacy statement| Terms of Service| Advertising| Contact us| Help| Sitemap|
Processed in 0.30874 second(s) , Gzip On .

© 2016 Powered by mzan.com design MATCHINFO