nyaload

Журнал Пушыстого

Журнал Пушыстого

Previous Entry Share Next Entry
Попадание точки в многоугольник. Количество заныриваний не равно количеству выныриваний.
nyaload
_winnie
Есть точка O, есть многоугольник A[n] на плоскости. Надо выяснить, попадает ли точка в многоугольник.

Основа алгоритма описывается очень просто: давайте возмём бесконечный горизонтальный луч из точки направо. Если луч пересекает границу четное число раз (в том числе 0 раз) , то точка вне многоугольника. Если нечетное число раз - то внутри.



Затем программист пытается обработать хитрые случаи: луч проходит через вершину многоугольника, сторона многоугольника совпадает с лучом. В коде возникают какие-то if, какие-то EPSILON..., молитвы о том, что редкие случаи когда две точки на расстоянии EPSILON но луч между ними - не существуют, насильно изгонясь из сознания

Надо так: надо работать не со сторонами, а с вершинами. Для каждой вершины A[n] - вычисляем булевский флажок, находится точка над или под лучом, то есть просто сравниваем координату y вершины с нулем, flag[n] = A[n] < 0. Затем идем по парам точек (A[n-1],A[n]) и считаем пересечения когда этот флажок меняется с true на false и наоборот. При этом проверяем, что отрезок проходит справа от точки O, для этого сравниваем знак определителя |A[n-1]-A[n], A[n]-O| с 0. Определитель получается не из божественного откровения, а из выражения для x-координаты пересечения отрезка и луча.

Алгоритм устойчив к дрожанию вершины на любой epsilon, к случаю когда она попадает почти на луч или ровно на луч. Например, на картинке 1) ниже - будет засчитано или ровно два пересечения, или ровно ноль. Даже если точка попадает ровно на луч. Благодаря дискретной природе булевского значения в точке - оно или true, или false, и количество их переключений - строго дискретно. Или меняется сразу на 2, но никак на 1. На картинке 2) - будет зачитано ровно одно пересечение, или на отрезке AB или на отрезке BC. Ни в коем случае не на двух сразу, и ни в коем случае не будет случайно пропущен один из двух отрезков.



Функция целиком, 30 строчек, никакой лапши из epsilon и особых случаев:

bool IsPointInside(Array<const Point2f> polygon, Point2f point)
{
    if (polygon.size <= 1)
        return false;

    int intersections_num = 0;
    int prev = polygon.size - 1;
    bool prev_under = polygon[prev].y < point.y;

    for (int i = 0; i < polygon.size; ++i)
    {
        bool cur_under = polygon[i].y < point.y;

        Point2f a = polygon[prev] - point;
        Point2f b = polygon[i]  - point;

        float t = (a.x*(b.y - a.y) - a.y*(b.x - a.x));
        if (cur_under && !prev_under)
        {
            if (t > 0)
                intersections_num += 1;
        }
        if (!cur_under && prev_under)
        {
            if (t < 0)
                intersections_num += 1;
        }

        prev = i;        
        prev_under = cur_under;        
    }

    return (intersections_num&1) != 0;
}
_Winnie C++ Colorizer



Смотри ещё: • Угол между двумя векторами
Смотри ещё: • Пересечение двух отрезков

  • 1
Вот ссылка на алгоритм попроще..

http://steps3d.narod.ru/snippets.html#ptinpoly

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

Угу. Там и у меня одинаковая математика, только там ради того что бы засунуть в одну строчку - повторно делаются сравнения и лишнее деления. Насчет "проще" не согласен :)

На то что я изобрёл что-то новое - не претендую, просто слишком часто вижу кривой код.


Винни слишком молод ещё, чтобы скатываться в маразм и переизобретать известные алгоритмы. Но это не мешает ему о них писать.

в 57-ой школе я переизобрёл матан *^_^*

А я в детском саду переизобрел электродвигатель. Только я щетки не изобрел. У меня ротор был с постоянным магнитом и нужно было специальную ручку крутить, чтобы менять полярность магниитов статора. Короче, недодумал конструкцию.

(Deleted comment)
(Deleted comment)
Вроде работает, а что не так?

#include <stdio.h>


int ptInPoly ( 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;
}

int main()
{
	const int N = 4;

	float xs[4]= { 10, 30, 50, 30 };
	float ys[4]= { 10, 3, 10, 20 };
	
	
	for (int y = 0; y < 30; ++y)
	{
		for (int x = 0; x < 78; ++x)
		{
			if (ptInPoly(N, xs, ys, x, y))
				printf("*");
			else
				printf("_");
		}
		printf("\n");
	}
}

