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

72

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.

ingrese la descripción de la imagen aquí

Desarrollador
fuente

Respuestas:

45

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!

julien
fuente
11
+1 Muy buena respuesta! Me gustaría señalar que la rotación real de la nube es innecesaria. Primero, probablemente quisiste decir esto, solo los vértices del casco tienen que ser considerados. Segundo, en lugar de rotar, represente el lado actual como un par de vectores unitarios ortogonales. Tomar sus productos puntuales con las coordenadas del vértice del casco (que podría hacerse como una operación de matriz única) proporciona las coordenadas rotadas: no es necesaria la trigonometría, es rápida y perfectamente precisa.
whuber
2
Gracias por los enlaces. De hecho, la rotación solo para # de bordes hace que el método propuesto sea muy eficiente. Podría encontrar que el papel lo demuestra. Aunque marqué esto como la respuesta a la lealtad a la primera buena respuesta (no puedo elegir dos / más respuestas geniales :() Me gustaría recomendar fuertemente considerando la respuesta completa de Whuber a continuación. La eficiencia del método dado allí (¡evitando rotaciones!) Es increíble, y todo el procedimiento es solo unas pocas líneas de código. Para mí es fácilmente traducible a Python :)
Desarrollador el
¿Pueden actualizar el enlace de implementación de Java?
Myra
sí, ya está hecho!
julien
1
Tenga en cuenta que la extensión en 3D es un poco más complicada que eso. Cada cara del casco convexo 3D define una posible orientación de una cara del cuadro delimitador, pero no la orientación de las caras perpendiculares a ella. El problema de cómo rotar la caja en ese plano se convierte en el problema del rectángulo de límite mínimo 2D en el plano de esa cara. Para cada borde del casco convexo de la nube proyectada en un plano dado, puede dibujar un cuadro delimitador que le dará un volumen diferente en 3D.
Se
40

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

MBR

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.

whuber
fuente
+1 ¡Qué solución tan increíble! Tal idea surge solo después de tener largas experiencias. A partir de ahora tendré curiosidad por optimizar mis códigos existentes inspirándome con esta gran respuesta.
Desarrollador
Desearía poder votar esto dos veces. Estoy aprendiendo R y sus respuestas son una fuente continua de inspiración.
John Powell
1
@retrovius El rectángulo delimitador de un conjunto de puntos (rotados) está determinado por cuatro números: la coordenada x más pequeña, la coordenada x más grande, la coordenada y más pequeña y la coordenada y más grande. A eso se refieren los "extremos a lo largo de los bordes".
whuber
1
@retrovius El origen no juega ningún papel en estos cálculos, porque todo se basa en diferencias de coordenadas, excepto al final, donde el mejor rectángulo calculado en coordenadas rotadas simplemente se gira hacia atrás. Aunque es una idea inteligente usar un sistema de coordenadas en el que el origen esté cerca de los puntos (para minimizar la pérdida de precisión de coma flotante), el origen de lo contrario es irrelevante.
whuber
1
@Retrovius Puede interpretar esto en términos de una propiedad de rotaciones: a saber, la matriz de una rotación es ortogonal. Por lo tanto, un tipo de recurso sería un estudio de álgebra lineal (generalmente) o geometría analítica euclidiana (específicamente). Sin embargo, descubrí que la forma más fácil de lidiar con las rotaciones (y las traslaciones y reescalamientos) en el plano es ver los puntos como números complejos: las rotaciones se llevan a cabo simplemente multiplicando los valores por números de longitud unitaria.
Whuber
8

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.

ingrese la descripción de la imagen aquí

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
JesseBuesking
fuente
Hola JesseBuesking, ¿puedes generar rectángulos con esquinas de 90 grados? Su código funciona muy bien para obtener paralelogramos, pero se requieren esquinas de 90 grados en mi caso de uso específico. ¿Podría recomendar cómo se puede modificar su código para llegar a eso? ¡Gracias!
Nader Alexan el
@NaderAlexan Si está preguntando si puede manejar cuadrados, ¡entonces sí puede! Acabo de probarlo en un cuadrado de la unidad points = np.array([[0, 0], [0, 1], [1, 0], [1, 1]]), y la salida array([[1.00000000e+00, 6.12323400e-17], [0.00000000e+00, 0.00000000e+00], [6.12323400e-17, 1.00000000e+00], [1.00000000e+00, 1.00000000e+00]])es el cuadrado de la unidad (incluidos algunos errores de redondeo de coma flotante). Nota: un cuadrado es solo un rectángulo con lados iguales, así que supongo que si puede manejar un cuadrado, se generaliza a todos los rectángulos.
JesseBuesking
Gracias por su respuesta. Sí, funciona muy bien, pero estoy tratando de forzarlo a producir siempre un rectángulo (4 lados con ángulos de 90 grados para cada lado) sobre cualquier otro polígono de 4 lados, aunque en ciertos casos produce un rectángulo, no parece para ser una restricción constante, ¿sabe cómo modificar el código para agregar esta restricción? ¡Gracias!
Nader Alexan
Tal vez gis.stackexchange.com/a/22934/48041 pueda guiarlo hacia una solución, dado que su respuesta parece tener esta restricción. Una vez que encuentre una solución, debe contribuir ya que estoy seguro de que otros la encontrarán útil. ¡Buena suerte!
JesseBuesking
7

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.

ingrese la descripción de la imagen aquí


fuente
4

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.

David
fuente
¿Se ha eliminado tu respuesta?
Paul Richter
@PaulRichter aparentemente. La fuente estaba aquí github.com/dbworth/minimum-area-bounding-rectangle aunque
sehe
3

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):

ingrese la descripción de la imagen aquí

Usted golpear
fuente
¿Es posible extender esta implementación al espacio 3d (es decir, encontrar un cuadro de volumen mínimo que incluya todos los puntos dados en el espacio 3d)?
Sasha
0

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 .

trampa
fuente