Ist ein Punkt in einem Polygon

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