DGtal  0.9.2
polyhedralizer.cpp
1 
30 #include <iostream>
33 #include <vector>
34 #include <set>
35 #include <map>
36 #include <queue>
37 
38 #include "DGtal/base/Common.h"
39 #include "DGtal/helpers/StdDefs.h"
40 #include "ConfigExamples.h"
42 
44 #include "DGtal/io/readers/VolReader.h"
45 #include "DGtal/images/ImageSelector.h"
46 #include "DGtal/images/ImageContainerBySTLVector.h"
47 #include "DGtal/images/SimpleThresholdForegroundPredicate.h"
49 
50 #include "DGtal/io/Display3D.h"
51 #include "DGtal/io/viewers/Viewer3D.h"
52 #include "DGtal/io/DrawWithDisplay3DModifier.h"
53 #include "DGtal/images/imagesSetsUtils/SetFromImage.h"
54 
55 #include "DGtal/topology/DigitalSurface.h"
56 #include "DGtal/topology/helpers/Surfaces.h"
57 #include "DGtal/topology/ImplicitDigitalSurface.h"
58 
59 #include "DGtal/graph/BreadthFirstVisitor.h"
60 #include "DGtal/geometry/surfaces/COBANaivePlaneComputer.h"
61 #include "DGtal/geometry/surfaces/ChordNaivePlaneComputer.h"
62 #include "DGtal/geometry/surfaces/ChordGenericNaivePlaneComputer.h"
63 
64 #include "DGtal/math/linalg/SimpleMatrix.h"
65 #include "DGtal/math/linalg/EigenDecomposition.h"
66 
68 
70 using namespace std;
71 using namespace DGtal;
72 using namespace Z3i;
74 
75 template <typename T1, typename T2>
76 struct PairSorted2nd
77 {
78  typedef PairSorted2nd<T1,T2> Self;
79  inline PairSorted2nd( const T1& t1, const T2& t2 ) : first( t1 ), second( t2 ) {}
80  bool operator<( const Self& other ) const
81  {
82  return second < other.second;
83  }
84  T1 first;
85  T2 second;
86 };
87 
88 template <typename T1, typename T2, typename T3>
89 struct Triple
90 {
91  T1 first;
92  T2 second;
93  T3 third;
94  Triple( T1 t1 = T1(), T2 t2 = T2(), T3 t3 = T3() )
95  : first( t1 ), second( t2 ), third( t3 )
96  {}
97 };
98 
99 // Least-Square Fit of a plane N.x=mu on points [itB,itE). Returns mu.
100 template <typename RealVector,
101  typename ConstIterator>
102 double LSF( RealVector& N, ConstIterator itB, ConstIterator itE )
103 {
104  typedef typename RealVector::Component Component;
105  typedef SimpleMatrix<Component,3,3> Matrix;
106  Matrix A; A.clear();
107  unsigned int nb = 0;
108  RealVector G = RealVector::zero; // center of gravity
109  for ( ConstIterator it = itB; it != itE; ++it )
110  {
111  G += RealVector( (*it)[ 0 ], (*it)[ 1 ], (*it)[ 2 ] );
112  ++nb;
113  }
114  G /= nb;
115  for ( ConstIterator it = itB; it != itE; ++it )
116  {
117  RealVector p( (*it)[ 0 ], (*it)[ 1 ], (*it)[ 2 ] );
118  p -= G;
119  for ( Dimension i = 0; i < 3; ++i )
120  for ( Dimension j = 0; j < 3; ++j )
121  A.setComponent( i, j, A( i, j ) + p[ i ] * p[ j ] );
122  }
123  // A is Gram matrix. We look for V such that V^t A V / |V|^2 is
124  // minimal. It is thus the first eigenvalue.
125  Matrix V;
126  RealVector values;
128  N = V.column( 0 ); // first eigenvector;
129  double mu = 0.0;
130  for ( ConstIterator it = itB; it != itE; ++it )
131  mu += N.dot( *it );
132  return mu/(double)nb;
133 }
134 
135 
136 int main( int argc, char** argv )
137 {
138  QApplication application(argc,argv);
139  string inputFilename = argc > 1 ? argv[ 1 ] : examplesPath+"/samples/Al.100.vol";
140  int threshold = argc > 2 ? atoi( argv[ 2 ] ) : 0;
141  int widthNum = argc > 3 ? atoi( argv[ 3 ] ) : 2;
142  int widthDen = argc > 4 ? atoi( argv[ 4 ] ) : 1;
143 
145  trace.beginBlock( "Reading vol file into an image." );
147  Image image = VolReader<Image>::importVol(inputFilename);
149  DigitalObject digitalObject( image, threshold );
150  trace.endBlock();
152 
154  trace.beginBlock( "Construct the Khalimsky space from the image domain." );
155  KSpace ks;
156  bool space_ok = ks.init( image.domain().lowerBound(), image.domain().upperBound(), true );
157  if (!space_ok)
158  {
159  trace.error() << "Error in the Khamisky space construction."<<endl;
160  return 2;
161  }
162  trace.endBlock();
164 
166  typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
167  MySurfelAdjacency surfAdj( false ); // exterior in all directions.
169 
171  trace.beginBlock( "Extracting boundary by tracking the surface. " );
172  typedef KSpace::Surfel Surfel;
173  Surfel start_surfel = Surfaces<KSpace>::findABel( ks, digitalObject, 100000 );
176  typedef MyDigitalSurface::ConstIterator ConstIterator;
177  MyContainer container( ks, digitalObject, surfAdj, start_surfel );
178  MyDigitalSurface digSurf( container );
179  trace.info() << "Digital surface has " << digSurf.size() << " surfels."
180  << endl;
181  trace.endBlock();
183 
185  // First pass to find biggest planes.
186  trace.beginBlock( "Decomposition first pass. Computes all planes so as to sort vertices by the plane size." );
189  map<Surfel,unsigned int> v2size;
190  for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
191  v2size[ *it ] = 0;
192  int j = 0;
193  int nb = digSurf.size();
194  NaivePlaneComputer planeComputer;
195  vector<Point> layer;
196  vector<Surfel> layer_surfel;
197  for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
198  {
199  if ( ( (++j) % 50 == 0 ) || ( j == nb ) ) trace.progressBar( j, nb );
200  Surfel v = *it;
201  planeComputer.init( widthNum, widthDen );
202  // The visitor takes care of all the breadth-first traversal.
203  Visitor visitor( digSurf, v );
204  layer.clear();
205  layer_surfel.clear();
206  Visitor::Size currentSize = visitor.current().second;
207  while ( ! visitor.finished() )
208  {
209  Visitor::Node node = visitor.current();
210  v = node.first;
211  int axis = ks.sOrthDir( v );
212  Point p = ks.sCoords( ks.sDirectIncident( v, axis ) );
213  if ( node.second != currentSize )
214  {
215  bool isExtended = planeComputer.extend( layer.begin(), layer.end() );
216  if ( isExtended )
217  {
218  for ( vector<Surfel>::const_iterator it_layer = layer_surfel.begin(),
219  it_layer_end = layer_surfel.end(); it_layer != it_layer_end; ++it_layer )
220  {
221  ++v2size[ *it_layer ];
222  }
223  layer_surfel.clear();
224  layer.clear();
225  currentSize = node.second;
226  }
227  else
228  break;
229  }
230  layer_surfel.push_back( v );
231  layer.push_back( p );
232  visitor.expand();
233  }
234  }
235  // Prepare queue
236  typedef PairSorted2nd<Surfel,int> SurfelWeight;
237  priority_queue<SurfelWeight> Q;
238  for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
239  Q.push( SurfelWeight( *it, v2size[ *it ] ) );
240  trace.endBlock();
242 
244  // Segmentation into planes
245  trace.beginBlock( "Decomposition second pass. Visits vertices from the one with biggest plane to the one with smallest plane." );
246  typedef Triple<NaivePlaneComputer, Color, pair<RealVector,double> > RoundPlane;
247  set<Surfel> processedVertices;
248  vector<RoundPlane*> roundPlanes;
249  map<Surfel,RoundPlane*> v2plane;
250  j = 0;
251  while ( ! Q.empty() )
252  {
253  if ( ( (++j) % 50 == 0 ) || ( j == nb ) ) trace.progressBar( j, nb );
254  Surfel v = Q.top().first;
255  Q.pop();
256  if ( processedVertices.find( v ) != processedVertices.end() ) // already in set
257  continue; // process to next vertex
258 
259  RoundPlane* ptrRoundPlane = new RoundPlane;
260  roundPlanes.push_back( ptrRoundPlane ); // to delete them afterwards.
261  v2plane[ v ] = ptrRoundPlane;
262  ptrRoundPlane->first.init( widthNum, widthDen );
263  ptrRoundPlane->third = make_pair( RealVector::zero, 0.0 );
264  // The visitor takes care of all the breadth-first traversal.
265  Visitor visitor( digSurf, v );
266  layer.clear();
267  layer_surfel.clear();
268  Visitor::Size currentSize = visitor.current().second;
269  while ( ! visitor.finished() )
270  {
271  Visitor::Node node = visitor.current();
272  v = node.first;
273  Dimension axis = ks.sOrthDir( v );
274  Point p = ks.sCoords( ks.sDirectIncident( v, axis ) );
275  if ( node.second != currentSize )
276  {
277  bool isExtended = ptrRoundPlane->first.extend( layer.begin(), layer.end() );
278  if ( isExtended )
279  {
280  for ( vector<Surfel>::const_iterator it_layer = layer_surfel.begin(),
281  it_layer_end = layer_surfel.end(); it_layer != it_layer_end; ++it_layer )
282  {
283  Surfel s = *it_layer;
284  processedVertices.insert( s );
285  if ( v2plane.find( s ) == v2plane.end() )
286  v2plane[ s ] = ptrRoundPlane;
287  }
288  layer.clear();
289  layer_surfel.clear();
290  currentSize = node.second;
291  }
292  else break;
293  }
294  layer_surfel.push_back( v );
295  layer.push_back( p );
296  if ( processedVertices.find( v ) != processedVertices.end() )
297  // surfel is already in some plane.
298  visitor.ignore();
299  else
300  visitor.expand();
301  }
302  if ( visitor.finished() )
303  {
304  for ( vector<Surfel>::const_iterator it_layer = layer_surfel.begin(),
305  it_layer_end = layer_surfel.end(); it_layer != it_layer_end; ++it_layer )
306  {
307  Surfel s = *it_layer;
308  processedVertices.insert( s );
309  if ( v2plane.find( s ) == v2plane.end() )
310  v2plane[ s ] = ptrRoundPlane;
311  }
312  }
313  // Assign random color for each plane.
314  ptrRoundPlane->second = Color( rand() % 192 + 64, rand() % 192 + 64, rand() % 192 + 64, 255 );
315  }
316  trace.endBlock();
318 
320  for ( vector<RoundPlane*>::iterator
321  it = roundPlanes.begin(), itE = roundPlanes.end();
322  it != itE; ++it )
323  {
324  NaivePlaneComputer& computer = (*it)->first;
325  RealVector normal;
326  double mu = LSF( normal, computer.begin(), computer.end() );
327  (*it)->third = make_pair( normal, mu );
328  }
330 
332  map<Surfel, RealPoint> coordinates;
333  for ( map<Surfel,RoundPlane*>::const_iterator
334  it = v2plane.begin(), itE = v2plane.end();
335  it != itE; ++it )
336  {
337  Surfel v = it->first;
338  RoundPlane* rplane = it->second;
339  Point p = ks.sKCoords( v );
340  RealPoint rp( (double)p[ 0 ]/2.0, (double)p[ 1 ]/2.0, (double)p[ 2 ]/2.0 );
341  double mu = rplane->third.second;
342  RealVector normal = rplane->third.first;
343  double lambda = mu - rp.dot( normal );
344  coordinates[ v ] = rp + lambda*normal;
345  }
346  typedef vector<Surfel> SurfelRange;
347  map<Surfel, RealPoint> new_coordinates;
348  for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
349  {
350  Surfel s = *it;
351  SurfelRange neighbors;
352  back_insert_iterator<SurfelRange> writeIt = back_inserter( neighbors );
353  digSurf.writeNeighbors( writeIt, *it );
354  RealPoint x = RealPoint::zero;
355  for ( SurfelRange::const_iterator its = neighbors.begin(), itsE = neighbors.end();
356  its != itsE; ++its )
357  x += coordinates[ *its ];
358  new_coordinates[ s ] = x / neighbors.size();
359  }
361 
363  typedef unsigned int Number;
364  typedef Mesh<RealPoint> MyMesh;
365  typedef MyMesh::MeshFace MeshFace;
366  typedef MyDigitalSurface::FaceSet FaceSet;
367  typedef MyDigitalSurface::VertexRange VertexRange;
368  map<Surfel, Number> index; // Numbers all vertices.
369  Number nbv = 0;
370  MyMesh polyhedron( true );
371  // Insert all projected surfels as vertices of the polyhedral surface.
372  for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
373  {
374  polyhedron.addVertex( new_coordinates[ *it ] );
375  index[ *it ] = nbv++;
376  }
377  // Define faces of the mesh. Outputs closed faces.
378  FaceSet faces = digSurf.allClosedFaces();
379  for ( typename FaceSet::const_iterator itf = faces.begin(), itf_end = faces.end();
380  itf != itf_end; ++itf )
381  {
382  MeshFace mface( itf->nbVertices );
383  VertexRange vtcs = digSurf.verticesAroundFace( *itf );
384  int i = 0;
385  for ( typename VertexRange::const_iterator itv = vtcs.begin(), itv_end = vtcs.end();
386  itv != itv_end; ++itv )
387  {
388  mface[ i++ ] = index[ *itv ];
389  }
390  polyhedron.addFace( mface, Color( 255, 243, 150, 255 ) );
391  }
393 
395  typedef Viewer3D<Space,KSpace> MyViewer3D;
396  MyViewer3D viewer( ks );
397  viewer.show();
398  bool isOK = polyhedron >> "test.off";
399  bool isOK2 = polyhedron >> "test.obj";
400  viewer << polyhedron;
401  viewer << MyViewer3D::updateDisplay;
402  application.exec();
404 
406  for ( vector<RoundPlane*>::iterator
407  it = roundPlanes.begin(), itE = roundPlanes.end();
408  it != itE; ++it )
409  delete *it;
411 
412  if (isOK && isOK2)
413  return 0;
414  else
415  return 1;
416 }
void beginBlock(const std::string &keyword="")
void init(Dimension axis, InternalInteger diameter, InternalInteger widthNumerator=NumberTraits< InternalInteger >::ONE, InternalInteger widthDenominator=NumberTraits< InternalInteger >::ONE)
void progressBar(const double currentValue, const double maximalValue)
std::vector< Vertex > VertexRange
The range of vertices is defined as a vector.
Aim: implements association bewteen points lying in a digital domain and values.
Definition: Image.h:69
Trace trace
Definition: Common.h:130
DigitalSurfaceContainer::SurfelConstIterator ConstIterator
Aim: A utility class for constructing surfaces (i.e. set of (n-1)-cells).
Definition: Surfaces.h:78
DGtal::uint32_t Dimension
Definition: Common.h:113
KhalimskySpaceND< 2, Integer > KSpace
Definition: StdDefs.h:77
Aim: This class is defined to represent a surface mesh through a set of vertices and faces...
Definition: Mesh.h:91
STL namespace.
double endBlock()
Aim: This class provides methods to compute the eigen decomposition of a matrix. Its objective is to ...
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
ConstIterator begin() const
std::set< Face > FaceSet
The set of faces is defined as set.
Aim: implements basic MxN Matrix services (M,N>=1).
Definition: SimpleMatrix.h:75
Aim: Define a simple Foreground predicate thresholding image values given a single thresold...
bool extend(const Point &p)
Aim: implements methods to read a "Vol" file format.
Definition: VolReader.h:88
ConstIterator end() const
DGtal is the top-level namespace which contains all DGtal functions and types.
Aim: A class that recognizes pieces of digital planes of given axis width. When the width is 1...
std::pair< Vertex, Data > Node
FIXME.
std::ostream & info()
Aim: This class is useful to perform a breadth-first exploration of a graph given a starting point or...
Structure representing an RGB triple with alpha component.
Definition: Color.h:66
const Domain & domain() const
Definition: Image.h:192
Space::RealPoint RealPoint
Definition: StdDefs.h:97
Space::RealVector RealVector
Definition: StdDefs.h:98
std::ostream & error()
Component dot(const Self &v) const
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as the boundary of an impl...