Point in polygon

How fast can be that, which we call slow today!

Goal of the document:

Целта на този документ е да съхранява знанието натрупано за намиране дали една точка лежи в полигон. Не само синтезираната информация да се съхранява тук, но всякакви идеи и размисли свързани с темата.

Author

Daniel Sapoundjiev

Please, don't copy this content across the web, put links to it, cause its often revised! Thanks!

Date

26 January 2009

What really deeply means a point to lay in polygon

Ако нарисувате един квадрат и по средата му нарисувате точка, и попитате някой дали тази точка лежи в квадрата той ще ви отговори мигновено. Да, със сигурност е вътре в квадрата. Но ако след това го попитате, защо е така, да обясни по каква логика определя дали точката е в квадрата, може да получите разнообразни отговори и дори да всеете объркване. Това отново свежда нещата към това, че човек не винаги може да обясни каква логика използва мозъкът му в някой ситуации, но това е друга тема, да не се откланяме. След известно проучване и размисли сякаш се избистря следната логика. Човек просто проверява дали всички линии на полигона заобграждат точката и ако се опита от дадена отдалечена позиция, която е очевидно, че е извън полигона, да достигне до точката, то ще трябва да премине през някое от ребрата на полигона.

Хората обикновено ограждат нещо за да го защитят от достъп, или пък да не допуснат нещо да избяга от заграждението. Но също така може да се използва заграждение, за да се посочи, кое точно пространство е собственост на някой, също така кое точно пространство ще се ползва за нещо. Тука сякаш напипах по-високо ниво на абстракция относно загражденията, тъй като никога досега не съм разглеждал кройката на дреха като заграждение. Но дефакто се загражда една част от плата, после се изрязва и се шие.

Polygon as fence

Можем да си представим, че полигонът е ограда, която загражда определена площ земя. А всяка една точка от полигона е колче от оградата.

Are we inside the fence

Нека се опитаме да разберем дали ние сме вътре в ограденото пространство. Как бихме могли да разберем това.

Crossing the fence

Единият начин за който се сещам е, например когато идваме от друго място. Тогава за да се озовем в заграденото пространство е необходимо да прескочим оградата. След като един път сме я прескочили, то от тук нататък както и да си променяме позицията, то ние сме вътре в ограждението. Разбира се под промяна на позицията имам в предвид да не прескачаме оградата отново. Прескочим ли я, когато сме били вътре, то ще се озовем извън нея. Тогава отново можем да си променяме позицията, и докато не прескачаме оградата ще бъдем извън ограждението.

Да, по този начин можем лесно да уточним дали сме вътре или извън ограждение.

Surrounding by fence

Но на практика, при повечето практически задачи, често не разполагаме с информация дали сме били извън оградата, дали сме прескочили оградата. А по скоро ни поставят на някое място заедно с една ограда и трябва от това положение нататък да се ориентираме. Вътре ли сме или не.

В този случай единственото което виждаме около нас са колчетата от оградата и самата ограда. Гледайки тези елементи, трябва да се ориентираме дали сме вътре или вън. Начина това да се определи е да проследим с поглед цялата ограда, започвайки от едно произволно избрано колче и завършвайки отново в него. Следейки оградата ще се завъртаме към всяко следващо колче, и по този начин ще обходим цялата ограда. Като разбира се обхождането на колчетата става последователно. Целта на това е да разберем дали оградата ще ни заобгради. Ако оградата ни заобгражда, то ние ще направим едно пълно завъртане около оста си. Ако не ни заобгражда то колкото и да се завъртаме после ще се завъртаме обратно в същата посока така, че общото завъртане ще е 0 градуса.

Which is fast is fast

От казаното до тук стана ясно, че за да определим дали точка се намира в полигона, то трябва да обходим всички точки на полигона, и да измерим ъгъла, който образува проверяваната точка с тях.

Но тъй като в практиката, често се налага при наличието на един полигон, да се проверяват различни точки дали лежат в него, то би било добре да се оптимизира полигона така, че да не се налага всеки път да се обхождат всичките му точки.

Но как би станало това. Ами ако знаем, например, че поредица от 100 точки от полигона имат Х по-голям от този на проверяваната точка. То тогава можем спокойно да вземем ъгъла между първата и последната точка от тази поредица.

