Punkt in Polygon

Für allgemeine Fragen zur Programmierung mit PureBasic.
Benutzeravatar
dige
Beiträge: 1182
Registriert: 08.09.2004 08:53

Punkt in Polygon

Beitrag von dige »

Hallo Zusammen,

da wir hier ja einige Grafikcracks haben, hoffe ich, ihr könnt mir bei folgender Aufgabenstellung helfen:

Ich habe eine Liste mit Punkten (Geokoordinaten im WGS85 Format) und mehrere Flächen (Umrisse, Polygone)
und möchte prüfen, welcher Punkt sich in welchem Polygon befindet.

Wie würdet ihr da rangehen?

Ciao Dige
"Papa, mein Wecker funktioniert nicht! Der weckert immer zu früh."
Benutzeravatar
NicTheQuick
Ein Admin
Beiträge: 8675
Registriert: 29.08.2004 20:20
Computerausstattung: Ryzen 7 5800X, 32 GB DDR4-3200
Ubuntu 22.04.3 LTS
GeForce RTX 3080 Ti
Wohnort: Saarbrücken
Kontaktdaten:

Re: Punkt in Polygon

Beitrag von NicTheQuick »

Wir hatten sowas schon mal hier im Forum. Falls die Polygone konvex sind, wird es einfach. Falls nicht, müssen härtere Geschütze ran.

Edit: Vielleicht hilft dir ja schon das da: Polygon Include - Vereinigung,Schnitt,Differenz - Testphase
Bild
DarkDragon
Beiträge: 6267
Registriert: 29.08.2004 08:37
Computerausstattung: Hoffentlich bald keine mehr
Kontaktdaten:

Re: Punkt in Polygon

Beitrag von DarkDragon »

Da es hier um geografische Daten (Länderumrisse usw.) geht vermute ich, dass es konkave, zweidimensionale Polygone sind. Wenn man es nicht allzu genau braucht, würde ich die Polygone auf ein Bild malen, jedes mit einer anderen Farbe, dann prüfen welche Farbe unter dem Punkt ist. Dann ist die Genauigkeit von der Bildauflösung abhängig.

Wenn man es genau braucht verwendet man Trapezoidal Maps, die nichts weiter sind als eine Aufteilung des Raums aller Polygone in Trapeze, indem man die Vertizen vertikal erweitert bis zur umgebenden Bounding Box.
Angenommen es gäbe einen Algorithmus mit imaginärer Laufzeit O(i * n), dann gilt O((i * n)^2) = O(-1 * n^2) d.h. wenn man diesen Algorithmus verschachtelt ist er fertig, bevor er angefangen hat.
Benutzeravatar
dige
Beiträge: 1182
Registriert: 08.09.2004 08:53

Re: Punkt in Polygon

Beitrag von dige »

Danke NicTheQuick und DarkDragon für die Tipps!

Durch einen Beitrag von Infratec bin ich gerade auf die Funktion

Code: Alles auswählen

Ergebnis = IsInsidePath(x.d, y.d [, KoordinatenSystem])
aufmerksam geworden.

Wenn ich die Geokoordinaten alle in Pixelkoordinaten umrechnen würde,
wäre dann diese Funktion geeignet um so einen Test zu machen? Mit der
Vector Lib habe ich leider noch keine Erfahrung.

Vielleicht geht es auch ganz anders? Die Fläche (Polygon) hole ich mir von OpenStreetMap:
https://nominatim.openstreetmap.org/sea ... und%20GmbH

Ich will nun prüfen welche Postleitzahlen sich in dieser Fläche befinden.
Die PLZ hole ich mir hier: https://github.com/mdornseif/pyGeoDb/bl ... db_plz.txt

Ciao Dige
"Papa, mein Wecker funktioniert nicht! Der weckert immer zu früh."
Benutzeravatar
alter Mann
Beiträge: 201
Registriert: 29.08.2008 09:13
Wohnort: hinterm Mond

Re: Punkt in Polygon

Beitrag von alter Mann »

...ist zwar schon etwas älter...

viewtopic.php?f=8&t=18209
Win11 64Bit / PB 6.0
ccode_new
Beiträge: 1214
Registriert: 27.11.2016 18:13
Wohnort: Erzgebirge

