Задача о принадлежности точки фигуре

Есть набор N точек (x1,y1) (x2,y2) .... (xn,yn), которыми задан (не обязательно) выпуклый многоугольник.

Как быстрее определить, принадлежит ли этому многоугольнику заданная точка (x,y)?

(нужно на Java для генерации распределения мобов по локациям на игровом сервере )

Может (условно) разбить выпуклый многоугольник на треугольники с некой общей вершиной внутри? Если не совсем выпуклый, то сперва разбить на выпуклые, а потом на треугольники.

Меня интересует конкретный алгоритм Поскольку задача не настолько насущная, чтобы садиться за его написание собственноручно, но и не настолько ненужная, чтобы её игнорировать

Таки пришлось самому решать.

Придумал такой алгоритм.

Рассматриваем в цикле все отрезки.

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

Если число пересечений с отрезками чётное (или нулевое) - то точка фигуре не принадлежит. Если нечётное - то принадлежит.

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

Упрощения:

• Если ΔYн и ΔYк имеют один знак, то с отрезком не пересекаемся.
• Если знаки и ΔXн и ΔXк > 0 - пересекаемся
• Если знаки и ΔXн и ΔXк < 0 - не пересекаемся

Остаётся один случай, когда начальная и конечная точки отрезка лежат по разные стороны от прямой, на которой находится наш луч, и по разные стороны от прямой, нормальной к нашему лучу.

Тогда считаем ΔXнп = Xп - Xн точки пересечения прямой с лучом и отрезка:
ΔXнп = ΔYнп*(ΔX/ΔY), где ΔYнп = Yл0 - Yн

Если ΔXнп > ΔXн, то не пересекается. Если меньше - пересекается.

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

Вот, если кому нужен рабочий (кажется) пример на Java. Редкий случай, когда программа заработала сразу (понятно, синтаксические ошибки не считаются)

1. //package net.sf.l2j.gameserver.model;
2.  
3. import java.util.ArrayList;
4.  
5. public class L2Area
6. {
7.         protected class Point
8.         {
9.                 protected double x, y;
10.                 Point (double _x, double _y) {x=_x; y=_y;}
11.         }
12.  
13.         private ArrayList<Point> _points;
14.  
15.         L2Area()
16.         {
17.                 _points = new ArrayList<Point>();
18.         }
19.  
20.         public void add(double x, double y)
21.         {
22.                 _points.add(new Point(x,y));
23.         }
24.  
25.         public void print()
26.         {
27.                 for(Point p : _points)
28.                         System.out.println("("+p.x+","+p.y+")");
29.         }
30.  
31.         public boolean isIntersect(double x, double y, Point p1, Point p2)
32.         {
33.                 double dy1 = p1.y - y;
34.                 double dy2 = p2.y - y;
35.  
36.                 if(Math.signum(dy1) == Math.signum(dy2))
37.                         return false;
38.                
39.                 double dx1 = p1.x - x;
40.                 double dx2 = p2.x - x;
41.  
42.                 if(dx1 >= 0 && dx2 >= 0)
43.                         return true;
44.                
45.                 if(dx1 < 0 && dx2 < 0)
46.                         return false;
47.  
48.                 double dx0 = dy1 * (p1.x-p2.x)/(p1.y-p2.y);
49.  
50.                 return dx0 <= dx1;
51.         }
52.  
53.         public boolean isInside(double x, double y)
54.         {
55.                 int intersect_count = 0;
56.                 for(int i=0; i<_points.size(); i++)
57.                 {
58.                         Point p1 = _points.get(i>0 ? i-1 : _points.size()-1);
59.                         Point p2 = _points.get(i);
60.  
61.                         System.out.println("("+p1.x+","+p1.y+")-("+p2.x+","+p2.y+") => "+isIntersect(x,y,p1,p2));
62.                        
63.                         if(isIntersect(x,y,p1,p2))
64.                                 intersect_count++;
65.                 }
66.  
67.                 return intersect_count%2 == 1;
68.         }
69.  
70.         public void test(L2Area a)
71.         {
72.  
73.         a.add(0,5);
74.         a.add(10,0);
75.         a.add(0,0);
76.  
77.                 System.out.println(a.isInside(1,2.5));
78.         }
79.  
80.     public static void main(String[] args)
81.     {
82.         L2Area a = new L2Area();
83.  
84.         a.test(a);
85.         }
86. }

Рома! Как не стыдно! Это же школьная задачка! Решается очень просто: последовательно суммируешь все углы из твоей точки до каждой вершины (т.е., если назвать твою точку 0, угол 1-0-2, 2-0-3 ... n-0-1), учитывая знак (скажем, против часовой - плюс, по - минус). Если точка внутри фигуры, сумма углов будет 2пи, иначе 0.

Рабочий (кажется ) пример на MATLAB:

1. function r = belong(A, p)
2. % A - массив 2хN, описывающий фигуру; p - точка (вектор длиной 2)
3. A(1,:) = A(1,:) - p(1); % вычитаем из фигуры координаты точки, чтобы работать с векторами.
4. A(2,:) = A(2,:) - p(2); % Теперь определяем как бы принадлежность начала координат
5. A = [A A(:,1)]; % добавляем в конец начальную точку
6. r = 0;
7. for i = 1:size(A,2)-1 % идем по всем точкам
8.     r = r + acos( dot(A(:,i), A(:,i+1)) / (norm(A(:,i))*norm(A(:,i+1))) ) ...
9.         * sign(A(1,i)*A(2,i+1) - A(2,i)*A(1,i+1));
10. % тут есть хитрость: чтобы определять углы >90°, сам угол вычисляется через скалярное произведение (dot)
11. % но для определения знака поворота используется знак векторного произведения.
12. % norm - норма (евклидова длина) вектора
13. % hint: скалярное произведение можно вычислить как x1*x2+y1*y2, если встроенного dot product нет
14. end
15. r = r > 100*eps; % eps - встроенная константа ~2E-16
16. % формально надо порядка N*eps, но вообще неважно, тут либо около 0 будет, либо 2*pi.

P.S. Я не делал проверку на сингулярность, которая возникает, когда тестовая точка совпадает с одной из вершин. Тут еще надо уточнить, чем считать такой случай Сам Матлаб нормально работает с делением на 0, возвращая NaN, так что добавить лишь придется if isnan(r) ... после цикла. Если точка лежит на ребре, в сумме будет pi - это тоже можно обработать при желании.

Zeus> Там вся тригонометрия - арккосинус

Неужели такие ЧЛ бывают? Да ещё от компьютеров??

>Зато сам алгоритм на порядок проще и легче проверяется (даже чисто математически).

Да ну, и чем же он проще? Только тем, что строчек меньше? Так у меня 90% кода - объявления классов. Это ж дело потом расширяться будет сильно, факторизовать нужно изначально...

Zeus>> Там вся тригонометрия - арккосинус
Balancer> Неужели такие ЧЛ бывают? Да ещё от компьютеров??

Хм, это комплимент или наоборот?

Balancer> Да ну, и чем же он проще? Только тем, что строчек меньше? Так у меня 90% кода - объявления классов.

Ну а у меня - 70% комментариев Нет, конечно, не напорядок, да и сам твой алгоритм известен (во всяком случае я припоминаю подсчет пересечений, просто само у меня это не всплыло); но с углами проще и по формулировке, и по обоснованию, и по отладке (меньше особенных ситуаций, на которые надо обращать внимание). Да и просто никакой "доморощенности", ничего не изобретается, нет логических ветвлений...

Нет, тригонометрию нам не нужно. До тех пор, пока плавучка не станет более экономной, чем целочисленка. А этого не случится никогда

Вообще, выбор более простого алгоритма (кстати, всё равно момент спорный. ИМХО, мой алгоритм нагляднее и проще - но это уже могут быть особенности информационного метаболизма ) перед более эффективным почти всегда есть признак непрофессионализма

Может, ты и линии будешь по тригонометрии рисовать, а не по Брэзенхему? Или сортировать пузырьком? Или множители за пределы цикла не будешь выносить?

Zeus> Рома! Как не стыдно!

А если точка в центре кольца?

Интересно, насколько живуч миф о больших потерях в вычислениях с ПТ...

Я не уверен, что acos тут будет выгоднее, но то что потери от него относительно невелики - это точно. На fpu первого пента acos - это где-то 100-300 тактов... С тех пор ситуация с fpu, вроде, только лучшела, а штрафы за плохой переход - только росли...

Да и ява...

Татарин> Я не уверен, что acos тут будет выгоднее, но то что потери от него относительно невелики - это точно. На fpu первого пента acos - это где-то 100-300 тактов...

Угу. А у нас массив на 500Мб координатных данных. Хранить всё это в плавучке - это значит потратить уже не 500Мб, а гиг с четвертью. Значит - храним в целочисленке. Да и глупо хранить изначально целочисленные координаты в плавучке. Значит, при каждой проверке - эти целочисленные данные постоянно надо грузить в сопроцессор. Потом - выгружать. Это - море штрафов, паразитные ожидания... Да и сами вычисления. Пусть acos хоть за 10 тактов считается. Всё равно это на порядок с лишним медленнее сравнений/вычитаний. А с паразитными расходами - там и два порядка набежит. И вдвое более широкий поток данных.

Короче, чтобы посчтитать на коленке в матлабе - тригонометрия годится.

Но не для нагруженного сервера...

...

Зато, глядя на этот топик, я теперь понимаю, почему нынешний софт так тормозит

Переведите кто-нить Зевсовский пример во что-то более низкоуровневое, а то я Матлаба не знаю, а лепить алгоритм с нуля ломает Измерим реальную скорость того и другого метода. Я бенч-оболочку уже слепил Целочисленный вариант проверяет у меня 10млн. точек на попадение в треугольник, занимающий 50% тестовой области за 6.8сек.

Ром, а как ты проверяешь, если луч пересекается с вершиной? Это ведь может быть случай и пересечения с контуром, и непересечения - в зависимости от того, острый угол или нет...

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

И по поводу Зевсового алгоритма...

На самом деле надо просто не через суммирование углов считать, а через суммирование синус-косинусов. По сути получится сложная сумма разниц координат... вполне целочисленно.

Хотя я бы считал через расчет попадания точки в отрезок :

Сначала нормируем все отрезки - то есть приводим их к виду Xн<Xк меняя концы.

Потом фильтруем по X на попадание точки в полосу отрезка.
Потом по Y - на попадание в прямоугольную область с диагональю отрезка.

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

Далее, надо рассматривать соседние отрезки. Если от концов отрезка соседние отрезки идут в разные стороны по вертикали - то засчитываем пересечение. Если нет - то не засчитываем - луч только 'чиркает' по фигуре и ее не пересекает...

Хотел повозмущаться попунктно, но потом увидел предложение забенчить.
ИМХО, показательнее и справедливее некуда.

Чуть попозже. В чем пишем? Ява?

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

Тригонометрической - нет.
Вот считать соседние отрезки - да, придется.