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

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

skrat
1#
skrat Published in 2018-02-12 07:24:02Z
 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. 
Dsm
2#
 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.
 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.