使用solve/optim基于两个点(已知半径)确定圆心

Question

我有一对积分,我想找到一个由这两点决定的已知r的圆圈.我将在模拟和可能的空间中使用它x并且y有边界(比如一个-200,200的盒子).

众所周知,半径的平方是

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

我现在想解决这个非线性方程组,得到两个潜在的圆心.我试过用包BB.这是我的微弱尝试,只给出了一点.我想得到的是两个可能的要点.任何指向正确方向的指针都将在第一时间获得免费啤酒.

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
                )

以图形方式,您可以使用以下代码看到这一点,但您需要一些额外的包.

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)

附录

谢谢大家的宝贵意见和代码.我提供了海报的答案,他们用代码称赞他们的答案.

    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

Answer 1

以下代码将为您提供两个所需圆圈的中心点.现在没时间对此进行评论或将结果转换为Spatial*对象,但这应该会给你一个良好的开端.

首先,这是一个引入点名称的ASCII艺术图.k并且K是已知的点,B是水平绘制的一个点k,C1并且C2是你所追求的圆圈的中心:

        C2





                            K


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






                                       C1

现在的代码:

# 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)

Answer 2

这是其他人都提到的解决该问题的基本几何方法。我使用 polyroot 来获取所得二次方程的根，但您可以轻松地直接使用二次方程。

# 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)