Нахождение минимальной площади прямоугольника для заданных точек?

Как вы видите на рисунке, вопрос заключается в следующем:

Как найти прямоугольник минимальной площади (MAR), расположенный в заданных точках?

и подтверждающий вопрос:

Есть ли аналитическое решение проблемы?

(Развитие вопроса будет состоять в том, чтобы поместить прямоугольник (3D) в кластер точек в трехмерном облаке точек.)

В качестве первого этапа я предлагаю найти выпуклую оболочку для точек, которые решают проблему (удаляя эти точки, которые не участвуют в решении) для: подгонки MAR к многоугольнику. Требуемый метод предоставит X ( центр прямоугольника ), D ( два измерения ) и A ( угол ).

Мое предложение для решения:

• Найти центр тяжести многоугольника (см. Поиск центра геометрии объекта? )
• [S] Установите простой подогнанный прямоугольник, т. Е. Параллельно осям X и Y
  • Вы можете использовать minmaxфункцию для X и Y заданных точек (например, вершин многоугольника)
• Сохраните область подогнанного прямоугольника
• Поворот многоугольника вокруг центроида, например, на 1 градус
• Повторите от [S] до полного вращения
• Сообщите угол минимальной площади как результат

Это выглядит многообещающе, однако существуют следующие проблемы:

• выбор хорошего разрешения для изменения угла может быть сложной задачей,
• высокая стоимость вычислений,
• решение не аналитическое, а экспериментальное.

Да, есть аналитическое решение этой проблемы. Алгоритм, который вы ищете, в обобщении полигонов известен как «наименьший окружающий прямоугольник».

Алгоритм, который вы описываете, хорош, но для решения перечисленных проблем вы можете использовать тот факт, что ориентация MAR такая же, как и у одного из краев выпуклой оболочки облака точек . Так что вам просто нужно проверить ориентацию выпуклых краев корпуса. Вам следует:

• Вычислить выпуклую оболочку облака.
• Для каждого края выпуклой оболочки:
  • вычислить ориентацию края (с помощью arctan),
  • поверните выпуклый корпус, используя эту ориентацию, чтобы легко вычислить площадь ограничивающего прямоугольника с минимумом / максимумом x / y повернутого выпуклого корпуса,
  • Сохраните ориентацию, соответствующую найденной минимальной площади,
• Вернуть прямоугольник, соответствующий минимальной найденной площади.

Пример реализации в Java доступен там .

В 3D применяется то же самое, за исключением:

• Выпуклая оболочка будет объемной,
• Проверенными ориентациями будут ориентации (в 3D) выпуклых граней корпуса.

Удачи!

В дополнение к отличному решению @ julien, здесь есть рабочая реализация R, которая может служить псевдокодом для руководства любой конкретной реализацией ГИС (или R, конечно, применяться напрямую ). Ввод представляет собой массив координат точки. Выход (значение mbr) представляет собой массив вершин минимального ограничивающего прямоугольника (с повторением первой, чтобы закрыть его). Обратите внимание на полное отсутствие каких-либо тригонометрических расчетов.

MBR <- function(p) {
  # Analyze the convex hull edges     
  a <- chull(p)                                   # Indexes of extremal points
  a <- c(a, a[1])                                 # Close the loop
  e <- p[a[-1],] - p[a[-length(a)], ]             # Edge directions
  norms <- sqrt(rowSums(e^2))                     # Edge lengths
  v <- e / norms                                  # Unit edge directions
  w <- cbind(-v[,2], v[,1])                       # Normal directions to the edges

  # Find the MBR
  vertices <- p[a, ]                              # Convex hull vertices
  x <- apply(vertices %*% t(v), 2, range)         # Extremes along edges
  y <- apply(vertices %*% t(w), 2, range)         # Extremes normal to edges
  areas <- (y[1,]-y[2,])*(x[1,]-x[2,])            # Areas
  k <- which.min(areas)                           # Index of the best edge (smallest area)

  # Form a rectangle from the extremes of the best edge
  cbind(x[c(1,2,2,1,1),k], y[c(1,1,2,2,1),k]) %*% rbind(v[k,], w[k,])
}

