Determine if 2 segments located in 2 matrices intersect

Question

I have two matrices stored as numpy arrays, each containing a number of segments. One segment per row. They look as such:

ts = np.array([[ 26.,  13.,  40.,   0.],
               [ 67.,  17.,  48.,  79.],
               [  8., 102.,  11., 117.]])
vs = np.array([[ 237.01219512,    7.52439024,  158.50929615,  119.50398406],
               [ 237.01219512,    7.52439024, -113.48853868,  481.39398281],
               [ 158.50929615,  119.50398406,  154.50068752,  102.91663802]])

Basically each storing the X1, Y1, X2, Y2 of the segments. Each matrix can store thousands of segments. I need to find which segments in vs intersect at least one segment in ts (and in theory they should only intersect one of them or none), flag it and remove them from vs.

For now, I have the following code to do so:

independent = np.ones_like(np.arange(vs.shape[0]), bool)
for index in range(len(vs)):
    v = vs[index]
    for t in ts:
        if intersect([t[0], t[1]], [t[2], t[3]], [v[0], v[1]], [v[2], v[3]]):
            independent[index] = False
            break

output = vs[independent]

The code for the intersect function is basically this one, marginally modified:

def intersect(a, b, c, d):
    def ccw(a, b, c):
        ax = a[0]
        ay = a[1]
        bx = b[0]
        by = b[1]
        cx = c[0]
        cy = c[1]
        return (cy - ay) * (bx - ax) > (by - ay) * (cx - ax)

    return ccw(a, c, d) != ccw(b, c, d) and ccw(a, b, c) != ccw(a, b, d)

My loop is pretty naive but doing the job. However, as the number of items in the matrices gets big, it's getting obviously pretty slow. I guess there should be a way to do this in a vectorized way without having to loop through all the elements but I can't see how.

My idea would be to call something like intersect(ts, vs) and get a list of segments that intersect, or their indices in the matrix. But I can't find a way to vectorize the function intersect. How can I do that? Is it even possible?

If it's of any help, ts corresponds to the edges of the segments that are part of the minimum spanning tree of a Delaunay triangulation and vs corresponds to the edges of the polygons making up the dual Voronoi diagram. So basically, I'm trying to identify the edges of the Voronoi diagram that intersect the minimum spanning tree of the Delaunay triangulation. My triangulation and diagram are computed using scipy.spatial.qhull.Delaunay and Voronoi respectively.

Thank you.

EDIT:

I came up with an updated version of my code. It feels faster but doesn't look very pretty and I don't think it's done a very pythonic way:

def intersect2(t, v):
    t0 = np.vstack(t[:, 0])
    t1 = np.vstack(t[:, 1])
    t2 = np.vstack(t[:, 2])
    t3 = np.vstack(t[:, 3])
    v = v.transpose()
    v0 = v[0]
    v1 = v[1]
    v2 = v[2]
    v3 = v[3]

    def ccw(ax, ay, bx, by, cx, cy):
        return (cy - ay) * (bx - ax) > (by - ay) * (cx - ax)

    return np.logical_and(ccw(t0, t1, v0, v1, v2, v3) != ccw(t2, t3, v0, v1, v2, v3),
                          ccw(t0, t1, t2, t3, v0, v1) != ccw(t0, t1, t2, t3, v2, v3))


def intersection_mask(t, v):
    inter = intersect2(t, v)
    those_that_intersect = np.array(np.where(inter))[1]
    mask = np.ones_like(np.arange(v.shape[0]), bool)
    mask[those_that_intersect] = False
    return v[mask]

I can now call intersection_mask(t, v) and it will flag those edges of the Voronoi diagram intersecting the minimum spanning tree and remove them. But again, I don't think it's really pythonic and although it seems faster than the first version above, it already starts to slow down for a few hundreds of segments per matrix.

Determine if 2 segments located in 2 matrices intersect

EDIT:

0 Answers0