За да можем да се възползваме от това обаче, ще е необходимо да знаем за тези поредици, от нарастващи, или намаляващи позиции на Х координатата. Същото е валидно и за У координата, тя също може да ни спести обикалянето на много точки.

А може дори и още по-бързо, обхождане по метода на б-дърветата.

Here gets more exciting

Ами когато сме седнали на оградата, кой може да ни каже, вътре ли сме или вън. А от различни ситуации стигам до извода, че някой път, когато си на оградата се смята, че по-скоро си вътре, а друг път изискването е да се приеме, че си вънка. Това веднага навежда нещата в посока, че трябва алгоритъма да може да определя дали точката лежи в полигон, в зависимост от това дали е върху връх или ребро от полигона.

Разглеждания дотук алгоритъм позволява лесно да се включи и тази опция към него. Това става като се следи за получаването на определени ъгли при сумирането. И то по-точно ъгли като 180, -180 градуса.

More thoughts and passions:

I made a search in internet for these two words fence and enclosure and here is some information in wikipedia about them

http://en.wikipedia.org/wiki/Fence
http://en.wikipedia.org/wiki/Enclosure

 

And a little code in Pascal

Методът isPointIn обикаля всички точки от полигона, за да определи дали дадена точка лежи в него А методът isPointIn2 пропуска голяма част от тях.

На вас се пада честта да изпробвате разнообразни полигони и да ми пишете относно впечатленията си.

Написах тази статия за да съм ви полезен на вас, бъдете ми полезни и вие, като ми пишете своите впечатления. Така ще я подобря! На е-mail daniel@microdor.com

Питайте за това, което не ви е ясно, ще го опиша по-подробно за вас и за хората, които го четат след вас

unit poly;

interface

uses
Contnrs;

type

TPolyPoint = class
public
x, y: Double;
Next: TPolyPoint;
Index: Integer;
end;

TPolyPointDirection = class
public
first: TPolyPoint;
last: TPolyPoint;
isInc: Boolean;
end;

TPoly = class
private
fPoints: TObjectList;
fXDirections: TObjectList;
fYDirections: TObjectList;
protected
function getPoint(idx: Integer): TPolyPoint; virtual;
public
constructor Create(); virtual;
destructor Destroy(); override;
function isPointIn(aX, aY: Double): Boolean;
function isPointIn2(aX, aY: Double): Boolean;
procedure add(aX, aY: Double);
procedure clear;
function count: Integer;
procedure buildDirections();
property Points[idx: Integer]: TPolyPoint read getPoint;
end;

implementation

uses
Math;

{ TPoly }

constructor TPoly.Create;
begin
fPoints := TObjectList.Create;
fXDirections := TObjectList.Create;
fYDirections := TObjectList.Create;
end;

destructor TPoly.Destroy;
begin
fPoints.Free;
fXDirections.Free;
fYDirections.Free;
inherited;
end;

procedure TPoly.add(aX, aY: Double);
var
newPoint: TPolyPoint;
cnt: Integer;
begin
newPoint := TPolyPoint.Create();
newPoint.x := aX;
newPoint.y := aY;

cnt := count;
if cnt > 0 then begin
newPoint.Next := Points[0];
TPolyPoint(fPoints[cnt-1]).Next := newPoint;
end;
newPoint.Index := cnt;
fPoints.Add(newPoint);
end;

function TPoly.getPoint(idx: Integer): TPolyPoint;
begin
Result := TPolyPoint(fPoints[idx]);
end;

function TPoly.count: Integer;
begin
Result := fPoints.Count;
end;

procedure TPoly.clear;
begin
fPoints.Clear;
fXDirections.Clear;
fYDirections.Clear;
end;

procedure TPoly.buildDirections;
var
i: Integer;
pointDirectionX: TPolyPointDirection;
pointDirectionY: TPolyPointDirection;
begin
if count < 3 then
exit;

pointDirectionX := TPolyPointDirection.Create;
pointDirectionX.first := Points[0];
pointDirectionX.isInc := Points[0].x < Points[1].x;

pointDirectionY := TPolyPointDirection.Create;
pointDirectionY.first := Points[0];
pointDirectionY.isInc := Points[0].y < Points[1].y;

