¿Encontrar el rectángulo de área mínima para los puntos dados?

Como puede ver en la figura, la pregunta es:

¿Cómo encontrar el rectángulo de área mínima (MAR) ajustado en los puntos dados?

y una pregunta de apoyo es:

¿Hay alguna solución analítica para el problema?

(Un desarrollo de la pregunta será ajustar un cuadro (3D) a un grupo de puntos en una nube de puntos 3D).

Como primera etapa, propongo encontrar el casco convexo para los puntos que reforma el problema (eliminando esos puntos que no están involucrados en la solución) para: ajustar un MAR a un polígono. El método requerido proporcionará X ( centro del rectángulo ), D ( dos dimensiones ) y A ( ángulo ).

Mi propuesta de solución:

• Encuentre el centroide del polígono (vea ¿ Encontrar el centro de geometría del objeto? )
• [S] Ajustar un rectángulo ajustado simple, es decir, paralelo a los ejes X e Y
  • puede usar la minmaxfunción para X e Y de los puntos dados (por ejemplo, vértices de polígonos)
• Almacene el área del rectángulo ajustado
• Gire el polígono sobre el centroide por ejemplo, 1 grado
• Repita desde [S] hasta completar una rotación completa
• Informe el ángulo del área mínima como resultado

Me parece prometedor, sin embargo, existen los siguientes problemas:

• elegir una buena resolución para el cambio de ángulo podría ser un desafío,
• el costo de cálculo es alto,
• La solución no es analítica sino experimental.

Sí, hay una solución analítica para este problema. El algoritmo que está buscando se conoce en la generalización de polígonos como "el rectángulo circundante más pequeño".

El algoritmo que describe está bien, pero para resolver los problemas que ha enumerado, puede utilizar el hecho de que la orientación del MAR es la misma que la de uno de los bordes del casco convexo de la nube de puntos . Por lo tanto, solo necesita probar las orientaciones de los bordes convexos del casco. Debieras:

• Calcule el casco convexo de la nube.
• Para cada borde del casco convexo:
  • calcular la orientación del borde (con arctan),
  • gire el casco convexo utilizando esta orientación para calcular fácilmente el área del rectángulo delimitador con un mínimo / máximo de x / y del casco convexo girado,
  • Almacene la orientación correspondiente al área mínima encontrada,
• Devuelve el rectángulo correspondiente al área mínima encontrada.

Un ejemplo de implementación en java está disponible allí .

En 3D, lo mismo se aplica, excepto:

• El casco convexo será un volumen,
• Las orientaciones probadas serán las orientaciones (en 3D) de las caras convexas del casco.

¡Buena suerte!

Para complementar la gran solución de @ julien, aquí hay una implementación funcional Rque podría servir como pseudocódigo para guiar cualquier implementación específica de SIG (o aplicarse directamente R, por supuesto). La entrada es una matriz de coordenadas de puntos. La salida (el valor de mbr) es una matriz de los vértices del rectángulo límite mínimo (con el primero repetido para cerrarlo). Tenga en cuenta la ausencia total de cualquier cálculo trigonométrico.

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,])
}

Aquí hay un ejemplo de su uso:

# 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

El tiempo está limitado por la velocidad del algoritmo de casco convexo, porque el número de vértices en el casco es casi siempre mucho menor que el total. La mayoría de los algoritmos de casco convexo son asintóticamente O (n * log (n)) para n puntos: puede calcular casi tan rápido como puede leer las coordenadas.

Acabo de implementar esto yo mismo y publiqué mi respuesta en StackOverflow , pero pensé que soltaría mi versión aquí para que otros la vean:

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

Aquí hay cuatro ejemplos diferentes de esto en acción. Para cada ejemplo, generé 4 puntos aleatorios y encontré el cuadro delimitador.

También es relativamente rápido para estas muestras en 4 puntos:

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

Hay una herramienta en Whitebox GAT ( http://www.uoguelph.ca/~hydrogeo/Whitebox/ ) llamada Minimum Bounding Box para resolver este problema exacto. También hay una herramienta mínima de casco convexo allí también. Varias de las herramientas en la caja de herramientas Forma de parche, por ejemplo, orientación y alargamiento de parche, se basan en encontrar el cuadro de límite mínimo.

Encontré este hilo mientras buscaba una solución de Python para un rectángulo delimitador de área mínima.

Aquí está mi implementación , para la cual los resultados fueron verificados con Matlab.

El código de prueba se incluye para polígonos simples, y lo estoy usando para encontrar el cuadro de límite mínimo 2D y las direcciones de los ejes para un PointCloud 3D.

Gracias la respuesta de @ whuber. Es una gran solución, pero lenta para grandes nubes de puntos. Encontré que la convhullnfunción en el paquete R geometryes mucho más rápida (138 s frente a 0.03 s para 200000 puntos). Pegué mis códigos aquí para que cualquiera sea interesante para una solución más rápida.

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)   

Dos métodos obtienen la misma respuesta (ejemplo para 2000 puntos):

Simplemente recomiendo la función incorporada de OpenCV minAreaRect, que encuentra un rectángulo girado del área mínima que encierra el conjunto de puntos 2D de entrada. Para ver cómo usar esta función, puede consultar este tutorial .