Ist ein Punkt in einem Polygon

[Amateurprofi]

Nicht immer, aber immer wieder taucht die Frage auf, wie man wohl feststellen kann, ob ein Punkt innerhalb eines Dreiecks oder innerhalb einer anderen Figur liegt.

Ich habe dieses Thema mal ganz allgemein betrachtet und eine kleine Funktion geschrieben, die prüft, ob sich ein Punkt innerhalb eines konvexen nicht überschlagenden Polygons befindet.
Ein Programm zum Testen ist im Anhang.
Vielleicht interessiert es den einen oder anderen.

Delphi-Quellcode:

{------------------------------------------------------------------------------}
{ PtInPolygon                                                                  }
{ Prüft, ob ein Punkt innerhalb eines Polygons liegt.                          }
{ Das Polygon muß konvex sein und darf keine Überschneidungen haben.           }
{------------------------------------------------------------------------------}
FUNCTION PtInPolygon(const points:array of TPoint; N:integer; const p:TPoint):boolean;
var i:integer; angles:single;
//------------------------------------------------------------------------------
FUNCTION Angle(const p1,p2:TPoint):Extended;
var a2,b2,c2,cosa:single;
//--------------------------------------------------------------
begin
   a2:=Sqr(p2.x-p1.x)+Sqr(p2.y-p1.y);
   b2:=Sqr(p2.x-p.x)+Sqr(p2.y-p.y);
   c2:=Sqr(p1.x-p.x)+Sqr(p1.y-p.y);
   cosa:=(b2+c2-a2)/(2*Sqrt(b2)*Sqrt(c2));
   if cosa<=-1 then result:=pi
      else if cosa>=1 then result:=0
         else Result:=ArcCos(cosa);
end;
//------------------------------------------------------------------------------
begin
   result:=(n>2) and (n-1<=High(points));
   if not result then exit;
   for i:=0 to n-1 do if (p.x=points[i].x) and (p.y=points[i].y) then exit;
   angles:=Angle(points[n-1],points[0]);
   for i:=0 to n-2 do angles:=angles+Angle(points[i],points[i+1]);
   result:=Abs(angles-Pi*2)<0.00001;
end;

[DP-Maintenance]

Dieses Thema wurde von "SirThornberry" von "Programmieren allgemein" nach "Open-Source" verschoben.

[EDatabaseError]

In dem angehängten Polygon erkennt er garkeinen Punkt der innen liegt.

Tobi

[Neutral General]

Zitat:

Zitat von Amateurprofi

Ich habe dieses Thema mal ganz allgemein betrachtet und eine kleine Funktion geschrieben, die prüft, ob sich ein Punkt innerhalb eines konvexen nicht überschlagenden Polygons befindet.

Kann es sein das du da was überlesen hast EDatabaseError ;)

Gruß
Neutral General

[Amateurprofi]

Zitat:

Zitat von EDatabaseError

In dem angehängten Polygon erkennt er garkeinen Punkt der innen liegt.

Tobi

Also funktioniert das ganz ausgezeichnet.

Bei "Deiner" Figur handelt es sich um ein nichtkonvexes, überschlagendes Polygon, also um genau das Gegenteil eines konvexexn, nichtüberschlagenden Polygons. Im Begleittext steht aber unmißverständlich, daß die Prüffunktion nur für konvexe, nichtüberschlagende Polygone gedacht ist.

[JasonDX]

Geil :D sehr coole Idee, wie dus geloest hast :thumb:
Nur eine ganz kleine Kleinigkeit zur Polygon-Bedingung: Ein konvexes Polygon kann sich IMHO nicht ueberschlagen, also koenntest du das theoretisch weglassen ;)

greetz
Mike

[Gandalfus]

passent zum Thema:
http://blubplayer.de/tipppointinpolygon.html

ein Programm zum testen hab ich auch noch:
http://www.blubplayer.de/PointinPolygon.zip


Wie genau ist ein konvexes Polygon definiert?

[Amateurprofi]

Zitat:

Zitat von Gandalfus

passent zum Thema:
http://blubplayer.de/tipppointinpolygon.html

Superschnell die Funktion.
Jedoch (wenn ich nicht grobe Fehler gemacht habe), leider nicht zuverlässig.

Die in anhängendem Bild rot markierten Punkte wurden nicht richtig erkannt.

Teilweise wird für Punkte außerhalb des Polygons True zurückgegeben, teilweise für Punkte innerhalb des Polygons false.

[Amateurprofi]

Zitat:

Zitat von JasonDX

Geil :D sehr coole Idee, wie dus geloest hast :thumb:
Nur eine ganz kleine Kleinigkeit zur Polygon-Bedingung: Ein konvexes Polygon kann sich IMHO nicht ueberschlagen, also koenntest du das theoretisch weglassen ;)

greetz
Mike

Ja, dem könnte ich mich anschließen.
Aber sieh Dir mal anhängendes Bild an.
Die Angaben, welche Punkte zum Polygon gehören stehen daneben.

Weiß nicht wie Fachleute das sehen.....
Die äußere Begrenzung ist ein Quadrat, m.E. konvex.
Durch die Senkrechten innerhalb des Quadrates wird es zu einem überschlagenden Polygon.
Wir haben also ein konvexes überschlagendes Polygon.
Sicherlich sehr konstruiert....

[EDatabaseError]

Zitat:

Zitat von Neutral General

Zitat:

Zitat von Amateurprofi

Ich habe dieses Thema mal ganz allgemein betrachtet und eine kleine Funktion geschrieben, die prüft, ob sich ein Punkt innerhalb eines konvexen nicht überschlagenden Polygons befindet.

Kann es sein das du da was überlesen hast EDatabaseError ;)

Gruß
Neutral General

:oops: kann ja mal vorkommen...

[Gandalfus]

Zitat:

Teilweise wird für Punkte außerhalb des Polygons True zurückgegeben, teilweise für Punkte innerhalb des Polygons false.

Bis auf den Rand habe ich Fehler behoben.

Delphi-Quellcode:

function PointInPolygon(const APolygon: TPointArray; const TestPoint: TPoint): Boolean;

function getNextPoint(const APolygon: TPointArray; const index: integer): TPoint;
begin
  if index=High(APolygon) then
    result := APolygon[0]
  else
    result := APolygon[index+1];
end;

function getPreviousPoint(const APolygon: TPointArray; const index: integer): TPoint;
begin
  if index=0 then
    result := APolygon[High(APolygon)]
  else
    result := APolygon[index-1];
end;

var
  i: integer;
  currentPolygonPoint: TPoint;
  nextPolygonPoint: TPoint;
  previousPolygonPoint: TPoint;
  crossLines: integer;
begin
  crossLines:=0;
  result := false;
  for i := 0 to High(APolygon) do               // check polygon vertices
  begin
    currentPolygonPoint := APolygon[i];
    nextPolygonPoint := getNextPoint(APolygon,i);


    if (TestPoint.y=currentPolygonPoint.y) and (TestPoint.x=currentPolygonPoint.x) or
       //Horizontale Linie:
       ((TestPoint.y=currentPolygonPoint.y) and (TestPoint.y=nextPolygonPoint.y) and
        (TestPoint.x<currentPolygonPoint.x) and (TestPoint.x>nextPolygonPoint.x)) or
       //vertikale Linie:
       ((TestPoint.x=currentPolygonPoint.x) and (TestPoint.x=nextPolygonPoint.x)and
        (TestPoint.y>currentPolygonPoint.y) and (TestPoint.y<nextPolygonPoint.y)) then
    begin  //Fall 1  oder Punkt auf seknkrehcter oder horiuontaler Linie
      crossLines := 1;
      break;
    end
    else
    begin
      if (TestPoint.y=currentPolygonPoint.y) and (TestPoint.x<currentPolygonPoint.x) then
      begin //Fall 21
        previousPolygonPoint := getPreviousPoint(APolygon,i);
        if ((previousPolygonPoint.y>TestPoint.y) and (nextPolygonPoint.y<TestPoint.y)) or
           ((previousPolygonPoint.y<TestPoint.y) and (nextPolygonPoint.y>TestPoint.y)) then
        begin
          inc(crossLines);
        end;
      end
      else
      begin
        if ((TestPoint.y > currentPolygonPoint.y) and (TestPoint.y < nextPolygonPoint.y)) or
           ((TestPoint.y > nextPolygonPoint.y) and (TestPoint.y < currentPolygonPoint.y)) then
        begin

          if (TestPoint.x<currentPolygonPoint.x) and (TestPoint.x<nextPolygonPoint.x)then
          begin //Fall 3
            inc(crossLines);
          end;

          if ((TestPoint.x >= nextPolygonPoint.x) and (TestPoint.x < currentPolygonPoint.x)) or
             ((TestPoint.x >= currentPolygonPoint.x) and (TestPoint.x < nextPolygonPoint.x)) then
          begin //Fall 4
            if currentPolygonPoint.y < nextPolygonPoint.y then
            begin // Punkt1 ist oben links}
              if (nextPolygonPoint.x-currentPolygonPoint.x)/(nextPolygonPoint.y-currentPolygonPoint.y) >=
                 (TestPoint.x-currentPolygonPoint.x)/(TestPoint.y-currentPolygonPoint.y) then
                inc(crossLines);
            end;

            if currentPolygonPoint.y > nextPolygonPoint.y then
            begin // Punkt1 ist unten links
              if (nextPolygonPoint.x-currentPolygonPoint.x)/(nextPolygonPoint.y-currentPolygonPoint.y) <=
                 (TestPoint.x-currentPolygonPoint.x)/(TestPoint.y-currentPolygonPoint.y) then
                inc(crossLines);
            end;
          end;
        end;

      end;
    end;

  end;
  if (crossLines mod 2) <> 0 then result := true;

end;

[Amateurprofi]

Ein toller Algorithmus, und vor allem funktioniert er auch bei nichtkonvexen und sogar bei überschlagenden Polygonen.

[DP-Maintenance]

Dieses Thema wurde von "JasonDX" von "Open-Source" nach "Neuen Beitrag zur Code-Library hinzufügen" verschoben.
Als Codestueck/einzelne Funktion passt es besser in die CodeLib als nach OpenSource...

[LoCrux]

Fettes Problem.. Werd mich auch mal daran versuchen....

Bei Erfolg --> GEBE LAUT...

[LoCrux]

OK.. Die einfachste Lösung (und auch die schnellste) ist gefunden... :lol:

Vorausgesetz, daß sich NUR das Polygon vollständig sichtbar auf einem CANVAS OBJEKT befindet.....

Delphi-Quellcode:

FUNCTION PointInPolygon(const aCanvas:TCanvas;const TestPoint: TPoint):Boolean;
BEGIN
  RESULT := aCanvas.Pixles[TestPoint.x,TestPoint.y]<>HINTERGRUNDFARBE
END;

und der funktioniert für jedes X-Beliebige Polygon.

:mrgreen:

Hat man keinen Darstellung dann gehts auch mit einem temporären Bitmap......

[JasonDX]

Zitat:

Zitat von LoCrux

OK.. Die einfachste Lösung (und auch die schnellste) ist gefunden... :lol:

Vorausgesetz, daß sich NUR das Polygon vollständig sichtbar auf einem CANVAS OBJEKT befindet.....

Delphi-Quellcode:

FUNCTION PointInPolygon(const aCanvas:TCanvas;const TestPoint: TPoint):Boolean;
BEGIN
  RESULT := aCanvas.Pixles[TestPoint.x,TestPoint.y]<>HINTERGRUNDFARBE
END;

und der funktioniert für jedes X-Beliebige Polygon.

:mrgreen:

Jap, aber nur unter folgenden Bedingungen:

• Das Polygon ist auf einem Canvas (wie gesagt)
• Das Polygon ist mit einer einheitlichen, eindeutigen Farbe gefuellt
• Diese Farbe ist uns bekannt

Zitat:

Zitat von LoCrux

Hat man keinen Darstellung dann gehts auch mit einem temporären Bitmap......

Ja, das wars das aber mit der schnellsten Loesung, denn Canvas-malereien sind relativ langsam :zwinker:

Zudem ist die Loesung mathematisch betrachtet nicht zu gebrauchen: Ein Canvas arbeitet lediglich mit Pixeln und Integers, was zu Rundungsfehlern fuehrt. Damit kann es durchaus passieren, dass ein Wert, der nahe einer Kante liegt, falsch angezeigt wird. ;)

greetz
Mike

[LoCrux]

HI JasonDX,

war mir soweit alles klar..

Vielleicht gefällt Dir die Lösung besser..

Delphi-Quellcode:

PROCEDURE AddPolyPoint(VAR APoly:TPolyPoints;pnt:TPoint);
VAR
  pos : INTEGER;
BEGIN
  pos := HIGH(APoly)+1;
  SETLENGTH(APoly,pos+1);
  APoly[pos] := pnt;
  FPolyNew := YES;
END;

FUNCTION CopyPolyRange(APoly:TPolyPoints;startindex,endindex:INTEGER):TPolyPoints;
VAR
  cnt   : INTEGER;
  tPoly : TPolyPoints;
BEGIN
  IF (startindex>=0) AND (HIGH(APoly)>=endindex)
  THEN FOR cnt := startindex TO EndIndex DO AddPolyPoint(tPoly,APoly[cnt]);
  RESULT := tPoly;
END;

FUNCTION PointInTriangle(Tri:TPolyPoints;pnt:TPoint):BOOLEAN;
VAR
  a,b,c : DOUBLE;
BEGIN
 IF (HIGH(Tri)<3)
 THEN BEGIN
   // NOW THE POINT SHOULD BE INDEXED WITH 3 IN THE TRI ARRAY
   Self.AddPolyPoint(Tri,pnt);
   // FIRST PLACE A POINT TO ZERO NOW MOVE THE OTHERS RESPECTIVLY
   Tri[1] := SubPoints(Tri[1],Tri[0]);
   Tri[2] := SubPoints(Tri[2],Tri[0]);
   Tri[3] := SubPoints(Tri[3],Tri[0]);
   Tri[0] := SubPoints(Tri[0],Tri[0]);
   // COMPUTE
   a := (Tri[1].x*Tri[2].y)-(Tri[2].x*Tri[1].y);
   b := (Tri[1].x*Tri[3].y)-(Tri[3].x*Tri[1].y);
   c := (Tri[2].x*Tri[3].y)-(Tri[3].x*Tri[2].y);
   // RESULT
   RESULT := ((a*b)>0) AND ((a*c)<0) AND ((a*(a-b+c))>0);  // IF TRUE THEN POINT IN TRIANGLE
 END;
END;

FUNCTION ClosePoly(APoly:TPolyPoints):TPolyPoints;
BEGIN
  IF HIGH(APoly)>2
  THEN BEGIN
    RESULT := APoly;
    AddPolyPoint(RESULT,APoly[0]);
  END;
END;

FUNCTION PointInPoly_A1(APoly:TPolyPoints;Point:TPoint;Depth:INTEGER):BOOLEAN;
VAR
  cnt        : INTEGER;
  HighPoly   : TPolyPoints;
  LowPoly    : TPolyPoints;
  ePoly      : TPolyPoints;
  tPoly      : TPolyPoints;
BEGIN
  IF HIGH(APoly)>2
  THEN BEGIN
    ePoly    := ClosePoly(APoly);
    LowPoly  := CopyPolyRange(ePoly,LOW(ePoly)   ,HIGH(ePoly) DIV 2);
    HighPoly := CopyPolyRange(ePoly,HIGH(LowPoly),HIGH(ePoly));
    RESULT   := PointInPoly_A1(LowPoly,Point,Depth+1);
    IF NOT(RESULT) THEN RESULT := PointInPoly_A1(HighPoly,Point,Depth+1);
  END
  ELSE IF (HIGH(APoly)=2) THEN RESULT := PointInTriangle(APoly,Point);
END;

Hänge mal das komplette Testprojekt+Source mit an....

ABER IST QUICK'N'DIRTY... UND GILT AUCH NUR FÜR SICH NICHT ÜBERSCHNEIDENDE POLYGONE.. MUSS AN DER TRINAGULIERUNG NOCH ARBEITEN...

[Amateurprofi]

Das Manko aller bisher gezeigten Versionen ist, daß bei Punkten die sehr nahe an oder auf einer Linie liegen oft nicht richtig erkannt wird, ob ein Punkt nun in dem Polygon liegt oder nicht.

@Mike:
Ich meine, es kommt bei solchen Funktionen nicht darauf an "mathematisch korrekt" zu arbeiten.
Die Funktion soll ja nicht für theoretische Ermittlungen dienen, sondern sie soll feststellen, ob ein bestimmter Punkt der z.B. angeklickt wurde, innerhalb einer auf dem Bildschirm dargestellten Polygonst liegt oder nicht.

Ich habe den von Gandalfus in #12 vorgestellten Algorithmus aufgegriffen und "etwas" optimiert".

Ich behaupte nicht daß die unten gezeigte Funktion immer und ausnahmslos nur richtige Ergebnisse liefert, jedoch bei allen von mir getesteten Polygonen wurden alle Punkte korrekt erkannt. Das trifft auch für nichtkonvexe und auch für sich überschlagende Polygone zu.
Getestet wurde unter Windows XP Prof. (Ich weiß nicht ob andere Windows-Versionen andere Algorithmen zum zeichen von Linien verwenden.).

Zugegeben, etwas lang, aber doch recht flink....

Ein Programm zum Testen und Spielen hängt an.

Delphi-Quellcode:

{------------------------------------------------------------------------------}
{ PtInPolygon                                                                  }
{ Diese Version wurde aus einem Vorschlag von Henning Brackmann, in der        }
{ Delphi-Praxis auch als "Gandalfus" bekannt, weiterentwickelt.                }
{ Wie funktioniert das ?:                                                      }
{ Für alle Punktepaare p1=points[0..n-1], p2=points[1..n mod n] wird geprüft   }
{ - Ob p auf der Linie liegt. Dann liegt p im Polygon.                         }
{ - Bei wie vielen Punktepaaren eines der Folgenden zutrifft:                  }
{   - p liegt vertikal zwischen p1 und p2 und links von der durch p1,p2        }
{     gebildeten Linie. "Zwischen" heißt unterhalb p1 und oberhalb p2          }
{     oder umgekehrt.                                                          }
{   - p liegt links von und auf gleicher Höhe wie p1 und vertikal zwischen     }
{     p2 und dem Vorgänger von p1                                              }
{   Ist die Anzahl positiver Ergebnisse ungerade, dann liegt p im Polygon.     }
{------------------------------------------------------------------------------}
FUNCTION PtInPolygon(const points:array of TPoint; n:integer; const p:TPoint):Boolean;
type TPixelPos=(ppLeft,ppRight,ppOnLine);
var i,cl,left:integer; cp,np,pp:TPoint; ppos:TPixelPos;
//------------------------------------------------------------------------------
FUNCTION GetFactor(err,divisor,bound:integer):integer;
asm
      push  ebx
      mov   ebx,edx
      add   ebx,ebx
      xor   edx,edx
      div   ebx
      or    edx,edx
      jz    @1
      mov   ecx,1
@1:   add   eax,ecx
      pop   ebx
