Introduction

Create a Voxel mesh from STL CAD data using the Signed Distance Function (SDF). The library used is Python's VTK. In addition, the code was entirely based on the following article.

https://qiita.com/torisan_piyopiyo/items/588e607c578b4b00903a

environment

Windows 10 home

Anaconda(Python 3.7.6)

VTK 8.2

Paraview 5.4.1 (used to display STL and voxel mesh)

CAD data (STL data)

Stanford Dragon.

Source code

import vtk

filename = "./dragon.stl"
#Read stl file
reader = vtk.vtkSTLReader()
reader.SetFileName(filename)
reader.Update()

#Extract Poly data
poly_data = reader.GetOutput()

#Determine the mesh size
cell_dims = [50, 50, 50]#x, y, z

#Create a mesh grid
#x_min:0 x_max:1, y_min:2,y_max:3,z_min:4,z_max:5
bounds = poly_data.GetBounds()
mesh_pitch = [(bounds[1] - bounds[0])/cell_dims[0],
              (bounds[3] - bounds[2])/cell_dims[1],
              (bounds[5] - bounds[4])/cell_dims[2]]
mins = [bounds[0], bounds[2], bounds[4]]

points = vtk.vtkPoints()
for ix in range(cell_dims[0]+1):
    for iy in range(cell_dims[1]+1):
        for iz in range(cell_dims[2]+1):
            x = ix * mesh_pitch[0] + mins[0]
            y = iy * mesh_pitch[1] + mins[1]
            z = iz * mesh_pitch[2] + mins[2]
            points.InsertNextPoint(x,y,z)
structured_base_mesh = vtk.vtkStructuredGrid()
structured_base_mesh.SetExtent(0, cell_dims[0], 0, cell_dims[1], 0, cell_dims[2])
structured_base_mesh.SetPoints(points)

#Convert structured grid data to unstructured grid data
append = vtk.vtkAppendFilter()
append.AddInputData(structured_base_mesh)
append.Update()
base_mesh = append.GetOutput()

#Find the coordinates of the center of Voxel
cell_centers = vtk.vtkCellCenters()
cell_centers.SetInputData(base_mesh)
cell_centers.Update()

#Create a Voxel mesh from STL
center_points = cell_centers.GetOutput().GetPoints()
cell_list = vtk.vtkIdList()
sdf = vtk.vtkImplicitPolyDataDistance()
sdf.SetInput(poly_data)
distance_sdf = vtk.vtkDoubleArray()
distance_sdf.SetName("sdf")
for idx in range(center_points.GetNumberOfPoints()):
    current_center = center_points.GetPoint(idx)
    distance = sdf.FunctionValue(current_center)
    distance_sdf.InsertNextValue(distance)
    if distance <= 0:
        cell_list.InsertNextId(idx)

base_mesh.GetCellData().SetScalars(distance_sdf)#Add SDF value

extract_cells = vtk.vtkExtractCells()
extract_cells.SetInputData(base_mesh)
extract_cells.SetCellList(cell_list)
extract_cells.Update()

#Output Voxel data as VTK file
writer = vtk.vtkXMLDataSetWriter()
writer.SetFileName("closed_voxel.vtu")
writer.SetInputData(extract_cells.GetOutput())
writer.Update()

#Output the grid data with the SDF value input as a VTK file
#writer = vtk.vtkXMLDataSetWriter()
#writer.SetFileName("all_voxel.vtu")
#writer.SetInputData(base_mesh)
#writer.Update()