Search code examples
rtriangulationleast-squares

Least Squares Triangulation in R


Which are the coordinates of an unknown point if they are given observations of the distances of 3 points with known coordinates?

eg:

x = c(30.0,10.0,50.0)
y = c(150.0,120.0,50.0)
distance = c("125.0 ± 0.5","133.5 ± 0.2","98.6 ± 0.2")
df = data.frame(x, y, distance) 

Is it possible to use R to calculate the coordinates of the unknown point using least squares with (a) indirect observations method and (b) condition equations?

How data in the form of "125.0 ± 0.5" can be imported in R?


Solution

  • This could be a start. I think there are ways to linearalize the equations, thus simplifying the procedure. But, this just takes the obvious approach of finding the intersection of the three circles, and uses nlm to minimize the squares. I don't deal with the errors at all here.

    ## Your data
    x = c(30.0,10.0,50.0)
    y = c(150.0,120.0,50.0)
    ## distance = c("125.0 ± 0.5","133.5 ± 0.2","98.6 ± 0.2")
    dists <- c(125, 133.5, 98.6)  # simplified
    
    ## Minimize this function:
    ## x: guesstimates
    ## centers: fixed points (x, y)
    ## b: distances
    f <- function(x, centers, b) {
        sqrt(sum((sqrt(colSums((x - centers)^2)) - b)^2))
    }
    
    ## Get estimate: initial guess of c(100, 100)
    centers <- matrix(c(x, y), nrow=2, byrow=TRUE)
    res <- nlm(f, c(100, 100), centers=centers, b=dists)
    
    ## Lets visualize to see if it worked
    circle <- function(x, y, d, col='black') {
        theta <- seq(0, 2*pi, length.out=100)
        data.frame(x=d*cos(theta)+x, y=d*sin(theta)+y, col, stringsAsFactors=FALSE)
    }
    cols <- colorRampPalette(c('blue', 'red'))(3)
    circs <- Map(circle, x, y, dists, cols)
    ps <- do.call(rbind, circs)
    plot(ps[1:2], type='n')
    grid()
    abline(h=0, v=0)
    points(x, y, col=cols, pch=16, cex=2)
    for (i in circs) points(i[1:2], col=i$col, type='l')
    
    ## Add the estimated point
    points(x=res$estimate[1], y=res$estimate[2], 
           col='firebrick', pch=8, cex=2)
    

    enter image description here