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) of basis coefficients but somehow using the formula from this link 导致我被零除。

所以...有没有完全不同的方法来计算这些系数或者这里有什么问题？

更新

我没有使用我自己的想法，而是尝试在评论中实现 here as suggested by 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

要构建递归金字塔以间隔计算系数，您必须从第一个真实（非重复）节点索引开始递归的顶层（计算的内循环）：

 For i := Test.Degree...

还要检查最后一个循环索引。

P.S。如果 constructor 和 destructor 只有 inherited.

，您可以从 class 描述和实现中删除它们

Answer 2

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

首先我觉得结不对。结似乎形成了一个 'ramp' 函数（夹线），虽然我无法计算出 '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 不能保证在循环外持续存在，所以你以后使用它不能保证成功（尽管它可能）

最后，同一个地方，你的范围判断在我看来是完全错误的。您应该寻找位于半开放线段中的参数，但您正在寻找靠近该线端点的参数。

把这两个放在一起[=13=]

   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的错误，但就是看不出来

无论如何，我希望这能帮助你更进一步。

B 样条曲线系数 - 除以零（DELPHI 中的代码）

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

delphi

bspline

cubic-bezier