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