Skip to content

Algorithmus zur Ermittlung des minimalen Flächenwinkels für gegebene Punkte, um die Länge der Haupt- und Nebenachse zu berechnen

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.

examples

(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.



Nutzen Sie unsere Suchmaschine

Suche
Generic filters

Schreibe einen Kommentar

Deine E-Mail-Adresse wird nicht veröffentlicht.