How to find all the intersection points between two contour sets in an effective way

Question

I wonder how best to find all the intersection points (before the rounding error) between the two sets of contour lines. What is the best method? Here is an example:

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

I need a little bit like:

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

gives

Answer 2

Work with @unutbu's answer and the intersection algorithm taken straight from numpy-and-line-intersections . I came up with this. It is slow due to the brute force method of finding intersections and loops inside loops inside loops. There may be a way to make cycles faster, but I'm not sure about the intersection algorithm.

 import matplotlib.pyplot as plt import numpy as np def find_intersections(A, B): #this function stolen from /questions/41518/numpy-and-line-intersections # 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() 

Edit: find_intersections above does not handle vertical and horizontal lines or overlapping line segments. Below is the linelineintersect function that handles these cases, but the whole process is still slow. I added an account, so you have an idea of ​​how long it should go. I also changed contour1 = plt.contour(Z1,colors='k') to contour1 = plt.contour(X,Y,Z1,colors='k') , so the axes and your intersections are in the actual X and Y coordinates of your mesh, not counting your mesh points.

I think the problem is that you just have a lot of data. One contour set has 24 lines with a total of 12818 segments, another set of contours has 29 lines with 11411 line segments. That there are many combinations of line segments to check (696 calls on linelineintersect ). You can get the result faster by using a coarser grid to calculate your functions (100 to 100 was much faster than your original 500 to 500). You may be able to execute some kind of linear sweep algorithm, but I don’t know how to implement such things. The line crossing problem has many applications in computer games and is known as collision detection - there must be some kind of graphics library that can quickly identify all the points of intersection.

 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(aT) 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(aT) 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 #eg 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() 

This graph is for a 50x50 grid with actual X and Y axes: