CODEWITHSUNDEEP
pythonnumpy
Matt
Matt

Reputation: 89

Intersection of line with rectangular prism (Python)

I am trying to calculate the intersection points of a line and a rectangular prism. The volume of the rectangular prism is represented by 6 planes.

s1 = np.array([[-0.25,0,0], [1,0,0]]); s2 = np.array([[0.25,0,0], [1,0,0]]);
s3 = np.array([[0,-0.15,0],[0,1,0]]); s4 = np.array([[0,0.15,0],[0,1,0]]); 
s5 = np.array([[0,0,0],[0,0,1]]); s6 = np.array([[0,0,0.3],[0,0,1]])
surfaces = np.array([s1,s2,s3,s4,s5,s6])

The line is represented by a point and direction.

point, direction = generateRandomUnitVelocity()
p = point; d = direction

To determine whether or not the line intersects a single plane, I run:

def intersectLinePlane(p0, u, p1, n):
   w = p0 - p1
   s = -1 * np.dot(n, w) / np.dot(n, u) # distance from p0 to intersection
   intersectionPoint = p0 + s*u
   return intersectionPoint

To determine where the line intercepts each plane of the volume, I iterate over each plane representing the volume. The line will intercept a side of the volume if it intercepts the plane within a range defined by the dimensions of the volume.

def volumeIntersectionPoints(p0,u,planes):
points = []
for plane in planes:
    point = intersectLinePlane(p0,u,plane[0],plane[1])
    #Define the dimensions of the rectangular prism 
    if ((abs(point[0]) <= 0.25) and (abs(point[1]) <= 0.15) and (0 <= point[2] <= 0.3)):
        points.append(point)
if points:
    return points

As far as I know, volumeIntersectionPoints should return two interception points if the line intercepts the volume as the line enters and exits the volume. This does not necessarily happen. Sometimes I get two interception points as expected, more frequently I'll get one.

#Run this 10000 to get some intersecting lines
for i in range (10000):
   point, direction = generateRandomUnitVelocity()
   p = point; d = direction
   points = volumeIntersectionPoints(p,d,surfaces)
   if points:
        print("Line direction:", d, "Intersection points:", points)
        print("Number of intersection points:", len(points))

Output:

Line direction: [ 0.39894237 -0.60698745 -0.68732178] Intersection points: [array([0.25      , 0.12644022, 0.23721007])]
Number of intersection points: 1
Line direction: [ 0.30552862 -0.72084221 -0.62212441] Intersection points: [array([-0.25      ,  0.1291224 ,  0.21573566]), array([-0.14405106, -0.12084588,  0.        ])]
Number of intersection points: 2
Line direction: [ 0.22272246 -0.14895558 -0.96343497] Intersection points: [array([ 0.13732939, -0.00621236,  0.        ])]
Number of intersection points: 1
Line direction: [ 0.00423704 -0.2063055  -0.97847846] Intersection points: [array([-0.18776991, -0.15      ,  0.0180488 ])]
Number of intersection points: 1
Line direction: [-0.52795903 -0.54028763 -0.65524693] Intersection points: [array([-0.20020341,  0.00703934,  0.        ])]
Number of intersection points: 1

Where am I going wrong here?

Upvotes: 0

Views: 440

Answers (1)

Tim Roberts
Tim Roberts

Reputation: 54733

You are being fooled by rounding errors. When the line does intersect the plane at x=0.25, often the point you compute is actually 0.2500001, and your test of <= 0.25 fails.

I added an epsilon = 0.000001, then widened the test to allow for that slop, and now every line gets 2 intersections.

Upvotes: 1

Related Questions