найти точку пересечения в R

У меня есть 2 вектора:

set.seed(1)
x1 = rnorm(100,0,1)
x2 = rnorm(100,1,1)

Я хочу построить их как линии, а затем найти точки пересечения линий, а также, если есть несколько точек пересечения, я хочу найти каждую из них.

введите описание изображения здесь

Я столкнулся с подобным вопросом и попытался решить эту проблему, используя spatstat, но мне не удалось преобразовать комбинированный фрейм данных, содержащий оба значения вектора, в psp object.


person SBS    schedule 11.12.2013    source источник
comment
Вы хотите найти все пересечения линий в plot(x1,x2, type='l') ?   -  person Stephen Henderson    schedule 11.12.2013
comment
Или вы имеете в виду пересечения plot(seq_along(x1), x1, type='l') и lines(seq_along(x2), x2, type='l', col="red")   -  person Stephen Henderson    schedule 11.12.2013
comment
Мне нужны координаты везде, где есть пересечение, я привел приведенные выше векторы в качестве игрушечных примеров, но мой фактический ряд нелинейный, уравнение которого не указано.   -  person SBS    schedule 11.12.2013
comment
Я имею в виду график (seq_along (x1), x1, тип = 'l') и строки (seq_along (x2), x2, тип = 'l', col = красный)   -  person SBS    schedule 11.12.2013
comment
Вы можете попробовать метод фиксированной точки Ньютона: en.wikipedia.org/wiki/Newton%27s_method   -  person lucas92    schedule 11.12.2013
comment
@ lucas92 lucas92, не могли бы вы предложить какой-нибудь пакет или функцию в R, которая могла бы это сделать?   -  person SBS    schedule 11.12.2013
comment
Я не программирую на r, но вам нужно найти пересечение двух функций. 1) Найдите две функции. 2) h(x) = f(x) - g(x), где f и g - две функции. 3) Применим метод Ньютона, используя функцию h(x). Ваш алгоритм должен попробовать несколько x0, если есть несколько точек пересечения. Я предлагаю вам начать тестирование метода с помощью двух линейных функций и посмотреть, сможете ли вы получить точку пересечения.   -  person lucas92    schedule 11.12.2013


Ответы (4)


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

set.seed(2)
x1 <- sample(1:10, 100, replace = TRUE)
x2 <- sample(1:10, 100, replace = TRUE)

# Find points where x1 is above x2.
above <- x1 > x2

# Points always intersect when above=TRUE, then FALSE or reverse
intersect.points <- which(diff(above) != 0)

# Find the slopes for each line segment.
x1.slopes <- x1[intersect.points+1] - x1[intersect.points]
x2.slopes <- x2[intersect.points+1] - x2[intersect.points]

# Find the intersection for each segment.
x.points <- intersect.points + ((x2[intersect.points] - x1[intersect.points]) / (x1.slopes-x2.slopes))
y.points <- x1[intersect.points] + (x1.slopes*(x.points-intersect.points))

# Joint points
joint.points <- which(x1 == x2)
x.points <- c(x.points, joint.points)
y.points <- c(y.points, x1[joint.points])

# Plot points
plot(x1,type='l')
lines(x2,type='l',col='red')
points(x.points,y.points,col='blue')

# Segment overlap
start.segment <- joint.points[-1][diff(joint.points) == 1] - 1
for (i in start.segment) lines(x = c(i, i+1), y = x1[c(i, i+1)], col = 'blue')

введите здесь описание изображения

person nograpes    schedule 11.12.2013
comment
@baptiste Я думаю, что вертикальные сегменты невозможны. Возможно, вы имели в виду горизонтальные сегменты, потому что да, горизонтальные перекрывающиеся сегменты могут создать проблему. На самом деле, также было бы проблемой, если бы точки пересекались точно в оцениваемой точке. Удобно, что вы можете проверить оба одновременно: проверить, равны ли точки в вашем векторе, добавить их к набору пересекающихся точек, а затем проверить, встречаются ли эти точки одна за другой. Может вы не это имели в виду? - person nograpes; 11.12.2013
comment
нет, я просто думал об общей проблеме пересечения сегментов, но здесь вы правы, ни один сегмент не является вертикальным (что вызвало бы проблемы с бесконечным наклоном). - person baptiste; 11.12.2013
comment
@nograpes, ваш ответ отлично сработал для класса timeSeries ... но в моем реальном приложении я имею дело с объектом xts, как с двумя временными рядами, и когда я пытаюсь построить точки с помощью xts, это требует x.points также быть объектом xts, тогда как x.points фактически возвращает числовые данные (номера индексов), которые я не могу преобразовать в объект xts, как мне с этим справиться? - person SBS; 13.12.2013
comment
@nograpes Я очень ценю ваше решение, оно простое, но есть один момент, который я до сих пор не понимаю; почему вы говорите above <- x1 > x2, а не above <- x1 >= x2, поскольку точки, которые уже имеют одинаковое значение, пересекаются? - person ExploreR; 15.01.2020
comment
Там есть ошибка; код не всегда будет работать, когда точки равны по определенному индексу. Но изменение на above <- x1 >= x2 не исправит ситуацию, вы должны явно обрабатывать эти случаи. - person nograpes; 16.01.2020
comment
Я исправил ошибку и добавил новую функцию, когда сегменты полностью перекрываются. В новом примере вы можете увидеть одну ситуацию по индексу 59, где сегменты перекрываются, и вы можете увидеть по первому индексу, что они совпадают по одному индексу. - person nograpes; 16.01.2020

Вот альтернативный код пересечения сегментов,

# segment-segment intersection code
# http://paulbourke.net/geometry/pointlineplane/
ssi <- function(x1, x2, x3, x4, y1, y2, y3, y4){

  denom <- ((y4 - y3)*(x2 - x1) - (x4 - x3)*(y2 - y1))
  denom[abs(denom) < 1e-10] <- NA # parallel lines

  ua <- ((x4 - x3)*(y1 - y3) - (y4 - y3)*(x1 - x3)) / denom
  ub <- ((x2 - x1)*(y1 - y3) - (y2 - y1)*(x1 - x3)) / denom

  x <- x1 + ua * (x2 - x1)
  y <- y1 + ua * (y2 - y1)
  inside <- (ua >= 0) & (ua <= 1) & (ub >= 0) & (ub <= 1)
  data.frame(x = ifelse(inside, x, NA), 
             y = ifelse(inside, y, NA))

}
# do it with two polylines (xy dataframes)
ssi_polyline <- function(l1, l2){
  n1 <- nrow(l1)
  n2 <- nrow(l2)
  stopifnot(n1==n2)
  x1 <- l1[-n1,1] ; y1 <- l1[-n1,2] 
  x2 <- l1[-1L,1] ; y2 <- l1[-1L,2] 
  x3 <- l2[-n2,1] ; y3 <- l2[-n2,2] 
  x4 <- l2[-1L,1] ; y4 <- l2[-1L,2] 
  ssi(x1, x2, x3, x4, y1, y2, y3, y4)
}
# do it with all columns of a matrix
ssi_matrix <- function(x, m){
  # pairwise combinations
  cn <- combn(ncol(m), 2)
  test_pair <- function(i){
    l1 <- cbind(x, m[,cn[1,i]])
    l2 <- cbind(x, m[,cn[2,i]])
    pts <- ssi_polyline(l1, l2)
    pts[complete.cases(pts),]
  }
  ints <- lapply(seq_len(ncol(cn)), test_pair)
  do.call(rbind, ints)

}
# testing the above
y1 = rnorm(100,0,1)
y2 = rnorm(100,1,1)
m = cbind(y1, y2)
x = 1:100
matplot(x, m, t="l", lty=1)
points(ssi_matrix(x, m))

введите описание изображения здесь

person baptiste    schedule 29.08.2015
comment
Предлагаемая настройка: использовать короткие логические операции, inside <- ua >0 $$ ua < 1 && ... для скорости. - person Carl Witthoft; 13.08.2020

Поздний ответ, но вот "пространственный" метод с использованием пакета SP и RGEOS. Для этого требуется, чтобы и x, и y были числовыми (или могли быть преобразованы в числовые). Проекция произвольная, но epsg:4269, похоже, работает хорошо:

