Rounding errors when computing intersection between two lines in python

Question

I'm attempting to create a Python script that identifies all the unique intersection points, eliminating duplicates, among all the lines passing through a subset of points of a regular grid in 2D or 3D cartesian space.

I successfully developed a script using the Sympy library, but its performance is inadequate due to excessive slowness. Subsequently, I came across a library called scikit-spatial that appears to simplify 3D object implementation. I crafted a script using this library, but I encountered precision issues resulting in some points being erroneously identified as distinct when they shouldn't be.

In my quest for a solution, I explored an alternative algorithm without relying on external libraries, adapted from this answer. Unfortunately, I'm facing the same precision issues. Other solutions I came across either proved too challenging to implement (given my limited coding ability) or were designed exclusively for 2D scenarios.

Here are a couple of code snippets that should reproduce the aforementioned issues. For simplicity, in these snippets, I'm only evaluating the unique intersection points between the lines generated from a 3x3 square grid, which should be 61 distinct points.

With skspatial:

import skspatial.objects as sks


def lineIntersection(p1, p2, q1, q2):
    v1, v2 = sks.Vector([p2[0]-p1[0], p2[1]-p1[1], p2[2]-p1[2]]), sks.Vector([q2[0]-q1[0], q2[1]-q1[1], q2[2]-q1[2]]) #direction vectors of the two lines
    if v1.is_parallel(v2): #I want to discard infinite intersection points of overlapping lines
        return ()
    L1 = sks.Line(point=sks.Point([p1[0], p1[1], p1[2]]), direction=v1) #Line defined via a point and its direction vector
    L2 = sks.Line(point=sks.Point([q1[0], q1[1], q1[2]]), direction=v2)
    intersection = tuple(L1.intersect_line(L2)) #intersection point
    return intersection

def grid_fnc(files, rows, cols = 0):
    grid = []
    for file in range(files):
        for row in range(rows):
            if cols != 0:
                for col in range(cols):
                    grid.append((file, row, col))
                continue
            grid.append((file, row, 0.0))
    return grid

def intersectionPoints_fnc(grid):
    intersectionPoints = []
    for p1 in grid: #start point of Line 1 (L1)
        for p2 in grid: #end point of L1
            if p1 == p2: #a line is defined by two distinct points
                continue
            for q1 in grid: #start point of L2
                for q2 in grid: #end point of L2
                    if q1 == q2:
                        continue
                    intersectionPoint = lineIntersection(p1, p2, q1, q2) #evaluating intersection between two lines
                    if intersectionPoint != ():
                        intersectionPoints.append(intersectionPoint)
    intersectionPoints = list(set(intersectionPoints)) #removing duplicate points
    return intersectionPoints

def main():
    grid = grid_fnc(3, 3)
    intersectionPoints = intersectionPoints_fnc(grid)
    print(len(intersectionPoints))
    print(intersectionPoints)


main()

Results (points in random order):

