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

63

Como se ve 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 reforman el problema (al eliminarlos 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] Encaja un rectángulo ajustado simple, es decir, paralelo a los ejes X e Y
    • puedes usar la función minmax para X e Y de los puntos dados (por ejemplo, los vértices de polígono)
  • Almacenar el área del rectángulo ajustado
  • Girar el polígono sobre el centroide, por ejemplo, 1 grado
  • Repita desde [S] hasta que se realice una rotación completa
  • Informe el ángulo del área mínima como el resultado

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

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

    
pregunta Developer 05.04.2012 - 10:03

7 respuestas

40

Sí, hay una solución analítica para este problema. El algoritmo que busca se conoce en 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 de la nube de puntos convexa casco . Así que solo necesita probar las orientaciones de los bordes del casco convexo. Deberías:

  • 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 there .

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!

    
respondido por el julien 05.04.2012 - 12:18
35

Para complementar la excelente solución de @ julien, aquí hay una implementación operativa en R , que podría servir como pseudocódigo para guiar cualquier implementación específica de GIS (o aplicarse directamente en R , de curso). 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 delimitador mínimo (con el primero repetido para cerrarlo). Tenga en cuenta la ausencia completa de cualquier cálculo trigonométrico.

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

Aquí hay un ejemplo de su uso:

# Create sample data
set.seed(23)
points <- 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, 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=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 los puntos n : puede calcular casi tan rápido como puede leer las coordenadas.

    
respondido por el whuber 05.04.2012 - 19:15
8

Simplemente implementé esto y publiqué mi respuesta en StackOverflow , pero pensé que dejaría mi versión aquí. para que otros 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 ello 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
    
respondido por el JesseBuesking 09.11.2015 - 22:59
7

Hay una herramienta en Whitebox GAT ( enlace ) llamada Caja de límites mínimos para resolver este problema exacto. También hay una herramienta de casco convexo mínimo allí también. Varias de las herramientas en la caja de herramientas Forma de parche, por ejemplo, La orientación y el alargamiento del parche se basan en encontrar el cuadro de límite mínimo.

    
respondido por el user21951 20.09.2013 - 15:41
4

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

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

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

    
respondido por el David 03.02.2013 - 19:38
3

Gracias @ la respuesta de whuber. Es una gran solución, pero lenta para la nube de puntos grandes. Encontré que la función convhulln en el paquete R geometry es mucho más rápida (138 s versus 0.03 s para 200000 puntos). He pegado mis códigos aquí para que cualquiera esté interesado en 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):

    
respondido por el Bangyou 21.12.2015 - 12:28
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 .

    
respondido por el pitfall 13.01.2018 - 06:25

Lea otras preguntas en las etiquetas