Комментарии 57
Может стоит добавить иллюстраций в статью? Без них читать не интересно, проще сразу скачать исходники.
+13
Если бы вы привели и другие алгоритмы проверки принадлежности, например разбивка невыпуклого многоугольника на выпуклые многоугольники или треугольники с последующей проверкой принадлежности к ним, то было бы интересней и полезней.
+1
Ваш способ на порядок быстрее ( если разбивка провелась заранее, что в жизни и делают для топологически статических полигонов)
+1
Алгоритм с разбиением хоть и возможен в принципе, но смысла в нем мало (на мой взгляд). А вот про алгоритм с подсчетом числа оборотов стоило упомянуть, так как он дает немного другой резльтат (для самопересекающихся многоугольников), что иногда и требуется. Очень нравится, как эти два алгоритма реализованы одновременно в Qt: qt.gitorious.org/qt/qt/blobs/4.7/src/gui/painting/qpolygon.cpp#line855
0
==
Она определяет пересекает ли луч выпущенный из точки c вдоль оси y в направлении уменьшения x отрезок (a,b)
==
Не очень понятна фраза «в направлении уменьшения x». Если луч идет вдоль оси y (параллелен ей), то x при этом не изменяется совсем. Или я просто неправильно себе все представляю?
Она определяет пересекает ли луч выпущенный из точки c вдоль оси y в направлении уменьшения x отрезок (a,b)
==
Не очень понятна фраза «в направлении уменьшения x». Если луч идет вдоль оси y (параллелен ей), то x при этом не изменяется совсем. Или я просто неправильно себе все представляю?
+1
Ваши функции поздно отсекают очевидные случаи. Во многих случаях умножать не надо вообще, а у вас всегда минимум два умножения. К примеру, рассмотрим случай, когда 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-битной машине не просто быстрее, а существенно быстрее :-)
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
Реабилитируюсь. Вот вам метод, который значимо быстрее всех на 32bit:
Экономия за счёт того, что во многих случаях int не превращаются в long вообще.
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
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). Ответа ноль можно избежать только если отрезок лежит точно на луче и обе точки находятся по одну сторону.
Как отрезок может лежать на луче и при этом его точки будут «по одну сторону» от луча?
Проверяем вырожденный случай, когда точка c лежит на прямой (a,b). Ответа ноль можно избежать только если отрезок лежит точно на луче и обе точки находятся по одну сторону.
Как отрезок может лежать на луче и при этом его точки будут «по одну сторону» от луча?
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, то итоговый алгоритм проверки будет правильным.
Но, если программа выдает результаты с точностью до замены 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 для всех сторон. Сейчас исправлю.
0
Замена ay < 0 ^ by < 0 на (ay^by) < 0 дает небольшое ускорение в 32-битном режиме, и почему-то небольшое замедление в 64-битном. Идея интересная.
0
Более правильный код:
На C#, 64 бита он обгоняет check7 в 2 раза. А если точки передавать как ref — получается еще вдвое быстрее.
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 — получается еще вдвое быстрее.
+1
На 32 битах выигрыш у check7 составляет 2.5% :)
0
На Java у меня получился серьезный выигрыш в 32-битном режиме (~25%) и такая же скорость на 64-битном по сравнению с check2 и check3. Отличный вариант по скорости!
0
Можно еще немного ускорить (3-6%), но не уверен, что на всех компиляторах это даст выигрыш:
Пока не могу придумать хороший метод найти int>=0? sign(long): -sign(long) — он помог бы сэкономить один условный переход.
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) — он помог бы сэкономить один условный переход.
+1
Есть! Нам ведь не нужно учитывать случай 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);
}
+3
У меня наблюдается небольшое замедление на тесте, где чередуются знаки по координате y. Видимо Java не может векторизовать его.
0
Недавно нам рассказали еще один очень интересный способ. Если приблизить многоуголник лемнискатой, а это можно сделать с любой точностью, то подставив в уравнение лемнискаты (многочлен в C) координаты точки, можно сразу сказать, лежит она внутри, или вне.
0
Можно еще поменять порядок проверок.
Это 64-битная версия. Она даёт выигрыш, если стороны многоугольника часто пересекают прямую y=middle.y. Для 32-битной последний return надо заменить на
Хорошая задачка, спасибо :)
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)) и линейную память.
+1
НЛО прилетело и опубликовало эту надпись здесь
Зарегистрируйтесь на Хабре, чтобы оставить комментарий
Проверка принадлежности точки невыпуклому многоугольнику