len(intersectionPoints) = 101
intersectionPoints = [(-1.0, 4.0, 0.0), (0.6666666666666666, 1.6666666666666667, 0.0), (1.3333333333333333, 1.6666666666666665, 0.0), (1.5, 1.5, 0.0), (0.8, 0.4, 0.0), (0.3999999999999999, 0.8, 0.0), (1.3333333333333333, 0.6666666666666666, 0.0), (0.5, 1.0, 0.0), (1.0, 1.0, 0.0), (0.6666666666666667, 1.6666666666666665, 0.0), (1.0, 0.5, 0.0), (1.3333333333333335, 1.3333333333333333, 0.0), (1.3333333333333335, 0.33333333333333337, 0.0), (4.0, 2.0, 0.0), (0.0, 0.0, 0.0), (0.3333333333333333, 1.3333333333333333, 0.0), (-1.0, 2.0, 0.0), (0.6666666666666666, 1.6666666666666665, 0.0), (-1.0, -2.0, 0.0), (2.0, 2.0, 0.0), (-2.0, -1.0, 0.0), (-1.0, 0.0, 0.0), (-2.0, 0.0, 0.0), (2.0, 4.0, 0.0), (0.6666666666666667, 1.3333333333333335, 0.0), (0.6666666666666667, 0.3333333333333333, 0.0), (1.3333333333333335, 1.6666666666666665, 0.0), (1.3333333333333335, 0.6666666666666666, 0.0), (0.0, 3.0, 0.0), (1.6666666666666665, 0.6666666666666667, 0.0), (4.0, -1.0, 0.0), (0.6666666666666667, 0.6666666666666667, 0.0), (0.6666666666666666, 1.3333333333333335, 0.0), (4.0, 0.0, 0.0), (0.4, 0.8, 0.0), (1.2, 0.3999999999999999, 0.0), (0.3333333333333333, 0.6666666666666666, 0.0), (1.5, 0.5, 0.0), (2.0, -1.0, 0.0), (2.0, -2.0, 0.0), (3.0, 4.0, 0.0), (1.3333333333333335, 0.6666666666666667, 0.0), (0.6666666666666666, 0.6666666666666667, 0.0), (1.3333333333333333, 0.3333333333333333, 0.0), (2.0, 0.0, 0.0), (0.33333333333333337, 0.6666666666666666, 0.0), (0.6666666666666667, 0.33333333333333337, 0.0), (1.3333333333333333, 1.3333333333333335, 0.0), (1.6666666666666667, 0.6666666666666667, 0.0), (0.3999999999999999, 1.2, 0.0), (1.6666666666666665, 1.3333333333333335, 0.0), (1.3333333333333333, 0.6666666666666667, 0.0), (0.5, 1.5, 0.0), (1.2, 0.4, 0.0), (1.0, 1.5, 0.0), (1.2, 1.6, 0.0), (0.0, 1.0, 0.0), (0.6666666666666666, 0.33333333333333337, 0.0), (3.0, 0.0, 0.0), (1.6666666666666665, 0.6666666666666666, 0.0), (0.8, 1.6, 0.0), (0.33333333333333337, 1.3333333333333335, 0.0), (1.6, 1.2, 0.0), (0.33333333333333337, 0.6666666666666667, 0.0), (1.3333333333333335, 0.3333333333333333, 0.0), (0.6666666666666666, 0.3333333333333333, 0.0), (-2.0, 3.0, 0.0), (1.3333333333333333, 1.6666666666666667, 0.0), (0.8, 0.3999999999999999, 0.0), (1.5, 1.0, 0.0), (3.0, 2.0, 0.0), (1.6666666666666665, 1.3333333333333333, 0.0), (0.0, 2.0, 0.0), (0.6666666666666667, 1.6666666666666667, 0.0), (1.6, 0.8, 0.0), (1.3333333333333335, 1.3333333333333335, 0.0), (1.0, 0.0, 0.0), (0.6666666666666667, 1.3333333333333333, 0.0), (1.6666666666666667, 1.3333333333333335, 0.0), (1.0, 2.0, 0.0), (0.3333333333333333, 1.3333333333333335, 0.0), (4.0, 3.0, 0.0), (0.33333333333333337, 1.3333333333333333, 0.0), (0.6666666666666666, 1.3333333333333333, 0.0), (0.4, 1.2, 0.0), (2.0, 1.0, 0.0), (1.6666666666666667, 0.6666666666666666, 0.0), (0.3333333333333333, 0.6666666666666667, 0.0), (2.0, 3.0, 0.0), (1.3333333333333335, 1.6666666666666667, 0.0), (3.0, -2.0, 0.0), (1.6666666666666667, 1.3333333333333333, 0.0), (0.0, -2.0, 0.0), (0.0, -1.0, 0.0), (0.5, 0.5, 0.0), (0.6666666666666667, 0.6666666666666666, 0.0), (1.3333333333333333, 1.3333333333333333, 0.0), (0.0, 4.0, 0.0), (1.3333333333333333, 0.33333333333333337, 0.0), (0.6666666666666666, 0.6666666666666666, 0.0), (-2.0, 2.0, 0.0)]

