Комментарии 57
Может стоит добавить иллюстраций в статью? Без них читать не интересно, проще сразу скачать исходники.
Да-да! Картинки нужны для привлечения внимания, факт!
Добавил картинку. Надеюсь я не перепутал никаких случаев…
Теперь стало ещё непонятнее.
Где собственно пересечения луча «вдоль y в сторону уменьшения x» с отрезками?
Если бы вы привели и другие алгоритмы проверки принадлежности, например разбивка невыпуклого многоугольника на выпуклые многоугольники или треугольники с последующей проверкой принадлежности к ним, то было бы интересней и полезней.
Ваш способ на порядок быстрее ( если разбивка провелась заранее, что в жизни и делают для топологически статических полигонов)
Алгоритм с разбиением хоть и возможен в принципе, но смысла в нем мало (на мой взгляд). А вот про алгоритм с подсчетом числа оборотов стоило упомянуть, так как он дает немного другой резльтат (для самопересекающихся многоугольников), что иногда и требуется. Очень нравится, как эти два алгоритма реализованы одновременно в Qt: qt.gitorious.org/qt/qt/blobs/4.7/src/gui/painting/qpolygon.cpp#line855
Ну почему же мало. Мы, например, в программе все многоугольники храним в виде набора выпуклых. Это ускоряет работу с геометрией в разы, а разбивка производится всего один раз.
==
Она определяет пересекает ли луч выпущенный из точки c вдоль оси y в направлении уменьшения x отрезок (a,b)
==
Не очень понятна фраза «в направлении уменьшения x». Если луч идет вдоль оси y (параллелен ей), то x при этом не изменяется совсем. Или я просто неправильно себе все представляю?
Моя ошибка. Исправил + привел в соответствие с программами.
Ваши функции поздно отсекают очевидные случаи. Во многих случаях умножать не надо вообще, а у вас всегда минимум два умножения. К примеру, рассмотрим случай, когда ay и by оба положительны или оба отрицательны (на больших многоугольниках и случайных расположениях точки middle это может выполняться для большинства рёбер). Тогда ни один if не сработает, и вы вернёте 1, но при этом выполните два абсолютно лишних 64-битных умножения.
Вот такой вариант потестируйте, если я нигде не ошибся:
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;
}
На 64-битной JVM работает медленнее всех на полностью случайном тесте. Сейчас попробую еще на 32-битной проверить.
У меня на 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-битной машине не просто быстрее, а существенно быстрее :-)
На 32-битной JVM ваш метод быстрее check3 и чуть медленнее (на уровне погрешности измерений) check2.
Ну правильно, в check2 как раз эти случаи сразу и отсекаются первым if. Я больше check3 критиковал, а на check2 не глянул сразу.
Реабилитируюсь. Вот вам метод, который значимо быстрее всех на 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 вообще.
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
Действительно, получается заметно лучше на 32-битах. Спасибо за этот вариант.
Да, видимо, в первом if должно быть строго меньше… Но в свете варианта от Mrrl проверять уже нет смысла :-)
А можете пояснить вот это:

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

Как отрезок может лежать на луче и при этом его точки будут «по одну сторону» от луча?
А в этом языке есть есть битовые операции? Если есть, можно попробовать такой код. Я погонял на 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;
        }
Ваш код даёт неверный результат на точках (1,0), (1,1), (0,0). По крайней мере, он не совпадает с результатом автора статьи (у вас -1, у автора 1)
Это очень странно — пример соответствует второму слева отрезку на картинке (который в верхней полуплоскости), а около него написано -1: точка лежит левее нижней вершины отрезка. Если семантика авторского кода отличается от примеров на картинках и -1 должно быть для точки, котрая лежит левее верхней вершины — то надо поменять знаки у ay и by. А учитывать точки, лежащие слева, к примеру, от a, просто нельзя — тогда для точек слева от верхней вершины треугольника у нас будет нечетное число пересечений. А они треугольнику, очевидно, не принадлежат.
Все-таки я на картинке перепутал знаки по сравнению с программой :(. Сейчас попытаюсь исправить картинку.
Но, если программа выдает результаты с точностью до замены y на -y, то итоговый алгоритм проверки будет правильным.
Вопрос — всё-таки, -1 дают точки, лежащие слева от верхней вершины, или слева от точки a?
-1 получается, когда отрезок пересекает положительную часть оси x + часть вырожденных случаев, когда одна из точек лежит на положительной части оси x (в этом случае -1 выдается если другая точка лежит строго ниже).
Вопрос в том, какая «одна из точек»? a, b, верхняя или нижняя? Если a или b, то код не годится для проверки принадлежности (если использовать его, вызывая Check(P[n],P[n+1],M) для всех n): он выдаст нечетное число пересечений для точек слева от любой вершины — хоть правой, хоть верхней, хоть нижней.
Вырожденные случаи, на которых выдается -1, это: (ay==0 && ax>0 && by<0) || (by==0 && bx>0 && ay<0).
То есть «слева от верхней точки». Тогда всё в порядке.
Да, у меня верхняя вершина треугольника выдаст 1 для всех сторон. Сейчас исправлю.
Замена ay < 0 ^ by < 0 на (ay^by) < 0 дает небольшое ускорение в 32-битном режиме, и почему-то небольшое замедление в 64-битном. Идея интересная.
Более правильный код:
        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 — получается еще вдвое быстрее.
На Java у меня получился серьезный выигрыш в 32-битном режиме (~25%) и такая же скорость на 64-битном по сравнению с check2 и check3. Отличный вариант по скорости!
Можно еще немного ускорить (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) — он помог бы сэкономить один условный переход.
Есть! Нам ведь не нужно учитывать случай 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);
        }
У меня наблюдается небольшое замедление на тесте, где чередуются знаки по координате y. Видимо Java не может векторизовать его.
На intel i5, 64 bit у меня вообще все результаты сравнялись — и check7, и разные версии check8. Разница меньше 10% и меняется от запуска к запуску.
Надо это еще на C испытать.
Недавно нам рассказали еще один очень интересный способ. Если приблизить многоуголник лемнискатой, а это можно сделать с любой точностью, то подставив в уравнение лемнискаты (многочлен в C) координаты точки, можно сразу сказать, лежит она внутри, или вне.
Наверное, можно воспользоваться и интегралом Шварца-Кристоффеля. Но тогда проще посчитать число обходов.
Можно еще поменять порядок проверок.

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

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