python：長軸と短軸の長さを計算するために与えられた点の最小面積矩形を見つけるアルゴリズムを実装するのを手助けします

Question

私は凸包から派生した点のセット（地理座標値の黒い点） （赤色）のポリゴン（青色）になります。図は、以下を参照してください

私は以下の手順に従ってメジャーとマイナーの軸の長さを計算する必要が

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

（R-プロジェクト内やJavaでこのpost書き込みを形成する）、または以下のthis example procedure

1. 雲の凸包を計算します。
2. 凸包の各辺について： 2a。エッジ方向、 2bを計算します。この向きを使用して凸包を回転させて、回転した凸包の最小/最大x/yで境界矩形領域を簡単に計算します（ 2c）。見つかった最小領域に対応する方向を格納します。
3. 見つかった最小領域に対応する矩形を返します。その後

我々は、（画像のy軸に外接する四角形の向きを表す）角度シータを知っています。シータ+のYI罪シータ

B（XI、YI）COS（XI、YI）= XI *

• ：の最小値と最大値は、すべての境界点上及びBが見出さ あります=シータ
COS XI *罪シータ+のYI

、値（a_max - a_min）及び（b_maxの - b_min）は、それぞれ、方向シータの境界矩形の の長さと幅を定義しました。

Answer 1

点集合の凸包内のn個の点の時計回りのリストが与えられた場合、それは最小面積を囲む矩形を見つけるためのO（n）操作です。

次のpythonプログラムでは、通常のO（n）アルゴリズムと同様の手法を使用して最大値を計算します（最大値を計算するには、O（n log n）時間内に凸包を見つけるにはactivestate.com recipe 66527を参照してください）。凸多角形の直径。つまり、与えられたベースラインに対して左端、反対端、右端の3つのインデックス（iL、iP、iR）を維持します。各インデックスは、最大でnポイント進んでいます。プログラムからのサンプル出力を（追加ヘッダで）次示されている：例えば

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 

、I = 10のエントリは、ポイント4は、ポイント10から11までのベースラインに対して、点0は左端であることを示し反対側にあり、ポイント7は右端にあり、187.451台の面積が得られます。

コードでは、mostfar()を使用して各インデックスを進めることに注意してください。 mx, myのパラメータをmostfar()に設定すると、テストする極端なことが分かります。一例として、mx,my = -1,0でmostfar()は、-rx（ここで、rxは点の回転したx）を最大にしようとします。したがって、最も左の点が見つかります。ただし、浮動小数点数が多い場合は、丸め誤差が問題となり、誤って指数を進ませない可能性がありますので、不正確な算術演算でif mx*rx + my*ry >= bestを実行した場合は、εの余裕が必要です。

コードを以下に示します。船体のデータは上記の質問から得られ、無関係な大きなオフセットと同一の小数点以下の桁は省略されています。

#!/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) 

注：が最小面積外接矩形の長さ及び幅を得るために、以下に示すように上記のコードを変更します。これは、第二及び第三の数字は長方形の長さおよび幅を示し、4つの整数はその側面上に位置する点のインデックス番号を与えている

Min rectangle: 187.451 18.037 10.393 10 0 4 7 

等出力線を生成します。

# 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

私はrecipe to compute convex hullsを見つけました。

"完全な解決策"（1つの機能で全体を処理する）について言えば、ArcGISプログラムの一部であるarcpyしか見つかりませんでした。それはあなたが探しているように見えるMinimumBoundingGeometry_management機能を提供します。しかし、オープンソースではありません。残念ながら、PythonのオープンソースのGISライブラリが欠けています。

Answer 3

私はちょうど自分自身はこれを実装したので、私は、私はビューに他の人のために、ここで私のバージョンをドロップしたい考え出し：ここ

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 

はアクションでそれ4つの異なる例があります。各例では、4つのランダムな点が生成され、境界ボックスが見つかりました。

（@heltonbikerによって編集） プロットするための簡単なコード：

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

（エンド編集）

それは4ポイントで、これらのサンプルのためにあまりにも比較的速いです：

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

Link to the same answer over on gis.stackexchange私の自身の参照のために。

Answer 4

githubには既にこれを行うモジュールがあります。 https://github.com/BebeSparkelSparkel/MinimumBoundingBox

ポイントクラウドを挿入するだけです。

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

次の方法でメジャーとマイナーの軸の長さを得ることができます。また、地域、長方形センター、長方形角、およびコーナーポイントを返し

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

。

Answer 5

OpenCVがこれです。

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

7.b.：これを参照してください。

cv2.minAreaRect(cnt) 

で回転四角形は

あなたは長さと幅の長方形のと同様にその角度を得ることができます。また、描画する場合は、角を計算することもできます。