私はnumpyで座標の点群を持っている。点の数が多い場合、点が点群の凸包にあるかどうかを調べたい。
pyhullを試してみましたが、点が ConvexHull
にあるかどうかを調べる方法がわかりません:
hull = ConvexHull(np.array([(1, 2), (3, 4), (3, 6)]))
for s in hull.simplices:
s.in_simplex(np.array([2, 3]))
raises LinAlgError:配列は正方形でなければなりません。
scipyだけでできる簡単な解決策を紹介しよう:
def in_hull(p, hull):
"""
Test if points in `p` are in `hull`
`p` should be a `NxK` coordinates of `N` points in `K` dimensions
`hull` is either a scipy.spatial.Delaunay object or the `MxK` array of the
coordinates of `M` points in `K`dimensions for which Delaunay triangulation
will be computed
"""
from scipy.spatial import Delaunay
if not isinstance(hull,Delaunay):
hull = Delaunay(hull)
return hull.find_simplex(p)>=0
これはブール値の配列を返し, True
の値は与えられた凸包内にある点を示す.これは次のように使うことができる.
tested = np.random.rand(20,3)
cloud = np.random.rand(50,3)
print in_hull(tested,cloud)
matplotlib がインストールされていれば、最初の関数を呼び出して結果をプロットする以下の関数も利用できます。Nx2`の配列で与えられる2次元データのみ:
def plot_in_hull(p, hull):
"""
plot relative to `in_hull` for 2d data
"""
import matplotlib.pyplot as plt
from matplotlib.collections import PolyCollection, LineCollection
from scipy.spatial import Delaunay
if not isinstance(hull,Delaunay):
hull = Delaunay(hull)
# plot triangulation
poly = PolyCollection(hull.points[hull.vertices], facecolors='w', edgecolors='b')
plt.clf()
plt.title('in hull')
plt.gca().add_collection(poly)
plt.plot(hull.points[:,0], hull.points[:,1], 'o', hold=1)
# plot the convex hull
edges = set()
edge_points = []
def add_edge(i, j):
"""Add a line between the i-th and j-th points, if not in the list already"""
if (i, j) in edges or (j, i) in edges:
# already added
return
edges.add( (i, j) )
edge_points.append(hull.points[ [i, j] ])
for ia, ib in hull.convex_hull:
add_edge(ia, ib)
lines = LineCollection(edge_points, color='g')
plt.gca().add_collection(lines)
plt.show()
# plot tested points `p` - black are inside hull, red outside
inside = in_hull(p,hull)
plt.plot(p[ inside,0],p[ inside,1],'.k')
plt.plot(p[-inside,0],p[-inside,1],'.r')
2Dの点群を使っているようなので、凸ポリゴンのポリゴン内点群検定用の包含検定を紹介します。
Scipyの凸包アルゴリズムでは、2次元以上の凸包を求めることができますが、これは2次元の点群に対して必要以上に複雑です。したがって、こちらのような別のアルゴリズムを使用することをお勧めします。なぜなら、凸包のポリゴン内点検定に本当に必要なのは、時計回りの順序で並んだ凸包点のリストと、ポリゴンの内側にある点だからである。
この手法の時間的性能は以下の通りである:
ここでNは点群内の点数、hは点群凸包内の点数である。
scipyを使い続けたいのであれば、凸包を作る必要がある。
>>> from scipy.spatial import ConvexHull
>>> points = np.random.rand(30, 2) # 30 random points in 2-D
>>> hull = ConvexHull(points)
を作成し、外皮上の点のリストを作成します。以下は外皮をプロットするコードです。
>>> import matplotlib.pyplot as plt
>>> plt.plot(points[:,0], points[:,1], 'o')
>>> for simplex in hull.simplices:
>>> plt.plot(points[simplex,0], points[simplex,1], 'k-')
このコードから始めて、私は船体上の点のリストを計算するために次のことを提案します。
pts_hull = [(points[simplex,0], points[simplex,1])
for simplex in hull.simplices]
(私は試していませんが)。
また、船体を計算し、x,y点を返す独自のコードを用意することもできます。
元のデータセットの点が外皮の上にあるかどうかを知りたいのであれば、これで完了だ。
ある点が外皮の内側にあるか外側にあるかを知りたい場合、もう少し工夫が必要です。あなたがしなければならないことは
船体の2つの単純点を結ぶすべての辺について、その点が船体の上にあるか下にあるかを決定する。
点がすべての線より下にあるか、すべての線より上にあれば、それは船体の外側にある。
スピードアップとして、ある点がある線より上にあり、別の線より下にある場合、その点は船体の内側にある。