Re: Punkt in Polygon

Beitrag von ccode_new »

Hallo!

siehe:
Link von Nic

->siehe:
;+ Tests the X- and Y-coordinates whether its inside, outside or on the border of a valid polygon.
Procedure PolygonPointTest(*Polygon.Polygon, X.d, Y.d)

; Quelle: https://de.wikipedia.org/wiki/Punkt-in- ... ach_Jordan
Aber ob das in die richtige Richtung führt ?

Sind die Polygone Ländern gleichzusetzen ? Sind die Polygone auf einem eigenen Bild (z.B. Landkarte) ? Haben diese Polygone eine einheitliche Abgrenzung (auch Farbe) ?
Betriebssysteme: div. Windows, Linux, Unix - Systeme

no Keyboard, press any key
no mouse, you need a cat
Michael Vogel
Beiträge: 71
Registriert: 16.03.2006 11:20

Re: Punkt in Polygon

Beitrag von Michael Vogel »

Der Thread ist vielleicht nicht ganz taufrisch, aber nachdem ich auf der Suche nach einer Funktion IsPointInsidePolygon bin, bin ich auch auf diesen Beitrag gestoßen:

https://stackoverflow.com/questions/275 ... in-polygon

Solange keine größeren Regionen der Erde untersucht werden, reicht die eukldische Geometrie - denn richtig rund läuft auf der Erde ohnehin nichts /:->


Übrigens, der folgende Code läuft auch noch nicht richtig rund...

Code: Alles auswählen


Global mousePosition.Point, mouseState, polygonPoints = 7
Global Dim polygonPoints.Point(10)


polygonPoints(0)\x = 100
polygonPoints(0)\y = 200
polygonPoints(1)\x = 350
polygonPoints(1)\y = 250
polygonPoints(2)\x = 200
polygonPoints(2)\y = 100
polygonPoints(3)\x = 550
polygonPoints(3)\y = 250
polygonPoints(4)\x = 400
polygonPoints(4)\y = 300
polygonPoints(5)\x = 350
polygonPoints(5)\y = 350
polygonPoints(6)\x = 300
polygonPoints(6)\y = 500
polygonPoints(7)\x = 250
polygonPoints(7)\y = 350


