DGtal 1.4.2
Loading...
Searching...
No Matches
polyhedralizer.cpp
Go to the documentation of this file.
1
55#include <iostream>
56#include <vector>
57#include <set>
58#include <map>
59#include <queue>
60
61#include "DGtal/base/Common.h"
62#include "DGtal/helpers/StdDefs.h"
63#include "ConfigExamples.h"
65
67#include "DGtal/io/readers/VolReader.h"
68#include "DGtal/images/ImageSelector.h"
69#include "DGtal/images/ImageContainerBySTLVector.h"
70#include "DGtal/images/SimpleThresholdForegroundPredicate.h"
72
73#include "DGtal/io/Display3D.h"
74#include "DGtal/io/viewers/Viewer3D.h"
75#include "DGtal/io/DrawWithDisplay3DModifier.h"
76#include "DGtal/images/imagesSetsUtils/SetFromImage.h"
77
78#include "DGtal/topology/DigitalSurface.h"
79#include "DGtal/topology/helpers/Surfaces.h"
80#include "DGtal/topology/ImplicitDigitalSurface.h"
81
82#include "DGtal/graph/BreadthFirstVisitor.h"
83#include "DGtal/geometry/surfaces/COBANaivePlaneComputer.h"
84#include "DGtal/geometry/surfaces/ChordNaivePlaneComputer.h"
85#include "DGtal/geometry/surfaces/ChordGenericNaivePlaneComputer.h"
86
87#include "DGtal/math/linalg/SimpleMatrix.h"
88#include "DGtal/math/linalg/EigenDecomposition.h"
89
91
93using namespace std;
94using namespace DGtal;
95using namespace Z3i;
97
98template <typename T1, typename T2>
99struct PairSorted2nd
100{
101 typedef PairSorted2nd<T1,T2> Self;
102 inline PairSorted2nd( const T1& t1, const T2& t2 ) : first( t1 ), second( t2 ) {}
103 bool operator<( const Self& other ) const
104 {
105 return second < other.second;
106 }
107 T1 first;
108 T2 second;
109};
110
111template <typename T1, typename T2, typename T3>
112struct Triple
113{
114 T1 first;
115 T2 second;
116 T3 third;
117 Triple( T1 t1 = T1(), T2 t2 = T2(), T3 t3 = T3() )
118 : first( t1 ), second( t2 ), third( t3 )
119 {}
120};
121
122// Least-Square Fit of a plane N.x=mu on points [itB,itE). Returns mu.
123template <typename RealVector,
124 typename ConstIterator>
125double LSF( RealVector& N, ConstIterator itB, ConstIterator itE )
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}
157
158
159int main( int argc, char** argv )
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}
Aim: This class is useful to perform a breadth-first exploration of a graph given a starting point or...
const Node & current() const
std::pair< Vertex, Data > Node
FIXME.
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:68
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
std::vector< Vertex > VertexRange
The range of vertices is defined as a vector.
FaceSet allClosedFaces() const
std::set< Face > FaceSet
The set of faces is defined as set.
ConstIterator begin() const
ConstIterator end() const
VertexRange verticesAroundFace(const Face &f) const
DigitalSurfaceContainer::SurfelConstIterator ConstIterator
void writeNeighbors(OutputIterator &it, const Vertex &v) const
static void getEigenDecomposition(const Matrix &matrix, Matrix &eigenVectors, Vector &eigenValues)
Compute both eigen vectors and eigen values from an input matrix.
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...
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
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
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).
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...
void beginBlock(const std::string &keyword="")
std::ostream & error()
std::ostream & info()
void progressBar(const double currentValue, const double maximalValue)
double endBlock()
Aim: Define a simple Foreground predicate thresholding image values given a single thresold....
DigitalSurface< MyDigitalSurfaceContainer > MyDigitalSurface
COBANaivePlaneComputer< Z3, InternalInteger > NaivePlaneComputer
MyDigitalSurface::ConstIterator ConstIterator
BreadthFirstVisitor< MyDigitalSurface > Visitor
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
Definition Common.h:136
Trace trace
Definition Common.h:153
STL namespace.
double LSF(RealVector &N, ConstIterator itB, ConstIterator itE)
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
static ImageContainer importVol(const std::string &filename, const Functor &aFunctor=Functor())
int main()
Definition testBits.cpp:56
MyPointD Point
ImageContainerBySTLVector< Domain, Value > Image
TriMesh::VertexRange VertexRange