Encontrando um retângulo de área mínima para determinados pontos?

Como você vê na figura, a pergunta é:

Como encontrar o retângulo de área mínima (MAR) ajustado nos pontos indicados?

e uma pergunta de apoio é:

Existe alguma solução analítica para o problema?

(O desenvolvimento da pergunta será ajustar uma caixa (3D) a um cluster de pontos em uma nuvem de pontos 3D.)

Como primeira etapa, proponho encontrar o casco convexo para os pontos que reformam o problema (removendo esses pontos não estão envolvidos na solução) para: ajustar uma MAR a um polígono. O método necessário fornecerá X ( centro do retângulo ), D ( duas dimensões ) e A ( ângulo ).

Minha proposta de solução:

• Encontre o centróide do polígono (consulte Localização do centro da geometria do objeto? )
• [S] Coloque um retângulo simples, isto é, paralelo aos eixos X e Y
  • você pode usar a minmaxfunção para X e Y dos pontos dados (por exemplo, vértices do polígono)
• Armazene a área do retângulo ajustado
• Gire o polígono em torno do centróide, por exemplo, 1 grau
• Repita de [S] até fazer uma rotação completa
• Relate o ângulo da área mínima como resultado

Parece-me promissor, no entanto, existem os seguintes problemas:

• escolher uma boa resolução para a mudança de ângulo pode ser desafiador,
• o custo da computação é alto,
• a solução não é analítica, mas experimental.

Sim, existe uma solução analítica para esse problema. O algoritmo que você está procurando é conhecido na generalização de polígonos como "menor retângulo circundante".

O algoritmo que você descreve é ​​bom, mas para resolver os problemas listados, você pode usar o fato de que a orientação do MAR é a mesma de uma das bordas do casco convexo da nuvem de pontos . Então, você só precisa testar as orientações das arestas convexas do casco. Você deve:

• Calcule o casco convexo da nuvem.
• Para cada aresta do casco convexo:
  • calcular a orientação da aresta (com arctan),
  • gire o casco convexo usando esta orientação para calcular facilmente a área do retângulo delimitador com min / max de x / y do casco convexo girado,
  • Armazene a orientação correspondente à área mínima encontrada,
• Retorne o retângulo correspondente à área mínima encontrada.

Um exemplo de implementação em java está disponível lá .

Em 3D, o mesmo se aplica, exceto:

• O casco convexo será um volume,
• As orientações testadas serão as orientações (em 3D) das faces convexas do casco.

Boa sorte!

Para complementar a ótima solução da @ julien, aqui está uma implementação funcional R, que poderia servir como pseudocódigo para orientar qualquer implementação específica de GIS (ou ser aplicada diretamente R, é claro). A entrada é uma matriz de coordenadas de ponto. Saída (o valor de mbr) é uma matriz dos vértices do retângulo delimitador mínimo (com o primeiro repetido para fechá-lo). Observe a ausência completa de quaisquer 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,])
}

Aqui está um exemplo de seu 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

O tempo é limitado pela velocidade do algoritmo convexo do casco, porque o número de vértices no casco é quase sempre muito menor que o total. A maioria dos algoritmos de casco convexo são assintoticamente O (n * log (n)) para n pontos: você pode calcular quase tão rápido quanto pode ler as coordenadas.

Acabei de implementar isso e postei minha resposta no StackOverflow , mas achei que deixaria minha versão aqui para que outras pessoas visualizassem:

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

Aqui estão quatro exemplos diferentes em ação. Para cada exemplo, gerei 4 pontos aleatórios e encontrei a caixa delimitadora.

Também é relativamente rápido para essas amostras em 4 pontos:

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

Há uma ferramenta no Whitebox GAT ( http://www.uoguelph.ca/~hydrogeo/Whitebox/ ) chamada Caixa Limite Mínima para resolver esse problema exato. Também existe uma ferramenta mínima de casco convexo. Várias das ferramentas da caixa de ferramentas Patch Shape, por exemplo, orientação e alongamento do patch, são baseadas na localização da caixa delimitadora mínima.

Me deparei com esse segmento enquanto procurava uma solução Python para um retângulo delimitador de área mínima.

Aqui está minha implementação , para a qual os resultados foram verificados com o Matlab.

O código de teste está incluído para polígonos simples e eu estou usando-o para encontrar a caixa delimitadora 2D mínima e as direções dos eixos para um PointCloud 3D.

Obrigado @ resposta whuber. É uma ótima solução, mas lenta para grandes nuvens de pontos. Achei que a convhullnfunção no pacote R geometryé muito mais rápida (138 s vs 0,03 s para 200000 pontos). Eu colei meus códigos aqui, pois qualquer um é interessante para uma solução mais 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)   

Dois métodos obtêm a mesma resposta (exemplo para 2000 pontos):

Simplesmente recomendo a função interna do OpenCV minAreaRect, que encontra um retângulo girado da área mínima que envolve o conjunto de pontos 2D de entrada. Para ver como usar esta função, pode-se consultar este tutorial .