DGtal  1.3.beta
tutorial-examples/polyhedralizer.cpp

Example of tutorial 2: making a polyhedron from a digital object

See also
Tutorial: making a polyhedron from a digital object
The polyhedral surface approaching Al Capone digital object, for width=3/1.
#include <iostream>
#include <vector>
#include <set>
#include <map>
#include <queue>
#include "DGtal/base/Common.h"
#include "DGtal/helpers/StdDefs.h"
#include "ConfigExamples.h"
#include "DGtal/io/readers/VolReader.h"
#include "DGtal/images/ImageSelector.h"
#include "DGtal/images/ImageContainerBySTLVector.h"
#include "DGtal/images/SimpleThresholdForegroundPredicate.h"
#include "DGtal/io/Display3D.h"
#include "DGtal/io/viewers/Viewer3D.h"
#include "DGtal/io/DrawWithDisplay3DModifier.h"
#include "DGtal/images/imagesSetsUtils/SetFromImage.h"
#include "DGtal/topology/DigitalSurface.h"
#include "DGtal/topology/helpers/Surfaces.h"
#include "DGtal/topology/ImplicitDigitalSurface.h"
#include "DGtal/graph/BreadthFirstVisitor.h"
#include "DGtal/geometry/surfaces/COBANaivePlaneComputer.h"
#include "DGtal/geometry/surfaces/ChordNaivePlaneComputer.h"
#include "DGtal/geometry/surfaces/ChordGenericNaivePlaneComputer.h"
#include "DGtal/math/linalg/SimpleMatrix.h"
#include "DGtal/math/linalg/EigenDecomposition.h"
using namespace std;
using namespace DGtal;
using namespace Z3i;
template <typename T1, typename T2>
struct PairSorted2nd
{
typedef PairSorted2nd<T1,T2> Self;
inline 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;
};
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 )
{}
};
// 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;
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;
}
int main( int argc, char** argv )
{
QApplication application(argc,argv);
string inputFilename = argc > 1 ? argv[ 1 ] : examplesPath+"/samples/Al.100.vol";
int threshold = argc > 2 ? atoi( argv[ 2 ] ) : 0;
int widthNum = argc > 3 ? atoi( argv[ 3 ] ) : 2;
int widthDen = argc > 4 ? atoi( argv[ 4 ] ) : 1;
trace.beginBlock( "Reading vol file into an image." );
Image image = VolReader<Image>::importVol(inputFilename);
DigitalObject digitalObject( image, threshold );
trace.beginBlock( "Construct the Khalimsky space from the image domain." );
KSpace ks;
bool space_ok = ks.init( image.domain().lowerBound(), image.domain().upperBound(), true );
if (!space_ok)
{
trace.error() << "Error in the Khamisky space construction."<<endl;
return 2;
}
typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
MySurfelAdjacency surfAdj( false ); // exterior in all directions.
trace.beginBlock( "Extracting boundary by tracking the surface. " );
Surfel start_surfel = Surfaces<KSpace>::findABel( ks, digitalObject, 100000 );
MyContainer container( ks, digitalObject, surfAdj, start_surfel );
MyDigitalSurface digSurf( container );
trace.info() << "Digital surface has " << digSurf.size() << " surfels."
<< endl;
// First pass to find biggest planes.
trace.beginBlock( "Decomposition first pass. Computes all planes so as to sort vertices by the plane size." );
map<Surfel,unsigned int> v2size;
for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
v2size[ *it ] = 0;
int j = 0;
int nb = digSurf.size();
NaivePlaneComputer planeComputer;
vector<Point> layer;
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;
planeComputer.init( widthNum, widthDen );
// The visitor takes care of all the breadth-first traversal.
Visitor visitor( digSurf, v );
layer.clear();
layer_surfel.clear();
Visitor::Size currentSize = visitor.current().second;
while ( ! visitor.finished() )
{
Visitor::Node node = visitor.current();
v = node.first;
int axis = ks.sOrthDir( v );
Point p = ks.sCoords( ks.sDirectIncident( v, axis ) );
if ( node.second != currentSize )
{
bool isExtended = planeComputer.extend( layer.begin(), layer.end() );
if ( isExtended )
{
for ( vector<Surfel>::const_iterator it_layer = layer_surfel.begin(),
it_layer_end = layer_surfel.end(); it_layer != it_layer_end; ++it_layer )
{
++v2size[ *it_layer ];
}
layer_surfel.clear();
layer.clear();
currentSize = node.second;
}
else
break;
}
layer_surfel.push_back( v );
layer.push_back( p );
visitor.expand();
}
}
// Prepare queue
typedef PairSorted2nd<Surfel,int> SurfelWeight;
priority_queue<SurfelWeight> Q;
for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
Q.push( SurfelWeight( *it, v2size[ *it ] ) );
// Segmentation into planes
trace.beginBlock( "Decomposition second pass. Visits vertices from the one with biggest plane to the one with smallest plane." );
typedef Triple<NaivePlaneComputer, Color, pair<RealVector,double> > RoundPlane;
set<Surfel> processedVertices;
vector<RoundPlane*> roundPlanes;
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 ( processedVertices.find( v ) != processedVertices.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->first.init( widthNum, widthDen );
ptrRoundPlane->third = make_pair( RealVector::zero, 0.0 );
// The visitor takes care of all the breadth-first traversal.
Visitor visitor( digSurf, v );
layer.clear();
layer_surfel.clear();
Visitor::Size currentSize = visitor.current().second;
while ( ! visitor.finished() )
{
Visitor::Node node = visitor.current();
v = node.first;
Dimension axis = ks.sOrthDir( v );
Point p = ks.sCoords( ks.sDirectIncident( v, axis ) );
if ( node.second != currentSize )
{
bool isExtended = ptrRoundPlane->first.extend( layer.begin(), layer.end() );
if ( isExtended )
{
for ( vector<Surfel>::const_iterator it_layer = layer_surfel.begin(),
it_layer_end = layer_surfel.end(); it_layer != it_layer_end; ++it_layer )
{
Surfel s = *it_layer;
processedVertices.insert( s );
if ( v2plane.find( s ) == v2plane.end() )
v2plane[ s ] = ptrRoundPlane;
}
layer.clear();
layer_surfel.clear();
currentSize = node.second;
}
else break;
}
layer_surfel.push_back( v );
layer.push_back( p );
if ( processedVertices.find( v ) != processedVertices.end() )
// surfel is already in some plane.
visitor.ignore();
else
visitor.expand();
}
if ( visitor.finished() )
{
for ( vector<Surfel>::const_iterator it_layer = layer_surfel.begin(),
it_layer_end = layer_surfel.end(); it_layer != it_layer_end; ++it_layer )
{
Surfel s = *it_layer;
processedVertices.insert( s );
if ( v2plane.find( s ) == v2plane.end() )
v2plane[ s ] = ptrRoundPlane;
}
}
// Assign random color for each plane.
ptrRoundPlane->second = Color( rand() % 192 + 64, rand() % 192 + 64, rand() % 192 + 64, 255 );
}
for ( 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 );
}
map<Surfel, RealPoint> coordinates;
for ( 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 vector<Surfel> SurfelRange;
map<Surfel, RealPoint> new_coordinates;
for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
{
Surfel s = *it;
SurfelRange neighbors;
back_insert_iterator<SurfelRange> writeIt = back_inserter( neighbors );
digSurf.writeNeighbors( writeIt, *it );
for ( SurfelRange::const_iterator its = neighbors.begin(), itsE = neighbors.end();
its != itsE; ++its )
x += coordinates[ *its ];
new_coordinates[ s ] = x / neighbors.size();
}
typedef unsigned int Number;
typedef Mesh<RealPoint> MyMesh;
typedef MyMesh::MeshFace MeshFace;
typedef MyDigitalSurface::FaceSet FaceSet;
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++;
}
// Define faces of the mesh. Outputs closed faces.
FaceSet faces = digSurf.allClosedFaces();
for ( typename FaceSet::const_iterator itf = faces.begin(), itf_end = faces.end();
itf != itf_end; ++itf )
{
MeshFace mface( itf->nbVertices );
VertexRange vtcs = digSurf.verticesAroundFace( *itf );
int i = 0;
for ( typename VertexRange::const_iterator itv = vtcs.begin(), itv_end = vtcs.end();
itv != itv_end; ++itv )
{
mface[ i++ ] = index[ *itv ];
}
polyhedron.addFace( mface, Color( 255, 243, 150, 255 ) );
}
typedef Viewer3D<Space,KSpace> MyViewer3D;
MyViewer3D viewer( ks );
viewer.show();
bool isOK = polyhedron >> "test.off";
bool isOK2 = polyhedron >> "test.obj";
viewer << polyhedron;
viewer << MyViewer3D::updateDisplay;
application.exec();
for ( vector<RoundPlane*>::iterator
it = roundPlanes.begin(), itE = roundPlanes.end();
it != itE; ++it )
delete *it;
if (isOK && isOK2)
return 0;
else
return 1;
}
operator<
bool operator<(const VertexSize &vs1, const VertexSize &vs2)
Definition: greedy-plane-segmentation-ex2.cpp:115
DGtal::Surfaces
Aim: A utility class for constructing surfaces (i.e. set of (n-1)-cells).
Definition: Surfaces.h:78
ConstIterator
MyDigitalSurface::ConstIterator ConstIterator
Definition: greedy-plane-segmentation-ex2.cpp:93
DGtal::PointVector::zero
static Self zero
Static const for zero PointVector.
Definition: PointVector.h:1595
DGtal::DigitalSurface::ConstIterator
DigitalSurfaceContainer::SurfelConstIterator ConstIterator
Definition: DigitalSurface.h:163
NaivePlaneComputer
COBANaivePlaneComputer< Z3, InternalInteger > NaivePlaneComputer
Definition: greedy-plane-segmentation-ex2.cpp:88
DGtal::COBANaivePlaneComputer::init
void init(Dimension axis, InternalInteger diameter, InternalInteger widthNumerator=NumberTraits< InternalInteger >::ONE, InternalInteger widthDenominator=NumberTraits< InternalInteger >::ONE)
DGtal::Trace::endBlock
double endBlock()
DGtal::DigitalSurface::VertexRange
std::vector< Vertex > VertexRange
The range of vertices is defined as a vector.
Definition: DigitalSurface.h:305
DGtal::COBANaivePlaneComputer::extend
bool extend(const Point &p)
DGtal::ImageContainerBySTLVector
Definition: ImageContainerBySTLVector.h:126
DGtal::DigitalSurface
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
Definition: DigitalSurface.h:139
DGtal::SurfelAdjacency< KSpace::dimension >
DGtal::Trace::error
std::ostream & error()
DGtal::Color
Structure representing an RGB triple with alpha component.
Definition: Color.h:66
DGtal::KhalimskySpaceND::init
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
RealVector
Z3i::RealVector RealVector
Definition: sphereCotangentLaplaceOperator.cpp:71
DGtal::EigenDecomposition
Aim: This class provides methods to compute the eigen decomposition of a matrix. Its objective is to ...
Definition: EigenDecomposition.h:86
DGtal::trace
Trace trace
Definition: Common.h:154
DGtal::COBANaivePlaneComputer::end
ConstIterator end() const
DGtal::ImplicitDigitalSurface
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as the boundary of an impl...
Definition: ImplicitDigitalSurface.h:71
DGtal::BreadthFirstVisitor
Aim: This class is useful to perform a breadth-first exploration of a graph given a starting point or...
Definition: BreadthFirstVisitor.h:94
DGtal::Dimension
DGtal::uint32_t Dimension
Definition: Common.h:137
DGtal::BreadthFirstVisitor::Node
std::pair< Vertex, Data > Node
FIXME.
Definition: BreadthFirstVisitor.h:115
LSF
double LSF(RealVector &N, ConstIterator itB, ConstIterator itE)
Definition: polyhedralizer.cpp:125
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
Visitor
BreadthFirstVisitor< MyDigitalSurface > Visitor
Definition: greedy-plane-segmentation-ex2.cpp:97
DGtal::KhalimskySpaceND::sCoords
Point sCoords(const SCell &c) const
Return its digital coordinates.
DGtal::SignedKhalimskyCell
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
Definition: KhalimskySpaceND.h:208
VertexRange
TriMesh::VertexRange VertexRange
Definition: testTriangulatedSurface.cpp:53
DGtal::DigitalSurface::FaceSet
std::set< Face > FaceSet
The set of faces is defined as set.
Definition: DigitalSurface.h:307
DGtal::PointVector::Component
TEuclideanRing Component
Type for Vector elements.
Definition: PointVector.h:614
DGtal::Trace::info
std::ostream & info()
DGtal::Trace::progressBar
void progressBar(const double currentValue, const double maximalValue)
Surfel
KSpace::SCell Surfel
Definition: testArithmeticalDSSComputerOnSurfels.cpp:50
DGtal::KhalimskySpaceND::sDirectIncident
SCell sDirectIncident(const SCell &p, Dimension k) const
Return the direct incident cell of [p] along [k] (the incident cell along [k])
DGtal::Viewer3D< Space, KSpace >
DGtal::SimpleMatrix
Aim: implements basic MxN Matrix services (M,N>=1).
Definition: SimpleMatrix.h:75
DGtal::Mesh
Aim: This class is defined to represent a surface mesh through a set of vertices and faces....
Definition: Mesh.h:91
Image
ImageContainerBySTLVector< Domain, Value > Image
Definition: testSimpleRandomAccessRangeFromPoint.cpp:45
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
index
unsigned int index(DGtal::uint32_t n, unsigned int b)
Definition: testBits.cpp:44
DGtal::KhalimskySpaceND::sKCoords
const Point & sKCoords(const SCell &c) const
Return its Khalimsky coordinates.
DGtal::ChordGenericNaivePlaneComputer
Aim: A class that recognizes pieces of digital planes of given axis width. When the width is 1,...
Definition: ChordGenericNaivePlaneComputer.h:118
DGtal::SimpleMatrix::clear
void clear()
main
int main(int argc, char **argv)
Definition: testArithmeticDSS-benchmark.cpp:147
DGtal::COBANaivePlaneComputer::begin
ConstIterator begin() const
DGtal::KhalimskySpaceND::sOrthDir
Dimension sOrthDir(const SCell &s) const
Given a signed surfel [s], returns its orthogonal direction (ie, the coordinate where the surfel is c...
DGtal::Image
Aim: implements association bewteen points lying in a digital domain and values.
Definition: Image.h:69
DGtal::PointVector::dot
auto dot(const PointVector< dim, OtherComponent, OtherStorage > &v) const -> decltype(DGtal::dotProduct(*this, v))
Dot product with a PointVector.
DGtal::PointVector
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:165
DGtal::VolReader
Aim: implements methods to read a "Vol" file format.
Definition: VolReader.h:89
DGtal::BreadthFirstVisitor::Size
Graph::Size Size
Definition: BreadthFirstVisitor.h:101
DGtal::COBANaivePlaneComputer< Z3, InternalInteger >
Point
MyPointD Point
Definition: testClone2.cpp:383
MyDigitalSurface
DigitalSurface< MyDigitalSurfaceContainer > MyDigitalSurface
Definition: greedy-plane-segmentation-ex2.cpp:92
DGtal::KhalimskySpaceND
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
Definition: KhalimskySpaceND.h:64
DGtal::functors::SimpleThresholdForegroundPredicate
Aim: Define a simple Foreground predicate thresholding image values given a single thresold....
Definition: SimpleThresholdForegroundPredicate.h:65