As you can see, there are repeated points treated differently due to rounding errors, like (1.6666666666666667, 0.6666666666666667, 0.0) and (1.6666666666666665, 0.6666666666666666, 0.0).

Without external libraries (fastest method by far):

def dot(u, v):
    return u[0] * v[0] + u[1] * v[1] + u[2] * v[2]

def norm2(v):
    return v[0] * v[0] + v[1] * v[1] + v[2] * v[2]

def cross(b, c):
    return [b[1] * c[2] - c[1] * b[2], b[2] * c[0] - c[2] * b[0], b[0] * c[1] - c[0] * b[1]]

def lineIntersection(a, b): #a=L1(p1, p2) b=L2(q1, q2)
    da = [a[1][0] - a[0][0], a[1][1] - a[0][1], a[1][2] - a[0][2]]
    db = [b[1][0] - b[0][0], b[1][1] - b[0][1], b[1][2] - b[0][2]]
    dc = [b[0][0] - a[0][0], b[0][1] - a[0][1], b[0][2] - a[0][2]]

    if dot(dc, cross(da, db)) != 0.0:
        return ()
    
    if norm2(cross(da, db)) == 0.0: #not interested in parallel lines
        return ()
        
    s = dot(cross(dc, db), cross(da, db)) / norm2(cross(da, db))
    
    x = a[0][0] + da[0] * s
    y = a[0][1] + da[1] * s
    z = a[0][2] + da[2] * s
    ip = (x, y, z) #intersection point
    return ip
    
def grid_fnc(files, rows, cols = 0):
    grid = []
    for file in range(files):
        for row in range(rows):
            if cols != 0:
                for col in range(cols):
                    grid.append((file, row, col))
                continue
            grid.append((file, row, 0.0))
    return grid

def intersectionPoints_fnc(grid):
    intersectionPoints = []
    for p1 in grid:
        for p2 in grid:
            if p1 == p2:
                continue
            for q1 in grid:
                for q2 in grid:
                    if q1 == q2:
                        continue
                    intersectionPoint = lineIntersection((p1, p2), (q1, q2))
                    if intersectionPoint != ():
                        intersectionPoints.append(intersectionPoint)
    intersectionPoints = list(set(intersectionPoints))
    return intersectionPoints

def main():
    grid = grid_fnc(3, 3)
    intersectionPoints = intersectionPoints_fnc(grid)
    print(len(intersectionPoints))
    print(intersectionPoints)

main()

Results have the same kind of issue:

