[Introduction] [Download] [Usage] [Disclaimer and Copyright]

A crucial part of OOF3D is the calculation of the volume of intersection of a finite element with a set of voxels (three dimensional pixels). This calculation must be fast, and it must be robust, so that small numerical errors in the position of an intersection point do not lead to large errors in the result. Because it is easy to write code that is not robust, we are here providing the OOF3D code in a separate module that can be easily used elsewhere.

These routines are called VSB, which stands for "Voxel Set Boundary" because "Intersection of Voxel Set Boundary with Convex Polyhedron", or IOVSBWCP, is too long. They are based on Powell and Abel's r3d algorithm.

This page provides information on how to download and use VSB. It doesn't explain the philosophy behind it. For details on how VSB works, see our paper and the comments in the code.

VSB requires a C++ compiler supporting the C++11 standard. It is a template library that must be compiled with your C++ program. It is not a stand-alone program.

Click to download vsb-1.0.2.tgz (35783 bytes, md5 checksum 1430d792087b3b68753005c092f8cb02).

Unpack the file with tar -xf vsb-1.0.2.tgz. This will create a directory called VSB containing C++ files.

Quick Summary

Include vsb.h. Create an instance of the VoxelSetBdy class. Its arguments include an image and an voxel value. Then call VoxelSetBdy::clippedVolume(). The arguments to clippedVolume are the sides of a polyhedron, and its return value is the volume of the intersection of the polyhedron with all voxels in the image with the given voxel value.

Details

At a minimum, you must include VSB/vsb.h in your C++ code, call the methods described below, and compile and link to vsb.C. (Rename vsb.C to vsb.cxx if that's your convention.)

Almost all of the classes and functions in VSB are templates. The template arguments are C++ classes that you must provide, representing coordinates, images, voxel values, etc. If you're writing a program that involves polyhedra and voxels you probably already have these classes. VSB assumes that these classes have certain minimal abilities, spelled out below and in VSB/vsb.h.

Template Arguments

These are the types that you must provide as template arguments to the VSB classes and functions.

• COORD: A point in 3-space, with double-valued components. The following member and non-member functions must be supported:

  • Coord(double x, double y, double z) constructor.
  • double Coord::operator[](int i) const retrieve a component (i = 0, 1, or 2).
  • double& Coord::operator[](int) retrieve or set a component.
  • COORD& COORD::operator+=(const COORD&) addition in place.
  • COORD operator+(const COORD&, const COORD&) vector addition.
  • COORD operator*(const COORD&, double) scalar multiplication.
  • COORD operator/(const COORD&, double) scalar division.
  • double dot(const COORD&, const COORD&) the three dimensional dot product.
  • COORD cross(const COORD&, const COORD&) the three dimensional cross product.

• ICOORD: Just like COORD, but with integer components. It must support the following functions:

  • ICOORD(int x, int y, int z) constructor.
  • int ICOORD::operator[](int i) const retrieve a component (i = 0, 1, or 2).
  • int& ICOORD::operator[](int) retrieve or set a component.
  • bool operator==(const ICOORD&, const ICOORD&) equality comparison.
  • bool operator!=(const ICOORD&, const ICOORD&) inequality comparison.
  • ICOORD operator+(const ICOORD&, const ICOORD&) vector addition.
  • ICOORD operator-(const ICOORD&, const ICOORD&) vector subtraction.

• IMAGEVAL: The value of a voxel (eg, a color or a segmentation index). It must have an equality operator:

  • bool operator==(const IMAGEVAL&, const IMAGEVAL&)

• IMAGE: A 3D array of voxels. It must allow a voxel value to be extracted, given its position as an ICOORD:

  const IMAGEVAL &IMAGE::operator[](const ICOORD&) const

The VoxelSetBdy class

VoxelSetBdy is the main class in VSB. It represents a set of voxels, which is stored in terms of the edges of boundary of the set. Its clippedVolume method computes the volume remaining after the voxel set is clipped by a set of planes. Clipping by all of the planes of a convex polyhedron produces the intersection of the voxel set and the polyhedron.

One VoxelSetBdy object can be clipped multiple times by multiple sets of planes. Clipping does not change the VoxelSetBdy.

The constructor is:

VoxelSetBdy<COORD, ICOORD, IMAGE, IMAGEVAL>(
	    const IMAGE &image,
	    const IMAGEVAL &imageVal,
	    double voxelVolume,
	    const std::vector<ICRectPrism<ICOORD>> &subregions)

• image: an IMAGE object containing the voxels
• imageval: an IMAGEVAL. The voxel set consists of all voxels with this value.
• voxelVolume: the volume of a single voxel. Set it to 1.0 if you want results in terms of the number of voxels.
• subregions: If a polyhedron is much smaller than the image, it's inefficient for VSB to spend time clipping voxels that are far from the polyhedron. It's better to divide the image into subregions and only examine voxels inside the nearby subregions. subregions is a vector of ICRectPrism objects representing rectilinear regions which must fill the image with no overlaps. ICRectPrism is defined in VSB/cprism.h

VoxelSetBdy methods:

•  double VoxelSetBdy::volume() const

  computes the volume of the voxel set.

• double VoxelSetBdy::clippedVolume(
  	      const CRectPrism<COORD> &bbox,
  	      const std::vector<VoxelSetBdy::Plane> &planes
  	      ) const 

  computes the volume of the voxel set after clipping by the given planes. bbox is the rectangular bounding box containing the region bounded by the planes. If the bounding box is too large, the computation may take longer than necessary, but if it's too small, the results may be incorrect. VoxelSetBdy::Plane is a typedef for VSBPlane<COORD>. CRectPrism is defined in VSB/cprism.h and VSBPlane is defined in VSB/cplane.h.

The VSBPlane class

VSBPlane is a templated class, defined in VSB/cplane.h that represents an oriented plane. The single template parameter is a COORD. The constructor is

VSBPlane(const COORD &normal, double offset)

The ICRectPrism and CRectPrism classes

ICRectPrism and CRectPrism, defined in VSB/cplane.h, represent rectangular prisms (3D rectangles). The difference between them is that the corners of an ICRectPrism are at integer coordinates, while the corners of a CRectPrism are floating point numbers. Both are template classes. The template parameter for ICRectPrism is ICOORD. The template parameter for CRectPrism is COORD.

The constructors are:

• CRectPrism() creates an unitialized prism. You should call another method to set its bounds before using it.
• CRectPrism(const COORD& a, const COORD& b) and CRectPrism(const COORD* a, const COORD* b) both create a prism with opposite corners at a and b.
• CRectPrism(const std::vector<COORD>& points) creates a prism containing all of the given points.

Other useful methods for setting up a prism are:

• CRectPrism::swallow(const POINT& point) expands a prism so that it includes the given point, which may be either a COORD or an ICOORD (or any other type that acts like a COORD or ICOORD).
• CRectPrism::swallowPrism(const PRISM& other) expands a prism so that it includes the given other prism, which may be a CRectPrism or ICRectPrism.
• CRectPrism::restrict(const PRISM& other) shrinks the prism to its intersection with the given other prism, which may be a CRectPrism or ICRectPrism.

See VSB/cprism.h for other methods, including retrieving corners and bounds, determining if a prism contains a given point, and whether or not two prisms intersect.

This software was produced by NIST, an agency of the U.S. government, and by statute is not subject to copyright in the United States. Recipients of this software assume all responsibilities associated with its operation, modification and maintenance. However, to facilitate maintenance we ask that before distributing modified versions of this software, you first contact the authors at oof_manager@list.nist.gov.