Пользователь
0,0
рейтинг
6 декабря 2012 в 10:35

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

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

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


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

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


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

1 and -1 cases

Оптимизируем функцию 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'ов разрешим выражения вида vi in S, где vi — одно из 5 signum выражений выше, а S — подмножество {-1, 0, 1}.
Теперь можно запустить перебор программ. В качестве входа функция перебора принимает множество разрешенных значений для каждого из 5 переменных и набор номеров тестов, которые этим условиям удовлетворяют. В качестве результата, функция перебора выдает дерево функции с минимальной сложностью, которая проходит все тесты. Для ускорения перебора, все результаты запоминаются и не рассчитываются повторно.
Количество вариантов оказалось достаточно мало, чтобы полный перебор программ завершался за несколько секунд. В результате, в качестве оптимальной получилась программа со сложностью 27. После того, как я разрешил возвращать кроме констант выражения вида vi и -vi, сложность упала до 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). Ответа ноль можно избежать только если точки лежат по одну сторону от оси y. Если они обе лежат на оси y, то случай когда 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
@halyavin
карма
46,0
рейтинг 0,0
Реклама помогает поддерживать и развивать наши сервисы

Подробнее
Спецпроект

Самое читаемое Разработка

Комментарии (57)

  • +13
    Может стоит добавить иллюстраций в статью? Без них читать не интересно, проще сразу скачать исходники.
    • +5
      Да-да! Картинки нужны для привлечения внимания, факт!
      • +1
        Добавил картинку. Надеюсь я не перепутал никаких случаев…
        • 0
          Теперь стало ещё непонятнее.
          Где собственно пересечения луча «вдоль y в сторону уменьшения x» с отрезками?
  • +1
    Если бы вы привели и другие алгоритмы проверки принадлежности, например разбивка невыпуклого многоугольника на выпуклые многоугольники или треугольники с последующей проверкой принадлежности к ним, то было бы интересней и полезней.
    • +1
      Ваш способ на порядок быстрее ( если разбивка провелась заранее, что в жизни и делают для топологически статических полигонов)
    • 0
      Алгоритм с разбиением хоть и возможен в принципе, но смысла в нем мало (на мой взгляд). А вот про алгоритм с подсчетом числа оборотов стоило упомянуть, так как он дает немного другой резльтат (для самопересекающихся многоугольников), что иногда и требуется. Очень нравится, как эти два алгоритма реализованы одновременно в Qt: qt.gitorious.org/qt/qt/blobs/4.7/src/gui/painting/qpolygon.cpp#line855
      • 0
        Ну почему же мало. Мы, например, в программе все многоугольники храним в виде набора выпуклых. Это ускоряет работу с геометрией в разы, а разбивка производится всего один раз.
  • +1
    ==
    Она определяет пересекает ли луч выпущенный из точки c вдоль оси y в направлении уменьшения x отрезок (a,b)
    ==
    Не очень понятна фраза «в направлении уменьшения x». Если луч идет вдоль оси y (параллелен ей), то x при этом не изменяется совсем. Или я просто неправильно себе все представляю?
    • 0
      Моя ошибка. Исправил + привел в соответствие с программами.
  • 0
    Ваши функции поздно отсекают очевидные случаи. Во многих случаях умножать не надо вообще, а у вас всегда минимум два умножения. К примеру, рассмотрим случай, когда ay и by оба положительны или оба отрицательны (на больших многоугольниках и случайных расположениях точки middle это может выполняться для большинства рёбер). Тогда ни один if не сработает, и вы вернёте 1, но при этом выполните два абсолютно лишних 64-битных умножения.
  • 0
    Вот такой вариант потестируйте, если я нигде не ошибся:
    public static int check3_mod(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 < 0 ^ by < 0)
        {
            int s = Long.signum(ax * by - ay * bx);
            if (by < 0)
                return s;
            return -s;
        }
        if ((ay == 0 || by == 0) && ax * bx <= 0 && (ax * by - ay * bx) == 0)
            return 0;     
        return 1;
    }
    • 0
      На 64-битной JVM работает медленнее всех на полностью случайном тесте. Сейчас попробую еще на 32-битной проверить.
      • 0
        У меня на 32bit
        rand
        Method 11 = 239818469
        Method 9, bad first = 316188510
        Method 9, bad last = 321210656
        Method 9, bad first &| = 312058097
        Method 9, bad first <=0 = 310681106
        Method 9, bad first min max = 300902768
        No extra mult = 251040286
        randpos
        Method 11 = 149962635
        Method 9, bad first = 279716835
        Method 9, bad last = 293585358
        Method 9, bad first &| = 278011033
        Method 9, bad first <=0 = 279741699
        Method 9, bad first min max = 278360518
        No extra mult = 138307522

        Похоже, 64-битные умножения на 64-битной машине не просто быстрее, а существенно быстрее :-)
    • 0
      На 32-битной JVM ваш метод быстрее check3 и чуть медленнее (на уровне погрешности измерений) check2.
      • 0
        Ну правильно, в check2 как раз эти случаи сразу и отсекаются первым if. Я больше check3 критиковал, а на check2 не глянул сразу.
  • 0
    Реабилитируюсь. Вот вам метод, который значимо быстрее всех на 32bit:
        public static int check7(Point a, Point b, Point middle)
        {
            if( ( a.y > middle.y ) ^ ( b.y <= middle.y ) )
                return 1;
            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)
            {
                if (ax * bx <= 0)
                    return 0;
                return 1;
            }
            if (ay < 0)
                return -s;
            if (by < 0)
                return s;
            return 1;
        }

    Экономия за счёт того, что во многих случаях int не превращаются в long вообще.
    • 0
      dry run
      Method 11 = 291993256
      Method 9, bad first = 304996864
      Method 9, bad last = 323442504
      Method 9, bad first &| = 317873082
      Method 9, bad first <=0 = 318828510
      Method 9, bad first min max = 346162152
      No extra mult = 98121688
      y=0
      Method 11 = 288414309
      Method 9, bad first = 305893626
      Method 9, bad last = 323014797
      Method 9, bad first &| = 320161362
      Method 9, bad first <=0 = 319378859
      Method 9, bad first min max = 346463028
      No extra mult = 96438793
      x=0
      Method 11 = 148143968
      Method 9, bad first = 293686208
      Method 9, bad last = 309175608
      Method 9, bad first &| = 306875315
      Method 9, bad first <=0 = 301019263
      Method 9, bad first min max = 303745308
      No extra mult = 97887302
      y+-
      Method 11 = 269566384
      Method 9, bad first = 288282729
      Method 9, bad last = 292221497
      Method 9, bad first &| = 288991757
      Method 9, bad first <=0 = 282111273
      Method 9, bad first min max = 274667311
      No extra mult = 256442927
      y+0
      Method 11 = 276192924
      Method 9, bad first = 276658905
      Method 9, bad last = 282801864
      Method 9, bad first &| = 281113661
      Method 9, bad first <=0 = 272635768
      Method 9, bad first min max = 276043184
      No extra mult = 254776235
      rand
      Method 11 = 240899053
      Method 9, bad first = 318759227
      Method 9, bad last = 320744396
      Method 9, bad first &| = 312402274
      Method 9, bad first <=0 = 308858528
      Method 9, bad first min max = 300334819
      No extra mult = 204082362
      randpos
      Method 11 = 147626583
      Method 9, bad first = 279747845
      Method 9, bad last = 294360316
      Method 9, bad first &| = 278517241
      Method 9, bad first <=0 = 274241279
      Method 9, bad first min max = 278710563
      No extra mult = 97374387
      max
      Method 11 = 288414309
      Method 9, bad first = 318759227
      Method 9, bad last = 323014797
      Method 9, bad first &| = 320161362
      Method 9, bad first <=0 = 319378859
      Method 9, bad first min max = 346463028
      No extra mult = 256442927
      • 0
        Действительно, получается заметно лучше на 32-битах. Спасибо за этот вариант.
      • 0
        У меня получается различие на тесте ax=-1,ay=-1,bx=0,by=0.
        • 0
          А что должно быть?
          • 0
            Должно выдавать 0, но выдает 1 на первом if'e.
            • 0
              Да, видимо, в первом if должно быть строго меньше… Но в свете варианта от Mrrl проверять уже нет смысла :-)
  • 0
    А можете пояснить вот это:

    Проверяем вырожденный случай, когда точка c лежит на прямой (a,b). Ответа ноль можно избежать только если отрезок лежит точно на луче и обе точки находятся по одну сторону.

    Как отрезок может лежать на луче и при этом его точки будут «по одну сторону» от луча?
    • 0
      Моя ошибка. Исправил в статье.
  • 0
    А в этом языке есть есть битовые операции? Если есть, можно попробовать такой код. Я погонял на C# на небольших примерах, вроде бы работает.

            public static int check8(Point a,Point b,Point middle) {
                int ax = middle.x-a.x;
                int ay = middle.y-a.y;
                int bx = middle.x-b.x;
                int by = middle.y-b.y;
    
                if((ay|by)==0) return (((ax^bx)|((-ax)^(-bx)))>>31)+1;
                if((ay^by)>=0) return 1;
                long m=(long)ax*by-(long)ay*bx;
                return m==0 ? 0 : (((int)(m>>32)^by)>>30)|1;
            }
    
    • 0
      Ваш код даёт неверный результат на точках (1,0), (1,1), (0,0). По крайней мере, он не совпадает с результатом автора статьи (у вас -1, у автора 1)
      • 0
        Это очень странно — пример соответствует второму слева отрезку на картинке (который в верхней полуплоскости), а около него написано -1: точка лежит левее нижней вершины отрезка. Если семантика авторского кода отличается от примеров на картинках и -1 должно быть для точки, котрая лежит левее верхней вершины — то надо поменять знаки у ay и by. А учитывать точки, лежащие слева, к примеру, от a, просто нельзя — тогда для точек слева от верхней вершины треугольника у нас будет нечетное число пересечений. А они треугольнику, очевидно, не принадлежат.
        • 0
          Все-таки я на картинке перепутал знаки по сравнению с программой :(. Сейчас попытаюсь исправить картинку.
          Но, если программа выдает результаты с точностью до замены y на -y, то итоговый алгоритм проверки будет правильным.
          • 0
            Вопрос — всё-таки, -1 дают точки, лежащие слева от верхней вершины, или слева от точки a?
            • 0
              -1 получается, когда отрезок пересекает положительную часть оси x + часть вырожденных случаев, когда одна из точек лежит на положительной части оси x (в этом случае -1 выдается если другая точка лежит строго ниже).
              • 0
                Вопрос в том, какая «одна из точек»? a, b, верхняя или нижняя? Если a или b, то код не годится для проверки принадлежности (если использовать его, вызывая Check(P[n],P[n+1],M) для всех n): он выдаст нечетное число пересечений для точек слева от любой вершины — хоть правой, хоть верхней, хоть нижней.
                • 0
                  Вырожденные случаи, на которых выдается -1, это: (ay==0 && ax>0 && by<0) || (by==0 && bx>0 && ay<0).
                  • 0
                    То есть «слева от верхней точки». Тогда всё в порядке.
      • 0
        Да, у меня верхняя вершина треугольника выдаст 1 для всех сторон. Сейчас исправлю.
        • 0
          В текущем коде получается не ноль при ax=-1, ay=1, bx=0, by=0.
    • 0
      Замена ay < 0 ^ by < 0 на (ay^by) < 0 дает небольшое ускорение в 32-битном режиме, и почему-то небольшое замедление в 64-битном. Идея интересная.
  • +1
    Более правильный код:
            public static int check8a(Point a,Point b,Point middle) {
                int ax = a.x-middle.x;
                int ay = a.y-middle.y;
                int bx = b.x-middle.x;
                int by = b.y-middle.y;
    
                if((ax|ay)==0 || (bx|by)==0) return 0;
                if((ay|by)==0) return ((ax^bx)>>31)+1;
                if((ay^by)>=0) return 1;
                long m=(long)ax*by-(long)ay*bx;
                return m==0 ? 0 : (((int)(m>>32)^by)>>30)|1;
            }
    

    На C#, 64 бита он обгоняет check7 в 2 раза. А если точки передавать как ref — получается еще вдвое быстрее.
    • 0
      На 32 битах выигрыш у check7 составляет 2.5% :)
    • 0
      На Java у меня получился серьезный выигрыш в 32-битном режиме (~25%) и такая же скорость на 64-битном по сравнению с check2 и check3. Отличный вариант по скорости!
      • +1
        Можно еще немного ускорить (3-6%), но не уверен, что на всех компиляторах это даст выигрыш:

                public static int check8(Point a,Point b,Point middle) {
                    int ax = a.x-middle.x;
                    int ay = a.y-middle.y;
                    int bx = b.x-middle.x;
                    int by = b.y-middle.y;
        
                    if((ax|ay)==0 || (bx|by)==0) return 0;
                    if((ay|by)==0) return ((ax^bx)>>31)+1;
                    if(ay>=0) {
                        if(by>=0) return 1;
                        long m=(long)ax*by-(long)ay*bx;
                        return (int)(m>>63)|(int)((ulong)(-m)>>63);
                    } else {
                        if(by<0) return 1;
                        long m=(long)ay*bx-(long)ax*by;
                        return (int)(m>>63)|(int)((ulong)(-m)>>63);
                    }
                }
         

        Пока не могу придумать хороший метод найти int>=0? sign(long): -sign(long) — он помог бы сэкономить один условный переход.
        • +3
          Есть! Нам ведь не нужно учитывать случай m=-2^63… Правда, выигрыш достигается только на 32 битах.
                  public static int check8(Point a,Point b,Point middle) {
                      int ax = a.x-middle.x;
                      int ay = a.y-middle.y;
                      int bx = b.x-middle.x;
                      int by = b.y-middle.y;
          
                      if((ax|ay)==0 || (bx|by)==0) return 0;
                      if((ay|by)==0) return ((ax^bx)>>31)+1;
                      if((ay^by)>=0) return 1;
                      long m=(long)ax*by-(long)ay*bx;
                      return (((int)(m>>32)^ay)>>31)-(((int)((-m)>>32)^ay)>>31);
                  }
          
          • 0
            У меня наблюдается небольшое замедление на тесте, где чередуются знаки по координате y. Видимо Java не может векторизовать его.
            • 0
              На intel i5, 64 bit у меня вообще все результаты сравнялись — и check7, и разные версии check8. Разница меньше 10% и меняется от запуска к запуску.
              Надо это еще на C испытать.
  • 0
    Недавно нам рассказали еще один очень интересный способ. Если приблизить многоуголник лемнискатой, а это можно сделать с любой точностью, то подставив в уравнение лемнискаты (многочлен в C) координаты точки, можно сразу сказать, лежит она внутри, или вне.
    • 0
      Наверное, можно воспользоваться и интегралом Шварца-Кристоффеля. Но тогда проще посчитать число обходов.
  • +1
    Можно еще поменять порядок проверок.

            public static int check9(Point a,Point b,Point middle) {
                int ax = a.x-middle.x;
                int ay = a.y-middle.y;
                int bx = b.x-middle.x;
                int by = b.y-middle.y;
    
                if((ay^by)>=0) {
                    if((ax|ay)==0 || (bx|by)==0) return 0;
                    if((ay|by)!=0) return 1;
                    return ((ax^bx)>>31)+1;
                }
                long m=(long)ax*by-(long)ay*bx;
                return (int)((m^ay)>>63)-(int)(((-m)^ay)>>63);
            }
    

    Это 64-битная версия. Она даёт выигрыш, если стороны многоугольника часто пересекают прямую y=middle.y. Для 32-битной последний return надо заменить на

    return (((int)(m>>32)^ay)>>31)-(((int)((-m)>>32)^ay)>>31);
    

    Хорошая задачка, спасибо :)
  • НЛО прилетело и опубликовало эту надпись здесь
    • +1
      За O(log(N)) операций (на каждую проверку) и O(N^2) памяти. Режем плоскость на горизонтальные полосы линиями, проходящими через вершины, а каждая полоса режется на трапеции и треугольники сторонами многоугольника. Дальше — два бинарных поиска.
    • +1
      А если нужно проверить принадлежность K точек невыпуклому N-угольнику, причем все точки известны сразу (т.е. оффлайн-обработка), то это решается за O((K+N)*log(K+N)) и линейную память.
  • НЛО прилетело и опубликовало эту надпись здесь
    • 0
      В данном случае лучше (ay^by)<0. А вообще, я тоже согласен, что некрасиво. Как и любая конверсия булевского типа к int.
      Хотя, после того, как в ассемблере появились команды копирования флагов в регистры, жизнь компилятора стала проще. И такой код (ay < 0 ^ by < 0) может даже оказаться эффективным.
      • НЛО прилетело и опубликовало эту надпись здесь
        • 0
          Нет, «исключающее или».
          • НЛО прилетело и опубликовало эту надпись здесь
            • 0
              В этом топике слово «float» впервые встретилось в вашем комментарии — до сих пор вся работа шла с int или long.
              • НЛО прилетело и опубликовало эту надпись здесь

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