На мой взгляд, есть два подхода, в зависимости от того, насколько надежным вам является решение. Если вам не нужно решение с высокой точностью, вы можете легко написать функцию самостоятельно, используя трассировку лучей ( пересечение лучей и линий и точка в многоугольнике ).
Под «неточным решением» я подразумеваю то, что вы соглашаетесь с тем, что l ie указывает точно на линию или в узлах многоугольника и т.д. c. будет очень, очень и очень редко неправильно классифицироваться. На компьютерах x64 этой проблемой можно пренебречь, если вы не построите алгоритм триангуляции высокой точности или что-то подобное (чего вы, вероятно, не сделали бы в R).
В таком случае это действительно просто, но как упоминалось ранее, такая функция иногда может дать сбой из-за арифметики с плавающей запятой c точности fl aws. Дополнительную информацию по этому вопросу можно найти в книгах по вычислительной геометрии и в документах CGAL .
Функция basi c может выглядеть так:
section.ray.intersection <- function(p,d,a,b)
# this function tests if a ray r=[p,d] intersects the segment [a,b]
# p - point c(x,y)
# d - ray direction c(x,y)
# a - one of the nodes of the section [a,b] (not ordered with respect to computational geometry standards)
# b - remaining node of the section [a,b]
{
v1 = p-a
v2 = b-a
v3 = c(-d[2],d[1])
t1 = (v2[1]*v1[2] - v2[2]*v1[1])/(v2 %*% v3)
t2 = (v1 %*% v3)/(v2 %*% v3)
return(t1 >=0. & (t2 >= 0. & t2 <= 1.0))
}
in.polygon <- function(px,py,Ax,Ay)
# basic function testing if the point is in the polygon
# px - x coordinate of the point
# py - y coordinate of the point
# Ax - a vector of x coordinates of the polygon
# Ay - a vector of y coordinates of the polygon
{
n <- length(Ax) # gent the number of nodes of the polygon
if (n != length(Ay)) return(NA) # can be ommitted if Ax and Ay are always valid
d = runif(2) # get the random vector of the direction
is = 0 # number of intersections
for (i in 2:n) # go trough all edges
is <- is + section.ray.intersection(c(px,py),d,c(Ax[i-1],Ay[i-1]),c(Ax[i],Ay[i])) # and count TRUEs
return( !(is %% 2==0)) # if the number of intersections is odd the point is inside
}
И пример использования:
# libs
library(ggplot2)
library(purrr)
library(dplyr)
# polygon
Ax <- c(0,1,2,0.25,0)
Ay <- c(0,0,1,1.00,0)
# points
N = 1000 # number of points
p <- data.frame(x = rnorm(N,0.5,.4), y= rnorm(N,0.5,.4)) # data frame od point coordinates
# apply the function to points - this is a "for loop" equivalent
p$is.inside <- apply(p, 1, function(X){return(in.polygon(X[1],X[2],Ax,Ay))})
# make a plot
ggplot() + geom_path(data=data.frame(x=Ax,y=Ay),aes(x,y)) + geom_point(data=p,aes(x,y,col = is.inside))
Однако, если вам нужно немного более точное решение, вы можете провести несколько тестов, а затем использовать некоторую методологию голосования. Например, предположим, что во время одного теста риск неправильной классификации составляет примерно p = 0,001 (что на самом деле намного ниже), тогда вероятность того, что во всех трех тестах (со случайными векторами направления) он будет классифицирован неправильно, будет только ppp = 0,000000001. Измененная версия функции может выглядеть так:
in.polygon.robust <- function(px,py,Ax,Ay)
{
n <- length(Ax)
if (n != length(Ay)) return(NA) # can be ommitted if Ax and Ay are always valid
d1 <- runif(2) # direction 1
d2 <- runif(2) # direction 2
d3 <- runif(2) # direction 3
is1 <- 0 # number of intersections in direction 1
is2 <- 0 # number of intersections in direction 2
is3 <- 0 # number of intersections in direction 3
for (i in 2:n) # go trough all edges
{
P <- c(px,py)
Ap1 <- c(Ax[i-1],Ay[i-1])
Ap2 <- c(Ax[i],Ay[i])
# count intersections
is1 <- is1 + section.ray.intersection(P,d1,Ap1,Ap2)
is2 <- is2 + section.ray.intersection(P,d2,Ap1,Ap2)
is3 <- is3 + section.ray.intersection(P,d3,Ap1,Ap2)
}
r <- (is1 %% 2==0) + (is2 %% 2==0) + (is3 %% 2==0) # sum the even outcomes
return( !(r>=2)) # return voting outcome
}
Однако, если этого все еще недостаточно, у вас есть возможность использовать, например, CGAL через R Cpp или ( R-привязки ), или если лицензия вам не подходит, вам нужно написать свой собственный надежный алгоритм, что действительно сложно.
Другое дело - проблема эффективности. Если вы используете sh, чтобы применить этот тест к множеству многоугольников (особенно описанных с большим количеством ребер), вам нужно использовать некоторые специфические c структуры данных, чтобы сделать все приемлемо эффективным. Более подробную информацию об этом можно найти в этой книге или в этой книге . И напишите это все на C или C ++.