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

# 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.
Processed in 0.30874 second(s) , Gzip On .