Table of Contents
- Set up your project
- Reading the volume file and making the digital object
- Wrapping a surface around the boundary of the digital object
- Visualising the surface
- Recognizing planar zones on the object
- Finding the best decomposition into planes
- Decomposing the object boundary into planes
- Making a polyhedral surface
- Going further
- 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:
- readers: to read an ".obj" file into memory (class VolReader)
- images: to store the ".obj" file into memory (class ImageContainerBySTLVector)
- predicate: to represent the characteristic function of the digital object (class functors::SimpleThresholdForegroundPredicate)
- cellular space: the cellular digital space in which the boundary of the digital object is well defined (class KhalimskySpaceND)
- digital surface: to connect all boundary elements of the digital object as a graph (class DigitalSurface)
- graph visitor: to define neighborhoods on the digital object boundary (class BreadthFirstVisitor)
- digital plane computer: to recognize points of the boundary that belongs to the same plane (class ChordGenericNaivePlaneComputer)
- umbrellas: to decompose the digital surface into cells
- viewers: to see intermediate and final results of our algorithm (class Viewer3D)
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).
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.
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.
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.
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.
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 non-null digital surface !
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.
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.
Recognizing planar zones on the object
We need two different elements to decompose the boundary of our object into planar zones.
- The first one is a way to determine if a set of points is planar. We will use the class ChordGenericNaivePlaneComputer to carry out this task, but you may also have a look at page Plane recognition and plane width computation to have a broader view on this topic. Such an object recognizes arbitrary pieces of planes with a given axis width. Hence it recognizes naive planes when the axis width is 1.
- The second element is a way to visit the points onto the object boundary in a reasonnable way, so that planar pieces are not too elongated. We adopt a simple approach to produce as "round" as possible planar pieces. The digital surface is a model of graph (concept CUndirectedSimpleGraph) and we may perform several breadth-first traversal onto it with the class BreadthFirstVisitor. Again, page Basic graph definitions and concepts gives a more complete overview.
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.
- 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.
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.
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
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.
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).
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.
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...
