File failed to load: https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.3/config/TeX-MML-AM_CHTML/MathJax.js
DGtal 2.0.0
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/viewers/PolyscopeViewer.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.

Data Structures

struct  PairSorted2nd< T1, T2 >
 [polyhedralizer-typedefs] More...
struct  Triple< T1, T2, T3 >

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 123 of file polyhedralizer.cpp.

124{
125 typedef typename RealVector::Component Component;
126 typedef SimpleMatrix<Component,3,3> Matrix;
127 Matrix A; A.clear();
128 unsigned int nb = 0;
129 RealVector G = RealVector::zero; // center of gravity
130 for ( ConstIterator it = itB; it != itE; ++it )
131 {
132 G += RealVector( (*it)[ 0 ], (*it)[ 1 ], (*it)[ 2 ] );
133 ++nb;
134 }
135 G /= nb;
136 for ( ConstIterator it = itB; it != itE; ++it )
137 {
138 RealVector p( (*it)[ 0 ], (*it)[ 1 ], (*it)[ 2 ] );
139 p -= G;
140 for ( Dimension i = 0; i < 3; ++i )
141 for ( Dimension j = 0; j < 3; ++j )
142 A.setComponent( i, j, A( i, j ) + p[ i ] * p[ j ] );
143 }
144 // A is Gram matrix. We look for V such that V^t A V / |V|^2 is
145 // minimal. It is thus the first eigenvalue.
146 Matrix V;
147 RealVector values;
149 N = V.column( 0 ); // first eigenvector;
150 double mu = 0.0;
151 for ( ConstIterator it = itB; it != itE; ++it )
152 mu += N.dot( *it );
153 return mu/(double)nb;
154}
static void getEigenDecomposition(const Matrix &matrix, Matrix &eigenVectors, Vector &eigenValues)
Compute both eigen vectors and eigen values from an input matrix.
auto dot(const PointVector< dim, OtherComponent, OtherStorage > &v) const -> decltype(DGtal::dotProduct(*this, v))
Dot product with a PointVector.
Aim: implements basic MxN Matrix services (M,N>=1).
Space::RealVector RealVector
Definition StdDefs.h:171
DGtal::uint32_t Dimension
Definition Common.h:119

References DGtal::PointVector< dim, TEuclideanRing, TContainer >::dot(), DGtal::EigenDecomposition< TN, TComponent, TMatrix >::getEigenDecomposition(), and DGtal::PointVector< dim, double >::zero.

Referenced by main().

â—† 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 157 of file polyhedralizer.cpp.

