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;

[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).
µ = (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 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 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 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 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:
      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