DGtal  0.9.4.1
Functions
polyhedralizer.cpp File Reference
#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"
Include dependency graph for polyhedralizer.cpp:

Go to the source code of this file.

Functions

template<typename RealVector , typename ConstIterator >
double LSF (RealVector &N, ConstIterator itB, ConstIterator itE)
 
int main (int argc, char **argv)
 

Detailed Description

This program is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.

Author
Jacques-Olivier Lachaud (jacqu.nosp@m.es-o.nosp@m.livie.nosp@m.r.la.nosp@m.chaud.nosp@m.@uni.nosp@m.v-sav.nosp@m.oie..nosp@m.fr ) Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
Date
2014/06/21

An example file named polyhedralizer.

This file is part of the DGtal library.

Definition in file polyhedralizer.cpp.

Function Documentation

◆ LSF()

template<typename RealVector , typename ConstIterator >
double LSF ( RealVector N,
ConstIterator  itB,
ConstIterator  itE 
)
Examples:
examples/tutorial-examples/polyhedralizer.cpp, and tutorial-examples/polyhedralizer.cpp.

Definition at line 125 of file polyhedralizer.cpp.

References DGtal::SimpleMatrix< TComponent, TM, TN >::clear(), DGtal::PointVector< dim, TEuclideanRing, TContainer >::dot(), and DGtal::PointVector< dim, TEuclideanRing, TContainer >::zero.

Referenced by main().

126 {
127  typedef typename RealVector::Component Component;
128  typedef SimpleMatrix<Component,3,3> Matrix;
129  Matrix A; A.clear();
130  unsigned int nb = 0;
131  RealVector G = RealVector::zero; // center of gravity
132  for ( ConstIterator it = itB; it != itE; ++it )
133  {
134  G += RealVector( (*it)[ 0 ], (*it)[ 1 ], (*it)[ 2 ] );
135  ++nb;
136  }
137  G /= nb;
138  for ( ConstIterator it = itB; it != itE; ++it )
139  {
140  RealVector p( (*it)[ 0 ], (*it)[ 1 ], (*it)[ 2 ] );
141  p -= G;
142  for ( Dimension i = 0; i < 3; ++i )
143  for ( Dimension j = 0; j < 3; ++j )
144  A.setComponent( i, j, A( i, j ) + p[ i ] * p[ j ] );
145  }
146  // A is Gram matrix. We look for V such that V^t A V / |V|^2 is
147  // minimal. It is thus the first eigenvalue.
148  Matrix V;
149  RealVector values;
151  N = V.column( 0 ); // first eigenvector;
152  double mu = 0.0;
153  for ( ConstIterator it = itB; it != itE; ++it )
154  mu += N.dot( *it );
155  return mu/(double)nb;
156 }
MyDigitalSurface::ConstIterator ConstIterator
Component dot(const Self &v) const
DGtal::uint32_t Dimension
Definition: Common.h:120
Aim: This class provides methods to compute the eigen decomposition of a matrix. Its objective is to ...
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:141
Aim: implements basic MxN Matrix services (M,N>=1).
Definition: SimpleMatrix.h:75
static Self zero
Static const for zero PointVector.
Definition: PointVector.h:798
TEuclideanRing Component
Type for Vector elements.
Definition: PointVector.h:155

◆ main()

int main ( int  argc,
char **  argv 
)

[polyhedralizer-readVol]

[polyhedralizer-readVol]

[polyhedralizer-KSpace]

[polyhedralizer-KSpace]

[polyhedralizer-SurfelAdjacency]

[polyhedralizer-SurfelAdjacency]

[polyhedralizer-ExtractingSurface]

[polyhedralizer-ExtractingSurface]

[polyhedralizer-ComputingPlaneSize]

[polyhedralizer-ComputingPlaneSize]

[polyhedralizer-segment]

[polyhedralizer-segment]

[polyhedralizer-lsf]

[polyhedralizer-lsf]

[polyhedralizer-projection]

[polyhedralizer-projection]

[polyhedralizer-MakeMesh]

[polyhedralizer-MakeMesh]

[polyhedralizer-visualization]

[polyhedralizer-visualization]

[polyhedralizer-freeMemory]

[polyhedralizer-freeMemory]

Definition at line 159 of file polyhedralizer.cpp.

