How to find an intersection of curve and circle?

Question

I have a curve, derived from empirical data, and I can obtain a reasonable model of it. I need to identify a point (x, y) where the curve intersects a circle of known center and radius. The following code illustrates the question.

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

To obtain the point of intersection we can use the optim function in r to do so:

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

Thus the two functions should evaluate to the same value at x=1.208466..

To make it even more precise, you can use the optim function:

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

Now you can evaluate:

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

From the above, you can tell that solution 2 is very close..

Plotting this point on the graph will be challenging since the draw.circle function draws circles in proportionality with the zxes.. Thus changing everytime depending on how big the plot region is.

If you were to write your own circle function:

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

Then you can do:

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

It's straightforward to find the intersection using functions from the sf package.

Calculate the circle values (inspired by this answer and 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)

Convert the predicted values and the circle values to "simple features" (LINESTRING), and find their intersection (a POINT):

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

Add the intersection to your plot:

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