Berekenen of coordinaten binnen een polygoon liggen *

[.Alex]

Goedemiddag,

Ik ben aan het nadenken over het volgende maar ik kom er maar niet uit:

1 - Ik stop een aantal coördinaten in een applicatie. (min 3 coördinaten, maximaal een heleboel) Dit zijn de coördinaten van een vlak dat dus in alle vormen en maten kan zijn.
2 - Nu stop ik nog een coördinaat in de applicatie en wil ik weten of deze in het vlak valt.

Het los ik dit wiskunde technisch op? Hoe bereken ik of het tweede coördinaat zich in het vlak bevindt?

Alvast enorm bedankt voor jullie input!

edit:
Ik hoop trouwens dat dit in het goede forum staat...

[ Voor 5% gewijzigd door .Alex op 25-06-2007 12:27 ]

[soulrider]

elk vlak kan omschreven worden door een functie.

mbv je opgegeven punten ga je die functie moeten zoeken, en nadien kijken of het extra punt daar ook aan voldoet.
(met 3 punten heb je genoeg voor een plat vlak, punten voor een lijn - dus als het enkel platte vlakken zijn ga je ook die extra coördinaten al moeten checken)


dat eerste is geen simpele opdracht. - maar na een opfrissing van je hogere wiskunde wrs wel mogelijk. (zekers niet als je vlak geen plat vlak is (wat ev. gekandelt is rond verschillende assen) )

edit: je wilt dus weten of een punt ligt in het gebied afgebakend door de lijnen die door de opgegeven punten gaan??

