İki kontur seti arasındaki tüm kesişim noktalarını verimli bir şekilde nasıl bulursunuz

Question

İki küme çizgisi çizgisi arasındaki tüm kesişim noktalarını (dönüş hatası) bulmak için en iyi yolu merak ediyorum. Wich en iyi yöntem mi? İşte örnek:

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

istediğim bazı benzer:

intersection_points = intersect(contour1,contour2) 
print intersection_points 
[(x1,y1),...,(xn,yn)] 

Answer 1

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

verim

Answer 2

@unutbu'nun yanıtından ve numpy-and-line-intersections dan koparılan bir kesişim algoritmasından yola çıkarak bunu ortaya çıkardım. Kesişmeleri ve ilmeklerin ilmek içindeki ilmekleri bulmak için kaba kuvvet türünden dolayı yavaştır. Döngüleri daha hızlı yapmanın bir yolu olabilir ama kesişim algoritmasından emin değilim.

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

Düzenleme:find_intersections yukarıda dikey ve yatay çizgiler veya üst üste gelen çizgi parçası işlemez. Aşağıdaki linelineintersect işlevi, bu durumları ele alır, ancak tüm süreç hala yavaştır. Bir sayımı ekledim, böylece ne kadar süreceğine dair bir fikrin var. Ayrıca contour1 = plt.contour(Z1,colors='k')'u contour1 = plt.contour(X,Y,Z1,colors='k') olarak değiştirdim, böylece eksenler ve kesişim noktalarınız ızgara noktalarınızın sayısından ziyade ızgaranızın gerçek X ve Y koordinatlarında.

Sanırım sorun sadece çok fazla veriye sahip olmanız. Bir kontur setinde toplam 12818 hat segmenti olan 24 hat vardır, diğer kontur seti 11411 hat segmenti ile 29 satıra sahiptir. Bu kontrol etmek için çok sayıda hat segmenti kombinasyonu (linelineintersect'a 696 çağrı). Fonksiyonlarınızı hesaplamak için daha kaba bir ızgara kullanarak daha hızlı sonuç alabilirsiniz (100 ile 100, 500 ile 500 arasında çok daha hızlıydı). Bir çeşit çizgi süpürme algoritması yapabilirsin ama böyle şeyleri nasıl uygulayacağımı bilmiyorum. Çizgi kesişim problemi bilgisayar oyunlarında birçok uygulamaya sahiptir ve çarpışma tespiti olarak bilinir - kesişme noktalarını hızlı bir şekilde belirleyebilecek bazı grafik kütüphaneleri olmalıdır.

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

Bu grafik, mevcut X 50x50 ızgara ve Y eksenleri: