DGtal  0.9.4beta
Tutorial: making a polyhedron from a digital object

Table of Contents

Author(s) of this documentation:
Jacques-Olivier Lachaud

This tutorial shows you how to construct a polyhedron (polyhedral surface) from a digital object. It uses several ingredients of DGtal:

Set up your project

We will write a new C++ source file (say makePolyhedron.cpp) that uses the DGtal library. You can follow for instance How to use DGtal in my own project ?. Then you can already start with the following skeleton that includes the main required headers (Qt will be used for displaying results).

#include <iostream>
#include <vector>
#include <set>
#include <map>
#include <queue>
#include <QtGui/qapplication.h>
#include "DGtal/base/Common.h"
#include "DGtal/helpers/StdDefs.h"

We will also play with some .vol file. There are a few in the directory DGtal/examples/samples. We pick below the Al Capone volume file. The program will also be parameterized by the threshold value on the volume file and the integers wp and wq, such that wp/wq is the tolerance used for decomposing the volume boundary into planes.

using namespace std;
using namespace DGtal;
using namespace Z3i;
int main( int argc, char** argv )
{
QApplication application(argc,argv);
string inputFilename = argc > 1 ? argv[ 1 ] : "DGtal/examples/samples/Al.100.vol";
int threshold = argc > 2 ? atoi( argv[ 2 ] ) : 0; // volume is (val > threshold)
int wp = argc > 3 ? atoi( argv[ 3 ] ) : 1;
int wq = argc > 4 ? atoi( argv[ 4 ] ) : 1;
...
return 0;
}

Reading the volume file and making the digital object

We store our volume as an image, i.e. a raster of value. We therefore use the class ImageContainerBySTLVector, parameterized by the usual Domain of Z3 (defined in namespace Z3i) and int.

typedef ImageContainerBySTLVector< Domain, int> Image;
Image image = VolReader<Image>::importVol(inputFilename);

However an image is not a digital object. We define a digital object implicitly as a Predicate : Point -> bool. Hence, we define the digital object as a threshold of the image with class functors::SimpleThresholdForegroundPredicate.

typedef functors::SimpleThresholdForegroundPredicate<Image> DigitalObject;
DigitalObject digitalObject( image, threshold );

Of course, class functors::SimpleThresholdForegroundPredicate is a wrapper around the image and only references the image.

Wrapping a surface around the boundary of the digital object

Up to now, the digital object is just a characteristic function, without any structure nor boundary. Digital surfaces are defined as inter voxel elemnts between face adjacent voxels, such that one is inside and the other is outside the digital object. It is thus natural to enrich the usual digital space \( \mathbb{Z}^d \) with lower dimensional cells. That is exactly the purpose of class KhalimskySpaceND. That kind of space allows to represent cells of arbitary dimension (see moduleCellularTopology). By convention, voxels are cells of maximal dimension d, surfels are d-1 dimensional cells, etc. Note also that surfels are oriented. It is because a surfel indicates not only a boundary element, but also in which direction is the inside.

KSpace ks; // defined in namespace Z3i, alias of KhalimskySpaceND<3,int>
ks.init( image.domain().lowerBound(), image.domain().upperBound(), true ); // can be checked if ok.

There are several kinds of digital surfaces (see also Digital surfaces). To address this diversity, the class DigitalSurface is templated by a model of CDigitalSurfaceContainer. Here we wish to define a digital surface as the boundary of a characteristic function (a model of concepts::CPointPredicate). We choose in our case the container ImplicitDigitalSurface, which does exactly this but around one boundary component only. You must provide one starting surfel for the implicitly defined surface, which indicates the component of the boundary that it is representing. Use Surfaces::findABel to get a first surfel element. Part of the code is given below

typedef KSpace::Surfel Surfel;
typedef ImplicitDigitalSurface< KSpace, DigitalObject > MyContainer;
typedef DigitalSurface< MyContainer > MyDigitalSurface;
MyContainer container( ks, digitalObject, surfAdj, start_surfel );
MyDigitalSurface digSurf( container );
trace.info() << "Digital surface has " << digSurf.size() << " surfels."
<< std::endl;

If everything is fine, you should have a non-null digital surface !

polyhedralizer-al-1.png
The original Al Capone digital object.

Visualising the surface

Just to check that everything is fine, we can visualize the surface with class Viewer3D, which is based on libQGLViewer (which in turn uses Qt and OpenGL). The principle is to instantiate the viewer and then to use output streams to feed the viewer with cells, colors, etc. The following piece of code displays the first surfel of your surface.