for i := 2 to count - 1 do begin

if pointDirectionX.isInc then begin
if (Points[i-1].x > Points[i].x) then begin
pointDirectionX.last := Points[i-1];
fXDirections.Add(pointDirectionX);

pointDirectionX := TPolyPointDirection.Create;
pointDirectionX.first := Points[i-1];
pointDirectionX.isInc := False;
end;
end
else begin
if (Points[i-1].x < Points[i].x) then begin
pointDirectionX.last := Points[i-1];
fXDirections.Add(pointDirectionX);

pointDirectionX := TPolyPointDirection.Create;
pointDirectionX.first := Points[i-1];
pointDirectionX.isInc := True;
end;
end;

if pointDirectionY.isInc then begin
if (Points[i-1].y > Points[i].y) then begin
pointDirectionY.last := Points[i-1];
fYDirections.Add(pointDirectionY);

pointDirectionY := TPolyPointDirection.Create;
pointDirectionY.first := Points[i-1];
pointDirectionY.isInc := False;
end;
end
else begin
if (Points[i-1].y < Points[i].y) then begin
pointDirectionY.last := Points[i-1];
fYDirections.Add(pointDirectionY);

pointDirectionY := TPolyPointDirection.Create;
pointDirectionY.first := Points[i-1];
pointDirectionY.isInc := True;
end;
end;

end;

if TPolyPointDirection(fXDirections[0]).isInc = pointDirectionX.isInc then
begin
TPolyPointDirection(fXDirections[0]).first := pointDirectionX.first;
pointDirectionX.Free;
end
else begin
pointDirectionX.last := Points[0];
fXDirections.Add(pointDirectionX);
end;

if TPolyPointDirection(fYDirections[0]).isInc = pointDirectionY.isInc then
begin
TPolyPointDirection(fYDirections[0]).first := pointDirectionY.first;
pointDirectionY.Free;
end
else begin
pointDirectionY.last := Points[0];
fYDirections.Add(pointDirectionY);
end;

end;

function TPoly.isPointIn(aX, aY: Double): Boolean;
var
i: Integer;
x1, y1, x2, y2, dp, cp, TotalSum: Double;
vPoint: TPolyPoint;
begin
if count < 3 then begin
Result := false;
exit;
end;

TotalSum := 0;
vPoint := Points[0];
x1 := (vPoint.X - aX);
y1 := (vPoint.Y - aY);
for i := 1 to count - 1 do begin
vPoint := points[i];
x2 := (vPoint.X - aX);
y2 := (vPoint.Y - aY);
cp := (x1*y2) - (y1*x2); // Cross-product
dp := (x1*x2) + (y1*y2); // Dot-product
TotalSum := TotalSum + ArcTan2(cp, dp);
x1 := x2;
y1 := y2;
end;

vPoint := points[0];
x2 := (vPoint.X - aX);
y2 := (vPoint.Y - aY);
cp := (x1*y2) - (y1*x2); // Cross-product
dp := (x1*x2) + (y1*y2); // Dot-product
TotalSum := TotalSum + ArcTan2(cp, dp);

Result := Abs(TotalSum) > 1;
end;

function TPoly.isPointIn2(aX, aY: Double): Boolean;
var
i: Integer;
vDirect: TPolyPointDirection;
x1, y1, x2, y2, dp, cp, TotalSum: Double;
vPoint: TPolyPoint;
found: Boolean;
index, indexStart, indexEnd: Integer;
begin
if count < 3 then begin
Result := false;
exit;
end;

TotalSum := 0;
if (fXDirections.Count <= fYDirections.Count) then begin
for i := 0 to fXDirections.Count - 1 do begin
vDirect := TPolyPointDirection(fXDirections[i]);

