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

void computeEstimation
( const po::variables_map& vm,     //< command-line parameters
  const KSpace& K,                 //< cellular grid space
  const ImplicitShape& shape,      //< implicit shape "ground truth"
  const Surface& surface,          //< digital surface approximating shape
  Estimator& estimator )           //< an initialized estimator
{
  typedef typename Surface::ConstIterator ConstIterator;
  typedef typename Surface::Surfel Surfel;
  typedef typename Estimator::Quantity Quantity;
  typedef double Scalar;
  typedef DepthFirstVisitor< Surface > Visitor;
  typedef GraphVisitorRange< Visitor > VisitorRange;
  typedef typename VisitorRange::ConstIterator VisitorConstIterator;
  
  std::string fname = vm[ "output" ].as<std::string>();
  string nameEstimator = vm[ "estimator" ].as<string>();

  trace.beginBlock( "Computing " + nameEstimator + " estimations." );
  CountedPtr<VisitorRange> range( new VisitorRange( new Visitor( surface, *(surface.begin()) )) );
  std::vector<Quantity> n_estimations;
  estimator.eval( range->begin(), range->end(), std::back_inserter( n_estimations ) );
  trace.info() << "- nb estimations  = " << n_estimations.size() << std::endl;
  trace.endBlock();

  trace.beginBlock( "Computing areas." );
  range = CountedPtr<VisitorRange>( new VisitorRange( new Visitor( surface, *(surface.begin()) )) );
  double area_est   = 0.0; // normal integration with absolute value.
  unsigned int i = 0;
  for ( typename VisitorRange::ConstIterator it = range->begin(), itE = range->end(); 
        it != itE; ++it, ++i )
    {
      Surfel s = *it;
      Dimension k = K.sOrthDir( s );
      area_est  += abs( n_estimations[ i ][ k ] );
    }  
  double h = vm["gridstep"].as<double>();
  trace.info() << setprecision(10) << "- Area_est     " << ( area_est * h * h )   << std::endl;
  std::ostringstream area_sstr;
  area_sstr << fname << "-" << nameEstimator << "-area-" << h << ".txt"; 
  std::ofstream area_output( area_sstr.str().c_str() );
  area_output << "# Area estimation by digital surface integration." << std::endl;
  area_output << "# X: " << nameEstimator << std::endl;
  area_output << "# h Area[X] nb_surf" << std::endl;
  area_output << setprecision(10) << h
              << " " << ( area_est * h * h )
              << " " << i << std::endl;
  area_output.close();
  trace.endBlock();
}

