2 セットの等高線の間のすべての交点 (丸め誤差) を見つける最良の方法について疑問に思っています。最良の方法はどれですか? 次に例を示します。
import matplotlib.pyplot as plt
import numpy as np
x = np.linspace(-1,1,500)
X,Y = np.meshgrid(x,x)
Z1 = np.abs(np.sin(2*X**2+Y))
Z2 = np.abs(np.cos(2*Y**2+X**2))
plt.contour(Z1,colors='k')
plt.contour(Z2,colors='r')
plt.show()
次のようなものが欲しい：
intersection_points = intersect(contour1,contour2)
print intersection_points
[(x1,y1),...,(xn,yn)]
import collections
import matplotlib.pyplot as plt
import numpy as np
import scipy.spatial as spatial
import scipy.spatial.distance as dist
import scipy.cluster.hierarchy as hier
def intersection(points1, points2, eps):
tree = spatial.KDTree(points1)
distances, indices = tree.query(points2, k=1, distance_upper_bound=eps)
intersection_points = tree.data[indices[np.isfinite(distances)]]
return intersection_points
def cluster(points, cluster_size):
dists = dist.pdist(points, metric='sqeuclidean')
linkage_matrix = hier.linkage(dists, 'average')
groups = hier.fcluster(linkage_matrix, cluster_size, criterion='distance')
return np.array([points[cluster].mean(axis=0)
for cluster in clusterlists(groups)])
def contour_points(contour, steps=1):
return np.row_stack([path.interpolated(steps).vertices
for linecol in contour.collections
for path in linecol.get_paths()])
def clusterlists(T):
'''
http://stackoverflow.com/a/2913071/190597 (denis)
T = [2, 1, 1, 1, 2, 2, 2, 2, 2, 1]
Returns [[0, 4, 5, 6, 7, 8], [1, 2, 3, 9]]
'''
groups = collections.defaultdict(list)
for i, elt in enumerate(T):
groups[elt].append(i)
return sorted(groups.values(), key=len, reverse=True)
# every intersection point must be within eps of a point on the other
# contour path
eps = 1.0
# cluster together intersection points so that the original points in each flat
# cluster have a cophenetic_distance < cluster_size
cluster_size = 100
x = np.linspace(-1, 1, 500)
X, Y = np.meshgrid(x, x)
Z1 = np.abs(np.sin(2 * X ** 2 + Y))
Z2 = np.abs(np.cos(2 * Y ** 2 + X ** 2))
contour1 = plt.contour(Z1, colors='k')
contour2 = plt.contour(Z2, colors='r')
points1 = contour_points(contour1)
points2 = contour_points(contour2)
intersection_points = intersection(points1, points2, eps)
intersection_points = cluster(intersection_points, cluster_size)
plt.scatter(intersection_points[:, 0], intersection_points[:, 1], s=20)
plt.show()
収量
@unutbu の回答と、numpy-and-line-intersectionsから直接抽出された交差アルゴリズムから作業して、これを思いつきました。ループ内のループ内の交差点とループ内のループを見つけるためのブルートフォースのような方法のために遅いです。ループを高速化する方法があるかもしれませんが、交差アルゴリズムについてはわかりません。
import matplotlib.pyplot as plt
import numpy as np
def find_intersections(A, B):
#this function stolen from https://stackoverflow.com/questions/3252194/numpy-and-line-intersections#answer-9110966
# min, max and all for arrays
amin = lambda x1, x2: np.where(x1<x2, x1, x2)
amax = lambda x1, x2: np.where(x1>x2, x1, x2)
aall = lambda abools: np.dstack(abools).all(axis=2)
slope = lambda line: (lambda d: d[:,1]/d[:,0])(np.diff(line, axis=0))
x11, x21 = np.meshgrid(A[:-1, 0], B[:-1, 0])
x12, x22 = np.meshgrid(A[1:, 0], B[1:, 0])
y11, y21 = np.meshgrid(A[:-1, 1], B[:-1, 1])
y12, y22 = np.meshgrid(A[1:, 1], B[1:, 1])
m1, m2 = np.meshgrid(slope(A), slope(B))
m1inv, m2inv = 1/m1, 1/m2
yi = (m1*(x21-x11-m2inv*y21) + y11)/(1 - m1*m2inv)
xi = (yi - y21)*m2inv + x21
xconds = (amin(x11, x12) < xi, xi <= amax(x11, x12),
amin(x21, x22) < xi, xi <= amax(x21, x22) )
yconds = (amin(y11, y12) < yi, yi <= amax(y11, y12),
amin(y21, y22) < yi, yi <= amax(y21, y22) )
return xi[aall(xconds)], yi[aall(yconds)]
x = np.linspace(-1,1,500)
X,Y = np.meshgrid(x,x)
Z1 = np.abs(np.sin(2*X**2+Y))
Z2 = np.abs(np.cos(2*Y**2+X**2))
contour1 = plt.contour(Z1,colors='k')
contour2 = plt.contour(Z2,colors='r')
xi = np.array([])
yi = np.array([])
for linecol in contour1.collections:
for path in linecol.get_paths():
for linecol2 in contour2.collections:
for path2 in linecol2.get_paths():
xinter, yinter = find_intersections(path.vertices, path2.vertices)
xi = np.append(xi, xinter)
yi = np.append(yi, yinter)
plt.scatter(xi, yi, s=20)
plt.show()
編集： find_intersections上記は、垂直線と水平線、または重なっている線分を処理しません。以下のlinelineintersect関数はこれらのケースを処理しますが、プロセス全体はまだ遅いです。カウントを追加したので、あとどれくらいかかるかがわかります。また、軸と交点がグリッド ポイントの数ではなく、グリッドの実際の X および Y 座標にあるように 変更contour1 = plt.contour(Z1,colors='k')しました。contour1 = plt.contour(X,Y,Z1,colors='k')
問題は、大量のデータがあることだと思います。1 つのコンター セットには合計 12818 のライン セグメントを持つ 24 本のラインがあり、もう 1 つのコンター セットには 11411 のライン セグメントを持つ 29 本のラインがあります。これは、チェックする線分の組み合わせがたくさんあります ( への 696 回の呼び出しlinelineintersect)。より粗いグリッドを使用して関数を計算すると、より迅速に結果を得ることができます (100 x 100 は、元の 500 x 500 よりもはるかに高速でした)。ある種のライン スイープ アルゴリズムを実行できるかもしれませんが、そのような実装方法はわかりません。線の交差の問題は、コンピューター ゲームで多くの用途があり、衝突検出として知られています。すべての交点をすばやく特定できるグラフィックス ライブラリが存在するはずです。
import numpy as np
from numpy.lib.stride_tricks import as_strided
import matplotlib.pyplot as plt
import itertools
def unique(a, atol=1e-08):
"""Find unique rows in 2d array
Parameters
----------
a : 2d ndarray, float
array to find unique rows in
atol : float, optional
tolerance to check uniqueness
Returns
-------
out : 2d ndarray, float
array of unique values
Notes
-----
Adapted to include tolerance from code at https://stackoverflow.com/questions/8560440/removing-duplicate-columns-and-rows-from-a-numpy-2d-array#answer-8564438
"""
if np.issubdtype(a.dtype, float):
order = np.lexsort(a.T)
a = a[order]
diff = np.diff(a, axis=0)
np.abs(diff,out = diff)
ui = np.ones(len(a), 'bool')
ui[1:] = (diff > atol).any(axis=1)
return a[ui]
else:
order = np.lexsort(a.T)
a = a[order]
diff = np.diff(a, axis=0)
ui = np.ones(len(a), 'bool')
ui[1:] = (diff != 0).any(axis=1)
return a[ui]
def linelineintersect(a, b, atol=1e-08):
"""Find all intersection points of two lines defined by series of x,y pairs
Intersection points are unordered
Colinear lines that overlap intersect at any end points that fall within the overlap
Parameters
----------
a, b : ndarray
2 column ndarray of x,y values defining a two dimensional line. 1st
column is x values, 2nd column is x values.
Notes
-----
An example of 3 segment line is: a = numpy.array([[0,0],[5.0,3.0],[4.0,1]])
Function faster when there are no overlapping line segments
"""
def x_y_from_line(line):
"""1st and 2nd column of array"""
return (line[:, 0],line[:, 1])
def meshgrid_as_strided(x, y, mask=None):
"""numpy.meshgrid without copying data (using as_strided)"""
if mask is None:
return (as_strided(x, strides=(0, x.strides[0]), shape=(y.size, x.size)),
as_strided(y, strides=(y.strides[0],0), shape=(y.size, x.size)))
else:
return (np.ma.array(as_strided(x, strides=(0, x.strides[0]), shape=(y.size, x.size)), mask=mask),
np.ma.array(as_strided(y, strides=(y.strides[0],0), shape=(y.size, x.size)), mask=mask))
#In the following the indices i, j represents the pairing of the ith segment of b and the jth segment of a
#e.g. if ignore[i,j]==True then the ith segment of b and the jth segment of a cannot intersect
ignore = np.zeros([b.shape[0]-1, a.shape[0]-1], dtype=bool)
x11, x21 = meshgrid_as_strided(a[:-1, 0], b[:-1, 0], mask=ignore)
x12, x22 = meshgrid_as_strided(a[1: , 0], b[1: , 0], mask=ignore)
y11, y21 = meshgrid_as_strided(a[:-1, 1], b[:-1, 1], mask=ignore)
y12, y22 = meshgrid_as_strided(a[1: , 1], b[1: , 1], mask=ignore)
#ignore segements with non-overlappiong bounding boxes
ignore[np.ma.maximum(x11, x12) < np.ma.minimum(x21, x22)] = True
ignore[np.ma.minimum(x11, x12) > np.ma.maximum(x21, x22)] = True
ignore[np.ma.maximum(y11, y12) < np.ma.minimum(y21, y22)] = True
ignore[np.ma.minimum(y11, y12) > np.ma.maximum(y21, y22)] = True
#find intersection of segments, ignoring impossible line segment pairs when new info becomes available
denom_ = np.empty(ignore.shape, dtype=float)
denom = np.ma.array(denom_, mask=ignore)
denom_[:, :] = ((y22 - y21) * (x12 - x11)) - ((x22 - x21) * (y12 - y11))
ua_ = np.empty(ignore.shape, dtype=float)
ua = np.ma.array(ua_, mask=ignore)
ua_[:, :] = (((x22 - x21) * (y11 - y21)) - ((y22 - y21) * (x11 - x21)))
ua_[:, :] /= denom
ignore[ua < 0] = True
ignore[ua > 1] = True
ub_ = np.empty(ignore.shape, dtype=float)
ub = np.ma.array(ub_, mask=ignore)
ub_[:, :] = (((x12 - x11) * (y11 - y21)) - ((y12 - y11) * (x11 - x21)))
ub_[:, :] /= denom
ignore[ub < 0] = True
ignore[ub > 1] = True
nans_ = np.zeros(ignore.shape, dtype = bool)
nans = np.ma.array(nans_, mask = ignore)
nans_[:,:] = np.isnan(ua)
if not np.ma.any(nans):
#remove duplicate cases where intersection happens on an endpoint
# ignore[np.ma.where((ua[:, :-1] == 1) & (ua[:, 1:] == 0))] = True
# ignore[np.ma.where((ub[:-1, :] == 1) & (ub[1:, :] == 0))] = True
ignore[np.ma.where((ua[:, :-1] < 1.0 + atol) & (ua[:, :-1] > 1.0 - atol) & (ua[:, 1:] < atol) & (ua[:, 1:] > -atol))] = True
ignore[np.ma.where((ub[:-1, :] < 1 + atol) & (ub[:-1, :] > 1 - atol) & (ub[1:, :] < atol) & (ub[1:, :] > -atol))] = True
xi = x11 + ua * (x12 - x11)
yi = y11 + ua * (y12 - y11)
return np.ma.compressed(xi), np.ma.compressed(yi)
else:
n_nans = np.ma.sum(nans)
n_standard = np.ma.count(x11) - n_nans
#I initially tried using a set to get unique points but had problems with floating point equality
#create n by 2 array to hold all possible intersection points, check later for uniqueness
points = np.empty([n_standard + 2 * n_nans, 2], dtype = float) #each colinear segment pair has two intersections, hence the extra n_colinear points
#add standard intersection points
xi = x11 + ua * (x12 - x11)
yi = y11 + ua * (y12 - y11)
points[:n_standard,0] = np.ma.compressed(xi[~nans])
points[:n_standard,1] = np.ma.compressed(yi[~nans])
ignore[~nans]=True
#now add the appropriate intersection points for the colinear overlapping segments
#colinear overlapping segments intersect on the two internal points of the four points describing a straight line
#ua and ub have already serverd their purpose. Reuse them here to mean:
#ua is relative position of 1st b segment point along segment a
#ub is relative position of 2nd b segment point along segment a
use_x = x12 != x11 #avoid vertical lines diviiding by zero
use_y = ~use_x
ua[use_x] = (x21[use_x]-x11[use_x]) / (x12[use_x]-x11[use_x])
ua[use_y] = (y21[use_y]-y11[use_y]) / (y12[use_y]-y11[use_y])
ub[use_x] = (x22[use_x]-x11[use_x]) / (x12[use_x]-x11[use_x])
ub[use_y] = (y22[use_y]-y11[use_y]) / (y12[use_y]-y11[use_y])
np.ma.clip(ua, a_min=0,a_max = 1, out = ua)
np.ma.clip(ub, a_min=0,a_max = 1, out = ub)
xi = x11 + ua * (x12 - x11)
yi = y11 + ua * (y12 - y11)
points[n_standard:n_standard + n_nans,0] = np.ma.compressed(xi)
points[n_standard:n_standard + n_nans,1] = np.ma.compressed(yi)
xi = x11 + ub * (x12 - x11)
yi = y11 + ub * (y12 - y11)
points[n_standard+n_nans:,0] = np.ma.compressed(xi)
points[n_standard+n_nans:,1] = np.ma.compressed(yi)
points_unique = unique(points)
return points_unique[:,0], points_unique[:,1]
x = np.linspace(-1,1,500)
X,Y = np.meshgrid(x,x)
Z1 = np.abs(np.sin(2*X**2+Y))
Z2 = np.abs(np.cos(2*Y**2+X**2))
contour1 = plt.contour(X,Y,Z1,colors='k')
contour2 = plt.contour(X,Y,Z2,colors='r')
xi = np.array([])
yi = np.array([])
i=0
ncombos = (sum([len(x.get_paths()) for x in contour1.collections]) *
sum([len(x.get_paths()) for x in contour2.collections]))
for linecol1, linecol2 in itertools.product(contour1.collections, contour2.collections):
for path1, path2 in itertools.product(linecol1.get_paths(),linecol2.get_paths()):
i += 1
print('line combo %5d of %5d' % (i, ncombos))
xinter, yinter = linelineintersect(path1.vertices, path2.vertices)
xi = np.append(xi, xinter)
yi = np.append(yi, yinter)
plt.scatter(xi, yi, s=20)
plt.show()
このプロットは、実際の X 軸と Y 軸を持つ 50x50 グリッド用です。