References DGtal::DigitalSurface< TDigitalSurfaceContainer >::allClosedFaces(), DGtal::COBANaivePlaneComputer< TSpace, TInternalInteger >::begin(), DGtal::DigitalSurface< TDigitalSurfaceContainer >::begin(), DGtal::Trace::beginBlock(), DGtal::BreadthFirstVisitor< TGraph, TMarkSet >::current(), DGtal::PointVector< dim, TEuclideanRing, TContainer >::dot(), DGtal::COBANaivePlaneComputer< TSpace, TInternalInteger >::end(), DGtal::DigitalSurface< TDigitalSurfaceContainer >::end(), DGtal::Trace::endBlock(), DGtal::Trace::error(), DGtal::BreadthFirstVisitor< TGraph, TMarkSet >::expand(), DGtal::COBANaivePlaneComputer< TSpace, TInternalInteger >::extend(), DGtal::BreadthFirstVisitor< TGraph, TMarkSet >::finished(), DGtal::BreadthFirstVisitor< TGraph, TMarkSet >::ignore(), index(), DGtal::Trace::info(), DGtal::COBANaivePlaneComputer< TSpace, TInternalInteger >::init(), DGtal::KhalimskySpaceND< dim, TInteger >::init(), LSF(), DGtal::Trace::progressBar(), DGtal::KhalimskySpaceND< dim, TInteger >::sCoords(), DGtal::KhalimskySpaceND< dim, TInteger >::sDirectIncident(), DGtal::DigitalSurface< TDigitalSurfaceContainer >::size(), DGtal::KhalimskySpaceND< dim, TInteger >::sKCoords(), DGtal::KhalimskySpaceND< dim, TInteger >::sOrthDir(), DGtal::trace, DGtal::DigitalSurface< TDigitalSurfaceContainer >::verticesAroundFace(), DGtal::DigitalSurface< TDigitalSurfaceContainer >::writeNeighbors(), DGtal::PointVector< 3, double >::zero, and DGtal::PointVector< dim, TEuclideanRing, TContainer >::zero.

