Алгоритм определения попадания точки в контур на основе комплексного анализа

Привет всем Хабра людям. Хочу представить уважаемым читателям пример, когда сухая и далекая от жизни в нашем понимании высшая математика дала не плохой практический результат.

Сначала немного воспоминаний

Было это в бытность мою студентом одного из технических Вузов в 90-е, курсе наверно втором. Попал я как-то на олимпиаду по программированию. И вот на этой самой олимпиаде и было задача: задать координаты треугольника, тестовой точки на плоскости, и определить принадлежит ли эта точка области треугольника. В общем, плевая задачка, но тогда я ее так и не решил. Но после задумался – над более общей задачей – принадлежность полигону. Повторюсь – была середина 90 –х, интернета не было, книжек по компьютерной геометрии не было, а были лекции по вышке и лаборатория 286 –х с турбо паскалем. И вот так совпали звезды, что как раз в то время когда я размышлял над проблемой, на вышке нам читали теорию комплексного переменного. И одна формула (о ней ниже) упала на благодатную почву. Алгоритм был придуман и реализован на паскале (к сожалению мой полутора гиговый винт погиб и унес в небытие этот код и кучу других моих юношеских наработок). После института я попал работать в один НИИ. Там мне пришлось заниматься разработкой ГИС для нужд работников института и собственной одной из задачей было определение попадания объектов в контур. Алгоритм был переписан на С++ и отлично зарекомендовал себя в работе.

Задача для алгоритма

Дано:

Г- замкнутая ломаная (далее полигон) на плоскости, заданная координатами своих вершин (xi,yi), и координата тестовой точки (x0,y0)

Определить:

принадлежит ли точка области D, ограниченной полигоном.

Вывод формул для последующего написания алгоритма ни в коем случае не претендует на математическую полноту и точность, а лишь демонстрирует инженерный (потребительский подход) к Царице полей наук.

И так начнем

Интегральная формула Коши

Пояснение с рабоче-крестьянской инженерной точки зрения:
— граница Г наш заданный контур,
— z0 -тестируемая точка
— f(z) — комплексная функция от комплексного аргумента нигде в контуре не обращается в бесконечность.

Те есть, чтобы установить принадлежность точки контуру, нам необходимо вычислить интеграл и сравнить его со значением функции в данной точки. Если они совпадают, то точка лежит в контуре. Замечание: интегральная теорема коши гласит, что если точка не лежит в контуре, те подынтегральное выражение нигде не обращается в бесконечность, то интеграл равен нулю. Это упрощает дело – нужно лишь вычислить интеграл и проверить его на равенство нулю: равен нулю точка не контура, отличен — лежит в контуре.
Займемся вычислением интеграла. За f(z) примем простую функцию 1. Не нарушая общности можно за z0 принять точку 0 (всегда можно сдвинуть координаты).


Избавляемся от мнимой единицы в знаменателе подынтегральной части и расщепим интеграл на действительную и мнимую части:



Получилось два криволинейных интеграла II рода.
Вычислим первый


Выполнятся условие не зависимости интеграла от пути, следовательно, первый интеграл равен нулю и его вычислять не нужно.

С мнимой частью такой фокус не проходит. Вспоминаем, что наша граница состоит из отрезков прямых, получаем:

Где Гi- это отрезок (xi,yi)- (xi+1,y i+1)
Вычислим i-ый интеграл. Для этого запишем уравнение i-го отрезка в параметрическом виде


Подставим в интеграл


и после громоздких и нудных преобразований получим следующую прельстивую формулу:


Окончательно получаем

Алгоритм на C++:

template <class T>
bool pt_in_polygon(const T &test,const std::vector &polygon)
{
if (polygon.size()<3) return false;

std::vector::const_iterator end=polygon.end();

T last_pt=polygon.back();

last_pt.x-=test.x;
last_pt.y-=test.y;

double sum=0.0;

for(
std::vector::const_iterator iter=polygon.begin();
iter!=end;
++iter
)
{
T cur_pt=*iter;
cur_pt.x-=test.x;
cur_pt.y-=test.y;

double del= last_pt.x*cur_pt.y-cur_pt.x*last_pt.y;
double xy= cur_pt.x*last_pt.x+cur_pt.y*last_pt.y;

sum+=
(
atan((last_pt.x*last_pt.x+last_pt.y*last_pt.y - xy)/del)+
atan((cur_pt.x*cur_pt.x+cur_pt.y*cur_pt.y- xy )/del)
);

last_pt=cur_pt;

}

return fabs(sum)>eps;

}



T – тип точки, например:
struct PointD
{
double x,y;
};


Пример

Пример работы алгоритма написан с применением самой на мой взгляд великой библиотеки 2D графики:Anti-Grain Geometry (AGG) .

Управление:
клик левой кнопкой – добавление новой точки контура
правой кнопкой - замыкание контура
левой с зажатым Shift-ом – перенос тестовой точки

PS

Господа, кому интересно, привожу более быстрый алгоритм. Уже не мой.
Отдельное и огромное спасибо forgotten за статейку.
template bool pt_in_polygon2(const T &test,const std::vector &polygon)
{

static const int q_patt[2][2]= { {0,1}, {3,2} };

if (polygon.size()<3) return false;

std::vector::const_iterator end=polygon.end();
T pred_pt=polygon.back();
pred_pt.x-=test.x;
pred_pt.y-=test.y;

int pred_q=q_patt[pred_pt.y<0][pred_pt.x<0];

int w=0;

for(std::vector::const_iterator iter=polygon.begin();iter!=end;++iter)
{
T cur_pt = *iter;

cur_pt.x-=test.x;
cur_pt.y-=test.y;

int q=q_patt[cur_pt.y<0][cur_pt.x<0];

switch (q-pred_q)
{
case -3:++w;break;
case 3:--w;break;
case -2:if(pred_pt.x*cur_pt.y>=pred_pt.y*cur_pt.x) ++w;break;
case 2:if(!(pred_pt.x*cur_pt.y>=pred_pt.y*cur_pt.x)) --w;break;
}

pred_pt = cur_pt;
pred_q = q;

}

return w!=0;

}