25

ポリゴン(赤)の凸包(青)から派生したポイントのセット(地理座標値の黒い点)があります。図を参照してください:ここに画像の説明を入力してください

[(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でこの投稿を作成)またはこの手順例に従って、長軸と短軸の長さを計算する必要があります

ここに画像の説明を入力してください

  1. 雲の凸包を計算します。
  2. 凸包の各エッジの場合:2a。エッジの向きを計算します、2b。回転した凸包の最小/最大x/y、2cで境界長方形領域を簡単に計算するために、この方向を使用して凸包を回転させます。見つかった最小領域に対応する方向を保存します。
  3. 見つかった最小領域に対応する長方形を返します。

その後、角度シータ(画像のy軸に対する外接する長方形の方向を表します)がわかります。すべての境界点にわたるabの最小値と最大値が見つかります。

  • a(xi、yi)= xi * cos Theta + yi sin Theta
  • b(xi、yi)= xi * sin Theta + yi cos Theta

値(a_max-a_min)と(b_max-b_min)は、方向シータの外接する長方形の長さと幅をそれぞれ定義しました。

ここに画像の説明を入力してください

4

6 に答える 6

24

私はこれを自分で実装したので、他の人が見ることができるように自分のバージョンをここにドロップすると思いました。

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

私自身の参照のためにgis.stackexchangeで同じ答えにリンクしてください。

于 2015-11-09T21:57:19.940 に答える
9

一連の点の凸包内のn個の点の時計回りの順序のリストが与えられた場合、最小面積の囲み長方形を見つけるのはO(n)演算です。(凸包の検出については、O(n log n)時間で、activestate.comレシピ66527を参照するか、tixxit.netの非常にコンパクトなグラハムスキャンコードを参照してください。)

次のPythonプログラムは、凸多角形の最大直径を計算するために、通常のO(n)アルゴリズムと同様の手法を使用します。つまり、特定のベースラインに対して左端、反対側、および右端のポイントに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エントリは、ポイント10から11までのベースラインに対して、ポイント0が左端、ポイント4が反対側、ポイント7が右端であり、187.451ユニットの面積を生成することを示します。

コードはmostfar()各インデックスを進めるために使用することに注意してください。テストする極端なものを指示するパラメーターmx, mymostfar()例として、を使用するとmx,my = -1,0mostfar()-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)

注:最小領域を囲む長方形の長さと幅を取得するには、以下に示すように上記のコードを変更します。これにより、次のような出力行が生成されます

Min rectangle:  187.451   18.037   10.393   10    0    4    7

ここで、2番目と3番目の数字は長方形の長さと幅を示し、4つの整数は長方形の側面にある点のインデックス番号を示します。

# 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
于 2012-11-24T20:06:33.127 に答える
9

すでに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)

また、面積、長方形の中心、長方形の角度、およびコーナーポイントも返します。

于 2016-11-11T01:49:51.517 に答える
6

OpenCVにはこれがあります。これを参照してください:

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

7.b. 回転した長方形

cv2.minAreaRect(cnt)

長方形の長さと幅、およびその角度を取得できます。描画したい場合は、コーナーを計算することもできます。

于 2017-05-13T06:42:26.707 に答える
1

凸包を計算するレシピを見つけました。

「完全なソリューション」(すべてを実行する1つの関数)について話している場合、プログラムarcpyの一部であるものだけが見つかりました。ArcGISそれはMinimumBoundingGeometry_managementあなたが探しているもののように見える機能を提供します。しかし、それはオープンソースではありません。残念ながら、PythonオープンソースGISライブラリはありません。

于 2012-11-24T18:36:10.727 に答える
0

2013年2月に最初に投稿されました:

上記のコード例は堅牢ではありません。実際のデータ(多くの点の凸包)でテストしたところ、ほぼ正しい結果が得られました。ただし、単純な4〜6面のポリゴンの場合は機能しません。

Pythonで記述された自己完結型のソリューションは次のとおりです: https ://github.com/dbworth/minimum-area-bounding-rectangle/

凸包は、qhull(QuickHull)アルゴリズムの2D実装を使用して検出されます。解決策は、ポリゴンのすべてのエッジの角度を計算し、回転の第1象限(90度)に固有の角度でのみ操作することです。最小面積の外接長方形を見つけた後、中心点と角点を含むすべてのデータを出力します。簡単なテストプログラムが提供されています。回答はMatlabを使用して検証されました。

このソリューションは、バウンディングボックスが入力ポリゴンと少なくとも1つのエッジを共有するという仮定に依存していることに注意してください。これは私のアプリケーションに適していましたが、著者がこれが当てはまらないいくつかのまれな解決策があることを示した論文があると聞きました。この結果が重要な場合は、調査する必要があります。

于 2020-02-01T00:03:17.390 に答える