Nach Rücksprache mit Spezialisten auf diesem Gebiet, Programmierern verschiedener Branchen und Professoren haben wir die Antwort auf das Dilemma gefunden und lassen sie in diesem Beitrag widerspiegeln.
Lösung:
Ich habe dies gerade selbst implementiert, also dachte ich mir, ich stelle meine Version hier ein, damit andere sie sehen können:
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
Hier sind vier verschiedene Beispiele in Aktion. Für jedes Beispiel habe ich 4 zufällige Punkte erzeugt und die Bounding Box gefunden.
(edit by @heltonbiker)
Ein einfacher Code zum Plotten:
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()
(end edit)
Es geht auch relativ schnell für diese Beispiele an 4 Punkten:
>>> %timeit minimum_bounding_rectangle(a)
1000 loops, best of 3: 245 µs per loop
Link zur gleichen Antwort drüben auf gis.stackexchange für meine eigene Referenz.
Bei einer im Uhrzeigersinn geordneten Liste von n Punkten in der konvexen Hülle eines Punktesatzes ist es eine O(n)-Operation, das umschließende Rechteck mit minimalem Flächeninhalt zu finden. (Für die Suche nach einer konvexen Hülle in O(n log n)-Zeit siehe Rezept 66527 von activestate.com oder den recht kompakten Graham-Scan-Code auf tixxit.net).
Das folgende Python-Programm verwendet Techniken, die denen des üblichen O(n)-Algorithmus zur Berechnung des maximalen Durchmessers eines konvexen Polygons ähneln. Das bedeutet, dass es drei Indizes (iL, iP, iR) für den äußersten linken, den gegenüberliegenden und den äußersten rechten Punkt relativ zu einer gegebenen Grundlinie verwaltet. Jeder Index geht durch höchstens n Punkte. Das folgende Beispiel zeigt die Ausgabe des Programms (mit einer zusätzlichen Kopfzeile):
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
Der Eintrag i=10 zeigt zum Beispiel an, dass relativ zur Grundlinie von Punkt 10 bis 11 der Punkt 0 ganz links, der Punkt 4 gegenüber und der Punkt 7 ganz rechts ist, was eine Fläche von 187,451 Einheiten ergibt.
Beachten Sie, dass der Code mostfar() verwendet, um jeden Index weiterzuschalten. Die mx, my Parameter zu mostfar() geben an, nach welchem Extremwert getestet werden soll; zum Beispiel mit mx,my = -1,0, mostfar() versucht, -rx zu maximieren (wobei rx das gedrehte x eines Punktes ist) und findet so den Punkt ganz links. Beachten Sie, dass ein Epsilon-Abzug wahrscheinlich verwendet werden sollte, wenn if mx*rx + my*ry >= best in ungenauer Arithmetik ausgeführt wird: Wenn eine Hülle zahlreiche Punkte hat, können Rundungsfehler ein Problem darstellen und dazu führen, dass die Methode fälschlicherweise einen Index nicht weiterführt.
Der Code ist unten dargestellt. Die Daten der Hülle stammen aus der obigen Frage, wobei irrelevante große Offsets und identische Nachkommastellen weggelassen wurden.
#!/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)
Anmerkung: Um die Länge und Breite des umschließenden Rechtecks mit minimaler Fläche zu erhalten, ändern Sie den obigen Code wie unten gezeigt. Dies ergibt eine Ausgabezeile wie
Min rectangle: 187.451 18.037 10.393 10 0 4 7
in der die zweite und dritte Zahl die Länge und Breite des Rechtecks angeben und die vier ganzen Zahlen die Indexnummern der Punkte, die auf den Seiten des Rechtecks liegen.
# 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
Es gibt bereits ein Modul auf github, das dies tut.
https://github.com/BebeSparkelSparkel/MinimumBoundingBox
Sie müssen nur noch Ihre Punktwolke in dieses Modul einfügen.
from MinimumBoundingBox import minimum_bounding_box
points = ( (1,2), (5,4), (-1,-3) )
bounding_box = minimum_bounding_box(points) # returns namedtuple
Die Längen der Haupt- und Nebenachsen erhält man durch:
minor = min(bounding_box.length_parallel, bounding_box.length_orthogonal)
major = max(bounding_box.length_parallel, bounding_box.length_orthogonal)
Außerdem werden Fläche, Rechteckmitte, Rechteckwinkel und Eckpunkte zurückgegeben.
Wenn Sie können, können Sie einen Aufsatz darüber schreiben, was Sie von dieser Aussage halten.