如何找到曲线和圆的交点？

Question

我有一条曲线，是由实测数据推导而来的，并且我可以得到一个合理的模型。我需要找到曲线与已知圆心和半径的圆相交的点(x,y)。以下代码展示了这个问题。

x <- c(0.05, 0.20, 0.35, 0.50, 0.65, 0.80, 0.95, 
   1.10, 1.25, 1.40, 1.55, 1.70, 1.85, 2.00, 
   2.15, 2.30, 2.45, 2.60, 2.75, 2.90, 3.05)

y <- c(1.52, 1.44, 1.38, 1.31, 1.23, 1.15, 1.06,
   0.96, 0.86, 0.76, 0.68, 0.61, 0.54, 0.47, 
   0.41, 0.36, 0.32, 0.29, 0.27, 0.26, 0.26)

fit <- loess(y ~ x, control = loess.control(surface = "direct"))
newx <- data.frame(x = seq(0, 3, 0.01))
fitline <- predict(fit, newdata = newx)
est <- data.frame(newx, fitline)

plot(x, y, type = "o",lwd = 2)
lines(est, col = "blue", lwd = 2)

library(plotrix)
draw.circle(x = 3, y = 0, radius = 2, nv = 1000, lty = 1, lwd = 1)

- Mercurial

---

你是想找到最接近与圆相交的x值，还是尽可能精确地逼近解决方案f(x, y) = circle(x, y)？ - MichaelChirico

请注意，plotrix::draw.circle()绘制的圆的x和y尺寸取决于绘图中x和y轴的缩放；在您的情况下，圆在x方向上的最大值将是（1,5）[x半径为2]，但在y方向上，最大值似乎只有约1.2。在解决问题时，可以安全地假设x和y轴比例确实是匹配的单位吗？（例如，请参见MASS::eqscplot） - Ben Bolker

我想尽可能精确地近似交点坐标。思考这些坐标的一种方式是将其视为沿x轴形成的直角三角形的顶点。谢谢！ - Mercurial

@Ben Bolker 谢谢！这确实非常有用。我通常使用ggplot绘图，并且必须考虑您在那里描述的问题。知道如何使用draw.circle()解决它真是太好了。 - Mercurial

Answer 1

为了获得交点，可以使用R中的optim函数来实现：

circle=function(x){
  if(4<(x-3)^2) return(NA)# Ensure it is limited within the radius
  sqrt(4-(x-3)^2)
}
fun=function(x)predict(fit,data.frame(x=x))  
g=function(x)(circle(x)-fun(x))# We need to set this to zero. Ie solve this
sol1=optimise(function(x)abs(g(x)),c(1,5))$min
 [1] 1.208466

因此，这两个函数在x=1.208466处应该评估为相同的值。

为了使其更加精确，您可以使用optim函数：

sol2= optim(1,function(x)abs(g(x)),g,method="Brent",upper=5,lower=1)$par
 [1] 1.208473

现在您可以进行评估：

circle(sol1)
[1] 0.889047
fun(sol1)
        1 
0.8890654 
circle(sol2)
[1] 0.889061
fun(sol2)
       1 
0.889061 

从上面可以看出，解决方案2非常接近...

在图表上绘制这个点将是具有挑战性的，因为 draw.circle 函数会根据zxes的比例绘制圆圈...因此每次都会根据绘图区域的大小而改变。

如果您要编写自己的圆函数：

circleplot=function(x,y,r){
  theta=seq(0,2*pi,length.out = 150)
  cbind(x+r*cos(theta),y+r*sin(theta))
}

然后你可以这样做：

plot(x, y, type = "o",lwd = 2)
lines(est, col = "blue", lwd = 2)
lines(circleplot(3,0,2))
abline(v=sol2,col=2) 
points(sol2,fun(sol2),col=2,pch=16)

- onyambu

---

太好了！谢谢。有一个问题，您如何更改circle()函数以包括圆的一组(x,y)坐标（而不是像我的示例中那样的(3,0)）？ - Mercurial

@Onyambu也许我可以重新表达我的问题。您能解释一下在circle()函数中这个表达式的术语吗：sqrt(4-(x-3)^2)。谢谢。 - Mercurial

y=sqrt(r^2-(x-h)^2)+k 其中 (h,k) 是你的中心位置 - Onyambu

我想问你能否再次帮助我。你提供的解决方案 nleqslv(1,g)[[1]] 对一些数据有效，但对其他数据则失败了。例如，如果将半径更改为0.2，则显然没有可用的解决方案。你能否重新审视一下并帮助我找出问题所在？谢谢！ - Mercurial

我认为，对于如何使用它进行研究可能会非常有帮助。我完全理解您想学习新功能甚至使用它们的愿望...但这超出了这个问题的范围...我相信如果您进行研究，就能够理解该函数的用途...我的目标是尝试编写一个R基础解决方案来解决您的问题，希望它确实给大家提供了一些见解...很高兴能够帮助。 - Onyambu

Answer 2

使用sf包中的函数很容易找到交叉点。

计算圆的值（灵感来自这个答案和@Onyambu所做的）

circ <- function(xc = 0, yc = 0, r = 1, n = 100){
  v <- seq(0, 2 * pi, len = n)
  cbind(x = xc + r * cos(v),
        y = yc + r * sin(v))
}

m <- circ(xc = 3, yc = 0, r = 2)

将预测值和圆值转换为“简单要素” (LINESTRING)，并找到它们的交点 (一个 POINT):

library(sf)
int <- st_intersection(st_linestring(as.matrix(est)),
                       st_linestring(m))
int
# POINT (1.2091 0.8886608)

将交点添加到您的图中：

plot(x, y, type = "o", lwd = 2)
lines(est, col = "blue", lwd = 2)
lines(m)
points(int[1], int[2], col = "red", pch = 19)

- Henrik

---

谢谢！这个解决方案对我的数据非常有效。我还有一个问题，就是关于这种方法背后的理论，是否真的如此，即只要圆和直线之间存在交点，函数st_intersection()就总能找到解决方案？或者说，尽管存在交点，函数仍可能失败的潜在情况吗？谢谢！ - Mercurial

我不知道。还可以查看vignette。 我认为PostGIS手册中的（更详细的）描述在这里也是有效的。另请参见此处。 - Henrik

非常感谢！我刚刚运行了我的数据分析。目前为止一切都很好。现在，我会检查和重新检查：）。我欠你一个人情。 - Mercurial