Creating voxels in Maya, or how to work with MFnMesh::allIntersections. Part 1

Hello everyone,
I would like to share with you some interesting stuff.

Let's fill an object with voxels. In order to do it we need to determine if voxel point is inside the mesh. And it's true, if a ray started at this point in any direction, intersects mesh triangles and number of intersections are even.

Therefore, to fill up the mesh with voxels we need to generate a 3D grid based on object's bounding box, and check each grid point, how many times it intersects the mesh. A key Maya API command in this task is allIntersections  in the class MFnMesh.
Okey, let's write some code. I will comment all important parts right in the script.

import maya.cmds as cmds
import maya.OpenMaya as om
import time

DENSITY = 10        # density of voxels
BORDER = 0.01       # shift border voxels inside a mesh for this amount,
                    # to avoid misses of border voxels due to math on floats

# Float range generator
def frange(start, stop, step):
    r = start
    while not r > stop:
        yield r
        r += step


def main():
    if not cmds.ls(sl = True):
        cmds.warning('Select mesh to generate voxels.')
        return

    # Get MFnMesh function set for selected mesh
    obj = cmds.ls(sl = True)[0]
    selection = om.MSelectionList()
    selection.add(obj)
    obj_mdagpath = om.MDagPath()
    selection.getDagPath(0, obj_mdagpath)
    obj_mfnMesh = om.MFnMesh(obj_mdagpath)

    # Get bounding box for selected mesh and calculate grid step
    bbox = cmds.xform (obj, q = True, ws = True, boundingBox = True)
    grp = cmds.group(em = True, name = 'grp')
    xlen = bbox[3] - bbox[0]
    ylen = bbox[4] - bbox[1]
    zlen = bbox[5] - bbox[2]
    step = round(((bbox[3] - bbox[0]) - 2* xlen * BORDER) /(DENSITY -1), 3)
    print 'step:', step

    # Init all required flags of MFnMesh.allIntersections()
    rayDirection =  om.MFloatVector(bbox[3]*2, bbox[4]*2, bbox[5]*2) # outer ray point
    faceIds = None
    triIds =None
    idsSorted = False
    space = om.MSpace.kWorld
    maxParam = 999999
    testBothDirections = False
    accelParams = None
    sortHits = True
    hitPoints =  om.MFloatPointArray()
    hitRayParams = om.MFloatArray()
    hitFaces = om.MIntArray()

    voxelCount = 0

    # it's not necessary but we want to know how many time computation takes.
    start_time = time.time()

    # Iterate through 3D grid based on mesh bounding box
    for x in frange(bbox[0] + xlen*BORDER, bbox[3], step):
        for y in frange(bbox[1]+ ylen*BORDER, bbox[4], step):
            for z in frange(bbox[2] + zlen*BORDER, bbox[5], step):
                raySource = om.MFloatPoint(x, y, z)
                # Don't forget to clear hitPoints variable or get unpredictable result
                hitPoints.clear()
                obj_mfnMesh.allIntersections(raySource,
                                             rayDirection,
                                             faceIds,
                                             triIds,
                                             idsSorted,
                                             space,
                                             maxParam,
                                             testBothDirections,
                                             accelParams,
                                             sortHits,
                                             hitPoints,
                                             hitRayParams,
                                             hitFaces,
                                             None,None,None,
                                             0.000001)

                # if number of intersections is odd, point of grid is inside mesh volume
                if hitPoints.length()%2 ==1:
                    loc = cmds.spaceLocator()[0]
                    cmds.setAttr(loc + '.scale',0.05,0.05,0.05)
                    cmds.xform(loc, ws = True, t = [x,y,z])
                    cmds.parent(loc, grp)
                    voxelCount +=1
    print 'Number of created voxel primitives:', voxelCount
    print 'Script completed in ', round(time.time() - start_time, 4), "seconds"

How to use it: Select any mesh object, and run the script.
Remark: Nothing special. Only one thing confused me, and it took time to understand what was wrong. In each loop of the cycle before usage in the mfnMesh.allIntersections method, we need to clear array hitPoints, or it has wrong values.
In Part2 we are going to create Voxel node using C# and Maya .NET API.

Thank you. See you soon in Part 2.