int main( int argc, char** argv )
{
    //! [greedy-plane-segmentation-ex3-parseCommandLine]
    trace.info() << "Segments the surface at given threshold within given volume into digital planes of rational width num/den." << std::endl;
    // Setting default options: ----------------------------------------------
    // input file used:
    string inputFilename =   examplesPath + "samples/Al.100.vol" ;
    trace.info() << "input file used " << inputFilename << std::endl;
    // parameter threshold
    unsigned int threshold = 0;
    trace.info() << "the value that defines the isosurface in the image (an integer between 0 and 255)= " << threshold<< std::endl;
    // parameter widthNum
    unsigned int widthNum = 1;
    trace.info() << "the numerator of the rational width (a non-null integer) =" << widthNum<< std::endl;
    // parameter widthDen
    unsigned int widthDen = 1;
    trace.info() << "the denominator of the rational width (a non-null integer)= " << widthDen<< std::endl;
    //! [greedy-plane-segmentation-ex3-parseCommandLine]

    //! [greedy-plane-segmentation-ex3-loadVolume]
    QApplication application(argc,argv);
    typedef ImageSelector < Domain, int>::Type Image;
    Image image = VolReader<Image>::importVol(inputFilename);
    DigitalSet set3d (image.domain());
    SetFromImage<DigitalSet>::append<Image>(set3d, image, threshold,255);
    //! [greedy-plane-segmentation-ex3-loadVolume]

    //! [greedy-plane-segmentation-ex3-makeSurface]
    trace.beginBlock( "Set up digital surface." );
    // We initializes the cellular grid space used for defining the
    // digital surface.
    KSpace ks;
    bool ok = ks.init( set3d.domain().lowerBound(),
                       set3d.domain().upperBound(), true );
    if ( ! ok ) std::cerr << "[KSpace.init] Failed." << std::endl;
    SurfelAdjacency<KSpace::dimension> surfAdj( true ); // interior in all directions.
    MyDigitalSurfaceContainer* ptrSurfContainer =
        new MyDigitalSurfaceContainer( ks, set3d, surfAdj );
    MyDigitalSurface digSurf( ptrSurfContainer ); // acquired
    trace.endBlock();
    //! [greedy-plane-segmentation-ex3-makeSurface]

    //! [greedy-plane-segmentation-ex3-segment]
    Point p;
    Dimension axis;
    unsigned int j = 0;
    unsigned int nb = digSurf.size();

    // First pass to find biggest planes.
    trace.beginBlock( "1) Segmentation first pass. Computes all planes so as to sort vertices by the plane size." );
    std::map<Vertex,unsigned int> v2size;
    NaivePlaneComputer planeComputer;
    std::priority_queue<VertexSize> Q;
    std::vector<Point> layer;
    for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
    {
        if ( ( (++j) % 50 == 0 ) || ( j == nb ) ) trace.progressBar( j, nb );
        Vertex v = *it;
        axis = ks.sOrthDir( v );
        planeComputer.init( axis, 500, widthNum, widthDen );
        // The visitor takes care of all the breadth-first traversal.
        Visitor visitor( digSurf, v );
        layer.clear();
        Visitor::Size currentSize = visitor.current().second;
        while ( ! visitor.finished() )
        {
            Visitor::Node node = visitor.current();
            v = node.first;
            axis = ks.sOrthDir( v );
            p = ks.sCoords( ks.sDirectIncident( v, axis ) );
            if ( node.second != currentSize )
            {
                // std::cerr << "Layer " << currentSize << ", size=" << layer.size() << std::endl;
                bool isExtended = planeComputer.extend( layer.begin(), layer.end() );
                if ( ! isExtended )
                    break;
                layer.clear();
                currentSize = node.second;
            }
            layer.push_back( p );
            visitor.expand();
        }
        // std::cerr << v << " -> " << planeComputer.size() << std::endl;
        Q.push( VertexSize( v, planeComputer.size() ) );
    }
    trace.endBlock();

    trace.beginBlock( "2) Segmentation second pass. Visits vertices from the one with biggest plane to the one with smallest plane." );
    std::set<Vertex> processedVertices;
    std::map<Vertex,SegmentedPlane*> v2plane;
    std::vector<SegmentedPlane*> segmentedPlanes;
    j = 0;
    while ( ! Q.empty() )
    {
        if ( ( (++j) % 50 == 0 ) || ( j == nb ) ) trace.progressBar( j, nb );
        Vertex v = Q.top().v;
        Q.pop();
        if ( processedVertices.find( v ) != processedVertices.end() ) // already in set
            continue; // process to next vertex

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

//.........这里部分代码省略.........
    trace.info() << "File extension not recognized: "<< extension << std::endl;
    return 0;
  }

  if(extension=="vol" || extension=="pgm3d" || extension=="pgm3D"
#ifdef WITH_ITK
     || extension =="dcm"
#endif
     ){
    trace.beginBlock( "Loading image into memory." );
#ifdef WITH_ITK
    int dicomMin = vm["dicomMin"].as<int>();
    int dicomMax = vm["dicomMax"].as<int>();
    typedef DGtal::functors::Rescaling<int ,unsigned char > RescalFCT;
    Image image = extension == "dcm" ? DicomReader< Image,  RescalFCT  >::importDicom( inputFilename,
										       RescalFCT(dicomMin,
												 dicomMax,
												 0, 255) ) :
      GenericReader<Image>::import( inputFilename );
#else
    Image image = GenericReader<Image>::import (inputFilename );
#endif
    trace.info() << "Image loaded: "<<image<< std::endl;
    trace.endBlock();

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

    trace.beginBlock( "Wrapping a digital set around image. " );
    typedef functors::IntervalForegroundPredicate<Image> ThresholdedImage;
    ThresholdedImage thresholdedImage( image, thresholdMin, thresholdMax );
    trace.endBlock();

    trace.beginBlock( "Extracting boundary by scanning the space. " );
    typedef KSpace::SurfelSet SurfelSet;
    typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
    typedef DigitalSurface< MySetOfSurfels > MyDigitalSurface;
    MySurfelAdjacency surfAdj( true ); // interior in all directions.
    MySetOfSurfels theSetOfSurfels( ks, surfAdj );
    Surfaces<KSpace>::sMakeBoundary( theSetOfSurfels.surfelSet(),
				     ks, thresholdedImage,
				     domain.lowerBound(),
				     domain.upperBound() );
    MyDigitalSurface digSurf( theSetOfSurfels );
    trace.info() << "Digital surface has " << digSurf.size() << " surfels."
		 << std::endl;
    trace.endBlock();

    trace.beginBlock( "Exporting everything." );
    Board3D<Space,KSpace> board(ks);

    board << SetMode3D(  ks.unsigns( *digSurf.begin() ).className(), "Basic" );

    typedef MyDigitalSurface::ConstIterator ConstIterator;
    if ( mode == "BDRY" )
      for ( ConstIterator it = digSurf.begin(), itE = digSurf.end(); it != itE; ++it )
        board << ks.unsigns( *it );
    else if ( mode == "INNER" )
      for ( ConstIterator it = digSurf.begin(), itE = digSurf.end(); it != itE; ++it )
        board << ks.sCoords( ks.sDirectIncident( *it, ks.sOrthDir( *it ) ) );
    else if ( mode == "OUTER" )
      for ( ConstIterator it = digSurf.begin(), itE = digSurf.end(); it != itE; ++it )
        board << ks.sCoords( ks.sIndirectIncident( *it, ks.sOrthDir( *it ) ) );
    else  if (mode == "CLOSURE")
    {
        std::set<KSpace::Cell> container;
        for ( ConstIterator it = digSurf.begin(), itE = digSurf.end(); it != itE; ++it )
        {
          container.insert( ks.unsigns( *it ) );
          KSpace::SCells oneNeig = ks.sLowerIncident(*it);
          //Processing linels
          for(KSpace::SCells::ConstIterator itt = oneNeig.begin(), ittend = oneNeig.end(); itt != ittend; ++itt)
          {
            container.insert( ks.unsigns( *itt) );
            KSpace::SCells oneNeig2 = ks.sLowerIncident(*itt);
            //Processing pointels
            for(KSpace::SCells::ConstIterator ittt = oneNeig2.begin(), itttend = oneNeig2.end(); ittt != itttend; ++ittt)
              container.insert( ks.unsigns(*ittt) );
          }
        }
      trace.info()<< "Exporting "<< container.size() << " cells"<<std::endl;
      for(auto cell: container)
        board << cell;
    }
    
    string outputFilename = vm["output"].as<std::string>();

    board.saveOBJ(outputFilename, normalization);
    trace.endBlock();
  }
  return 0;
}

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 << ") "
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  int thresholdMin = vm["thresholdMin"].as<int>();
  int thresholdMax = vm["thresholdMax"].as<int>();
  string mode = vm["mode"].as<string>();

  QApplication application(argc,argv);

  string extension = inputFilename.substr(inputFilename.find_last_of(".") + 1);
  if(extension!="vol" && extension != "p3d" && extension != "pgm3D" && extension != "pgm3d" && extension != "sdp" && extension != "pgm"
#ifdef WITH_ITK
     && extension !="dcm"
#endif
     ){
    trace.info() << "File extension not recognized: "<< extension << std::endl;
    return 0;
  }

  if(extension=="vol" || extension=="pgm3d" || extension=="pgm3D"
#ifdef WITH_ITK
     || extension =="dcm"
#endif
     ){
    trace.beginBlock( "Loading image into memory." );
#ifdef WITH_ITK
    int dicomMin = vm["dicomMin"].as<int>();
    int dicomMax = vm["dicomMax"].as<int>();
    typedef DGtal::functors::Rescaling<int ,unsigned char > RescalFCT;
    Image image = extension == "dcm" ? DicomReader< Image,  RescalFCT  >::importDicom( inputFilename,
                           RescalFCT(dicomMin,
                         dicomMax,
                         0, 255) ) :
      GenericReader<Image>::import( inputFilename );
#else
    Image image = GenericReader<Image>::import (inputFilename );
#endif
    trace.info() << "Image loaded: "<<image<< std::endl;
    trace.endBlock();

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

    //! [3dVolBoundaryViewer-Set3D]
    trace.beginBlock( "Wrapping a digital set around image. " );
    typedef functors::IntervalForegroundPredicate<Image> ThresholdedImage;
    ThresholdedImage thresholdedImage( image, thresholdMin, thresholdMax );
    trace.endBlock();
    //! [3dVolBoundaryViewer-Set3D]

    //! [3dVolBoundaryViewer-ExtractingSurface]
    trace.beginBlock( "Extracting boundary by scanning the space. " );
    typedef KSpace::SurfelSet SurfelSet;
    typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
    typedef DigitalSurface< MySetOfSurfels > MyDigitalSurface;
    MySurfelAdjacency surfAdj( true ); // interior in all directions.
    MySetOfSurfels theSetOfSurfels( ks, surfAdj );
    Surfaces<KSpace>::sMakeBoundary( theSetOfSurfels.surfelSet(),
             ks, thresholdedImage,
             domain.lowerBound(),
             domain.upperBound() );
    MyDigitalSurface digSurf( theSetOfSurfels );
    trace.info() << "Digital surface has " << digSurf.size() << " surfels."
     << std::endl;
    trace.endBlock();
    //! [3dVolBoundaryViewer-ExtractingSurface]

    //! [3dVolBoundaryViewer-ViewingSurface]
    trace.beginBlock( "Displaying everything. " );
    Viewer3D<Space,KSpace> viewer(ks);
    viewer.setWindowTitle("Simple boundary of volume Viewer");
    viewer.show();
    typedef MyDigitalSurface::ConstIterator ConstIterator;
    if ( mode == "BDRY" ){
      viewer << SetMode3D(ks.unsigns( *(digSurf.begin()) ).className(), "Basic");
      for ( ConstIterator it = digSurf.begin(), itE = digSurf.end(); it != itE; ++it )
  viewer << ks.unsigns( *it );
    }else if ( mode == "INNER" )
      for ( ConstIterator it = digSurf.begin(), itE = digSurf.end(); it != itE; ++it )
  viewer << ks.sCoords( ks.sDirectIncident( *it, ks.sOrthDir( *it ) ) );
    else if ( mode == "OUTER" )
      for ( ConstIterator it = digSurf.begin(), itE = digSurf.end(); it != itE; ++it )
  viewer << ks.sCoords( ks.sIndirectIncident( *it, ks.sOrthDir( *it ) ) );
    else{
      trace.error() << "Warning display mode (" << mode << ") not implemented." << std::endl;
      trace.error() << "The display will be empty." << std::endl;
    }
    viewer << Viewer3D<>::updateDisplay;
    trace.endBlock();
    return application.exec();
  }
  return 0;
}

