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

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( )
{

  KSpace K;
  Point plow(0,0,0);
  Point pup(1,1,0);
  Domain domain( plow, pup );
  K.init( plow, pup, true );

  Board3D<Space, KSpace> board(K);

  // drawing cells of dimension 3

  SCell v2 = K.sSpel( Point( 1, 0, 0 ), KSpace::POS ); // +v
  SCell v3 = K.sSpel( Point( 0, 1, 0 ), KSpace::POS ); // +v
  SCell v4 = K.sSpel( Point( 1, 1, 0 ), KSpace::NEG ); // +v
  //! [SetKSIllustrationMode3D]
  SCell v = K.sSpel( Point( 0, 0, 0 ), KSpace::POS ); // +v
  board << SetMode3D( v.className(), "Illustration" );
  //! [SetKSIllustrationMode3D]

  board << v << v2 << v3;
  board.saveOBJ("board3D-2-ks.obj");

  Board3D<Space, KSpace> board2(K);

  // Surfel of Voxel (0,0)
  //! [KSIllustrationModeTransformed]
  SCell sx = K.sIncident( v, 0, true ); // surfel further along x
  DGtal::TransformedPrism tsx (sx, v);
  //! [KSIllustrationModeTransformed]

  SCell sy = K.sIncident( v, 1, true ); // surfel further along y
  SCell sz = K.sIncident( v, 2, true ); // surfel further along z
  SCell sxn = K.sIncident( v, 0, false ); // surfel further along x
  SCell syn = K.sIncident( v, 1, false ); // surfel further along y
  SCell szn = K.sIncident( v, 2, false ); // surfel further along z

  // Resizing and shifting the surfel towords its associated voxel (v).

  DGtal::TransformedPrism tsy (sy, v);
  DGtal::TransformedPrism tsz (sz, v);
  DGtal::TransformedPrism tsxn (sxn, v);
  DGtal::TransformedPrism tsyn (syn, v);
  DGtal::TransformedPrism tszn (szn, v);

  board2 << tsx << tsy << tsz << tsxn << tsyn << tszn;


  // Surfel of Voxel (1,0)
  SCell sx2 = K.sIncident( v2, 0, true ); // surfel further along x
  SCell sy2 = K.sIncident( v2, 1, true ); // surfel further along y
  SCell sz2 = K.sIncident( v2, 2, true ); // surfel further along z
  SCell sxn2 = K.sIncident( v2, 0, false ); // surfel further along x
  SCell syn2 = K.sIncident( v2, 1, false ); // surfel further along y
  SCell szn2 = K.sIncident( v2, 2, false ); // surfel further along z

  // Resizing and shifting the surfel towords its associated voxel (v2).
  DGtal::TransformedPrism tsx2 (sx2, v2);
  DGtal::TransformedPrism tsy2 (sy2, v2);
  DGtal::TransformedPrism tsz2 (sz2, v2);
  DGtal::TransformedPrism tsxn2 (sxn2, v2);
  DGtal::TransformedPrism tsyn2 (syn2, v2);
  DGtal::TransformedPrism tszn2 (szn2, v2);

  board2 << tsx2 << tsy2 << tsz2 << tsxn2 << tsyn2 << tszn2;


  // Surfel of Voxel (0,1)
  SCell sx3 = K.sIncident( v3, 0, true ); // surfel further along x
  SCell sy3 = K.sIncident( v3, 1, true ); // surfel further along y
  SCell sz3 = K.sIncident( v3, 2, true ); // surfel further along z
  SCell sxn3 = K.sIncident( v3, 0, false ); // surfel further along x
  SCell syn3 = K.sIncident( v3, 1, false ); // surfel further along y
  SCell szn3 = K.sIncident( v3, 2, false ); // surfel further along z

  // Shifting the surfel to its associated voxel (v3).
  DGtal::TransformedPrism tsx3 (sx3, v3);
  DGtal::TransformedPrism tsy3 (sy3, v3);
  DGtal::TransformedPrism tsz3 (sz3, v3);
  DGtal::TransformedPrism tsxn3 (sxn3, v3);
  DGtal::TransformedPrism tsyn3 (syn3, v3);
  DGtal::TransformedPrism tszn3 (szn3, v3);


  board2 << tsx3 << tsy3 << tsz3 << tsxn3 << tsyn3 << tszn3;

  std::cout << "save obj" << std::endl;
  board2.saveOBJ("board3D-2bis-ks.obj");


}

