halyavin Dec 6 2012 at 10:35

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

5 min

37K

Sport programming*Algorithms*

+47

Comments 57

s1im Dec 6 2012 at 11:32

Может стоит добавить иллюстраций в статью? Без них читать не интересно, проще сразу скачать исходники.

+13

sergof Dec 6 2012 at 12:25

Да-да! Картинки нужны для привлечения внимания, факт!

halyavin Dec 6 2012 at 12:53

Добавил картинку. Надеюсь я не перепутал никаких случаев…

qmax Dec 6 2012 at 13:19

Теперь стало ещё непонятнее.
Где собственно пересечения луча «вдоль y в сторону уменьшения x» с отрезками?

KvanTTT Dec 6 2012 at 12:27

Если бы вы привели и другие алгоритмы проверки принадлежности, например разбивка невыпуклого многоугольника на выпуклые многоугольники или треугольники с последующей проверкой принадлежности к ним, то было бы интересней и полезней.

PapaBubaDiop Dec 6 2012 at 12:43

Ваш способ на порядок быстрее ( если разбивка провелась заранее, что в жизни и делают для топологически статических полигонов)

kxyu Dec 6 2012 at 13:00

Алгоритм с разбиением хоть и возможен в принципе, но смысла в нем мало (на мой взгляд). А вот про алгоритм с подсчетом числа оборотов стоило упомянуть, так как он дает немного другой резльтат (для самопересекающихся многоугольников), что иногда и требуется. Очень нравится, как эти два алгоритма реализованы одновременно в Qt: qt.gitorious.org/qt/qt/blobs/4.7/src/gui/painting/qpolygon.cpp#line855

Woodroof Dec 8 2012 at 09:32

Ну почему же мало. Мы, например, в программе все многоугольники храним в виде набора выпуклых. Это ускоряет работу с геометрией в разы, а разбивка производится всего один раз.

tvolf Dec 6 2012 at 12:30

==
Она определяет пересекает ли луч выпущенный из точки c вдоль оси y в направлении уменьшения x отрезок (a,b)
==
Не очень понятна фраза «в направлении уменьшения x». Если луч идет вдоль оси y (параллелен ей), то x при этом не изменяется совсем. Или я просто неправильно себе все представляю?

halyavin Dec 6 2012 at 12:55

Моя ошибка. Исправил + привел в соответствие с программами.

tagir_valeev Dec 6 2012 at 12:56

Ваши функции поздно отсекают очевидные случаи. Во многих случаях умножать не надо вообще, а у вас всегда минимум два умножения. К примеру, рассмотрим случай, когда ay и by оба положительны или оба отрицательны (на больших многоугольниках и случайных расположениях точки middle это может выполняться для большинства рёбер). Тогда ни один if не сработает, и вы вернёте 1, но при этом выполните два абсолютно лишних 64-битных умножения.

tagir_valeev Dec 6 2012 at 13:01

Вот такой вариант потестируйте, если я нигде не ошибся:

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;
}

halyavin Dec 6 2012 at 13:14

На 64-битной JVM работает медленнее всех на полностью случайном тесте. Сейчас попробую еще на 32-битной проверить.

tagir_valeev Dec 6 2012 at 13:19

У меня на 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-битной машине не просто быстрее, а существенно быстрее :-)

halyavin Dec 6 2012 at 13:21

На 32-битной JVM ваш метод быстрее check3 и чуть медленнее (на уровне погрешности измерений) check2.

tagir_valeev Dec 6 2012 at 13:26

Ну правильно, в check2 как раз эти случаи сразу и отсекаются первым if. Я больше check3 критиковал, а на check2 не глянул сразу.

tagir_valeev Dec 6 2012 at 13:41

Реабилитируюсь. Вот вам метод, который значимо быстрее всех на 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 вообще.

tagir_valeev Dec 6 2012 at 13:45

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

halyavin Dec 6 2012 at 15:56

Действительно, получается заметно лучше на 32-битах. Спасибо за этот вариант.

halyavin Dec 6 2012 at 16:04

У меня получается различие на тесте ax=-1,ay=-1,bx=0,by=0.

Mrrl Dec 6 2012 at 16:08

А что должно быть?

halyavin Dec 6 2012 at 16:54

Должно выдавать 0, но выдает 1 на первом if'e.

tagir_valeev Dec 6 2012 at 18:16

Да, видимо, в первом if должно быть строго меньше… Но в свете варианта от Mrrl проверять уже нет смысла :-)

guyfawkes Dec 6 2012 at 14:29

А можете пояснить вот это:

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

