3D line segment box intersection

Question

Looking for code to detect an intersection between a 3D segment (not a line/ray) and a 3D box (not necessarily a cube, but always axis-aligned). The boxes are voxels so they have regular spacing.

Already have the code to find a segment/plane intersection. Ideally, I'd like to find an efficient solution to adapt this for a rectangle, repeat for each face of the 3d box, and then iterate for tens of thousands of segments and boxes.

seg_start = array([x1,y1,z1])
seg_end = array([x2,y2,z2])
plane_point = array([x3,y3,z3])
plane_normal = array([x4,y4,z4])
u = seg_end - seg_start
w = seg_start - plane_point
D = dot(plane_normal,u)
N = -dot(plane_normal,w)
sI = N / D
if sI >= 0 and sI <= 1:
    return 1

Answer 1

Just to have this archived with a sufficient background (the answers above seem to not include some collision cases):

This question looks to be quite the same, with a good answer.

In this thesis the problem is explained with pictures and descriptions on pages 15-19.

link went bad over time

Answer 2

Because the box is axis-aligned, all you need to do is to check for interval intersection in each coordinate.

Here is an example in python, with some tests. Note that it is generic for N dimensions, and it is the same algorithm for box-box intersection:

def are_intervals_intersecting(a0, a1, b0, b1):
    '''
    @param a0: float
    @param a1: float
    @param b0: float
    @param b1: float
    '''
    if (a1 < a0):
        a1, a0 = a0, a1

    if (b1 < b0):
        b1, b0 = b0, b1

    # 6 conditions:

    # 1)
    #        a0 ---------- a1                              a0 < b0 and a1 < b0
    #                             b0 ---------- b1         (no intersection)

    # 2)
    #               a0 ---------- a1
    #                      b0 ---------- b1                (intersection)

    # 3)
    #               a0 ------------------------ a1
    #                      b0 ---------- b1                (intersection)

    # 4)
    #                      a0 ---------- a1         
    #               b0 ------------------------ b1         (intersection)

    # 5)
    #                             a0 ---------- a1         (intersection)
    #                      b0 ---------- b1

    # 6)
    #                                    a0 ---------- a1  b0 < a0 and b1 < a0         
    #               b0 ---------- b1                       (no intersection)

    if b0 < a0:
        # conditions 4, 5 and 6
        return a0 < b1 # conditions 4 and 5
    else:
        # conditions 1, 2 and 3
        return b0 < a1 # conditions 2 and 3


def is_segment_intersecting_box(P0, P1, B0, B1):
    '''
    @param P0: tuple(float)
    @param P1: tuple(float)
    @param B0: tuple(float)
    @param B1: tuple(float)
    '''
    for i in xrange(len(P0)):
        if not are_intervals_intersecting(P0[i], P1[i], B0[i], B1[i]):
            return False
    return True


if __name__ == '__main__':
    assert not is_segment_intersecting_box(
        (0.0, 0.0, 0.0), (1.0, 1.0, 1.0), (2.0, 2.0, 2.0), (3.0, 3.0, 3.0))

    assert not is_segment_intersecting_box(
        (0.0, 0.0, 0.0), (4.0, 1.0, 1.0), (2.0, 2.0, 2.0), (3.0, 3.0, 3.0))

    assert not is_segment_intersecting_box(
        (1.5, 1.5, 0.0), (4.0, 2.5, 1.0), (2.0, 2.0, 2.0), (3.0, 3.0, 3.0))

    assert is_segment_intersecting_box(
        (1.5, 1.5, 0.0), (4.0, 2.5, 2.5), (2.0, 2.0, 2.0), (3.0, 3.0, 3.0))

    assert is_segment_intersecting_box(
        (1.5, 1.5, 1.5), (2.5, 2.5, 2.5), (2.0, 2.0, 2.0), (3.0, 3.0, 3.0))

    assert is_segment_intersecting_box(
        (2.5, 2.5, 2.5), (2.6, 2.6, 2.6), (2.0, 2.0, 2.0), (3.0, 3.0, 3.0))

    assert is_segment_intersecting_box(
        (2.5, 2.5), (2.5, 3.5), (2.0, 2.0), (3.0, 3.0))

    print 'ok'

Answer 3

First off, you probably meant and rather than or in your if-condition, otherwise it'll just always return true. Second, if you're just testing whether there is an intersection or not, you can do it faster (without the float-division):

• You can determine which "side" of each plane any given point is on using vector math:
  side = dot(my_point - plane_point, plane_normal)
  Now if side is positive, my_point is "in front of" the plane (i.e. it's on the side the normal is pointing towards); if negative, it's "behind" the plane. If side is zero, your point lies on the plane.

• You can check whether your segment intersects an (infinite) plane by just testing to see if the start point and end point are on different sides:

  start_side = dot(seg_start - plane_point, plane_normal)
  end_side = dot(seg_end - plane_point, plane_normal)
  return start_side * end_side
  #if < 0, both points lie on different sides, hence intersection
  #if = 0, at least one point lies on the plane
  #if > 0, both points lie on the same side, i.e. no intersection
  

• You can use the "side" check to do the axis-aligned-cuboid intersection too (actually, this will work for any parallelpiped):

  • Treat your box as a set of six planes
  • Make sure the plane normals are all pointing either "outwards" or "inwards" from the box. I'll assume you're using "outwards"
  • For any point to lie inside your box, it has to be "behind" all six planes. If it isn't, it lies outside the box.
  • For any segment to intersect the box, one point has to lie outside it and one inside.
  • That's all!

edit: The last point is actually incorrect; as you say, voxels can be intersected even if both endpoints lie outside. So it's not the whole solution - actually, you can't really do this without calculating the intersection point. You can, however, still use the "side test" as an early-reject mechanism, so as to cut down on the number of full calculations you need to do: if both points are on the same side of any one of the six planes, there can be no intersection.

As far as your specific case goes, it seems like you're trying to find all intersecting voxels for some given line segment? In that case, you'd probably be better served using something like Bresenham's to explicitly calculate the path, instead of testing for intersections against all of the voxels...