Wie kann man viele Punkte schnell vergleichen?

[Medium]

Im 3D-Bereich wird zum "culling" (Ablehnen von Polygonen/Vertices) gerne das Separating Axis Theorem ausgenutzt. Gerade mit achsenparallelen Rechtecken kann man da einiges optimieren, und das kann extrem flott gezüchtet werden. (Wird immerhin zum Frustum-Culling teils Bild für Bild in Spielen gemacht, und sogar in 3D mit schiefen Ebenen.) Es gibt die Wikipediaseite auch auf Deutsch, ich fand die Englische jetzt nur weit aus anschaulicher.
Das dürfte mehr oder minder State of the Art sein, es sei denn es gibt für 2D eine noch flottere Klamotte die mir dann unbekannt wäre

[Uwe Raabe]

Wenn ich das recht verstehe, betrachtest du lediglich die X-Werte der Punkte mit den jeweiligen Fenstergrenzen. Die Y-Werte interessieren dabei noch gar nicht.

Die Punkte scheinen in X aufsteigend sortiert zu sein, sonst funktioniert dein Algorithmus nämlich nicht.

Wenn man nun die Fenstergrenzen ebenfalls sortiert und mit Index-Arrays arbeitet, kann man mit zwei Durchläufen die Punktbereiche in den Fenstern ermitteln.

Hier der passende Code:

Delphi-Quellcode:

{ Punkte gibt die X-Werte der Punkte an, FensterLo die Untergrenzen, FensterHi die Obergrenzen der Fenster.
  insideLo und insideHi geben die Indizes der Punkte an, die gerade noch in dem jeweiligen Fenster liegen }
procedure FindeFenster(const Punkte, FensterLo, FensterHi: array of Double; var insideLo, insideHi: TIntegerDynArray);
var
  CntPunkte: Integer;
  CntFenster: Integer;
  orderLo: TIntegerDynArray;
  orderHi: TIntegerDynArray;
  I: Integer;
  idx: Integer;
  rng: Integer;
begin
  CntPunkte := Length(Punkte);
  CntFenster := Length(FensterLo);

  Assert(CntPunkte > 0);
  Assert(CntFenster > 0);
  Assert(CntFenster = Length(FensterHi));

  SetLength(insideLo, CntFenster);
  SetLength(insideHi, CntFenster);

  { orderLo enthält eine Liste von Indizes auf FensterLo, so daß die Werte aufsteigend sortiert sind.
    orderHi enthält eine Liste von Indizes auf FensterHi, so daß die Werte aufsteigend sortiert sind.
    Läßt sich mit einem angepassten QuickSort leicht implementieren. }

  orderLo := GetSortOrder(FensterLo);
  orderHi := GetSortOrder(FensterHi);

  { erster Durchlauf: untere Grenzen ermitteln }
  idx := 0;
  rng := orderLo[idx];
  for I := 0 to Length(Punkte) - 1 do begin
    if Punkte[I] >= FensterLo[rng] then begin
      { untere Grenze des nächsten Fensters überschritten }
      insideLo[rng] := I;
      { die folgenden Fenster prüfen bis wir eins finden, in dem wir noch nicht sind }
      Inc(idx);
      while idx < CntFenster do begin
        rng := orderLo[idx];
        if Punkte[I] < FensterLo[rng] then
          Break;
        insideLo[rng] := I;
        Inc(idx);
      end;
      { schon alle Fenster überprüft? }
      if idx >= CntFenster then
        Break;
      rng := orderLo[idx]
    end;
  end;

  { zweiter Durchlauf: obere Grenzen ermitteln }
  idx := CntFenster - 1;
  rng := orderHi[idx];
  for I := CntPunkte - 1 downto 0 do begin
    if Punkte[I] <= FensterHi[rng] then begin
      { obere Grenze des nächsten Fensters unterschritten }
      insideHi[rng] := I;
      { die folgenden Fenster prüfen bis wir eins finden, in dem wir noch nicht sind }
      Dec(idx);
      while idx >= 0 do begin
        rng := orderHi[idx];
        if Punkte[I] > FensterHi[rng] then
          Break;
        insideHi[rng] := I;
        Dec(idx);
      end;
      { schon alle Fenster überprüft? }
      if idx < 0 then
        Break;
      rng := orderHi[idx];
    end;
  end;
