DGtal  0.9.4beta
Digital surfaces

Table of Contents

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

Part of the Topology package.

This part of the manual describes how to define digital surfaces, closed or open. A lot of the ideas, concepts, algorithms, documentation and code is a backport from ImaGene. [39]

The following programs are related to this documentation: volToOFF.cpp, volMarchingCubes.cpp, volScanBoundary.cpp, volTrackBoundary.cpp, frontierAndBoundary.cpp, volBreadthFirstTraversal.cpp, trackImplicitPolynomialSurfaceToOFF

Introduction to digital surfaces

Possible definitions for digital surfaces

Although continuous surfaces are well defined as n-manifolds, a digital or discrete analog of surface is more tricky to define. Several methods have been explored to define consistent digital surfaces. We mention three approaches here.

We focus here on the last method to define digital surfaces.

Digital surface as a set of n-1-cells

Formally, the elements of the digital space \( Z^n \) are called spels, and often pixels in 2D and voxels in 3D. A surfel is a couple (u,v) of face-adjacent voxels. A digital surface is a set of surfels. It is obvious that a spel is a n-cell in the cellular grid decomposition of the space, while a surfel is clearly some oriented n-1-cell which is incident to the two n-cells u and v (see Cellular grid space and topology, unoriented and oriented cells, incidence).

Let s be some surfel. It is incident to two oriented voxels. By convention, its interior pixel is the one that is positively oriented.

We assume from now on that you have instantiated some cellular space K of type KSpace (see dgtal_ctopo_sec4), for instance with the following lines:

KSpace K;
Point low( -10, -10, -10 );
Point high( 10, 10, 10 );
bool space_ok = K.init( low, high, true );

A surfel is an oriented \(n-1\)-cell, i.e. some SCell. Surfel may be obtained from spels by incidence relation. Reciprocally, you can use incidence to get spels.

typedef KSpace::SCell Surfel; // or typedef KSpace::Surfel Surfel;
typedef KSpace::SCell Spel; // or typedef KSpace::Surfel Surfel;
Spel v = K.sSpel( Point( 0, 0, 0 ), POS ); // +v
Surfel sx = K.sIncident( v, 0, true ); // surfel further along x
Surfel sy = K.sIncident( v, 1, true ); // surfel further along y
Surfel sz = K.sIncident( v, 2, true ); // surfel further along z
Spel qx = K.sDirectIncident( s, 0 ); // positive coboundary of s
Spel qy = K.sDirectIncident( s, 1 ); // positive coboundary of s
Spel qz = K.sDirectIncident( s, 2 ); // positive coboundary of s
ASSERT( v == qx && v == qy && v == qz ); // same as qx, qy, qz

The direct orientation to some cell s is the one that gives a positively oriented cell in the boundary or coboundary of s.

You may now for instance define the digital surface that lies in the boundary of some digital shape \( S \subset Z^n \) as the set of oriented surfels between spels of S and spels not in S. Algebraically, S is the formal of its positively oriented spels, and its boundary is obtained by applying the boundary operator on S.

