C++ KSpace::sDirs方法代码示例

//-----------------------------------------------------------------------------
// Testing LightImplicitDigitalSurface
//-----------------------------------------------------------------------------
bool testLightImplicitDigitalSurface()
{
  using namespace Z3i;
  typedef ImplicitDigitalEllipse3<Point> ImplicitDigitalEllipse;
  typedef LightImplicitDigitalSurface<KSpace,ImplicitDigitalEllipse> Boundary;
  typedef Boundary::SurfelConstIterator ConstIterator;
  typedef Boundary::Tracker Tracker;
  typedef Boundary::Surfel Surfel;

  unsigned int nbok = 0;
  unsigned int nb = 0;
  trace.beginBlock ( "Testing block ... LightImplicitDigitalSurface" );
  Point p1( -10, -10, -10 );
  Point p2( 10, 10, 10 );
  KSpace K;
  nbok += K.init( p1, p2, true ) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
	       << "K.init() is ok" << std::endl;
  ImplicitDigitalEllipse ellipse( 6.0, 4.5, 3.4 );
  Surfel bel = Surfaces<KSpace>::findABel( K, ellipse, 10000 );
  Boundary boundary( K, ellipse, 
                     SurfelAdjacency<KSpace::dimension>( true ), bel );
  unsigned int nbsurfels = 0;
  for ( ConstIterator it = boundary.begin(), it_end = boundary.end();
        it != it_end; ++it )
    {
      ++nbsurfels;
    }
  trace.info() << nbsurfels << " surfels found." << std::endl;
  trace.beginBlock ( "Checks if adjacent surfels are part of the surface." );

  for ( ConstIterator it = boundary.begin(), it_end = boundary.end();
        it != it_end; ++it )
    {
      Tracker* ptrTracker = boundary.newTracker( *it );
      Surfel s = ptrTracker->current();
      Dimension trackDir = * K.sDirs( s );
      Surfel s1, s2;
      // unsigned int m1 = 
      ptrTracker->adjacent( s1, trackDir, true ); 
      // unsigned int m2 = 
      ptrTracker->adjacent( s2, trackDir, false ); 
      // trace.info() << "s = " << s << std::endl;
      // trace.info() << "s1 = " << s1 << " m1 = " << m1 << std::endl;
      // trace.info() << "s2 = " << s2 << " m2 = " << m2 << std::endl;
      nb++, nbok += boundary.isInside( s1 ) ? 1 : 0;
      // trace.info() << "(" << nbok << "/" << nb << ") "
      //              << "boundary.isInside( s1 )" << std::endl;
      nb++, nbok += boundary.isInside( s2 ) ? 1 : 0;
      // trace.info() << "(" << nbok << "/" << nb << ") "
      //              << "boundary.isInside( s2 )" << std::endl;
      delete ptrTracker;
    }
  trace.info() << "(" << nbok << "/" << nb << ") isInside tests." << std::endl;
  trace.endBlock();
  trace.endBlock();
  return nbok == nb;
}

/**
 * Example of a test. To be completed.
 *
 */