end;
//------------------------------------------------------------------------------
FUNCTION GetPosition(const p,cp,np:TPoint):TPixelPos;
var dx,dy,incx,incy,e,x,y,bound:integer;
begin
   // Prüft, ob p links von, auf oder rechts von der schrägen Linie cp --> np
   // liegt und gibt ppLeft,ppRight oder ppOnLine zurück.
   // Es wird vorausgesetzt, daß p sowohl vertikal wie auch horizontal
   // zwischen cp und np liegt.
   // Der Algorithmus ist von einem modifizierten Bresenham-Algorithmus
   // abgeleitet.
   // -----------------------------------------------------------------------
   // X und Y Differenzen ermitteln
   dx:=np.x-cp.x;
   dy:=np.y-cp.y;
   // Veränderungswerte für x und y festlegen
   if dx=0 then incx:=0
      else if dx>0 then incx:=1
         else begin
            incx:=-1;
            dx:=-dx;
         end;
   if dy=0 then incy:=0
      else if dy>0 then incy:=1
         else begin
            incy:=-1;
            dy:=-dy;
         end;
   // Schwellenwert festlegen
   if incx=incy then begin
      if incx>0 then bound:=0 else bound:=1;
   end else if incx>0 then begin
      if dx<=dy then bound:=0 else bound:=1;
   end else begin
      if dx>=dy then bound:=0 else bound:=1;
   end;
   // x Position des Pixels auf der Linie cp --> np  ermitteln, der auf der
   // gleichen y-Position liegt, wie p
   if dx<dy then begin
      // Dann wird y bei jedem Schritt um incy verändert. x dagegen wird
      // nur gelegentlich verändert
      // Die vertikale Distanz cp,p gibt direkt an, wie weit das Pixel, zu dem
      // die x-Position gesucht wird, auf der schrägen Linie von cp entfert ist.
      e:=Abs(p.y-cp.y)*(dx shl 1)-dy;
      if e<0 then x:=cp.x else x:=cp.x+incx*GetFactor(e,dy,bound);
      if p.x<x then result:=ppLeft
         else if p.x>x then result:=ppRight
            else result:=ppOnline;
   end else if dx=dy then begin
      // Dann ist das eine 45 Grad Linie. d.h. x und y werden bei jedem
      // Schritt incx bzw. incy verändert.
      x:=cp.x+Abs(p.y-cp.y)*incx;
      if p.x<x then result:=ppLeft
         else if p.x>x then result:=ppRight
            else result:=ppOnline;
   end else begin
      // Dann wird x bei jedem Schritt um incx verändert. y wird nur
      // gelegentlich verändert.
      // Anders als bei dx=dy und dx<dy ist in diesem Falle die vertikale
      // Distanz p-cp, nicht identisch mit der Entfernung des gesuchten Pixels
      // auf der schrägen Linie. Deshalb wird hier nicht aus der vertikalen
      // Distanz p-cp die horizontale Position des Pixels auf der Linie
      // ermittelt, sondern aus der horizontalen Distanz p-cp die
      // vertikale Position, aus der dann abgeleitet wird, ob p auf, links von
      // oder rechts von der Linie liegt
      e:=Abs(p.x-cp.x)*(dy shl 1)-dx;
      if e<0 then y:=cp.y else y:=cp.y+incy*GetFactor(e,dx,bound);
      if p.y=y then begin
         result:=ppOnline
      end else if incx=incy then begin // Linie läuft "\"
         if p.y<y then result:=ppRight else result:=ppLeft;
      end else begin // Linie läuft "/"
         if p.y<y then result:=ppLeft else result:=ppRight;
      end;
   end;
end;
//------------------------------------------------------------------------------
begin
   dec(n);
   result:=n<=High(points);
   if not result then exit;
   left:=maxint;
   cl:=0;
   for i := 0 to n do begin
      // current point und next point holen
      cp:=points[i];
      if i=n then np:=points[0] else np:=points[i+1];
      // Weitesten links liegenden Punkt definieren
      if cp.x<left then left:=cp.x;
      // Prüfen, ob p unter, über oder rechts von cp UND np liegt
      if (p.x>cp.x) and (p.x>np.x) then continue; // p liegt rechts von der Linie
      if (p.y>cp.y) and (p.y>np.y) then continue; // p liegt unterhalb der Linie
      if (p.y<cp.y) and (p.y<np.y) then continue; // p liegt oberhalb der Linie
      // Prüfen ob cp, np eine waagerechte, senkrechte oder schräge Linie bilden
      if cp.y=np.y then begin // cp, np bilden waagerechte linie
         if not ((p.x<cp.x) and (p.x<np.x)) then exit; // p liegt auf der Linie
         if cp.x<np.x then inc(cl);
      end else if cp.x=np.x then begin // cp, np bilden senkrechte linie
         if p.x=cp.x then exit; // p liegt auf der Linie
         if p.y<>np.y then inc(cl); // p liegt links von der Linie
      end else begin // cp, np bilden schräge linie
         if p.y=cp.y then begin
            if p.x<cp.x then begin
               if (p.x>np.x) and (Abs(cp.x-np.x)>Abs(cp.y-np.y)) then
                  if GetPosition(p,cp,np)=ppOnline then exit;
               if i=0 then pp:=points[n] else pp:=points[i-1];
               if ((p.y<pp.y) and (p.y>np.y)) or ((p.y>pp.y) and (p.y<np.y)) then inc(cl);
            end else if p.x>cp.x then begin
               if (p.x<np.x) and (Abs(cp.x-np.x)>Abs(cp.y-np.y)) then
                  if GetPosition(p,cp,np)=ppOnLine then exit;
            end else begin
               exit; // p liegt auf dem Eckpunkt
            end;
         end else if (p.x<cp.x) and (p.x<np.x) then begin
            if (p.y<>np.y) then inc(cl)
         end else begin
            ppos:=GetPosition(p,cp,np);
            if ppos=ppOnline then exit  // p liegt auf der Linie
               else if ppos=ppLeft then
                  if p.y<>np.y then inc(cl); // p liegt links von Linie
         end;
      end;
   end;
   result:=(p.x>=left) and (Odd(cl));
end;

[LoCrux]

@Amateurprofi

Hi.. und nun die Schelte... wie konntest Du dieses Thema nur aufgreifen?

Ich bin süchtig geworden nach diesem Problem.. das ist echt ne harte Nuss...
Bin gerade dabei so sämtliche verfügbare Literatur durchzustöbern...
Aber ich knack das... auch für sich überschlagende Polys (hoffentlich)...

irgendwie, irgendwo, irgendwann [NENA]

Merci An Dich....

LoCrux

[Amateurprofi]

Zitat:

Zitat von LoCrux

Hi.. und nun die Schelte... wie konntest Du dieses Thema nur aufgreifen?
Aber ich knack das... auch für sich überschlagende Polys (hoffentlich)...
Merci An Dich....
LoCrux

Ja, so ging es mir.
Falls du es überlesen hast : Die obige Version arbeitet für alle Polygone, konvexe, nichtkonvexe, und für sich üerschlagende.
Und um es nicht in Vergessenheit geraten zu lassen: Die Basis dafür wurde von Gandalfus veröffentlicht.....

[DerManu]

Hi erstmal! (hab mich für diesen post angemeldet)
ich finde Gandalfus' lösung des polygonproblems wirklich sehr elegant und konnte mir eine geometrisch analytische umsetzung einfach nicht verkneifen (Gandalfus, evtl. kennen wir uns sogar, ich schrieb vor einigen jahren als "Crosskiller" im delphisource forum).

Die obig gezeigte implementierung sieht mir recht umständlich aus, meine arbeitet mit einfachen schnitten zweier geraden als zentrale elemente.
Die mathematische überlegung dahinter:
man hat die zwei geraden in parameterform gegeben

Delphi-Quellcode:

g: x = a + µ*r
h: x = b + m*v

Nun setze ich beide geraden gleich und erhalte (weil wir ja zweidimensional arbeiten) zwei gleichungen (eine für die x eine für die y koordinate). Nun löse ich die (II) gleichung nach m in abhängigkeit von µ auf. Dies setze ich in die (I) ein und nach einigen mathematischen umformungen erhält man schließlich ein eindeutiges µ in abhängigkeit vom aufpunkt der beiden geraden (a,b) und den richtungsvektoren der beiden geraden (r,v).

Delphi-Quellcode:

µ = (vx*ax - vy*bx - vx*ay + vx*by) / (vx*ry - vy*rx)

zu beachten ist natürlich, dass bei parallelen geraden der divisor 0 wird und der schnittpunkt somit gegen unendlich geht.
Wollte man nun noch den genauen schnittpunkt ermitteln, so müsste man einfach den errechneten skalar µ in die erste gleichung (x = a+µr) einsetzen und man erhielte den exakten schnittpunkt. Uns genügt aber zu prüfen, ob der schnittpunkt der geraden im bereich der strecke [a+0*r] --> [a+1*r] liegt (damit nicht ein schnitt auftritt, wo unsere polygonwand bereits "zuende" ist). Man prüft nun also ob µ im intervall [0;1] liegt. Ist dies der fall, haben wir einen gültigen schnitt und "merken" ihn uns (inc(crosscounter) im code).

Diese schnittprüfung wenden wir systematisch auf alle wände des polygons an, die wir vorher in aufpunkt b und richtungsvektor v zerlegt haben. als zweite gerade nehmen wir, nach gandalfus'schem lösungsprinzip ;), eine gerade mit unserem zu prüfenden punkt als aufpunkt und einem richtugsvektor, welcher nach rechts zeigt ([1|0]).
Bleibt nun ein problem: der algorithmus sieht nun die letztere gerade nicht als strahl nach rechts sondern - ihr habt's vermutet - als eine gerade, berechnet also auch schnitte "links" vom punkt. Also die ganze geschichte nochmal umgekehrt: wir schneiden so, dass wir den skalar m erhalten, welcher die lage des schnittpunktes auf unserem [1|0] richtungsvektor angibt. Ist m < 0 so heisst das, dass wir links von unserem betrachteten aufpunkt sind, also der schnittpunkt nicht gültig ist.
zusammenfassend: schnittpunkt, wenn m > 0 und µ element von [0;1].