len(intersectionPoints) = 101
intersectionPoints = [(-1.0, 4.0, 0.0), (0.6666666666666666, 1.6666666666666667, 0.0), (1.3333333333333333, 1.6666666666666665, 0.0), (1.5, 1.5, 0.0), (0.8, 0.4, 0.0), (0.3999999999999999, 0.8, 0.0), (1.3333333333333333, 0.6666666666666666, 0.0), (0.5, 1.0, 0.0), (1.0, 1.0, 0.0), (0.6666666666666667, 1.6666666666666665, 0.0), (1.0, 0.5, 0.0), (1.3333333333333335, 1.3333333333333333, 0.0), (1.3333333333333335, 0.33333333333333337, 0.0), (4.0, 2.0, 0.0), (0.0, 0.0, 0.0), (0.3333333333333333, 1.3333333333333333, 0.0), (-1.0, 2.0, 0.0), (0.6666666666666666, 1.6666666666666665, 0.0), (-1.0, -2.0, 0.0), (2.0, 2.0, 0.0), (-2.0, -1.0, 0.0), (-1.0, 0.0, 0.0), (-2.0, 0.0, 0.0), (2.0, 4.0, 0.0), (0.6666666666666667, 1.3333333333333335, 0.0), (0.6666666666666667, 0.3333333333333333, 0.0), (1.3333333333333335, 1.6666666666666665, 0.0), (1.3333333333333335, 0.6666666666666666, 0.0), (0.0, 3.0, 0.0), (1.6666666666666665, 0.6666666666666667, 0.0), (4.0, -1.0, 0.0), (0.6666666666666667, 0.6666666666666667, 0.0), (0.6666666666666666, 1.3333333333333335, 0.0), (4.0, 0.0, 0.0), (0.4, 0.8, 0.0), (1.2, 0.3999999999999999, 0.0), (0.3333333333333333, 0.6666666666666666, 0.0), (1.5, 0.5, 0.0), (2.0, -1.0, 0.0), (2.0, -2.0, 0.0), (3.0, 4.0, 0.0), (1.3333333333333335, 0.6666666666666667, 0.0), (0.6666666666666666, 0.6666666666666667, 0.0), (1.3333333333333333, 0.3333333333333333, 0.0), (2.0, 0.0, 0.0), (0.33333333333333337, 0.6666666666666666, 0.0), (0.6666666666666667, 0.33333333333333337, 0.0), (1.3333333333333333, 1.3333333333333335, 0.0), (1.6666666666666667, 0.6666666666666667, 0.0), (0.3999999999999999, 1.2, 0.0), (1.6666666666666665, 1.3333333333333335, 0.0), (1.3333333333333333, 0.6666666666666667, 0.0), (0.5, 1.5, 0.0), (1.2, 0.4, 0.0), (1.0, 1.5, 0.0), (1.2, 1.6, 0.0), (0.0, 1.0, 0.0), (0.6666666666666666, 0.33333333333333337, 0.0), (3.0, 0.0, 0.0), (1.6666666666666665, 0.6666666666666666, 0.0), (0.8, 1.6, 0.0), (0.33333333333333337, 1.3333333333333335, 0.0), (1.6, 1.2, 0.0), (0.33333333333333337, 0.6666666666666667, 0.0), (1.3333333333333335, 0.3333333333333333, 0.0), (0.6666666666666666, 0.3333333333333333, 0.0), (-2.0, 3.0, 0.0), (1.3333333333333333, 1.6666666666666667, 0.0), (0.8, 0.3999999999999999, 0.0), (1.5, 1.0, 0.0), (3.0, 2.0, 0.0), (1.6666666666666665, 1.3333333333333333, 0.0), (0.0, 2.0, 0.0), (0.6666666666666667, 1.6666666666666667, 0.0), (1.6, 0.8, 0.0), (1.3333333333333335, 1.3333333333333335, 0.0), (1.0, 0.0, 0.0), (0.6666666666666667, 1.3333333333333333, 0.0), (1.6666666666666667, 1.3333333333333335, 0.0), (1.0, 2.0, 0.0), (0.3333333333333333, 1.3333333333333335, 0.0), (4.0, 3.0, 0.0), (0.33333333333333337, 1.3333333333333333, 0.0), (0.6666666666666666, 1.3333333333333333, 0.0), (0.4, 1.2, 0.0), (2.0, 1.0, 0.0), (1.6666666666666667, 0.6666666666666666, 0.0), (0.3333333333333333, 0.6666666666666667, 0.0), (2.0, 3.0, 0.0), (1.3333333333333335, 1.6666666666666667, 0.0), (3.0, -2.0, 0.0), (1.6666666666666667, 1.3333333333333333, 0.0), (0.0, -2.0, 0.0), (0.0, -1.0, 0.0), (0.5, 0.5, 0.0), (0.6666666666666667, 0.6666666666666666, 0.0), (1.3333333333333333, 1.3333333333333333, 0.0), (0.0, 4.0, 0.0), (1.3333333333333333, 0.33333333333333337, 0.0), (0.6666666666666666, 0.6666666666666666, 0.0), (-2.0, 2.0, 0.0)]

Expected results (points in random order):

len(intersectionPoints) = 61
intersectionPoints = [(-1.00000000000000, 4.00000000000000, 0), (0.500000000000000, 1.50000000000000, 0), (4.00000000000000, -1.00000000000000, 0), (4.00000000000000, 3.00000000000000, 0), (-2.00000000000000, 0, 0), (1.00000000000000, 1.50000000000000, 0), (1.00000000000000, 0, 0), (2.00000000000000, 1.00000000000000, 0), (0, 1.00000000000000, 0), (0.333333333333333, 1.33333333333333, 0), (3.00000000000000, -2.00000000000000, 0), (0.666666666666667, 1.66666666666667, 0), (-1.00000000000000, -2.00000000000000, 0), (3.00000000000000, 4.00000000000000, 0), (0.666666666666667, 0.666666666666667, 0), (1.33333333333333, 1.66666666666667, 0), (0.500000000000000, 0.500000000000000, 0), (2.00000000000000, 4.00000000000000, 0), (1.33333333333333, 0.666666666666667, 0), (0.400000000000000, 0.800000000000000, 0), (3.00000000000000, 0, 0), (-1.00000000000000, 0, 0), (0, -2.00000000000000, 0), (0, 4.00000000000000, 0), (1.00000000000000, 0.500000000000000, 0), (2.00000000000000, 0, 0), (-2.00000000000000, 2.00000000000000, 0), (2.00000000000000, -2.00000000000000, 0), (0.666666666666667, 1.33333333333333, 0), (0.800000000000000, 1.60000000000000, 0), (0.400000000000000, 1.20000000000000, 0), (0.333333333333333, 0.666666666666667, 0), (1.60000000000000, 0.800000000000000, 0), (1.66666666666667, 1.33333333333333, 0), (-2.00000000000000, -1.00000000000000, 0), (-2.00000000000000, 3.00000000000000, 0), (0, 0, 0), (1.00000000000000, 2.00000000000000, 0), (0.800000000000000, 0.400000000000000, 0), (-1.00000000000000, 2.00000000000000, 0), (1.50000000000000, 1.00000000000000, 0), (1.33333333333333, 0.333333333333333, 0), (3.00000000000000, 2.00000000000000, 0), (1.60000000000000, 1.20000000000000, 0), (4.00000000000000, 0, 0), (1.66666666666667, 0.666666666666667, 0), (0.666666666666667, 0.333333333333333, 0), (1.50000000000000, 1.50000000000000, 0), (1.20000000000000, 0.400000000000000, 0), (2.00000000000000, 2.00000000000000, 0), (2.00000000000000, -1.00000000000000, 0), (0, 2.00000000000000, 0), (0.500000000000000, 1.00000000000000, 0), (1.33333333333333, 1.33333333333333, 0), (4.00000000000000, 2.00000000000000, 0), (2.00000000000000, 3.00000000000000, 0), (0, -1.00000000000000, 0), (0, 3.00000000000000, 0), (1.20000000000000, 1.60000000000000, 0), (1.50000000000000, 0.500000000000000, 0), (1.00000000000000, 1.00000000000000, 0)]

How can I address these issues? Can you suggest an alternative algorithm or approach that would achieve the desired result without compromising too much on performance (Maybe symmetries could be exploited in some ways to cut down on computation time)? Thank you!

Can you please provide the results you get and expected results. — Yuri Ginsburg
– Yuri Ginsburg, Commented Nov 26, 2023 at 3:47

aerobiomat · Accepted Answer · 2023-11-27 13:33:41Z

1

You should consider using Numpy to represent points and to operate with them.

The problem with floating point calculations is that supposedly equivalent results may vary slightly. So the only possibility is to round them or to use a tolerance when comparing floating their values.

This is your exact algorithm "numpyfied". The only important change is that points are rounded before finding unique ones.

DECIMALS = 6            # Expected precision
EPS = 10**-DECIMALS

def line_intersection(a, b):  # a=L1(p1, p2) b=L2(q1, q2)
    da = a[1] - a[0]
    db = b[1] - b[0]
    dc = b[0] - a[0]

    x = np.cross(da, db)
    if np.abs(np.dot(dc, x)) > EPS:
        return None

    x2 = np.inner(x, x)
    if np.abs(x2) < EPS:  # not interested in parallel lines
        return None

    s = np.dot(np.cross(dc, db), x) / x2
    ip = a[0] + da * s
    return ip

def grid_fnc(files, rows, cols=0):
    if cols == 0:
        cols = 1
    return np.array(np.meshgrid(np.arange(files),
                                np.arange(rows),
                                np.arange(cols))).T.reshape(-1, 3)

def intersectionPoints_fnc(grid):
    intersectionPoints = []
    for i, p1 in enumerate(grid):
        for p2 in grid[i+1:]:
            for j, q1 in enumerate(grid):
                for q2 in grid[j+1:]:
                    intersectionPoint = line_intersection((p1, p2), (q1, q2))
                    if intersectionPoint is not None:
                        intersectionPoints.append(intersectionPoint)
    intersectionPoints = np.unique(np.round(intersectionPoints, decimals=DECIMALS), axis=0)
    return intersectionPoints

grid = grid_fnc(3, 3)
print(grid)
intersectionPoints = intersectionPoints_fnc(grid)
print(len(intersectionPoints))
print(intersectionPoints)

The intersections are compatible with your expected results.

However, there remains an important problem: Loops, and especially nested loops tend to be very slow in Python.

Update: vectorized version

from numpy.core.umath_tests import inner1d

DECIMALS = 6            # Expected precision

def line_intersection(a, b):  # a=L1(p1, p2) b=L2(q1, q2)
    da = a[1] - a[0]
    db = b[1] - b[0]
    dc = b[0] - a[0]

    x = np.cross(da, db)
    x2 = inner1d(x, x)
    s = inner1d(np.cross(dc, db), x) / x2
    ip = (a[0] + da * s[..., None]).reshape(-1, 3)
    valid = np.isfinite(ip).any(axis=-1)
    return ip[valid]

def grid(files, rows, cols=0):
    if cols == 0:
        cols = 1
    return np.array(np.meshgrid(np.arange(files),
                                np.arange(rows),
                                np.arange(cols))).T.reshape(-1, 3)

def intersection_points(grid):
    i1, i2 = np.triu_indices(len(grid), k=1)
    points = line_intersection((grid[i1], grid[i2]), (grid[i1, None], grid[i2, None]))
    return np.unique(np.round(points, decimals=DECIMALS), axis=0)

grid = grid(3, 3)
with np.errstate(all='ignore'):
    intersectionPoints = intersection_points(grid)
print(len(intersectionPoints))
print(intersectionPoints)

edited Nov 27, 2023 at 13:33

answered Nov 27, 2023 at 8:55

aerobiomat

3,4371 gold badge18 silver badges22 bronze badges

Sign up to request clarification or add additional context in comments.

4 Comments

Edoardo Serra Over a year ago

Thank you! I'll try them this evening. I thought of using the fractions module that in this case would give accurate answers since all the lines have rational coefficients, but I guess it wouldn't work in more general cases with irrational coefficients involved. I'm still trying to properly understand how floating-point arithmetics works.

Edoardo Serra Over a year ago

Regarding the vectorized version, it's WAY faster than the other but I think there might be a problem somewhere. I've changed inner1d with np.einsum('...i,...i', x, x), as inner1d is deprecated, and still got the correct result for a 3x3 grid. However, I don't get correct results for bigger grids, like the 4x4x4 one, in which I get a len of 688315 whereas it should be of "only" 13307.

aerobiomat Over a year ago

There is an error in valid = np.isfinite(ip).any(axis=-1). It should be all() instead of any(). This may be the problem or not. I'll check it if I find some time.

Edoardo Serra Over a year ago

I've tried it but still got the wrong result. My guess is there is something that goes wrong when the column points are not all 0s.

Collectives™ on Stack Overflow

Rounding errors when computing intersection between two lines in python

1 Answer 1

4 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

4 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related