Как найти наименьшую описанную окружность неправильного многоугольника в проекте R?

Question

Мне было интересно, как найти наименьшую описанную окружность неправильного многоугольника. Я работал с пространственными многоугольниками в R.

Я хочу воспроизвести некоторые метрики fragstats в векторном режиме, потому что у меня были трудности с пакетом «landscapemetrics» для огромного количества данных. В частности, c я хотел бы реализовать круг (http://www.umass.edu/landeco/research/fragstats/documents/Metrics/Shape%20Metrics/Metrics/P11%20-%20CIRCLE.htm). Пока мне не удалось найти формулу или сценарий для наименьшего описанного круга.

Все ваши комментарии приветствуются.

Чем вы

Answer 1

Как я уже упоминал в комментарии, я не знаю существующего R-кода для этого, но поиск методом грубой силы должен быть достаточно быстрым, если у вас не слишком много точек, которые должны быть в круге. Я только что написал это. Функция center() основана на коде из Википедии для рисования круга вокруг треугольника; circumcircle() - это нужная функция, найденная методом перебора всех кругов, проходящих через 2 или 3 точки в наборе. На моем ноутбуке обработка 100 точек занимает около 4 секунд. Если у вас несколько больших наборов, вы, вероятно, сможете получить приемлемые результаты, переводя на C ++, но это скорость роста n^4, поэтому вам понадобится лучшее решение для действительно большого набора.

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)

Это показывает результат