如何根据截距和斜率（包括原点）找到线之间的曲面？

Question

我正在寻找一种方法来可视化多条直线之间的表面，这些直线是通过截距和斜率在数据框中定义的。我正在寻找的曲面是包含原点 (0, 0) 的曲面。

\n

线的数量可能会有所不同（即使在下面的简化示例中我只有 6 条），并且其中一些可能是多余的（即它们没有包围我正在寻找的表面，因为其他线更具约束性）。

\n

让我们采用这个简单的数据框：

\n

df <- data.frame("Line" = c("A", "B", "C", "D", "E", "F"),\n                 "Intercept" = c(4, 3, -2.5, -1.5, -5, -.5),\n                 "Slope" = c(-1, 1, 2.4, -.6, -.8, .6))\n

\n

绘制这些线ggplot2：

\n

ggplot(data = df) +\n  geom_vline(xintercept = 0) +\n  geom_hline(yintercept = 0) +\n  geom_abline(mapping = aes(intercept = Intercept, slope = Slope),\n              colour = "red") +\n  coord_cartesian(xlim = c(-6, 6), ylim = c(-6, 6))\n

\n

给我以下输出：

\n

\n

基本上我想找到包围原点 (0, 0) 的线之间的交点，忽略多余的线（在本例中为左下角，截距 = -5 且斜率 = -0.8）。然后，这 5 个交点将用于绘制凸包。

\n

我的主要问题在于找到约束线的交点（下面的绿色点）以便能够找到蓝色表面。

\n

\n

问题：关于如何在 R 中处理这个问题，有什么建议吗？最好是以可以扩展到更大数据帧（包括更多约束和冗余行）的方式？

\n

附加问题： geom_abline()没有类似于 的群体geom_line()审美，可用于识别该行。ggplot2有谁知道根据斜率和截距（或直线的两个用户定义点）绘制直线的解决方法？

\n

预先感谢您的任何建议或（部分）潜在解决方案！

\n