(als voorbeeld: als je de coordinaten van de 4 hoeken van een vierkant opgeeft, wil je weten of het apart opgegeven punt binnen of buiten dat vierkant ligt ?
Dan kan je ook gebruik maken van de wiskunde:
elke zijlijn is wiskundig formuleerbaar, en het vlak aan 1 der zijdes is verkrijgbaar door de = te vervangen door <= of => afh. van welke zijde je moet hebben.
(bijvoorbeeld de lijn: x=0 (de y-as): het vlak links van ervan is: x<=0, rechts ervan x >= 0)

door zo van al die zijlijnen de voorwaarden op te stellen (uit te zoeken) weet je onmiddelijk of de coordinaten van dat extra punt voldoet aan al die voorwaarden, indien niet ligt ie niet in dat afgebakend vlak.

bv: vierkant met punten 0,0; 0,2 ; 2,0 en 2,2 is omschrijfbaar als:
x>=0, x<=2, y>=0 en y<=2 elk punt (bv 1,1) dat voldoet aan die 4 voorwaarden zit in dat vierkant

[ Voor 75% gewijzigd door soulrider op 25-06-2007 14:59 ]

[RobIII]

Je wil dus berekenen of een punt x,y (,z?) in een bepaalde polygoon valt? Met de juiste termen googled namelijk een stuk makkelijker Voor 3D heb je hier wat leesmateriaal om te beginnen

[ Voor 38% gewijzigd door RobIII op 25-06-2007 12:41 ]

[.oisyn]

Wat bedoel je met 'in' het vlak? Ik heb het idee dat je eigenlijk polygoon bedoelt en niet vlak - een vlak is oneindig, een polygoon is de subset van punten op een vlak dat binnen de lijnen van de polygoon ligt.

Hoe zijn je punten die het polygoon definieren gesorteerd? Zijn deze altijd met de klok mee danwel tegen de klok in? Het simpelste is dan gewoon hoeken bij elkaar optellen - als P het nieuwe punt is, en V_i de punten van je polygoon met i van 0 t/m N-1, bereken dan de som van de hoeken van lijnstuk P.V_i met lijnstuk P.V_i+1. Als je op 2pi (of -2pi, afhankelijk van hoe je het berekent en hoe je punten gesorteerd zijn) dan ligt het punt binnen de polygoon, anders erbuiten.

Als je echt gewoon vlak bedoelt moet je de normaal van het vlak berekenen (uitwendig product of "kruisproduct" tussen 2 lijnstukken gedefinieerd door 3 willekeurige punten op dat vlak) en de afstand tot het nulpunt (inwendig product of "dotproduct" van de normaal en een willekeurig punt op het vlak). Een ander punt P ligt dan op het vlak als het inwendig product van dat punt met de normaal dezelfde afstand oplevert als de afstand van het vlak tot het nulpunt dat je eerder berekend hebt (het verschil in die afstanden is dan ook eigenlijk de afstand van dat punt tot het vlak - als die 0 is dan ligt het natuurlijk op het vlak).

[ Voor 30% gewijzigd door .oisyn op 25-06-2007 12:41 ]

[.Alex]

Heb nog niet alles gelezen, heb een meeting over 5 minuten. RobIII verwoordt precies wat ik wil. Het gaat om 2d. Om het wat concreter te maken zou je het kunnen vergelijken met GPS coördinaten. Ik voer een coördinaat in en de applicatie moet aangeven of deze een van tevoren aangegeven vlak valt.

Ik kom er na de meeting op terug

[Janoz]

Als het puur het vlak is heb je de volgende formule:

ax+by+cy=d

Waarbij (a,b,c) de normaal op het vlak is en d is de afstand van het vlak tot de oorsprong (in aantal keren de eerder genoemde normaal).

Kijken of een punt op dat vlak ligt is gewoon een kwestie van invullen.


Ja, dan is het natuurlijk een ander verhaal. Daarom is terminologie ook zo belangrijk. Een vlak is geen polygoon.

Om van een berg punten een polygoon te maken moet je in eerste instantie bepalen wat nu eigenlijk je vlak is. Liggen al die punten op de rand, of is het een wolk waarvan je de convex hull(*) wilt. Als je eenmaal je polygoon hebt is het vervolgens van belang of deze convex is (in dat geval is het een simpel 'links van elke edge' algoritme) of niet.

*convex hull van een set punten is alsof je voor elk punt een spijkertje in de plank slaat en er vervolgens een elastiekje omheen doet. Wiskundige omschrijving is 'Voor elke twee punten die binnen het polygoon liggen ligt de lijn tussen die twee punten ook helemaal in het polygoon'.

[ Voor 59% gewijzigd door RobIII op 25-06-2007 14:12 . Reden: Linkje naar wikipedia toegevoegd :+ ]

[MrBucket]

Een makkelijk te implementeren methode is door een ray* te schieten in een willekeurige richting vanuit je query point Q, en te kijken hoevaak hij de edges van je polygon snijdt. Als dit een even aantal is, dan ligt Q buiten je polygon, als dit een oneven aantal is, dan ligt Q in je polygon. (Logisch, toch?).

Dit ray schieten kun je implementeren door een "lang genoeg" lijnsegment S te nemen waarvan Q een van de eindpunten is, en het andere eindpunt P bijvoorbeeld buiten de bounding box van {Q en je polygon} ligt. Vervolgens loop je door alle lijnsegmenten uit je polygon en hou je bij hoeveel er daarvan snijden met S.

Let wel even op met de 'degenerate case' wanneer S een hoekpunt van je polygon snijdt. In dat geval is het niet duidelijk of dit als 1 of als 2 snijpunten moet tellen - om je een hoop hoofdbrekens te sparen is de makkelijkste oplossing: kies een willekeurig ander eindpunt P en probeer het nog een keer.

* een Ray is een lijnsegment met maar 1 eindpunt; de andere kant loopt tot in het oneindige door.

[Anoniem: 1572]

Als het daadwerkelijk om gps coordinaten gaat maakt dit het probleem een stuk ingewikkelder. Een polygoon is dan namelijk niet 2d omdat hij op een bol ligt. Hierdoor gaat de standaard wiskunde niet op, een van de axioma's van wiskunde houd namelijk niet meer: 2 parallele lijnen kunnen elkaar wel degelijk kruizen. (kijk maar naar de lengte graden).

[Teun_2]

Wat je wil is nog niet zo simpel. Ik heb daar over een dik uur een examen over.
Wat je eerst moet doen is je polgoon opdelen in allemaal 3-hoeken. (Zoeken op google naar triangulatie).
Vervolgens kan je redelijk eenvoudig per driehoek zien of een punt in die driehoek ligt.

Straks, na mijn examen computer graphics heb ik misschien even wat meer tijd.

[Dido]

Wat Zoutvat zei was ook het eerste wat in me opkwam, eigenlijk. Het bepalen van al je driehoeken mag niet moeilijk zijn (simpelweg alle mogelijkheden aflopen moet te doen zijn) en bepalen of je punt binnen een driehoek ligt is ook niet het einde van de wereld, lijkt me.

Het enige wat me wat tegenstaat is dat het wel een brute-force benadering is die veel rekentijd gaat kosten als je veel punten gebruikt...

[Haan]

In bijvoorbeeld Java heb je daar toch een mooie kant-en-klare functie voor?

[Janoz]

@dido & zoutvat: Je hoeft niet helemaal te trianguleren. Opdelen in convex polygonen is genoeg. Dat algoritme is ook een stuk efficienter dan triangulatie algoritmen.

En als we toch gaan pochen met examens: * Janoz is afgestudeerd in Scientific computing and visualisation waar Computer Graphics een groot deel van uitmaakt.

[MrBucket]

Janoz schreef op dinsdag 26 juni 2007 @ 08:50:
@dido & zoutvat: Je hoeft niet helemaal te trianguleren. Opdelen in convex polygonen is genoeg. Dat algoritme is ook een stuk efficienter dan triangulatie algoritmen.

En als we toch gaan pochen met examens: * Janoz is afgestudeerd in Scientific computing and visualisation waar Computer Graphics een groot deel van uitmaakt.

Dan zou je moeten weten dat het opdelen in monotone polygonen zowel makkelijker te implementeren als efficienter is, en dat je deze ook in lineaire tijd kunt trianguleren

Maar on-topic: een triangulatie maken van een polygon, alleen om te kijken of een enkel query-punt erin ligt is behoorlijk overkill, en niet te vergeten, redelijk complex om te implementeren als je geen library hebt met geometrische functies.

[Janoz]

Het trianguleren van monotone polygonen (wat polygonen zijn die als het ware convex in 1 richting zijn) is inderdaad lineair, maar het splitsen in monotone polygonen is maar iets makkelijker dan het splitsen in convex polygonen. Dit verlies wordt weer ongedaan gemaakt doordat een triangulatie uiteindelijk veel meer edges krijgt die je moet gebruiken bij het controleren of een punt binnen 1 van de geconstrueerde polygonen ligt.

Maar inderdaad, het zal niet veel schelen.

[MrBucket]

Janoz schreef op dinsdag 26 juni 2007 @ 09:33:
Het trianguleren van monotone polygonen (wat polygonen zijn die als het ware convex in 1 richting zijn) is inderdaad lineair, maar het splitsen in monotone polygonen is maar iets makkelijker dan het splitsen in convex polygonen. Dit verlies wordt weer ongedaan gemaakt doordat een triangulatie uiteindelijk veel meer edges krijgt die je moet gebruiken bij het controleren of een punt binnen 1 van de geconstrueerde polygonen ligt.

Maar inderdaad, het zal niet veel schelen.

Hmm, ik was in de veronderstelling dat een polygon opdelen in convexe stukken nog niet zo eenvoudig was. Kun je een globale omschrijving van een algoritme geven, of heb je misschien een link naar een algoritme om dit te doen?

[Janoz]

Door simpel langs de edge te wandelen kun je telkens stukken van het polygoon afknippen. Vergelijk het met het algoritme om een convex hull te maken. Uitgaande van een polygoon definitie met de klok mee. Je begint met twee punten. Het volgende punt van het resultaat is het eerst volgende punt dat rechts van de lijn ligt. Alles wat je overslaat maak je een nieuw polygoon van en zet je in een queue om straks nog even te controleren. Als je edge gedefinieerd is als een linked list is het gewoon een kwestie van de list doorknippen in twee apparte lists waarvan je bij elk de uiteindes weer aan elkaar plakt.

Dat is iig een algoritme wat ik zo even uit mijn hoofd kan bedenken. Er zullen er vast nog wel anderen zijn.

[Soultaker]

.oisyn schreef op maandag 25 juni 2007 @ 12:37:
Als je echt gewoon vlak bedoelt moet je de normaal van het vlak berekenen (uitwendig product of "kruisproduct" tussen 2 lijnstukken gedefinieerd door 3 willekeurige punten op dat vlak) en de afstand tot het nulpunt (inwendig product of "dotproduct" van de normaal en een willekeurig punt op het vlak). [..]

Deze methode klinkt interessant, maar ik volg 'm niet helemaal. Het inwendig project van twee vectoren die haaks op elkaar staan is toch altijd 0?

MrBucket schreef op maandag 25 juni 2007 @ 20:50:
Een makkelijk te implementeren methode is door een ray* te schieten in een willekeurige richting vanuit je query point Q, en te kijken hoevaak hij de edges van je polygon snijdt. Als dit een even aantal is, dan ligt Q buiten je polygon, als dit een oneven aantal is, dan ligt Q in je polygon.

Dit is inderaad één standaardmethode; de andere is bepalen of het winding number (hoe vaak gaat de polygon om Q heen?) nul is. Hier staat de methode wel goed uitgelegd, hoewel het geen lichte kost is: Fast Winding Number Inclusion of a Point in a Polygon

[.oisyn]

Soultaker schreef op dinsdag 26 juni 2007 @ 10:24:

Deze methode klinkt interessant, maar ik volg 'm niet helemaal. Het inwendig project van twee vectoren die haaks op elkaar staan is toch altijd 0?

Klopt, maar hoe brengt dat jou in de war dan?

Soultaker schreef op dinsdag 26 juni 2007 @ 10:24:
Dit is inderaad één standaardmethode; de andere is bepalen of het winding number (hoe vaak gaat de polygon om Q heen?) nul is.

Hey, de methode die ik als eerst in de draad noemde is óók een standaard-methode hoor . Én makkelijk te implementeren:

C++:

bool pointInPolygon(const vector & p, const vector * pVerts, uint32 numVerts)
{
    float angle = 0.f;
    for (int i = 0, lastI = numVerts - 1; i < numVerts; lastI = i++)
    {
        vector e0 = pVerts[lastI] - p;
        vector e1 = pVerts[i] - p;
        float z = e0.x * e1.y - e0.y * e1.x;
        angle += asin(z / (length(e0) * length(e1));
    }

    return abs(angle) > 0.1f;
}

Deze werkt ongeacht de winding order van de polygoon (CW of CCW). Als het punt erbuiten staat zou je op een totale hoek van 0 uit moeten komen. Als het punt erbinnen staat kom je op 2pi of -2pi uit, afhankelijk van de winding order en de handedness van je coordinatenstelsel.

[ Voor 72% gewijzigd door .oisyn op 26-06-2007 11:50 ]

[.Alex]

Wow! Ik wil jullie om te beginnen allemaal bedanken voor jullie reacties en moet eerlijk toegeven dat ik meer dan de helft niet kon volgen (vooral vanwege de terminologie) maar goed, dat geeft niet. Ik geloof dat ik er op zich wel uit kom met de 'cutting edges'-methode. Die is relatief simpel te programmeren en het mag nou eenmaal niet teveel tijd kosten :-)

Waarom zou ik een andere functie gebruiken eigenlijk? Volgens mij voldoet deze prima. (Ja, schiet maar )

Op zich heb ik nu een (tijdelijke) oplossing, maar mochten jullie willen doorgaan met discussiëren in dit topic dan volg ik de discussie met alle plezier (Gratis kant-en-klare leermomentjes zeg maar)

edit:
En ja, ik heb hier gelezen wat de verschillen zijn tussen de "Winding Numbers" en "Cutting Edges".
http://www.geometryalgori...m_0103/algorithm_0103.htm

Prima linkje!

[ Voor 13% gewijzigd door .Alex op 26-06-2007 11:39 ]

[.oisyn]

Lees mijn edit nog even

[Janoz]

Bedenk dat dergelijke algoritmen niet zomaar toepasbaar zijn op gps coordinaten. Er wordt bij dergelijke algoritmen namelijk van een plat vlak uitgegaan. Dit betekent bij het gerbuik van gps coordinaten dat je even goed moet letten op het 'wrappen' van je raster (360 graden = 0 graden, Wanneer je in het overgangs gebied zit kun je natuurlijk gewoon bij het ene deel 360 optellen, bij de polen heb je echter een groter probleem) en dat de randen van je polygonen op de aarde uiteindelijk niet recht lopen.

[Soultaker]

Janoz schreef op dinsdag 26 juni 2007 @ 11:52:
Bedenk dat dergelijke algoritmen niet zomaar toepasbaar zijn op gps coordinaten. Er wordt bij dergelijke algoritmen namelijk van een plat vlak uitgegaan.

Volgens mij werken de meeste genoemde algoritmen prima met coördinaten die uit willekeurige 2D ruimtes afkomstig zijn, zolang ze maar ingebed kunnen worden op zijn met een plat vlak, en in dat vlak de edges recht zijn. Dat geldt dus niet voor een hele bol, maar wel voor een deel van een bol (bijvoorbeeld: een gebied in Nederland).

.oisyn schreef op dinsdag 26 juni 2007 @ 11:35:
Klopt, maar hoe brengt dat jou in de war dan?

Ik volg niet helemaal welke vectoren je nu neemt en vermenigvuldigt; ik zal er later nog eens naar kijken en als ik het dan nog niet snap kom ik er op terug.

Hey, de methode die ik als eerst in de draad noemde is óók een standaard-methode hoor . Én makkelijk te implementeren

Dat is inderdaad een heel handige implementatie; maar is gewoon het bepalen van een winding number toch? (Wel makkelijk te coden zo, al zal het niet het snelste zijn met al die calls van asin().)

edit:
Moet dat geen acos() zijn trouwens?

[.oisyn]

Anoniem: 1572 schreef op dinsdag 26 juni 2007 @ 06:01:
Als het daadwerkelijk om gps coordinaten gaat maakt dit het probleem een stuk ingewikkelder. Een polygoon is dan namelijk niet 2d omdat hij op een bol ligt. Hierdoor gaat de standaard wiskunde niet op, een van de axioma's van wiskunde houd namelijk niet meer: 2 parallele lijnen kunnen elkaar wel degelijk kruizen. (kijk maar naar de lengte graden).

Dat is alleen een probleem als een van de polen in je polygoon liggen, of het te testen punt ligt op een van de polen. Het daadwerkelijke probleem is dat de edges van de polygoon de kortste cirkelboog van het ene punt op de bol naar het andere volgen, en dat is niet hetzelfde als een rechte lijn in GPS coordinaten. Als je de pool-cases kunt negeren, dan is dit probleem met een simpele Mercator projectie op te lossen zodat lijnen wel recht lopen.

Zoutvat schreef op dinsdag 26 juni 2007 @ 07:46:
triangulatie

Janoz schreef op dinsdag 26 juni 2007 @ 08:50:
Opdelen in convex polygonen is genoeg.

MrBucket schreef op dinsdag 26 juni 2007 @ 09:08:
monotone polygonen

Wat nodeloos ingewikkeld allemaal . Je hoeft helemaal niet op te delen, en je kunt de cutting edges methode gewoon gebruiken met het daadwerkelijke polygoon. Je moet het alleen niet in 2D doen met poly edges, maar in 3D met vlakken. Projecteer de GPS coordinaten op een unit sphere, en maak vlakken door 2 aanliggende coordinaten en de oorsprong. Vervolgens kun je gewoon het aantal snijpunten berekenen met een straal vanuit het te testen punt in een richting loodrecht op de richting naar de oorsprong.

[.oisyn]

Soultaker schreef op dinsdag 26 juni 2007 @ 12:05:
Ik volg niet helemaal welke vectoren je nu neemt en vermenigvuldigt; ik zal er later nog eens naar kijken en als ik het dan nog niet snap kom ik er op terug.

Afstand vlak tot oorsprong is N*P, waarbij N de normaal is en P een willekeurig punt op het vlak (dus de vector vanuit de oorsprong naar het punt). Als P haaks op N staat, dan impliceert dat idd dat de afstand tot de oorsprong 0 is. Want als je haaks op de normaal "beweegt", dan blijf je in het vlak, dus de oorsprong ligt dan ook in het vlak.

Dat is inderdaad een heel handige implementatie; maar is gewoon het bepalen van een winding number toch? (Wel makkelijk te coden zo, al zal het niet het snelste zijn met al die calls van asin().)

Met winding number bereken je toch het aantal snijpunten, of ben ik nou gek?
.edit: nee, my bad. Je hebt idd gelijk

Moet dat geen acos() zijn trouwens?

Nee asin.
P∙Q = cos(a)∙|P|∙|Q|
|P×Q| = sin(a)∙|P|∙|Q|

Ik gebruik hier het crossproduct (wat met 2d coordinaten gewoon een 'virtueel' z-coordinaat als antwoord geeft) omdat je daarmee het verschil tussen linksom en rechtsom kunt bepalen. Bij het inproduct kan dat niet.

[roy-t]

2D?

Stel A is het vlak, Z is het punt

Visual Basic:

A.left = 10
A.top = 10 'relatieve positie 
A.width =50
A.height = 50

Z.left = 20
Z.top = 20  'ligt er dus in

if Z.top < ( A.top + A.height ) then
    if Z.top > ( A.top ) then
       'controleer breedte
        if Z.left < ( A.left + A.width ) then
            if Z.left > (A.left ) then
             'Jaaa hij zit erin
             msgbox "Ik zit erin!" , vbinformation
            End if
        End if
      End if
End if
  

Als je ook 3d, of ingewikkelde vormen wilt hebben, zou ik even naar C# XNA documentatie kijken daar staan veel checks in, en als je toevallig C# gebruikt kom je dan een heel eind komen.

Edit: reageert te laat, en te kort door de bocht

[ Voor 3% gewijzigd door roy-t op 26-06-2007 12:45 ]

[.oisyn]

Edit: reageert te laat, en te kort door de bocht

No shit . Je test hier een punt in een rechthoek.

[Janoz]

Soultaker schreef op dinsdag 26 juni 2007 @ 12:05:
Volgens mij werken de meeste genoemde algoritmen prima met coördinaten die uit willekeurige 2D ruimtes afkomstig zijn, zolang ze maar ingebed kunnen worden op zijn met een plat vlak, en in dat vlak de edges recht zijn. Dat geldt dus niet voor een hele bol, maar wel voor een deel van een bol (bijvoorbeeld: een gebied in Nederland).

Als je graden gewoon als roosterpunten interpreteert, zullen de edges waarvan in het algoritme uitgegaan wordt dat ze recht lopen in het echt niet recht lopen. Dit hoeft inderdaad niet uit te maken, daarom geef ik ook aan dat ze niet zomaar toepasbaar zijn. Je moet even controleren of de mogelijke afwijkingen binnen een gewenste marges vallen.

[.oisyn]

Die edges zijn dus wel recht te krijgen middels een mercator projectie, zoals ik al zei.

p.x = GPS.long;
p.y = .5 * ln((1 + sin(GPS.lat)) / (1 - sin(GPS.lat)))

.edit: ik bedenk me net dat dit niet klopt, latitude-lijnen zijn niet geodetisch

[ Voor 31% gewijzigd door .oisyn op 26-06-2007 14:40 ]