int main( int argc, char** argv )
{
  QApplication application(argc,argv);
  string inputFilename = argc > 1 ? argv[ 1 ] : examplesPath+"/samples/Al.100.vol";
  int threshold = argc > 2 ? atoi( argv[ 2 ] ) : 0;
  int widthNum = argc > 3 ? atoi( argv[ 3 ] ) : 2;
  int widthDen = argc > 4 ? atoi( argv[ 4 ] ) : 1;

  //! [polyhedralizer-readVol]
  trace.beginBlock( "Reading vol file into an image." );
  typedef ImageContainerBySTLVector< Domain, int> Image;
  Image image = VolReader<Image>::importVol(inputFilename);
  typedef functors::SimpleThresholdForegroundPredicate<Image> DigitalObject;
  DigitalObject digitalObject( image, threshold );
  trace.endBlock();
  //! [polyhedralizer-readVol]

  //! [polyhedralizer-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."<<endl;
      return 2;
    }
  trace.endBlock();
  //! [polyhedralizer-KSpace]

  //! [polyhedralizer-SurfelAdjacency]
  typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
  MySurfelAdjacency surfAdj( false ); // exterior in all directions.
  //! [polyhedralizer-SurfelAdjacency]

  //! [polyhedralizer-ExtractingSurface]
  trace.beginBlock( "Extracting boundary by tracking the surface. " );
  typedef KSpace::Surfel Surfel;
  Surfel start_surfel = Surfaces<KSpace>::findABel( ks, digitalObject, 100000 );
  typedef ImplicitDigitalSurface< KSpace, DigitalObject > MyContainer;
  typedef DigitalSurface< MyContainer > MyDigitalSurface;
  typedef MyDigitalSurface::ConstIterator ConstIterator;
  MyContainer container( ks, digitalObject, surfAdj, start_surfel );
  MyDigitalSurface digSurf( container );
  trace.info() << "Digital surface has " << digSurf.size() << " surfels."
               << endl;
  trace.endBlock();
  //! [polyhedralizer-ExtractingSurface]

  //! [polyhedralizer-ComputingPlaneSize]
  // First pass to find biggest planes.
  trace.beginBlock( "Decomposition first pass. Computes all planes so as to sort vertices by the plane size." );
  typedef BreadthFirstVisitor<MyDigitalSurface> Visitor;
  typedef ChordGenericNaivePlaneComputer<Z3,Z3::Point, DGtal::int64_t> NaivePlaneComputer;
  map<Surfel,unsigned int> v2size;
  for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
    v2size[ *it ] = 0;
  int j = 0;
  int nb = digSurf.size();
  NaivePlaneComputer planeComputer;
  vector<Point> layer;
  vector<Surfel> layer_surfel;
  for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
    {
      if ( ( (++j) % 50 == 0 ) || ( j == nb ) ) trace.progressBar( j, nb );
      Surfel v = *it;
      planeComputer.init( widthNum, widthDen );
      // The visitor takes care of all the breadth-first traversal.
      Visitor visitor( digSurf, v );
      layer.clear();
      layer_surfel.clear();
      Visitor::Size currentSize = visitor.current().second;
      while ( ! visitor.finished() )
        {
          Visitor::Node node = visitor.current();
          v = node.first;
          int axis = ks.sOrthDir( v );
          Point p = ks.sCoords( ks.sDirectIncident( v, axis ) );
          if ( node.second != currentSize )
            {
              bool isExtended = planeComputer.extend( layer.begin(), layer.end() );
              if ( isExtended )
                {
                  for ( vector<Surfel>::const_iterator it_layer = layer_surfel.begin(),
                          it_layer_end = layer_surfel.end(); it_layer != it_layer_end; ++it_layer )
                    {
                      ++v2size[ *it_layer ];
                    }
                  layer_surfel.clear();
                  layer.clear();
                  currentSize = node.second;
                }
              else
                break;
            }
          layer_surfel.push_back( v );
          layer.push_back( p );
          visitor.expand();
        }
    }
  // Prepare queue
//.........这里部分代码省略.........