The world’s leading publication for data science, AI, and ML professionals.

How to Voxelize Meshes and Point Clouds in Python

Step-by-step tutorial on voxelization using Open3D, Trimesh, PyVista, and pyntcloud – extracting features and creating interactive visuals

This article goes through the steps of generating voxel representations of point clouds and meshes using four widely popular Python libraries – Open3D, Trimesh, PyVista, and pyntcloud. Voxelization is an important pre-processing step for a lot of 3D deep learning models. The article shows how to calculate voxel level features like colors, point density, and occupancy among others. Finally, a demonstration of how to create simple interactive voxelization and thresholding examples is also given. Who said that voxelization should be complicated?

Examples of voxelization of a point cloud in Open3D - different voxel sizes (left) and construction of a voxel grid (right) | Image by the author
Examples of voxelization of a point cloud in Open3D – different voxel sizes (left) and construction of a voxel grid (right) | Image by the author

Deep Learning for 3D data is becoming a more and more important part of machine learning and understanding the world around us. As new 3D data extraction hardware like depth cameras and LiDARs, is becoming commonplace in CCTVs, cameras, and smartphones more and more people are using the additional dimension that it provides. In addition, photogrammetry and Structure from Motion are becoming a normal part of the 3D reconstruction and modeling pipelines, and extracting and manipulating large 3D datasets is becoming a necessity. Unstructured data for 3D deep learning can have different representations:

  • Point clouds [3,4, 11]
  • Voxels and voxel grids [1, 8, 9]
  • Depth maps [2]
  • Graphs [5, 10]

These are far from all possible 3D data presentations, with others like parametric CAD models, multi-view images, volumes, etc. For a really good overview of some of these, you can read the informative article by Florent Paux on "How to represent 3D Data ?".

In this article, we will focus on representing 3D data as voxels. But first what is a voxel? The simplest comparison is that a voxel is a 3D pixel. Voxels are ordered into voxel grids, which can be seen as the 3D equivalent of the ordered structure of images. When a point cloud or mesh is turned into a voxel representation, it is intersected with a voxel grid. Points in the point cloud or mesh then fall in certain voxels. These voxels are left, while all others that do not intersect any points are either discarded or just zeroed out and what we are left with is a sculpted representation of the object. The voxelization can just be surface level or throughout the whole mesh/point cloud volume.

Example of voxelization where the angel statue on the left is transformed into a voxel representation on the right using PyVista | Image by the author
Example of voxelization where the angel statue on the left is transformed into a voxel representation on the right using PyVista | Image by the author

Building voxelized representations of meshes and point clouds is an important step in data preprocessing for many deep learning methods. Voxelization is also widely used to process point clouds – subsampling, feature extraction, and occupancy analysis, among others. Finally, generating voxel representation of meshes can be also useful for creating assets for games and simplifying surfaces for simulations.

