使用求解/优化根据两点（已知半径）确定圆心答案

【问题标题】：determine circle center based on two points (radius known) with solve/optim使用求解/优化根据两点（已知半径）确定圆心
【发布时间】：2012-08-29 04:38:17
【问题描述】：

我有一对点，我想找到一个由这两个点确定的已知 r 的圆。我将在模拟中使用它，x 和 y 的可能空间有边界（比如 -200、200 的框）。

It is known那个平方的半径是

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

附录

感谢大家提供宝贵的 cmets 和代码。我会为那些用代码称赞他们的答案的发帖者提供答案的时间。

    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

【问题讨论】：

点是否在圆周上？
可以手动求解方程组并使用公式
@James，是的，这些点位于圆周上的某个位置。我已经更新了显示结果的答案。

标签： r optimization geometry

【解决方案1】：

这是一个答案的骨架，如果我以后有时间我会充实它们。如果你跟着文字一起画，这应该很容易理解，对不起，我这台电脑上没有合适的软件为你画图。

抛开退化的情况，即点相同（无限解）或相距太远而无法位于所选半径的同一圆上（无解）。

标记点X 和Y 以及两个圆的未知中心点c1 和c2。 c1和c2位于XY的垂直平分线上；调用这条线c1c2，在这个阶段，我们不知道c1和c2的位置的所有细节是无关紧要的。

所以，计算线c1c2 的方程。它通过XY 的中点（称为此点Z）并且斜率等于XY 的负倒数。现在你有了c1c2 的等式（或者如果这些骨头上有肉，你就会得到）。

现在构造一个三角形，从一点到直线与其垂直平分线的交点和圆的中心点（比如XZc1）。您仍然不知道 c1 的确切位置，但这从未阻止任何人绘制几何图形。你有一个直角三角形的两条边长已知（XZ 和Xc1），所以很容易找到Zc1。对另一个三角形和圆心重复该过程。

当然，这种方式与 OP 最初的方式有很大不同，可能没有吸引力。

【讨论】：

【解决方案2】：

一些警告要摆脱，但这应该让你开始。可能存在性能问题，因此使用基本几何完全解决它可能是更好的方法。

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)

【讨论】：

【解决方案3】：

我希望你知道一些基本的几何图形，因为很遗憾我不会画它。

垂直平分线是穿过 A 和 B 的圆的每个中点所在的线。

现在你有了 AB 和 r 的中点，所以你可以用点 A、AB 的中点和未知的圆的中点画一个直角三角形。

现在使用毕达哥拉斯定理得到AB中点到圆中点的距离，从这里开始计算圆的位置应该不难，使用基本的sin/cos组合。

【讨论】：

【解决方案4】：

不需要数值方程求解。只是公式：

您知道，由于 A 点和 B 点都位于圆上，因此从每个点到给定中心的距离就是半径 r。
以两个已知点的弦为底，第三个点为圆心，形成一个等腰三角形。
将 A 和 B 之间的三角形一分为二，得到一个直角三角形。
http://mathworld.wolfram.com/IsoscelesTriangle.html 为您提供基于底长和半径的高度。
按照 AB 弦 (See this SO Answer) 的法线，计算从该点在每个方向上刚刚计算的高度的距离。

【讨论】：

【解决方案5】：

以下代码将为您提供两个所需圆的中心点。现在没有时间对此发表评论或将结果转换为 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)

【讨论】：

【解决方案6】：

这是其他人都提到的解决问题的基本几何方法。我使用 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)

【讨论】：

【解决方案7】：

按照@PhilH 的解决方案，仅在 R 中使用三角函数：

radius=40

在半径上画出原点

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

找到ab的中点c（这也是de两个圆心的中点）

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)

求弦的长度和角度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

计算从c到两个圆心的直线的长度和角度

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

计算并绘制两个中心d 和e 从垂直于ab 的位置和长度：

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

演出：

【讨论】：