bool testDigitalSetBoundary()
{
  unsigned int nbok = 0;
  unsigned int nb = 0;
  
  trace.beginBlock ( "Testing block ... DigitalSetBoundary" );
  using namespace Z2i;
  typedef DigitalSetBoundary<KSpace,DigitalSet> Boundary;
  typedef Boundary::SurfelConstIterator ConstIterator;
  typedef Boundary::Tracker Tracker;
  typedef Boundary::Surfel Surfel;
  Point p1( -10, -10 );
  Point p2( 10, 10 );
  Domain domain( p1, p2 );
  DigitalSet dig_set( domain );
  Shapes<Domain>::addNorm2Ball( dig_set, Point( 0, 0 ), 5 );
  Shapes<Domain>::removeNorm2Ball( dig_set, Point( 0, 0 ), 1 );
  KSpace K;
  nbok += K.init( domain.lowerBound(), domain.upperBound(), true ) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
	       << "K.init() is ok" << std::endl;
  Boundary boundary( K, dig_set );
  unsigned int nbsurfels = 0;
  for ( ConstIterator it = boundary.begin(), it_end = boundary.end();
        it != it_end; ++it )
    {
      ++nbsurfels;
    }
  trace.info() << nbsurfels << " surfels found." << std::endl;
  nb++, nbok += nbsurfels == ( 12 + 44 ) ? 1 : 0;
  trace.info() << "(" << nbok << "/" << nb << ") "
	       << "nbsurfels == (12 + 44 )" << std::endl;
  for ( ConstIterator it = boundary.begin(), it_end = boundary.end();
        it != it_end; ++it )
    {
      Tracker* ptrTracker = boundary.newTracker( *it );
      Surfel s = ptrTracker->current();
      Dimension trackDir = * K.sDirs( s );
      Surfel s1, s2;
      unsigned int m1 = ptrTracker->adjacent( s1, trackDir, true ); 
      unsigned int m2 = ptrTracker->adjacent( s2, trackDir, false ); 
      trace.info() << "s = " << s << std::endl;
      trace.info() << "s1 = " << s1 << " m1 = " << m1 << std::endl;
      trace.info() << "s2 = " << s2 << " m2 = " << m2 << std::endl;
      nb++, nbok += boundary.isInside( s1 ) ? 1 : 0;
      trace.info() << "(" << nbok << "/" << nb << ") "
                   << "boundary.isInside( s1 )" << std::endl;
      nb++, nbok += boundary.isInside( s2 ) ? 1 : 0;
      trace.info() << "(" << nbok << "/" << nb << ") "
                   << "boundary.isInside( s2 )" << std::endl;
      delete ptrTracker;
    }
  trace.endBlock();
  return nbok == nb;
}

bool testCellularGridSpaceND()
{
  typedef typename KSpace::Cell Cell;
  typedef typename KSpace::SCell SCell;
  typedef typename KSpace::Point Point;
  typedef typename KSpace::DirIterator DirIterator;
  typedef typename KSpace::Cells Cells;
  typedef typename KSpace::SCells SCells;
  unsigned int nbok = 0;
  unsigned int nb = 0;
  
  trace.beginBlock ( "Testing block KSpace instantiation and scan ..." );
  KSpace K;
  int xlow[ 4 ] = { -3, -2, -2, -1 };
  int xhigh[ 4 ] = { 5, 3, 2, 3 };
  Point low( xlow );
  Point high( xhigh ); 
  bool space_ok = K.init( low, high, true );
  nbok += space_ok ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
         << "K.init( low, high )" << std::endl;
  trace.info() << "K.dim()=" << K.dimension << endl;
  int spel[ 4 ] = { 1, 1, 1, 1 }; // pixel
  Point kp( spel );
  Cell center = K.uCell( kp );
  Cell c1 = K.uCell( kp );
  Cell clow = K.uCell( low, kp );
  Cell chigh = K.uCell( high, kp );
  trace.info() << c1 << clow << chigh 
         << " topo(c1)=" << K.uTopology( c1 ) << " dirs=";
  for ( DirIterator q = K.uDirs( clow ); q != 0; ++q )
    trace.info() << " " << *q;
  trace.info() << endl;
  Cell f = K.uFirst( c1 );
  Cell l = K.uLast( c1 );
  trace.info() << "Loop in " << clow << chigh << endl;
  c1 = f;
  unsigned int nbelems = 0;
  do {
    ++nbelems;
    // trace.info() << c1;
  } while ( K.uNext( c1, f, l ) );
  trace.info() << " -> " << nbelems << " elements." << endl;
  unsigned int exp_nbelems = 1;
  for ( Dimension i = 0; i < K.dimension; ++i )
    exp_nbelems *= K.size( i );
  nbok += nbelems == exp_nbelems ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
         << nbelems << " scanned elements == "
         << exp_nbelems << " space size."
         << std::endl;
  trace.endBlock();
  trace.beginBlock ( "Testing neighborhoods in KSpace..." );
  Cells N = K.uNeighborhood( center );
  nbok += N.size() == ( K.dimension*2 + 1 ) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
         << N.size() << "(neighborhood size) == "
         << ( K.dimension*2 + 1 ) << "(2*dim()+1)" << endl;
  Cells Np = K.uProperNeighborhood( center );
  nbok += Np.size() == ( K.dimension*2 ) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
         << Np.size() << "(proper neighborhood size) == "
         << ( K.dimension*2 ) << "(2*dim())" << endl;
  trace.endBlock();

  trace.beginBlock ( "Testing faces in KSpace..." );
  Cells Nf = K.uFaces( center );
  nbok += Nf.size() == ceil( std::pow( 3.0 ,(int) K.dimension ) - 1 ) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
         << Nf.size() << "(faces size) == "
         << floor( std::pow( 3.0, (int)K.dimension ) - 1 ) << "(3^dim()-1)" << endl;
  trace.endBlock();
  
  trace.beginBlock ( "Testing block Incidence in KSpace..." );
  SCell sspel = K.sCell( kp, K.POS );
  for ( DirIterator q1 = K.sDirs( sspel ); q1 != 0; ++q1 )
    for ( DirIterator q2 = K.sDirs( sspel ); q2 != 0; ++q2 )
      {
  if ( *q1 != *q2 )
    {
      SCell s0 = K.sIncident( sspel, *q1, true );
      SCell s1 = K.sIncident( sspel, *q2, true );
      SCell l10 = K.sIncident( s0, *q2, true );
      SCell l01 = K.sIncident( s1, *q1, true );
      trace.info() << "D+_" << *q2 << "(D+_" << *q1 << "(V))=" << l10 
       << " D+_" << *q1 << "(D+_" << *q2 << "(V))=" << l01
       << endl;
      nbok += l10 == K.sOpp( l01 ) ? 1 : 0; 
      nb++;
    }
      }
  trace.info() << "(" << nbok << "/" << nb << ") "
         << "anti-commutativity of incidence operators." << std::endl;
  trace.endBlock();

//.........这里部分代码省略.........

int main( int argc, char** argv )
{
  // for 3D display with Viewer3D
  QApplication application(argc,argv);
  
  //! [digitalSurfaceSlice-readVol]
  trace.beginBlock( "Reading vol file into an image." );
  typedef ImageSelector < Domain, int>::Type Image;
  std::string inputFilename = examplesPath + "samples/Al.100.vol"; 
  Image image = VolReader<Image>::importVol(inputFilename);
  DigitalSet set3d (image.domain());
  SetPredicate<DigitalSet> set3dPredicate( set3d );
  SetFromImage<DigitalSet>::append<Image>(set3d, image, 
                                          0, 1 );
  Viewer3D viewer;  
  viewer.show(); 
  trace.endBlock();
  //! [digitalSurfaceSlice-readVol]
  
  //! [digitalSurfaceSlice-KSpace]
  trace.beginBlock( "Construct the Khalimsky space from the image domain." );
  KSpace ks;
  bool space_ok = ks.init( image.domain().lowerBound(), 
                           image.domain().upperBound(), true );
  if (!space_ok)
    {
      trace.error() << "Error in the Khamisky space construction."<<std::endl;
      return 2;
    }
  trace.endBlock();
  //! [digitalSurfaceSlice-KSpace]

  //! [digitalSurfaceSlice-SurfelAdjacency]
  typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
  MySurfelAdjacency surfAdj( true ); // interior in all directions.
  //! [digitalSurfaceSlice-SurfelAdjacency]

  //! [digitalSurfaceSlice-ExtractingSurface]
  trace.beginBlock( "Extracting boundary by scanning the space. " );
  typedef KSpace::Surfel Surfel;
  typedef KSpace::SurfelSet SurfelSet;
  typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
  typedef DigitalSurface< MySetOfSurfels > MyDigitalSurface;
  MySetOfSurfels theSetOfSurfels( ks, surfAdj );
  Surfaces<KSpace>::sMakeBoundary( theSetOfSurfels.surfelSet(),
                                   ks, set3dPredicate,
                                   image.domain().lowerBound(), 
                                   image.domain().upperBound() );
  MyDigitalSurface digSurf( theSetOfSurfels );
  trace.info() << "Digital surface has " << digSurf.size() << " surfels."
               << std::endl;
  trace.endBlock();
  //! [digitalSurfaceSlice-ExtractingSurface]
  
  //! [digitalSurfaceSlice-ExtractingSlice]
  trace.beginBlock( "Extract slices." );
  typedef MyDigitalSurface::DigitalSurfaceTracker MyTracker;
  typedef DigitalSurface2DSlice< MyTracker > My2DSlice;
  //Extract an initial boundary cell
  Surfel surf = *digSurf.begin();
  MyTracker* tracker1 = digSurf.container().newTracker( surf );
  MyTracker* tracker2 = digSurf.container().newTracker( surf );
  // Extract the bondary contour associated to the initial surfel in
  // its first direction
  My2DSlice slice1( tracker1, *(ks.sDirs( surf )) ); 
  // Extract the bondary contour associated to the initial surfel in
  // its second direction
  My2DSlice slice2( tracker2, *++(ks.sDirs( surf )) ); 
  delete tracker1;
  delete tracker2;
  trace.endBlock();
  //! [digitalSurfaceSlice-ExtractingSlice]

  ASSERT( slice1.start() == slice1.begin() );
  ASSERT( slice1.cstart() == slice1.c() );
  ASSERT( *slice1.begin() == surf );
  ASSERT( *slice1.c() == surf );
  ASSERT( *slice1.start() == surf );
  ASSERT( *slice1.cstart() == surf );
  ASSERT( *slice1.rcstart() == surf );
  ASSERT( slice1.rcstart() == slice1.rc() );
  ASSERT( *slice1.rc() == surf );
  ASSERT( *(slice1.c()+1) == *(slice1.begin()+1) );
  ASSERT( *(slice1.rc()+1) == *(slice1.rbegin()) );

  //! [digitalSurfaceSlice-displayingAll]
  trace.beginBlock( "Display all with QGLViewer." );
  // Displaying all the surfels in transparent mode
  viewer << SetMode3D( surf.className(), "Transparent");
  for( MyDigitalSurface::ConstIterator it = theSetOfSurfels.begin(),
         it_end = theSetOfSurfels.end(); it != it_end; ++it )
    viewer<< *it;
  
  // Displaying First surfels cut with gradient colors.;
  GradientColorMap<int> cmap_grad( 0, slice1.size() );
  cmap_grad.addColor( Color( 50, 50, 255 ) );
  cmap_grad.addColor( Color( 255, 0, 0 ) );
  cmap_grad.addColor( Color( 255, 255, 10 ) );
  
  // Need to avoid surfel superposition (the surfel size in increased)
//.........这里部分代码省略.........

int main( int argc, char** argv )
{
  trace.info() <<  "exampleGridCurve3d: the type can be changed in example source code with  <gridcurve>, <inner>, <outer>, <incident> " << std::endl;

  string type = "gridcurve";

  //vol reading and digital set construction
  trace.beginBlock( "Reading vol file into an image." );
  typedef ImageSelector < Domain, int>::Type Image;
  std::string inputFilename = examplesPath + "samples/cat10.vol";
  Image image = VolReader<Image>::importVol(inputFilename);
  DigitalSet set3d (image.domain());
  setFromImage( image, DigitalSetInserter<DigitalSet>(set3d), 1);
  trace.info() << set3d.size() << " voxels." << std::endl;
  trace.endBlock();

  //Khalimsky space construction
  trace.beginBlock( "Construct the Khalimsky space from the image domain." );
  KSpace ks;
/*  bool space_ok = ks.init( image.domain().lowerBound(),
                           image.domain().upperBound(), true );
  if (!space_ok)
    {
      trace.error() << "Error in the Khamisky space construction."<<std::endl;
      return 2;
    }*/
  trace.endBlock();

  //digital surface construction
  typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
  MySurfelAdjacency surfAdj( true ); // interior in all directions.

  trace.beginBlock( "Extracting boundary by scanning the space. " );
  typedef KSpace::Surfel Surfel;
  typedef KSpace::SurfelSet SurfelSet;
  typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
  typedef DigitalSurface< MySetOfSurfels > MyDigitalSurface;
  MySetOfSurfels theSetOfSurfels( ks, surfAdj );
  Surfaces<KSpace>::sMakeBoundary( theSetOfSurfels.surfelSet(),
                                   ks, set3d,
                                   image.domain().lowerBound(),
                                   image.domain().upperBound() );
  MyDigitalSurface digSurf( theSetOfSurfels );
  trace.info() << digSurf.size() << " surfels." << std::endl;
  trace.endBlock();

  //slice retrieving
  trace.beginBlock( "Extracting slice and constructing a grid curve. " );
  typedef MyDigitalSurface::DigitalSurfaceTracker MyTracker;
  typedef DigitalSurface2DSlice< MyTracker > My2DSlice;

  //Extract an initial boundary cell
  Surfel surf = *digSurf.begin();
  MyTracker* tracker = digSurf.container().newTracker( surf );

  // Extract the bondary contour associated to the initial surfel in
  // its first direction
  My2DSlice slice( tracker, *(ks.sDirs( surf )) );
  delete tracker;

  //! [exampleGridCurve3d-Construction]
  GridCurve<KSpace> gc(ks);
  gc.initFromSCellsRange( slice.begin(), slice.end() );
  //! [exampleGridCurve3d-Construction]

  trace.endBlock();


  // for 3D display with Viewer3D
  QApplication application(argc,argv);
  trace.beginBlock( "Display all with QGLViewer." );
  Viewer3D<> viewer;
  viewer.show();
  // Displaying all the surfels in transparent mode
  viewer << SetMode3D( surf.className(), "Transparent");
  for( MyDigitalSurface::ConstIterator it = theSetOfSurfels.begin(),
         it_end = theSetOfSurfels.end(); it != it_end; ++it )
    viewer<< *it;


  // Displaying slice
  viewer << Viewer3D<>::shiftSurfelVisu;
  viewer << SetMode3D( surf.className(), "");
  viewer.setFillColor( Color( 50, 50, 255 ) );

  if (type == "gridcurve")
    {
      viewer  << gc;
    }
  else if (type == "inner")
    {
      viewer << gc.getInnerPointsRange();
    }
  else if (type == "outer")
    {
      viewer << gc.getOuterPointsRange();
    }
  else if (type == "incident")
    {
      viewer << gc.getIncidentPointsRange();
//.........这里部分代码省略.........

bool testUmbrellaComputer()
{
  using namespace Z3i;
  
  typedef Space::RealPoint RealPoint;
  typedef ImplicitBall<Space> EuclideanShape;
  typedef GaussDigitizer<Space,EuclideanShape> DigitalShape; 
  typedef GaussDigitizer<Space,EuclideanShape>::Domain Domain;
  typedef LightImplicitDigitalSurface<KSpace,DigitalShape> Boundary;
  typedef Boundary::SurfelConstIterator ConstIterator;
  //typedef Boundary::Tracker Tracker;
  typedef Boundary::Surfel Surfel;
  typedef Boundary::DigitalSurfaceTracker DigitalSurfaceTracker;
  typedef DigitalSurface<Boundary> MyDigitalSurface;
  typedef UmbrellaComputer<DigitalSurfaceTracker> MyUmbrellaComputer;

  unsigned int nbok = 0;
  unsigned int nb = 0;
  trace.beginBlock ( "Testing block ... UmbrellaComputer" );
  // Creating shape
  Point c( 0, 0, 0 );
  EuclideanShape ball( c, 2 ); // ball r=4
  DigitalShape shape;
  shape.attach( ball );
  shape.init( RealPoint( -10.0, -10.0, -10.0 ), 
	      RealPoint( 10.0, 10.0, 10.0 ), 1.0 );
  // Creating cellular grid space around.
  Domain domain = shape.getDomain();
  KSpace K;
  nbok += K.init( domain.lowerBound(), domain.upperBound(), true ) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
	       << "K.init() is ok" << std::endl;
  // Find start surfel on surface.
  Surfel bel = Surfaces<KSpace>::findABel( K, shape, 10000 );
  // Define surface container then surface itself.
  Boundary boundary( K, // cellular space
		     shape, // point predicate
                     SurfelAdjacency<KSpace::dimension>( true ), // adjacency
		     bel // starting surfel
		     );
  MyDigitalSurface digSurf( boundary ); // boundary is cloned

  // Get tracker on surface.
  DigitalSurfaceTracker* ptrTracker = boundary.newTracker( bel );
  MyUmbrellaComputer umbrella;
  KSpace::DirIterator dirIt = K.sDirs( bel );
  Dimension k = *dirIt;
  Dimension j = *(++dirIt);
  trace.beginBlock ( "Testing block ... forward umbrella" );
  umbrella.init( *ptrTracker, k, true, j );
  unsigned int nb_forward = 0;
  Surfel init_bel = bel;
  do {
    Point x = K.sKCoords( bel );
    trace.info() << x << std::endl;
    umbrella.next();
    ++nb_forward;
    bel = umbrella.surfel();
  } while ( bel != init_bel );
  trace.endBlock();
  trace.beginBlock ( "Testing block ... backward umbrella" );
  unsigned int nb_backward = 0;
  do {
    Point x = K.sKCoords( bel );
    trace.info() << x << std::endl;
    umbrella.previous();
    ++nb_backward;
    bel = umbrella.surfel();
  } while ( bel != init_bel );
  nb++, nbok += nb_forward == nb_backward ? 1 : 0;
  
  trace.info() << "(" << nbok << "/" << nb << ") "
               << " nb_forward(" << nb_forward
	       << ") == nb_backward(" << nb_backward << ")"
	       << std::endl;
  trace.endBlock();
  unsigned int nbsurfels = 0;
  for ( ConstIterator it = boundary.begin(), it_end = boundary.end();
        it != it_end; ++it )
    {
      ++nbsurfels;
    }
  trace.info() << nbsurfels << " surfels found." << std::endl;

  trace.endBlock();

  delete ptrTracker;
  return nbok == nb;
}

//.........这里部分代码省略.........
    {
      viewer.setFillColor( Color( 200, 200, 250 ) );
      Display3DFactory<Space,KSpace>::drawOrientedSurfelWithNormal( viewer, *it, n_estimations[ *it ], false );
    }
  viewer << Viewer3D<>::updateDisplay;
  trace.endBlock();

  //-----------------------------------------------------------------------------
  // Defining Discrete Calculus.
  typedef CubicalComplex< KSpace >                                 CComplex;
  typedef DiscreteExteriorCalculus<2,3, EigenLinearAlgebraBackend> Calculus;
  typedef Calculus::Index                                          Index;
  typedef Calculus::PrimalForm0                                    PrimalForm0;
  typedef Calculus::PrimalForm1                                    PrimalForm1;
  typedef Calculus::PrimalForm2                                    PrimalForm2;
  typedef Calculus::PrimalIdentity0                                PrimalIdentity0;
  typedef Calculus::PrimalIdentity1                                PrimalIdentity1;
  typedef Calculus::PrimalIdentity2                                PrimalIdentity2;
  trace.beginBlock( "Creating Discrete Exterior Calculus. " );
  Calculus calculus;
  calculus.initKSpace<Domain>( domain );
  const KSpace& Kc = calculus.myKSpace; // should not be used.
  // Use a cubical complex to find all lower incident cells easily.
  CComplex complex( K );
  for ( ConstIterator it = digSurf.begin(), itE = digSurf.end(); it != itE; ++it )
    complex.insertCell( 2, K.unsigns( *it ) );
  complex.close();
  for ( CComplex::CellMapIterator it = complex.begin( 0 ), itE = complex.end( 0 ); it != itE; ++it )
    calculus.insertSCell( K.signs( it->first, K.POS ) );
  
  for ( CComplex::CellMapIterator it = complex.begin( 1 ), itE = complex.end( 1 ); it != itE; ++it )
    {
      SCell     linel = K.signs( it->first, K.POS );
      Dimension k     = * K.sDirs( linel );
      bool      pos   = K.sDirect( linel, k );
      calculus.insertSCell( pos ? linel : K.sOpp( linel ) );
      // calculus.insertSCell( K.signs( it->first, K.POS ) );
    }

  // for ( CComplex::CellMapIterator it = complex.begin( 2 ), itE = complex.end( 2 ); it != itE; ++it )
  // calculus.insertSCell( K.signs( it->first, K.POS ) );
  for ( ConstIterator it = digSurf.begin(), itE = digSurf.end(); it != itE; ++it )
    {
      calculus.insertSCell( *it );
      // SCell     surfel = *it;
      // Dimension k      = K.sOrthDir( surfel );
      // bool      pos    = K.sDirect( surfel, k );
      // calculus.insertSCell( pos ? *it : K.sOpp( *it ) );
    }
  calculus.updateIndexes();
  trace.info() << calculus << endl;

  std::vector<PrimalForm2> g;
  g.reserve( 3 );
  g.push_back( PrimalForm2( calculus ) );
  g.push_back( PrimalForm2( calculus ) );
  g.push_back( PrimalForm2( calculus ) );
  Index nb2 = g[ 0 ].myContainer.rows();
  
  for ( Index index = 0; index < nb2; index++)
    {
      const Calculus::SCell& cell = g[ 0 ].getSCell( index );
      if ( theSetOfSurfels.isInside( cell ) ) 
        {
          const RealVector&      n    = n_estimations[ cell ];
          g[ 0 ].myContainer( index ) = n[ 0 ];