digital-surface-BdryOp-1s.png
Using the boundary operator to compute the boundary of a digital shape S.
digital-surface-BdryOp-2s.png
The boundary operator is linear with cells, thus we compute spel by spel.
digital-surface-BdryOp-3s.png
The boundary operator is linear with cells, thus we compute spel by spel.
digital-surface-BdryOp-4s.png
Opposite cells cancel each other.
digital-surface-BdryOp-5s.png
This is the resulting set of surfels (in blue and magenta.

Defining a digital surface as a set of surfels is not enough for geometry. Indeed we need to relate surfels together so as to have a topology on the digital surface. The first step is to transform the digital surface into a graph.

Digital surface as a graph: adding adjacencies between surfels

The idea here is to connect surfels that share some \(n-2\)-cells. The obtained adjacency relations are called bel adjacencies in the terminology of Herman, Udupa and others. Generally an \(n-2\)-cell is shared by at most two \(n-1\)-cells, except in "cross configuration", symbolized below:

   O | X    X: interior spels               O a X
   - + -    O: exterior spels               b + c
   X | O    - and |: surfels a,b,c,d        X d O
            +: n-2-cell 

A bel adjacency makes a deterministic choice to connect b to d and a to c when interior, and b to a and c to d when exterior. This choice has to be made along each possible pair of directions when going nD. In DGtal, this is encoded with the class SurfelAdjacency.

digital-surface-IntAdjacency.png
Interior adjacency on digital surfaces in 2D

The following snippet shows the interior surfel adjacency (i.e. (6,18)-surfaces).

typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
MySurfelAdjacency surfAdj( intAdjacency ); // interior in all directions.

Once the adjacency has been chosen, it is possible to determine what are the adjacent surfels to a given surfel. More precisely, any surfel (or \(n-1\)-cell) has exactly \(2n-2\) adjacent surfels (for closed surfaces). More precisely, it has exactly 2 adjacent surfels along directions different from the orthogonal direction of the surfel.

digital-surface-SurfaceTracking2.png
Any surfel in 3D has 4 adjacent surfels.

In fact, we can be even more precise. We can use orientation to orient the adjacencies consistently. Given two surfels sharing an \(n-2\)-cell, this cell is positively oriented in the boundary of one surfel and negatively oriented in the boundary of the other. We have thus that there are \(n-1\) adjacent cells in the direct orientation (positive) and \(n-1\) adjacent cells in the indirect orientation (negative). The direct adjacent surfels look like:

digital-surface-SurfaceTracking.png
Any surfel in 3D has 2 adjacent surfels in the direct orientation.

The class SurfelNeighborhood is the base class for computing adjacent surfels. You may use it as follow, but generally this class is hidden for you since you will have more high-level classes to move on surfaces.

// K is a KSpace, surfAdj is a SurfelAdjacency, surfel is some SCell.
SurfelNeighborhood<KSpace> sNeigh;
sNeigh.init( &K, &SAdj, surfel );
trace.info() << "surfel=" << surfel << endl;
trace.info() << "follower1(+)=" << sNeigh.follower1( 1, true ) << endl;
trace.info() << "follower2(+)=" << sNeigh.follower2( 1, true ) << endl;
trace.info() << "follower3(+)=" << sNeigh.follower3( 1, true ) << endl;
trace.info() << "follower1(-)=" << sNeigh.follower1( 1, false ) << endl;
trace.info() << "follower2(-)=" << sNeigh.follower2( 1, false ) << endl;
trace.info() << "follower3(-)=" << sNeigh.follower3( 1, false ) << endl;

Generally, methods SurfelNeighborhood::getAdjacentOnSpelSet, SurfelNeighborhood::getAdjacentOnDigitalSet, SurfelNeighborhood::getAdjacentOnPointPredicate, SurfelNeighborhood::getAdjacentOnSurfelPredicate are used to find adjacent surfels.

Since we have defined adjacencies between surfels on a digital surface, the digital surface forms a graph (at least in theory).

Tracking digital surfaces

This section shows how to extract the boundary of a digital set, given some predicate telling whether we are inside or outside the surface.

Constructing digital surfaces by scanning

Given a domain and a predicate telling whether we are inside or outside the object of interest, it is easy to determine the set of surfels by a simple scanning of the space. This is done for you by static methods Surfaces::uMakeBoundary and Surfaces::sMakeBoundary.

The following snippet shows how to get a set of surfels that is the boundary of some predicate on point by a simple scan of the whole domain (see volScanBoundary.cpp).

trace.beginBlock( "Extracting boundary by scanning the space. " );
KSpace::SCellSet boundary;
Surfaces<KSpace>::sMakeBoundary( boundary, ks, set3d,
image.domain().lowerBound(),
image.domain().upperBound() );

Constructing digital surfaces by tracking

In many circumstances, it is better to use the above mentioned graph structure of digital surfaces. For instance it may be used to find the surface just by searching it by adjacencies. This process is called tracking. This is done for you by static method Surfaces::trackBoundary.

The following snippet shows how to get a set of surfels that is the boundary of some predicate on point given only a starting surfel and by tracking (see volTrackBoundary.cpp).

trace.beginBlock( "Extracting boundary by tracking from an initial bel." );
KSpace::SCellSet boundary;
SCell bel = Surfaces<KSpace>::findABel( ks, set3d, 100000 );
surfAdj,
set3d, bel );

On the lobser.vol volume, volScanBoundary.cpp extracts 155068 surfels in 3866ms, while volTrackBoundary.cpp extracts 148364 surfels in 351ms.

# Commands
$ ./examples/topology/volScanBoundary ../examples/samples/lobster.vol 50 255
$ ./examples/topology/volTrackBoundary ../examples/samples/lobster.vol 50 255
volTrackBoundary-lobster.png
Digital surface that is the boundary of a (6,18)-connected component in image lobster.vol, extracted by tracking from an initial surfel in 351ms.

