In Python, how do I voxelize a 3D mesh

Question

I need help on getting started in Python (which I almost know nothing of) to voxelize a 3D mesh generated from Rhino. The data input will be a .OBJ file and so will the output. The ultimate purpose for this usage is to find the shortest distance between two points within a building. But that's for later. As for now, I need to first voxelize a 3D mesh. The voxelization primitive may just be a simple cube.

so far I can read from a OBJ file parser and out of the parsed obj with V, VT, VN, F prefixes stripped, and using those coordinates to find the bounding box of 3D object. What should be the proper way to voxelize the mesh?

import objParser
import math

inputFile = 'test.obj'
vList = []; vtList = []; vnList = []; fList = []

def parseOBJ(inputFile):
list = []
vList, vtList, vnList, fList = objParser.getObj(inputFile)
print 'in parseOBJ'
#print vList, vtList, vnList, fList
return vList, vtList, vnList, fList

def findBBox(vList):
   i = 0; j=0; x_min = float('inf'); x_max = float('-inf'); y_min = float('inf');
   y_max =  float('-inf'); z_min = float('inf'); z_max = float('-inf');
   xWidth = 0; yWidth = 0; zWidth =0

print 'in findBBox'
while i < len(vList):
        #find min and max x value
        if vList[i][j] < x_min:
            x_min = float(vList[i][j])
        elif vList[i][j] > x_max:
            x_max = float(vList[i][j])

        #find min and max y value
        if vList[i][j + 1] < y_min:
            y_min = float(vList[i][j + 1])
        elif vList[i][j + 1] > y_max:
            y_max = float(vList[i][j + 1])

        #find min and max x value
        if vList[i][j + 2] < z_min:
            z_min = vList[i][j + 2]
        elif vList[i][j + 2] > z_max:
            z_max = vList[i][j + 2]

        #incriment the counter int by 3 to go to the next set of (x, y, z)
        i += 3; j=0

xWidth = x_max - x_min
yWidth = y_max - y_min
zWidth = z_max - z_min
length = xWidth, yWidth, zWidth
volume = xWidth* yWidth* zWidth
print 'x_min, y_min, z_min : ', x_min, y_min, z_min
print 'x_max, y_max, z_max : ', x_max, y_max, z_max
print 'xWidth, yWidth, zWidth : ', xWidth, yWidth, zWidth
return length, volume

def init():
    list = parseOBJ(inputFile)
    findBBox(list[0])

print init()

Answer 1

I haven't used it but you can try this one: http://packages.python.org/glitter/api/examples.voxelization-module.html

Or this tool: http://www.patrickmin.com/binvox/

If you want to do this by yourself you have two main approaches:

• A 'real' voxelization, when you take into account the mesh interior - rather complex and you need much of CSG operations. I can't help you there.

• A 'fake' one - just voxelize every triangle of your mesh. This is much simpler, all you need to do is to check the intersection of a triangle and axis-aligned cube. Then you just do:

  for every triagle:
      for every cube:
          if triangle intersects cube:
              set cube = full
          else:
              set cube = empty
  

All you need to do is to implement a BoundingBox-Triangle intersection. Ofcourse you can optimize those for loops a bit :)

Answer 2

For every one that still has this issue. I have written a file the makes a voxel of a STL 3d file (you can save .obj files as stl). I used the method suggested by kolenda to solve this. This is for so calles 'fake' voxelization. Still working on filling this voxel in. I used numpy.stl to import the file and just standard packages for the rest of the file. I used the following links for my program.

For line plane collision

to check if point is in triangle

It is probably not the most efficiant but it works.

import numpy as np
import os
from stl import mesh
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt
import math


if not os.getcwd() == 'path_to_model':
    os.chdir('path_to model')


your_mesh = mesh.Mesh.from_file('file.stl') #the stl file you want to voxelize


## set the height of your mesh
for i in range(len(your_mesh.vectors)):
    for j in range(len(your_mesh.vectors[i])):
        for k in range(len(your_mesh.vectors[i][j])):
            your_mesh.vectors[i][j][k] *= 5

## define functions
def triangle_voxalize(triangle):
    trix = []
    triy = []
    triz= []
    triangle = list(triangle)

    #corners of triangle in array formats
    p0 = np.array(triangle[0])
    p1 = np.array(triangle[1])
    p2 = np.array(triangle[2])

    #vectors and the plane of the triangle
    v0 = p1 - p0
    v1 = p2 - p1
    v2 = p0 - p2
    v3 = p2 - p0
    plane = np.cross(v0,v3)

    #minimun and maximun coordinates of the triangle
    for i in range(3):
        trix.append(triangle[i][0])
        triy.append(triangle[i][1])
        triz.append(triangle[i][2])
    minx, maxx = int(np.floor(np.min(trix))), int(np.ceil(np.max(trix)))
    miny, maxy = int(np.floor(np.min(triy))), int(np.ceil(np.max(triy)))
    minz, maxz = int(np.floor(np.min(triz))), int(np.ceil(np.max(triz)))

    #safe the points that are inside triangle
    points = []

    #go through each point in the box of minimum and maximum x,y,z
    for x in range (minx,maxx+1):
        for y in range(miny,maxy+1):
            for z in range(minz,maxz+1):

                #check the diagnals of each voxel cube if they are inside triangle
                if LinePlaneCollision(triangle,plane,p0,[1,1,1],[x-0.5,y-0.5,z-0.5],[x,y,z]):
                    points.append([x,y,z])
                elif LinePlaneCollision(triangle,plane,p0,[-1,-1,1],[x+0.5,y+0.5,z-0.5],[x,y,z]):
                    points.append([x,y,z])
                elif LinePlaneCollision(triangle,plane,p0,[-1,1,1],[x+0.5,y-0.5,z-0.5],[x,y,z]):
                    points.append([x,y,z])
                elif LinePlaneCollision(triangle,plane,p0,[1,-1,1],[x-0.5,y+0.5,z-0.5],[x,y,z]):
                    points.append([x,y,z])

                #check edge cases and if the triangle is completly inside the box
                elif intersect(triangle,[x,y,z],v0,p0):
                    points.append([x,y,z])
                elif intersect(triangle,[x,y,z],v1,p1):
                    points.append([x,y,z])
                elif intersect(triangle,[x,y,z],v2,p2):
                    points.append([x,y,z])

    #return the points that are inside the triangle
    return(points)

#check if the point is on the triangle border
def intersect(triangle,point,vector,origin):
    x,y,z = point[0],point[1],point[2]
    origin = np.array(origin)

    #check the x faces of the voxel point
    for xcube in range(x,x+2):
        xcube -= 0.5
        if LinePlaneCollision(triangle,[1,0,0], [xcube,y,z], vector, origin,[x,y,z]):
            return(True)

    #same for y and z
    for ycube in range(y,y+2):
        ycube -= 0.5
        if LinePlaneCollision(triangle,[0,1,0], [x,ycube,z], vector, origin,[x,y,z]):
            return(True)
    for zcube in range(z,z+2):
        zcube -= 0.5
        if LinePlaneCollision(triangle,[0,0,1], [x,y,zcube], vector, origin,[x,y,z]):
            return(True)

    #check if the point is inside the triangle (in case the whole tri is in the voxel point)
    if origin[0] <= x+0.5 and origin[0] >= x-0.5:
        if origin[1] <= y+0.5 and origin[1] >= y-0.5:
            if origin[2] <= z+0.5 and origin[2] >= z-0.5:
                return(True)

    return(False)

# I modified this file to suit my needs:
# https://gist.github.com/TimSC/8c25ca941d614bf48ebba6b473747d72
#check if the cube diagnals cross the triangle in the cube
def LinePlaneCollision(triangle,planeNormal, planePoint, rayDirection, rayPoint,point, epsilon=1e-6):
    planeNormal = np.array(planeNormal)
    planePoint = np.array(planePoint)
    rayDirection = np.array(rayDirection)
    rayPoint = np.array(rayPoint)


    ndotu = planeNormal.dot(rayDirection)
    if abs(ndotu) < epsilon:
        return(False)

    w = rayPoint - planePoint
    si = -planeNormal.dot(w) / ndotu
    Psi = w + si * rayDirection + planePoint

    #check if they cross inside the voxel cube
    if np.abs(Psi[0]-point[0]) <= 0.5 and np.abs(Psi[1]-point[1]) <= 0.5 and np.abs(Psi[2]-point[2]) <= 0.5:
        #check if the point is inside the triangle and not only on the plane
        if PointInTriangle(Psi, triangle):
            return (True)
    return (False)

# I used the following file for the next 2 functions, I converted them to python. Read the article. It explains everything way better than I can.
# https://blackpawn.com/texts/pointinpoly#:~:text=A%20common%20way%20to%20check,triangle%2C%20otherwise%20it%20is%20not.
#check if point is inside triangle
def SameSide(p1,p2, a,b):
    cp1 = np.cross(b-a, p1-a)
    cp2 = np.cross(b-a, p2-a)
    if np.dot(cp1, cp2) >= 0:
        return (True)
    return (False)

def PointInTriangle(p, triangle):
    a = triangle[0]
    b = triangle[1]
    c = triangle[2]
    if SameSide(p,a, b,c) and SameSide(p,b, a,c) and SameSide(p,c, a,b):
        return (True)
    return (False)

##
my_mesh = your_mesh.vectors.copy() #shorten the name

voxel = []
for i in range (len(my_mesh)): # go though each triangle and voxelize it
    new_voxel = triangle_voxalize(my_mesh[i])
    for j in new_voxel:
        if j not in voxel: #if the point is new add it to the voxel
            voxel.append(j)

##
print(len(voxel)) #number of points in the voxel

#split in x,y,z points
x_points = []
y_points = []
z_points = []
for a in range (len(voxel)):
    x_points.append(voxel[a][0])
    y_points.append(voxel[a][1])
    z_points.append(voxel[a][2])

## plot the voxel
ax = plt.axes(projection="3d")
ax.scatter3D(x_points, y_points, z_points)
plt.xlabel('x')
plt.ylabel('y')
plt.show()

## plot 1 layer of the voxel
for a in range (len(z_points)):
    if z_points[a] == 300:
        plt.scatter(x_points[a],y_points[a])

plt.show()

Answer 3

You can achieve this with pymadcad The class PositionMap has a very optimized method to rasterize points, lines and triangles into voxels.

This is only filling the voxels on the surface of the mesh, and not filling the interior. But using those functions the interior will always be separated from the exterior ;)

here what it should look:

from madcad.hashing import PositionMap
from madcad.io import read

# load the obj file in a madcad Mesh
mymesh = read('mymesh.obj')
# choose the cell size of the voxel
size = 1

voxel = set()    # this is a sparse voxel
hasher = PositionMap(size)   # ugly object creation, just to use one of its methods
for face in mymesh.faces:
    voxel.update(hasher.keysfor(mymesh.facepoints(face)))

Yeah, it's not that pretty, beacause even though the functions exist, madcad has not (yet) any idiomatic way to do it.

explanations

• instead of a 3d array of booleans (which is very costly and overkill to store mostly zeroes), we use a set that stores only locations of non-zero cells
• the method hasher.keysfor generates all the locations of the voxel's cells reached by the input primitive (triangles here)
• a set is a hashtable and allow to check very efficiently wheter a an element is inside or not (eg. can check very efficiently if a cell is occupied by the primitive)
• So the loop is just updating the voxel with the cells occupied by the primitives, the occupied cells are then the cells on the surface of the mesh

an other way to perform pathfinding

It seems that you are trying to use some pathfinding to find your shortest distance inside this building. voxelize the area is a good approach. There is an other one if you need to try it:

A pathfinding (like A* or dikstra) usually works on a graph. It can be connected cells in a voxel, or any kind of irregularly spaced cells.

So an other solution would be to triangulate your 3d space between your building walls by generating tetraedrons. Instead of propagating you algorithm from voxel cell to voxel cell, you would propagate from tetraedron to tetraedron. Tetraedrons do have the advantage of being faster to generate from random meshes, but also to not need a memory propotionate to the volume of the area, and do not loose any obstacle precision (tetraedrons have adaptative sizes).

You can have a look here at what it looks like.