(Deleted comment)
Код написан не очень красиво/понятно. Поэтому за '&&' и за '||' не сразу видно, деления на ноль тут не будет. Дело в том, что если (a < x <= b), то деление (x - a)/(b-a) вполне корректно.

(Deleted comment)
в интернетах численно-стабильных решений этой задачи много меньше, чем решений вообще, т.ч. пост не бесполезен. к тому же автор на новизну идеи и не претендует, насколько я понимаю.

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

А с вертикальными и горизонтальными отрезками тут все ОК?

Да. Если у концов отрезка одинаковые y, то он будет или целиком под лучом, или над лучом, не влияя на четность количества пересечений.

Спасибо.

Определитель получается из уравнения полуплоскости :)

20 лет спустя :)

(Anonymous)
Отличный и очень полезный пост, спасибо!

Может никто уже и не прочитает это, но все-таки решусь написать.

Мне кажется, что этот конкретный код работает с лучом, направленным влево, а не вправо.
Нам, конечно же, все равно куда направлять - на результат это никак не влияет.
Однако же текст и иллюстрации говорят о луче вправо, а код работает с лучом влево :)

Собственно, это и вызвало у меня небольшой диссонанс и желание написать комментарий.
Просто вдруг я чего-то не "догнал" :)

Конкретный пример, на котором можно рассмотреть - такое ребро (для простоты допуская, что дан point = (0; 0)):
(5; -5) -> (5; 1)
Пересечение не будет засчитано, т.к.
t = (5 * (1 - (-5)) - (-5) * (5 - 5)) = 5 * 6 = 30
и выходит, что t > 0, в то время как истинно второе условие (!cur_under && prev_under).

Для зеркального же случая (-5; 5) -> (-5; -1) пересечение засчитано будет, т.к. t = 30 по-прежнему, но теперь истинно первое условие (cur_under && !prev_under).

Спасибо, дружище за старый, но нужный пост =)

А всегда л работает?

(Anonymous)
Здравствуйте.

Проведу небольшую эксгумацию двухлетней, но все еще актуальной темы.

А что если луч от точки, лежащей за границей полигона войдет в него через вершину? Да потом еще и пересечет вершину, которая нарушает выпуклость?
Как отреагирует алгоритм?

PS Я долго искал решение, но остановлся на подсчете суммы углов между отрезками, соедияющими точку и вершины. Если 0, то снаружи, 180-внутри.

мой код.