Вот пример его использования:

# Create sample data
set.seed(23)
p <- matrix(rnorm(20*2), ncol=2)                 # Random (normally distributed) points
mbr <- MBR(points)

# Plot the hull, the MBR, and the points
limits <- apply(mbr, 2, range) # Plotting limits
plot(p[(function(x) c(x, x[1]))(chull(p)), ], 
     type="l", asp=1, bty="n", xaxt="n", yaxt="n",
     col="Gray", pch=20, 
     xlab="", ylab="",
     xlim=limits[,1], ylim=limits[,2])                # The hull
lines(mbr, col="Blue", lwd=3)                         # The MBR
points(points, pch=19)                                # The points

Время ограничено скоростью алгоритма выпуклой оболочки, потому что количество вершин в корпусе почти всегда намного меньше, чем общее. Большинство алгоритмов выпуклой оболочки асимптотически O (n * log (n)) для n точек: вы можете вычислить почти так же быстро, как вы можете прочитать координаты.

Я сам реализовал это и опубликовал свой ответ в StackOverflow , но я решил оставить свою версию здесь для просмотра другими:

import numpy as np
from scipy.spatial import ConvexHull

def minimum_bounding_rectangle(points):
    """
    Find the smallest bounding rectangle for a set of points.
    Returns a set of points representing the corners of the bounding box.

    :param points: an nx2 matrix of coordinates
    :rval: an nx2 matrix of coordinates
    """
    from scipy.ndimage.interpolation import rotate
    pi2 = np.pi/2.

    # get the convex hull for the points
    hull_points = points[ConvexHull(points).vertices]

    # calculate edge angles
    edges = np.zeros((len(hull_points)-1, 2))
    edges = hull_points[1:] - hull_points[:-1]

    angles = np.zeros((len(edges)))
    angles = np.arctan2(edges[:, 1], edges[:, 0])

    angles = np.abs(np.mod(angles, pi2))
    angles = np.unique(angles)

    # find rotation matrices
    # XXX both work
    rotations = np.vstack([
        np.cos(angles),
        np.cos(angles-pi2),
        np.cos(angles+pi2),
        np.cos(angles)]).T
#     rotations = np.vstack([
#         np.cos(angles),
#         -np.sin(angles),
#         np.sin(angles),
#         np.cos(angles)]).T
    rotations = rotations.reshape((-1, 2, 2))

    # apply rotations to the hull
    rot_points = np.dot(rotations, hull_points.T)

    # find the bounding points
    min_x = np.nanmin(rot_points[:, 0], axis=1)
    max_x = np.nanmax(rot_points[:, 0], axis=1)
    min_y = np.nanmin(rot_points[:, 1], axis=1)
    max_y = np.nanmax(rot_points[:, 1], axis=1)

    # find the box with the best area
    areas = (max_x - min_x) * (max_y - min_y)
    best_idx = np.argmin(areas)

    # return the best box
    x1 = max_x[best_idx]
    x2 = min_x[best_idx]
    y1 = max_y[best_idx]
    y2 = min_y[best_idx]
    r = rotations[best_idx]

    rval = np.zeros((4, 2))
    rval[0] = np.dot([x1, y2], r)
    rval[1] = np.dot([x2, y2], r)
    rval[2] = np.dot([x2, y1], r)
    rval[3] = np.dot([x1, y1], r)

    return rval

Вот четыре различных примера этого в действии. Для каждого примера я сгенерировал 4 случайных точки и нашел ограничивающий прямоугольник.

Это сравнительно быстро для этих образцов на 4 балла:

>>> %timeit minimum_bounding_rectangle(a)
1000 loops, best of 3: 245 µs per loop

