15

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)

enter image description here

Henrik
  • 65,555
  • 14
  • 143
  • 159
Mercurial
  • 195
  • 6
  • 1
    Are you trying to find the closest value of `x` to intersecting the circle, or would you like to approximate as closely as possible the solution to `f(x, y) = circle(x, y)` – MichaelChirico Apr 29 '18 at 15:08
  • 1
    you should be aware that `plotrix::draw.circle()` draws a circle whose x and y dimensions depend on the scaling of the x and y axes in the plots; in your case the maxima of the circle in the x direction will be (1,5) [x-radius of 2], but in the y dimension it looks like the maximum is only about 1.2. Is it safe to assume when solving the problem that x and y axis scales are indeed in matching units? (e.g. see `MASS::eqscplot`) – Ben Bolker Apr 29 '18 at 15:21
  • I would like to approximate the intersection coordinates as closely as possible. One way to think about these coordinates is as a vertex of a right-angle triangle formed along the x-axis. Thanks! – Mercurial Apr 29 '18 at 15:31
  • @Ben Bolker Thank you! This was very useful, indeed. I typically plot with ggplot, and had to account for the problem you describe there. It's great to know how to do it with draw.circle() – Mercurial May 01 '18 at 14:38

2 Answers2

10

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)

enter image description here

Onyambu
  • 67,392
  • 3
  • 24
  • 53
  • instead of the optim you can also do `library(nleqslv);nleqslv(1,g)[[1]]` – Onyambu Apr 29 '18 at 17:24
  • This is fantastic! Thank you. One question, how would you change the circle() function to include a set of (x, y) coordinates for the circle (other than (3, 0), as in my example? – Mercurial Apr 29 '18 at 18:07
  • the `(3,0)` is the centre of the circle. Its not the x,y coordinates... if you run `circleplot(3,0,2)` it will give you x and y coordinates for a circle centered at 3,0 with radius 2 – Onyambu Apr 29 '18 at 18:11
  • Thank you so much for your help! I found your comments very helpful. – Mercurial Apr 29 '18 at 20:42
  • Thanks for your answer. My question about coordinates was about the ```circle()``` function above, not ```circleplot()```. How would you formulate your ```circle()``` function to include (x,y) coordinates at the center, e.g., (x=3, y=1.5). Thanks! – Mercurial Apr 30 '18 at 00:28
  • Am so sorry I do not understand your question. Buy I hope you did find the best solution to your problem – Onyambu Apr 30 '18 at 00:49
  • @ Onyambu perhaps I can rephrase my question. Could you explain the terms of this expression: ```sqrt(4-(x-3)^2)``` in the ```circle()``` function. Thanks. – Mercurial Apr 30 '18 at 01:05
  • `y=sqrt(r^2-(x-h)^2)+k` whre `(h,k)` is your centre – Onyambu Apr 30 '18 at 01:12
  • I wonder if I could ask your help again. The solution you offered ```nleqslv(1,g)[[1]]``` works for some of the data, but fails for others. For example, if you change the radius to 0.2, there is apparently no solution available. Could you take another look and help me figure it out? Thanks! – Mercurial May 02 '18 at 12:56
  • Yes, I have been using the solution provided by Henrik and it works great. However, I'd like to be able to use the nleqslv method, because, it should work, in principle, as you suggested. It's just frustrating not to be able to use it for some of the data. Thanks! – Mercurial May 02 '18 at 13:31
  • I believe maybe doing a research on how touse it will be very helpful. I completely understand how you would like to learn new functions and even use them... but that is beyond the scope of this question... i believe if you do a research, you will be able to understand the use of the function... my aim was to try and write a R base solution to your problem and i hope it did give an insight to all... am happy to help – Onyambu May 02 '18 at 13:57
  • you can use `uniroot` from base R since this solves for a unique value.. eg `uniroot(g,c(1,2))[[1]]` – Onyambu May 02 '18 at 23:49
  • Thank you so much. This is very helpful. However, for values of R < 1, e.g., 0.6, the circle equation produces ```NaN``` so the ```uniroot()``` will return an error. Do you know why that is? Thanks. I feel like I'm getting closer to the solution, thanks to your help. – Mercurial May 03 '18 at 13:34
  • But if you have a circle with radius 2 and centered at 3...the lowest it can go is 3-2=1.. how cas 0.6 be within the circle?? – Onyambu May 03 '18 at 14:14
  • yes, of course. But let's try to find the point of intersection between the curve and the circle, centered at 3, but with the radius of 0.6. What would the ```circle()``` function look like then? – Mercurial May 03 '18 at 14:26
  • Just do circle=function(x,r,cent){ if(r^2<(x-cent[1])^2) return(NA)# Ensure it is limited within the radius sqrt(r^2-(x-cent[1])^2)+cent[2] }. Now you can put any radius and centre – Onyambu May 03 '18 at 14:31
  • That's it. You solved it! Thank you so much. I was making the mistake of feeding unrealistic values to the function. When I apply my actual data, rather than the unrealistic test data, to your solution, it seems to work. – Mercurial May 03 '18 at 15:13
8

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)

enter image description here

Henrik
  • 65,555
  • 14
  • 143
  • 159
  • Thank you! This solution works very well for my data. I another question, as far as the theory behind the method is concerned, is it true that the function st_intersection() will always find a solution as long as there exists an intersection between the circle and the line? Or are there potential cases where the function would fail, despite there being an intersection? Thanks! – Mercurial Apr 29 '18 at 19:48
  • Not that I know of. Also check [the vignette](https://cran.r-project.org/web/packages/sf/vignettes/sf3.html#operations-returning-a-geometry). I _think_ the (more thorough) description in [the PostGIS manual](http://postgis.net/docs/manual-2.1/ST_Intersection.html) is valid here as well. See also [here](http://workshops.boundlessgeo.com/postgis-intro/geometry_returning.html#st-intersection) – Henrik Apr 29 '18 at 20:32
  • 1
    Thank you so much! I just ran my data analysis. So far so good. Now, I check, and re-check :). I owe you one. – Mercurial Apr 29 '18 at 20:40