delphi - B 样条曲线系数 - 除以零(DELPHI 中的代码)

Question

我试图在我的代码中实现以下递归公式

但令我惊讶的是，在 DELPHI 中实现此操作后，由于被零除而出现错误。我 98% 确定我的结向量计算正确，这在某种程度上意味着不应该被零除。我 70% 确定递归公式已正确实现，因此我在这里发布我的代码:

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.

基本上我不知道如何继续。我需要基本系数的矩阵 N ( see this link ) 的值，但以某种方式使用 this link 中的公式让我除以零。

那么...是否有完全不同的方法来计算这些系数，或者这里有什么问题？

更新

我没有使用自己的想法，而是尝试实现 here 中的算法。根据Dsm的建议在评论中。结果，不再被零除，但结果无论如何都是完全出乎意料的。

对于 n + 1 = 10 个样条阶数为 3 的随机拟合点，基础矩阵 N(参见 link )是奇异的 - 如附图所示。

我希望矩阵是 band matrix 。无论如何，这是我更新的代码:

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

这似乎不是完整的答案，但它可能会对您有所帮助，并且结果更接近您的期望，但正如我所说，并不完全存在。

首先，我觉得这些结看起来不太对劲。这些结似乎形成了一个“斜坡”函数(夹紧线)，虽然我无法计算出“m”是否有任何特定值，但我希望该函数是连续的，而你的不是。使其连续可以提供更好的结果，例如

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;

为了方便起见，我引入了 iL = Length( Fitpoints ) - 这并不重要。

我发现的第二个问题更多的是编程问题。在 GetBasis 例程中，您可以通过中断 for 循环来计算 k。问题在于 k 不能保证在循环之外持续存在，因此不能保证您稍后使用它会成功(尽管可能会成功)

最后，在同一个地方，我认为你的范围确定是完全错误的。您应该寻找位于半开线段中的参数，但实际上您正在寻找靠近该线端点的参数。

把这两个放在一起

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

应替换为

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

其中 k1 是一个整数。

我不禁感觉某个地方有一个加1的错误，但我无法发现它。

无论如何，我希望这可以帮助您走得更远。