In the case you know beforehands that your surface is closed (i.e. without boundary), you should use preferentially Surfaces::trackClosedBoundary. This technique only follows direct adjacent surfels and extracts the whole surface when it is closed. This process can be illustrated as follows:

suivi-parcours-largeur.png
Tracking a digital surface by following adjacencies.
suivi-artzy.png
Tracking a digital surface by following only direct adjacencies.

Other constructions by tracking: 2D contours, surface connected components, 2D slices in 3D surface

You can have a look at page Helpers for digital surfaces to see how to construct directly a 2D contour in \( Z^2 \), how to get the set of surface components, how to get 2D contours as slices onto a 3D surface.

You may also use the template class DigitalSurface2DSlice, which, given a tracker on a digital surface, computes and store a 2d slice on the associated surface. You may after visit the slice with iterators or circulators.

// Tracker is a type of tracker on some digital surface.
typedef DigitalSurface2DSlice<Tracker> MySlice;
Tracker* ptrTracker = ... ; // some pointer on a tracker.
MySlice slicex( ptrTracker, 0 ); // slice containing x-axis
MySlice slicey( ptrTracker, 1 ); // slice containing y-axis
MySlice slicez( ptrTracker, 2 ); // slice containing z-axis
// One of them is invalid in 3D.
See also
digitalSurfaceSlice.cpp

High-level classes for digital surfaces

Digital surfaces arise in many different contexts:

Since there are so many digital surfaces, it is necessary to provide a mechanism to handle them generically. The class DigitalSurface will be the common proxy to hide models of concepts::CDigitalSurfaceContainer.

A common architecture for digital surfaces

Since there may be several types of digital surfaces, there are several container classes for them. All of them follow the concept concepts::CDigitalSurfaceContainer. Essentially, a model of this class should provide methods begin() and end() to visit all the surfels, and a Tracker which allows to move by adjacencies on the surface. A Tracker should be a model of concepts::CDigitalSurfaceTracker. The architecture is sumed up below:

diag-digital-surface-1.png
Class architecture of digital surfaces.

Models of digital surface containers

For now, the implemented digital surface container are (realization of concepts::CDigitalSurfaceContainer):

Depending of what is your digital surface, you should choose your container accordingly:

Light versions of containers should be preferred in mostly two cases:

  1. The digital surface is big and you do not wish to keep it in memory. This is probably the case when tracking some implicit function and outputing it in some stream.
  2. The digital surface is likely to evolve (i.e. the predicate do not give the same response at a point/surfel depending on when it is called).

Relating a digital surface to some container

The following snippet shows how to create a digital surface associated to a LightImplicitDigitalSurface. The LightImplicitDigitalSurface is connected to a shape described by a predicate on point (a simple digital set). The full code is in volBreadthFirstTraversal.cpp.

trace.beginBlock( "Set up digital surface." );
typedef LightImplicitDigitalSurface<KSpace, DigitalSet >
MyDigitalSurfaceContainer;
typedef DigitalSurface<MyDigitalSurfaceContainer> MyDigitalSurface;
SCell bel = Surfaces<KSpace>::findABel( ks, set3d, 100000 );
MyDigitalSurfaceContainer* ptrSurfContainer =
new MyDigitalSurfaceContainer( ks, set3d, surfAdj, bel );
MyDigitalSurface digSurf( ptrSurfContainer ); // acquired

Once this is done, the object digSurf is a proxy to your container. Here the dynamically allocated object is acquired by the digital surface, which will take care of its deletion.

Note
This creation of the container and the connection to a digital surface takes almost no time, since the chosen container (Light...) is lazy: the whole surface has not been extracted yet.

A digital surface is a graph, example of breadth-first traversal

Any DigitalSurface is a model of concepts::CUndirectedSimpleGraph (a refinement of concepts::CUndirectedSimpleLocalGraph). Essentially, you have methods DigitalSurface::begin() and DigitalSurface::end() to visit all the vertices, and for each vertex, you obtain its neighbors with DigitalSurface::writeNeighbors. You may thus for instance use the class BreadthFirstVisitor to do a breadth-first traversal of the digital surface.

We may extract the whole surface by doing a breadth-first traversal on the graph. The snippet below shows how to do it on any kind of digital surface, whatever the container.

trace.beginBlock( "Extracting boundary by tracking from an initial bel." );
typedef BreadthFirstVisitor<MyDigitalSurface> MyBreadthFirstVisitor;
typedef MyBreadthFirstVisitor::Node MyNode;
typedef MyBreadthFirstVisitor::Size MySize;
MyBreadthFirstVisitor visitor( digSurf, bel );
unsigned long nbSurfels = 0;
MyNode node;
while ( ! visitor.finished() )
{
node = visitor.current();
++nbSurfels;
visitor.expand();
}
MySize maxDist = node.second;

Here, we only use it to get the maximal distance to the starting bel. A second pass could be done to display cells with a color that depends on the distance to the starting surfel.

trace.beginBlock( "Displaying surface in Viewer3D." );
QApplication application(argc,argv);
Viewer3D<> viewer;
viewer.show();
HueShadeColorMap<MySize,1> hueShade( 0, maxDist );
MyBreadthFirstVisitor visitor2( digSurf, bel );
viewer << CustomColors3D( Color::Black, Color::White )
<< ks.unsigns( bel );
visitor2.expand();
while ( ! visitor2.finished() )
{
node = visitor2.current();
Color c = hueShade( node.second );
viewer << CustomColors3D( Color::Red, c )
<< ks.unsigns( node.first );
visitor2.expand();
}
viewer << Viewer3D<>::updateDisplay;
trace.info() << "nb surfels = " << nbSurfels << std::endl;
return application.exec();

We may call it as follows

# Commands
$ ./examples/topology/volBreadthFirstTraversal ../examples/samples/lobster.vol 50 255
$ ./examples/topology/volBreadthFirstTraversal ../examples/samples/Al.100.vol 0 1
$ ./examples/topology/volBreadthFirstTraversal ../examples/samples/cat10.vol 1 255

to get these pictures:

digital-surface-bfv-all.png
Examples of breadth-first traversal on a digital surface.
Note
In fact, the specific container LightImplicitDigitalSurface performs itself a breadth-first traversal to extract its vertices by itself. Therefere, you may use a simple begin() and end() to get the vertices in this order, when your container is a LightImplicitDigitalSurface.
Todo:
The concepts concepts::CUndirectedSimpleLocalGraph and concepts::CUndirectedSimpleGraph are susceptible to evolve to meet other standards.

Boundary and frontiers examples of digital surface

Surfels of a digital surface can also be defined by a predicate telling whether or not some surfel belongs to it. Given an image of labels, the class functors::BoundaryPredicate and functors::FrontierPredicate are such predicates:

  1. functors::BoundaryPredicate ( KSpace, Image, l ) is a predicate returning true for a surfel s iff the interior spel of s has label l while the exterior spel of s has label different from l.
  2. functors::FrontierPredicate ( KSpace, Image, l1, l2 ) is a predicate returning true for a surfel s iff the interior spel of s has label l1 while the exterior spel of s has label l2.

Using container ExplicitDigitalSurface, we can very simply define digital surfaces that are connected boundary of frontiers in some image. The following snippet is an excerpt from frontierAndBoundary.cpp.

trace.beginBlock( "Set up digital surface." );
typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
MySurfelAdjacency surfAdj( true ); // interior in all directions.
typedef functors::FrontierPredicate<KSpace, Image> FSurfelPredicate;
typedef ExplicitDigitalSurface<KSpace,FSurfelPredicate> FrontierContainer;
typedef DigitalSurface<FrontierContainer> Frontier;
typedef functors::BoundaryPredicate<KSpace, Image> BSurfelPredicate;
typedef ExplicitDigitalSurface<KSpace,BSurfelPredicate> BoundaryContainer;
typedef DigitalSurface<BoundaryContainer> Boundary;
// frontier between label 1 and 0 (connected part containing bel10)
SCell vox1 = K.sSpel( c1 + Point( radius1, 0, 0 ), K.POS );
SCell bel10 = K.sIncident( vox1, 0, true );
FSurfelPredicate surfPredicate10( K, image, 1, 0 );
Frontier frontier10 =
new FrontierContainer( K, surfPredicate10, surfAdj, bel10 );
// frontier between label 2 and 0 (connected part containing bel20)
SCell vox2 = K.sSpel( c2 - Point( radius2, 0, 0 ), K.POS );
SCell bel20 = K.sIncident( vox2, 0, false );
FSurfelPredicate surfPredicate20( K, image, 2, 0 );
Frontier frontier20 =
new FrontierContainer( K, surfPredicate20, surfAdj, bel20 );
// boundary of label 3 (connected part containing bel32)
SCell vox3 = K.sSpel( c1 - Point( radius1, 0, 0 ), K.POS );
SCell bel32 = K.sIncident( vox3, 0, false );
BSurfelPredicate surfPredicate3( K, image, 3 );
Boundary boundary3 =
new BoundaryContainer( K, surfPredicate3, surfAdj, bel32 );