bool testCellDrawOnBoard()
{
  typedef typename KSpace::Integer Integer;
  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;
  typedef SpaceND<2, Integer> Z2;
  typedef HyperRectDomain<Z2> Domain;
  unsigned int nbok = 0;
  unsigned int nb = 0;
  trace.beginBlock ( "Testing cell draw on digital board ..." );
  KSpace K;
  int xlow[ 4 ] = { -3, -3 };
  int xhigh[ 4 ] = { 5, 3 };
  Point low( xlow );
  Point high( xhigh ); 
  bool space_ok = K.init( low, high, true );
  Domain domain( low, high );
  Board2D board;
  board.setUnit( LibBoard::Board::UCentimeter );
  board << SetMode( domain.className(), "Paving" )
  << domain;
  int spel[ 2 ] = { 1, 1 }; // pixel 0,0
  Point kp( spel );
  Cell uspel = K.uCell( kp );
  board << uspel 
  << low << high
  << K.uIncident( uspel, 0, false )
  << K.uIncident( uspel, 1, false );
  int spel2[ 2 ] = { 5, 1 }; // pixel 2,0
  Point kp2( spel2 );
  SCell sspel2 = K.sCell( kp2, K.POS );
  board << CustomStyle( sspel2.className(), 
      new CustomPen( Color( 200, 0, 0 ), 
               Color( 255, 100, 100 ),
               2.0, 
               Board2D::Shape::SolidStyle ) )
  << sspel2 
      << K.sIncident( sspel2, 0, K.sDirect( sspel2, 0 ) )
  << K.sIncident( sspel2, 1, K.sDirect( sspel2, 0 ) );
  board.saveEPS( "cells-1.eps" );
  board.saveSVG( "cells-1.svg" );
  trace.endBlock();
  board.clear();
  board << domain;
  SCell slinel0 = K.sIncident( sspel2, 0, K.sDirect( sspel2, 0 ) );
  SCell spointel01 = K.sIncident( slinel0, 1, K.sDirect( slinel0, 1 ) );
  board << CustomStyle( sspel2.className(), 
      new CustomColors( Color( 200, 0, 0 ), 
            Color( 255, 100, 100 ) ) )
  << sspel2
  << CustomStyle( slinel0.className(), 
      new CustomColors( Color( 0, 200, 0 ), 
            Color( 100, 255, 100 ) ) )
  << slinel0
  << CustomStyle( spointel01.className(), 
      new CustomColors( Color( 0, 0, 200 ), 
            Color( 100, 100, 255 ) ) )
  << spointel01;
  board.saveEPS( "cells-3.eps" );
  board.saveSVG( "cells-3.svg" );
  
  return ((space_ok) && (nbok == nb));
}

