determine circle center based on two points (radius known) with solve/optim

Question

I have a pair of points and I would like to find a circles of known r that are determined by these two points. I will be using this in a simulation and possible space for x and y have boundaries (say a box of -200, 200).

It is known that square of radius is

(x-x1)^2 + (y-y1)^2 = r^2
(x-x2)^2 + (y-y2)^2 = r^2

I would now like to solve this non-linear system of equations to get two potential circle centers. I tried using package BB. Here is my feeble attempt which gives only one point. What I would like to get is both possible points. Any pointers into right direction will be met with complimentary beer on first possible occasion.

library(BB)
known.pair <- structure(c(-46.9531139599816, -62.1874917150412, 25.9011462171242, 
16.7441676243879), .Dim = c(2L, 2L), .Dimnames = list(NULL, c("x", 
"y")))

getPoints <- function(ps, r, tr) {
    # get parameters
     x <- ps[1]
     y <- ps[2]
                
     # known coordinates of two points
     x1 <- tr[1, 1]
     y1 <- tr[1, 2]
     x2 <- tr[2, 1]
     y2 <- tr[2, 2]
                
     out <- rep(NA, 2)
     out[1] <- (x-x1)^2 + (y-y1)^2 - r^2
     out[2] <- (x-x2)^2 + (y-y2)^2 - r^2
     out
}

slvd <- BBsolve(par = c(0, 0),
                fn = getPoints,
                method = "L-BFGS-B",
                tr = known.pair,
                r = 40
                )

Graphically you can see this with the following code, but you will need some extra packages.

library(sp)
library(rgeos)
plot(0,0, xlim = c(-200, 200), ylim = c(-200, 200), type = "n", asp = 1)
points(known.pair)
found.pt <- SpatialPoints(matrix(slvd$par, nrow = 1))
plot(gBuffer(found.pt, width = 40), add = T)

enter image description here

ADDENDUM

Thank you all for your valuable comments and code. I provide timings for answers by posters who complimented their answers with code.

    test replications elapsed relative user.self sys.self user.child sys.child
4   alex          100    0.00       NA      0.00        0         NA        NA
2  dason          100    0.01       NA      0.02        0         NA        NA
3   josh          100    0.01       NA      0.02        0         NA        NA
1 roland          100    0.15       NA      0.14        0         NA        NA

It is possible to solve equation system by hands and use formulas — MBo, Sep 04 '12 at 13:49
@James, yes the points lie somewhere on the circumference. I've updated my answer that shows the result. — Roman Luštrik, Sep 04 '12 at 13:50

Josh O'Brien · Answer 1 · 2012-09-04T17:34:16.123

The following code will get you the points at the centers of the two desired circles. No time right now to comment this up or convert the results to Spatial* objects, but this should give you a good start.

First, here's an ASCII-art diagram to introduce point names. k and K are the known points, B is a point on the horizontal drawn through k, and C1 and C2 are the centers of the circles you are after:

        C2





                            K


                  k----------------------B






                                       C1

Now the code:

# Example inputs
r <- 40
known.pair <- structure(c(-46.9531139599816, -62.1874917150412, 
25.9011462171242, 16.7441676243879), .Dim = c(2L, 2L), 
.Dimnames = list(NULL, c("x", "y")))

## Distance and angle (/_KkB) between the two known points
d1 <- sqrt(sum(diff(known.pair)^2))
theta1 <- atan(do.call("/", as.list(rev(diff(known.pair)))))

## Calculate magnitude of /_KkC1 and /_KkC2
theta2 <- acos((d1/2)/r)

## Find center of one circle (using /_BkC1)
dx1 <- cos(theta1 + theta2)*r
dy1 <- sin(theta1 + theta2)*r
p1 <- known.pair[2,] + c(dx1, dy1)

## Find center of other circle (using /_BkC2)
dx2 <- cos(theta1 - theta2)*r
dy2 <- sin(theta1 - theta2)*r
p2 <- known.pair[2,] + c(dx2, dy2)

## Showing that it worked
library(sp)
library(rgeos)
plot(0,0, xlim = c(-200, 200), ylim = c(-200, 200), type = "n", asp = 1)
points(known.pair)
found.pt <- SpatialPoints(matrix(slvd$par, nrow = 1))
points(p1[1], p1[2], col="blue", pch=16)
points(p2[1], p2[2], col="green", pch=16)

enter image description here

score 4 · Accepted Answer · answered Sep 04 '12 at 15:53

This is the basic geometric way of going about solving it that everybody else is mentioning. I use polyroot to get the roots of the resulting quadratic equation but you could easily just use the quadratic equation directly.

# x is a vector containing the two x coordinates
# y is a vector containing the two y coordinates
# R is a scalar for the desired radius
findCenter <- function(x, y, R){
    dy <- diff(y)
    dx <- diff(x)
    # The radius needs to be at least as large as half the distance
    # between the two points of interest
    minrad <- (1/2)*sqrt(dx^2 + dy^2)
    if(R < minrad){
        stop("Specified radius can't be achieved with this data")
    }

    # I used a parametric equation to create the line going through
    # the mean of the two points that is perpendicular to the line
    # connecting the two points
    # 
    # f(k) = ((x1+x2)/2, (y1+y2)/2) + k*(y2-y1, x1-x2)
    # That is the vector equation for our line.  Then we can
    # for any given value of k calculate the radius of the circle
    # since we have the center and a value for a point on the
    # edge of the circle.  Squaring the radius, subtracting R^2,
    # and equating to 0 gives us the value of t to get a circle
    # with the desired radius.  The following are the coefficients
    # we get from doing that
    A <- (dy^2 + dx^2)
    B <- 0
    C <- (1/4)*(dx^2 + dy^2) - R^2

    # We could just solve the quadratic equation but eh... polyroot is good enough
    k <- as.numeric(polyroot(c(C, B, A)))

    # Now we just plug our solution in to get the centers
    # of the circles that meet our specifications
    mn <- c(mean(x), mean(y))
    ans <- rbind(mn + k[1]*c(dy, -dx),
                 mn + k[2]*c(dy, -dx))

    colnames(ans) = c("x", "y")

    ans
}

