python: aiuto per implementare un algoritmo per trovare il rettangolo dell'area minima per i punti dati al fine di calcolare la lunghezza dell'asse maggiore e minore

Question

Ho un insieme di punti (punti neri nel valore di coordinate geografiche) derivato dallo scafo convesso (blu) di un poligono (rosso). vedi Figura:

[(560023.44957588764,6362057.3904932579), 
(560023.44957588764,6362060.3904932579), 
(560024.44957588764,6362063.3904932579), 
(560026.94957588764,6362068.3904932579), 
(560028.44957588764,6362069.8904932579), 
(560034.94957588764,6362071.8904932579), 
(560036.44957588764,6362071.8904932579), 
(560037.44957588764,6362070.3904932579), 
(560037.44957588764,6362064.8904932579), 
(560036.44957588764,6362063.3904932579), 
(560034.94957588764,6362061.3904932579), 
(560026.94957588764,6362057.8904932579), 
(560025.44957588764,6362057.3904932579), 
(560023.44957588764,6362057.3904932579)] 

devo calcolare la maggiore e minore lunghezza dell'asse seguente modo (formano questo post scrittura in R-progetto e nella Java) o alla this example procedure

1. Calcola lo scafo convesso del cloud.
2. Per ciascun lato dello scafo convesso: 2a. calcolare l'orientamento del bordo, 2b. ruotare lo scafo convesso usando questo orientamento per calcolare facilmente l'area del rettangolo di delimitazione con min/max di x/y dello scafo convesso ruotato, 2c. Memorizza l'orientamento corrispondente all'area minima trovata,
3. Restituisce il rettangolo corrispondente all'area minima trovata.

Dopo che sappiamo The dell'angolo Theta (rappresentato l'orientamento del rettangolo rispetto alla asse y dell'immagine). Il minimo e il massimo di un e b su tutti i punti di confine sono trovati:

• un (xi, yi) = xi * cos Theta + yi peccato Theta
• b (xi, yi) = xi * sin theta + yi cos Theta

I valori (a_max - a_min) e (b_max - b_min) definisce la lunghezza e larghezza, rispettivamente, del rettangolo di una direzione Theta.

Answer 1

Dato un elenco ordinato di n punti nello scafo convesso di un insieme di punti, è un'operazione di tipo O (n) trovare il rettangolo che racchiude l'area minima. (Per la ricerca di convesso scafo, a O (n log n), vedi activestate.com recipe 66527 o vedere il abbastanza compatto Graham scan code at tixxit.net.)

Il seguente programma Python usa tecniche simili a quelle dei soliti O (n) algoritmo per la massima calcolo diametro di un poligono convesso. Cioè, mantiene tre indici (iL, iP, iR) ai punti più a sinistra, opposti e più a destra rispetto a una data linea di riferimento. Ogni indice avanza al massimo n punti.output di esempio del programma è indicata successivo (con un'intestazione aggiunta):

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 

Ad esempio, il i = 10 entrata indica che rispetto alla linea di base dal punto 10 a 11, punto 0 è più a sinistra, il punto 4 è di fronte, e il punto 7 è più a destra, producendo un'area di 187.451 unità.

Si noti che il codice utilizza mostfar() per avanzare ogni indice. I parametri mx, my su mostfar() indicano quale estremo è necessario testare; ad esempio, con mx,my = -1,0, mostfar() proverà a massimizzare -rx (dove rx è la x ruotata di un punto), trovando così il punto più a sinistra. Si noti che probabilmente una tolleranza epsilon dovrebbe essere utilizzata quando lo if mx*rx + my*ry >= best viene eseguito in aritmetica inesatta: quando uno scafo ha numerosi punti, l'errore di arrotondamento potrebbe essere un problema e fare in modo che il metodo non anticipi in modo errato un indice.

Il codice è mostrato di seguito. I dati dello scafo sono presi dalla domanda di cui sopra, con offset grandi irrilevanti e posizioni decimali identiche eluite.

#!/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: Per ottenere la lunghezza e la larghezza della minima area racchiude rettangolo, modificare il codice precedente come illustrato di seguito. Ciò produrrà una linea di uscita come

Min rectangle: 187.451 18.037 10.393 10 0 4 7 

in cui i secondi e terzi numeri indicano la lunghezza e la larghezza del rettangolo, e le quattro interi invia numeri indice di punti che giacciono su lati di esso.

# 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 

Answer 2

ho trovato recipe to compute convex hulls.

Se stiamo parlando di "soluzioni complete" (una funzione per fare cose intere), ho trovato solo arcpy che fa parte del programma ArcGIS. Fornisce la funzione MinimumBoundingGeometry_management che sembra quella che stai cercando. Ma non è open source. Sfortunatamente, mancano le librerie GIS open source Python.

Answer 3

Ho appena implementato questo me stesso, quindi ho pensato di cadere la mia versione qui per gli altri a vista:

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 

Qui ci sono quattro diversi esempi di esso in azione. Per ogni esempio, ho generato 4 punti casuali e ho trovato il riquadro di delimitazione.

(modifica di @heltonbiker) Un codice semplice per il tracciato:

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

(fine edit)

E 'relativamente veloce troppo per questi campioni su 4 punti:

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

Link to the same answer over on gis.stackexchange per riferimento personale.

Answer 4

C'è un modulo che lo fa già su github. https://github.com/BebeSparkelSparkel/MinimumBoundingBox

Tutto quello che devi fare è inserire la tua nuvola di punti in esso.

from MinimumBoundingBox import minimum_bounding_box 
points = ((1,2), (5,4), (-1,-3)) 
bounding_box = minimum_bounding_box(points) # returns namedtuple 

è possibile ottenere maggiori e minori lunghezze asse da:

minor = min(bounding_box.length_parallel, bounding_box.length_orthogonal) 
major = max(bounding_box.length_parallel, bounding_box.length_orthogonal) 

Inoltre restituisce zona, centro del rettangolo, angolo di rettangolo e punti d'angolo.

Answer 5

OpenCV ha questo. Vedere questo:

http://docs.opencv.org/trunk/dd/d49/tutorial_py_contour_features.html

7.b. Ruotato rettangolo

Con

cv2.minAreaRect(cnt) 

È possibile ottenere la lunghezza e la larghezza del rettangolo, così come il suo angolo. Puoi anche calcolare gli angoli se desideri disegnarli.