end;

Die Funktion GetSortOrder ruft einen angepassten Quicksort auf:

Delphi-Quellcode:

procedure QuickSort(const Values: array of Double; var Index: array of Integer; L, R: Integer);
var
  I, J, T: Integer;
  P: Double;
begin
  repeat
    I := L;
    J := R;
    P := Values[Index[(L + R) shr 1]];
    repeat
      while (Values[Index[I]] < P) do
        Inc(I);
      while (Values[Index[J]] > P) do
        Dec(J);
      if I <= J then
      begin
        if I <> J then
        begin
          T := Index[I];
          Index[I] := Index[J];
          Index[J] := T;
        end;
        Inc(I);
        Dec(J);
      end;
    until I > J;
    if L < J then
      QuickSort(Values, Index, L, J);
    L := I;
  until I >= R;
end;

function GetSortOrder(const Values: array of Double): TIntegerDynArray;
var
  I: Integer;
begin
  SetLength(result, Length(Values));
  for I := 0 to Length(result) - 1 do
    result[I] := I;

  QuickSort(Values, result, 0, Length(result) - 1);
end;

[Medium]

Öhm, sortiert und nur für eine Achse ginge das doch noch VIEL einfacher! Die Fenstergrenzen müssen dafür auch nach gleicher Weise wie die Punkte sortiert sein, und du kommst mit nur einem Durchlauf aus.

Mal so ein Fetzen Pseudocode:

Delphi-Quellcode:

type
  TWindow = class;
  TEdge = class
  public
    Value: Float;
    IsLeftEdge: Boolean; // im Konstruktor zuweisen, wird noch wichtig :)
    Window: TWindow; // im Konstruktor zuweisen, quasi der Parent, also Rückbezug
  end;
  TWindow = class
  public
    Left, Right: TEdge;
    LeftPointLeft: Float;
    LeftPointRight: Float;
    RightPointLeft: Float;
    RightPointRight: Float;
  end;

var
  Windows: array of TWindow;
  Edges: array of TEdge; // referenzen auf alle Edges aus dem Windows-Array sortiert nach Value
  CurrentEdgeIndex: Integer;
begin
  for i := 0 to maxPoints do
  begin
    if Points[i].X > Edges[CurrentEdgeIndex].Value then
    begin
      if Edges[CurrentEdgeIndex].IsLeftEdge then
      begin
        Edges[CurrentEdgeIndex].Window.LeftPointLeft := Points[i-1].Value;
        Edges[CurrentEdgeIndex].Window.LeftPointRight := Points[i].Value;
      end
      else
      begin
        Edges[CurrentEdgeIndex].Window.RightPointLeft := Points[i-1].Value;
        Edges[CurrentEdgeIndex].Window.RightPointRight := Points[i].Value;
      end;
      inc(CurrentEdgeIndex);
    end;
  end;
end;

Und schwupps stehen im TWindows-Array alle Punkte, die links und rechts der X-Fenstergrenzen liegen, in O(n).

Man könnte natürlich auch die Punkte in die TEdge Instanzen werfen, und sich dieses IsLeftEdge damit sparen:

Delphi-Quellcode:

type
  TWindow = class;
  TEdge = class
  public
    Value: Float;
    LeftValue, RightValue: Float;
  end;
  TWindow = class
  public
    Left, Right: TEdge;
  end;

var
  Windows: array of TWindow;
  Edges: array of TEdge; // referenzen auf alle Edges aus dem Windows-Array sortiert nach Value
  CurrentEdgeIndex: Integer;
