How to detect self-intersections in a spatial path in Python?

Question

Let's say I have a path in a multidimensional space, given by an array of points along the path. E.g.

p[0] = [  0,  0,  0,   0,...]
p[1] = [  0,  0,  0,0.01,...]

etc.

What's an efficient algorithm for detecting if the path has any self-intersections? It's entirely possible that the path intersects at a point in between the points I have stored.

For example, in 2D :

p[0] = [  0,   0]
p[1] = [  1,   0]
p[2] = [  1,   1]
p[3] = [0.5,-0.5]

This path doesn't have any identical points, but has a self-intersection .

Edit:

I'm currently testing if any points fall within cylinders created at each pair of points.

Here's the code I've been playing around with:

import numpy as np

def pairs(l):
    for i in range(0,len(l),2):
        yield (i,l[i:i+2])

def in_cyl(h0,h1,tol,p):
    # efficient cylinder test from https://www.flipcode.com/archives/Fast_Point-In-Cylinder_Test.shtml
    l = np.linalg.norm(h0-h1)
    dx = h1-h0
    for point in p:
        pdx = point-h0
        dot = np.dot(pdx,dx)
        if dot < 0 or dot > l**2:
            # outside end caps
            continue
        else:
            # point lies within end caps. Find dist from point to cyl axis
            dsq = np.dot(pdx,pdx) - dot**2/l**2
            if (dsq > tol**2):
                # outside cyl
                continue
            else:
                # inside cyl
                return True
    return False

def self_intersect(p,tol):
    for i,(h0,h1) in pairs(p):
        if in_cyl(h0,h1,tol,np.vstack([p[:i],p[i+2:]])):
            return True
    return False

# 50-dimensional test. Intersections should be very rare
dim = 50
test_points = np.random.randn(2500,(50))
print(self_intersect(test_points,0.1))

On my machine, this is fairly slow.

%timeit self_intersect(np.random.randn(2000,(10)),0.1)
10 s ± 94.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Answer 1

Self-intersection search can be implemented via BVH traverse with complexity of O(log(n)) + O(n*log(n)) to build the tree (usually BVH is built once and used many times for different purposes).

For this case traverse is quite simple, let me show example for binary BVH:

1. First we define task - two nodes of BVH that intersects
2. Start with task(root node, root node)
3. If task.nodeA = task.nodeB:
   • add Task(task.nodeA.left,task.nodeA.left)
   • add Task(task.nodeA.right,task.nodeA.right)
   • add Task(task.nodeA.left,task.nodeA.right)
4. Otherwise if task.nodeA intersects task.nodeB:
   • if task.nodeA and task.nodeB are leaves - check intersection (be careful with lines that share point, they should not be counted)
   • otherwise split the node with less volume (if it is not leaf, then split other one) and add Task(task.smallNode,task.splitedNode.left) and Task(task.smallNode,task.splitedNode.right)

You can find c++ implementation for 2d polyline in opensource library MeshLib (it uses AABB tree BVH structure inside). It also has python modules where you can run following code:

from meshlib import mrmeshpy
from meshlib import mrmeshnumpy
import numpy as np

# mrmeshpy uses float32 for vertex coordinates
# however, you could also use float64
points = np.ndarray(shape=(4,2), dtype=np.float32, buffer=np.array([[0.0,0.0],[1.0,0.0],[1.0,1.0],[0.5,-0.5]], dtype=np.float32))

# create polyline structure from numpy points array
polyline = mrmeshnumpy.polyline2FromPoints(points)

# find indices of self-intersecting polyline edges
selfInters = mrmeshpy.findSelfCollidingEdges(polyline)
assert (selfInters.size() == 1)

# get line segments of that edges
segm1 = polyline.edgeSegment(mrmeshpy.EdgeId(selfInters[0].aUndirEdge))
segm2 = polyline.edgeSegment(mrmeshpy.EdgeId(selfInters[0].bUndirEdge))

# find coordinates of intersection
inter = mrmeshpy.intersection(segm1,segm2) # inter == [0.666667,0.0]
assert (abs(inter.x - 2.0/3.0) < 0.001 and inter.y == 0)