DGtalTools  1.2.0
3dImplicitSurfaceExtractorByThickening.cpp
1 
32 #include <iostream>
33 #include <map>
34 #include "CLI11.hpp"
35 
36 #include "DGtal/base/Common.h"
37 #include "DGtal/base/CountedPtr.h"
38 #include "DGtal/helpers/StdDefs.h"
39 #include "DGtal/math/MPolynomial.h"
40 #include "DGtal/io/readers/MPolynomialReader.h"
41 #include "DGtal/io/DrawWithDisplay3DModifier.h"
42 #include "DGtal/io/viewers/Viewer3D.h"
43 #include "DGtal/topology/KhalimskySpaceND.h"
44 #include "DGtal/topology/CubicalComplex.h"
45 #include "DGtal/topology/CubicalComplexFunctions.h"
46 #include "DGtal/topology/SetOfSurfels.h"
47 #include "DGtal/topology/DigitalSurface.h"
48 #include "DGtal/topology/helpers/Surfaces.h"
49 #include "DGtal/shapes/GaussDigitizer.h"
50 #include "DGtal/shapes/Mesh.h"
51 #include "DGtal/shapes/implicit/ImplicitPolynomial3Shape.h"
53 
54 using namespace std;
55 using namespace DGtal;
56 
57 
58 
117 template <typename CellOutputIterator, typename DigitalSurface>
118 void
119 analyzeSurface( CellOutputIterator itSure, CellOutputIterator itUnsure, DigitalSurface surface )
120 {
121  typedef typename DigitalSurface::KSpace KSpace;
122  typedef typename DigitalSurface::Surfel Surfel;
123  typedef typename DigitalSurface::ConstIterator ConstIterator;
124  typedef typename DigitalSurface::DigitalSurfaceContainer Container;
125  typedef typename DigitalSurface::DigitalSurfaceTracker Tracker;
126  typedef typename KSpace::Integer Integer;
127  typedef typename KSpace::Cell Cell;
128  const Dimension t = KSpace::dimension - 1;
129  const Container& C = surface.container();
130  const KSpace& K = surface.container().space();
131  Surfel s2;
132  for ( ConstIterator it = surface.begin(), itE = surface.end(); it != itE; ++it )
133  {
134  Surfel s = *it;
135  Cell is = K.unsigns( s );
136  Integer s_xt = K.sKCoord( s, t );
137  if ( s_xt == 0 ) *itUnsure++ = is;
138  else if ( s_xt == -1 )
139  {
140  CountedPtr<Tracker> tracker( C.newTracker( s ) );
141  if ( tracker->adjacent( s2, t, true ) != 0 )
142  {
143  Integer s2_xt = K.sKCoord( s2, t );
144  Cell ic = K.uIncident( is, t, true );
145  if ( s2_xt > 0 ) *itSure++ = ic;
146  else *itUnsure++ = ic;
147  }
148  }
149  else if ( s_xt == 1 )
150  {
151  CountedPtr<Tracker> tracker( C.newTracker( s ) );
152  if ( tracker->adjacent( s2, t, false ) != 0 )
153  {
154  Integer s2_xt = K.sKCoord( s2, t );
155  Cell ic = K.uIncident( is, t, false );
156  if ( s2_xt < 0 ) *itSure++ = ic;
157  else *itUnsure++ = ic;
158  }
159  }
160  else cout << " " << s_xt;
161  }
162 }
163 
172 template <typename ImplicitSurface, typename RealPoint>
173 RealPoint projectNewton( const ImplicitSurface & is,
174  RealPoint x,
175  typename RealPoint::Coordinate epsilon,
176  unsigned int max_iter )
177 {
178  typedef typename RealPoint::Coordinate Scalar;
179  RealPoint gx;
180  Scalar f, g2;
181  Scalar eps2 = epsilon * epsilon;
182  while ( max_iter-- != 0 )
183  {
184  f = is( x );
185  if ( abs( f ) < epsilon ) return x;
186  gx = is.gradient( x );
187  g2 = gx.dot( gx );
188  if ( g2 < eps2 ) return x;
189  x -= (f/g2) * gx;
190  }
191  return x;
192 }
193 
194 
195 template <typename CubicalComplex3, typename ImplicitShape3,
196  typename ImplicitDigitalShape3>
197 void projectComplex( std::vector< typename ImplicitShape3::RealPoint >& points,
198  const CubicalComplex3& complex3,
199  const ImplicitShape3& shape,
200  const ImplicitDigitalShape3& dshape,
201  double epsilon,
202  unsigned int max_iter,
203  double max_distance )
204 {
205  typedef typename CubicalComplex3::Cell Cell3;
206  typedef typename CubicalComplex3::Point Point3;
207  typedef typename CubicalComplex3::CellMapConstIterator CellMapConstIterator;
208  typedef typename ImplicitShape3::RealPoint RealPoint3;
209  typedef typename ImplicitShape3::Ring Ring;
210  points.clear();
211  for ( CellMapConstIterator it = complex3.begin( 0 ), itE = complex3.end( 0 ); it != itE; ++it )
212  {
213  Cell3 cell = it->first;
214  Point3 dp = complex3.space().uKCoords( cell ) - Point3::diagonal( 1 );
215  RealPoint3 p = dshape->embed( dp ) / 2.0;
216  RealPoint3 q = projectNewton( shape, p, epsilon, max_iter );
217  double d = (p-q).norm();
218  if ( d > max_distance ) q = p + (q-p)*( max_distance / d );
219  points.push_back( q );
220  }
221 }
222 
223 template <typename CubicalComplex3, typename ImplicitShape3,
224  typename ImplicitDigitalShape3>
225 typename ImplicitShape3::Ring
226 getValue( const CubicalComplex3& complex3,
227  const typename CubicalComplex3::Cell& cell,
228  const ImplicitShape3& shape,
229  const ImplicitDigitalShape3& dshape )
230 {
231  typedef typename CubicalComplex3::Cell Cell3;
232  typedef typename CubicalComplex3::Cells Cells3;
233  typedef typename CubicalComplex3::Point Point3;
234  typedef typename ImplicitShape3::RealPoint RealPoint3;
235  typedef typename ImplicitShape3::Ring Ring;
236 
237  Point3 dp = complex3.space().uKCoords( cell ) - Point3::diagonal( 1 );
238  RealPoint3 p = dshape->embed( dp ) / 2.0;
239  Ring v = shape( p );
240  return v;
241 }
242 
243 
244 template <typename CubicalComplex3, typename ImplicitShape3,
245  typename ImplicitDigitalShape3>
246 void doNotProjectComplex( std::vector< typename ImplicitShape3::RealPoint >& points,
247  const CubicalComplex3& complex3,
248  const ImplicitShape3& shape,
249  const ImplicitDigitalShape3& dshape )
250 {
251  typedef typename CubicalComplex3::Cell Cell3;
252  typedef typename CubicalComplex3::Point Point3;
253  typedef typename CubicalComplex3::CellMapConstIterator CellMapConstIterator;
254  typedef typename ImplicitShape3::RealPoint RealPoint3;
255  typedef typename ImplicitShape3::Ring Ring;
256  points.clear();
257  for ( CellMapConstIterator it = complex3.begin( 0 ), itE = complex3.end( 0 ); it != itE; ++it )
258  {
259  Cell3 cell = it->first;
260  Point3 dp = complex3.space().uKCoords( cell );
261  RealPoint3 p = dshape->embed( dp ) / 2.0;
262  points.push_back( p );
263  }
264 }
265 
266 
267 
268 int main( int argc, char** argv )
269 {
270  typedef int Integer;
271  typedef SpaceND<3,Integer> Space3;
272  typedef KhalimskySpaceND<3,Integer> KSpace3;
273  typedef KSpace3::Cell Cell3;
274  typedef std::map<Cell3, CubicalCellData> Map3;
275  // typedef boost::unordered_map<Cell3, CubicalCellData> Map3;
276  typedef CubicalComplex< KSpace3, Map3 > CC3;
277  typedef Space3::Point Point3;
278  typedef Space3::RealPoint RealPoint3;
279  typedef Space3::RealVector RealVector3;
280  typedef RealPoint3::Coordinate Ring;
281  typedef Ring Scalar;
282  typedef MPolynomial<3, Ring> Polynomial3;
283  typedef MPolynomialReader<3, Ring> Polynomial3Reader;
284  typedef ImplicitPolynomial3Shape<Space3> ImplicitShape3;
285  typedef GaussDigitizer< Space3, ImplicitShape3 >
286  ImplicitDigitalShape3;
287  typedef ImplicitDigitalShape3::Domain Domain3;
288  typedef CC3::CellMapIterator CellMapIterator;
289  typedef CC3::CellMapConstIterator CellMapConstIterator;
290  typedef CC3::Cells Cells3;
291 
292  // parse command line using CLI ----------------------------------------------
293  CLI::App app;
294  string poly_str;
295  std::string outputFileName {"result.raw"};
296  DGtal::int64_t rescaleInputMin {0};
297  DGtal::int64_t rescaleInputMax {255};
298  Ring min_x {-10.0};
299  Ring max_x {10.0};
300  Ring h {1.0};
301  Ring t {1e-2};
302  std::string project {"Newton"};
303  double epsilon {1e-6};
304  unsigned int max_iter {500};
305  std::string view {"Normal"};
306 
307  app.description( "Computes the zero level set of the given polynomial. Usage: 3dImplicitSurfaceExtractorByThickening -p <polynomial> [options]\n Example:\n 3dImplicitSurfaceExtractorByThickening -p \"x^2-y*z^2\" -g 0.1 -a -2 -A 2 -v Singular\n - whitney : x^2-y*z^2 \n - 4lines : x*y*(y-x)*(y-z*x) \n - cone : z^2-x^2-y^2 \n - simonU : x^2-z*y^2+x^4+y^4 \n - cayley3 : 4*(x^2 + y^2 + z^2) + 16*x*y*z - 1 \n - crixxi : -0.9*(y^2+z^2-1)^2-(x^2+y^2-1)^3 \n Some other examples (more difficult): \n 3dImplicitSurfaceExtractorByThickening -a -2 -A 2 -p \"((y^2+z^2-1)^2-(x^2+y^2-1)^3)*(y*(x-1)^2-z*(x+1))^2\" -g 0.025 -e 1e-6 -n 50000 -v Singular -t 0.5 -P Newton \n 3dImplicitSurfaceExtractorByThickening -a -2 -A 2 -p \"(x^5-4*z^3*y^2)*((x+y)^2-(z-x)^3)\" -g 0.025 -e 1e-6 -n 50000 -v Singular -t 0.05 -P Newton ");
308  app.add_option("-p,--polynomial,1", poly_str, "the implicit polynomial whose zero-level defines the shape of interest." )
309  ->required();
310  app.add_option("--minAABB,-a",min_x, "the min value of the AABB bounding box (domain)" , true);
311  app.add_option("--maxAABB,-A",max_x, "the max value of the AABB bounding box (domain)" , true);
312  app.add_option("--gridstep,-g", h, "the gridstep that defines the digitization (often called h).", true);
313  app.add_option("--thickness,-t",t, "the thickening parameter for the implicit surface." , true);
314  app.add_option("--project,-P", project, "defines the projection: either No or Newton.", true)
315  -> check(CLI::IsMember({"No", "Newton"}));
316  app.add_option("--epsilon,-e", epsilon, "the maximum precision relative to the implicit surface in the Newton approximation of F=0.", true);
317  app.add_option("--max_iter,-n", max_iter, "the maximum number of iteration in the Newton approximation of F=0.", true );
318  app.add_option("--view,-v", view, "specifies if the surface is viewed as is (Normal) or if places close to singularities are highlighted (Singular), or if unsure places should not be displayed (Hide).",true )
319  -> check(CLI::IsMember({"Singular", "Normal", "Hide"}));
320 
321  app.get_formatter()->column_width(40);
322  CLI11_PARSE(app, argc, argv);
323  // END parse command line using CLI ----------------------------------------------
324 
325 
326  //-------------- read polynomial and creating 3d implicit fct -----------------
327  trace.beginBlock( "Reading polynomial and creating 3D implicit function" );
328  Polynomial3 poly;
329  Polynomial3Reader reader;
330  string::const_iterator iter = reader.read( poly, poly_str.begin(), poly_str.end() );
331  if ( iter != poly_str.end() )
332  {
333  trace.error() << "ERROR reading polynomial: I read only <" << poly_str.substr( 0, iter - poly_str.begin() )
334  << ">, and I built P=" << poly << std::endl;
335  return 2;
336  }
337  // Creating implicit shape and storing it with smart pointer for automatic deallocation.
338  CountedPtr<ImplicitShape3> shape( new ImplicitShape3( poly ) );
339 
340  RealPoint3 p1( min_x, min_x, min_x );
341  RealPoint3 p2( max_x, max_x, max_x );
342  // Creating digitized shape and storing it with smart pointer for automatic deallocation.
343  CountedPtr<ImplicitDigitalShape3> dshape( new ImplicitDigitalShape3() );
344  dshape->attach( *shape );
345  dshape->init( p1, p2, RealVector3( h, h, h ) );
346  Domain3 domain3 = dshape->getDomain();
347  KSpace3 K3;
348  K3.init( domain3.lowerBound(), domain3.upperBound(), true );
349  trace.info() << "- domain is " << domain3 << std::endl;
350  trace.endBlock();
351 
352  //-------------- read polynomial and creating 3d implicit fct -----------------
353  trace.beginBlock( "Extracting thickened isosurface [-t,+t] of 3D polynomial. " );
354  CubicalCellData unsure_data( 0 );
355  CubicalCellData sure_data( CC3::FIXED );
356  CC3 complex3( K3 );
357  for ( Domain3::ConstIterator it = domain3.begin(), itE = domain3.end(); it != itE; ++it )
358  {
359  Cell3 spel = K3.uSpel( *it );
360  // RealPoint3 px = dshape->embed( *it );
361  RealPoint3 px = dshape->embed( K3.uKCoords( spel ) - Point3::diagonal( 1 ) ) / 2.0;
362  Ring s = (*shape)( px );
363  if ( (-t <= s ) && ( s <= t ) ) complex3.insertCell( spel, unsure_data );
364  }
365  trace.info() << "- K[-t,+t] = " << complex3 << endl;
366  complex3.close();
367  trace.info() << "- Cl K[-t,+t] = " << complex3 << endl;
368  std::vector<Cell3> separating_cells;
369  std::back_insert_iterator< std::vector<Cell3> > outItSurface( separating_cells );
370  Surfaces<KSpace3>::uWriteBoundary( outItSurface,
371  K3, *dshape, domain3.lowerBound(), domain3.upperBound() );
372  trace.info() << "- separating S = " << separating_cells.size() << " 2-cells." << endl;
373  complex3.insertCells( separating_cells.begin(), separating_cells.end(), sure_data );
374  for ( std::vector<Cell3>::const_iterator it = separating_cells.begin(), itE = separating_cells.end(); it != itE; ++it )
375  {
376  Cells3 bdry = K3.uFaces( *it );
377  for ( Cells3::const_iterator itBdry = bdry.begin(), itBdryE = bdry.end(); itBdry != itBdryE; ++itBdry )
378  complex3.insertCell( *itBdry, sure_data );
379  }
380  separating_cells.clear();
381  trace.info() << "- Cl K[-t,+t] + Cl S = " << complex3 << endl;
382  trace.endBlock();
383 
384  //-------------- Get boundary and inner cells --------------------------------
385  trace.beginBlock( "Get boundary and inner cells. " );
386  std::vector<Cell3> inner;
387  std::vector<Cell3> bdry;
388  functions::filterCellsWithinBounds
389  ( complex3, K3.uKCoords( K3.lowerCell() ), K3.uKCoords( K3.upperCell() ),
390  std::back_inserter( bdry ), std::back_inserter( inner ) );
391  trace.info() << "- there are " << inner.size() << " inner cells." << endl;
392  trace.info() << "- there are " << bdry.size() << " boundary cells." << endl;
393  trace.endBlock();
394 
395  //-------------- Compute priority function -----------------------------------
396  trace.beginBlock( "Compute priority function. " );
397  Dimension d = complex3.dim();
398  for ( Dimension i = 0; i <= d; ++i )
399  {
400  for ( CellMapIterator it = complex3.begin( i ), itE = complex3.end( i ); it != itE; ++it )
401  {
402  Cell3 cell = it->first;
403  Ring v = getValue( complex3, cell, *shape, dshape );
404  v = abs( 10000.0*v );
405  if ( v > 10000000.0 ) v = 10000000.0;
406  it->second.data &= ~CC3::VALUE;
407  it->second.data |= (DGtal::uint32_t) floor( v );
408  // std::cout << " " << it->second.data;
409  }
410  }
411  trace.endBlock();
412 
413  //-------------- Collapse boundary -------------------------------------------
414  trace.beginBlock( "Collapse boundary. " );
415  typename CC3::DefaultCellMapIteratorPriority priority;
416  CC3 bdry_complex3( K3 );
417  for ( std::vector<Cell3>::const_iterator it = bdry.begin(), itE = bdry.end(); it != itE; ++it )
418  {
419  Cell3 cell = *it;
420  Dimension d = K3.uDim( cell );
421  CellMapConstIterator cmIt = complex3.findCell( d, cell );
422  bdry_complex3.insertCell( d, cell, cmIt->second );
423  }
424  trace.info() << "- [before collapse] K_bdry =" << bdry_complex3 << endl;
425  functions::collapse( bdry_complex3, bdry.begin(), bdry.end(), priority, true, true, false );
426  trace.info() << "- [after collapse] K_bdry =" << bdry_complex3 << endl;
427  for ( std::vector<Cell3>::const_iterator it = bdry.begin(), itE = bdry.end(); it != itE; ++it )
428  {
429  Cell3 cell = *it;
430  Dimension d = K3.uDim( cell );
431  CellMapConstIterator cmIt = bdry_complex3.findCell( d, cell );
432  if ( cmIt != bdry_complex3.end( d ) ) {
433  CellMapIterator cmIt2 = complex3.findCell( d, cell );
434  cmIt2->second = sure_data;
435  }
436  }
437  trace.endBlock();
438 
439  //-------------- Collapse all -------------------------------------------
440  trace.beginBlock( "Collapse all. " );
441  std::copy( bdry.begin(), bdry.end(), std::back_inserter( inner ) );
442  functions::collapse( complex3, inner.begin(), inner.end(), priority, true, true, true );
443  trace.info() << "- K = " << complex3 << endl;
444  trace.endBlock();
445 
446  //-------------- Project complex onto surface --------------------------------
447  trace.beginBlock( "Project complex onto surface. " );
448  std::vector<RealPoint3> points;
449  if ( project == "Newton" )
450  projectComplex( points, complex3, *shape, dshape, epsilon, max_iter, h * sqrt(3.0));
451  else
452  doNotProjectComplex( points, complex3, *shape, dshape );
453  trace.endBlock();
454 
455  //-------------- Create Mesh -------------------------------------------
456  trace.beginBlock( "Create Mesh. " );
457  bool highlight = ( view == "Singular" );
458  bool hide = ( view == "Hide" );
459  Mesh<RealPoint3> mesh( true );
460  std::map<Cell3,unsigned int> indices;
461  int idx = 0;
462  for ( CellMapConstIterator it = complex3.begin( 0 ), itE = complex3.end( 0 ); it != itE; ++it, ++idx )
463  {
464  Cell3 cell = it->first;
465  indices[ cell ] = idx;
466  mesh.addVertex( points[ idx ] );
467  }
468  for ( CellMapConstIterator it = complex3.begin( 2 ), itE = complex3.end( 2 ); it != itE; ++it )
469  {
470  Cell3 cell = it->first;
471  bool fixed = it->second.data & CC3::FIXED;
472  Cells3 bdry = complex3.cellBoundary( cell, true );
473  std::vector<unsigned int> face_idx;
474  for ( Cells3::const_iterator itC = bdry.begin(), itCE = bdry.end(); itC != itCE; ++itC )
475  {
476  if ( complex3.dim( *itC ) == 0 )
477  face_idx.push_back( indices[ *itC ] );
478  }
479  if ( ( ! fixed ) && hide ) continue;
480  Color color = highlight
481  ? ( fixed ? Color::White : Color(128,255,128) )
482  : Color::White;
483  RealVector3 diag03 = points[ face_idx[ 0 ] ] - points[ face_idx[ 3 ] ];
484  RealVector3 diag12 = points[ face_idx[ 1 ] ] - points[ face_idx[ 2 ] ];
485  if ( diag03.dot( diag03 ) <= diag12.dot( diag12 ) )
486  {
487  mesh.addTriangularFace( face_idx[ 0 ], face_idx[ 1 ], face_idx[ 3 ], color );
488  mesh.addTriangularFace( face_idx[ 0 ], face_idx[ 3 ], face_idx[ 2 ], color );
489  }
490  else
491  {
492  mesh.addTriangularFace( face_idx[ 0 ], face_idx[ 1 ], face_idx[ 2 ], color );
493  mesh.addTriangularFace( face_idx[ 1 ], face_idx[ 3 ], face_idx[ 2 ], color );
494  }
495  //mesh.addQuadFace( face_idx[ 0 ], face_idx[ 1 ], face_idx[ 3 ], face_idx[ 2 ], color );
496  }
497  trace.endBlock();
498 
499  //-------------- View surface -------------------------------------------
500  QApplication application(argc,argv);
501  Viewer3D<Space3,KSpace3> viewer( K3 );
502  viewer.setWindowTitle("Implicit surface viewer by thickening");
503  viewer.show();
504  viewer << mesh;
505  // Display lines that are not in the mesh.
506  for ( CellMapConstIterator it = complex3.begin( 1 ), itE = complex3.end( 1 ); it != itE; ++it )
507  {
508  Cell3 cell = it->first;
509  bool fixed = it->second.data & CC3::FIXED;
510  std::vector<Cell3> dummy;
511  std::back_insert_iterator< std::vector<Cell3> > outIt( dummy );
512  complex3.directCoFaces( outIt, cell );
513  if ( ! dummy.empty() ) continue;
514 
515  Cells3 bdry = complex3.cellBoundary( cell, true );
516  Cell3 v0 = *(bdry.begin() );
517  Cell3 v1 = *(bdry.begin() + 1);
518  if ( ( ! fixed ) && hide ) continue;
519  Color color = highlight
520  ? ( fixed ? Color::White : Color(128,255,128) )
521  : Color::White;
522  viewer.setLineColor( color );
523  viewer.addLine( points[ indices[ v0 ] ], points[ indices[ v1 ] ], h/2.0 );
524  }
525  viewer << Viewer3D<Space3,KSpace3>::updateDisplay;
526  return application.exec();
527 }
Definition: ATu0v1.h:57