rgisspatial-index

How to find the smallest circumcircle of an irregular polygon on R project?


I was wondering about how to find the smallest circumcircle of an irregular polygon. I've worked with spatial polygons in R.

I want to reproduce some of the fragstats metrics in a vector mode because I had hard times with the package 'landscapemetrics' for a huge amount of data. In specific I would like to implement the circle (http://www.umass.edu/landeco/research/fragstats/documents/Metrics/Shape%20Metrics/Metrics/P11%20-%20CIRCLE.htm). So far, I could not find the formula or script for the smallest circumcircle.

All your comments are more than welcome.

Than you


Solution

  • As I mentioned in a comment, I don't know of existing R code for this, but a brute force search should be fast enough if you don't have too many points that need to be in the circle. I just wrote this one. The center() function is based on code from Wikipedia for drawing a circle around a triangle; circumcircle() is the function you want, found by brute force search through all circles that pass through 2 or 3 points in the set. On my laptop it takes about 4 seconds to handle 100 points. If you have somewhat bigger sets, you can probably get tolerable results by translating to C++, but it's an n^4 growth rate, so you'll need a better solution for a really large set.

    center <- function(D) {
      if (NROW(D) == 0)
        matrix(numeric(), ncol = 2)
      else if (NROW(D) == 1)
        D
      else if (NROW(D) == 2) {
        (D[1,] + D[2,])/2
      } else if (NROW(D) == 3) {
        B <- D[2,] - D[1,]
        C <- D[3,] - D[1,]
        Dprime <- 2*(B[1]*C[2] - B[2]*C[1])
        if (Dprime == 0) {
          drop <- which.max(c(sum((B-C)^2), sum(C^2), sum(B^2)))
          center(D[-drop,])
        } else 
          c((C[2]*sum(B^2) - B[2]*sum(C^2))/Dprime,
            (B[1]*sum(C^2) - C[1]*sum(B^2))/Dprime) + D[1,]
      } else 
        center(circumcircle(D))
    }
    
    radius <- function(D, U = center(D))
      sqrt(sum((D[1,] - U)^2))
    
    circumcircle <- function(P) {
      n <- NROW(P)
      if (n < 3) 
        return(P)
      P <- P[sample(n),]
      bestset <- NULL
      bestrsq <- Inf
      # Brute force search
      for (i in 1:(n-1)) {
        for (j in (i+1):n) {
          D <- P[c(i,j),]
          U <- center(D)
          rsq <- sum((D[1,] - U)^2)
          if (rsq >= bestrsq)
            next
          failed <- FALSE
          for (k in (1:n)[-j][-i]) {
            Pk <- P[k,,drop = FALSE]
            if (sum((Pk - U)^2) > rsq) {
              failed <- TRUE
              break
            }
          }
          if (!failed) {
            bestset <- c(i,j)
            bestrsq <- rsq
          }
        }
      }
      # Look for the best 3 point set
      for (i in 1:(n-2)) {
        for (j in (i+1):(n-1)) {
          for (l in (j+1):n) {
            D <- P[c(i,j,l),]
            U <- center(D)
            rsq <- sum((D[1,] - U)^2)
            if (rsq >= bestrsq)
              next
            failed <- FALSE
            for (k in (1:n)[-l][-j][-i]) {
              Pk <- P[k,,drop = FALSE]
              if (sum((Pk - U)^2) > rsq) {
                failed <- TRUE
                break
              }
            }
            if (!failed) {
              bestset <- c(i,j,l)
              bestrsq <- rsq
            }
          }
        }  
      }
      P[bestset,]
    }
    
    showP <- function(P, ...) {
      plot(P, asp = 1, type = "n", ...)
      text(P, labels = seq_len(nrow(P)))
    }
    
    showD <- function(D) {
      U <- center(D)
      r <- radius(D, U)
      theta <- seq(0, 2*pi, len = 100)
      lines(U[1] + r*cos(theta), U[2] + r*sin(theta))
    }
    
    n <- 100
    P <- cbind(rnorm(n), rnorm(n))
    D <- circumcircle(P)
    showP(P)
    showD(D)
    

    This shows the output

    screenshot