library(sp)
library(rgeos)
# dummy x data
x1 = rnorm(100,0,1)
x2 = rnorm(100,1,1)

#dummy y data 
y1 <- seq(1, 100, 1)
y2 <- seq(1, 100, 1) 

# convert to a sp object (spatial lines)
l1 <- Line(matrix(c(x1, y1), nc = 2, byrow = F))
l2 <- Line(matrix(c(x2, y2), nc = 2, byrow = F))
ll1 <- Lines(list(l1), ID = "1")
ll2 <- Lines(list(l2), ID = "1")
sl1 <- SpatialLines(list(ll1), proj4string = CRS("+init=epsg:4269"))
sl2 <- SpatialLines(list(ll2), proj4string = CRS("+init=epsg:4269"))

# Calculate locations where spatial lines intersect
int.pts <- gIntersection(sl1, sl2, byid = TRUE)
int.coords <- int.pts@coords

# Plot line data and points of intersection
plot(x1, y1, type = "l")
lines(x2, y2, type = "l", col = "red")
points(int.coords[,1], int.coords[,2], pch = 20, col = "blue")
person ken    schedule 20.04.2017

Мне нужно было пересечение для другого приложения, и я обнаружил, что ответ nograpes работает неправильно:

# another example
x=seq(-4,6,length.out=10)
x1=dnorm(x, 0, 1)
x2=dnorm(x,2,2)

# Find points where x1 is above x2.
above <- x1 > x2

# Points always intersect when above=TRUE, then FALSE or reverse
intersect.points <- which(diff(above) != 0)

# Find the slopes for each line segment.
x1.slopes <- x1[intersect.points+1] - x1[intersect.points]
x2.slopes <- x2[intersect.points+1] - x2[intersect.points]

# Find the intersection for each segment.
x.points <- x[intersect.points] + ((x2[intersect.points] - x1[intersect.points]) / (x1.slopes-x2.slopes))
y.points <- x1[intersect.points] + (x1.slopes*(x.points-x[intersect.points]))

# Joint points
joint.points <- which(x1 == x2)
x.points <- c(x.points, joint.points)
y.points <- c(y.points, x1[joint.points])

# Plot points 
# length(x); length(x1)
plot(x, x1,type='l')
lines(x, x2,type='l',col='red')
points(x.points,y.points,col='blue')

Для бинормального распределения точки отключены

Для этих бинормальных распределений точки не совпадают, в данном случае особенно правая точка пересечения. Это происходит, когда значения на оси x не являются последовательными целыми числами и, следовательно, не имеют разницы в 1 для последовательных точек. Я заменил intersect.points на x[intersect.points], чего недостаточно. Жаль, так как метод относительно прост по сравнению с другими методами. Метод, предоставленный baptiste, работает намного лучше:

m = cbind(x1, x2)
matplot(x, m, t="l", lty=1)
points(ssi_matrix(x, m))

Следуя той же идее, более общая реализация, учитывающая различия последовательных значений x != 1:

intersect.2lines <- function (x, y1, y2){
  above = y1 > y2
  intersect.points <- which(diff(above) != 0) 
  
  y1.diff <- y1[intersect.points+1] - y1[intersect.points]
  y2.diff <- y2[intersect.points+1] - y2[intersect.points]
  x.diff <- x[intersect.points+1]-x[intersect.points]
  
  slope1 = y1.diff/x.diff
  slope2 = y2.diff/x.diff
  intercept1 = y1[intersect.points]-slope1*x[intersect.points]
  intercept2 = y2[intersect.points]-slope2*x[intersect.points]
  x.points = ifelse(slope1==slope2, NA, 
                  (intercept2-intercept1)/(slope1-slope2))
  y.points = ifelse(slope1==slope2, NA,
                  slope1*x.points+intercept1)
  # Joint points
  joint.points <- which(y1 == y2)
  x.points <- c(x.points, x[joint.points])
  y.points <- c(y.points, y1[joint.points])
  return(data.frame(x.points=x.points, y.points=y.points))
}

Это реализация формулы, приведенной в Википедии, с учетом двухстрочных уравнений Line пересечение линий

Результаты теперь идентичны результатам, полученным методом баптиста.

person tak101    schedule 28.03.2021