In the article, we will explore the voxelization capabilities of four Python libraries – Open3D, Trimesh, PyVista, and pyntcloud. I have selected these libraries as they provide relatively straightforward voxelization capabilities, together with built-in solutions for the analysis of the created voxels. Other smaller libraries like Py[Voxelizer](https://github.com/archi-coder/Voxelizer), Voxelizer, simple-3dviz, and DepthVisualizer also provide voxelization functions, but they are deemed too limited. Another library that provides voxelization functionality is Vedo, but after initial tests, I have seen that they are directed towards 3D Volume data, which limits their utilization. If you want to read more on how to create beautiful visualizations using these libraries and many more, you can look at my previous articles on "Python Libraries for Mesh, Point Cloud, and Data Visualization" (Part 1 and 2). The only library we will use that has not been previously discussed is pyntcloud, so we will go into more detail about how to install and set it up.

Python Libraries for Mesh, Point Cloud, and Data Visualization (Part 1)

Python Libraries for Mesh, Point Cloud, and Data Visualization (Part 2)

To demonstrate the voxelization on both point clouds and meshes, I have provided two objects. First, a bunny statue point cloud in .txt format, which contains the X, Y, and Z coordinates of each point, together with their R, G, and B colors, and finally the Nx, Ny, and Nz normals. Second, a rooster statue mesh in a .obj format, together with a .mat file and a texture in .jpg format. Both objects were created using Structure from Motion photogrammetry and are free to use in commercial, non-commercial, public, and private projects. These objects are a part of a larger dataset [6] and have been used in the development of methods for noise detection and inspection [7]. To follow the tutorials, in addition to the used libraries and their dependencies, you will also need NumPy and SciPy. All the code is available at the GitHub repository HERE


Voxelization using Open3D

Open3D result of voxelization with different sizes of voxel grids | Image by the author
Open3D result of voxelization with different sizes of voxel grids | Image by the author

Open3D is one of the most feature-rich Python libraries for 3D analysis, mesh and point cloud manipulation, and visualization. It contains excellent tools for generating voxels from both point clouds and meshes. The library also has support for a large array of 3D primitives like spheres, cubes, cylinders, etc. which makes it easy to generate 3D meshes from the voxels. Finally, the voxel grids can be also used to subsample point clouds and generate segmentations. For a step-by-step tutorial on how to install and use Open3D, you can take a look at my previous article "Python Libraries for Mesh, Point Cloud, and Data Visualization (Part 1)". We will be using a lot of the same code for visualization and callbacks.

Again I would like to point out an excellent tutorial by Florent Poux, Ph.D. – "How to Automate Voxel Modelling of 3D Point Cloud with Python", where he gives a good overview and use cases for the generation of voxels from point clouds using Open3D. In my article, I build upon the information provided in the tutorial.

For a better understanding of the voxelization process, we will build two examples using Open3D – one that visualizes multiple voxel grid levels and how they affect the voxelized model and another one that shows the voxelization process.

We start by first loading the bunny statue point cloud using NumPy’s loadtxt() function. We separate the array into points, colors, and normals and create an Open3D point cloud by first transforming the arrays into vectors using theVector3dVector function. From here creating a voxel grid is trivial by calling VoxelGrid.create_from_point_cloud() and giving it the point cloud and voxel size parameter. If we are starting from a mesh we can do the same thing with create_from_triangle_mesh(). If we are using the voxel grid to manipulate and analyze the point cloud/mesh in Open3D we can stop here. We can also quickly visualize the voxel grid by calling by adding it as a geometry to a Visualizer object. The code for this is given below.

We would also to continue processing the voxel grid and transform it into a mesh because we would like to create an animation with it, and as of the time of writing voxel grids cannot be transformed. We use the method proposed by Florent Poux, Ph.D. We create an empty triangle mesh by calling geometry.TriangleMesh() , we then extract all voxels from the grid by calling voxel_grid.get_voxels() and getting the true size of the created voxels. This way we can create an Open3D box primitive at each voxel position, scale the boxes to the size of the voxels and position them at the voxel centers by calling voxel_grid.get_voxel_center_coordinate(). We finally add the thus created, scaled, and positioned box to the TriangleMesh object.

To demonstrate how different voxel grid sizes change the final voxel output and to make the example a little bit more interesting we will animate a rotating transition between the voxel outputs. To do this we register a new callback to the Visualizer by calling register_animation_callback(). In the callback function, we generate the voxel grid and create a voxel mesh from it we rotate the thus created mesh 360 degrees, destroy the mesh and repeat with a voxel grid of different sizes. We also utilize a small class to hold all variables that will change in the callback update loop. The code for the full rotation visualization is given below.

To get a better understanding of the process of creating the voxels and the voxel mesh, we will visualize it through animation. The code for this visualization is very similar to the one that we have already shown, but this time we visualize the process of creating the voxel mesh, by adding each new voxel as part of the callback function loop. To have the outline of the bunny statue, we also use the voxel function get_voxel_center_coordinates to extract the voxel centers and use them to generate a point cloud. This way we generate in a roundabout way a subsampled and uniformly sampled version of the input point cloud. The code for the whole visualization is given below.

Visual representation of generating the voxel mesh from the input point cloud in Open3D| Image by the author
Visual representation of generating the voxel mesh from the input point cloud in Open3D| Image by the author

Voxelization using Pyntcloud

Pyntcloud result of voxelization, voxels colored based on density | Image by the author
Pyntcloud result of voxelization, voxels colored based on density | Image by the author

Pyntcloud is a lightweight and powerful Python 3 library that is directed toward the analysis and pre-processing of point clouds and meshes. It contains tools for voxelization, feature extraction, mesh smoothing and decimation, point cloud subsampling, outlier detection, and many more. It is directed use is as part of jupyter notebooks and contains visualization bindings for threejs, pythreejs, PyVista, and matplotlib. It can be easily integrated into the workflows of other larger libraries like Open3D and PyVista.

The voxelization capabilities of pyntcloud are quite robust, with voxel representations, already containing a number of pre-calculated features like voxel to points correspondences, voxel grid segment information, the density of points in each voxel, and binary mask, among others.

Pyntcloud works with Python 3 on Linux, Mac, and Windows. It requires NumPy and Pandas as pre-requisites. The library is in constant development and gets updated monthly. It can be installed using either pip or Anaconda. We normally install each library in a new Anaconda environment to keep everything neat, but in this case, we will install pyntcloud together with Open3d to leverage the visualization capabilities of the former. We do this because voxel and mesh visualizations in pyntcloud require threejs or pythreejs, which in turn require either jupyter notebooks or starting a local server. These are outside of the scope of the article, but I will still give examples of how to visualize the voxels directly in pyntcloud.

conda create -n pyntcloud_env python=3.8
conda activate pyntcloud_env
pip install open3d
pip install pyntcloud
OR
conda install pyntcloud -c conda-forge

Once we have installed everything and no errors are present, we can test both libraries by first calling import open3d as o3d together with print(o3d.__version__) and then from pyntcloud import PyntCloud together with print(PyntCloud.__version__). If both libraries are imported we can progress with the tutorial.

The I/O module of pyntcloud is quite simplified for better or worse. It uses PyntCloud.from_file() to load both meshes and point clouds. The function decides what internal functions to call based on the extension of the file that needs to be loaded. For a list of all the supported file types of pyntcloud, you can read their I/O page HERE. For ASCII files containing point clouds, pyntcloud directly calls Pandas and requires additional input parameters like a separator, names of columns, is there a header part, etc. In the case of this tutorial, we import the bunny statue point cloud and need to explicitly specify the names of all the columns. Once we have loaded the point cloud, we need to call .add_structure("voxelgrid", grid_size_x, grid_size_y, grid_size_z) to create a voxel grid on top of the point cloud. The output of this function is the voxel grid ids where points from the point cloud intersect the grid. We then use these ids to create the occupied voxel representation of the point cloud by calling name_of_point_cloud.structures[voxelgrid_ids]. From here on we can use this object to extract information for the voxels and visualize everything. The code for importing the point cloud and extracting the voxels is given below.

Once we have created the voxel representation of the point cloud, we can directly visualize it using pyntcloud by calling voxelgrid.plot(cmap="hsv", backend="threejs"). Different backends can be called here, but for meshes and voxels only threejs and pythreejs are supported. To be able to easier visualize everything without the need of a running server we will transfer the created voxel grid to Open3D and visualize it using the technique shown in the Open3D part. This way we also show how easy it is to combine these two libraries, leveraging their strengths.

We first create an empty Open3D geometry.TriangleMesh(), which will be filled with voxels. We then loop through all voxels on the pyntcloud voxel grid. We can get the current voxel id from the voxelgrid.voxel_n array and we can use this id to get its voxel center by calling voxelgrid.voxel_centers[voxel_id]. To get the density value that we will use as a color for the voxels we use the X, Y, and Z positions of the voxel in the voxel grid, that we already extracted. From here everything is the same as in the example of Open3D. We color the box primitive, scale it and translate it to the voxel center. Finally, we initialize an Open3D visualizer and add the created voxel mesh as a geometry to it. The code for the full pyntcloud/Open3D visualization is given below.

Voxel meshes using cones (top left), cylinders (top right), octahedrons (bottom left), and spheres (bottom right) | Image by the author
Voxel meshes using cones (top left), cylinders (top right), octahedrons (bottom left), and spheres (bottom right) | Image by the author

To have some more fun with this we can substitute the box primitives when creating the voxel mesh with other primitives. Open3D has a large number of primitives – spheres, cylinders, toruses, octahedrons, cones, etc. We create a small selection function that takes the name of the primitive and returns the 3D object. Then the only thing that we change is primitive = o3d.geometry.TriangleMesh.createbox() to the new function call primitive = choose_primitive("cylinder"). The code for the function is given below.

Voxelization using Trimesh

Process for creating a color voxelized mesh in Trimesh - input mesh with texture (left), mesh with color information mapped to each vertex (center), and colored voxelized mesh (right) | Image by the author
Process for creating a color voxelized mesh in Trimesh – input mesh with texture (left), mesh with color information mapped to each vertex (center), and colored voxelized mesh (right) | Image by the author

Trimesh is one of the most widely-used Python libraries for working with meshes. A large number of other libraries and stand-alone projects use it or are directly built on top of it. It provides a large array of point cloud and mesh analysis and transformation tools. For installation instructions for Trimesh and examples of visualization use, you can look at the article "Python Libraries for Mesh, Point Cloud, and Data Visualization (Part 1)", where we show how to work on both meshes and point datasets.

For voxelization, Trimesh has a number of methods for creating voxel grids, hollowing them inside, as well as filling them using morphological operations like closing and binary dilation. It also has direct methods for extracting the surface or the fill of a voxel mesh by callingtrimesh.voxel.morphology.surface and trimesh.voxel.morphology.fill respectively. The library also gives tools for quick analysis of the voxelized mesh, like extracting bounds, checking if a voxel space in the grid is empty or filled, and calculating the scale and volume of the occupied voxels. A function for performing marching cubes calculation is also readily available.

One functionality that is quite tricky with Trimesh is extracting a mesh’s texture color information and giving it as colors to the voxels representation. To do this we first extract the texture information and map it as color information to each vertex of the mesh. This is done by calling name_of_the_imported_mesh.visual.to_color().vertex_colors. We can then transform this into a NumPy array and freely access it for each mesh vertex. We can then create a voxel grid from the mesh by directly calling name_of_the_imported_mesh.voxelized().hollow() where we can specify the resolution of the voxel grid. We also transform the voxelization into a hollow representation, as we will need to calculate distances between the mesh and the voxel representation.

As Trimesh does not contain functionality for direct mapping of mesh colors to the voxels of the voxel representation we first need to call trimesh.proximity.ProximityQuery(name_of_the_mesh).vertex(voxel_grid_name.points)to calculate the distance between each vertex of the mesh and the already created voxel representation points. We then need to go through each point, extract the color of the closest mesh vertex, and map it into a color array of size (X, Y, Z, 4), where X, Y, and Z are the sizes of the generated voxel grid and 3 is the R, G, B, and A components of the color. Once we have the mapping we can call voxel_grid_name.as_boxes(colors = mapped_color_array), which makes the voxel grid into a Trimesh mesh object with voxel/box colors equal to the initial mesh ones. One thing that readers will see when running the code is that the voxelized mesh has holes in it. This most probably happened because the rooster statue mesh is not perfectly watertight because of the Structure from Motion process. The code for the voxelization process is given below.

Voxelization using PyVista

PyVista multi-plot visualization of mesh (top left), voxelized representation (top right), cones representation instead of cubes (bottom left), and per voxel distance visualization to the mesh (bottom right) | Image by the author
PyVista multi-plot visualization of mesh (top left), voxelized representation (top right), cones representation instead of cubes (bottom left), and per voxel distance visualization to the mesh (bottom right) | Image by the author

PyVista is a fully-featured mesh visualization and analysis library, built on top of VTK. It can be used to create complex applications containing multiple windows, GUI elements, and interactions through the mouse, keyboard, and GUI. The library uses filtering routines to transform point clouds and meshes and to extract any required information. It can be used as part of physical simulations, computer graphics, GIS, architectural visualizations, ray tracing, and mathematical computation showcasing. For installation instructions for PyVista and examples of visualization use, you can look at the article "Python Libraries for Mesh, Point Cloud, and Data Visualization (Part 1)", where we show how to work on both meshes and point datasets.

As PyVista is built on top of VTK, it contains a straightforward interface for transforming meshes into UniformGrid objects for building voxel representations. An important point to mention is that PyVista can be used for voxelization of only meshes. If a point cloud is given, it first needs to be triangulated into a mesh. For simpler point cloud surfaces this can be done directly in PyVista, using delaunay_3d triangulation. Another important point is that PyVista expects watertight meshes by default for voxelization. If the mesh has holes, overlapping edges, or other incorrect geometry, trying to voxelize it will result in an error. The voxelization is done by calling pyvista.voxelize(input_mesh). if this method throws an error that the used mesh is not watertight, there are two ways to proceed – either fix the mesh in another software like Blender or add the parameter check_surface=False. In case of small imperfections in the mesh, the second way will yield good results without any noise or artifacts. Be warned that in case there are large holes or geometric imperfections skipping the surface check may result in incorrect voxelization and noise.

Once voxelize is called the result is a uniform grid object that contains all the voxels. Additional information can then be extracted from the representation like voxel centers, coordinates of the X, Y, Z bounds of each voxel, convex hull of the voxels containing parts of the mesh, etc. A useful analysis of the voxel representation is given in the examples of PyVista. Once we have the voxels, we can call compute_implicit_distance to calculate the distance between the voxels and the mesh. This can be useful for inspecting if there have been any inconsistencies with the voxelization process, as well as used to filter out voxels inside of the mesh to create a hollow voxel representation. Another useful transformation is using different primitive objects than cubes to represent the voxels. This can be directly set when calling voxelize or can be created later by calling voxelized_object.glyph(geom = type_of_primitive()). Finally, for this part, we will also visualize everything in separate sub-windows and link them. This is done in the same way as matplotlib, where the number of windows is given when setting up the plotter pyvista.Plotter(shape=(num_horiz_windows,num_vertical_windows)). Once this is set then each sub-window can be invoked by calling pyvista.subplot(horizontal_num, vertical_num). The code for producing four sub-window visualizations of the input mesh, voxelized representation with cubes, voxelized representation with cones, and finally distances between mesh and voxels are given below.

PyVista thresholding voxel mesh based on point density in each voxel, using GUI widgets | Image by the author
PyVista thresholding voxel mesh based on point density in each voxel, using GUI widgets | Image by the author

We will use the voxelization context to show how to also create a simple GUI in PyVista and add interactivity to the generating process. First, a short example demonstrating how to generate point density for each voxel. This was chosen as this is a useful analysis step, which is not present in PyVista out of the box. For this, the KDTree functionality of SciPy is required. We use it to get all neighboring mesh vertices to each voxel center. SciPy has the function query_ball_point, which makes a sphere with a center at a point and a specified radius and gets all neighboring points to it. Of course, as we are searching for points inside of a voxel cube shape, using a sphere will result in imperfect results. For this tutorial, we will use this simplification. To calculate the radius of the sphere, so it has the same volume as the voxel cube we use the known cube side x to first find the volume of the cube:

Next, we get the volume of the sphere:

From there if want the two volumes to be equal, we can calculate the radius as :

This can be written in Python in a small helper function as seen below.

Once we have all vertices in the neighborhood of each voxel center we count them and divide them by the maximum count to normalize them between 0 and 1 for visualization. Finally, we add these values as a field to the voxels. The code is separated into another function given below.

PyVista contains a ready-to-use built-in visualization tool for thresholding meshes and point clouds called add_mesh_threshold. This directly adds the object to the scene and creates a slider widget connected to the most recent active field or a specified one. The code for the full visualization is given below.

PyVista custom GUI widget for interactively changing voxelization voxel size | Image by the author
PyVista custom GUI widget for interactively changing voxelization voxel size | Image by the author

Finally, we will create a custom slider widget to demonstrate how easy it is in PyVista. Moving the slider will change the size of the voxels, with which the mesh is voxelized, and recalculate the whole thing. For a more neat package, we will add everything to a class. To create a slider widget we just need to call add_slider_widget(name_of_callback_function, [min, max], start_value, title_of_widget, how_often_should_the_widget_call_callback). The last input determines if the callback function will be invoked when we click the mouse button, release it, or every time the mouse is moved. In our case, we use the last option. In the callback function, we do the voxelization calculation and add it to the Plotter. Here it’s important to specify the name of the mesh object in the plotter. This way PyVista knows that we are updating the same mesh and not creating a new one every time we move the slider. The code for the whole GUI interactivity is given below.

Conclusion

Voxelization can be one of the strongest tools in your arsenal when working with point clouds and meshes. Using Python to generate and work with voxels can be straightforward. In this article, we explored how to create a voxel grid, extract features from voxels, and change the primitives from cubes to spheres, cylinders, cones, etc. We also showed how the voxel building process is done, as well as how to create interactive voxelization and thresholding tools.

Once you have a voxel representation, it can be easy to plug it into existing deep learning or asset creation pipelines for achieving the required results. In the next articles, we will explore more ways to extract features and useful information from point clouds and meshes through adjacency analysis, PCA feature extraction, RANSAC implementations, and many more.

If you want to read more about extracting features from point clouds and meshes, you look at some of my articles on 3D surface inspection and noise detection [11] and [7]. You can find the articles, plus my other research on my Page, and if you see something interesting or just want to chat feel free to drop me a message. Stay tuned for more!

References

  1. Qi, C. R., Su, H., Nießner, M., Dai, A., Yan, M., & Guibas, L. J. (2016). Volumetric and multi-view cnns for object classification on 3d data. In Proceedings of the IEEE conference on computer vision and pattern recognition (pp. 5648–5656); https://openaccess.thecvf.com/content_cvpr_2016/html/Qi_Volumetric_and_Multi-View_CVPR_2016_paper.html
  2. Su, H., Maji, S., Kalogerakis, E., & Learned-Miller, E. (2015). Multi-view convolutional neural networks for 3d shape recognition. In Proceedings of the IEEE international conference on computer vision (pp. 945–953); https://www.cv-foundation.org/openaccess/content_iccv_2015/html/Su_Multi-View_Convolutional_Neural_ICCV_2015_paper.html
  3. Qi, C. R., Su, H., Mo, K., & Guibas, L. J. (2017). Pointnet: Deep learning on point sets for 3d classification and segmentation. In Proceedings of the IEEE conference on computer vision and pattern recognition (pp. 652–660); https://openaccess.thecvf.com/content_cvpr_2017/html/Qi_PointNet_Deep_Learning_CVPR_2017_paper.html
  4. Qi, C. R., Yi, L., Su, H., & Guibas, L. J. (2017). Pointnet++: Deep hierarchical feature learning on point sets in a metric space. Advances in neural information processing systems, 30; https://www.cs.toronto.edu/~bonner/courses/2022s/csc2547/papers/point_nets/pointnet++,_qi,_nips2017.pdf
  5. Wang, Y., Sun, Y., Liu, Z., Sarma, S. E., Bronstein, M. M., & Solomon, J. M. (2019). Dynamic graph cnn for learning on point clouds. Acm Transactions On Graphics (tog), 38(5), 1–12; https://dl.acm.org/doi/abs/10.1145/3326362
  6. Nikolov, I.; Madsen, C. (2020), "GGG – Rough or Noisy? Metrics for Noise Detection in SfM Reconstructions", Mendeley Data, V2; https://doi.org/10.17632/xtv5y29xvz.2
  7. Nikolov, I., & Madsen, C. (2020). Rough or Noisy? Metrics for Noise Estimation in SfM Reconstructions. Sensors, 20(19), 5725; https://doi.org/10.3390/s20195725
  8. Poux, F., & Billen, R. (2019). Voxel-based 3D point cloud semantic segmentation: unsupervised geometric and relationship featuring vs deep learning methods. ISPRS International Journal of Geo-Information. 8(5), 213; https://doi.org/10.3390/ijgi8050213
  9. Xu, J., Zhang, R., Dou, J., Zhu, Y., Sun, J., & Pu, S. (2021). Rpvnet: A deep and efficient range-point-voxel fusion network for lidar point cloud segmentation. In Proceedings of the IEEE/CVF International Conference on Computer Vision (pp. 16024–16033); https://openaccess.thecvf.com/content/ICCV2021/html/Xu_RPVNet_A_Deep_and_Efficient_Range-Point-Voxel_Fusion_Network_for_LiDAR_ICCV_2021_paper.html
  10. Haurum, J. B., Allahham, M. M., Lynge, M. S., Henriksen, K. S., Nikolov, I. A., & Moeslund, T. B. (2021, February). Sewer Defect Classification using Synthetic Point Clouds. In VISIGRAPP (5: VISAPP) (pp. 891–900); https://www.scitepress.org/Papers/2021/102079/102079.pdf
  11. Nikolov, I., & Madsen, C. B. (2021). Quantifying Wind Turbine Blade Surface Roughness using Sandpaper Grit Sizes: An Initial Exploration. In 16th International Conference on Computer Vision Theory and Application (pp. 801–808). SCITEPRESS Digital Library; https://doi.org/10.5220/0010283908010808

Related Articles