在Python中，如何将一个三维网格模型进行体素化处理？

Question

我需要帮助开始学习Python（我几乎一无所知），以将由Rhino生成的3D网格进行体素化。数据输入将是一个.OBJ文件，输出也是如此。这个用途的最终目的是在建筑物内找到两点之间的最短距离。但那是以后的事情。现在，我需要先对3D网格进行体素化。体素化原语可能只是一个简单的立方体。

到目前为止，我可以从OBJ文件解析器中读取数据，并从解析后的OBJ中去除V、VT、VN、F前缀，使用这些坐标来找到3D对象的边界框。应该采取什么合适的方式来对网格进行体素化呢？

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()

- FongYu

Answer 1

我虽然没有使用过，但你可以试试这个：http://packages.python.org/glitter/api/examples.voxelization-module.html

或者尝试这个工具：http://www.patrickmin.com/binvox/

如果你想自己完成这个任务，有两种主要方法：

• 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
  

您只需要实现一个边界框和三角形相交的函数即可。当然，您可以对那些循环进行一些优化 :)

- kolenda

Answer 2

对于仍然存在此问题的每个人。我已经编写了一个文件，可以将STL 3D文件制作成体素（您可以将.obj文件保存为stl）。我使用了kolenda建议的方法来解决这个问题。这是用于所谓的“假”体素化。仍在努力填充这个体素。

我使用numpy.stl导入文件，其余文件使用标准包。我使用了以下链接来进行编程。 用于线平面碰撞 检查点是否在三角形中 它可能不是最有效的，但它可以工作。

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()

- Wycher

Answer 3

您可以使用pymadcad来实现这一点。

PositionMap类具有非常优化的方法，可以将点、线和三角形光栅化为体素。

这仅填充网格表面上的体素，而不是填充内部。但是使用这些函数，内部始终与外部分离；）

这是它应该的样子：

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)))

是的，这并不太美观，因为尽管函数存在，但madcad还没有（尚未）任何惯用的方法来执行它。

解释

• 我们使用一个set而不是一个3D布尔数组（存储大多数零非常昂贵和过度），该set仅存储非零单元格的位置
• 方法hasher.keysfor生成由输入基元（这里是三角形）到达的体素单元格的所有位置
• set是一个哈希表，可以非常高效地检查元素是否在内部（例如，可以非常高效地检查单元格是否被基元占据）
• 因此，循环只是使用由基元占据的单元格更新体素，占据的单元格然后是网格表面上的单元格

执行路径查找的另一种方法

看起来您正在尝试使用一些路径查找来找到建筑物内的最短距离。对区域进行体素化是一个好方法。如果需要尝试其他方法，则还有另一种方法：

路径规划（如A*或Dijkstra算法）通常在图形上运作。它可以是体素中连接的单元格，也可以是任何类型的不规则间隔单元格。

因此，另一种解决方案是通过生成四面体来三角化建筑墙壁之间的3D空间。您不需要从体素单元格到体素单元格传播算法，而是需要从四面体到四面体进行传播。 四面体具有从随机网格更快生成的优点，而且不需要与区域体积成比例的内存，并且不会失去任何障碍物精度（四面体具有自适应大小）。

您可以在这里看看它的样子。

- jimy byerley

Answer 4

对于仍然感兴趣的人，您可以尝试voxeltool，它似乎是专门为此目的而制作的。

从文档中可以看到：

体素工具为Rhino提供了轻量级的体素几何。它允许您快速从网格、Brep、曲线和点生成和操作体素化几何，并在体素网格之间提供布尔运算。它可以将体素网格转换为实体网格外壳。

- Thiago Petrone