determinar el centro del círculo basado en dos puntos (radio conocido) con resolver/optim

Question

Tengo un par de puntos y me gustaría encontrar un círculo de r conocido que esté determinado por estos dos puntos. Voy a usar esto en una simulación y espacio posible para x y y tienen límites (por ejemplo, un cuadro de -200, 200).

It is known ese cuadrado del radio es

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

ahora me gustaría resolver este sistema no lineal de ecuaciones para obtener dos centros potenciales círculo. Intenté usar el paquete BB. Aquí está mi intento débil que da solo un punto. Lo que me gustaría obtener son ambos puntos posibles. Cualquier puntero en la dirección correcta se encontrará con cerveza de cortesía en la primera ocasión posible.

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 
       ) 

Gráficamente, puede ver esto con el siguiente código, pero necesitará algunos paquetes adicionales.

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) 

ADENDA

Gracias a todos por sus valiosos comentarios y código. Proporciono tiempos para las respuestas de los carteles que felicitaron sus respuestas con el código.

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

Esta es la forma geométrica básica de resolverlo que todos los demás mencionan. Yo uso polyroot para obtener las raíces de la ecuación cuadrática resultante, pero puedes usar la ecuación cuadrática directamente.

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

Answer 2

Aquí están los huesos de una respuesta, si tengo tiempo más tarde los voy a procesar. Esto debería ser lo suficientemente fácil de seguir si dibujas junto con las palabras, lo siento, no tengo el software adecuado en esta computadora para dibujar la imagen para ti.

Deje a un lado los casos degenerados donde los puntos son idénticos (soluciones infinitas) o demasiado separados para estar en el mismo círculo del radio elegido (sin soluciones).

Etiquete los puntos X y Y y los puntos centrales desconocidos de los 2 círculos c1 y c2. c1 y c2 se encuentran en la bisectriz perpendicular de XY; llame a esta línea c1c2, en este momento es indiferente que no conozcamos todos los detalles de las ubicaciones de c1 y c2.

Por lo tanto, averiguar la ecuación de la línea c1c2. Pasa por el punto medio de XY (llame a este punto Z) y tiene una pendiente igual al recíproco negativo de XY. Ahora tiene la ecuación de c1c2 (o lo haría si hubiera algo de carne en estos huesos).

Ahora construye el triángulo desde un punto hasta la intersección de la línea y su bisectriz perpendicular y el punto central de un círculo (digamos XZc1). Aún no sabes exactamente dónde está c1, pero eso nunca detuvo a nadie dibujando la geometría. Tiene un triángulo rectángulo con dos longitudes de lado conocidas (XZ y Xc1), por lo que es fácil encontrar Zc1. Repita el proceso para el otro triángulo y centro del círculo.

Por supuesto, este enfoque es bastante diferente del enfoque inicial de OP y puede no ser atractivo.

Answer 3

Algunas advertencias para deshacerse, pero esto debería comenzar. Puede haber un problema de rendimiento, por lo que resolverlo completamente con geometría básica podría ser un mejor enfoque.

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) 

Answer 4

Espero que sepas algo de geometría básica, porque no puedo dibujarla desafortunadamente.

La bisectriz perpendicular es la línea donde se ubica cada punto medio de un círculo que cruza tanto A como B.

Ahora tiene el medio de AB y r, por lo que puede dibujar un triángulo rectángulo con el punto A, el centro de AB y el punto medio desconocido del círculo.

Ahora use el teorema de Pitágoras para obtener la distancia desde el punto medio de AB hasta el punto medio del círculo, y calcular la posición del círculo no debería ser difícil desde aquí, usando combinaciones básicas de sin/cos.

Answer 5

No es necesario resolver ecuaciones numéricas. Solo fórmulas:

1. Sabe que dado que ambos puntos A y B se encuentran en el círculo, la distancia de cada uno a un centro dado es el radio r.
2. Forma un triángulo isósceles con el acorde de los dos puntos conocidos en la base y el tercer punto en el centro del círculo.
3. Biseca el triángulo a mitad de camino entre A y B, que te da un triángulo en ángulo recto.
4. http://mathworld.wolfram.com/IsoscelesTriangle.html le da la altura en términos de la longitud de la base y el radio.
5. Siga la normal al acorde AB (See this SO Answer) para una distancia de la altura recién calculada en cada dirección desde el punto.

Answer 6

El siguiente código obtendrá los puntos en los centros de los dos círculos deseados. No hay tiempo ahora para comentar esto o convertir los resultados a objetos Spatial*, pero esto debería darle un buen comienzo.

Primero, aquí hay un diagrama ASCII-art para introducir nombres de puntos. k y K son los puntos conocidos, B es un punto sobre la horizontal trazada a través k y C1 y C2 son los centros de los círculos que está buscando:

 C2 





          K 


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






             C1 

Ahora el código:

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

Siguiendo @ solución de PhilH, simplemente utilizando la trigonometría en R:

radius=40 

Dibujar los puntos originales en el radio

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

Encuentra el punto medio c de ab (que es también el punto medio de de los dos centros del círculo)

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) 

Encuentra la longitud y el ángulo de la cuerda 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 

Calcular la longitud y el ángulo de la línea de c a los centros de los dos círculos

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

calcular y trazar la posición de los dos centros d y e desde la perpendicular a ab y la longitud:

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: