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