DGtal
0.9.3beta

This tutorial shows you how to construct a polyhedron (polyhedral surface) from a digital object. It uses several ingredients of DGtal:
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).
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.
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
.
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.
Of course, class functors::SimpleThresholdForegroundPredicate is a wrapper around the image and only references the image.
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.
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
If everything is fine, you should have a nonnull digital 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.
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.
We need two different elements to decompose the boundary of our object into planar zones.
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
.
Ok, this is not an easy task, the smallest decomposition is even NPhard 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.
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.
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.
Then we put all the surfels into the priority queue Q
.
The best candidates for planar decomposition are now waiting in order in the queue.
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
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.
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.
This is the kind of result you should obtain.
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 leastsquare 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 leastsquare fit of a plane N.X=mu
to an arbitrary range of points [itB,itE).
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
.
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.
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.
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).
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:
This gives the following results for different digital objects and plane width.
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...