更新：为了使多边形围绕点 (a,b) 而不是原点 (0, 0) 居中，我修改了原始代码（特别是 @AllanCameron 的 \xc3\xacnnermost()` 函数，如下所示：

\n

innermost <- function(slopes, intercepts, a, b) {\n\n  meetings <- function(slopes, intercepts) {\n    meets_at <- function(i1, s1, i2, s2) {\n      ifelse(s1 - s2 == 0, NA, (i2 - i1)/(s1 - s2))\n    }\n    xvals <- outer(seq_along(slopes), seq_along(slopes), function(i, j) {\n      meets_at(intercepts[i], slopes[i], intercepts[j], slopes[j])\n    })\n\n    yvals <- outer(seq_along(slopes), seq_along(slopes), function(i, j) {\n      intercepts + slopes *\n        meets_at(intercepts[i], slopes[i], intercepts[j], slopes[j])\n    })\n\n    cbind(x = xvals[lower.tri(xvals)], y = yvals[lower.tri(yvals)])\n\n  }\n\n  xy <- meetings(slopes, intercepts)\n  xy[,1] <- xy[,1] - a\n  xy[,2] <- xy[,2] - b\n\n  is_cut <- function(x, y, slopes, intercepts, a, b) {\n    d <- sqrt(x^2 + y^2)\n    slope <- y / x\n    xvals <- (intercepts + slopes*a - b) / (slope - slopes)\n    yvals <- xvals * slopes + (intercepts + slopes*a - b)\n    ds <- sqrt(xvals^2 + yvals^2)\n    any(d - ds > 1e-6 & sign(xvals) == sign(x) & sign(yvals) == sign(y))\n  }\n\n  xy <- xy[sapply(seq(nrow(xy)), function(i) {\n    !is_cut(xy[i, 1], xy[i, 2], slopes, intercepts, a, b)\n  }),]\n\n  xy <- xy[order(atan2(xy[,2], xy[,1])),]\n  \n  xy[,1] <- xy[,1] + a\n  xy[,2] <- xy[,2] + b\n  \n  as.data.frame(rbind(xy, xy[1,]))\n}\n

\n

Answer 1

这是一个只需要一点几何和代数的解决方案，仅使用基本 R。我们可以定义一个函数，innermost查找内部多边形的 x、y 坐标，并以逆时针顺序将它们作为数据框返回。这允许您通过执行以下操作来创建 ggplot：

ggplot(data = df) +
  geom_vline(xintercept = 0) +
  geom_hline(yintercept = 0) +
  geom_abline(mapping = aes(intercept = Intercept, slope = Slope),
              colour = "red") +
  geom_polygon(data = innermost(df$Slope, df$Intercept), aes(x, y),
               fill = "#99d9ea") +
  geom_point(data = innermost(df$Slope, df$Intercept), aes(x, y), 
             color = 'green3') +
  coord_cartesian(xlim = c(-6, 6), ylim = c(-6, 6))

该函数innermost定义如下：

innermost <- function(slopes, intercepts) {
  
  meetings <- function(slopes, intercepts) {
    meets_at <- function(i1, s1, i2, s2) {
      ifelse(s1 - s2 == 0, NA, (i2 - i1)/(s1 - s2))
    }    
    xvals <- outer(seq_along(slopes), seq_along(slopes), function(i, j) {
      meets_at(intercepts[i], slopes[i], intercepts[j], slopes[j])
    })
    
    yvals <- outer(seq_along(slopes), seq_along(slopes), function(i, j) {
      intercepts + slopes *
        meets_at(intercepts[i], slopes[i], intercepts[j], slopes[j])
    })
    
    cbind(x = xvals[lower.tri(xvals)], y = yvals[lower.tri(yvals)])
    
  xy <- meetings(slopes, intercepts)
  is_cut <- function(x, y, slopes, intercepts) {
    d <- sqrt(x^2 + y^2)
    slope <- y / x
    xvals <- intercepts / (slope - slopes)
    yvals <- xvals * slopes + intercepts
    ds <- sqrt(xvals^2 + yvals^2)
    any(d - ds > 1e-6 & sign(xvals) == sign(x) & sign(yvals) == sign(y))
  }
  xy <- xy[sapply(seq(nrow(xy)), function(i) {
    !is_cut(xy[i, 1], xy[i, 2], slopes, intercepts)
  }),]
  xy <- xy[order(atan2(xy[,2], xy[,1])),]
  as.data.frame(rbind(xy, xy[1,]))
}

解释

首先，获得两条直线的交点很简单。直线的公式由下式给出，其中m是斜率，c是截距。所以如果我们有两条不同的直线由方程给出和那么他们必须在哪里见面

重新排列我们得到

或者

这意味着两条线相交

也就是说，如果第一行有intercept1and ，第二行有and ，那么我们可以使用这个简单的函数找到交汇点的 x 值：slope1intercept2slope2

meets_at <- function(intercept1, slope1, intercept2, slope2) {
  ifelse(slope1 - slope2 == 0, NA, (intercept2 - intercept1)/(slope1 - slope2))
}

请注意，如果线是平行的，即slope1 - slope2 == 0，它们将没有唯一的交点，并且此函数应返回NA

我们可以在所有线对上使用此函数来获取所有交点outer：

meetings <- function(slopes, intercepts) {
  
  xvals <- outer(seq_along(slopes), seq_along(slopes), function(i, j) {
            meets_at(intercepts[i], slopes[i], intercepts[j], slopes[j])
           })
  
  yvals <- outer(seq_along(slopes), seq_along(slopes), function(i, j) {
            intercepts + slopes *
            meets_at(intercepts[i], slopes[i], intercepts[j], slopes[j])
          })
  
  cbind(x = xvals[lower.tri(xvals)], y = yvals[lower.tri(yvals)])
}

例如，如果我们绘制这些点，我们将看到绘制的线的所有交点：

plot(seq(-6, 6), seq(-6, 6), type = 'n')
for(i in seq(nrow(df))) abline(a = df$Intercept[i], b = df$Slope[i])
xy <- meetings(df$Slope, df$Intercept)
points(xy[,1], xy[,2], col = 'red')

这仍然给我们留下了仅找到形成所需多边形轮廓的最内部点的问题。为此，请注意当我们从原点 (0, 0) 到上图中的每个交点绘制一条线时会发生什么：

xy <- xy[abs(xy[,1]) < 6 & abs(xy[,2] < 6),] # Remove intersections outside plot

for(i in seq(nrow(df))) {
  segments(x0 = 0, y0 = 0, x1 = xy[,1], y1 = xy[,2], col = '#0000ff20')
}

请注意，从原点到最里面的顶点（我们要保留的顶点）的蓝线不与任何其他线交叉。换句话说，我们要丢弃的顶点是那些无法在不与另一条线相交的情况下绘制到原点的直线的顶点。

因此，我们可以计算出是否有任何线穿过连接原点和交点的线段，并仅保留那些未交叉的线。

我们还需要按圆周顺序排列最终点。这是通过计算 x 轴与绘制到该点的直线之间的角度的反正切来完成的，该角度正好为atan2(y, x)。这给了我们一个介于 -pi 和 pi 之间的数字，可以从 9 点钟位置开始逆时针对点进行排序：

innermost <- function(slopes, intercepts) {
  xy <- meetings(slopes, intercepts)
  is_cut <- function(x, y, slopes, intercepts) {
    d <- sqrt(x^2 + y^2)
    slope <- y / x
    xvals <- intercepts / (slope - slopes)
    yvals <- xvals * slopes + intercepts
    ds <- sqrt(xvals^2 + yvals^2)
    any(d - ds > 1e-6 & sign(xvals) == sign(x) & sign(yvals) == sign(y))
  }
  xy <- xy[sapply(seq(nrow(xy)), function(i) {
    !is_cut(xy[i, 1], xy[i, 2], slopes, intercepts)
  }),]
  xy <- xy[order(atan2(xy[,2], xy[,1])),]
  as.data.frame(rbind(xy, xy[1,]))
}

我们可以使用上面的函数来找到由一组线创建的最里面的多边形。在 R 基础上绘图，我们可以这样做：

plot(seq(-6, 6), seq(-6, 6), type = 'n')
for(i in seq(nrow(df))) abline(a = df$Intercept[i], b = df$Slope[i])
xy <- innermost(df$Slope, df$Intercept)
points(xy$x, xy$y, col = 'red')
polygon(xy$x, xy$y, col = 'gray')