Running it give the pictures:

digital-surface-intersecting-balls.png
Frontiers and boundary of two intersecting balls. The first ball is labelled 1 (red), the second 2 (yellow), their intersection 3 (orange). The frontier between 1 and 0 is displayed in red, the frontier between 2 and 0 is displayed in yellow, the boundary of region 3 is displayed in orange.

The digital surface graph is a combinatorial manifold

Umbrellas and faces

We have shown before that a digital surface has a dual structure that is the graph whose vertices are n-1-cells and whose arcs are (almost) n-2-cells. We can go further and use n-3-cells to define faces on this graph. This is related to the concept of umbrellas in 3D (see [57]). More precisely, given a start surfel, an incident n-2-cell (the separator) and an incident n-3-cell (the pivot), one can turn around the pivot progressively to get a face. This gives precisely in 3D the dual to a digital surface that is a kind of marching-cube surface (see [42], [43]).

The main class for computing umbrellas is the class UmbrellaComputer. Turning around the pivot means moving the face and the separator once (in the track direction), such that the pivot is the same (ie +p), the track and separator directions being updated. Repeating this process a sufficient number of times brings the umbrella back in its original position, except in the case when the DigitalSurface has a boundary touching the pivot. You may use it to compute umbrellas given a tracker on your surface.

You may look at file testUmbrellaComputer.cpp to see how to use this class directly.

Vertices, arcs and faces on a digital surface

However, it is simpler to use directly the methods defined in DigitalSurface. It offers three types: DigitalSurface::Vertex, DigitalSurface::Arc, DigitalSurface::Face to get the combinatorial structure of the surface. An arc is an oriented edge, a face is a sequence of edges that is closed when the pivot is not on the boundary and open otherwise.

The following code lists the faces of a DigitalSurface object, and for each face it lists the coordinates of its vertices.

// MyDigitalSurface is a DigitalSurface realization.
// digSurf is any object of type MyDigitalSurface.
// K is a cellular space.
MyDigitalSurface::FaceSet all_faces = digSurf.allFaces();
for ( MyDigitalSurface::FaceSet::const_iterator it = all_faces.begin(),
it_end = all_faces.end(); it != it_end; ++it )
{
std::cerr << " face=" << K.sKCoords( digSurf.pivot( *it ) ) << ":";
std::cerr << "(" << it->nbVertices << ")" << (it->isClosed() ? "C": "O");
MyDigitalSurface::VertexRange vtx = digSurf.verticesAroundFace( *it );
for ( unsigned int i = 0; i < vtx.size(); ++i )
{
std::cerr << " " << K.sKCoords( vtx[ i ] );
}
std::cerr << std::endl;
}

DigitalSurface provides you with many methods to get faces from edges or vertices and reciprocally.

Application to export surface in OFF format

In 3D, you may use faces to transform an arbitrary digital surface into a polygonal manifold (closed or open). You can for instance directly the method DigitalSurface::exportSurfaceAs3DOFF to do so, or look at its code to see how it works.

The following snippets of file volToOFF.cpp show how to extract all surfels in an image and then how to export the surface in OFF format. The output is a surface that is very much the classical marching-cube surface, except that vertices lies in the middle of the edge.

trace.beginBlock( "Extracting boundary by scanning the space. " );
typedef KSpace::SurfelSet SurfelSet;
typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
typedef DigitalSurface< MySetOfSurfels > MyDigitalSurface;
MySetOfSurfels theSetOfSurfels( K, surfAdj );
Surfaces<KSpace>::sMakeBoundary( theSetOfSurfels.surfelSet(),
K, set3d,
image.domain().lowerBound(),
image.domain().upperBound() );
MyDigitalSurface digSurf( theSetOfSurfels );
trace.info() << "Digital surface has " << digSurf.size() << " surfels."
<< std::endl;
trace.beginBlock( "Making OFF surface <marching-cube.off>. " );
ofstream out( "marching-cube.off" );
if ( out.good() )
digSurf.exportSurfaceAs3DOFF( out );
out.close();
digital-surface-mc-cat10.png
Marching-cube surface of cat10.vol file.
digital-surface-mc-lobster.png
Marching-cube surface of lobster.vol file.