Procedure InsidePolygon(*point.point,dots,Array poly.point(1))

	; int pnpoly(int npol, float *xp, float *yp, float x, float y)
	; {
	;    int i, j, c = 0;
	;    For (i = 0, j = npol-1; i < npol; j = i++) {
	;      If ((((yp[i] <= y) && (y < yp[j])) ||
	;         ((yp[j] <= y) && (y < yp[i]))) &&
	;         (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
	;        c = !c;
	;    }
	;    Return c;
	; }

	Protected i,j,c

	While i<dots
		j=(i+dots-1)%dots
		Debug Str(i)+" "+Str(j)
			If ( ( (poly(i)\y <= *point\y) And (*point\y < poly(j)\y) ) Or
				( (poly(j)\y <= *point\y) And (*point\y < poly(i)\y) ) ) And
				( *point\x < (poly(j)\x-poly(i)\x) * (*point\y-poly(i)\y) / (poly(j)\y-poly(i)\y) + poly(i)\x)
				c!1
			EndIf
		i+1
	Wend

	ProcedureReturn c

EndProcedure

Procedure PointInPolygon(*polygon.Point, sides, *mousePosition.Point)

	Static filled
	result = #False
	

	If StartVectorDrawing(CanvasVectorOutput(0))
		
		VectorSourceColor($FFFFFFFF)
		FillVectorOutput()
		
		VectorSourceColor($FF000000)
		
		result=insidePolygon(*mousePosition,polygonPoints,polygonPoints())

		MovePathCursor(polygonPoints(0)\x,polygonPoints(0)\y)
		For i = 1 To polygonPoints-1
			AddPathLine(polygonPoints(i)\x,polygonPoints(i)\y)
		Next i
		ClosePath()
		
		If result
			FillPath()
		Else
			StrokePath(5)
		EndIf
		
		For i = 0 To polygonPoints-1
			MovePathCursor(polygonPoints(i)\x,polygonPoints(i)\y)
			DrawVectorText("("+Str(i)+") "+Str(polygonPoints(i)\x)+" | "+Str(polygonPoints(i)\y))
		Next i
		
		MovePathCursor(10,50)
		DrawVectorText(Str(*mousePosition\x)+" | "+Str(*mousePosition\y))

		FillPath()

		StopVectorDrawing()
		Delay(1)
	EndIf

	ProcedureReturn result

EndProcedure


polygonPoints(polygonPoints)=polygonPoints(0)


wFlags = #PB_Window_SystemMenu | #PB_Window_ScreenCentered
If OpenWindow(0, 0, 0, 600, 630, "PointInPolygon()", wFlags)
	CanvasGadget(0, 0, 30, 600, 600, #PB_Canvas_Container)
	TextGadget(1, 0, 0, 600, 30, "", #PB_Text_Center)

	Repeat
		event = WaitWindowEvent()
		Select event
		Case #PB_Event_CloseWindow
			appQuit = 1
		Case #PB_Event_Gadget
			If EventType() = #PB_EventType_MouseMove
				mousePosition\x = GetGadgetAttribute(0, #PB_Canvas_MouseX)
				mousePosition\y = GetGadgetAttribute(0, #PB_Canvas_MouseY)
				If mouseState <> PointInPolygon(@polygonPoints(), polygonPoints, @mousePosition)
					mouseState ! 1
					If mousestate
						txt$ = "Mouse is in..."
					Else
						txt$ = "Mouse is out..."
					EndIf
					SetGadgetText(1, Str(mousePosition\x)+" | "+Str(mousePosition\y)+" >> "+txt$)
				EndIf
			EndIf
		EndSelect
	Until appQuit
EndIf
ccode_new
Beiträge: 1214
Registriert: 27.11.2016 18:13
Wohnort: Erzgebirge

Re: Punkt in Polygon

Beitrag von ccode_new »

@Michael Vogel
Michael Vogel hat geschrieben:Übrigens, der folgende Code läuft auch noch nicht richtig rund...
Warum?

Das läuft doch schon ganz schön rund.

Code: Alles auswählen

Structure TPoint
  x.d
  y.d
EndStructure

Structure TCircle
  p.TPoint
  r.d
EndStructure


Procedure EinfacherKreis(*A.TPoint, *B.TPoint, *C.TPoint, *Kreis.TCircle)
  Protected.TPoint p1, p2
  Protected.d vB, vC, vD
  p1\x = *B\x - *A\x : p1\y = *B\y - *A\y
  p2\x = *C\x - *A\x : p2\y = *C\y - *A\y
  vB = p1\x * p1\x + p1\y * p1\y
  vC = p2\x * p2\x + p2\y * p2\y
  vD = p1\x * p2\y - p1\y * p2\x
  *Kreis\p\x = (p2\y * vB - p1\y * vC) / (2 * vD)
  *Kreis\p\y = (p1\x * vC - p2\x * vB) / (2 * vD)
  *Kreis\p\x + *A\x : *Kreis\p\y + *A\y
  *Kreis\r = Sqr(Pow(*Kreis\p\x - *A\x, 2) + Pow(*Kreis\p\y - *A\y, 2))
EndProcedure

Global mousePosition.TPoint, mouseState, polygonPoints = 360
Global Dim polygonPoints.TPoint(360)
Global.i Kreis.TCircle

polygonPoints(0)\x = 100
polygonPoints(0)\y = 200
polygonPoints(1)\x = 350
polygonPoints(1)\y = 250
polygonPoints(2)\x = 200
polygonPoints(2)\y = 100


Procedure InsidePolygon(*point.TPoint, dots, Array poly.TPoint(1))
  
  ; int pnpoly(int npol, float *xp, float *yp, float x, float y)
  ; {
  ;    int i, j, c = 0;
  ;    For (i = 0, j = npol-1; i < npol; j = i++) {
  ;      If ((((yp[i] <= y) && (y < yp[j])) ||
  ;         ((yp[j] <= y) && (y < yp[i]))) &&
  ;         (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
  ;        c = !c;
  ;    }
  ;    Return c;
  ; }
  
  Protected i=0,j,c
  
  While i<dots
    j=(i+dots-1)%dots
    ;Debug Str(i)+" "+Str(j)
    If ( ( (poly(i)\y <= *point\y) And (*point\y < poly(j)\y) ) Or
         ( (poly(j)\y <= *point\y) And (*point\y < poly(i)\y) ) ) And
       ( *point\x < (poly(j)\x-poly(i)\x) * (*point\y-poly(i)\y) / (poly(j)\y-poly(i)\y) + poly(i)\x)
      c!1
    EndIf
    i+1
  Wend
  
  ProcedureReturn c
  
EndProcedure

Procedure PointInPolygon(*polygon.TPoint, sides, *mousePosition.TPoint)
  result = #False
  
  If StartDrawing(CanvasOutput(1))
    
    result=insidePolygon(*mousePosition,polygonPoints,polygonPoints())
    
    For i = 0 To polygonPoints-1
      Plot(polygonPoints(i)\x, polygonPoints(i)\y, RGB(0,0,0))
    Next i
    
    If result
      FillArea(WindowWidth(0)/2, WindowHeight(0)/2, RGB(0,0,0), RGB(255,0,0))
    Else
      FillArea(WindowWidth(0)/2, WindowHeight(0)/2, RGB(255,255,255), RGB(255,255,255))
    EndIf
    
    DrawText(0, 0, Str(*mousePosition\x)+" | "+Str(*mousePosition\y))
    
    StopDrawing()
    Delay(1)
  EndIf
  
  ProcedureReturn result
  
EndProcedure


polygonPoints(polygonPoints)=polygonPoints(0)


wFlags = #PB_Window_SystemMenu | #PB_Window_ScreenCentered
If OpenWindow(0, 0, 0, 600, 630, "PointInPolygon()", wFlags)
  TextGadget(0, 0, 0, 600, 30, "", #PB_Text_Center)
  CanvasGadget(1, 0, 30, 600, 600, #PB_Canvas_Container)
  ;Kreis erstellen
  EinfacherKreis(polygonPoints(0), polygonPoints(1), polygonPoints(2), @Kreis)
  ;Debug "Radius: " + Kreis\r
  
  For i = 0 To 359
    polygonPoints(i)\x = ((Kreis\r) * Sin(i)) + Kreis\p\x
    polygonPoints(i)\y = ((Kreis\r) * Cos(i)) + Kreis\p\y
  Next
  
  Repeat
    event = WaitWindowEvent()
    Select event
      Case #PB_Event_CloseWindow
        appQuit = 1
      Case #PB_Event_Gadget
        If EventGadget() = 1 And EventType() = #PB_EventType_MouseMove
          mousePosition\x = GetGadgetAttribute(1, #PB_Canvas_MouseX)
          mousePosition\y = GetGadgetAttribute(1, #PB_Canvas_MouseY)
          If mouseState <> PointInPolygon(@polygonPoints(), polygonPoints, @mousePosition)
            mouseState ! 1
            If mousestate
              txt$ = "Mouse is in..."
            Else
              txt$ = "Mouse is out..."
            EndIf
            SetGadgetText(0, StrD(mousePosition\x)+" | "+StrD(mousePosition\y)+" >> "+txt$)
          EndIf
        EndIf
    EndSelect
  Until appQuit
EndIf
Betriebssysteme: div. Windows, Linux, Unix - Systeme

no Keyboard, press any key
no mouse, you need a cat
Benutzeravatar
chi
Beiträge: 90
Registriert: 17.05.2007 09:30
Wohnort: Linz - Austria

Re: Punkt in Polygon

Beitrag von chi »

Unter Windows könnte man auch CreatePolygonRgn und PtInRegion nutzen...

Code: Alles auswählen

Dim poly.POINT(2)

poly(0)\x = 10
poly(0)\y = 10
poly(1)\x = 100
poly(1)\y = 10
poly(2)\x = 50
poly(2)\y = 50

pRgn = CreatePolygonRgn_(@poly(), 3, #WINDING)

Debug PtInRegion_(pRgn, 50, 25)
Debug PtInRegion_(pRgn, 5, 5)

DeleteObject_(pRgn)
Antworten