und abschließend: wenn wir eine ungerade anzahl an schnittpunkten finden, wissen wir, dass der punkt innerhalb des polygons liegt.

Hier nun der code. Beachtet, dass ich auf meinem zettel vor mir und im code das obige m als lambda bezeichnet habe. Leider gibts kein lambda in ascii codes ;(

Delphi-Quellcode:

type
  TVector = record
    X,Y: single;
  end;

function Makevector(ix,iy:single): TVector;
begin
  result.x := ix; result.y := iy;
end;

function GetStraightIntersect(a1,a2,r1,r2: TVector): single;
var divisor: single;
begin
  divisor := r2.x*r1.y-r2.y*r1.x;
  if not iszero(divisor) then 
    result := (r2.y*a1.x-r2.y*a2.x-r2.x*a1.y+r2.x*a2.y)/divisor
  else 
    result := INFINITY;
end;

function PtInPolygon(p:TVector; poly: array of TVector): boolean;
var i, crosscounter: integer;
    pv,polyv: TVector;
    mu, lambda: single;
begin
  crosscounter := 0;
  pv.x := 1; pv.y := 0;
  for i := 0 to high(poly) do
  begin
    if i = high(poly) then
      polyv := makevector(poly[0].x-poly[i].x,poly[0].y-poly[i].y)
    else
      polyv := makevector(poly[i+1].x-poly[i].x,poly[i+1].y-poly[i].y);
    mu := GetStraightIntersect(poly[i],p,polyv,pv);
    lambda := GetStraightIntersect(p,poly[i],pv,polyv);
    if inrange(mu,0,1) and (lambda > 0) then
      inc(crosscounter);
  end;
  result := not (crosscounter mod 2 = 0);
end;

Hab den code nun an einigen dingern getestet und es funktioniert gut.
Wer langeweile hat kann ja mal die vorgestellten lösungen benchmarken und sehen ob die analytische methode hier wenn's hart auf hart kommt dann nicht doch zurückstecken muss. wird ja doch einiges multipliziert und dividiert...
Ganz fein wäre natürlich eine komplette übersetzung in ASM ;)

Mit freundlichen Grüßen
Manu

[Amateurprofi]

Zitat:

Zitat von DerManu

Hab den code nun an einigen dingern getestet und es funktioniert gut.
Wer langeweile hat kann ja mal die vorgestellten lösungen benchmarken und sehen ob die analytische methode hier wenn's hart auf hart kommt dann nicht doch zurückstecken muss. wird ja doch einiges multipliziert und dividiert...
Ganz fein wäre natürlich eine komplette übersetzung in ASM ;)

Mit freundlichen Grüßen
Manu

Manu,
auch eine schöne Lösung.
aber :
Polygondaten liegen i.d.R. nicht als Singles vor sondern als Integers.
Die Umstellung auf Singles sollte deshalb innerhalb der Funktion erfolgen.
und:
Scheint nicht richtig zu arbeiten.
Die im anhängenden Bild rot markierten Punkte werden falsch erkannt.

[DerManu]

Zitat:

Zitat von Amateurprofi

Polygondaten liegen i.d.R. nicht als Singles vor sondern als Integers.
Die Umstellung auf Singles sollte deshalb innerhalb der Funktion erfolgen.

Für prüfungen auf z.b. der 2D programmoberfläche mit mauscursor stimmt das zweifellos. Meine PtInPoly-prüfung jedoch findet zukünftig in einer "zoombaren" umgebung statt, in der u.a. punkte mit single-genauigkeit dargestellt/gesetzt/usw. werden können. Da wird dann eine so zu sagen mathematische genauigkeit gebraucht.

Zitat:

Zitat von Amateurprofi

und:
Scheint nicht richtig zu arbeiten.
Die im anhängenden Bild rot markierten Punkte werden falsch erkannt.

Jap, da hast du recht. Wenn man das so fies hinbiegt, dass der strahl genau auf einen poly-punkt schießt, geht's schief ;)
Ich denke eine relativ sinnvolle abhilfe wäre, den strahl per zufall in eine richtung zu schicken und zu prüfen, ob einer der polygonpunkte auf ihm liegt, also die geradengleichung erfüllbar ist. Falls ja, einfach einen neuen vektor zufällig generieren.

Für meine einsatzzwecke genügt das bisherige aber, da es praktisch wohl nie vorkommen wird, dass der benutzer es schafft mit der maus auf alle 7 stellen des singles exakt auf die höhe eines punktes des polygons zu kommen (hoffentlich ;)).

Vielen Dank für die Kritik,
Manu

[Amateurprofi]

Zitat:

Zitat von DerManu

Vielen Dank für die Kritik,
Manu

Das war nicht als Kritik gemeint - nur als Hinweis. Ich weiß mittlerweile, wie schwer die Probleme bei bestimmten Konstellationen in den Griff zu bekommen sind.