158{
159 string inputFilename = argc > 1 ? argv[ 1 ] : examplesPath+"/samples/Al.100.vol";
160 int threshold = argc > 2 ? atoi( argv[ 2 ] ) : 0;
161 int widthNum = argc > 3 ? atoi( argv[ 3 ] ) : 2;
162 int widthDen = argc > 4 ? atoi( argv[ 4 ] ) : 1;
163
165 trace.beginBlock( "Reading vol file into an image." );
167 Image image = VolReader<Image>::importVol(inputFilename);
169 DigitalObject digitalObject( image, threshold );
170 trace.endBlock();
172
174 trace.beginBlock( "Construct the Khalimsky space from the image domain." );
175 KSpace ks;
176 bool space_ok = ks.init( image.domain().lowerBound(), image.domain().upperBound(), true );
177 if (!space_ok)
178 {
179 trace.error() << "Error in the Khamisky space construction."<<endl;
180 return 2;
181 }
182 trace.endBlock();
184
186 typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
187 MySurfelAdjacency surfAdj( false ); // exterior in all directions.
189
191 trace.beginBlock( "Extracting boundary by tracking the surface. " );
192 typedef KSpace::Surfel Surfel;
193 Surfel start_surfel = Surfaces<KSpace>::findABel( ks, digitalObject, 100000 );
197 MyContainer container( ks, digitalObject, surfAdj, start_surfel );
198 MyDigitalSurface digSurf( container );
199 trace.info() << "Digital surface has " << digSurf.size() << " surfels."
200 << endl;
201 trace.endBlock();
203
205 // First pass to find biggest planes.
206 trace.beginBlock( "Decomposition first pass. Computes all planes so as to sort vertices by the plane size." );
209 map<Surfel,unsigned int> v2size;
210 for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
211 v2size[ *it ] = 0;
212 int j = 0;
213 int nb = digSurf.size();
214 NaivePlaneComputer planeComputer;
215 vector<Point> layer;
216 vector<Surfel> layer_surfel;
217 for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
218 {
219 if ( ( (++j) % 50 == 0 ) || ( j == nb ) ) trace.progressBar( j, nb );
220 Surfel v = *it;
221 planeComputer.init( widthNum, widthDen );
222 // The visitor takes care of all the breadth-first traversal.
223 Visitor visitor( digSurf, v );
224 layer.clear();
225 layer_surfel.clear();
226 Visitor::Size currentSize = visitor.current().second;
227 while ( ! visitor.finished() )
228 {
229 Visitor::Node node = visitor.current();
230 v = node.first;
231 int axis = ks.sOrthDir( v );
232 Point p = ks.sCoords( ks.sDirectIncident( v, axis ) );
233 if ( node.second != currentSize )
234 {
235 bool isExtended = planeComputer.extend( layer.begin(), layer.end() );
236 if ( isExtended )
237 {
238 for ( vector<Surfel>::const_iterator it_layer = layer_surfel.begin(),
239 it_layer_end = layer_surfel.end(); it_layer != it_layer_end; ++it_layer )
240 {
241 ++v2size[ *it_layer ];
242 }
243 layer_surfel.clear();
244 layer.clear();
245 currentSize = node.second;
246 }
247 else
248 break;
249 }
250 layer_surfel.push_back( v );
251 layer.push_back( p );
252 visitor.expand();
253 }
254 }
255 // Prepare queue
256 typedef PairSorted2nd<Surfel,int> SurfelWeight;
257 priority_queue<SurfelWeight> Q;
258 for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
259 Q.push( SurfelWeight( *it, v2size[ *it ] ) );
260 trace.endBlock();
262
264 // Segmentation into planes
265 trace.beginBlock( "Decomposition second pass. Visits vertices from the one with biggest plane to the one with smallest plane." );
267 set<Surfel> processedVertices;
268 vector<RoundPlane*> roundPlanes;
269 map<Surfel,RoundPlane*> v2plane;
270 j = 0;
271 while ( ! Q.empty() )
272 {
273 if ( ( (++j) % 50 == 0 ) || ( j == nb ) ) trace.progressBar( j, nb );
274 Surfel v = Q.top().first;
275 Q.pop();
276 if ( processedVertices.find( v ) != processedVertices.end() ) // already in set
277 continue; // process to next vertex
278
279 RoundPlane* ptrRoundPlane = new RoundPlane;
280 roundPlanes.push_back( ptrRoundPlane ); // to delete them afterwards.
281 v2plane[ v ] = ptrRoundPlane;
282 ptrRoundPlane->first.init( widthNum, widthDen );
283 ptrRoundPlane->third = make_pair( RealVector::zero, 0.0 );
284 // The visitor takes care of all the breadth-first traversal.
285 Visitor visitor( digSurf, v );
286 layer.clear();
287 layer_surfel.clear();
288 Visitor::Size currentSize = visitor.current().second;
289 while ( ! visitor.finished() )
290 {
291 Visitor::Node node = visitor.current();
292 v = node.first;
293 Dimension axis = ks.sOrthDir( v );
294 Point p = ks.sCoords( ks.sDirectIncident( v, axis ) );
295 if ( node.second != currentSize )
296 {
297 bool isExtended = ptrRoundPlane->first.extend( layer.begin(), layer.end() );
298 if ( isExtended )
299 {
300 for ( vector<Surfel>::const_iterator it_layer = layer_surfel.begin(),
301 it_layer_end = layer_surfel.end(); it_layer != it_layer_end; ++it_layer )
302 {
303 Surfel s = *it_layer;
304 processedVertices.insert( s );
305 if ( v2plane.find( s ) == v2plane.end() )
306 v2plane[ s ] = ptrRoundPlane;
307 }
308 layer.clear();
309 layer_surfel.clear();
310 currentSize = node.second;
311 }
312 else break;
313 }
314 layer_surfel.push_back( v );
315 layer.push_back( p );
316 if ( processedVertices.find( v ) != processedVertices.end() )
317 // surfel is already in some plane.
318 visitor.ignore();
319 else
320 visitor.expand();
321 }
322 if ( visitor.finished() )
323 {
324 for ( vector<Surfel>::const_iterator it_layer = layer_surfel.begin(),
325 it_layer_end = layer_surfel.end(); it_layer != it_layer_end; ++it_layer )
326 {
327 Surfel s = *it_layer;
328 processedVertices.insert( s );
329 if ( v2plane.find( s ) == v2plane.end() )
330 v2plane[ s ] = ptrRoundPlane;
331 }
332 }
333 // Assign random color for each plane.
334 ptrRoundPlane->second = Color( rand() % 192 + 64, rand() % 192 + 64, rand() % 192 + 64, 255 );
335 }
336 trace.endBlock();
338
340 for ( vector<RoundPlane*>::iterator
341 it = roundPlanes.begin(), itE = roundPlanes.end();
342 it != itE; ++it )
343 {
344 NaivePlaneComputer& computer = (*it)->first;
345 RealVector normal;
346 double mu = LSF( normal, computer.begin(), computer.end() );
347 (*it)->third = make_pair( normal, mu );
348 }
350
352 map<Surfel, RealPoint> coordinates;
353 for ( map<Surfel,RoundPlane*>::const_iterator
354 it = v2plane.begin(), itE = v2plane.end();
355 it != itE; ++it )
356 {
357 Surfel v = it->first;
358 RoundPlane* rplane = it->second;
359 Point p = ks.sKCoords( v );
360 RealPoint rp( (double)p[ 0 ]/2.0, (double)p[ 1 ]/2.0, (double)p[ 2 ]/2.0 );
361 double mu = rplane->third.second;
362 RealVector normal = rplane->third.first;
363 double lambda = mu - rp.dot( normal );
364 coordinates[ v ] = rp + lambda*normal;
365 }
366 typedef vector<Surfel> SurfelRange;
367 map<Surfel, RealPoint> new_coordinates;
368 for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
369 {
370 Surfel s = *it;
371 SurfelRange neighbors;
372 back_insert_iterator<SurfelRange> writeIt = back_inserter( neighbors );
373 digSurf.writeNeighbors( writeIt, *it );
375 for ( SurfelRange::const_iterator its = neighbors.begin(), itsE = neighbors.end();
376 its != itsE; ++its )
377 x += coordinates[ *its ];
378 new_coordinates[ s ] = x / neighbors.size();
379 }
381
383 typedef unsigned int Number;
384 typedef Mesh<RealPoint> MyMesh;
385 typedef MyMesh::MeshFace MeshFace;
386 typedef MyDigitalSurface::FaceSet FaceSet;
388 map<Surfel, Number> index; // Numbers all vertices.
389 Number nbv = 0;
390 MyMesh polyhedron( true );
391 // Insert all projected surfels as vertices of the polyhedral surface.
392 for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
393 {
394 polyhedron.addVertex( new_coordinates[ *it ] );
395 index[ *it ] = nbv++;
396 }
397 // Define faces of the mesh. Outputs closed faces.
398 FaceSet faces = digSurf.allClosedFaces();
399 for ( typename FaceSet::const_iterator itf = faces.begin(), itf_end = faces.end();
400 itf != itf_end; ++itf )
401 {
402 MeshFace mface( itf->nbVertices );
403 VertexRange vtcs = digSurf.verticesAroundFace( *itf );
404 int i = 0;
405 for ( typename VertexRange::const_iterator itv = vtcs.begin(), itv_end = vtcs.end();
406 itv != itv_end; ++itv )
407 {
408 mface[ i++ ] = index[ *itv ];
409 }
410 polyhedron.addFace( mface, Color( 255, 243, 150, 255 ) );
411 }
413
415 typedef PolyscopeViewer<Space,KSpace> MyViewer3D;
416 MyViewer3D viewer( ks );
417 // bool isOK = polyhedron >> "test.off";
418 // bool isOK2 = polyhedron >> "test.obj";
419 viewer << polyhedron;
420 viewer.show();
422
424 for ( vector<RoundPlane*>::iterator
425 it = roundPlanes.begin(), itE = roundPlanes.end();
426 it != itE; ++it )
427 delete *it;
429 return true;
430 }
Aim: This class is useful to perform a breadth-first exploration of a graph given a starting point or...
void init(Dimension axis, InternalInteger diameter, InternalInteger widthNumerator=NumberTraits< InternalInteger >::ONE, InternalInteger widthDenominator=NumberTraits< InternalInteger >::ONE)
ConstIterator end() const
ConstIterator begin() const
bool extend(const Point &p)
Aim: A class that recognizes pieces of digital planes of given axis width. When the width is 1,...
Structure representing an RGB triple with alpha component.
Definition Color.h:77
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
DigitalSurfaceContainer::SurfelConstIterator ConstIterator
Aim: implements association bewteen points lying in a digital domain and values.
Definition Image.h:70
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as the boundary of an impl...
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
Dimension sOrthDir(const SCell &s) const
Given a signed surfel [s], returns its orthogonal direction (ie, the coordinate where the surfel is c...
Point sCoords(const SCell &c) const
Return its digital coordinates.
const Point & sKCoords(const SCell &c) const
Return its Khalimsky coordinates.
SCell sDirectIncident(const SCell &p, Dimension k) const
Return the direct incident cell of [p] along [k] (the incident cell along [k])
Aim: This class is defined to represent a surface mesh through a set of vertices and faces....
Definition Mesh.h:92
static Dimension size()
static SCell findABel(const KSpace &K, const PointPredicate &pp, unsigned int nbtries=1000)
Aim: Represent adjacencies between surfel elements, telling if it follows an interior to exterior ord...
Aim: Define a simple Foreground predicate thresholding image values given a single thresold....
DigitalSurface< MyDigitalSurfaceContainer > MyDigitalSurface
COBANaivePlaneComputer< Z3, InternalInteger > NaivePlaneComputer
BreadthFirstVisitor< MyDigitalSurface > Visitor
KhalimskySpaceND< 3, Integer > KSpace
Definition StdDefs.h:146
Space::RealPoint RealPoint
Definition StdDefs.h:170
Trace trace
double LSF(RealVector &N, ConstIterator itB, ConstIterator itE)
static ImageContainer importVol(const std::string &filename, const Functor &aFunctor=Functor())
[polyhedralizer-typedefs]
unsigned int index(DGtal::uint32_t n, unsigned int b)
Definition testBits.cpp:44
Image image(domain)
TriMesh::VertexRange VertexRange

References DGtal::DigitalSurface< TDigitalSurfaceContainer >::allClosedFaces(), DGtal::COBANaivePlaneComputer< TSpace, TInternalInteger >::begin(), DGtal::DigitalSurface< TDigitalSurfaceContainer >::begin(), DGtal::BreadthFirstVisitor< TGraph, TMarkSet >::current(), DGtal::PointVector< dim, TEuclideanRing, TContainer >::dot(), DGtal::COBANaivePlaneComputer< TSpace, TInternalInteger >::end(), DGtal::DigitalSurface< TDigitalSurfaceContainer >::end(), DGtal::BreadthFirstVisitor< TGraph, TMarkSet >::expand(), DGtal::COBANaivePlaneComputer< TSpace, TInternalInteger >::extend(), DGtal::Surfaces< TKSpace >::findABel(), DGtal::BreadthFirstVisitor< TGraph, TMarkSet >::finished(), DGtal::BreadthFirstVisitor< TGraph, TMarkSet >::ignore(), image(), DGtal::VolReader< TImageContainer, TFunctor >::importVol(), index(), DGtal::COBANaivePlaneComputer< TSpace, TInternalInteger >::init(), DGtal::KhalimskySpaceND< dim, TInteger >::init(), LSF(), DGtal::KhalimskySpaceND< dim, TInteger >::sCoords(), DGtal::KhalimskySpaceND< dim, TInteger >::sDirectIncident(), DGtal::DigitalSurface< TDigitalSurfaceContainer >::size(), DGtal::PointVector< dim, TEuclideanRing, TContainer >::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, double >::zero.