bool testSurfelAdjacency()
{
  typedef typename KSpace::Integer Integer;
  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, -3, -3, -3 };
  int xhigh[ 4 ] = { 5, 3, 3, 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;
  trace.endBlock();
  trace.beginBlock ( "Testing surfel adjacency ..." );
  SurfelAdjacency<KSpace::dimension> SAdj( true );
  for ( Dimension i = 0; i < K.dimension; ++i )
    for ( Dimension j = 0; j < K.dimension; ++j )
      if ( i != j )
  trace.info() << "(" << i << "," << j << ")=" 
         << ( SAdj.getAdjacency( i, j ) ? "i2e" : "e2i" );
  trace.info() << endl;
  trace.endBlock();

  int spel[ 4 ] = { 1, 1, 1, 1 }; // pixel
  Point kp( spel );
  SCell sspel = K.sCell( kp, K.POS );
  trace.beginBlock ( "Testing surfel directness ..." );
  for ( Dimension k = 0; k < K.dimension; ++k )
    {
      SCell surfel = K.sIncident( sspel, k, true );
      SCell innerspel = K.sDirectIncident( surfel, K.sOrthDir( surfel ) );
      trace.info() << "spel=" << sspel << " surfel=" << surfel
       << " innerspel=" << innerspel << endl;
      nbok += sspel == innerspel ? 1 : 0; 
      nb++;
      trace.info() << "(" << nbok << "/" << nb << ") "
       << "spel == innerspel" << std::endl;
      surfel = K.sIncident( sspel, k, false );
      innerspel = K.sDirectIncident( surfel, K.sOrthDir( surfel ) );
      trace.info() << "spel=" << sspel << " surfel=" << surfel
       << " innerspel=" << innerspel << endl;
      nbok += sspel == innerspel ? 1 : 0; 
      nb++;
      trace.info() << "(" << nbok << "/" << nb << ") "
       << "spel == innerspel" << std::endl;
    }
  trace.endBlock();

  SurfelNeighborhood<KSpace> SN;
  trace.beginBlock ( "Testing surfel neighborhood ..." );
  SCell surfel = K.sIncident( sspel, 0, false );
  SN.init( &K, &SAdj, surfel );
  trace.info() << "surfel      =" << surfel << endl;
  trace.info() << "follower1(+)=" << SN.follower1( 1, true ) << endl;
  trace.info() << "follower2(+)=" << SN.follower2( 1, true ) << endl;
  trace.info() << "follower3(+)=" << SN.follower3( 1, true ) << endl;
  trace.info() << "follower1(-)=" << SN.follower1( 1, false ) << endl;
  trace.info() << "follower2(-)=" << SN.follower2( 1, false ) << endl;
  trace.info() << "follower3(-)=" << SN.follower3( 1, false ) << endl;
  trace.endBlock();

  trace.beginBlock ( "Testing surface tracking ..." );
  typedef SpaceND< KSpace::dimension, Integer > Space;
  typedef HyperRectDomain<Space> Domain;
  typedef typename DigitalSetSelector< Domain, BIG_DS+HIGH_BEL_DS >::Type DigitalSet;
  Domain domain( low, high );
  DigitalSet shape_set( domain );
  SetPredicate<DigitalSet> shape_set_predicate( shape_set );
  int center[ 4 ] = { 1, 0, 0, 0 }; // pixel
  Point pcenter( center );
  Shapes<Domain>::addNorm1Ball( shape_set, pcenter, 1 );
  trace.info() << "surfel      = " << surfel << endl;
  SCell other1, other2;
  SN.getAdjacentOnDigitalSet( other1, shape_set, 1, K.sDirect( surfel, 1 ) );
  SN.getAdjacentOnDigitalSet( other2, shape_set, 1, !K.sDirect( surfel, 1 ) );
  trace.info() << "directNext  = " << other1 << endl;
  trace.info() << "indirectNext= " << other2 << endl;
  std::set<SCell> bdry;

  // surfel = Surfaces<KSpace>::findABel( K, shape_set );

  Surfaces<KSpace>::trackBoundary( bdry,
           K, SAdj, shape_set_predicate, surfel );
  trace.info() << "tracking finished, size=" << bdry.size() 
         << ", should be " << 2*K.dimension*(2*K.dimension-1) << endl;
  nbok += bdry.size() == ( 2*K.dimension*(2*K.dimension-1) ) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
//.........这里部分代码省略.........

//-----------------------------------------------------------------------------
// Testing LightExplicitDigitalSurface
//-----------------------------------------------------------------------------
bool testLightExplicitDigitalSurface()
{
  using namespace Z3i;
  typedef ImageContainerBySTLVector<Domain,DGtal::uint8_t> Image;
  typedef FrontierPredicate<KSpace, Image> SurfelPredicate;
  typedef LightExplicitDigitalSurface<KSpace,SurfelPredicate> Frontier;
  typedef Frontier::SurfelConstIterator ConstIterator;
  typedef Frontier::SCell SCell;

  unsigned int nbok = 0;
  unsigned int nb = 0;
  trace.beginBlock ( "Testing block ... LightExplicitDigitalSurface" );
  Point p1( -5, -5, -5 );
  Point p2( 5, 5, 5 );
  KSpace K;
  nbok += K.init( p1, p2, true ) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
	       << "K.init() is ok" << std::endl;
  Image image( Domain(p1, p2) );
  fillImage3D( image, p1, p2, 0 );
  fillImage3D( image, Point(-2,-2,-2 ), Point( 2, 2, 2 ), 1 );
  fillImage3D( image, Point( 0, 0,-2 ), Point( 0, 0, 2 ), 2 );
  fillImage3D( image, Point(-1,-1, 2 ), Point( 1, 1, 2 ), 2 );
  {
    SCell vox2  = K.sSpel( Point( 0, 0, 2 ), K.POS );
    SCell bel20 = K.sIncident( vox2, 2, true );
    SurfelPredicate surfPredicate( K, image, 2, 0 );
    Frontier frontier20( K, surfPredicate,
                         SurfelAdjacency<KSpace::dimension>( true ), 
                         bel20 );
    unsigned int nbsurfels = 0;
    for ( ConstIterator it = frontier20.begin(), it_end = frontier20.end();
          it != it_end; ++it )
      {
        ++nbsurfels;
      }
    trace.info() << nbsurfels << " surfels found." << std::endl;
    nb++, nbok += nbsurfels == 9 ? 1 : 0;
    trace.info() << "(" << nbok << "/" << nb << ") "
                 << "frontier20: nbsurfels == 9" << std::endl;
  }
  {
    SCell vox1  = K.sSpel( Point( 2, 0, 0 ), K.POS );
    SCell bel10 = K.sIncident( vox1, 0, true );
    SurfelPredicate surfPredicate( K, image, 1, 0 );
    Frontier frontier10( K, surfPredicate,
                         SurfelAdjacency<KSpace::dimension>( true ), 
                         bel10 );
    unsigned int nbsurfels = 0;
    for ( ConstIterator it = frontier10.begin(), it_end = frontier10.end();
          it != it_end; ++it )
      {
        ++nbsurfels;
      }
    trace.info() << nbsurfels << " surfels found." << std::endl;
    nb++, nbok += nbsurfels == 140 ? 1 : 0; // 4*25(sides) + 16(top) + 24(bot)
    trace.info() << "(" << nbok << "/" << nb << ") "
                 << "frontier10: nbsurfels == 140" << std::endl;
  }
  {
    SCell vox1  = K.sSpel( Point( 1, 0, 0 ), K.POS );
    SCell bel12 = K.sIncident( vox1, 0, false );
    SurfelPredicate surfPredicate( K, image, 1, 2 );
    Frontier frontier12( K, surfPredicate,
                         SurfelAdjacency<KSpace::dimension>( true ), 
                         bel12 );
    unsigned int nbsurfels = 0;
    for ( ConstIterator it = frontier12.begin(), it_end = frontier12.end();
          it != it_end; ++it )
      {
        ++nbsurfels;
      }
    trace.info() << nbsurfels << " surfels found." << std::endl;
    nb++, nbok += nbsurfels == 36 ? 1 : 0; // 8+12(top) + 16(axis) 
    trace.info() << "(" << nbok << "/" << nb << ") "
                 << "frontier12: nbsurfels == 36" << std::endl;
  }
  {
    typedef BoundaryPredicate<KSpace, Image> SecondSurfelPredicate;
    typedef LightExplicitDigitalSurface<KSpace,SecondSurfelPredicate> Boundary;
    typedef Boundary::SurfelConstIterator LEConstIterator;
    //typedef Boundary::Tracker Tracker;
    //typedef Boundary::SCell SCell;
    //typedef Boundary::Surfel Surfel;
    SCell vox1  = K.sSpel( Point( 1, 0, 0 ), K.POS );
    SCell bel1x = K.sIncident( vox1, 0, false );
    SecondSurfelPredicate surfPredicate( K, image, 1 );
    Boundary boundary1x( K, surfPredicate,
                         SurfelAdjacency<KSpace::dimension>( true ), 
                         bel1x );
    unsigned int nbsurfels = 0;
    for ( LEConstIterator it = boundary1x.begin(), it_end = boundary1x.end();
          it != it_end; ++it )
      {
        ++nbsurfels;
      }
//.........这里部分代码省略.........

int main( int argc, char** argv )
{
  //! [frontierAndBoundary-LabelledImage]
  typedef Space::RealPoint RealPoint;
  typedef ImplicitBall<Space> EuclideanShape;
  typedef GaussDigitizer<Space,EuclideanShape> DigitalShape;
  typedef ImageContainerBySTLVector<Domain,DGtal::uint8_t> Image;
  Point c1( 2, 0, 0 );
  int radius1 = 6;
  EuclideanShape ball1( c1, radius1 ); // ball r=6
  DigitalShape shape1;
  shape1.attach( ball1 );
  shape1.init( RealPoint( -10.0, -10.0, -10.0 ),
               RealPoint( 10.0, 10.0, 10.0 ), 1.0 );
  Point c2( -2, 0, 0 );
  int radius2 = 5;
  EuclideanShape ball2( c2, radius2 ); // ball r=6
  DigitalShape shape2;
  shape2.attach( ball2 );
  shape2.init( RealPoint( -10.0, -10.0, -10.0 ),
               RealPoint( 10.0, 10.0, 10.0 ), 1.0 );
  Domain domain = shape1.getDomain();
  Image image( domain ); // p1, p2 );
  std::cerr << std::endl;
  for ( Domain::ConstIterator it = domain.begin(), it_end = domain.end();
        it != it_end; ++it )
    {
      DGtal::uint8_t label = shape1( *it ) ? 1 : 0;
      label += shape2( *it ) ? 2 : 0;
      image.setValue( *it, label );
      std::cerr << (int) image( *it );
    }
  std::cerr << std::endl;
  //! [frontierAndBoundary-LabelledImage]

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

  //! [frontierAndBoundary-SetUpDigitalSurface]
  trace.beginBlock( "Set up digital surface." );
  typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
  MySurfelAdjacency surfAdj( true ); // interior in all directions.
  typedef functors::FrontierPredicate<KSpace, Image> FSurfelPredicate;
  typedef ExplicitDigitalSurface<KSpace,FSurfelPredicate> FrontierContainer;
  typedef DigitalSurface<FrontierContainer> Frontier;
  typedef functors::BoundaryPredicate<KSpace, Image> BSurfelPredicate;
  typedef ExplicitDigitalSurface<KSpace,BSurfelPredicate> BoundaryContainer;
  typedef DigitalSurface<BoundaryContainer> Boundary;
  // frontier between label 1 and 0 (connected part containing bel10)
  SCell vox1  = K.sSpel( c1 + Point( radius1, 0, 0 ), K.POS );
  SCell bel10 = K.sIncident( vox1, 0, true );
  FSurfelPredicate surfPredicate10( K, image, 1, 0 );
  Frontier frontier10 =
    new FrontierContainer( K, surfPredicate10, surfAdj, bel10 );
  // frontier between label 2 and 0 (connected part containing bel20)
  SCell vox2  = K.sSpel( c2 - Point( radius2, 0, 0 ), K.POS );
  SCell bel20 = K.sIncident( vox2, 0, false );
  FSurfelPredicate surfPredicate20( K, image, 2, 0 );
  Frontier frontier20 =
    new FrontierContainer( K, surfPredicate20, surfAdj, bel20 );
  // boundary of label 3 (connected part containing bel32)
  SCell vox3  = K.sSpel( c1 - Point( radius1, 0, 0 ), K.POS );
  SCell bel32 = K.sIncident( vox3, 0, false );
  BSurfelPredicate surfPredicate3( K, image, 3 );
  Boundary boundary3 =
    new BoundaryContainer( K, surfPredicate3, surfAdj, bel32 );
  trace.endBlock();
  //! [frontierAndBoundary-SetUpDigitalSurface]

  //! [volBreadthFirstTraversal-DisplayingSurface]
  trace.beginBlock( "Displaying surface in Viewer3D." );
  QApplication application(argc,argv);
  Viewer3D<> viewer;
  viewer.show();
  viewer << SetMode3D( domain.className(), "BoundingBox" )
         << domain;
  Cell dummy;
  // Display frontier between 1 and 0.
  unsigned int nbSurfels10 = 0;
  viewer << CustomColors3D( Color::Red, Color::Red );
  for ( Frontier::ConstIterator
          it = frontier10.begin(), it_end = frontier10.end();
        it != it_end; ++it, ++nbSurfels10 )
    viewer << *it;
  // Display frontier between 2 and 0.
  unsigned int nbSurfels20 = 0;
  viewer << CustomColors3D( Color::Yellow, Color::Yellow );
  for ( Frontier::ConstIterator
          it = frontier20.begin(), it_end = frontier20.end();
        it != it_end; ++it, ++nbSurfels20 )
    viewer << *it;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
          const PrimalForm1 former_v = v;
          trace.info() << "Building matrix tu_tA_A_u + BB + Mw2" << endl;
          PrimalIdentity1 V_Id1 = BB;
          for ( unsigned int d = 0; d < 3; ++d )
            {
              const PrimalIdentity1 A_u = diag( calculus, dual_h1 * dual_D0 * primal_h2 * u[ d ] );
              V_Id1.myContainer += square( calculus, A_u ).myContainer;
            }
          trace.info() << "Prefactoring matrix tu_tA_A_u + BB + Mw2" << endl;
          solver_v.compute( V_Id1 );
          trace.info() << "Solving (tu_tA_A_u + BB + Mw2) v = 1/(4eps) * l" << endl;
          v = solver_v.solve( (1.0/eps) * l_sur_4 );
          trace.info() << "  => " << ( solver_v.isValid() ? "OK" : "ERROR" ) 
                       << " " << solver_v.myLinearAlgebraSolver.info() << endl;
          trace.info() << "-------------------------------------------------------------------------------" << endl;
          trace.endBlock();

          for ( Index index = 0; index < nb2; index++)
            {
              double n2 = 0.0;
              for ( unsigned int d = 0; d < 3; ++d )
                n2 += u[ d ].myContainer( index ) * u[ d ].myContainer( index );
              double norm = sqrt( n2 );
              for ( unsigned int d = 0; d < 3; ++d )
                u[ d ].myContainer( index ) /= norm;
            }

          trace.beginBlock("Checking v, computing norms");
          double m1 = 1.0;
          double m2 = 0.0;
          double ma = 0.0;
          for ( Index index = 0; index < nb1; index++)
            {
              double val = v.myContainer( index );
              m1 = std::min( m1, val );
              m2 = std::max( m2, val );
              ma += val;
            }
          trace.info() << "1-form v: min=" << m1 << " avg=" << ( ma / nb1 ) << " max=" << m2 << std::endl;
          for ( Index index = 0; index < nb1; index++)
            v.myContainer( index ) = std::min( std::max(v.myContainer( index ), 0.0) , 1.0 );
          double n_infty = 0.0;
          double n_2 = 0.0;
          double n_1 = 0.0;
          for ( Index index = 0; index < nb1; index++)
            {
              n_infty = std::max( n_infty, fabs( v.myContainer( index ) - former_v.myContainer( index ) ) );
              n_2    += ( v.myContainer( index ) - former_v.myContainer( index ) )
                * ( v.myContainer( index ) - former_v.myContainer( index ) );
              n_1    += fabs( v.myContainer( index ) - former_v.myContainer( index ) );
            }
          n_1 /= v.myContainer.rows();
          n_2 = sqrt( n_2 / v.myContainer.rows() );
          
          trace.info() << "Variation |v^k+1 - v^k|_oo = " << n_infty << endl;
          trace.info() << "Variation |v^k+1 - v^k|_2 = " << n_2 << endl;
          trace.info() << "Variation |v^k+1 - v^k|_1 = " << n_1 << endl;
          trace.endBlock();
          if ( n_infty < 1e-4 ) break;
        } // for ( int i = 0; i < n; ++i )
    }
  trace.endBlock();

  //-----------------------------------------------------------------------------
  // Displaying regularized normals
  trace.beginBlock( "Displaying regularized normals. " );
  Viewer3D<Space,KSpace> viewerR( K );
  viewerR.setWindowTitle("Regularized normals");
  viewerR.show();
  viewerR << SetMode3D(K.unsigns( *(digSurf.begin()) ).className(), "Basic");
  viewerR.setFillColor( Color( 200, 200, 250 ) );
  for ( Index index = 0; index < nb2; index++)
    {
      const SCell& cell    = u[ 0 ].getSCell( index );
      // const RealVector& n  = n_estimations[ cell ];
      RealVector nr        = RealVector( u[ 0 ].myContainer( index ), 
                                         u[ 1 ].myContainer( index ), 
                                         u[ 2 ].myContainer( index ) );
      nr /= nr.norm();
      if ( theSetOfSurfels.isInside( cell ) ) 
        Display3DFactory<Space,KSpace>::drawOrientedSurfelWithNormal( viewerR, cell, nr, false );
      else
        Display3DFactory<Space,KSpace>::drawOrientedSurfelWithNormal( viewerR, K.sOpp( cell ), nr, false );
    }
  viewerR.setLineColor( Color( 255, 0, 0 ) );
  for ( Index index = 0; index < nb1; index++)
    {
      const SCell& cell    = v.getSCell( index );
      Dimension    k       = * K.sDirs( cell ); 
      const SCell  p0      = K.sIncident( cell, k, true );
      const SCell  p1      = K.sIncident( cell, k, false );
      if ( v.myContainer( index ) >= 0.5 ) continue;
      viewerR.addLine( embedder.embed( p0 ), embedder.embed( p1 ), (0.5 - v.myContainer( index ))/ 5.0 );
    }
  viewerR << Viewer3D<>::updateDisplay;
  trace.endBlock();

  return application.exec();

}