begin
  for i := 0 to maxPoints do
  begin
    if Points[i].X > Edges[CurrentEdgeIndex].Value then
    begin
      Edges[CurrentEdgeIndex].LeftValue := Points[i-1].Value;
      Edges[CurrentEdgeIndex].RightValue := Points[i].Value;
      inc(CurrentEdgeIndex);
    end;
  end;
end;

Das spart evtl. noch (sehr minimal) Zeit in der Schleife (ein Vergleich und ein paar Dereferenzierungen, uhu), und die Werte stehen ggf. an sinnvollerer Stelle, je nach dem wie du es nachher brauchst.

[Uwe Raabe]

@Medium: Es gibt da neben ein paar Grenzfällen, in denen der Code einen ERangeError auslöst (z.B. erster Punkt liegt innerhalb des ersten Fensters, letzter Punkt liegt außerhalb des letzten Fensters), noch das Problem, daß er fehlerhafte Resultate liefert. Mit den gegebenen Werten ergänzt um ein paar weitere Fenster mit Spezialfällen (drittes Fenster umschließt erstes, kompletter Bereich, genau ein Punkt, nur ein Punkt im Bereich) bekomme ich folgende Ergebnisse:

Code:

Punkte: 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 

erwartet: -------------------------------------------------------------
0.7-1.8: 1.0, 1.5
1.4-3.7: 1.5, 2.0, 2.5, 3.0, 3.5
0.6-2.1: 1.0, 1.5, 2.0
0.0-5.0: 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0
2.0-2.0: 2.0
2.9-3.1: 3.0

erhalten: -------------------------------------------------------------
0.7-1.8: 1.5 , 2.0 
1.4-3.7: 
0.6-2.1: 1.0 , 1.5 , 2.0 , 2.5 , 3.0 , 3.5 
0.0-5.0: 
2.0-2.0: 3.0 
2.9-3.1: 4.5

[himitsu]

Hab jetzt nur mal auf die Schnelle das Find der TStringList versucht anzupassen und hoffe es stimmt soweit auch.

Delphi-Quellcode:

type TMyArray = array of record  X, Y: Integer;  end;

function FindLowerXIndex(const Arr: TMyArray; LowX: Integer): Integer;
var
  L, H, I: Integer;
begin
  L := 0;
  H := High(Arr);
  while L <= H do
  begin
    I := (L + H) shr 1;
    if Arr[I].X < LowX then L := I + 1 else H := I - 1;
  end;
  Result := L;
end;

Das Array muß natürlich nach X sortiert sein.

[Medium]

Zitat von Uwe Raabe:

@Medium: Es gibt da neben ein paar Grenzfällen, in denen der Code einen ERangeError auslöst (z.B. erster Punkt liegt innerhalb des ersten Fensters, letzter Punkt liegt außerhalb des letzten Fensters)

Dessen war ich mir bewusst, war aber gestern Nacht zu faul explizit drauf einzugehen - ich nahm einfach mal an, dass das offensichtlich sein, dass wenn solche Fälle existieren, sie einer Spezialbehandlung bedürfen

Zitat:

noch das Problem, daß er fehlerhafte Resultate liefert. Mit den gegebenen Werten ergänzt um ein paar weitere Fenster mit Spezialfällen (drittes Fenster umschließt erstes, kompletter Bereich, genau ein Punkt, nur ein Punkt im Bereich) bekomme ich folgende Ergebnisse:

Wie gesagt, das hab ich durchaus erwartet. Sobald Fenster Edges teilen (bzw. identische haben), ein Fenster 0 Punkte umschließt oder Edges ausserhalb des Bereiches liegen, müsste ein wenig ge-if-then-else betrieben werden. Meine Intention war auch nicht, einen fix und fertigen Algo einzustellen, sondern nur fix die spontane Idee zu illustrieren

Schöne Grüße,
Medium