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