typedef Viewer3D<Space,KSpace> MyViewer3D;
MyViewer3D viewer( ks );
viewer.show();
viewer << CustomColors3D( Color::Black, Color::White );
viewer << *( digSurf.begin() ); // surfel pointed by the begin iterator.
viewer << MyViewer3D::updateDisplay;
application.exec();

Update this code to display the whole surface, just by iterating on the range provided by the digital surface. To speed up a little bit the display, you may add the following line (after show() ), which forces a simplified display for surfels.

viewer << SetMode3D( start_surfel.className(), "Basic" );

Recognizing planar zones on the object

We need two different elements to decompose the boundary of our object into planar zones.

typedef ChordGenericNaivePlaneComputer<Z3,Z3::Point, DGtal::int64_t> NaivePlaneComputer;
typedef BreadthFirstVisitor<MyDigitalSurface> Visitor;

A breadth first visitor, given a starting vertex on a graph, visits in sequence vertices at increasing (topological) distance. The current Visitor::Node holds a pair <Vertex,Size> which represents the current vertex in the traversal and its distance to the starting vertex. You may choose at each step whether you can continue normally to the next vertex (BreadthFirstVisitor::expand) or if you continue to the next vertex but you block the expansion of the visitor through the current vertex (BreadthFirstVisitor::ignore). We use the visitor layer by layer (i.e. vertices with same distance). Here is a skeleton code for extracting a "round" and planar neighborhood aroud surfel start_surfel. We call this plane the round plane around start_surfel.

NaivePlaneComputer planeComputer;
std::vector<Point> layer;
planeComputer.init( wp, wq );
// The visitor takes care of all the breadth-first traversal.
Visitor visitor( digSurf, start_surfel );
Visitor::Size currentSize = visitor.current().second;
layer.clear();
while ( ! visitor.finished() )
{
Visitor::Node node = visitor.current();
Surfel v = node.first;
int axis = ks.sOrthDir( v );
Point p = ks.sCoords( ks.sDirectIncident( v, axis ) ); // point interior to surfel
if ( node.second != currentSize )
{ // tries to extend the plane with the whole layer.
bool isExtended = planeComputer.extend( layer.begin(), layer.end() );
if ( isExtended )
{ // ... do your stuff
...
layer.clear();
currentSize = node.second;
}
else
break; // there is no plane that includes also the next layer.
}
layer.push_back( p );
visitor.expand();
}
trace.info() << "Vertex " << start_surfel << " has plane " << planeComputer.primitive() << std::endl;
Note
A surfel has not really integer coordinates that may be used readily for a digital plane recognition algorithm. Our choice here is to see the surfel as if it had the coordinates of its interior voxel. You could use other embeddings: (1) use the four pointels of each surfel, (2) use the doubled half-integer centroid coordinates of the surfel.

Finding the best decomposition into planes

Ok, this is not an easy task, the smallest decomposition is even NP-hard in the general case. We use the following heuristic, which is reasonnably costly. First, we associate an integer to each surfel, and we initialize it to zero.

// Initialisation
std::map<Surfel,unsigned int> v2size;
for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
v2size[ *it ] = 0;

Secondly, for each surfel, we compute its round plane. Then we increment the counter of each surfel that belongs to this round plane. You must reuse the code computing the round plane (Recognizing planar zones on the object), but you must not only store the points of the layer bu also the surfels of the layer. This gives something like that.

int j = 0;
int nb = digSurf.size();
std::vector<Point> layer;
std::vector<Surfel> layer_surfel;
for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
{
if ( ( (++j) % 50 == 0 ) || ( j == nb ) ) trace.progressBar( j, nb );
Surfel v = *it;
Visitor visitor( digSurf, v );
layer.clear();
layer_surfel.clear();
// compute round plane here and update v2size accordingly.
...
}

Note that the execution of this part may take some time for big objects. Once this is done, the surfels that belongs to a lot of planes have a high value in v2size. And generally, surfels at the center of big planar regions have a high value. We will thus put the surfels in a priority queue, such that the first surfels to pop out are the one with high value in v2size: they are indeed excellent candidates for finding a good decomposition. To do this, you can use the following generic class for storing pairs and ordering them according to the second member.

template <typename T1, typename T2>
struct PairSorted2nd
{
typedef PairSorted2nd<T1,T2> Self;
PairSorted2nd( const T1& t1, const T2& t2 ) : first( t1 ), second( t2 ) {}
bool operator<( const Self& other ) const { return second < other.second; }
T1 first;
T2 second;
};

Then we put all the surfels into the priority queue Q.