bool DXF_Area::Is_Inside(Point Point1)
{
if(Lines->Count<3) return false; //check polygon

if(Point1.X>Max_X && Point1.Y>Max_Y) return false; //check limits
if(Point1.X
[Error: Irreparable invalid markup ('<min_x [...] dxf_line^>') in entry. Owner must fix manually. Raw contents below.]

Здравствуйте.

Проведу небольшую эксгумацию двухлетней, но все еще актуальной темы.

А что если луч от точки, лежащей за границей полигона войдет в него через вершину? Да потом еще и пересечет вершину, которая нарушает выпуклость?
Как отреагирует алгоритм?

PS Я долго искал решение, но остановлся на подсчете суммы углов между отрезками, соедияющими точку и вершины. Если 0, то снаружи, 180-внутри.

мой код.

bool DXF_Area::Is_Inside(Point Point1)
{
if(Lines->Count<3) return false; //check polygon

if(Point1.X>Max_X && Point1.Y>Max_Y) return false; //check limits
if(Point1.X<Min_X && Point1.Y<Min_Y) return false;

double Result_Angle=0;
double Angle1,Angle2;
double Current_Dalta;

DXF_Line^ Checking_Line=gcnew DXF_Line(); //line for gey angle

Checking_Line->X1=Point1.X; //set first point
Checking_Line->Y1=Point1.Y;

Checking_Line->X2=Lines[Lines->Count-1]->X2; //set first point to end of polygon
Checking_Line->Y2=Lines[Lines->Count-1]->Y2;
Angle1=Checking_Line->Angle; //get angle
if(Angle1<0) Angle1+=(Math::PI*2); //angle shoud be positive

for(int i=0;i!=Lines->Count;i++)
{
Checking_Line->X2=Lines[i]->X1; //set next point for check
Checking_Line->Y2=Lines[i]->Y1;

if(Checking_Line->X1==Checking_Line->X2 && Checking_Line->Y1==Checking_Line->Y2) return true; //if point is on one of polygon points return 3

Angle2=Checking_Line->Angle; //get new angle
if(Angle2<0) Angle2+=(Math::PI*2); //angle shoud be positive

if(Math::Abs(Math::Abs(Angle1-Angle2)-Math::PI)<ANGLE_ACCURACY) return true; //if point is on this intercest return 2

Current_Dalta=Angle1-Angle2; //calculating delta

if(Current_Dalta>Math::PI) Current_Dalta-=Math::PI; //take shitrest angle
if(Current_Dalta<-Math::PI) Current_Dalta+=Math::PI;

Result_Angle+=Current_Dalta; //store delta in result
Angle1=Angle2; //store point for next step
}

Checking_Line->X2=Lines[Lines->Count-1]->X2; //set last point
Checking_Line->Y2=Lines[Lines->Count-1]->Y2;
Angle2=Checking_Line->Angle; //get angle
if(Angle2<0) Angle2+=(Math::PI*2); //angle shoud be positive

Current_Dalta=Angle1-Angle2; //calculating delta

if(Current_Dalta>Math::PI) Current_Dalta-=Math::PI; //take shitrest angle
if(Current_Dalta<-Math::PI) Current_Dalta+=Math::PI;

Result_Angle+=Current_Dalta; //store delta in result

if(Math::Abs(Result_Angle)<ANGLE_ACCURACY) return false; //point is outside return 0

if(Math::Abs(Result_Angle-Math::PI)<ANGLE_ACCURACY) return true; //point is inside return 1

throw 1; //unknown result throw?
}

Re: А всегда л работает?

То, что луч проходит через вершину - корректно учтено, так как даже в таком случае логическое выражение (точка ниже луча) считается одинаково, и пересечение засчитываются только с тем отрезком, для которого это логическоне значение поменялось.

Наверное можно и через угол, только лучше вермишель из 'if' внутри цикла заменить на atan2(axby - bxay, axbx + ayby) (см. http://users.livejournal.com/_winnie/237714.html )
Замечу ещё, что решение через угол будет давать другой результат, если вместо полигона - спираль, которая делает 2 оборота вокруг точки.



По коду - сильно удивлён использованию специфичных расширений Microsoft-компилятора в математическом коде (gcnew), удивлён что точки - это ОБЪЕКТЫ с доступом через указатель (а не просто пара double), ну и throw 1. Хорошо бы прочитать книжку про C++, например Страуструпа, и по чистому C.




Edited at 2014-01-13 09:20 am (UTC)

Re: А всегда ли работает?

Перед предыдущим постом забыл залогиться. Сори.

Хорошо. А как считается прохождение луча через вершину, которая нарушает выпуклость? считатся 2 раза или 1? Поменяетися ли четность?
Я подставил Ваш код к себе в программу и он работает. Но вопрос стоит тот же.

У меня нет спиралей, как и пересечений. Работаю с реальными объектами. Но если проверять на кратность пи то должно работать.


По коду:
gcnew-можно убрать и сделать так DXF_Line^ Checking_Line. Тогда Checking_Line->X1 будет Checking_Line.X1. Есть принципиальная разница?

ОБЪЕКТЫ- для получения угла и простоты писания и использования.Тоесть отрезок это 2 точки а угол свойство. Для рисования просто передаю Point в функцию рисования. Внутри и есть пара double. Угол это property функция.

throw 1- из за недописанности кода. Пока все сыро.

Спасибо за критику.

Edited at 2014-01-13 01:12 pm (UTC)

Re: А всегда ли работает?

Четность корректно меняется. Или переключается туда и обратно 2 раза, или переключается 0 раз (вообще не переключается), вот как раз этот случай изображен слева:
http://dobrokot.ru/pics/2010-06-15__22-24-30_3kb.png



Edited at 2014-01-13 02:40 pm (UTC)

Учет попадание точки на контур

(Anonymous)
По идее чтобы точка на контуре тоже засчитывалась, надо везде < и > заменить на <= и >=. Но всё равно точки на гранях не учитываются, не соображу почему...,?

Re: Учет попадание точки на контур

Случай "точки прямо на контуре" в данном случае не обрабатывается, так как в компьютере не существует контура в математическом смысле, из-за неточности хранения числа в float/double.

Для этого нужна логика не "да/нет", а "да/нет/точка на расстоянии миллиметра от контура".


Edited at 2015-04-22 09:40 pm (UTC)

Re: Учет попадание точки на контур

(Anonymous)
Сейчас, если я правильно понимаю, точка считается не в контуре, если она не в области контура. А можно поменять логику на такую - "если точка не снаружи контура, то она в контуре"? Тогда будет учтено попадание на границу.

Re: Учет попадание точки на контур

Об этом есть смысл думать, если вместо double/float используется точная арифметика. Если используется хранение чисел в double/float, то замена < на <= ни в каком месте не поможет: если в результате ошибок округления получилось 1.9999998, то неважно сравнение с 2 записано как <= или <.

Если же используется точная арифметика, то во-первых надо поменять тип t, во-вторых думать про сравнения (t > 0) и (t < 0).

1) Есть подозрение, что нестрогое неравенство нужно только в одном месте, а не в обоих.

2) В каком именнно месте - зависит, обходим мы многоугольник по часовой стрелке, или против (если пойти в обратную сторону, то знаки всех переменных поменяются, и от этого поменяются местами "внешность" и "внутренность").

Re: Учет попадание точки на контур

(Anonymous)
Да, если данные из реального мира, то вероятность что точка попадет на границу очень мала. А теперь представим такую ситуацию - два квадрата, один большой и один маленький
{
{0.0, 0.0},
{0.0, 100.0},
{100.0, 100.0},
{100.0, 0.0}
}
{
{100.0, 25.0},
{100.0, 75.0},
{150.0, 75.0},
{150.0, 25.0}
}
Они касаются. Но если я возьму середину границы первого прямоугольника - точку {100.0,50.0}, алгоритм выдаст что она снаружи второго прямоугольника -> т.е. прямоугольники не пересекаются и необъединяемы и это не зависит от типа, хоть возьми long double, хоть int.
Предложенный вами алгоритм работает хорошо, но в предложенной мной задаче его логика не подходит. Я не говорю что надо как-то там отдельно учитывать границу, просто тут может быть два варианта как учитывать точку, каждый из которых применим в отдельных случаях:
1. Точка не попадает в область контура => она снаружи (это то, что у вас)
2. Точка не попадает в область "наружи" => она в контуре.

Re: Учет попадание точки на контур

Скорее всего надо смотреть на задачу целиком. Математическая она или практическая (про 2d-3d-моделирование реальных/игровых объектов).

Если у соседних многоугольников вершины - общие, то можно решить задачу "отнести точку к одному или другому многоугольников без пропусков и наложений", и решать её без epsilon-допусков.

Если у нас все многоугольнки - это квадраты и прямоугольники с вертикальными-горизонтальными линиями, и поэтому точное решение теоретически возможно даже в float/double числах - то нужно решать её по-другому, а не через обобщенную функцию для любых многоугольников.

Если у соседних многоугольников вершины не совпадают ( http://en.wikipedia.org/wiki/T-vertices ) - то задача в принципе не решаема без введения допусков вида "щель тоньше 10 микрометров не считается щелью".






Edited at 2015-05-07 02:28 pm (UTC)

Re: Учет попадание точки на контур

Если предположить, что у нас арифметика - точная, то алгоритм наверное надо изменить так, чтобы для многоугольников против часовой стрелки контур был внутри:

        if (cur_under && !prev_under)
        {
            if (t >= 0) // изменение только тут
                intersections_num += 1;
        }
        if (!cur_under && prev_under)
        {
            if (t < 0) // а тут как раньше.
                intersections_num += 1;
        }


Для многоугольников по часовой стрелке - ровно наоборот.

Ещё, нужно дополнительно обработать случай горизонтальных отрезков (сейчас для горизонтальных отрезков никогда не выполняется условие prev_under != cur_under)

Только это может не решить проблему, так как арифметика неточная, и баг для искусственных квадратов [0,100]x[0,100] починится, а для наклонных квадратов - нет.

Я очень скептично отношусь к тому, что в какой-то задаче полезно точное решение задачи про точку на контуре, но это голословное заявление.

Можно решать полезную задачу "отнестить точку ровно к одному треугольнику" для полигональных моделей, где к стороне каждого треугольника модели примыкает сторона другого треугольника, но если у этих сторон не совпадают в точности две вершины - то задача снова ломается от неточности float.


Edited at 2015-05-07 03:00 pm (UTC)

  • 1
?

Log in