python algorithm geometry categorization

Find if segment comes within distance of another one

I have a bunch of segments (the data I have are the 2 points that makes the segment [x1,y1] and [x2, y2]) and would like to categorize them according to their position. If a segment is close enough to another one, then I want to put them together. If I had to describe it in sentence: I would like to find all neighbor segments with a distance of 5px away from any point of the segment.

With rules similar to the drawn picture:

I was looking at different algorithms but most of them deal with intersections and predict whether the lines will intersect. That does not really work for me as I don't want to continue the lines to infinity, I just want to know if they come within 5px of each other.

Does anyone knows how I can approach this problem (and relatively fast)? Do you know any frameworks that could help? (I was looking at the nearest neighbors but I cannot find any framework that deals with geometry instead of points).

Thanks!

Solution

(I have revised my previous answer. That answer hat some shortcomings. I think my new answer shows a simpler and more robust solution.)

You have two segments, S with points S₀ and S₁ and T with poins T₀ and T₁. A collision is detected when these segments are less than a distance of r apart at one point.

For the segment S, you get the direction vector Δs, the segment length s and the normalized direction vector u.

Δs = S₁ − S₀
s = |Δs|
u = Δs / s

The unit vector u and the point S₀ can describe a transformation of any point P to a point P′:

P′_x = (P_x − S_0x) · u_x + (P_y − S_0y) · u_y
P′_y = − (P_x − S_0x) · u_y + (P_y − S_0y) · u_x

In this transformation the points of the segment S are:

S′₀ = (0, 0)
S′₁ = (s, 0)

For the transformed points T′₀ and T′₁, the y components can be interpreted as signed distance to S. Now several tests can be performed:

The first test is whether T′₀ or T′₁ are within a distance of r of the segment S or within a radius of r of either S₀′ or S₁′. If so, we have a hit.
The next test is whether the two lines intersect. That can only happen if the signs of T′_0y or T′_1y are different. If so, we have a hit.
For the last test, we reverse the first test by transforming S to S′′ in a system where T is aligned to the x-axis. Then test whether one of the transformed points S′′₀ or S′′₁ are within r of T′′. If so, we have a hit, otherwise it's a miss.

Python code is below. I've also updated my JS Fiddle.

Notes:

The longitudinal variable a and the distance d in my old answer were in effect the same as the x′ and y′ here. I think the transformation is simpler.
This solution tests only (1) whether the points of T are within a distance of r from S, (2) whether the lines intersect and (3) whether the points of S are within a distance of r from T. The case of collinear line segments is caught by the tests (1) and (3).
The code below does not handle zero-length segments (S₀ = S₁ or T₀ = T₁) explicitly, but returning a non-null vector as norm of a null vector seems to do the trick – tests (1) and (3) catch these cases.

Python code:

import math

class Point:
    """ A point P(x, y) in 2D space
    """

    def __init__(self, x, y):
        self.x = float(x)
        self.y = float(y)

class Vector:
    """ A vector v(x, y) in 2D space
    """

    def __init__(self, x, y):
        self.x = x
        self.y = y
    
    def mag(self):
        """ magnitude of the vector
        """
        
        return math.hypot(self.x, self.y)
    
    def norm(self):
        """ return the normalized vector or (0, 0)
        """
    
        a = self.mag()
        
        if a*a < 1.0e-16:
            return Vector(1, 0)
        
        return Vector(self.x / a, self.y / a)
    


def diff(p, q):
    """ difference vector (q - p)
    """

    return Vector(q.x - p.x, q.y - p.y)

def within(p, dx, r):
    """ Is p within r of point (dx, 0)?
    """

    x = p.x - dx
    y = p.y
    
    return x*x + y*y <= r*r

def rot(p, u):
    """ Rotate point p to a coordinate system aligned with u.
    """

    return Point(p.x * u.x + p.y * u.y,
                -p.x * u.y + p.y * u.x)

def collision(s, t, r):
    """ Do the line segments s and t collide with a radius r
    """

    ds = diff(s[0], s[1])
    ss = ds.mag()    
    u = ds.norm()
    
    a0 = rot(diff(s[0], t[0]), u)
    a1 = rot(diff(s[0], t[1]), u)

    # Test T0 and T1 against S
    
    if -r <= a0.y <= r and -r <= a0.x <= ss + r:
        if a0.x < 0: return within(a0, 0, r)
        if a0.x > ss: return within(a0, ss, r)
        return True    
    
    if -r <= a1.y <= r and -r <= a1.x <= ss + r:
        if a1.x < 0: return within(a1, 0, r)
        if a1.x > ss: return within(a1, ss, r)
        return True

    # Test intersection
    
    if a0.y * a1.y < -0.9 * r * r:
        a = -a0.y * (a1.x - a0.x) / (a1.y - a0.y) + a0.x
        
        if 0 <= a <= ss: return True

    # Test S0 and S1 against T
    
    dt = diff(t[0], t[1])
    tt = dt.mag()    
    v = dt.norm()
    
    b0 = rot(diff(t[0], s[0]), v)
    b1 = rot(diff(t[0], s[1]), v)

    if 0 <= b0.x <= tt and -r <= b0.y <= r: return True
    if 0 <= b1.x <= tt and -r <= b1.y <= r: return True
       
    return False