Solución:
Acabo de implementar esto yo mismo, así que 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.
(editado por @heltonbiker) Un código simple para trazar:
import matplotlib.pyplot as plt
for n in range(10):
points = np.random.rand(4,2)
plt.scatter(points[:,0], points[:,1])
bbox = minimum_bounding_rectangle(points)
plt.fill(bbox[:,0], bbox[:,1], alpha=0.2)
plt.axis('equal')
plt.show()
(fin de editar)
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
Enlace a la misma respuesta en gis.stackexchange para mi propia referencia.
Dada una lista ordenada en el sentido de las agujas del reloj de n puntos en el casco convexo de un conjunto de puntos, es una operación O (n) encontrar el rectángulo circundante de área mínima. (Para encontrar cascos convexos, en tiempo O (n log n), vea la receta 66527 de activestate.com o vea el código de escaneo de Graham bastante compacto en tixxit.net).
El siguiente programa de Python utiliza técnicas similares a las del algoritmo O (n) habitual para calcular el diámetro máximo de un polígono convexo. Es decir, mantiene tres índices (iL, iP, iR) en los puntos más a la izquierda, opuestos y más a la derecha en relación con una línea de base determinada. Cada índice avanza como máximo n puntos. La salida de muestra del programa se muestra a continuación (con un encabezado agregado):
i iL iP iR Area
0 6 8 0 203.000
1 6 8 0 211.875
2 6 8 0 205.800
3 6 10 0 206.250
4 7 12 0 190.362
5 8 0 1 203.000
6 10 0 4 201.385
7 0 1 6 203.000
8 0 3 6 205.827
9 0 3 6 205.640
10 0 4 7 187.451
11 0 4 7 189.750
12 1 6 8 203.000
Por ejemplo, la entrada i = 10 indica que en relación con la línea de base del punto 10 al 11, el punto 0 está en el extremo izquierdo, el punto 4 es el opuesto y el punto 7 está en el extremo derecho, lo que arroja un área de 187.451 unidades.
Tenga en cuenta que el código utiliza mostfar() para avanzar cada índice. los mx, my parámetros a mostfar() dígale qué extremo probar; como ejemplo, con mx,my = -1,0, mostfar() intentará maximizar -rx (donde rx es la x rotada de un punto), encontrando así el punto más a la izquierda. Tenga en cuenta que probablemente debería utilizarse una asignación épsilon cuando if mx*rx + my*ry >= best se realiza en aritmética inexacta: cuando un casco tiene numerosos puntos, el error de redondeo puede ser un problema y hacer que el método no avance incorrectamente un índice.
El código se muestra a continuación. Los datos del casco se toman de la pregunta anterior, con grandes desplazamientos irrelevantes y lugares decimales idénticos omitidos.
#!/usr/bin/python
import math
hull = [(23.45, 57.39), (23.45, 60.39), (24.45, 63.39),
(26.95, 68.39), (28.45, 69.89), (34.95, 71.89),
(36.45, 71.89), (37.45, 70.39), (37.45, 64.89),
(36.45, 63.39), (34.95, 61.39), (26.95, 57.89),
(25.45, 57.39), (23.45, 57.39)]
def mostfar(j, n, s, c, mx, my): # advance j to extreme point
xn, yn = hull[j][0], hull[j][1]
rx, ry = xn*c - yn*s, xn*s + yn*c
best = mx*rx + my*ry
while True:
x, y = rx, ry
xn, yn = hull[(j+1)%n][0], hull[(j+1)%n][1]
rx, ry = xn*c - yn*s, xn*s + yn*c
if mx*rx + my*ry >= best:
j = (j+1)%n
best = mx*rx + my*ry
else:
return (x, y, j)
n = len(hull)
iL = iR = iP = 1 # indexes left, right, opposite
pi = 4*math.atan(1)
for i in range(n-1):
dx = hull[i+1][0] - hull[i][0]
dy = hull[i+1][1] - hull[i][1]
theta = pi-math.atan2(dy, dx)
s, c = math.sin(theta), math.cos(theta)
yC = hull[i][0]*s + hull[i][1]*c
xP, yP, iP = mostfar(iP, n, s, c, 0, 1)
if i==0: iR = iP
xR, yR, iR = mostfar(iR, n, s, c, 1, 0)
xL, yL, iL = mostfar(iL, n, s, c, -1, 0)
area = (yP-yC)*(xR-xL)
print ' {:2d} {:2d} {:2d} {:2d} {:9.3f}'.format(i, iL, iP, iR, area)
Nota: Para obtener la longitud y el ancho del rectángulo circundante de área mínima, modifique el código anterior como se muestra a continuación. Esto producirá una línea de salida como
Min rectangle: 187.451 18.037 10.393 10 0 4 7
en el que los números segundo y tercero indican la longitud y el ancho del rectángulo, y los cuatro números enteros dan números de índice de los puntos que se encuentran a los lados del mismo.
# add after pi = ... line:
minRect = (1e33, 0, 0, 0, 0, 0, 0) # area, dx, dy, i, iL, iP, iR
# add after area = ... line:
if area < minRect[0]:
minRect = (area, xR-xL, yP-yC, i, iL, iP, iR)
# add after print ... line:
print 'Min rectangle:', minRect
# or instead of that print, add:
print 'Min rectangle: ',
for x in ['{:3d} '.format(x) if isinstance(x, int) else '{:7.3f} '.format(x) for x in minRect]:
print x,
print
Hay un módulo que ya hace esto en github. https://github.com/BebeSparkelSparkel/MinimumBoundingBox
Todo lo que necesita hacer es insertar su nube de puntos en él.
from MinimumBoundingBox import minimum_bounding_box
points = ( (1,2), (5,4), (-1,-3) )
bounding_box = minimum_bounding_box(points) # returns namedtuple
Puede obtener las longitudes de los ejes mayor y menor mediante:
minor = min(bounding_box.length_parallel, bounding_box.length_orthogonal)
major = max(bounding_box.length_parallel, bounding_box.length_orthogonal)
También devuelve el área, el centro del rectángulo, el ángulo del rectángulo y los puntos de las esquinas.