findCenter(c(-2, 0), c(1, 1), 3)

Alex Brown · Answer 3 · 2012-09-04T18:17:15.400

Following @PhilH's solution, just using trigonometry in R:

radius=40

Draw the original points on the radius

plot(known.pair,xlim=100*c(-1,1),ylim=100*c(-1,1),asp=1,pch=c("a","b"),cex=0.8)

Find the midpoint c of ab (which is also the midpoint of de the two circle centers)

AB.bisect=known.pair[2,,drop=F]/2+known.pair[1,,drop=F]/2
C=AB.bisect
points(AB.bisect,pch="c",cex=0.5)

Find the length and angle of the chord ab

AB.vector=known.pair[2,,drop=F]-known.pair[1,,drop=F]
AB.len=sqrt(sum(AB.vector^2))
AB.angle=atan2(AB.vector[2],AB.vector[1])
names(AB.angle)<-NULL

Calculate the length and angle of the line from c to the centers of the two circles

CD.len=sqrt(diff(c(AB.len/2,radius)^2))
CD.angle=AB.angle-pi/2

Calculate and plot the position of the two centers d and e from the perpendicular to ab and the length:

center1=C+CD.len*c(x=cos(CD.angle),y=sin(CD.angle))
center2=C-CD.len*c(x=cos(CD.angle),y=sin(CD.angle))
points(center1[1],center1[2],col="blue",cex=0.8,pch="d")
points(center2[1],center2[2],col="blue",cex=0.8,pch="e")

Shows:

enter image description here

score 3 · Answer 4 · edited May 23 '17 at 12:01

No numerical equation solving required. Just formulae:

You know that since both points A and B lie on the circle, the distance from each to a given centre is the radius r.
Form an isosceles triangle with the chord of the two known points at the base and the third point at the circle centre.
Bisect the triangle midway between A and B, giving you a right-angle triangle.
http://mathworld.wolfram.com/IsoscelesTriangle.html gives you the height in terms of the base length and the radius.
Follow the normal to the AB chord (See this SO Answer) for a distance of the height just calculated in each direction from the point.

score 1 · Answer 5 · answered Sep 04 '12 at 14:18

Here are the bones of an answer, if I have time later I'll flesh them out. This should be easy enough to follow if you draw along with the words, sorry I don't have the right software on this computer to draw the picture for you.

Leave aside degenerate cases where the points are identical (infinite solutions) or too far apart to lie on the same circle of your chosen radius (no solutions).

Label the points X and Y and the unknown centre points of the 2 circles c1 and c2. c1 and c2 lie on the perpendicular bisector of XY; call this line c1c2, at this stage it's immaterial that we don't know all the details of the locations of c1 and c2.

So, figure out the equation of line c1c2. It passes through the half-way point of XY (call this point Z) and has slope equal to the negative reciprocal of XY. Now you have the equation of c1c2 (or you would if there was any flesh on these bones).

Now construct the triangle from one point to the intersection of the line and its perpendicular bisector and the centre point of a circle (say XZc1). You still don't know exactly where c1 is but that never stopped anyone sketching the geometry. You have a right triangle with two side lengths known (XZ and Xc1), so it's easy to find Zc1. Repeat the process for the other triangle and circle centre.

Of course, this approach is quite different from OP's initial approach and may not appeal.

Roland · Answer 6 · 2012-09-04T14:40:27.517

Some warnings to get rid of, but this should get you started. There might be a performance issue, so solving it completely with basic geometry could be a better approach.

known.pair <- structure(c(-46.9531139599816, -62.1874917150412, 25.9011462171242, 
                          16.7441676243879), .Dim = c(2L, 2L), .Dimnames = list(NULL, c("x", 
                                                                                        "y")))

findCenter <- function(p,r) {
  yplus <- function(y) {
    ((p[1,1]+sqrt(r^2-(y-p[1,2])^2)-p[2,1])^2+(y-p[2,2])^2-r^2)^2
  }


 yp <- optimize(yplus,interval=c(min(p[,2]-r),max(p[,2]+r)))$minimum
 xp <- p[1,1]+sqrt(r^2-(yp-p[1,2])^2)  
 cp <- c(xp,yp)
 names(cp)<-c("x","y") 

 yminus <- function(y) {
   ((p[1,1]-sqrt(r^2-(y-p[1,2])^2)-p[2,1])^2+(y-p[2,2])^2-r^2)^2
 }


 ym <- optimize(yminus,interval=c(min(p[,2]-r),max(p[,2]+r)))$minimum
 xm <- p[1,1]-sqrt(r^2-(ym-p[1,2])^2)  
 cm <- c(xm,ym)
 names(cm)<-c("x","y")


 list(c1=cp,c2=cm)
}

cent <- findCenter(known.pair,40)

score 0 · Answer 7 · answered Sep 04 '12 at 14:31

I hope you know some basic geometry, because I cannot draw it unfortunately.

The perpendicular bisector is the line where every middle point of a circle which crosses both A and B lays.

Now you have the middle of AB and r, so you can draw a right triangle with the point A, the middle of AB and the unknown middle point of the circle.

Now use the pythagoras' theorem to get the distance from the middle point of AB to the middle point of the circle, and calculate the position of the circle shouldn't be hard from here, using basic sin/cos combinations.

determine circle center based on two points (radius known) with solve/optim

7 Answers7