[mimi]

Der Letzte Beitrag ist zwar etwas her, aber ich denke das macht nix oder ?

Ich habe eine Frage, arbeiten den jetzt die Letzte "DerManu" Funktion genau ?

Ich möchte sie gerne in einem 2D Spiel was ich mit Canvas mache benutzen. Das habe ich mir so vorgestellt:
Alle Objekte die ich nutzte, haben ein Polygone-Array. Damit ich prüfen kann, ob ein bestimmter Punkt da drin liegt.

[DerDan]

Hallo Polygonprofis,


ich bin auf der Suche nach einer Methode auch für Polygone, die
aus Segmenten und Kreisbögen bestehen, zu ermitteln ob ein gegebener Punkt
innerhalb liegt.
Eine Aufteilung der Bögen in viele kleine Segmente wollt ich mir aber sparen,
da

1.) evtl. zu ungenau
2.) die Berechnung sicher länger dauert, wenn man die Polygone immer erst zerlegen muss.

hat jemand eine Idee?


mfg

DerDan

[mr.winkle]

Zu Amateurprofis Methode:
bei mir tritt eine AV an folgender Stelle auf:

Delphi-Quellcode:

      if (p.x>cp.x) and (p.x>np.x) then continue; // p liegt rechts von der Linie

woran kann das liegen?

[Medium]

Vom Prinzip her müsste doch folgendes eigentlich völlig ausreichen. Ich hab es bei mir mehrfach im Einsatz, und es geht mit konvexen, konkaven, überschlagenen und sogar mit Polygonen mit Loch (so lange man sich an die Konvention hält, dass das Loch-Polygon dem äusseren Polygon entgegen orientiert ist).

Delphi-Quellcode:

function PtInPoly(const points: array of TPoint; p: TPoint): Boolean;
var
  i     : Integer;
  angle : Double;

  function Angle2D(x1, y1, x2, y2: Double): Double;
  var
    dtheta, theta1, theta2: Double;
  begin
    theta1 := ArcTan2(y1, x1);
    theta2 := ArcTan2(y2, x2);
    dtheta := theta2 - theta1;
    while (dtheta > PI) do
       dtheta := dtheta - 2*pi;
    while (dtheta < -PI) do
       dtheta := dtheta + 2*pi;
    result := dtheta;
  end;

begin
  angle := 0;
  for i := 0 to High(points) do
  begin
    p1.X := points[i].X - p.X;
    p1.Y := points[i].Y - p.Y;
    p2.X := points[(i+1) mod Length(points)].X - p.X;
    p2.Y := points[(i+1) mod Length(points)].Y - p.Y;
    angle := angle + Angle2D(p1.X, p1.Y, p2.X, p2.Y);
  end;
  result := Abs(angle) < pi;
end;

Bislang hatte ich mit keinen Exceptions oder Polygonen zu tun, die nicht richtig erkannt wurden.

[Mongfice]

Moin,
ich weiß, das Thema ist schon ewig alt, aber da ich keine endgültige Lösung sehe, frag ich einfach nochmal weiter.

Ich arbeite im Moment an einem Projekt indem es um Polygon-Überdeckung geht.
Ich bekomme die Eckdaten der Polygone aus einer Datenbank, jedes Polygon wird mit 4 Ecken angegeben, die aber nicht notwendigerweise verschieden sein müssen.

Im Endeffekt sollen auf einer Zeichenfläche dann alle Bildpunkte entsprechend der Anzahl sich dort überschneidender Polygone gefärbt werden.

Ich hab jetzt eigentlich alle hier vorgestellten Lösungen durchprobiert und die beste Lösung für meine Zwecke scheint mir die Lösung von Amateurprofi zu sein, da sie die geringste Fehlerquote hat (im angehängten Beispiel 17/1422900). Leider liegen die Fehler wie man im angehängten Bild sieht denkbar ungünstig - und tritt bei praktisch allen meinen Polygonen auf. Die entstehenden Fehlerlinie liegt zwar mal oben und mal unten (jenachdem in welche Richtung das Polygon "gekippt" ist), ist aber fast immer vorhanden und richtet sich in ihrer Größe anscheinend nach dem Winkel.

Ich hab schon etwas rumprobiert, aber leider liegen die anderen Lösungen mit Fehlerquoten von ca. 200-230 Punkten doch zu weit daneben und auch eine Anpassung der Lösung von Amateuerprofi ist mir leider bisher nicht gelungen.

Hat jemand von euch da vielleicht eine gute Idee oder eine Lösung?

Da ich auch Polygone mit nur genau einem Punkt habe (haben könnte), bräuchte ich eigentlich eine gänzlich exakte Lösung - aber die will ich eigentlich gar nicht erwarten *gg*

Achja, falls jemand nachbauen möchte, die Koordinaten für den Screenshot waren (in dieser Reihenfolge) (199,295),(179,295),(196,499),(214,499).

Gruß
Mongfice