Как отрезок может лежать на луче и при этом его точки будут «по одну сторону» от луча?

halyavin Dec 6 2012 at 15:57

Моя ошибка. Исправил в статье.

Mrrl Dec 6 2012 at 14:49

А в этом языке есть есть битовые операции? Если есть, можно попробовать такой код. Я погонял на 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;
        }

tagir_valeev Dec 6 2012 at 15:11

Ваш код даёт неверный результат на точках (1,0), (1,1), (0,0). По крайней мере, он не совпадает с результатом автора статьи (у вас -1, у автора 1)

Mrrl Dec 6 2012 at 15:34

Это очень странно — пример соответствует второму слева отрезку на картинке (который в верхней полуплоскости), а около него написано -1: точка лежит левее нижней вершины отрезка. Если семантика авторского кода отличается от примеров на картинках и -1 должно быть для точки, котрая лежит левее верхней вершины — то надо поменять знаки у ay и by. А учитывать точки, лежащие слева, к примеру, от a, просто нельзя — тогда для точек слева от верхней вершины треугольника у нас будет нечетное число пересечений. А они треугольнику, очевидно, не принадлежат.

halyavin Dec 6 2012 at 16:08

Все-таки я на картинке перепутал знаки по сравнению с программой :(. Сейчас попытаюсь исправить картинку.
Но, если программа выдает результаты с точностью до замены y на -y, то итоговый алгоритм проверки будет правильным.

Mrrl Dec 6 2012 at 16:11

Вопрос — всё-таки, -1 дают точки, лежащие слева от верхней вершины, или слева от точки a?

halyavin Dec 6 2012 at 16:22

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

Mrrl Dec 6 2012 at 16:30

Вопрос в том, какая «одна из точек»? a, b, верхняя или нижняя? Если a или b, то код не годится для проверки принадлежности (если использовать его, вызывая Check(P[n],P[n+1],M) для всех n): он выдаст нечетное число пересечений для точек слева от любой вершины — хоть правой, хоть верхней, хоть нижней.

halyavin Dec 6 2012 at 16:36

Вырожденные случаи, на которых выдается -1, это: (ay==0 && ax>0 && by<0) || (by==0 && bx>0 && ay<0).

Mrrl Dec 6 2012 at 16:41

То есть «слева от верхней точки». Тогда всё в порядке.

Mrrl Dec 6 2012 at 16:10

Да, у меня верхняя вершина треугольника выдаст 1 для всех сторон. Сейчас исправлю.

halyavin Dec 6 2012 at 16:16

В текущем коде получается не ноль при ax=-1, ay=1, bx=0, by=0.

halyavin Dec 6 2012 at 15:50

Замена ay < 0 ^ by < 0 на (ay^by) < 0 дает небольшое ускорение в 32-битном режиме, и почему-то небольшое замедление в 64-битном. Идея интересная.

Mrrl Dec 6 2012 at 16:32

Более правильный код:

        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 — получается еще вдвое быстрее.

Mrrl Dec 6 2012 at 16:45

На 32 битах выигрыш у check7 составляет 2.5% :)

halyavin Dec 6 2012 at 16:52

На Java у меня получился серьезный выигрыш в 32-битном режиме (~25%) и такая же скорость на 64-битном по сравнению с check2 и check3. Отличный вариант по скорости!

Mrrl Dec 6 2012 at 17:11

Можно еще немного ускорить (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) — он помог бы сэкономить один условный переход.

Mrrl Dec 6 2012 at 17:53

Есть! Нам ведь не нужно учитывать случай 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);
        }

halyavin Dec 6 2012 at 20:12

У меня наблюдается небольшое замедление на тесте, где чередуются знаки по координате y. Видимо Java не может векторизовать его.

Mrrl Dec 6 2012 at 20:22

На intel i5, 64 bit у меня вообще все результаты сравнялись — и check7, и разные версии check8. Разница меньше 10% и меняется от запуска к запуску.
Надо это еще на C испытать.

Laytlas Dec 6 2012 at 23:23

Недавно нам рассказали еще один очень интересный способ. Если приблизить многоуголник лемнискатой, а это можно сделать с любой точностью, то подставив в уравнение лемнискаты (многочлен в C) координаты точки, можно сразу сказать, лежит она внутри, или вне.

Mrrl Dec 6 2012 at 23:47

Наверное, можно воспользоваться и интегралом Шварца-Кристоффеля. Но тогда проще посчитать число обходов.

Mrrl Dec 7 2012 at 13:33

Можно еще поменять порядок проверок.

        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);

Хорошая задачка, спасибо :)