typedef PairSorted2nd<Surfel,int> WeightedSurfel;
std::priority_queue<WeightedSurfel> Q;
for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
Q.push( WeightedSurfel( *it, v2size[ *it ] ) );

The best candidates for planar decomposition are now waiting in order in the queue.

Decomposing the object boundary into planes

The idea is to pop the surfels out of the queue in order. If the popped surfel has not been marked before, then we build its round plane. All the unmarked surfels in the round plane are associated to this round plane, and are in turn marked. We process until the queue is empty.

In order to store the round planes (and other information like its color for display), you shoud define a structure like

template <typename T1, typename T2, typename T3>
struct Triple
{
T1 first;
T2 second;
T3 third;
Triple( T1 t1 = T1(), T2 t2 = T2(), T3 t3 = T3() )
: first( t1 ), second( t2 ), third( t3 )
{}
};
// ...
typedef Triple< NaivePlaneComputer, Color, std::pair<RealVector,double> > RoundPlane;

Since a round plane will be shared by many surfels, it is more convenient to allocate them dynamically. A CountedPtr could also be used, but this is enough for a tutorial.

std::set<Surfel> markedSurfels;
std::vector<RoundPlane*> roundPlanes;
std::map<Surfel,RoundPlane*> v2plane;
j = 0;
while ( ! Q.empty() )
{
if ( ( (++j) % 50 == 0 ) || ( j == nb ) ) trace.progressBar( j, nb );
Surfel v = Q.top().first;
Q.pop();
if ( markedSurfels.find( v ) != markedSurfels.end() ) // already in set
continue; // process to next vertex
RoundPlane* ptrRoundPlane = new RoundPlane;
roundPlanes.push_back( ptrRoundPlane ); // to delete them afterwards.
v2plane[ v ] = ptrRoundPlane;
ptrRoundPlane->init( wp, wq );
ptrRoundPlane->second = Color( rand() % 192 + 64, rand() % 192 + 64, rand() % 192 + 64, 255 );
// same extraction of round plane as above
// Do not forget to associate ptrRoundPlane to each unmarked visited surfels of the round plane.
// Do not forget to mark all surfels of the round plane.
}

At the end of your program, you should of course call delete on each allocated round plane. Since each round plane has a (random) color, we may update our visualisation to display for each surfel the color of its round plane.

MyViewer3D viewer( ks );
viewer.show();
viewer << SetMode3D( start_surfel.className(), "Basic" );
for ( std::map<Surfel,RoundPlane*>::const_iterator
it = v2plane.begin(), itE = v2plane.end();
it != itE; ++it )
{
Surfel v = it->first;
RoundPlane* rplane = it->second;
viewer << CustomColors3D( rplane->second, rplane->second );
viewer << v;
}
viewer << MyViewer3D::updateDisplay;
application.exec();

This is the kind of result you should obtain.

polyhedralizer-al-round.png
The round planes of Al Capone digital object, for width=3/2.

Making a polyhedral surface

Now that meaningul planes have been extracted, there are many ways to build a polyhedral surface. We adopt a simple approach, which gives correct results. First of all, we must associate a real plane to each (digital) round plane. The method ChordGenericNaivePlaneComputer::primitive returns the parallel strip that encloses the digital points. We could for instance take the middle plane of the strip as real plane. Another way is to fit by least-square a plane to the points of the round plane. For our purpose, it gives nicer results. That's what we do here. The function below realizes the least-square fit of a plane N.X=mu to an arbitrary range of points [itB,itE).

// Least-Square Fit of a plane N.x=mu on points [itB,itE). Returns mu.
template <typename RealVector,
typename ConstIterator>
double LSF( RealVector& N, ConstIterator itB, ConstIterator itE )
{
typedef typename RealVector::Component Component;
typedef SimpleMatrix<Component,3,3> Matrix;
Matrix A; A.clear();
unsigned int nb = 0;
RealVector G = RealVector::zero; // center of gravity
for ( ConstIterator it = itB; it != itE; ++it )
{
G += RealVector( (*it)[ 0 ], (*it)[ 1 ], (*it)[ 2 ] );
++nb;
}
G /= nb;
for ( ConstIterator it = itB; it != itE; ++it )
{
RealVector p( (*it)[ 0 ], (*it)[ 1 ], (*it)[ 2 ] );
p -= G;
for ( Dimension i = 0; i < 3; ++i )
for ( Dimension j = 0; j < 3; ++j )
A.setComponent( i, j, A( i, j ) + p[ i ] * p[ j ] );
}
// A is Gram matrix. We look for V such that V^t A V / |V|^2 is
// minimal. It is thus the first eigenvalue.
Matrix V;
RealVector values;
N = V.column( 0 ); // first eigenvector;
double mu = 0.0;
for ( ConstIterator it = itB; it != itE; ++it )
mu += N.dot( *it );
return mu/(double)nb;
}

You must now compute the LSF plane for each round plane. This is easily done by iterating over the vector roundPlanes. We store the normal vector N and the offset mu into the RoundPlane.

for ( std::vector<RoundPlane*>::iterator
it = roundPlanes.begin(), itE = roundPlanes.end();
it != itE; ++it )
{
NaivePlaneComputer& computer = (*it)->first;
RealVector normal;
double mu = LSF( normal, computer.begin(), computer.end() );
(*it)->third = make_pair( normal, mu );
}

Then each surfel is first projected onto the real plane associated to its round plane. In a second sequence, we average the coordinates of the surfels with the coordinates of their neighbors. The effect is that surfels inside a round plane are still projected onto their plane while surfels just on the rim of a round plane are averaged with the neighboring plane.

std::map<Surfel, RealPoint> coordinates;
for ( std::map<Surfel,RoundPlane*>::const_iterator
it = v2plane.begin(), itE = v2plane.end();
it != itE; ++it )
{
Surfel v = it->first;
RoundPlane* rplane = it->second;
Point p = ks.sKCoords( v );
RealPoint rp( (double)p[ 0 ]/2.0, (double)p[ 1 ]/2.0, (double)p[ 2 ]/2.0 );
double mu = rplane->third.second;
RealVector normal = rplane->third.first;
double lambda = mu - rp.dot( normal );
coordinates[ v ] = rp + lambda*normal;
}
typedef std::vector<Surfel> SurfelRange;
std::map<Surfel, RealPoint> new_coordinates;
for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
{
Surfel s = *it;
SurfelRange neighbors;
std::back_insert_iterator<SurfelRange> writeIt = std::back_inserter( neighbors );
digSurf.writeNeighbors( writeIt, *it );
RealPoint x = RealPoint::zero;
for ( SurfelRange::const_iterator its = neighbors.begin(), itsE = neighbors.end();
its != itsE; ++its )
x += coordinates[ *its ];
new_coordinates[ s ] = x / neighbors.size();
}

It remains to build the polyhedral surface approximating the digital boundary. We use the class Mesh (from package shapes), which let us construct a polyhedral surface by simply adding vertices and then adding faces as vectors of vertex indices. Class Mesh is devoted to i/o. You can output it on a viewer to display it or use the stream operator to export it as some OBJ or OFF file.

typedef unsigned int Number;
typedef Mesh<RealPoint> MyMesh;
typedef MyMesh::MeshFace MeshFace;
typedef MyDigitalSurface::FaceSet FaceSet;
typedef MyDigitalSurface::VertexRange VertexRange;
std::map<Surfel, Number> index; // Numbers all vertices.
Number nbv = 0;
MyMesh polyhedron( true );
// Insert all projected surfels as vertices of the polyhedral surface.
for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
{
polyhedron.addVertex( new_coordinates[ *it ] );
index[ *it ] = nbv++;
}

Since surfels are mapped as vertices into the mesh, we need to build faces on the mesh such that their vertices are surfels. It is in fact a kind of dual surface to the digitized boundary. The class DigitalSurface is also able to list the polygonal faces that joins surfels, through method DigitalSurface::allClosedFaces. These faces form umbrellas of surfels around some pointel (see the works of Fran├žon).

FaceSet faces = digSurf.allClosedFaces();

It remains just to iterate over all the faces, and for each face, to use method Mesh::addFace to add a Mesh::MeshFace with the correct vertex indices. Have a look at method DigitalSurface::verticesAroundFace in order to get the surfels around a given face.

Once the polyhedron is build, you may view it and export it as follows:

viewer.clear();
viewer.show();
viewer << polyhedron;
viewer << MyViewer3D::updateDisplay;
application.exec();
bool isOK = polyhedron >> "test.off";
bool isOK2 = polyhedron >> "test.obj";

This gives the following results for different digital objects and plane width.

polyhedral-al-finished-w3.png
The polyhedral surface approaching Al Capone digital object, for width=3/1.
polyhedral-bunny-finished-w2.png
The polyhedral surface approaching bunny-256 dataset, for width=2/1.
polyhedral-ssphere-ctmv.png
A nicer viewer for the polyhedral surface approaching sharp-sphere-129, for width=2/1.

Going further

It is clear that a better polyhedralizer would not keep vertices inside a planar zone. It is not too difficult (although not straightforward) to go further along this path. It remains to identify vertices as surfels touching several round planes, and then to build the faces manually. Have fun...