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

Solución:

Sí, existe una solución analítica para este problema. El algoritmo que está buscando se conoce en la generalización de polígonos como “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. Entonces, solo necesita probar las orientaciones de los bordes convexos del casco. Debería:

• Calcule el casco convexo de la nube.
• Para cada borde del casco convexo:
  • calcular la orientación del borde (con arctan),
  • rotar el casco convexo usando 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 rotado,
  • Almacene la orientación correspondiente al área mínima encontrada,
• Devuelve el rectángulo correspondiente al área mínima encontrada.

Allí se encuentra disponible un ejemplo de implementación en Java.

En 3D, se aplica lo mismo, 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 de trabajo en R, que podría servir como pseudocódigo para guiar cualquier implementación específica de SIG (o aplicarse directamente en R, por supuesto). La entrada es una matriz de coordenadas de puntos. Salida (el valor de mbr) es una matriz de los vértices del rectángulo delimitador mínimo (con el primero repetido para cerrarlo). Tenga en cuenta la ausencia total de cálculos trigonométricos.

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

A continuación, se muestra un ejemplo de su uso:

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

# 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(p, pch=19)                                     # The points

El tiempo está limitado por la velocidad del algoritmo del 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 norte 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 dejarí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 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

---

---

Utiliza Nuestro Buscador