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

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;
}

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

  trace.beginBlock ( "Testing direct Incidence in KSpace..." );
  for ( DirIterator q1 = K.sDirs( sspel ); q1 != 0; ++q1 )
    for ( DirIterator q2 = K.sDirs( sspel ); q2 != 0; ++q2 )
      {
  if ( *q1 != *q2 )
    {
      SCell s0 = K.sDirectIncident( sspel, *q1 );
      SCell l10 = K.sDirectIncident( s0, *q2 );
      SCell s1 = K.sDirectIncident( sspel, *q2 );
      SCell l01 = K.sDirectIncident( s1, *q1 );
      trace.info() << "Dd_" << *q2 << "(Dd_" << *q1 << "(V))=" << l10 
       << " Dd_" << *q1 << "(Dd_" << *q2 << "(V))=" << l01
       << endl;
      nbok += l10 != l01 ? 1 : 0; 
      nbok += K.sSign( s0 ) == K.POS ? 1 : 0;
      nbok += K.sSign( s1 ) == K.POS ? 1 : 0;
      nbok += K.sSign( l10 ) == K.POS ? 1 : 0;
      nbok += K.sSign( l01 ) == K.POS ? 1 : 0;
      nbok += s0 == K.sIncident( sspel, *q1, K.sDirect( sspel, *q1 ) )
        ? 1 : 0;
      nbok += s1 == K.sIncident( sspel, *q2, K.sDirect( sspel, *q2 ) )
        ? 1 : 0;
      nbok += l10 == K.sIncident( s0, *q2, K.sDirect( s0, *q2 ) )
        ? 1 : 0;
      nbok += l01 == K.sIncident( s1, *q1, K.sDirect( s1, *q1 ) )
        ? 1 : 0;
      nb += 9;
    }
      }
  trace.info() << "(" << nbok << "/" << nb << ") "
         << "correctness of direct and indirect orientations." << std::endl;
  
  trace.endBlock();
  
  
  return 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 << ") "
//.........这里部分代码省略.........

int main( int argc, char** argv )
{
  // parse command line ----------------------------------------------
  po::options_description general_opt("Allowed options are: ");
  general_opt.add_options()
    ("help,h", "display this message")
    ("scale,s", po::value<double>()->default_value(1.0), "set the scale of the maximal level. (default 1.0)")
    ("colorMap,c", "define the heightmap color with a pre-defined colormap (GradientColorMap)")
    ("colorTextureImage,t",po::value<std::string>(),  "define the heightmap color from a given color image (32 bits image).")
    ("input-file,i", po::value<std::string>(), "2d input image representing the height map (given as grayscape image cast into 8 bits)." );

  bool parseOK=true;
  po::variables_map vm;
  try{
    po::store(po::parse_command_line(argc, argv, general_opt), vm);
  }catch(const std::exception& ex){
    parseOK=false;
    trace.info()<< "Error checking program options: "<< ex.what()<< endl;
  }
  po::notify(vm);
  if( !parseOK || vm.count("help")||argc<=1)
    {
      std::cout << "Usage: " << argv[0] << " [input-file]\n"
                << "Displays 2D image as heightmap by using QGLviewer. "<< endl
                << general_opt << "\n" <<
                "Exemple of use:  visualisation/3dHeightMapViewer -i ${DGtal}/examples/samples/church.pgm -s 0.2" << std::endl;

      return 0;
    }

  if(! vm.count("input-file"))
    {
      trace.error() << " The file name was defined" << endl;
      return 0;
    }


  string inputFilename = vm["input-file"].as<std::string>();
  double scale = vm["scale"].as<double>();


  typedef DGtal::ImageContainerBySTLVector<Z2i::Domain, unsigned char> Image2DG ;
  typedef DGtal::ImageContainerBySTLVector<Z2i::Domain, unsigned int> Image2DCol ;

  Image2DG image = GenericReader<Image2DG>::import( inputFilename );
  Image2DCol imageTexture(image.domain());
  int maxHeight = (int)(std::numeric_limits<Image2DG::Value>::max()*scale);
  trace.info()<< "Max height from scale:" << maxHeight << std::endl;

  if(vm.count("colorTextureImage")){
    imageTexture =  GenericReader<Image2DCol>::import( vm["colorTextureImage"].as<std::string>() );
  }

  QApplication application(argc,argv);
  Z2i::RealPoint plow (image.domain().lowerBound()[0]-0.5,
                       image.domain().lowerBound()[1]-0.5);
  
  Z2i::RealPoint pup (image.domain().upperBound()[0]+0.5,
                      image.domain().upperBound()[1]+0.5);
  
  Viewer3DImageSpec<> viewer(plow, pup) ;
  viewer.setWindowTitle("Height Map Viewer");
  viewer.show();


  KSpace K;
  K.init(Z3i::Point(0,0,0),Z3i::Point(image.domain().upperBound()[0], image.domain().upperBound()[1], maxHeight+1), true);
  std::set<KSpace::SCell> boundVect;
  Image3DPredicatFrom2DImage<Image2DG, Z3i::Point> image3Dpredicate(image, scale);
  trace.info() << "Constructing boundary... ";
  Surfaces<KSpace>::sMakeBoundary (boundVect, K, image3Dpredicate, Z3i::Point(0,0,0),
                                   Z3i::Point(image.domain().upperBound()[0], image.domain().upperBound()[1], maxHeight+1)  );
  trace.info() << "[done]"<< std::endl;

  viewer << SetMode3D((*(boundVect.begin())).className(), "Basic" );
  GradientColorMap<Image2DG::Value,CMAP_JET>  gradientShade( 0, std::numeric_limits<Image2DG::Value>::max());
  GrayscaleColorMap<Image2DG::Value>  grayShade(0, std::numeric_limits<Image2DG::Value>::max());


  for(std::set<KSpace::SCell>::const_iterator it = boundVect.begin();
      it!= boundVect.end(); it++){
    Z3i::Point pt = K.sCoords(K.sDirectIncident( *it, 2 ));
    functors::Projector<SpaceND<2,int> > proj;
    Image2DG::Value val = image(proj(pt));
    if(vm.count("colorMap")){
      viewer.setFillColor(gradientShade(val));
    }else if (vm.count("colorTextureImage")) {
      viewer.setFillColor(Color(imageTexture(proj(pt))));
    }else{
      viewer.setFillColor(grayShade(val));
    }
  viewer << *it;
  }

  viewer << Viewer3D<>::updateDisplay;
  return application.exec();
}

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