Проверка принадлежности точки невыпуклому многоугольнику

Проверить принадлежность точки невыпуклому многоугольнику за линейное время совсем не сложно. Один из самых распространенных методов — выпустить луч и посчитать число точек пересечения. Однако, при этом нужно аккуратно рассматривать случаи, когда точки многоугольника попадают на луч. Отсюда естественно возникает вопрос, как рассмотреть эти случаи проще всего?

Формулировка задачи

Будем решать вариант задачи, в котором вершины многоугольника и проверяемая точка заданы целочисленными координатами до 10⁸. На вершины многоугольника, для отсутствия простоты, не наложим никаких ограничений. Разрешим совпадающие вершины, вершины на сторонах, самокасания и нулевую площадь. Запретим самопересечения, иначе не возможно определить внутренность получившейся ломаной. Программа должна выдавать 0, если точка лежит на границе многоугольника, в противном случае нужно выдавать 1, если точка лежит вне многоугольника и -1, если точка лежит внутри многоугольника.

Общая схема решения

Пусть x₁, x₂, …, x_n — вершины многоугольника (они могут быть заданы как по, так и против часовой стрелки), y — точка, которую нам нужно проверить. Если n=1, то достаточно проверить равенство y и x₁. В противном случае посчитаем выражение F(y)=f(x₁,x₂,y)·f(x₂,x₃,y)···f(x_n-1,x_n,y)·f(x_n,x₁,y). Функция f(a,b,c) принимает значения 1, 0, -1. Она определяет пересекает ли луч выпущенный из точки c вдоль оси y в направлении уменьшения x отрезок (a,b) с учетом всех крайних случаев. В случае пересечения она будет возвращать -1. Если точка c лежит на отрезке (a,b) она вернет 0, в противном случае она будет возвращать 1.

Оптимизируем функцию f

Вся сложность проверки принадлежности точки многоугольнику сосредоточена в функции f. Поэтому у меня возникла идея перебрать ее реализации и найти среди них самую короткую. Понятно, что первый шаг в функции f — это вычисление координат a и b относительно точки c. Занесем эти координаты в локальные переменные ax, ay, bx и by. Для проверки реализаций, я построил все возможные отрезки (a,b) с координатами не больше 5 и нашел для них правильные ответы. Кроме того, выражения, которые можно использовать в f, я ограничил величинами signum(ax * by - ay * bx), signum(ax), signum(bx), signum(ay), signum(by). Таким образом, функция f должна выдать ответ, основываясь на значении этих 5 величин, которые могут принимать значения -1, 0 или 1.
В начале, я попытался выразить ответ как знак линейной комбинации этих величин. К сожалению, из этого ничего не вышло. Следующая идея — использовать цепочку вложенных if’ов. В качестве сложности программы выберем количество if’ов + количество return’ов. В качестве условий if’ов разрешим выражения вида v_i in S, где v_i — одно из 5 signum выражений выше, а S — подмножество {-1, 0, 1}.
Теперь можно запустить перебор программ. В качестве входа функция перебора принимает множество разрешенных значений для каждого из 5 переменных и набор номеров тестов, которые этим условиям удовлетворяют. В качестве результата, функция перебора выдает дерево функции с минимальной сложностью, которая проходит все тесты. Для ускорения перебора, все результаты запоминаются и не рассчитываются повторно.
Количество вариантов оказалось достаточно мало, чтобы полный перебор программ завершался за несколько секунд. В результате, в качестве оптимальной получилась программа со сложностью 27. После того, как я разрешил возвращать кроме констант выражения вида v_i и —v_i, сложность упала до 18. Это все еще слишком много. Внимательное изучение программы позволило обнаружить xor-like структуру при проверке знаков ax и bx. Поэтому я добавил переменную signum(ax * bx) и сложность упала до 13. Эксперимент с добавлением переменной signum(ay * by) уменьшил сложность до 11. В результате получилась следующая функция:

