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;
}
|