Voxelize a watertight mesh as a Numpy NxNxN array. All the mesh interior is filled, no spatial acceleration structures, as vectorized as I could make it. Works only if all the triangles are within the cube 0x0x0 => NxNxN

voxelize.py

def get_z_line_triangles_intersections_generator(triangles):
    # Precompute everything that is independent from the ray origin
    Z = numpy.array([0., 0., 1.])

    tri_AB = triangles[:, 1] - triangles[:, 0]
    tri_AC = triangles[:, 2] - triangles[:, 0]
    tri_H = numpy.cross(Z, tri_AC)
    tri_a = numpy.sum(tri_AB * tri_H, axis = 1) # dot product along axis 0
            
    # Remove all triangles on a plane parallel to the ray
    is_invisible = numpy.isclose(tri_a, 0)

    tri_A = numpy.delete(triangles[:, 0], is_invisible, axis = 0)    
    tri_AB = numpy.delete(tri_AB, is_invisible, axis = 0)
    tri_AC = numpy.delete(tri_AC, is_invisible, axis = 0)
    tri_H  = numpy.delete(tri_H, is_invisible, axis = 0)
    tri_a  = numpy.delete(tri_a, is_invisible, axis = 0)
    
    # Intersection routine, accelerated with Numba
    def get_intersections(x, y):
        O = numpy.array([x, y, 0.], dtype = 'float64')

        tri_S = O - tri_A
        tri_u = numpy.sum(tri_S * tri_H, axis = 1) / tri_a
        tri_Q = numpy.cross(tri_S, tri_AB)
        tri_v = tri_Q[:,2] / tri_a
        tri_d = numpy.sum(tri_AC * tri_Q, axis = 1) / tri_a
        
        return tri_d[numpy.logical_and(numpy.logical_and(tri_u >= 0, tri_u <= 1.), numpy.logical_and(tri_v >= 0., tri_v + tri_u <= 1.))]
        
    return get_intersections

def voxelize_mesh(triangles, N):    
    ret = numpy.zeros((N, N, N), dtype = bool)

    get_intersections = get_z_line_triangles_intersections_generator(triangles)
    
    # Ray casting loop
    for i in range(N):
        for j in range(N):
            # Get all intersections all intersections along Z axis
            z_list = get_intersections(i + .5, j + .5)
            if z_list.shape[0] < 2:
                continue
                
            z_list.sort()
            z_list = numpy.floor(z_list)
            
            # Rasterize along the Z axis           
            for z0, z1 in zip(z_list[:-1:2], z_list[1::2]):                    
                # Skip segments out of bound
                z0, z1 = int(z0), int(z1)
                if z1 < 0:
                    continue

                # Stop if all the volume have been traversed
                if z0 > N - 1:
                    break
                    
                # Rasterize the segment
                ret[i, j, max(0, z0):min(N, z1 + 1)] = True
                                
    # Job done
    return ret