2012-09-04 14 views
8

Mam parę punktów i chciałbym znaleźć kółka o znanym r, które są określone przez te dwa punkty. Będę używał tego w symulacji i możliwej przestrzeni dla x i y mają granice (powiedzmy pudełko z -200, 200).wyznaczyć środek okręgu na podstawie dwóch punktów (znany promień) z rozwiązać/optyma

It is known że kwadrat promienia jest

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

Chciałbym teraz, aby rozwiązać ten nieliniowy układ równań, aby uzyskać dwa potencjalne centra koła. Próbowałem użyć pakietu BB. Oto moja słaba próba, która daje tylko jeden punkt. To, co chciałbym uzyskać, to oba możliwe punkty. Wszelkie wskazówki w dobrym kierunku spotkają się z bezpłatnym piwem przy pierwszej możliwej okazji.

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 
       ) 

Graficznie można to zobaczyć za pomocą następującego kodu, ale potrzebne będą dodatkowe pakiety.

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

DODATEK

Dziękuję wszystkim za cenne uwagi i kodu. Dostaję czas na odpowiedzi przez plakaty, które uzupełniły swoje odpowiedzi kodem.

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 
+1

wykonaj punkty leżą na obwodzie? – James

+0

Możliwe jest rozwiązanie układów równania rękami i używanie wzorów – MBo

+0

@ James, tak, punkty leżą gdzieś na obwodzie. Zaktualizowałem swoją odpowiedź, która pokazuje wynik. –

Odpowiedz

4

Jest to podstawowy geometryczny sposób rozwiązywania problemu, o którym wspominają wszyscy inni. Używam polyroot, aby uzyskać korzenie wynikowego równania kwadratowego, ale możesz łatwo użyć równania kwadratowego bezpośrednio.

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

Oto kości odpowiedzi, jeśli będę miał czas później, wypalę je. To powinno być łatwe do naśladowania, jeśli narysujesz razem ze słowami, przepraszam, nie mam odpowiedniego oprogramowania na tym komputerze, aby narysować obraz dla ciebie.

Pomijaj zdegenerowane przypadki, w których punkty są identyczne (nieskończone rozwiązania) lub zbyt daleko, aby leżeć na tym samym okręgu o wybranym promieniu (brak rozwiązań).

Oznacz punkty X i Y oraz nieznane punkty środkowe 2 okręgów c1 i c2. c1 i c2 leżą na dwusiecznej prostopadłej z XY; nazwij tę linię c1c2, na tym etapie nie ma znaczenia, że ​​nie znamy wszystkich szczegółów lokalizacji c1 i c2.

Wyobraźmy sobie równanie linii c1c2. Przechodzi przez punkt pół-drogi XY (nazywa ten punkt Z) i ma nachylenie równe ujemnemu odwrotności XY. Teraz masz równanie c1c2 (lub gdybyś miał jakiekolwiek ciało na tych kościach).

Teraz skonstruuj trójkąt od jednego punktu do przecięcia linii i jej prostopadłej dwusiecznej oraz środkowego punktu okręgu (na przykład XZc1). Nadal nie wiesz dokładnie, gdzie jest c1, ale nigdy nie zatrzymał nikogo, kto szkicuje geometrię. Masz trójkąt prostokątny o dwóch znanych długościach boków (XZ i Xc1), więc łatwo jest znaleźć Zc1. Powtórz ten proces dla drugiego trójkąta i środka okręgu.

Oczywiście podejście to różni się od podejścia początkowego OP i nie może się odwoływać.

1

Kilka ostrzeżeń, których należy się pozbyć, ale to powinno zacząć. Możliwe, że wystąpił problem z wydajnością, więc rozwiązanie go w całości za pomocą podstawowej geometrii może być lepszym rozwiązaniem.

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

Mam nadzieję, że znasz podstawową geometrię, ponieważ niestety nie potrafię jej narysować.

Dwusieczna prostopadła to linia, w której znajduje się każdy środkowy punkt koła, który przecina zarówno A, jak i B.

Teraz masz środek AB i r, więc możesz narysować trójkąt prostokątny z punktem A, środkiem AB i nieznanym środkowym punktem koła.

Teraz użyj twierdzenia Pitagorasa, aby uzyskać odległość od środkowego punktu AB do środkowego punktu koła, i oblicz położenie koła nie powinno być tutaj trudne, używając podstawowych kombinacji sin/cos.

3

Nie wymaga się rozwiązywania równań numerycznych. Tylko formuły:

  1. Wiesz, że ponieważ oba punkty A i B leżą na okręgu, odległość od każdego do danego środka jest promieniem r.
  2. Tworzy trójkąt równoramienny z cięciwą dwóch znanych punktów w podstawie i trzecim punktem w środku okręgu.
  3. Podziel pół trójkąta między A i B, podając trójkąt prostokątny.
  4. http://mathworld.wolfram.com/IsoscelesTriangle.html podaje wysokość pod względem długości podstawy i promienia.
  5. Podążaj normalnie do akordu AB (See this SO Answer), aby uzyskać odległość o wysokości obliczonej w każdym kierunku od punktu.
9

Poniższy kod wyświetli punkty w centrach dwóch żądanych okręgów. Nie ma teraz czasu, aby to skomentować lub przekonwertować wyniki na obiekty Spatial*, ale powinno to dać dobry początek.

Po pierwsze, oto schemat ASCII-art, aby wprowadzić nazwy punktów. k i K są znane punkty, B jest punktem na linii poziomej poprowadzonej przez k i C1 i C2 są centra kręgach szukasz:

 C2 





          K 


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






             C1 

teraz kodu:

# 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

4

Po @ rozwiązania PhilH, tuż za pomocą trygonometrii w R:

radius=40 

Draw oryginalnych punktów na promieniu

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

Znajdź punkt środkowy c z ab (który jest również midpoint of de dwóch ośrodków okręgu)

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) 

Znajdź długość i kąt cięciwy 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 

Obliczenie długości i kąta linii z c do centrów dwóch kół

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

obliczenia i wykreślenia pozycję dwa centra: d i e od prostopadłej do ab i długości:

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

Pokazy:

enter image description here

Powiązane problemy