Как я могу определить, находится ли 2D точка внутри многоугольника?

Я пытаюсь создать быструю 2D точку внутри алгоритма многоугольника для использования при тестировании попаданий (например Polygon.contains(p:Point)). Предложения для эффективных методов будут оценены.

Для графики я бы предпочел не предпочитать целые числа. Многие системы используют целые числа для рисования пользовательского интерфейса (в конце концов, пиксели - это целые числа), но, например, macOS использует float для всего. macOS знает только точки, и точка может переводиться в один пиксель, но в зависимости от разрешения монитора она может переводиться во что-то другое. На экранах сетчатки половина точки (0,5 / 0,5) составляет пиксель. Тем не менее, я никогда не замечал, что пользовательские интерфейсы macOS значительно медленнее, чем другие пользовательские интерфейсы. Ведь 3D API (OpenGL или Direct3D) также работают с плавающей точкой, а современные графические библиотеки очень часто используют преимущества ускорения GPU.

Теперь вы сказали, что скорость - ваша главная задача, хорошо, давайте перейдем к скорости. Прежде чем запускать какой-либо сложный алгоритм, сначала выполните простой тест. Создайте ограничивающий прямоугольник с осями вокруг вашего многоугольника. Это очень легко, быстро и уже может обезопасить вас от многих расчетов. Как это работает? Выполните итерацию по всем точкам многоугольника и найдите минимальные / максимальные значения X и Y.

Например, у вас есть очки (9/1), (4/3), (2/7), (8/2), (3/6). Это означает, что Xmin равно 2, Xmax равно 9, Ymin равно 1 и Ymax равно 7. Точка за пределами прямоугольника с двумя ребрами (2/1) и (9/7) не может быть внутри многоугольника.

// p is your point, p.x is the x coord, p.y is the y coord
if (p.x < Xmin || p.x > Xmax || p.y < Ymin || p.y > Ymax) {
    // Definitely not within the polygon!
}

Это первый тест для любой точки. Как видите, этот тест очень быстрый, но он также очень грубый. Для обработки точек, которые находятся внутри ограничительного прямоугольника, нам нужен более сложный алгоритм. Есть несколько способов, как это можно рассчитать. Какой метод работает, также зависит от того, может ли многоугольник иметь отверстия или всегда будет сплошным. Вот примеры твердых (одна выпуклая, одна вогнутая):

А вот один с дыркой:

Зеленый имеет отверстие в середине!

Самый простой алгоритм, который может обрабатывать все три приведенных выше случая и который все еще довольно быстр, называется приведением лучей . Идея алгоритма довольно проста: нарисуйте виртуальный луч из любой точки за пределами многоугольника к своей точке и посчитайте, как часто он попадает на сторону многоугольника. Если число совпадений четное, оно находится за пределами многоугольника, если оно нечетное, оно внутри.

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

У вас все еще есть ограничительная рамка выше, помните? Просто выберите точку за пределами ограничительной рамки и используйте ее в качестве отправной точки для вашего луча. Например, точка (Xmin - e/p.y)находится вне многоугольника точно.

А что есть e? Ну, e(на самом деле эпсилон) дает ограничивающей рамке немного отступов . Как я уже сказал, трассировка лучей завершится неудачей, если мы начнем слишком близко к линии многоугольника. Поскольку ограничивающий прямоугольник может равняться многоугольнику (если многоугольник является выровненным по оси прямоугольником, ограничивающий прямоугольник равен самому многоугольнику!), Нам нужно добавить отступы, чтобы сделать это безопасно, вот и все. Насколько большой вы должны выбрать e? Не слишком большой Это зависит от масштаба системы координат, который вы используете для рисования. Если ширина шага вашего пикселя равна 1,0, просто выберите 1,0 (но 0,1 также сработало бы)

Теперь, когда у нас есть луч с его начальной и конечной координатами, проблема смещается с « точки внутри многоугольника » на « как часто луч пересекает сторону многоугольника ». Поэтому мы не можем просто работать с точками многоугольника, как раньше, теперь нам нужны фактические стороны. Сторона всегда определяется двумя точками.

side 1: (X1/Y1)-(X2/Y2)
side 2: (X2/Y2)-(X3/Y3)
side 3: (X3/Y3)-(X4/Y4)
:

Вам нужно проверить луч со всех сторон. Считайте, что луч - вектор, а каждая сторона - вектор. Луч должен поразить каждую сторону ровно один раз или никогда. Он не может попасть в одну и ту же сторону дважды. Две линии в 2D-пространстве всегда будут пересекаться ровно один раз, если они не параллельны, и в этом случае они никогда не пересекаются. Однако, поскольку векторы имеют ограниченную длину, два вектора могут не быть параллельными и никогда не пересекаться, потому что они слишком короткие, чтобы когда-либо встречаться.

// Test the ray against all sides
int intersections = 0;
for (side = 0; side < numberOfSides; side++) {
    // Test if current side intersects with ray.
    // If yes, intersections++;
}
if ((intersections & 1) == 1) {
    // Inside of polygon
} else {
    // Outside of polygon
}

Пока все хорошо, но как проверить, пересекаются ли два вектора? Вот некоторый код на C (не проверенный), который должен сделать свое дело:

#define NO 0
#define YES 1
#define COLLINEAR 2

int areIntersecting(
    float v1x1, float v1y1, float v1x2, float v1y2,
    float v2x1, float v2y1, float v2x2, float v2y2
) {
    float d1, d2;
    float a1, a2, b1, b2, c1, c2;

    // Convert vector 1 to a line (line 1) of infinite length.
    // We want the line in linear equation standard form: A*x + B*y + C = 0
    // See: http://en.wikipedia.org/wiki/Linear_equation
    a1 = v1y2 - v1y1;
    b1 = v1x1 - v1x2;
    c1 = (v1x2 * v1y1) - (v1x1 * v1y2);

    // Every point (x,y), that solves the equation above, is on the line,
    // every point that does not solve it, is not. The equation will have a
    // positive result if it is on one side of the line and a negative one 
    // if is on the other side of it. We insert (x1,y1) and (x2,y2) of vector
    // 2 into the equation above.
    d1 = (a1 * v2x1) + (b1 * v2y1) + c1;
    d2 = (a1 * v2x2) + (b1 * v2y2) + c1;

    // If d1 and d2 both have the same sign, they are both on the same side
    // of our line 1 and in that case no intersection is possible. Careful, 
    // 0 is a special case, that's why we don't test ">=" and "<=", 
    // but "<" and ">".
    if (d1 > 0 && d2 > 0) return NO;
    if (d1 < 0 && d2 < 0) return NO;

    // The fact that vector 2 intersected the infinite line 1 above doesn't 
    // mean it also intersects the vector 1. Vector 1 is only a subset of that
    // infinite line 1, so it may have intersected that line before the vector
    // started or after it ended. To know for sure, we have to repeat the
    // the same test the other way round. We start by calculating the 
    // infinite line 2 in linear equation standard form.
    a2 = v2y2 - v2y1;
    b2 = v2x1 - v2x2;
    c2 = (v2x2 * v2y1) - (v2x1 * v2y2);

    // Calculate d1 and d2 again, this time using points of vector 1.
    d1 = (a2 * v1x1) + (b2 * v1y1) + c2;
    d2 = (a2 * v1x2) + (b2 * v1y2) + c2;

    // Again, if both have the same sign (and neither one is 0),
    // no intersection is possible.
    if (d1 > 0 && d2 > 0) return NO;
    if (d1 < 0 && d2 < 0) return NO;

    // If we get here, only two possibilities are left. Either the two
    // vectors intersect in exactly one point or they are collinear, which
    // means they intersect in any number of points from zero to infinite.
    if ((a1 * b2) - (a2 * b1) == 0.0f) return COLLINEAR;

    // If they are not collinear, they must intersect in exactly one point.
    return YES;
}

Входными значениями являются две конечные точки вектора 1 ( v1x1/v1y1и v1x2/v1y2) и вектора 2 ( v2x1/v2y1и v2x2/v2y2). Итак, у вас есть 2 вектора, 4 точки, 8 координат. YESи NOпонятно. YESувеличивает пересечения, NOничего не делает.

Как насчет Коллинеар? Это означает, что оба вектора лежат на одной и той же бесконечной линии, в зависимости от положения и длины, они вообще не пересекаются или пересекаются в бесконечном количестве точек. Я не совсем уверен, как справиться с этим делом, я бы не посчитал это пересечением в любом случае. Ну, во всяком случае, этот случай довольно редок на практике из-за ошибок округления с плавающей запятой; лучший код, вероятно, не будет проверяться, == 0.0fно вместо этого для чего-то вроде < epsilon, где epsilon довольно небольшое число.

Если вам нужно протестировать большее количество точек, вы, безусловно, можете немного ускорить все это, сохраняя в памяти стандартные формы сторон многоугольника в памяти, поэтому вам не придется каждый раз пересчитывать их. Это сэкономит вам два умножения с плавающей запятой и три вычитания с плавающей запятой при каждом тесте в обмен на сохранение в памяти трех значений с плавающей запятой на каждую сторону многоугольника. Это типичное соотношение памяти и времени вычислений.

И последнее, но не менее важное: если вы можете использовать 3D-оборудование для решения проблемы, есть интересная альтернатива. Просто пусть GPU сделает всю работу за вас. Создайте поверхность для рисования вне экрана. Заполните его полностью черным цветом. Теперь позвольте OpenGL или Direct3D нарисовать ваш многоугольник (или даже все ваши многоугольники, если вы просто хотите проверить, находится ли точка в каком-либо из них, но вам не важно, какая из них), и заполните многоугольники другим цвет, например белый. Чтобы проверить, находится ли точка внутри многоугольника, получите цвет этой точки с поверхности рисования. Это всего лишь O (1) выборка из памяти.

Конечно, этот метод можно использовать только в том случае, если ваша поверхность для рисования не обязательно должна быть огромной. Если он не может поместиться в память графического процессора, этот метод медленнее, чем в ЦП. Если он должен быть огромным, и ваш графический процессор поддерживает современные шейдеры, вы все равно можете использовать графический процессор, реализовав приведенное выше приведение лучей в качестве графического шейдера, что абсолютно возможно. Для большего количества полигонов или большого количества точек для тестирования это окупится, учитывая, что некоторые графические процессоры смогут тестировать от 64 до 256 точек параллельно. Тем не менее, обратите внимание, что передача данных из CPU в GPU и обратно всегда обходится дорого, поэтому для простого тестирования пары точек на пару простых многоугольников, где либо точки, либо многоугольники являются динамическими и будут часто меняться, подход с использованием графического процессора редко платит выкл.

Я думаю, что следующий фрагмент кода - лучшее решение (взято отсюда ):

int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy)
{
  int i, j, c = 0;
  for (i = 0, j = nvert-1; i < nvert; j = i++) {
    if ( ((verty[i]>testy) != (verty[j]>testy)) &&
     (testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
       c = !c;
  }
  return c;
}

аргументы

• nvert : количество вершин в многоугольнике. Повторять ли первую вершину в конце было обсуждено в статье, упомянутой выше.
• vertx, verty : Массивы, содержащие x- и y-координаты вершин многоугольника.
• testx, testy : координаты X и Y контрольной точки.

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

Идея этого довольно проста. Автор описывает это следующим образом:

Я запускаю полубесконечный луч по горизонтали (с увеличением x, фиксированный y) из контрольной точки и подсчитываю, сколько ребер он пересекает. На каждом пересечении луч переключается между внутренним и внешним. Это называется теоремой кривой Жордана.

Переменная c переключается с 0 на 1 и с 1 на 0 каждый раз, когда горизонтальный луч пересекает любое ребро. Так что в основном он отслеживает, является ли количество пересеченных ребер четным или нечетным. 0 означает четное, а 1 означает нечетное.

Вот версия ответа Nirg на языке C # от профессора RPI . Обратите внимание, что использование кода из этого источника RPI требует указания авторства.

Ограничительный флажок был добавлен вверху. Тем не менее, как указывает Джеймс Браун, основной код работает почти так же быстро, как и сам флажок ограничивающего прямоугольника, поэтому флажок ограничивающего прямоугольника может фактически замедлить всю операцию, в случае, если большинство проверяемых точек находятся внутри ограничивающего прямоугольника. , Таким образом, вы можете оставить флажок ограничивающим, или альтернативой может быть предварительный расчет ограничивающих прямоугольников ваших полигонов, если они не меняют форму слишком часто.

public bool IsPointInPolygon( Point p, Point[] polygon )
{
    double minX = polygon[ 0 ].X;
    double maxX = polygon[ 0 ].X;
    double minY = polygon[ 0 ].Y;
    double maxY = polygon[ 0 ].Y;
    for ( int i = 1 ; i < polygon.Length ; i++ )
    {
        Point q = polygon[ i ];
        minX = Math.Min( q.X, minX );
        maxX = Math.Max( q.X, maxX );
        minY = Math.Min( q.Y, minY );
        maxY = Math.Max( q.Y, maxY );
    }

    if ( p.X < minX || p.X > maxX || p.Y < minY || p.Y > maxY )
    {
        return false;
    }

    // https://wrf.ecse.rpi.edu/Research/Short_Notes/pnpoly.html
    bool inside = false;
    for ( int i = 0, j = polygon.Length - 1 ; i < polygon.Length ; j = i++ )
    {
        if ( ( polygon[ i ].Y > p.Y ) != ( polygon[ j ].Y > p.Y ) &&
             p.X < ( polygon[ j ].X - polygon[ i ].X ) * ( p.Y - polygon[ i ].Y ) / ( polygon[ j ].Y - polygon[ i ].Y ) + polygon[ i ].X )
        {
            inside = !inside;
        }
    }

    return inside;
}

Вот вариант ответа М. Каца на JavaScript, основанный на подходе Нирга:

function pointIsInPoly(p, polygon) {
    var isInside = false;
    var minX = polygon[0].x, maxX = polygon[0].x;
    var minY = polygon[0].y, maxY = polygon[0].y;
    for (var n = 1; n < polygon.length; n++) {
        var q = polygon[n];
        minX = Math.min(q.x, minX);
        maxX = Math.max(q.x, maxX);
        minY = Math.min(q.y, minY);
        maxY = Math.max(q.y, maxY);
    }

    if (p.x < minX || p.x > maxX || p.y < minY || p.y > maxY) {
        return false;
    }

    var i = 0, j = polygon.length - 1;
    for (i, j; i < polygon.length; j = i++) {
        if ( (polygon[i].y > p.y) != (polygon[j].y > p.y) &&
                p.x < (polygon[j].x - polygon[i].x) * (p.y - polygon[i].y) / (polygon[j].y - polygon[i].y) + polygon[i].x ) {
            isInside = !isInside;
        }
    }

    return isInside;
}

Вычислить ориентированную сумму углов между точкой р и каждым из вершин многоугольника. Если общий угол ориентирования составляет 360 градусов, точка находится внутри. Если сумма равна 0, точка находится за пределами.

Мне больше нравится этот метод, потому что он более надежный и менее зависит от точности чисел.

Методы, которые вычисляют равномерность количества пересечений, ограничены, потому что вы можете «поразить» вершину во время вычисления количества пересечений.

РЕДАКТИРОВАТЬ: Кстати, этот метод работает с вогнутыми и выпуклыми полигонами.

РЕДАКТИРОВАТЬ: я недавно нашел целую статью в Википедии по этой теме.

Этот вопрос так интересен. У меня есть другая работоспособная идея, отличная от других ответов на этот пост. Идея состоит в том, чтобы использовать сумму углов, чтобы решить, находится ли цель внутри или снаружи. Лучше известный как число обмоток .

Пусть х будет целевой точкой. Пусть массив [0, 1, .... n] будет всеми точками области. Соедините целевую точку с каждой граничной точкой линией. Если целевая точка находится внутри этой области. Сумма всех углов будет 360 градусов. Если не углы будут меньше 360.

Обратитесь к этому изображению, чтобы получить общее представление о идее:

Мой алгоритм предполагает, что по часовой стрелке это положительное направление. Вот потенциальный вклад:

[[-122.402015, 48.225216], [-117.032049, 48.999931], [-116.919132, 45.995175], [-124.079107, 46.267259], [-124.717175, 48.377557], [-122.92315, 47.047963], [-122.402015, 48.225216]]

Ниже приведен код Python, который реализует идею:

def isInside(self, border, target):
degree = 0
for i in range(len(border) - 1):
    a = border[i]
    b = border[i + 1]

    # calculate distance of vector
    A = getDistance(a[0], a[1], b[0], b[1]);
    B = getDistance(target[0], target[1], a[0], a[1])
    C = getDistance(target[0], target[1], b[0], b[1])

    # calculate direction of vector
    ta_x = a[0] - target[0]
    ta_y = a[1] - target[1]
    tb_x = b[0] - target[0]
    tb_y = b[1] - target[1]

    cross = tb_y * ta_x - tb_x * ta_y
    clockwise = cross < 0

    # calculate sum of angles
    if(clockwise):
        degree = degree + math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C)))
    else:
        degree = degree - math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C)))

if(abs(round(degree) - 360) <= 3):
    return True
return False

Статья Эрик Haines цитируется bobobobo действительно отлично. Особенно интересными являются таблицы, сравнивающие производительность алгоритмов; метод суммирования углов действительно плох по сравнению с другими. Также интересно то, что оптимизация, например использование поисковой сетки для дальнейшего разделения полигона на секторы «in» и «out», может сделать тест невероятно быстрым даже на полигонах с> 1000 сторонами.

Во всяком случае, это первые дни, но мой голос идет по методу «пересечений», который, по-моему, описывает Меки. Однако я обнаружил, что это наиболее кратко описано и кодифицировано Дэвидом Бурком . Мне нравится, что не требуется никакой реальной тригонометрии, она работает для выпуклых и вогнутых поверхностей и работает достаточно хорошо, так как количество сторон увеличивается.

Кстати, вот одна из таблиц производительности из статьи Эрика Хейнса для интереса, тестирование на случайных полигонах.

                       number of edges per polygon
                         3       4      10      100    1000
MacMartin               2.9     3.2     5.9     50.6    485
Crossings               3.1     3.4     6.8     60.0    624
Triangle Fan+edge sort  1.1     1.8     6.5     77.6    787
Triangle Fan            1.2     2.1     7.3     85.4    865
Barycentric             2.1     3.8    13.8    160.7   1665
Angle Summation        56.2    70.4   153.6   1403.8  14693

Grid (100x100)          1.5     1.5     1.6      2.1      9.8
Grid (20x20)            1.7     1.7     1.9      5.7     42.2
Bins (100)              1.8     1.9     2.7     15.1    117
Bins (20)               2.1     2.2     3.7     26.3    278

Свифт версия ответа по nirg :

extension CGPoint {
    func isInsidePolygon(vertices: [CGPoint]) -> Bool {
        guard !vertices.isEmpty else { return false }
        var j = vertices.last!, c = false
        for i in vertices {
            let a = (i.y > y) != (j.y > y)
            let b = (x < (j.x - i.x) * (y - i.y) / (j.y - i.y) + i.x)
            if a && b { c = !c }
            j = i
        }
        return c
    }
}

Очень нравится решение, опубликованное Nirg и отредактированное bobobobo. Я только сделал его дружественным к JavaScript и немного более разборчивым для моего использования:

function insidePoly(poly, pointx, pointy) {
    var i, j;
    var inside = false;
    for (i = 0, j = poly.length - 1; i < poly.length; j = i++) {
        if(((poly[i].y > pointy) != (poly[j].y > pointy)) && (pointx < (poly[j].x-poly[i].x) * (pointy-poly[i].y) / (poly[j].y-poly[i].y) + poly[i].x) ) inside = !inside;
    }
    return inside;
}

Я работал над этим, когда был исследователем под руководством Майкла Стоунбрейкера - вы знаете, профессора, который придумал Ingres , PostgreSQL и т. Д.

Мы поняли, что самым быстрым способом было сначала создать ограничивающую рамку, потому что это СУПЕР быстро. Если это вне ограничительной рамки, это снаружи. В противном случае, вы делаете тяжелую работу ...

Если вы хотите отличный алгоритм, посмотрите на исходный код PostgreSQL с открытым исходным кодом для гео работы ...

Я хочу отметить, что мы никогда не понимали право против левши (также выражаемое как проблема «внутри» и «снаружи» ...

ОБНОВИТЬ

Ссылка БКБ предоставила множество разумных алгоритмов. Я работал над проблемами наук о Земле и поэтому нуждался в решении, которое работает по широте / долготе, и у него есть особая проблема с ручностью - область внутри меньшей области или большей области? Ответ заключается в том, что «направление» вершин имеет значение - оно либо левостороннее, либо правостороннее, и таким образом вы можете указать любую область как «внутри» любого данного многоугольника. Таким образом, моя работа использовала решение три, перечисленное на этой странице.

Кроме того, в моей работе использовались отдельные функции для тестов «на линии».

... Так как кто-то спросил: мы выяснили, что тесты ограничивающего прямоугольника были наилучшими, когда число вершин превысило некоторое число - проведите очень быстрый тест, прежде чем выполнять более длинный тест, если это необходимо ... Ограничивающий прямоугольник создается простым взятием наибольший x, наименьший x, наибольший y и наименьший y и объединение их в четыре точки в квадрате ...

Еще один совет для тех, кто следует: мы выполнили все наши более сложные и «светорегуляторные» вычисления в пространстве сетки, все в положительных точках на плоскости, а затем снова проецировали обратно в «реальную» долготу / широту, избегая, таким образом, возможных ошибок обтекание, когда одна пересекает линию 180 долготы и при работе с полярными областями. Работал отлично!

Ответ Дэвида Сегонда в значительной степени является стандартным общим ответом, а Ричард Т - наиболее распространенная оптимизация, хотя есть и другие. Другие сильные оптимизации основаны на менее общих решениях. Например, если вы собираетесь проверять один и тот же многоугольник с множеством точек, триангуляция многоугольника может значительно ускорить процесс, поскольку существует ряд очень быстрых алгоритмов поиска TIN. Другой вариант: если многоугольник и точки находятся на ограниченной плоскости при низком разрешении, например, на экране, вы можете нарисовать многоугольник в отображаемом в память буфере отображения заданного цвета и проверить цвет данного пикселя, чтобы увидеть, лежит ли он в полигонах.

Как и многие оптимизации, они основаны на конкретных, а не общих случаях, и приносят выгоды, основанные на амортизированном времени, а не на единичном использовании.

Работая в этой области, я обнаружил, что вычислительная геометрия Джозефа О'Руркса в C 'ISBN 0-521-44034-3 очень помогает.

Тривиальным решением было бы разделить многоугольник на треугольники и нажать «проверить треугольники», как описано здесь

Если у вас многоугольник CONVEX, то может быть лучше. Посмотрите на многоугольник как на коллекцию бесконечных линий. Каждая строка делит пространство на две части. для каждой точки легко сказать, находится ли она на одной стороне или другой стороне линии. Если точка находится на одной стороне всех линий, она находится внутри многоугольника.

Я понимаю, что это старый, но вот алгоритм наложения лучей, реализованный в Какао, на случай, если кому-то будет интересно. Не уверен, что это самый эффективный способ сделать что-то, но он может помочь кому-то.

- (BOOL)shape:(NSBezierPath *)path containsPoint:(NSPoint)point
{
    NSBezierPath *currentPath = [path bezierPathByFlatteningPath];
    BOOL result;
    float aggregateX = 0; //I use these to calculate the centroid of the shape
    float aggregateY = 0;
    NSPoint firstPoint[1];
    [currentPath elementAtIndex:0 associatedPoints:firstPoint];
    float olderX = firstPoint[0].x;
    float olderY = firstPoint[0].y;
    NSPoint interPoint;
    int noOfIntersections = 0;

    for (int n = 0; n < [currentPath elementCount]; n++) {
        NSPoint points[1];
        [currentPath elementAtIndex:n associatedPoints:points];
        aggregateX += points[0].x;
        aggregateY += points[0].y;
    }

    for (int n = 0; n < [currentPath elementCount]; n++) {
        NSPoint points[1];

        [currentPath elementAtIndex:n associatedPoints:points];
        //line equations in Ax + By = C form
        float _A_FOO = (aggregateY/[currentPath elementCount]) - point.y;  
        float _B_FOO = point.x - (aggregateX/[currentPath elementCount]);
        float _C_FOO = (_A_FOO * point.x) + (_B_FOO * point.y);

        float _A_BAR = olderY - points[0].y;
        float _B_BAR = points[0].x - olderX;
        float _C_BAR = (_A_BAR * olderX) + (_B_BAR * olderY);

        float det = (_A_FOO * _B_BAR) - (_A_BAR * _B_FOO);
        if (det != 0) {
            //intersection points with the edges
            float xIntersectionPoint = ((_B_BAR * _C_FOO) - (_B_FOO * _C_BAR)) / det;
            float yIntersectionPoint = ((_A_FOO * _C_BAR) - (_A_BAR * _C_FOO)) / det;
            interPoint = NSMakePoint(xIntersectionPoint, yIntersectionPoint);
            if (olderX <= points[0].x) {
                //doesn't matter in which direction the ray goes, so I send it right-ward.
                if ((interPoint.x >= olderX && interPoint.x <= points[0].x) && (interPoint.x > point.x)) {  
                    noOfIntersections++;
                }
            } else {
                if ((interPoint.x >= points[0].x && interPoint.x <= olderX) && (interPoint.x > point.x)) {
                     noOfIntersections++;
                } 
            }
        }
        olderX = points[0].x;
        olderY = points[0].y;
    }
    if (noOfIntersections % 2 == 0) {
        result = FALSE;
    } else {
        result = TRUE;
    }
    return result;
}

Obj-C версия ответа nirg с примером метода для тестирования точек. Ответ Нирга хорошо сработал для меня.

- (BOOL)isPointInPolygon:(NSArray *)vertices point:(CGPoint)test {
    NSUInteger nvert = [vertices count];
    NSInteger i, j, c = 0;
    CGPoint verti, vertj;

    for (i = 0, j = nvert-1; i < nvert; j = i++) {
        verti = [(NSValue *)[vertices objectAtIndex:i] CGPointValue];
        vertj = [(NSValue *)[vertices objectAtIndex:j] CGPointValue];
        if (( (verti.y > test.y) != (vertj.y > test.y) ) &&
        ( test.x < ( vertj.x - verti.x ) * ( test.y - verti.y ) / ( vertj.y - verti.y ) + verti.x) )
            c = !c;
    }

    return (c ? YES : NO);
}

- (void)testPoint {

    NSArray *polygonVertices = [NSArray arrayWithObjects:
        [NSValue valueWithCGPoint:CGPointMake(13.5, 41.5)],
        [NSValue valueWithCGPoint:CGPointMake(42.5, 56.5)],
        [NSValue valueWithCGPoint:CGPointMake(39.5, 69.5)],
        [NSValue valueWithCGPoint:CGPointMake(42.5, 84.5)],
        [NSValue valueWithCGPoint:CGPointMake(13.5, 100.0)],
        [NSValue valueWithCGPoint:CGPointMake(6.0, 70.5)],
        nil
    ];

    CGPoint tappedPoint = CGPointMake(23.0, 70.0);

    if ([self isPointInPolygon:polygonVertices point:tappedPoint]) {
        NSLog(@"YES");
    } else {
        NSLog(@"NO");
    }
}

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

Основано на алгоритме симуляции простоты в http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html

Некоторые предикаты помощника:

exor(A,B):- \+A,B;A,\+B.
in_range(Coordinate,CA,CB) :- exor((CA>Coordinate),(CB>Coordinate)).

inside(false).
inside(_,[_|[]]).
inside(X:Y, [X1:Y1,X2:Y2|R]) :- in_range(Y,Y1,Y2), X > ( ((X2-X1)*(Y-Y1))/(Y2-Y1) +      X1),toggle_ray, inside(X:Y, [X2:Y2|R]); inside(X:Y, [X2:Y2|R]).

get_line(_,_,[]).
get_line([XA:YA,XB:YB],[X1:Y1,X2:Y2|R]):- [XA:YA,XB:YB]=[X1:Y1,X2:Y2]; get_line([XA:YA,XB:YB],[X2:Y2|R]).

Уравнение прямой с учетом 2 точек A и B (Линия (A, B)) имеет вид:

                    (YB-YA)
           Y - YA = ------- * (X - XA) 
                    (XB-YB) 

Важно, чтобы направление вращения линии было установлено по часовой стрелке для границ и против часовой стрелки для отверстий. Мы собираемся проверить, находится ли точка (X, Y), то есть проверенная точка, в левой полуплоскости нашей линии (это вопрос вкуса, это также может быть правая сторона, но также и направление границ линии должны быть изменены в этом случае), это проецирует луч из точки вправо (или влево) и подтверждает пересечение с линией. Мы решили проецировать луч в горизонтальном направлении (опять же, это вопрос вкуса, это также может быть сделано по вертикали с аналогичными ограничениями), поэтому мы имеем:

               (XB-XA)
           X < ------- * (Y - YA) + XA
               (YB-YA) 

Теперь нам нужно знать, находится ли точка только на левой (или правой) стороне сегмента линии, а не на всей плоскости, поэтому нам нужно ограничить поиск только этим сегментом, но это легко, поскольку находиться внутри сегмента только одна точка на линии может быть выше, чем Y на вертикальной оси. Поскольку это более сильное ограничение, оно должно быть первым, чтобы проверить, поэтому мы сначала выбираем только те строки, которые отвечают этому требованию, а затем проверяем его возможность. По теореме Кривая Джордана любой луч, проецируемый на многоугольник, должен пересекаться на четном числе линий. Итак, мы закончили, мы бросим луч вправо, а затем каждый раз, когда он пересекает линию, переключают свое состояние. Однако в нашей реализации мы собираемся проверить длину пакета решений, отвечающих заданным ограничениям, и принять решение об участии в нем. для каждой линии в многоугольнике это должно быть сделано.

is_left_half_plane(_,[],[],_).
is_left_half_plane(X:Y,[XA:YA,XB:YB], [[X1:Y1,X2:Y2]|R], Test) :- [XA:YA, XB:YB] =  [X1:Y1, X2:Y2], call(Test, X , (((XB - XA) * (Y - YA)) / (YB - YA) + XA)); 
                                                        is_left_half_plane(X:Y, [XA:YA, XB:YB], R, Test).

in_y_range_at_poly(Y,[XA:YA,XB:YB],Polygon) :- get_line([XA:YA,XB:YB],Polygon),  in_range(Y,YA,YB).
all_in_range(Coordinate,Polygon,Lines) :- aggregate(bag(Line),    in_y_range_at_poly(Coordinate,Line,Polygon), Lines).

traverses_ray(X:Y, Lines, Count) :- aggregate(bag(Line), is_left_half_plane(X:Y, Line, Lines, <), IntersectingLines), length(IntersectingLines, Count).

% This is the entry point predicate
inside_poly(X:Y,Polygon,Answer) :- all_in_range(Y,Polygon,Lines), traverses_ray(X:Y, Lines, Count), (1 is mod(Count,2)->Answer=inside;Answer=outside).

C # версия ответа nirg здесь: я просто поделюсь кодом. Это может сэкономить кому-то время.

public static bool IsPointInPolygon(IList<Point> polygon, Point testPoint) {
            bool result = false;
            int j = polygon.Count() - 1;
            for (int i = 0; i < polygon.Count(); i++) {
                if (polygon[i].Y < testPoint.Y && polygon[j].Y >= testPoint.Y || polygon[j].Y < testPoint.Y && polygon[i].Y >= testPoint.Y) {
                    if (polygon[i].X + (testPoint.Y - polygon[i].Y) / (polygon[j].Y - polygon[i].Y) * (polygon[j].X - polygon[i].X) < testPoint.X) {
                        result = !result;
                    }
                }
                j = i;
            }
            return result;
        }

Версия Java:

public class Geocode {
    private float latitude;
    private float longitude;

    public Geocode() {
    }

    public Geocode(float latitude, float longitude) {
        this.latitude = latitude;
        this.longitude = longitude;
    }

    public float getLatitude() {
        return latitude;
    }

    public void setLatitude(float latitude) {
        this.latitude = latitude;
    }

    public float getLongitude() {
        return longitude;
    }

    public void setLongitude(float longitude) {
        this.longitude = longitude;
    }
}

public class GeoPolygon {
    private ArrayList<Geocode> points;

    public GeoPolygon() {
        this.points = new ArrayList<Geocode>();
    }

    public GeoPolygon(ArrayList<Geocode> points) {
        this.points = points;
    }

    public GeoPolygon add(Geocode geo) {
        points.add(geo);
        return this;
    }

    public boolean inside(Geocode geo) {
        int i, j;
        boolean c = false;
        for (i = 0, j = points.size() - 1; i < points.size(); j = i++) {
            if (((points.get(i).getLongitude() > geo.getLongitude()) != (points.get(j).getLongitude() > geo.getLongitude())) &&
                    (geo.getLatitude() < (points.get(j).getLatitude() - points.get(i).getLatitude()) * (geo.getLongitude() - points.get(i).getLongitude()) / (points.get(j).getLongitude() - points.get(i).getLongitude()) + points.get(i).getLatitude()))
                c = !c;
        }
        return c;
    }

}

Порт .Net:

    static void Main(string[] args)
    {

        Console.Write("Hola");
        List<double> vertx = new List<double>();
        List<double> verty = new List<double>();

        int i, j, c = 0;

        vertx.Add(1);
        vertx.Add(2);
        vertx.Add(1);
        vertx.Add(4);
        vertx.Add(4);
        vertx.Add(1);

        verty.Add(1);
        verty.Add(2);
        verty.Add(4);
        verty.Add(4);
        verty.Add(1);
        verty.Add(1);

        int nvert = 6;  //Vértices del poligono

        double testx = 2;
        double testy = 5;


        for (i = 0, j = nvert - 1; i < nvert; j = i++)
        {
            if (((verty[i] > testy) != (verty[j] > testy)) &&
             (testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i]))
                c = 1;
        }
    }

VBA VERSION:

Примечание. Помните, что если полигон является областью на карте, то широта / долгота являются значениями Y / X, а не X / Y (широта = Y, долгота = X) из-за того, что, как я понимаю, являются историческими последствиями, когда Долгота не была измерением.

МОДУЛЬ КЛАССА: CPoint

Private pXValue As Double
Private pYValue As Double

'''''X Value Property'''''

Public Property Get X() As Double
    X = pXValue
End Property

Public Property Let X(Value As Double)
    pXValue = Value
End Property

'''''Y Value Property'''''

Public Property Get Y() As Double
    Y = pYValue
End Property

Public Property Let Y(Value As Double)
    pYValue = Value
End Property

МОДУЛЬ:

Public Function isPointInPolygon(p As CPoint, polygon() As CPoint) As Boolean

    Dim i As Integer
    Dim j As Integer
    Dim q As Object
    Dim minX As Double
    Dim maxX As Double
    Dim minY As Double
    Dim maxY As Double
    minX = polygon(0).X
    maxX = polygon(0).X
    minY = polygon(0).Y
    maxY = polygon(0).Y

    For i = 1 To UBound(polygon)
        Set q = polygon(i)
        minX = vbMin(q.X, minX)
        maxX = vbMax(q.X, maxX)
        minY = vbMin(q.Y, minY)
        maxY = vbMax(q.Y, maxY)
    Next i

    If p.X < minX Or p.X > maxX Or p.Y < minY Or p.Y > maxY Then
        isPointInPolygon = False
        Exit Function
    End If


    ' SOURCE: http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html

    isPointInPolygon = False
    i = 0
    j = UBound(polygon)

    Do While i < UBound(polygon) + 1
        If (polygon(i).Y > p.Y) Then
            If (polygon(j).Y < p.Y) Then
                If p.X < (polygon(j).X - polygon(i).X) * (p.Y - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then
                    isPointInPolygon = True
                    Exit Function
                End If
            End If
        ElseIf (polygon(i).Y < p.Y) Then
            If (polygon(j).Y > p.Y) Then
                If p.X < (polygon(j).X - polygon(i).X) * (p.Y - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then
                    isPointInPolygon = True
                    Exit Function
                End If
            End If
        End If
        j = i
        i = i + 1
    Loop   
End Function

Function vbMax(n1, n2) As Double
    vbMax = IIf(n1 > n2, n1, n2)
End Function

Function vbMin(n1, n2) As Double
    vbMin = IIf(n1 > n2, n2, n1)
End Function


Sub TestPointInPolygon()

    Dim i As Integer
    Dim InPolygon As Boolean

'   MARKER Object
    Dim p As CPoint
    Set p = New CPoint
    p.X = <ENTER X VALUE HERE>
    p.Y = <ENTER Y VALUE HERE>

'   POLYGON OBJECT
    Dim polygon() As CPoint
    ReDim polygon(<ENTER VALUE HERE>) 'Amount of vertices in polygon - 1
    For i = 0 To <ENTER VALUE HERE> 'Same value as above
       Set polygon(i) = New CPoint
       polygon(i).X = <ASSIGN X VALUE HERE> 'Source a list of values that can be looped through
       polgyon(i).Y = <ASSIGN Y VALUE HERE> 'Source a list of values that can be looped through
    Next i

    InPolygon = isPointInPolygon(p, polygon)
    MsgBox InPolygon

End Sub

Я сделал реализацию Python из nirg в C ++ код :

входные

• bounding_points: узлы, которые составляют многоугольник.

• bounding_box_positions: кандидат указывает на фильтр. (В моей реализации создано из ограничительной рамки.

  (Входы списки кортежей в формате: [(xcord, ycord), ...])

Возвращает

• Все точки, которые находятся внутри многоугольника.

def polygon_ray_casting(self, bounding_points, bounding_box_positions):
    # Arrays containing the x- and y-coordinates of the polygon's vertices.
    vertx = [point[0] for point in bounding_points]
    verty = [point[1] for point in bounding_points]
    # Number of vertices in the polygon
    nvert = len(bounding_points)
    # Points that are inside
    points_inside = []

    # For every candidate position within the bounding box
    for idx, pos in enumerate(bounding_box_positions):
        testx, testy = (pos[0], pos[1])
        c = 0
        for i in range(0, nvert):
            j = i - 1 if i != 0 else nvert - 1
            if( ((verty[i] > testy ) != (verty[j] > testy))   and
                    (testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i]) ):
                c += 1
        # If odd, that means that we are inside the polygon
        if c % 2 == 1: 
            points_inside.append(pos)


    return points_inside

Опять же, идея взята отсюда

Удивило, что никто не поднял этот вопрос раньше, кроме тех прагматиков, которым нужна база данных: MongoDB имеет отличную поддержку запросов Geo, включая этот.

То, что вы ищете, это:

db.neighborhoods.findOne ({geometry: {$ geoIntersects: {$ geometry: {тип: "Точка", координаты: ["долгота", "широта"]}}}})

Neighborhoodsэто коллекция, которая хранит один или несколько полигонов в стандартном формате GeoJson. Если запрос возвращает значение NULL, он не пересекается, иначе это так.

Очень хорошо задокументировано здесь: https://docs.mongodb.com/manual/tutorial/geospatial-tutorial/

Производительность для более чем 6000 точек, классифицированных в сетке из 330 нерегулярных многоугольников, была меньше одной минуты без какой-либо оптимизации вообще, включая время для обновления документов соответствующим полигоном.

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

/* math lib (defined below) */
static float dot_v2v2(const float a[2], const float b[2]);
static float angle_signed_v2v2(const float v1[2], const float v2[2]);
static void copy_v2_v2(float r[2], const float a[2]);

/* intersection function */
bool isect_point_poly_v2(const float pt[2], const float verts[][2], const unsigned int nr,
                         const bool use_holes)
{
    /* we do the angle rule, define that all added angles should be about zero or (2 * PI) */
    float angletot = 0.0;
    float fp1[2], fp2[2];
    unsigned int i;
    const float *p1, *p2;

    p1 = verts[nr - 1];

    /* first vector */
    fp1[0] = p1[0] - pt[0];
    fp1[1] = p1[1] - pt[1];

    for (i = 0; i < nr; i++) {
        p2 = verts[i];

        /* second vector */
        fp2[0] = p2[0] - pt[0];
        fp2[1] = p2[1] - pt[1];

        /* dot and angle and cross */
        angletot += angle_signed_v2v2(fp1, fp2);

        /* circulate */
        copy_v2_v2(fp1, fp2);
        p1 = p2;
    }

    angletot = fabsf(angletot);
    if (use_holes) {
        const float nested = floorf((angletot / (float)(M_PI * 2.0)) + 0.00001f);
        angletot -= nested * (float)(M_PI * 2.0);
        return (angletot > 4.0f) != ((int)nested % 2);
    }
    else {
        return (angletot > 4.0f);
    }
}

/* math lib */

static float dot_v2v2(const float a[2], const float b[2])
{
    return a[0] * b[0] + a[1] * b[1];
}

static float angle_signed_v2v2(const float v1[2], const float v2[2])
{
    const float perp_dot = (v1[1] * v2[0]) - (v1[0] * v2[1]);
    return atan2f(perp_dot, dot_v2v2(v1, v2));
}

static void copy_v2_v2(float r[2], const float a[2])
{
    r[0] = a[0];
    r[1] = a[1];
}

Примечание: это один из менее оптимальных методов, так как он включает в себя множество вызовов atan2f, но он может быть интересен разработчикам, читающим этот поток (в моих тестах он примерно в 23 раза медленнее, чем при использовании метода пересечения линий).

Для обнаружения попадания в полигон нам нужно проверить две вещи:

1. Если точка находится внутри полигона. (может быть выполнено с помощью алгоритма преобразования лучей)
2. Если точка находится на границе многоугольника (может быть выполнено тем же алгоритмом, который используется для обнаружения точки на полилинии (линии)).

Чтобы справиться со следующими частными случаями в алгоритме приведения лучей :

1. Луч перекрывает одну из сторон многоугольника.
2. Точка находится внутри многоугольника, и луч проходит через вершину многоугольника.
3. Точка находится за пределами многоугольника, и луч просто касается одного из углов многоугольника.

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

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

Или вы можете проверить, что сумма внутренних углов от вашей точки до каждой пары из двух последовательных вершин многоугольника до вашей контрольной точки равна 360, но у меня есть ощущение, что первый вариант быстрее, потому что он не включает деления или вычисления инверсии тригонометрических функций.

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

Вы также можете опубликовать вопрос в математическом сообществе. Могу поспорить, у них есть миллион способов сделать это

Если вы ищете библиотеку java-script, есть класс javascript google maps v3 для класса Polygon, чтобы определить, находится ли точка внутри него.

var polygon = new google.maps.Polygon([], "#000000", 1, 1, "#336699", 0.3);
var isWithinPolygon = polygon.containsLatLng(40, -90);

Google Extention Github

Когда используешь кварты(Qt 4.3+), можно использовать функцию QPolygon containsPoint

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

Для простых полигонов лучшим алгоритмом является алгоритм преобразования лучей (число пересечений). Для сложных полигонов этот алгоритм не обнаруживает точки, которые находятся внутри перекрывающихся областей. Так что для сложных полигонов вы должны использовать алгоритм числа обмоток.

Вот отличная статья с реализацией C обоих алгоритмов. Я попробовал их, и они хорошо работают.

http://geomalgorithms.com/a03-_inclusion.html

Scala-версия решения по nirg (предполагается, что предварительная проверка ограничивающего прямоугольника выполняется отдельно):

def inside(p: Point, polygon: Array[Point], bounds: Bounds): Boolean = {

  val length = polygon.length

  @tailrec
  def oddIntersections(i: Int, j: Int, tracker: Boolean): Boolean = {
    if (i == length)
      tracker
    else {
      val intersects = (polygon(i).y > p.y) != (polygon(j).y > p.y) && p.x < (polygon(j).x - polygon(i).x) * (p.y - polygon(i).y) / (polygon(j).y - polygon(i).y) + polygon(i).x
      oddIntersections(i + 1, i, if (intersects) !tracker else tracker)
    }
  }

  oddIntersections(0, length - 1, tracker = false)
}

Вот golang-версия ответа @nirg (вдохновлена ​​кодом C # @@ m-katz)

func isPointInPolygon(polygon []point, testp point) bool {
    minX := polygon[0].X
    maxX := polygon[0].X
    minY := polygon[0].Y
    maxY := polygon[0].Y

    for _, p := range polygon {
        minX = min(p.X, minX)
        maxX = max(p.X, maxX)
        minY = min(p.Y, minY)
        maxY = max(p.Y, maxY)
    }

    if testp.X < minX || testp.X > maxX || testp.Y < minY || testp.Y > maxY {
        return false
    }

    inside := false
    j := len(polygon) - 1
    for i := 0; i < len(polygon); i++ {
        if (polygon[i].Y > testp.Y) != (polygon[j].Y > testp.Y) && testp.X < (polygon[j].X-polygon[i].X)*(testp.Y-polygon[i].Y)/(polygon[j].Y-polygon[i].Y)+polygon[i].X {
            inside = !inside
        }
        j = i
    }

    return inside
}