DGtalTools  0.9.4
3dImplicitSurfaceExtractorBy4DExtension.cpp
1 
32 #include <iostream>
34 #include <map>
35 #include <boost/program_options/options_description.hpp>
36 #include <boost/program_options/parsers.hpp>
37 #include <boost/program_options/variables_map.hpp>
38 
39 #include "DGtal/base/Common.h"
40 #include "DGtal/base/CountedPtr.h"
41 #include "DGtal/helpers/StdDefs.h"
42 #include "DGtal/math/MPolynomial.h"
43 #include "DGtal/io/readers/MPolynomialReader.h"
44 #include "DGtal/io/DrawWithDisplay3DModifier.h"
45 #include "DGtal/io/viewers/Viewer3D.h"
46 #include "DGtal/topology/KhalimskySpaceND.h"
47 #include "DGtal/topology/CubicalComplex.h"
48 #include "DGtal/topology/CubicalComplexFunctions.h"
49 #include "DGtal/topology/SetOfSurfels.h"
50 #include "DGtal/topology/DigitalSurface.h"
51 #include "DGtal/topology/helpers/Surfaces.h"
52 #include "DGtal/shapes/GaussDigitizer.h"
53 #include "DGtal/shapes/Mesh.h"
54 #include "DGtal/shapes/implicit/ImplicitPolynomial3Shape.h"
55 
57 
58 using namespace std;
59 using namespace DGtal;
60 
130 template <typename TSpace3, typename TSpace4, typename TPolynomial3>
131 struct ImplicitSurface4DProductExtension {
132  typedef TSpace3 Space3;
133  typedef TSpace4 Space4;
134  typedef TPolynomial3 Polynomial3;
135  typedef typename Space3::RealPoint RealPoint3;
136  typedef typename Space3::RealVector RealVector3;
137  typedef typename Space3::Integer Integer;
138  typedef typename RealPoint3::Coordinate Ring;
139  typedef typename Space4::RealPoint RealPoint4;
140  typedef typename Space4::RealVector RealVector4;
141  typedef RealPoint4 RealPoint;
142  typedef RealVector4 RealVector;
143 
144  ImplicitSurface4DProductExtension( Clone<Polynomial3> poly )
145  : f( poly )
146  {
147  fx = derivative<0>( f );
148  fy = derivative<1>( f );
149  fz = derivative<2>( f );
150  fxx = derivative<0>( fx );
151  fxy = derivative<0>( fy );
152  fxz = derivative<0>( fz );
153  fyy = derivative<1>( fy );
154  fyz = derivative<1>( fz );
155  fzz = derivative<2>( fz );
156  }
157 
162  Ring operator()(const RealPoint &aPoint) const
163  {
164  Ring x = aPoint[ 0 ];
165  Ring y = aPoint[ 1 ];
166  Ring z = aPoint[ 2 ];
167  Ring t = aPoint[ 3 ];
168  RealVector3 grad( fx(x)(y)(z), fy(x)(y)(z), fz(x)(y)(z) );
169  return f(x)(y)(z) - t * grad.norm();
170  }
171 
176  bool isInside(const RealPoint &aPoint) const
177  {
178  return this->operator()( aPoint ) <= NumberTraits<Ring>::ZERO;
179  }
180 
187  Orientation orientation(const RealPoint &aPoint) const
188  {
189  Ring v = this->operator()( aPoint );
190  if ( v > NumberTraits<Ring>::ZERO ) return OUTSIDE;
191  else if ( v < NumberTraits<Ring>::ZERO ) return INSIDE;
192  else return ON;
193  }
194 
199  inline
200  RealVector gradient( const RealPoint &aPoint ) const
201  {
202  Ring x = aPoint[ 0 ];
203  Ring y = aPoint[ 1 ];
204  Ring z = aPoint[ 2 ];
205  Ring t = aPoint[ 3 ];
206  RealVector3 grad( fx(x)(y)(z), fy(x)(y)(z), fz(x)(y)(z) );
207  Ring ngrad = grad.norm();
208  grad /= ngrad;
209  RealVector d_gradx( fxx(x)(y)(z), fxy(x)(y)(z),fxz(x)(y)(z) );
210  RealVector d_grady( fxy(x)(y)(z), fyy(x)(y)(z),fyz(x)(y)(z) );
211  RealVector d_gradz( fxz(x)(y)(z), fyz(x)(y)(z),fzz(x)(y)(z) );
212  return RealVector( grad[ 0 ] - t * d_gradx.dot( grad ),
213  grad[ 1 ] - t * d_grady.dot( grad ),
214  grad[ 2 ] - t * d_gradz.dot( grad ),
215  -ngrad );
216 
217  }
218 
219 private:
220  Polynomial3 f;
221  Polynomial3 fx;
222  Polynomial3 fy;
223  Polynomial3 fz;
224  Polynomial3 fxx;
225  Polynomial3 fxy;
226  Polynomial3 fxz;
227  Polynomial3 fyy;
228  Polynomial3 fyz;
229  Polynomial3 fzz;
230 };
231 
236 template <typename TSpace3, typename TSpace4, typename TPolynomial3>
237 struct ImplicitSurface4DExtension {
238  typedef TSpace3 Space3;
239  typedef TSpace4 Space4;
240  typedef TPolynomial3 Polynomial3;
241  typedef typename Space3::RealPoint RealPoint3;
242  typedef typename Space3::RealVector RealVector3;
243  typedef typename Space3::Integer Integer;
244  typedef typename RealPoint3::Coordinate Ring;
245  typedef typename Space4::RealPoint RealPoint4;
246  typedef typename Space4::RealVector RealVector4;
247  typedef RealPoint4 RealPoint;
248  typedef RealVector4 RealVector;
249 
250  ImplicitSurface4DExtension( Clone<Polynomial3> poly )
251  : f( poly )
252  {
253  fx = derivative<0>( f );
254  fy = derivative<1>( f );
255  fz = derivative<2>( f );
256  fxx = derivative<0>( fx );
257  fxy = derivative<0>( fy );
258  fxz = derivative<0>( fz );
259  fyy = derivative<1>( fy );
260  fyz = derivative<1>( fz );
261  fzz = derivative<2>( fz );
262  }
263 
268  Ring operator()(const RealPoint &aPoint) const
269  {
270  Ring x = aPoint[ 0 ];
271  Ring y = aPoint[ 1 ];
272  Ring z = aPoint[ 2 ];
273  Ring t = aPoint[ 3 ];
274  // RealVector3 grad( fx(x)(y)(z), fy(x)(y)(z), fz(x)(y)(z) );
275  return f(x)(y)(z) - t;
276  }
277 
282  bool isInside(const RealPoint &aPoint) const
283  {
284  return this->operator()( aPoint ) <= NumberTraits<Ring>::ZERO;
285  }
286 
293  Orientation orientation(const RealPoint &aPoint) const
294  {
295  Ring v = this->operator()( aPoint );
296  if ( v > NumberTraits<Ring>::ZERO ) return OUTSIDE;
297  else if ( v < NumberTraits<Ring>::ZERO ) return INSIDE;
298  else return ON;
299  }
300 
305  inline
306  RealVector gradient( const RealPoint &aPoint ) const
307  {
308  Ring x = aPoint[ 0 ];
309  Ring y = aPoint[ 1 ];
310  Ring z = aPoint[ 2 ];
311  Ring t = aPoint[ 3 ];
312  return RealVector( fx(x)(y)(z), fy(x)(y)(z), fz(x)(y)(z), -1.0 );
313  }
314 
315 private:
316  Polynomial3 f;
317  Polynomial3 fx;
318  Polynomial3 fy;
319  Polynomial3 fz;
320  Polynomial3 fxx;
321  Polynomial3 fxy;
322  Polynomial3 fxz;
323  Polynomial3 fyy;
324  Polynomial3 fyz;
325  Polynomial3 fzz;
326 };
327 
328 
329 template <typename CellOutputIterator, typename DigitalSurface>
330 void
331 analyzeSurface( CellOutputIterator itSure, CellOutputIterator itUnsure,
332  DigitalSurface surface )
333 {
334  typedef typename DigitalSurface::KSpace KSpace;
335  typedef typename DigitalSurface::Surfel Surfel;
336  typedef typename DigitalSurface::ConstIterator ConstIterator;
337  typedef typename DigitalSurface::DigitalSurfaceContainer Container;
338  typedef typename DigitalSurface::DigitalSurfaceTracker Tracker;
339  typedef typename KSpace::Integer Integer;
340  typedef typename KSpace::Cell Cell;
341  const Dimension t = KSpace::dimension - 1;
342  const Container& C = surface.container();
343  const KSpace& K = surface.container().space();
344  Surfel s2;
345  for ( ConstIterator it = surface.begin(), itE = surface.end(); it != itE; ++it )
346  {
347  Surfel s = *it;
348  Cell is = K.unsigns( s );
349  Integer s_xt = K.sKCoord( s, t );
350  if ( s_xt == 0 ) *itUnsure++ = is;
351  else if ( s_xt == -1 )
352  {
353  CountedPtr<Tracker> tracker( C.newTracker( s ) );
354  if ( tracker->adjacent( s2, t, true ) != 0 )
355  {
356  Integer s2_xt = K.sKCoord( s2, t );
357  Cell ic = K.uIncident( is, t, true );
358  if ( s2_xt > 0 ) *itSure++ = ic;
359  else *itUnsure++ = ic;
360  }
361  }
362  else if ( s_xt == 1 )
363  {
364  CountedPtr<Tracker> tracker( C.newTracker( s ) );
365  if ( tracker->adjacent( s2, t, false ) != 0 )
366  {
367  Integer s2_xt = K.sKCoord( s2, t );
368  Cell ic = K.uIncident( is, t, false );
369  if ( s2_xt < 0 ) *itSure++ = ic;
370  else *itUnsure++ = ic;
371  }
372  }
373  else cout << " " << s_xt;
374  }
375 }
376 
385 template <typename ImplicitSurface, typename RealPoint>
386 RealPoint projectNewton( const ImplicitSurface & is,
387  RealPoint x,
388  typename RealPoint::Coordinate epsilon,
389  unsigned int max_iter )
390 {
391  typedef typename RealPoint::Coordinate Scalar;
392  RealPoint gx;
393  Scalar f, g2;
394  Scalar eps2 = epsilon * epsilon;
395  while ( max_iter-- != 0 )
396  {
397  f = is( x );
398  if ( abs( f ) < epsilon ) return x;
399  gx = is.gradient( x );
400  g2 = gx.dot( gx );
401  if ( g2 < eps2 ) return x;
402  x -= (f/g2) * gx;
403  }
404  return x;
405 }
406 
407 
408 template <typename CubicalComplex4, typename ImplicitShape4,
409  typename ImplicitDigitalShape4, typename ImplicitShape3>
410 void projectComplex( std::vector< typename ImplicitShape3::RealPoint >& points,
411  const CubicalComplex4& complex4,
412  const ImplicitShape4& shape,
413  const ImplicitDigitalShape4& dshape,
414  const ImplicitShape3& shape3,
415  double epsilon,
416  unsigned int max_iter )
417 {
418  typedef typename CubicalComplex4::Cell Cell4;
419  typedef typename CubicalComplex4::Point Point4;
420  typedef typename CubicalComplex4::CellMapConstIterator CellMapConstIterator;
421  typedef typename ImplicitShape4::RealPoint RealPoint4;
422  typedef typename ImplicitShape4::Ring Ring;
423  typedef typename ImplicitShape3::RealPoint RealPoint3;
424  points.clear();
425  for ( CellMapConstIterator it = complex4.begin( 0 ), itE = complex4.end( 0 ); it != itE; ++it )
426  {
427  Cell4 cell = it->first;
428  Point4 dp = complex4.space().uCoords( cell );
429  RealPoint4 p = dshape->embed( dp );
430  RealPoint3 p3 = RealPoint3( p[ 0 ], p[ 1 ], p[ 2 ] );
431  RealPoint3 q = projectNewton( shape3, p3, epsilon, max_iter );
432  points.push_back( q );
433  }
434 }
435 
436 template <typename CubicalComplex4, typename ImplicitShape4,
437  typename ImplicitDigitalShape4, typename ImplicitShape3>
438 void doNotProjectComplex( std::vector< typename ImplicitShape3::RealPoint >& points,
439  const CubicalComplex4& complex4,
440  const ImplicitShape4& shape,
441  const ImplicitDigitalShape4& dshape,
442  const ImplicitShape3& shape3 )
443 {
444  typedef typename CubicalComplex4::Cell Cell4;
445  typedef typename CubicalComplex4::Point Point4;
446  typedef typename CubicalComplex4::CellMapConstIterator CellMapConstIterator;
447  typedef typename ImplicitShape4::RealPoint RealPoint4;
448  typedef typename ImplicitShape4::Ring Ring;
449  typedef typename ImplicitShape3::RealPoint RealPoint3;
450  for ( CellMapConstIterator it = complex4.begin( 0 ), itE = complex4.end( 0 ); it != itE; ++it )
451  {
452  Cell4 cell = it->first;
453  Point4 dp = complex4.space().uCoords( cell );
454  RealPoint4 p = dshape->embed( dp );
455  RealPoint3 p3 = RealPoint3( p[ 0 ], p[ 1 ], p[ 2 ] );
456  points.push_back( p3 );
457  }
458 }
459 
460 
461 int main( int argc, char** argv )
462 {
463  typedef int Integer;
464  typedef SpaceND<3,Integer> Space3;
465  typedef KhalimskySpaceND<3,Integer> KSpace3;
466  typedef KSpace3::Cell Cell3;
467  typedef std::map<Cell3, CubicalCellData> Map3;
468  typedef CubicalComplex< KSpace3, Map3 > CC3;
469  typedef Space3::Point Point3;
470  typedef Space3::RealPoint RealPoint3;
471  typedef RealPoint3::Coordinate Ring;
472  typedef Ring Scalar;
473  typedef MPolynomial<3, Ring> Polynomial3;
474  typedef MPolynomialReader<3, Ring> Polynomial3Reader;
475  typedef ImplicitPolynomial3Shape<Space3> ImplicitShape3;
476  typedef SpaceND<4,Integer> Space4;
477  typedef KhalimskySpaceND<4,Integer> KSpace4;
478  typedef KSpace4::Cell Cell4;
479  typedef KSpace4::Cells Cells4;
480  typedef KSpace4::SCell SCell4;
481  typedef Space4::Point Point4;
482  typedef Space4::RealPoint RealPoint4;
483  typedef Space4::RealVector RealVector4;
484  typedef ImplicitSurface4DProductExtension<Space3,Space4,Polynomial3>
485  ImplicitShape4;
486  typedef GaussDigitizer< Space4, ImplicitShape4 >
487  ImplicitDigitalShape4;
488  typedef ImplicitDigitalShape4::Domain Domain4;
489  typedef std::map<Cell4, CubicalCellData> Map4;
490  typedef CubicalComplex< KSpace4, Map4 > CC4;
491  typedef CC4::CellMapIterator CellMapIterator;
492  typedef CC4::CellMapConstIterator CellMapConstIterator;
493  typedef CC4::Cells Cells4;
494 
495  //-------------- parse command line ----------------------------------------------
496  namespace po = boost::program_options;
497  po::options_description general_opt("Allowed options are: ");
498  general_opt.add_options()
499  ("help,h", "display this message")
500  ("polynomial,p", po::value<string>(), "the implicit polynomial whose zero-level defines the shape of interest." )
501  ("minAABB,a", po::value<double>()->default_value( -10.0 ), "the min value of the AABB bounding box (domain)" )
502  ("maxAABB,A", po::value<double>()->default_value( 10.0 ), "the max value of the AABB bounding box (domain)" )
503  ("gridstep,g", po::value< double >()->default_value( 1.0 ), "the gridstep that defines the digitization (often called h). " )
504  ("timestep,t", po::value< double >()->default_value( 0.000001 ), "the gridstep that defines the digitization in the 4th dimension (small is generally a good idea, default is 1e-6). " )
505  ("project,P", po::value< std::string >()->default_value( "Newton" ), "defines the projection: either No or Newton." )
506  ("epsilon,e", po::value< double >()->default_value( 0.0000001 ), "the maximum precision relative to the implicit surface." )
507  ("max_iter,n", po::value< unsigned int >()->default_value( 500 ), "the maximum number of iteration in the Newton approximation of F=0." )
508  ("view,v", po::value< std::string >()->default_value( "Normal" ), "specifies if the surface is viewed as is (Normal) or if places close to singularities are highlighted (Singular)." )
509  ;
510 
511  bool parseOK=true;
512  po::variables_map vm;
513  try {
514  po::store(po::parse_command_line(argc, argv, general_opt), vm);
515  } catch ( const exception& ex ) {
516  parseOK=false;
517  cerr << "Error checking program options: "<< ex.what()<< endl;
518  }
519  po::notify(vm);
520  if( ! parseOK || vm.count("help") || ! vm.count( "polynomial" ) )
521  {
522  if ( ! vm.count( "polynomial" ) )
523  cerr << "Need parameter --polynomial / -p" << endl;
524 
525  cerr << "Usage: " << argv[0] << " -p <polynomial> [options]\n"
526  << "Computes the zero level set of the given polynomial."
527  << endl
528  << general_opt << "\n";
529  cerr << "Example:\n" << endl
530  << argv[ 0 ] << " -p \"-0.9*(y^2+z^2-1)^2-(x^2+y^2-1)^3\" -g 0.06125 -a -2 -A 2 -v Singular -t 0.02" << endl
531  << " - whitney : x^2-y*z^2" << endl
532  << " - 4lines : x*y*(y-x)*(y-z*x)" << endl
533  << " - cone : z^2-x^2-y^2" << endl
534  << " - simonU : x^2-z*y^2+x^4+y^4" << endl
535  << " - cayley3 : 4*(x^2 + y^2 + z^2) + 16*x*y*z - 1" << endl
536  << " - crixxi : -0.9*(y^2+z^2-1)^2-(x^2+y^2-1)^3" << endl << endl;
537  return 0;
538  }
539 
540  //-------------- read polynomial and creating 4d implicit fct -----------------
541  trace.beginBlock( "Reading polynomial and creating 4D implicit function" );
542  string poly_str = vm[ "polynomial" ].as<string>();
543  Polynomial3 poly;
544  Polynomial3Reader reader;
545  string::const_iterator iter = reader.read( poly, poly_str.begin(), poly_str.end() );
546  if ( iter != poly_str.end() )
547  {
548  trace.error() << "ERROR reading polynomial: I read only <" << poly_str.substr( 0, iter - poly_str.begin() )
549  << ">, and I built P=" << poly << std::endl;
550  return 2;
551  }
552  // Creating implicit shape and storing it with smart pointer for automatic deallocation.
553  ImplicitShape3 shape3( poly );
554  CountedPtr<ImplicitShape4> shape( new ImplicitShape4( poly ) );
555 
556  Ring min_x = vm[ "minAABB" ].as<double>();
557  Ring max_x = vm[ "maxAABB" ].as<double>();
558  Ring h = vm[ "gridstep" ].as<double>();
559  Ring t = vm[ "timestep" ].as<double>();
560  RealPoint4 p1( min_x, min_x, min_x, -t*h );
561  RealPoint4 p2( max_x, max_x, max_x, 0 );
562  // Creating digitized shape and storing it with smart pointer for automatic deallocation.
563  CountedPtr<ImplicitDigitalShape4> dshape( new ImplicitDigitalShape4() );
564  dshape->attach( *shape );
565  dshape->init( p1, p2, RealVector4( h, h, h, t ) );
566  Domain4 domain4 = dshape->getDomain();
567  KSpace4 K4;
568  K4.init( domain4.lowerBound(), domain4.upperBound(), true );
569  trace.info() << "- domain is " << domain4 << std::endl;
570  trace.endBlock();
571 
572  //-------------- read polynomial and creating 4d implicit fct -----------------
573  trace.beginBlock( "Extracting isosurface 0 of 4D polynomial. " );
574  CC4 complex4( K4 );
575  std::vector<SCell4> surface_cells;
576  std::back_insert_iterator< std::vector<SCell4> > outItSurface( surface_cells );
577  Surfaces<KSpace4>::sWriteBoundary( outItSurface,
578  K4, *dshape, domain4.lowerBound(), domain4.upperBound() );
579  trace.info() << "- 4D surface has " << surface_cells.size() << " 3-cells." << endl;
580  typedef SetOfSurfels<KSpace4> SurfaceContainer;
581  SurfaceContainer* ptrF0 = new SurfaceContainer( K4, SurfelAdjacency< KSpace4::dimension >( false ) );
582  ptrF0->surfelSet().insert( surface_cells.begin(), surface_cells.end() );
583  DigitalSurface<SurfaceContainer> surface( ptrF0 );
584  std::vector<Cell4> sure_cells;
585  std::vector<Cell4> unsure_cells;
586  analyzeSurface( std::back_inserter( sure_cells ), std::back_inserter( unsure_cells ),
587  surface );
588  trace.info() << "- sure cells = " << sure_cells.size() << endl;
589  trace.info() << "- unsure cells = " << unsure_cells.size() << endl;
590  CubicalCellData unsure_data( 0 );
591  CubicalCellData sure_data( CC4::FIXED );
592  complex4.insertCells( sure_cells.begin(), sure_cells.end(), sure_data );
593  complex4.insertCells( unsure_cells.begin(), unsure_cells.end(), unsure_data );
594  surface_cells.clear();
595  trace.endBlock();
596 
597  //-------------- Get boundary and inner cells --------------------------------
598  trace.beginBlock( "Get boundary and inner cells. " );
599  std::vector<Cell4> inner;
600  std::vector<Cell4> bdry;
601  complex4.close();
602  trace.info() << "- K=" << complex4 << endl;
603  functions::filterCellsWithinBounds
604  ( complex4, K4.uKCoords( K4.lowerCell() ), K4.uKCoords( K4.upperCell() ),
605  std::back_inserter( bdry ), std::back_inserter( inner ) );
606  trace.info() << "- there are " << inner.size() << " inner cells." << endl;
607  trace.info() << "- there are " << bdry.size() << " boundary cells." << endl;
608  trace.endBlock();
609 
610  //-------------- Compute priority function -----------------------------------
611  trace.beginBlock( "Compute priority function. " );
612  Dimension d = complex4.dim();
613  for ( Dimension i = 0; i <= d; ++i )
614  {
615  for ( CellMapIterator it = complex4.begin( i ), itE = complex4.end( i ); it != itE; ++it )
616  {
617  Cell4 cell = it->first;
618  Point4 dp = K4.uCoords( cell );
619  RealPoint4 p = dshape->embed( dp );
620  Ring v = (*shape)( p );
621  v = abs( 1000.0*v );
622  if ( v > 1000000.0 ) v = 1000000.0;
623  it->second.data &= ~CC4::VALUE;
624  it->second.data |= (DGtal::uint32_t) floor( v );
625  // std::cout << " " << it->second.data;
626  }
627  }
628  trace.endBlock();
629 
630  //-------------- Collapse boundary -------------------------------------------
631  trace.beginBlock( "Collapse boundary. " );
632  typename CC4::DefaultCellMapIteratorPriority priority;
633  CC4 bdry_complex4( K4 );
634  for ( std::vector<Cell4>::const_iterator it = bdry.begin(), itE = bdry.end(); it != itE; ++it )
635  {
636  Cell4 cell = *it;
637  Dimension d = K4.uDim( cell );
638  CellMapConstIterator cmIt = complex4.findCell( d, cell );
639  bdry_complex4.insertCell( d, cell, cmIt->second );
640  }
641  //bdry_complex4.insertCells( bdry.begin(), bdry.end() );
642  trace.info() << "- [before collapse] K_bdry =" << bdry_complex4 << endl;
643  functions::collapse( bdry_complex4, bdry.begin(), bdry.end(), priority, true, true, true );
644  trace.info() << "- [after collapse] K_bdry =" << bdry_complex4 << endl;
645  for ( std::vector<Cell4>::const_iterator it = bdry.begin(), itE = bdry.end(); it != itE; ++it )
646  {
647  Cell4 cell = *it;
648  Dimension d = K4.uDim( cell );
649  CellMapConstIterator cmIt = bdry_complex4.findCell( d, cell );
650  if ( cmIt != bdry_complex4.end( d ) ) {
651  CellMapIterator cmIt2 = complex4.findCell( d, cell );
652  cmIt2->second = sure_data;
653  }
654  }
655  trace.endBlock();
656 
657  //-------------- Collapse all -------------------------------------------
658  trace.beginBlock( "Collapse all. " );
659  std::copy( bdry.begin(), bdry.end(), std::back_inserter( inner ) );
660  //typename CC4::DefaultCellMapIteratorPriority priority;
661  functions::collapse( complex4, inner.begin(), inner.end(), priority, true, true, true );
662  trace.info() << "- K=" << complex4 << endl;
663  trace.endBlock();
664 
665  //-------------- Project complex onto surface --------------------------------
666  trace.beginBlock( "Project complex onto surface. " );
667  std::string project = vm[ "project" ].as<std::string>();
668  double epsilon = vm[ "epsilon" ].as<double>();
669  unsigned int max_iter = vm[ "max_iter" ].as<unsigned int>();
670  std::vector<RealPoint3> points;
671  if ( project == "Newton" )
672  projectComplex( points, complex4, *shape, dshape, shape3, epsilon, max_iter );
673  else
674  doNotProjectComplex( points, complex4, *shape, dshape, shape3 );
675  trace.endBlock();
676 
677  //-------------- Create Mesh -------------------------------------------
678  trace.beginBlock( "Create Mesh. " );
679  std::string view = vm[ "view" ].as<std::string>();
680  bool highlight = ( view == "Singular" );
681  Mesh<RealPoint3> mesh( true );
682  std::map<Cell4,unsigned int> indices;
683  int idx = 0;
684  for ( CellMapConstIterator it = complex4.begin( 0 ), itE = complex4.end( 0 ); it != itE; ++it, ++idx )
685  {
686  Cell4 cell = it->first;
687  indices[ cell ] = idx;
688  mesh.addVertex( points[ idx ] );
689  }
690  for ( CellMapConstIterator it = complex4.begin( 2 ), itE = complex4.end( 2 ); it != itE; ++it, ++idx )
691  {
692  Cell4 cell = it->first;
693  bool fixed = it->second.data & CC4::FIXED;
694  Cells4 bdry = complex4.cellBoundary( cell, true );
695  std::vector<unsigned int> face_idx;
696  for ( Cells4::const_iterator itC = bdry.begin(), itCE = bdry.end(); itC != itCE; ++itC )
697  {
698  if ( complex4.dim( *itC ) == 0 )
699  face_idx.push_back( indices[ *itC ] );
700  }
701  Color color = highlight
702  ? ( fixed ? Color::White : Color(128,255,128) )
703  : Color::White;
704  mesh.addQuadFace( face_idx[ 0 ], face_idx[ 1 ], face_idx[ 3 ], face_idx[ 2 ], color );
705  }
706  trace.endBlock();
707 
708  //-------------- View surface -------------------------------------------
709  QApplication application(argc,argv);
710  Point4 low4 = K4.lowerBound();
711  Point4 up4 = K4.upperBound();
712  KSpace3 K3;
713  K3.init( Point3( low4[ 0 ], low4[ 1 ], low4[ 2 ] ),
714  Point3( up4 [ 0 ], up4 [ 1 ], up4 [ 2 ] ), true );
715  Viewer3D<Space3,KSpace3> viewer( K3 );
716  viewer.setWindowTitle("Implicit surface viewer by 4d extension");
717  viewer.show();
718  viewer << mesh;
719  viewer.setLineColor( highlight ? Color::Red : Color( 120, 120, 120 ) );
720  // Drawing lines
721  for ( CellMapConstIterator it = complex4.begin( 1 ), itE = complex4.end( 1 ); it != itE; ++it )
722  {
723  Cell4 cell = it->first;
724  std::vector<Cell4> dummy;
725  std::back_insert_iterator< std::vector<Cell4> > outIt1( dummy );
726  complex4.directCoFaces( outIt1, cell );
727  if ( ! dummy.empty() ) continue;
728  Cells4 vertices = complex4.cellBoundary( cell );
729  ASSERT( vertices.size() == 2 );
730  RealPoint3 p1 = points[ indices[ vertices.front() ] ];
731  RealPoint3 p2 = points[ indices[ vertices.back() ] ];
732  viewer.addLine( p1, p2, 0.05 );
733  }
734  viewer << Viewer3D<Space3,KSpace3>::updateDisplay;
735  return application.exec();
736 }
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::Z2i::Integer
DGtal::int32_t Integer
DGtal::uint32_t
boost::uint32_t uint32_t
DGtal::PointVector::Coordinate
Component Coordinate
DGtal::DigitalSurface::ConstIterator
DigitalSurfaceContainer::SurfelConstIterator ConstIterator
DGtal::Dimension
DGtal::uint32_t Dimension
std
STL namespace.
DGtal::Trace::endBlock
double endBlock()
DGtal::DigitalSurface::DigitalSurfaceTracker
DigitalSurfaceContainer::DigitalSurfaceTracker DigitalSurfaceTracker
DGtal::DigitalSurface::DigitalSurfaceContainer
TDigitalSurfaceContainer DigitalSurfaceContainer
DGtal::KhalimskySpaceND::init
bool init(const Point &lower, const Point &upper, bool isClosed)
DGtal::KhalimskySpaceND::uIncident
Cell uIncident(const Cell &c, Dimension k, bool up) const
DGtal::DigitalSurface::KSpace
DigitalSurfaceContainer::KSpace KSpace
DGtal::trace
Trace trace(traceWriterTerm)
DGtal::Trace::info
std::ostream & info()
DGtal::KhalimskySpaceND::sKCoord
Integer sKCoord(const SCell &c, Dimension k) const
boost::vertices
std::pair< typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_iterator, typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_iterator > vertices(const DGtal::DigitalSurface< TDigitalSurfaceContainer > &digSurf)
DGtal::DigitalSurface::container
const DigitalSurfaceContainer & container() const
DGtal::KhalimskySpaceND::unsigns
Cell unsigns(const SCell &p) const
DGtal::functions::abs
T abs(const T &a)
DGtal::Z2i::RealVector
Space::RealVector RealVector
DGtal::Trace::error
std::ostream & error()
DGtal::KhalimskySpaceND::dimension
static const constexpr Dimension dimension
DGtal::PointVector::dot
Component dot(const Self &v) const
DGtal::DigitalSurface::begin
ConstIterator begin() const
DGtal::DigitalSurface::Surfel
DigitalSurfaceContainer::Surfel Surfel
DGtal::DigitalSurface::end
ConstIterator end() const
Companion project associated to DGtal ; Generated on Mon Mar 26 2018 12:12:21 for DGtalTools by   doxygen 1.8.10