7

I'm wondering about the best way to find all the intersection points (to roundoff error) between two sets of contour lines. Which is the best method? Here is the 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()

enter image description here

I want some similar to:

intersection_points = intersect(contour1,contour2)
print intersection_points
[(x1,y1),...,(xn,yn)]
Jason Aller
  • 3,541
  • 28
  • 38
  • 38
Pablo
  • 2,443
  • 1
  • 20
  • 32
  • Are you looking for the intersection of say the z=0.75 contour from Z1 with the z=0.75 contour from Z2? Or do you want all the intersections of all the lines in your plot? – rtrwalker Jul 02 '13 at 03:08
  • I'm looking for the intersection of all common points between all the lines on the plot. All contours with all contours. Actually, for example, if you solve the problem for one red and all blacks, the problem is solved... – Pablo Jul 02 '13 at 03:13

2 Answers2

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

yields

enter image description here

unutbu
  • 842,883
  • 184
  • 1,785
  • 1,677
  • good aproximation! but in the middle you have a lot of wrong intersection points... My question is, intersection to "roundoff error", with good performance calculation. – Pablo Jul 02 '13 at 03:17
  • @Pablo: I've added some clustering code to pick just one point for each intersection. – unutbu Jul 02 '13 at 14:11
  • Nice!! With this method you're computing distances between points, but they are not true intersections... I think this method is prone to fail with a more sophisticated example. I'm feel that it is not the best way to attack this problem. I don't know, I'm apologize if I'm not understand something. I will think about your solution... Maybe your method could be developed in a more general way. Anycase, is a very good answer!! It's fast and give a good portrait of intersections... I will wait some days before the definitive answer. Thank you very much! Really, your method is amazing! – Pablo Jul 02 '13 at 15:14
  • Yes, there might be a better way. I look for to learning from it when the answer comes. – unutbu Jul 02 '13 at 16:16
4

Working from the answer of @unutbu and an intersection algorithm plucked straight from numpy-and-line-intersections I've come up with this. It is slow owing to the brute-force sort of way of finding the intersections and the loops within loops within loops. There might be a way to make the loops 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 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()

enter image description here

Edit: find_intersections above does not handle vertical and horizontal lines nor overlapping line segments. The linelineintersect function below does handle those cases but the whole process is still slow. I have added a count so you have some idea of how long ther is to 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 coords of your grid rather than a count of your grid points.

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

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

This plot is for 50x50 grid with the actual X and Y axes: enter image description here

Community
  • 1
  • 1
rtrwalker
  • 1,021
  • 6
  • 13
  • I think the intersection algorithm breaks down for vertical line segments and horizontal line segments – rtrwalker Jul 02 '13 at 07:33
  • Amazing result!!! But it's too slow... Anycase, it give the solution! I don't want to mark it as the answer because the performance is poor, but obviously, it's a possible answer! I will think about it... – Pablo Jul 02 '13 at 13:13