В Whitebox GAT есть инструмент ( http://www.uoguelph.ca/~hydrogeo/Whitebox/ ), который называется Minimum Bounding Box для решения именно этой проблемы. Там также есть инструмент с минимальным выпуклым корпусом. Некоторые из инструментов в наборе инструментов Patch Shape, например, ориентация и удлинение патча, основаны на поиске минимальной ограничительной рамки.

Я наткнулся на эту тему, когда искал решение Python для ограничивающего прямоугольника с минимальной площадью.

Вот моя реализация , результаты которой были проверены с помощью Matlab.

Тестовый код включен для простых полигонов, и я использую его, чтобы найти двухмерную минимальную ограничивающую рамку и направления осей для 3D PointCloud.

Спасибо, ответ @ whuber. Это отличное решение, но медленное для большого облака точек. Я обнаружил, что convhullnфункция в пакете R geometryработает намного быстрее (138 с против 0,03 с для 200000 баллов). Я вставил здесь свои коды для всех, кто заинтересован в более быстром решении.

library(alphahull)                                  # Exposes ashape()
MBR <- function(points) {
    # Analyze the convex hull edges                       
    a <- ashape(points, alpha=1000)                 # One way to get a convex hull...
    e <- a$edges[, 5:6] - a$edges[, 3:4]            # Edge directions
    norms <- apply(e, 1, function(x) sqrt(x %*% x)) # Edge lengths
    v <- diag(1/norms) %*% e                        # Unit edge directions
    w <- cbind(-v[,2], v[,1])                       # Normal directions to the edges

    # Find the MBR
    vertices <- (points) [a$alpha.extremes, 1:2]    # Convex hull vertices
    minmax <- function(x) c(min(x), max(x))         # Computes min and max
    x <- apply(vertices %*% t(v), 2, minmax)        # Extremes along edges
    y <- apply(vertices %*% t(w), 2, minmax)        # Extremes normal to edges
    areas <- (y[1,]-y[2,])*(x[1,]-x[2,])            # Areas
    k <- which.min(areas)                           # Index of the best edge (smallest area)

    # Form a rectangle from the extremes of the best edge
    cbind(x[c(1,2,2,1,1),k], y[c(1,1,2,2,1),k]) %*% rbind(v[k,], w[k,])
}

MBR2 <- function(points) {
    tryCatch({
        a2 <- geometry::convhulln(points, options = 'FA')

        e <- points[a2$hull[,2],] - points[a2$hull[,1],]            # Edge directions
        norms <- apply(e, 1, function(x) sqrt(x %*% x)) # Edge lengths

        v <- diag(1/norms) %*% as.matrix(e)                        # Unit edge directions


        w <- cbind(-v[,2], v[,1])                       # Normal directions to the edges

        # Find the MBR
        vertices <- as.matrix((points) [a2$hull, 1:2])    # Convex hull vertices
        minmax <- function(x) c(min(x), max(x))         # Computes min and max
        x <- apply(vertices %*% t(v), 2, minmax)        # Extremes along edges
        y <- apply(vertices %*% t(w), 2, minmax)        # Extremes normal to edges
        areas <- (y[1,]-y[2,])*(x[1,]-x[2,])            # Areas
        k <- which.min(areas)                           # Index of the best edge (smallest area)

        # Form a rectangle from the extremes of the best edge
        as.data.frame(cbind(x[c(1,2,2,1,1),k], y[c(1,1,2,2,1),k]) %*% rbind(v[k,], w[k,]))
    }, error = function(e) {
        assign('points', points, .GlobalEnv)
        stop(e)  
    })
}


# Create sample data
#set.seed(23)
points <- matrix(rnorm(200000*2), ncol=2)                 # Random (normally distributed) points
system.time(mbr <- MBR(points))
system.time(mmbr2 <- MBR2(points))


# Plot the hull, the MBR, and the points
limits <- apply(mbr, 2, function(x) c(min(x),max(x))) # Plotting limits
plot(ashape(points, alpha=1000), col="Gray", pch=20, 
     xlim=limits[,1], ylim=limits[,2])                # The hull
lines(mbr, col="Blue", lwd=10)                         # The MBR
lines(mbr2, col="red", lwd=3)                         # The MBR2
points(points, pch=19)   

Два метода получают один и тот же ответ (пример для 2000 баллов):

Я просто рекомендую встроенную функцию OpenCV minAreaRect, которая находит повернутый прямоугольник минимальной области, содержащей входной набор 2D-точек. Чтобы увидеть, как использовать эту функцию, можно обратиться к этому руководству .