160 {
161  QApplication application(argc,argv);
162  string inputFilename = argc > 1 ? argv[ 1 ] : examplesPath+"/samples/Al.100.vol";
163  int threshold = argc > 2 ? atoi( argv[ 2 ] ) : 0;
164  int widthNum = argc > 3 ? atoi( argv[ 3 ] ) : 2;
165  int widthDen = argc > 4 ? atoi( argv[ 4 ] ) : 1;
166 
168  trace.beginBlock( "Reading vol file into an image." );
170  Image image = VolReader<Image>::importVol(inputFilename);
172  DigitalObject digitalObject( image, threshold );
173  trace.endBlock();
175 
177  trace.beginBlock( "Construct the Khalimsky space from the image domain." );
178  KSpace ks;
179  bool space_ok = ks.init( image.domain().lowerBound(), image.domain().upperBound(), true );
180  if (!space_ok)
181  {
182  trace.error() << "Error in the Khamisky space construction."<<endl;
183  return 2;
184  }
185  trace.endBlock();
187 
189  typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
190  MySurfelAdjacency surfAdj( false ); // exterior in all directions.
192 
194  trace.beginBlock( "Extracting boundary by tracking the surface. " );
195  typedef KSpace::Surfel Surfel;
196  Surfel start_surfel = Surfaces<KSpace>::findABel( ks, digitalObject, 100000 );
200  MyContainer container( ks, digitalObject, surfAdj, start_surfel );
201  MyDigitalSurface digSurf( container );
202  trace.info() << "Digital surface has " << digSurf.size() << " surfels."
203  << endl;
204  trace.endBlock();
206 
208  // First pass to find biggest planes.
209  trace.beginBlock( "Decomposition first pass. Computes all planes so as to sort vertices by the plane size." );
212  map<Surfel,unsigned int> v2size;
213  for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
214  v2size[ *it ] = 0;
215  int j = 0;
216  int nb = digSurf.size();
217  NaivePlaneComputer planeComputer;
218  vector<Point> layer;
219  vector<Surfel> layer_surfel;
220  for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
221  {
222  if ( ( (++j) % 50 == 0 ) || ( j == nb ) ) trace.progressBar( j, nb );
223  Surfel v = *it;
224  planeComputer.init( widthNum, widthDen );
225  // The visitor takes care of all the breadth-first traversal.
226  Visitor visitor( digSurf, v );
227  layer.clear();
228  layer_surfel.clear();
229  Visitor::Size currentSize = visitor.current().second;
230  while ( ! visitor.finished() )
231  {
232  Visitor::Node node = visitor.current();
233  v = node.first;
234  int axis = ks.sOrthDir( v );
235  Point p = ks.sCoords( ks.sDirectIncident( v, axis ) );
236  if ( node.second != currentSize )
237  {
238  bool isExtended = planeComputer.extend( layer.begin(), layer.end() );
239  if ( isExtended )
240  {
241  for ( vector<Surfel>::const_iterator it_layer = layer_surfel.begin(),
242  it_layer_end = layer_surfel.end(); it_layer != it_layer_end; ++it_layer )
243  {
244  ++v2size[ *it_layer ];
245  }
246  layer_surfel.clear();
247  layer.clear();
248  currentSize = node.second;
249  }
250  else
251  break;
252  }
253  layer_surfel.push_back( v );
254  layer.push_back( p );
255  visitor.expand();
256  }
257  }
258  // Prepare queue
259  typedef PairSorted2nd<Surfel,int> SurfelWeight;
260  priority_queue<SurfelWeight> Q;
261  for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
262  Q.push( SurfelWeight( *it, v2size[ *it ] ) );
263  trace.endBlock();
265 
267  // Segmentation into planes
268  trace.beginBlock( "Decomposition second pass. Visits vertices from the one with biggest plane to the one with smallest plane." );
269  typedef Triple<NaivePlaneComputer, Color, pair<RealVector,double> > RoundPlane;
270  set<Surfel> processedVertices;
271  vector<RoundPlane*> roundPlanes;
272  map<Surfel,RoundPlane*> v2plane;
273  j = 0;
274  while ( ! Q.empty() )
275  {
276  if ( ( (++j) % 50 == 0 ) || ( j == nb ) ) trace.progressBar( j, nb );
277  Surfel v = Q.top().first;
278  Q.pop();
279  if ( processedVertices.find( v ) != processedVertices.end() ) // already in set
280  continue; // process to next vertex
281 
282  RoundPlane* ptrRoundPlane = new RoundPlane;
283  roundPlanes.push_back( ptrRoundPlane ); // to delete them afterwards.
284  v2plane[ v ] = ptrRoundPlane;
285  ptrRoundPlane->first.init( widthNum, widthDen );
286  ptrRoundPlane->third = make_pair( RealVector::zero, 0.0 );
287  // The visitor takes care of all the breadth-first traversal.
288  Visitor visitor( digSurf, v );
289  layer.clear();
290  layer_surfel.clear();
291  Visitor::Size currentSize = visitor.current().second;
292  while ( ! visitor.finished() )
293  {
294  Visitor::Node node = visitor.current();
295  v = node.first;
296  Dimension axis = ks.sOrthDir( v );
297  Point p = ks.sCoords( ks.sDirectIncident( v, axis ) );
298  if ( node.second != currentSize )
299  {
300  bool isExtended = ptrRoundPlane->first.extend( layer.begin(), layer.end() );
301  if ( isExtended )
302  {
303  for ( vector<Surfel>::const_iterator it_layer = layer_surfel.begin(),
304  it_layer_end = layer_surfel.end(); it_layer != it_layer_end; ++it_layer )
305  {
306  Surfel s = *it_layer;
307  processedVertices.insert( s );
308  if ( v2plane.find( s ) == v2plane.end() )
309  v2plane[ s ] = ptrRoundPlane;
310  }
311  layer.clear();
312  layer_surfel.clear();
313  currentSize = node.second;
314  }
315  else break;
316  }
317  layer_surfel.push_back( v );
318  layer.push_back( p );
319  if ( processedVertices.find( v ) != processedVertices.end() )
320  // surfel is already in some plane.
321  visitor.ignore();
322  else
323  visitor.expand();
324  }
325  if ( visitor.finished() )
326  {
327  for ( vector<Surfel>::const_iterator it_layer = layer_surfel.begin(),
328  it_layer_end = layer_surfel.end(); it_layer != it_layer_end; ++it_layer )
329  {
330  Surfel s = *it_layer;
331  processedVertices.insert( s );
332  if ( v2plane.find( s ) == v2plane.end() )
333  v2plane[ s ] = ptrRoundPlane;
334  }
335  }
336  // Assign random color for each plane.
337  ptrRoundPlane->second = Color( rand() % 192 + 64, rand() % 192 + 64, rand() % 192 + 64, 255 );
338  }
339  trace.endBlock();
341 
343  for ( vector<RoundPlane*>::iterator
344  it = roundPlanes.begin(), itE = roundPlanes.end();
345  it != itE; ++it )
346  {
347  NaivePlaneComputer& computer = (*it)->first;
348  RealVector normal;
349  double mu = LSF( normal, computer.begin(), computer.end() );
350  (*it)->third = make_pair( normal, mu );
351  }
353 
355  map<Surfel, RealPoint> coordinates;
356  for ( map<Surfel,RoundPlane*>::const_iterator
357  it = v2plane.begin(), itE = v2plane.end();
358  it != itE; ++it )
359  {
360  Surfel v = it->first;
361  RoundPlane* rplane = it->second;
362  Point p = ks.sKCoords( v );
363  RealPoint rp( (double)p[ 0 ]/2.0, (double)p[ 1 ]/2.0, (double)p[ 2 ]/2.0 );
364  double mu = rplane->third.second;
365  RealVector normal = rplane->third.first;
366  double lambda = mu - rp.dot( normal );
367  coordinates[ v ] = rp + lambda*normal;
368  }
369  typedef vector<Surfel> SurfelRange;
370  map<Surfel, RealPoint> new_coordinates;
371  for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
372  {
373  Surfel s = *it;
374  SurfelRange neighbors;
375  back_insert_iterator<SurfelRange> writeIt = back_inserter( neighbors );
376  digSurf.writeNeighbors( writeIt, *it );
378  for ( SurfelRange::const_iterator its = neighbors.begin(), itsE = neighbors.end();
379  its != itsE; ++its )
380  x += coordinates[ *its ];
381  new_coordinates[ s ] = x / neighbors.size();
382  }
384 
386  typedef unsigned int Number;
387  typedef Mesh<RealPoint> MyMesh;
388  typedef MyMesh::MeshFace MeshFace;
389  typedef MyDigitalSurface::FaceSet FaceSet;
391  map<Surfel, Number> index; // Numbers all vertices.
392  Number nbv = 0;
393  MyMesh polyhedron( true );
394  // Insert all projected surfels as vertices of the polyhedral surface.
395  for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
396  {
397  polyhedron.addVertex( new_coordinates[ *it ] );
398  index[ *it ] = nbv++;
399  }
400  // Define faces of the mesh. Outputs closed faces.
401  FaceSet faces = digSurf.allClosedFaces();
402  for ( typename FaceSet::const_iterator itf = faces.begin(), itf_end = faces.end();
403  itf != itf_end; ++itf )
404  {
405  MeshFace mface( itf->nbVertices );
406  VertexRange vtcs = digSurf.verticesAroundFace( *itf );
407  int i = 0;
408  for ( typename VertexRange::const_iterator itv = vtcs.begin(), itv_end = vtcs.end();
409  itv != itv_end; ++itv )
410  {
411  mface[ i++ ] = index[ *itv ];
412  }
413  polyhedron.addFace( mface, Color( 255, 243, 150, 255 ) );
414  }
416 
418  typedef Viewer3D<Space,KSpace> MyViewer3D;
419  MyViewer3D viewer( ks );
420  viewer.show();
421  bool isOK = polyhedron >> "test.off";
422  bool isOK2 = polyhedron >> "test.obj";
423  viewer << polyhedron;
424  viewer << MyViewer3D::updateDisplay;
425  application.exec();
427 
429  for ( vector<RoundPlane*>::iterator
430  it = roundPlanes.begin(), itE = roundPlanes.end();
431  it != itE; ++it )
432  delete *it;
434 
435  if (isOK && isOK2)
436  return 0;
437  else
438  return 1;
439 }
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
MyDigitalSurface::ConstIterator ConstIterator
unsigned int index(DGtal::uint32_t n, unsigned int b)
Definition: testBits.cpp:44
Trace trace
Definition: Common.h:137
TriMesh::VertexRange VertexRange
DigitalSurfaceContainer::SurfelConstIterator ConstIterator
COBANaivePlaneComputer< Z3, InternalInteger > NaivePlaneComputer
ConstIterator begin() const
Component dot(const Self &v) const
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:120
Aim: This class is defined to represent a surface mesh through a set of vertices and faces...
Definition: Mesh.h:91
double endBlock()
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
BreadthFirstVisitor< MyDigitalSurface > Visitor
ConstIterator end() const
std::set< Face > FaceSet
The set of faces is defined as set.
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:141
Dimension sOrthDir(const SCell &s) const
Aim: Define a simple Foreground predicate thresholding image values given a single thresold...
bool init(const Point &lower, const Point &upper, bool isClosed)
bool extend(const Point &p)
Aim: implements methods to read a "Vol" file format.
Definition: VolReader.h:89
const Point & sKCoords(const SCell &c) const
Point sCoords(const SCell &c) const
MyPointD Point
Definition: testClone2.cpp:383
SCell sDirectIncident(const SCell &p, Dimension k) const
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()
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value...
ImageContainerBySTLVector< Domain, Value > Image
Aim: This class is useful to perform a breadth-first exploration of a graph given a starting point or...
double LSF(RealVector &N, ConstIterator itB, ConstIterator itE)
Structure representing an RGB triple with alpha component.
Definition: Color.h:66
static Self zero
Static const for zero PointVector.
Definition: PointVector.h:798
std::ostream & error()
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex...
DigitalSurface< MyDigitalSurfaceContainer > MyDigitalSurface
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as the boundary of an impl...