determinare il centro del cerchio in base a due punti (raggio noto) con solve/optim

Question

Ho una coppia di punti e mi piacerebbe trovare un cerchi di r noti che sono determinati da questi due punti. Lo userò in una simulazione e uno spazio possibile per x e con limiti (diciamo una casella di -200, 200).

It is known quel quadrato del raggio è

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

Vorrei ora risolvere questo sistema non lineare di equazioni per ottenere due potenziali centri di cerchio. Ho provato a utilizzare il pacchetto BB. Ecco il mio debole tentativo che dà solo un punto. Quello che vorrei ottenere è entrambi i punti possibili. Ogni suggerimento nella giusta direzione sarà accolto con birra in omaggio alla prima occasione possibile.

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 
       ) 

Graficamente si può vedere questo con il seguente codice, ma avrete bisogno di alcuni pacchetti aggiuntivi.

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) 

ADDENDUM

Grazie a tutti per i vostri commenti e il codice di valore. Fornisco i tempi per le risposte dei poster che si sono complimentati con le loro risposte con il codice.

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

Questo è il modo geometrico di base per risolvere il problema che tutti gli altri stanno menzionando. Io uso polyroot per ottenere le radici dell'equazione quadratica risultante, ma si potrebbe facilmente utilizzare direttamente l'equazione quadratica.

# 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

Ecco le ossa di una risposta, se ho tempo dopo io li rimpolpare. Questo dovrebbe essere abbastanza facile da seguire se si disegna con le parole, mi dispiace non ho il software giusto su questo computer per disegnare l'immagine per voi.

Lasciare da parte casi degenerati in cui i punti sono identici (soluzioni infinite) o troppo distanti per trovarsi nello stesso cerchio del raggio scelto (nessuna soluzione).

Etichettare i punti X e Y e i punti centrali sconosciuti dei 2 cerchi c1 e c2. c1 e c2 si trovano sulla bisettrice perpendicolare di XY; chiama questa riga c1c2, in questa fase è irrilevante che non conosciamo tutti i dettagli delle posizioni di c1 e c2.

Quindi, capire l'equazione della linea c1c2. Passa attraverso il punto intermedio di XY (chiama questo punto Z) e ha una pendenza uguale al reciproco negativo di XY. Ora hai l'equazione di c1c2 (o lo faresti se ci fosse carne su queste ossa).

Ora costruisci il triangolo da un punto all'intersezione della linea e della sua bisettrice perpendicolare e il punto centrale di un cerchio (ad esempio XZc1). Non si sa ancora esattamente dove sia c1 ma non si è mai fermato nessuno a disegnare la geometria. Hai un triangolo rettangolo con due lunghezze laterali conosciute (XZ e Xc1), quindi è facile trovare Zc1. Ripeti il ​​processo per l'altro triangolo e centro del cerchio.

Ovviamente, questo approccio è molto diverso dall'approccio iniziale di OP e potrebbe non fare appello.

Answer 3

Alcuni avvisi per sbarazzarsi di, ma questo dovrebbe iniziare. Potrebbe esserci un problema di prestazioni, quindi risolverlo completamente con la geometria di base potrebbe essere un approccio migliore.

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

Spero che tu sappia un po 'di geometria di base, perché non riesco a disegnarlo purtroppo.

La bisettrice perpendicolare è la linea in cui ogni punto medio di un cerchio che attraversa sia A sia B pone.

Ora avete il centro di AB e r, quindi potete disegnare un triangolo rettangolo con il punto A, il centro di AB e il punto medio sconosciuto del cerchio.

Ora utilizzare il teorema di Pitagora per ottenere la distanza dal punto medio di AB al punto centrale del cerchio e calcolare la posizione del cerchio non dovrebbe essere difficile da qui, utilizzando le combinazioni sin/cos di base.

Answer 5

Nessuna risoluzione di equazione numerica richiesta. Solo formule:

1. Sapete che poiché entrambi i punti A e B si trovano sul cerchio, la distanza da ciascuno a un centro dato è il raggio r.
2. Forma un triangolo isoscele con l'accordo dei due punti noti alla base e il terzo punto al centro del cerchio.
3. Tagliare il triangolo a metà strada tra A e B, dando un triangolo ad angolo retto.
4. http://mathworld.wolfram.com/IsoscelesTriangle.html fornisce l'altezza in termini di lunghezza della base e raggio.
5. Seguire la normale sull'accordo AB (See this SO Answer) per una distanza dell'altezza appena calcolata in ciascuna direzione dal punto.

Answer 6

Il seguente codice consente di ottenere i punti al centro dei due cerchi desiderati. Non c'è tempo ora per commentare questo o convertire i risultati in oggetti Spatial*, ma questo dovrebbe darti un buon inizio.

Per prima cosa, ecco uno schema ASCII per introdurre i nomi dei punti. k e K sono i punti noti, B è un punto sul piano orizzontale che passa per k e C1 e C2 sono i centri dei cerchi che stai cercando:

 C2 





          K 


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






             C1 

Ora il codice:

# 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

seguito @ soluzione PhilH, semplicemente utilizzando la trigonometria in R:

radius=40 

Disegnare i punti originali sul raggio

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

Trova il punto medio c di ab (che è anche il punto centrale di de centri dei due cerchi)

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) 

Trovare la lunghezza e l'angolo della corda 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 

Calcolare la lunghezza e l'angolo della linea da c ai centri dei due cerchi

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

calcolare e tracciare la posizione dei i due centri d e e dalla perpendicolare a ab e la lunghezza:

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

Spettacoli: