一連の点の凸包内の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, 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)
注:最小領域を囲む長方形の長さと幅を取得するには、以下に示すように上記のコードを変更します。これにより、次のような出力行が生成されます
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