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

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)

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)

Answer 2

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

计算圆值（灵感来自this answer 和as done by @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)