public static int check2(Point a, Point b, Point middle) {     long ax = a.x - middle.x;     long ay = a.y - middle.y;     long bx = b.x - middle.x;     long by = b.y - middle.y;     if (ay * by > 0)         return 1;     int s = Long.signum(ax * by - ay * bx);     if (s == 0)     {         if (ax * bx <= 0)             return 0;         return 1;     }     if (ay < 0)         return -s;     if (by < 0)         return s;     return 1; }

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

1. В начале проверяем не лежит ли отрезок (a,b) строго по одну сторону от луча. В этом случае пересечения точно нет.
2. Проверяем вырожденный случай, когда точка c лежит на прямой (a,b). Ответа ноль можно избежать только если отрезок лежит точно на луче и обе точки находятся по одну сторону.
3. Xor-like код для проверки пересечения луча с отрезком. Знак у s зависит от того, какая точка ниже луча. Здесь можно поставить минус у s в любом из двух случаев; знак меньше, заменить на больше или меньше или равно.
4. Возвращаем 1 на случай, когда одна точка лежит на луче, а другая находится по некоторую фиксированную сторону от него. Такие пересечения не засчитываются, чтобы не посчитать одно пересечение два раза, если вершина многоугольника попала на луч.

Но и это не предел! В этой реализации ay умножается на by. Почему бы не избавиться от умножения, заменив его на что нибудь вроде ay < 0 ^ by < 0? В результате получится не совсем то же самое, но можно добавить эту переменную в перебор и посмотреть что получится. Очередной запуск алгоритма… и он выдает решение со сложностью 9! Новый вариант выглядит так:

public static int check3(Point a, Point b, Point middle) {     long ax = a.x - middle.x;     long ay = a.y - middle.y;     long bx = b.x - middle.x;     long by = b.y - middle.y;     int s = Long.signum(ax * by - ay * bx);         if (ay < 0 ^ by < 0)     {         if (by < 0)             return s;         return -s;     }     if (s == 0 && ax * bx <= 0)         return Long.signum(ay * by);          return 1; } 

Его можно еще немного облагородить. После проверки на ay < 0 ^ by < 0 величина Long.signum(ay * by) не может равняться -1. Поэтому последние строчки можно переписать так:

public static int check3(Point a, Point b, Point middle) {     long ax = a.x - middle.x;     long ay = a.y - middle.y;     long bx = b.x - middle.x;     long by = b.y - middle.y;     int s = Long.signum(ax * by - ay * bx);         if (ay < 0 ^ by < 0)     {         if (by < 0)             return s;         return -s;     }     if (s == 0 && (ay == 0 || by == 0) && ax * bx <= 0)         return 0;          return 1; } 

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

public static int check3(Point a, Point b, Point middle) {     long ax = a.x - middle.x;     long ay = a.y - middle.y;     long bx = b.x - middle.x;     long by = b.y - middle.y;     int s = Long.signum(ax * by - ay * bx);         if (s == 0 && (ay == 0 || by == 0) && ax * bx <= 0)         return 0;       if (ay < 0 ^ by < 0)     {         if (by < 0)             return s;         return -s;     }           return 1; } 

В новом алгоритме опять же есть большое число степеней свободы. В выражении ay < 0 ^ by < 0, вместо знака меньше, можно использовать знак больше или нестрогий знак. Знак минус у s опять же можно ставить у любого из двух случаев. Переменную by можно поменять во вложенном if’e на ay.
Сложнее всего в новом алгоритме понять случай нуля. В нем одновременно проверяются два независимых случая. В первом случае одна из точек отрезка (a,b) совпадает с точкой c, во втором случае отрезок параллелен оси x и проходит через точку c. Эти проверки можно выполнить и другими способами, но способ выше — самый короткий.

Тестирование на скорости работы в худшем случае показывает, что на 32-битной JVM функция check2 быстрее check3. На 64-битной JVM результат не ясный. Однако, разница в обоих случаях довольно мала — 5-10%. Попытки замены ax * bx <= 0 неравенствами почему-то не дают заметных улучшений ни в 32-битной, ни в 64-битной версии java машины.

Исходники для тех, кто хочет поэкспериментировать: PolyCheck.zip