if vDirect.isInc then begin
if (vDirect.first.x > aX) or (vDirect.last.x < aX) then begin
x1 := (vDirect.first.X - aX);
y1 := (vDirect.first.Y - aY);
x2 := (vDirect.last.X - aX);
y2 := (vDirect.last.Y - aY);
cp := (x1*y2) - (y1*x2); // Cross-product
dp := (x1*x2) + (y1*y2); // Dot-product
TotalSum := TotalSum + ArcTan2(cp, dp);
end
else begin
indexStart := vDirect.first.Index;
indexEnd := vDirect.last.Index;
found := false;
repeat
if (indexStart < indexEnd) then
index := (indexStart + indexEnd) div 2
else begin
index := (count + indexEnd + indexStart) div 2;
if index > count then
index := index - count;
end;
if index = indexStart then begin
x1 := (Points[indexStart].X - aX);
y1 := (Points[indexStart].Y - aY);
x2 := (Points[indexEnd].X - aX);
y2 := (Points[indexEnd].Y - aY);
found := true;
end
else
if (Points[index].x < aX) then begin
x1 := (Points[indexStart].X - aX);
y1 := (Points[indexStart].Y - aY);
x2 := (Points[index].X - aX);
y2 := (Points[index].Y - aY);
indexStart := index;
end
else begin
x1 := (Points[index].X - aX);
y1 := (Points[index].Y - aY);
x2 := (Points[indexEnd].X - aX);
y2 := (Points[indexEnd].Y - aY);
indexEnd := index;
end;
cp := (x1*y2) - (y1*x2); // Cross-product
dp := (x1*x2) + (y1*y2); // Dot-product
TotalSum := TotalSum + ArcTan2(cp, dp);
until found;
end
end
else begin
if (vDirect.first.x < aX) or (vDirect.last.x > aX) then begin
x1 := (vDirect.first.X - aX);
y1 := (vDirect.first.Y - aY);
x2 := (vDirect.last.X - aX);
y2 := (vDirect.last.Y - aY);
cp := (x1*y2) - (y1*x2); // Cross-product
dp := (x1*x2) + (y1*y2); // Dot-product
TotalSum := TotalSum + ArcTan2(cp, dp);
end
else begin
indexStart := vDirect.first.Index;
indexEnd := vDirect.last.Index;
found := false;
repeat
if (indexStart < indexEnd) then
index := (indexStart + indexEnd) div 2
else begin
index := (count + indexEnd + indexStart) div 2;
if index > count then
index := index - count;
end;
if index = indexStart then begin
x1 := (Points[indexStart].X - aX);
y1 := (Points[indexStart].Y - aY);
x2 := (Points[indexEnd].X - aX);
y2 := (Points[indexEnd].Y - aY);
found := true;
end
else
if (Points[index].x > aX) then begin
x1 := (Points[indexStart].X - aX);
y1 := (Points[indexStart].Y - aY);
x2 := (Points[index].X - aX);
y2 := (Points[index].Y - aY);
indexStart := index;
end
else begin
x1 := (Points[index].X - aX);
y1 := (Points[index].Y - aY);
x2 := (Points[indexEnd].X - aX);
y2 := (Points[indexEnd].Y - aY);
indexEnd := index;
end;
cp := (x1*y2) - (y1*x2); // Cross-product
dp := (x1*x2) + (y1*y2); // Dot-product
TotalSum := TotalSum + ArcTan2(cp, dp);
until found;
end;
end;
end;
end
else begin
for i := 0 to fYDirections.Count - 1 do begin
vDirect := TPolyPointDirection(fYDirections[i]);
if ((vDirect.first.y < aY) and (vDirect.last.y < aY)) or
((vDirect.first.y > aY) and (vDirect.last.y > aY)) then
begin
x1 := (vDirect.first.X - aX);
y1 := (vDirect.first.Y - aY);
x2 := (vDirect.last.X - aX);
y2 := (vDirect.last.Y - aY);
cp := (x1*y2) - (y1*x2); // Cross-product
dp := (x1*x2) + (y1*y2); // Dot-product
TotalSum := TotalSum + ArcTan2(cp, dp);
end
else begin
vPoint := vDirect.first;
x1 := (vPoint.X - aX);
y1 := (vPoint.Y - aY);
repeat
vPoint := vPoint.Next;
x2 := (vPoint.X - aX);
y2 := (vPoint.Y - aY);
cp := (x1*y2) - (y1*x2); // Cross-product
dp := (x1*x2) + (y1*y2); // Dot-product
TotalSum := TotalSum + ArcTan2(cp, dp);
x1 := x2;
y1 := y2;
until (vPoint <> vDirect.last);
end;
end;
end;
Result := Abs(TotalSum) > 1;
end;

end.

Back to main menu