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;

[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.

[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...

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.



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

[JasonDX]

Zitat von LoCrux:

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

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.

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 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

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 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.....