Voxelizer

The Voxelizer uses a robust, high-performance algorithm that converts a surface model in 3D into voxel data, as described in this paper, and can handle broken input fairly well. The reverse operation is performed by the SurfaceExtractor.

Voxelizer can be used with either Matrix, boost's multi_array, or with Blitz++'s Array through BlitzArray. Additionally, the Voxelizer depends on Geometry.

Alternatively, Voxelizer can output Octrees implemented as a 3-dimensional DTree. This representation usually requires much less memory and the voxelization is usually significantly faster.

It is also posisble to determine which fraction of the voxels is filled by the geometry. A seperate function exists that produces fractions in stead of a binary decision.

Note: MeshLab can be used to correct facet representations or to convert other datatypes to STL.

Interface

The Voxelizer can be partially parallelized by using OpenMP.

Note: the typename matrix_type is used to denote either one of the accepted matrix types.

  template <typename Tg, typename voxel_type>
  bool voxelize(const Geometry<Tg> &geometry,
		matrix_type<voxel_type, 3> &voxels,
		const value_type voxelSize = 1.0,
		const unsigned	 pad = 0u,
		const voxel_type inside = 1,
		const voxel_type outside = 0);                  
Voxelize a given geometry and store a representation in voxels in a 3D matrix.
The size of the geometry and the voxels, together with the padding, determine the size of the matrix. Parameters:
  • geometry The geometry to voxelize.
  • matrix The matrix that will contain the voxels. The matrix will be resized automatically.
  • voxelSize The size of voxels.
  • inValue The value given to elements of the matrix which are found to be inside the geometry.
  • outValue The value given to elements of the matrix which are found to be outside the geometry.
  • pad A number of extra cells around the outer dimensions of the geometry.
  template <typename Tg, typename Tf, typename voxel_type>
  bool fractionVoxelize(const Geometry<Tg>  &geometry,        
			Matrix_t<Tf, 3u, A> &voxelFractions,
			const double voxelSize,
			unsigned samples = 16u,
			unsigned pad = 0u)    
Voxelize a given geometry and store a representation in fractions of voxels in a 3D matrix.
The size of the geometry and the voxels, together with the padding, determine the size of the matrix. Parameters:
  • geometry The geometry to voxelize.
  • voxelFractions The matrix that will contain the fractions of the voxels. Type Tf must be a floating-point type. The matrix will be resized automatically.
  • voxelSize The size of voxels.
  • samples The number of samples along each dimension of the voxels.
  • pad A number of extra cells around the outer dimensions of the geometry.
  template <typename Tg, typename voxel_type>
  bool voxelize(const Geometry<Tg>  &geometry,
		DTree<voxel_type, 3u> &voxels,
		const double voxelSize,
		const voxel_type inside = 1, 
		const voxel_type outside = 0)    
Voxelize to an Octree.
  • geometry The geometry to voxelize.
  • voxels The Octree that will contain the voxels.
  • voxelSize The size of voxels.
  • inValue The value given to elements of the matrix which are found to be inside the geometry.
  • outValue The value given to elements of the matrix which are found to be outside the geometry.

Fraction Voxelizing Explained

Instead of producing a binary "in-or-out" decision, fraction voxelization will sample the space taken by a voxel and determine what fraction of it lies within the geometry. The complexity of the method is exponential in the precision of the result.

The output is in the range [0, 1], where 0 means "completely outside", and 1 means "completely inside". Values inbetween correspond to the fraction of the voxel that is filled by the geometry.

The advantage of the method is that it doesn't require much more memory than regular voxelization, nor complex analytical solutions.

A disadvantage of the method is the complexity. The function takes a parameter samples, which is the number of samples in each of the 3 dimensions X Y. Thus, if samples is 4, as used in the example section below, the algorithm will examine 4^3 sub-voxels per voxel, producing 64 possible values, which shows that the resulting precision is 6 bits. Each sample requies one "regular" voxelization, hence, the number of voxelizations equals samples^3; and the coresponding precision will be log_2(samples^3).
Conversely, the number of required number of voxelizations is exponential in the obtained precision. Therefore, anything more than 16 bit precision is probably not feasible.

Voxel-Center Distances

  template <typename Tg, typename voxel_type, typename Td>
  bool distances(const Geometry<Tg> &geometry,
	const DTree<voxel_type, 3u> &voxels,
	const double voxelSize,
	const std::vector<iPoint3D> &directions,
	std::map<std::size_t, /* DNode id */
		 std::vector<std::pair<int, Td> > > & distances,
	const voxel_type inside = 1)  
See text.
  template <template <typename Tm, std::size_t D, typename Aux> class Matrix_t,  
	  typename Tg, typename voxel_type, typename A, typename Td>
  bool distances(const Geometry<Tg> &geometry,
	const Matrix_t<voxel_type, 3u, A> &voxels,
	const double voxelSize,
	const std::vector<iPoint3D> &directions,
	std::map<iPoint3D, std::vector<std::pair<int, Td> > > & distances,
	const std::size_t pad = 0u, const voxel_type inside = 1)  
See text.

Given a 3D volume where 'geometry' is a facet representation and voxels an octree- or matrix representation of the same volume, the function distances calculates the distances from the center of each voxel to the nearest facet along the vectors in 'directions'.

Given two representations of the same 3D volume in facets and an octree by geometry and voxels respectively, calculate distances in voxel units from the center of voxels to the closest facet along the indicated directions. voxelSize must be the parameter passed to the voxelize() function that was used to transform geometry into voxels. The directions are vectors that indicate in which direction the distance to the facets is to be found. The direction vectors are understood to be contained within the voxel, and expressed relative to half the voxelsize, i.e. a vector [1,1,1] measures from the center to a corner. The direction vectors should therefor be composed out of values from the set [-1, 0, 1] only. The optional parameter inside is to indicate for which voxels the distances should be calculated. For matrices of voxels, the parameter pad should also be identical to the one passed to voxelize, or in both cases be omitted.

The calculated distances are in lattice units, i.e. a distance from one voxel center to one directly adjacent center is one. Note that the size of the lattice units varies in octrees. Due to diagonal directions, the distances are in the range [0 ... sqrt(3)]. If a distance is larger, than this is either to numerical error, or due to errors in the facet representation.

The data calculated will be placed in distances, which is a mapping from a location to a vector of direction-index-and-distance pairs. For voxel matrices, the location is indicated by an iPoint3D, for octrees, the location is indicated by the id of the node. The node can be retrieved from the tree using the retrieve() member function documented in DTree. For convenience, a member function index_trail() can be used to help deduce the location of the node in the tree.

The functions return true if for all directions a distance was found, and false if any were missed. Facet-representations are often inconsistent, or ill-formed, i.e. have missing facets, or have overlapping facets, et cetera. Although the voxelization procedure is somewhat robust to such errors in the input, the distances functions are not.

The current implementation is fairly basic, and not optimized for speed. It will iterate over the all voxels, all directions, and all facets, thus creating a 5-fold nested loop. Inside this loop is a matrix inversion. There is no parallelization at this point.

Note: The return value for the distances function is currently bogus.

Example

Simple Voxelization

Read an STL file and voxelize it:
#include <cvmlcpp/base/Matrix>
#include <cvmlcpp/volume/Geometry>
#include <cvmlcpp/volume/VolumeIO>
#include <cvmlcpp/volume/Voxelizer>

using namespace cvmlcpp;

int main(int argc, char **argv)
{
	Matrix<int, 3u> voxels;
	Geometry<float> geometry;

	readSTL(geometry, "cube.stl");

	voxelize(geometry, voxels);

	return 0;
}

Fraction Voxelization

Geometry<float> geometry;
readSTL(geometry, "cube.stl");

Matrix<float, 3u> fract;
unsigned pad = 1;
unsigned samples = 4;
double voxelSize = 0.1;
fractionVoxelize(geometry, fract, voxelSize, samples, pad);

double weight = 0.0;
for (unsigned x = 0u; x < 12; ++x)
for (unsigned y = 0u; y < 12; ++y)
for (unsigned z = 0u; z < 12; ++z)
	weight += fract[x][y][z];
assert(std::abs(weight - 1000.0) < 0.1);

Voxelization to Octree

Geometry<float> geometry;
readSTL(geometry, "cube.stl");

// Calculate a proper voxel size
const double nrVoxels = 1024;
double maxGeometrySize = 0.0;
for (unsigned d = 0; d < 3; ++d)
	maxGeometrySize = std::max(maxGeometrySize, double(geometry.max(d))-double(geometry.min(d)));
const double voxelSize = maxGeometrySize / nrVoxels;

DTree<short int, 3u> octree(0);
voxelize(geometry, octree, voxelSize);

std::cout << d << std::endl;

Distances

#include <cvmlcpp/base/StringTools>
#include <cvmlcpp/volume/DTree>
#include <cvmlcpp/volume/Geometry>
#include <cvmlcpp/volume/VolumeIO>
#include <cvmlcpp/volume/Voxelizer>

int main(int argc, char **argv)
{
	using namespace cvmlcpp;

	// Read facet data
	Geometry<float> g;
	if (argc >= 2)
		assert(readSTL(g, argv[1]));
	else
		assert(readSTL(g, "cube.stl"));

	// Determine voxelSize
	double voxelSize;
	if (argc == 3)
		voxelSize = std::atof(argv[2]);
	else
		voxelSize = 0.1;

	// Voxelize
	DTree<char, 3u> octree(0);
	assert(voxelize(g, octree, voxelSize));

	// D3Q15, I think
	std::vector<iPoint3D> directions;
	directions.push_back(iPoint3D( 0, 0, 0));

	directions.push_back(iPoint3D( 1, 0, 0));
	directions.push_back(iPoint3D( 0, 1, 0));
	directions.push_back(iPoint3D( 0, 0, 1));

	directions.push_back(iPoint3D(-1, 0, 0));
	directions.push_back(iPoint3D( 0,-1, 0));
	directions.push_back(iPoint3D( 0, 0,-1));

	directions.push_back(iPoint3D( 1, 1, 1));
	directions.push_back(iPoint3D( 1, 1,-1));
	directions.push_back(iPoint3D( 1,-1, 1));
	directions.push_back(iPoint3D(-1, 1, 1));
	directions.push_back(iPoint3D(-1,-1, 1));
	directions.push_back(iPoint3D(-1, 1,-1));
	directions.push_back(iPoint3D( 1,-1,-1));
	directions.push_back(iPoint3D(-1,-1,-1));

	// Calculate distances
	std::map<std::size_t, std::vector<std::pair<int, double> > > dists;
	assert(distances(g, octree, voxelSize, directions, dists));

	typedef std::map<std::size_t, std::vector<std::pair<int, double> > >::const_iterator map_const_iterator;
	for (map_const_iterator mi = dists.begin(); mi != dists.end(); ++mi)
	{
		const std::size_t node_id = mi->first;
		
		octree.retrieve(node_id);
		const std::vector<DTree<char, 3u>::index_t>
			trail = octree.retrieve(node_id).index_trail();
		std::cout << node_id << to_string(trail.begin(), trail.end()) << std::endl;
		typedef std::vector<std::pair<int, double> >::const_iterator vec_const_iterator;
		for (vec_const_iterator vi = mi->second.begin(); vi != mi->second.end(); ++vi)
		{
			const int direction_index = vi->first;
			const double distance     = vi->second;
			std::cout << "\t" << direction_index << " " 
				<< directions[direction_index].to_string() << " " 
				<< distance << std::endl;
			assert(direction_index > 0); // should skip (0, 0, 0) direction
		}
	}

	return 0;
}