C++ Domain::begin方法代码示例

bool Apply_AntiBruit (MetaImage & image){
  Domain domain = image.domain();
  MetaImage new_image = MetaImage(domain);
  bool change=false;
  for (Domain::Iterator it = domain.begin(); it != domain.end();it++){
    int noirs =0; 		//variable indiquant le nombre de noirs
    int blancs =0;		//variable indiquant le nombre de blancs
    bool couleur=image(*it)>0;
    couleur?blancs=1:noirs=1;
    vector<Point> voisins = Anti_Bruit_Make_Voisins(*it);
    for (unsigned int i = 0; i < voisins.size(); i++){
        if (domain.isInside(voisins[i])){
	    image(voisins[i])>0?blancs++:noirs++;
	    //noirs  += 1-image(voisins[i]);
	    //blancs += image(voisins[i]);
	}
    }
    if((blancs>noirs)^couleur)change=true;
    if (noirs > blancs) {new_image.setValue(*it,0);}
    else {new_image.setValue(*it,1);}
//    new_image.setValue(*it,(1+(noirs>blancs)) %2);
  }
  image=new_image;
  return change;
}

void ConfigManager::writeDomain(FILE *file, const String &name, const Domain &domain) {
	if (domain.empty())
		return;		// Don't bother writing empty domains.

	String comment;

	// Write domain comment (if any)
	comment = domain.getDomainComment();
	if (!comment.empty())
		fprintf(file, "%s", comment.c_str());

	// Write domain start
	fprintf(file, "[%s]\n", name.c_str());

	// Write all key/value pairs in this domain, including comments
	Domain::const_iterator x;
	for (x = domain.begin(); x != domain.end(); ++x) {
		const String &value = x->_value;
		if (!value.empty()) {
			// Write comment (if any)
			if (domain.hasKVComment(x->_key)) {
				comment = domain.getKVComment(x->_key);
				fprintf(file, "%s", comment.c_str());
			}
			// Write the key/value pair
			fprintf(file, "%s=%s\n", x->_key.c_str(), value.c_str());
		}
	}
	fprintf(file, "\n");
}

bool MetaImage::removeNoise()
{
    Domain mydomain = this->domain();
    int l,h;
    l=mydomain.upperBound()[0]+1;
    h=mydomain.upperBound()[1]+1;
    vector<vector<bool> > newcol(l,vector<bool>(h,false));
  
  bool change=false;
  for (Domain::Iterator it = mydomain.begin(); it != mydomain.end();it++){
    int noirs =0; 		//variable indiquant le nombre de noirs
    int blancs =0;		//variable indiquant le nombre de blancs
    bool couleur=(*this)(*it)>0;
    couleur?blancs=1:noirs=1;
    vector<Point> voisins = TheVoisins(*it);
    for (unsigned int i = 0; i < voisins.size(); i++){
        if (mydomain.isInside(voisins[i])){
	    (*this)(voisins[i])>0?blancs++:noirs++;
	    //noirs  += 1-image(voisins[i]);
	    //blancs += image(voisins[i]);
	}
    }
    if((blancs>noirs)^couleur)change=true;
    if (noirs < blancs) {newcol[(*it)[0]][(*it)[1]]=1;}
//    new_image.setValue(*it,(1+(noirs>blancs)) %2);
  }
  
  for(int i=0;i<l;i++)
    for(int j=0;j<h;j++)
      setValue(Point(i,j),newcol[i][j]);
  return change;

}

void MetaImage::Open()
{
    Domain mydomain = this->domain();
    /*int l,h;
    l=domain.upperBound()[0]+1;
    h=domain.upperBound()[1]+1;*/
    queue<Point> colorize;
  
  for (Domain::Iterator it = mydomain.begin(); it != mydomain.end();it++){
    if((*this)(*it)>0)
    {
      vector<Point> voisins = TheVoisins(*it);
      for(Point p:voisins)
	if (mydomain.isInside(p))
	  colorize.push(*it);
    }
  }
  
  while(!colorize.empty())
  {
      setValue(colorize.front(),1);
      colorize.pop();
  }

}

// Iterate other the whole domain to get the first point
Point getInitPoint(DigitalSet set, Domain dom) {
	for(Domain::ConstIterator it = dom.begin(); it != dom.end(); ++it) {	// Iterate lines per lines from the bottom to the up, and (inside a line) from the left to the right
		Point p = *it;
		if (set(p))
			return p;
	}
}

void ConfigManager::writeDomain(WriteStream &stream, const String &name, const Domain &domain) {
	if (domain.empty())
		return;     // Don't bother writing empty domains.

	// WORKAROUND: Fix for bug #1972625 "ALL: On-the-fly targets are
	// written to the config file": Do not save domains that came from
	// the command line
	if (domain.contains("id_came_from_command_line"))
		return;

	String comment;

	// Write domain comment (if any)
	comment = domain.getDomainComment();
	if (!comment.empty())
		stream.writeString(comment);

	// Write domain start
	stream.writeByte('[');
	stream.writeString(name);
	stream.writeByte(']');
#ifdef _WIN32
	stream.writeByte('\r');
	stream.writeByte('\n');
#else
	stream.writeByte('\n');
#endif

	// Write all key/value pairs in this domain, including comments
	Domain::const_iterator x;
	for (x = domain.begin(); x != domain.end(); ++x) {
		if (!x->_value.empty()) {
			// Write comment (if any)
			if (domain.hasKVComment(x->_key)) {
				comment = domain.getKVComment(x->_key);
				stream.writeString(comment);
			}
			// Write the key/value pair
			stream.writeString(x->_key);
			stream.writeByte('=');
			stream.writeString(x->_value);
#ifdef _WIN32
			stream.writeByte('\r');
			stream.writeByte('\n');
#else
			stream.writeByte('\n');
#endif
		}
	}
#ifdef _WIN32
	stream.writeByte('\r');
	stream.writeByte('\n');
#else
	stream.writeByte('\n');
#endif
}

vector<double> distancetransform(Image image, Point barycentre){
	typedef functors::SimpleThresholdForegroundPredicate<Image> PointPredicate;
  PointPredicate predicate(image,0);
  //! [DTPredicate]  

  //! [DTCompute]
  typedef  DistanceTransformation<Z2i::Space, PointPredicate, Z2i::L2Metric> DTL2;
  typedef  DistanceTransformation<Z2i::Space, PointPredicate, Z2i::L1Metric> DTL1;
 
 
  DTL2 dtL2(image.domain(), predicate, Z2i::l2Metric);
  DTL1 dtL1(image.domain(), predicate, Z2i::l1Metric);
  //! [DTCompute]

  Domain d = dtL2.domain();
  DTL2::Value maxv2=0;
  vector<Point> sup;
  for(DTL2::Domain::ConstIterator it = d.begin(), itend = d.end(); it != itend; ++it) {
		DTL2::Value v = dtL2(*it);
		if(v>maxv2){
			sup.clear();
			sup.push_back(*it);
			maxv2 = v;
      //cout << "maxv2 " << maxv2 << " at " << *it << endl;
		}
		else if(v==maxv2){
			sup.push_back(*it);
      //cout << "more sup for " << maxv2 << endl;
		}
	}

  // distance moyenne entre le barycentre et les points atteignant maxv2 
  // ( a priori un seul point )
  double dist = 0;
  for(int i=0; i<sup.size(); i++){
    dist += (sup[i]-barycentre).norm();
  }
  dist = dist/ sup.size();

 
  DTL1::Value maxv1=0;
  //We compute the maximum DT value on the L1 map
  for ( DTL1::ConstRange::ConstIterator it = dtL1.constRange().begin(), itend = dtL1.constRange().end();it != itend; ++it)
    if ( (*it) > maxv1)  maxv1 = (*it);

  vector<double> result;
  result.push_back( (double) maxv2 );
  result.push_back( sup[0][0]);
  result.push_back( sup[0][1]);
  result.push_back( (double) maxv1 );
  result.push_back( (double) dtL2(barycentre) );
  result.push_back( dist );


  return result;
}

void
displayPredicate( Viewer3D & viewer,
                  const Domain & domain, const Predicate & pred )
{
  for ( typename Domain::ConstIterator itB = domain.begin(), itE = domain.end();
        itB != itE; ++itB )
    {
      if ( pred( *itB ) )
        viewer << *itB;
    }
}

std::vector<typename Domain::Point> pointsInStandardPlane
( const Domain & domain,
  typename Domain::Integer a, 
  typename Domain::Integer b,
  typename Domain::Integer c,
  typename Domain::Integer mu )
{
  typedef typename Domain::Integer Integer;
  typedef typename Domain::Point Point;
  typedef typename Domain::ConstIterator ConstIterator;
  std::vector<Point> pts;
  Integer mup = mu + abs(a) + abs(b) + abs(c);
  for ( ConstIterator it = domain.begin(), itE = domain.end();
        it != itE; ++it )
    {
      Point p = *it;
      Integer r = a * p[ 0 ] + b * p[ 1 ] + c * p[ 2 ];
      if ( ( mu <= r ) && ( r < mup ) )
        pts.push_back( p );
    }
  return pts;
}

//.........这里部分代码省略.........
  viewer.show();
  viewer.setGLScale(sx, sy, sz);

#ifdef WITH_ITK
  int dicomMin = vm["dicomMin"].as<int>();
  int dicomMax = vm["dicomMax"].as<int>();
  typedef DGtal::functors::Rescaling<int ,unsigned char > RescalFCT;

  Image3D image = extension == "dcm" ? DicomReader< Image3D,  RescalFCT  >::importDicom( inputFilename,
                                                                                         RescalFCT(dicomMin,
                                                                                                   dicomMax,
                                                                                                   0, 255) ) :
    GenericReader<Image3D>::import( inputFilename );
#else
  Image3D image = GenericReader<Image3D>::import( inputFilename );
#endif
  Domain domain = image.domain();

  trace.info() << "Image loaded: "<<image<< std::endl;
  viewer.setVolImage(&image);
  viewer << Z3i::Point(512, 512, 0);
  // Used to display 3D surface
  Z3i::DigitalSet set3d(domain);



  viewer << Viewer3D<>::updateDisplay;
  if(vm.count("thresholdImage")){
    GradientColorMap<long> gradient( thresholdMin, thresholdMax);
    gradient.addColor(Color::Blue);
    gradient.addColor(Color::Green);
    gradient.addColor(Color::Yellow);
    gradient.addColor(Color::Red);
    for(Domain::ConstIterator it = domain.begin(), itend=domain.end(); it!=itend; ++it){
      unsigned char  val= image( (*it) );
      Color c= gradient(val);
      if(val<=thresholdMax && val >=thresholdMin){
  if(!vm.count("displayDigitalSurface")){
          viewer <<  CustomColors3D(Color((float)(c.red()), (float)(c.green()),(float)(c.blue()), transp),
                                    Color((float)(c.red()), (float)(c.green()),(float)(c.blue()), transp));
          viewer << *it;
  }
      }else{
  set3d.insert(*it);
      }
    }
  }

  if(vm.count("displaySDP")){
    if(vm.count("colorSDP")){
      std::vector<int> vcol= vm["colorSDP"].as<std::vector<int > >();
      if(vcol.size()<4){
        trace.error() << "Not enough parameter: color specification should contains four elements: red, green, blue and alpha values."  << std::endl;
        return 0;
      }
      Color c(vcol[0], vcol[1], vcol[2], vcol[3]);
      viewer << CustomColors3D(c, c);
    }

    vector<Z3i::Point> vectVoxels;
    if(vm.count("SDPindex")) {
      std::vector<unsigned int > vectIndex = vm["SDPindex"].as<std::vector<unsigned int > >();
      if(vectIndex.size()!=3){
        trace.error() << "you need to specify the three indexes of vertex." << std::endl;
        return 0;
      }

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")
    ("input-file,i", po::value<std::string>(), "vol file (.vol) , pgm3d (.p3d or .pgm3d, pgm (with 3 dims)) file or sdp (sequence of discrete points)" )
    ("thresholdMin,m",  po::value<int>()->default_value(0), "threshold min to define binary shape" ) 
    ("thresholdMax,M",  po::value<int>()->default_value(255), "threshold max to define binary shape" )
    ("numMaxVoxel,n",  po::value<int>()->default_value(500000), "set the maximal voxel number to be displayed." )
#ifdef WITH_ITK
    ("dicomMin", po::value<int>()->default_value(-1000), "set minimum density threshold on Hounsfield scale")
    ("dicomMax", po::value<int>()->default_value(3000), "set maximum density threshold on Hounsfield scale")
#endif    
    ("transparency,t",  po::value<uint>()->default_value(255), "transparency") ; 

  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"
                << "Display volume file as a voxel set by using QGLviewer"<< endl
                << general_opt << "\n";
      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>();
  int thresholdMin = vm["thresholdMin"].as<int>();
  int thresholdMax = vm["thresholdMax"].as<int>();
  unsigned char transp = vm["transparency"].as<uint>();
  
  bool limitDisplay=false;
  if(vm.count("numMaxVoxel")){
    limitDisplay=true;
  }
  unsigned int numDisplayedMax = vm["numMaxVoxel"].as<int>();
  
  
  QApplication application(argc,argv);
  Viewer3D<> viewer;
  viewer.setWindowTitle("simple Volume Viewer");
  viewer.show();
 
  typedef ImageSelector<Domain, unsigned char>::Type Image;
  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
){
    unsigned int numDisplayed=0;
    
#ifdef WITH_ITK
   int dicomMin = vm["dicomMin"].as<int>();
   int dicomMax = vm["dicomMax"].as<int>();
   typedef DGtal::RescalingFunctor<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;
    Domain domain = image.domain();
    GradientColorMap<long> gradient( thresholdMin, thresholdMax);
    gradient.addColor(Color::Blue);
    gradient.addColor(Color::Green);
    gradient.addColor(Color::Yellow);
    gradient.addColor(Color::Red);
    for(Domain::ConstIterator it = domain.begin(), itend=domain.end(); it!=itend; ++it){
      unsigned char  val= image( (*it) );     
      if(limitDisplay && numDisplayed > numDisplayedMax)
	break;
      Color c= gradient(val);
      if(val<=thresholdMax && val >=thresholdMin){
	viewer <<  CustomColors3D(Color((float)(c.red()), (float)(c.green()),(float)(c.blue()), transp),
//.........这里部分代码省略.........

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")
    ("input-file,i", po::value<std::string>(), "volume file" )
    ("thresholdMin,m",  po::value<int>()->default_value(0), "threshold min to define binary shape" ) 
    ("thresholdMax,M",  po::value<int>()->default_value(255), "threshold max to define binary shape" )
    ("transparency,t",  po::value<uint>()->default_value(255), "transparency") ; 
  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"
    << "Display volume file as a voxel set by using QGLviewer"
    << general_opt << "\n";
      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>();
  int thresholdMin = vm["thresholdMin"].as<int>();
  int thresholdMax = vm["thresholdMax"].as<int>();
  unsigned char transp = vm["transparency"].as<uint>();
 
  QApplication application(argc,argv);
  Viewer3D viewer;
  viewer.setWindowTitle("simple Volume Viewer");
  viewer.show();
 
  typedef ImageSelector<Domain, unsigned char>::Type Image;
  Image image = VolReader<Image>::importVol( inputFilename );

  trace.info() << "Image loaded: "<<image<< std::endl;

  Domain domain = image.domain();
  GradientColorMap<long> gradient( thresholdMin, thresholdMax);
  gradient.addColor(Color::Blue);
  gradient.addColor(Color::Green);
  gradient.addColor(Color::Yellow);
  gradient.addColor(Color::Red);
  for(Domain::ConstIterator it = domain.begin(), itend=domain.end(); it!=itend; ++it){
    unsigned char  val= image( (*it) );     
   
    Color c= gradient(val);
    if(val<=thresholdMax && val >=thresholdMin){
      viewer <<  CustomColors3D(Color((float)(c.red()), (float)(c.green()),(float)(c.blue()), transp),
        Color((float)(c.red()), (float)(c.green()),(float)(c.blue()), transp));     
      viewer << *it;     
    }     
  }
  viewer << Viewer3D::updateDisplay;
  return application.exec();
}

int main( int argc, char** argv )
{
  QApplication app( argc, argv );
  trace.beginBlock ( "Testing class IVViewer" );
  IVViewer ivv( argc, argv );

  // Load volumetric image
  typedef ImageSelector<Domain, unsigned char>::Type Image;
  std::string filename = testPath + "samples/cat10.vol";
  Image image = VolReader<Image>::importVol( filename );
  trace.info() << image <<endl;

  // make shape.
  Domain domain =  image.domain() ;
  DigitalSet shape_set( domain );
  for ( Domain::ConstIterator it = domain.begin(), itend = domain.end();
  it != itend;   
  ++it )
    {
      if ( image( *it ) != 0 )
  shape_set.insert( *it );
    }

  // find first surfel on the boundary (not nice).
  Point first;
  for ( Domain::ConstIterator it = domain.begin(), itend = domain.end();
  it != itend;   
  ++it )
    {
      if ( image( *it ) != 0 )
  {
    first = *it;
    break;
  }
    }

  // Builds Khalimsky space.
  typedef KhalimskySpaceND<3> KSpace;
  typedef KSpace::SCell SCell;
  KSpace K3;
  SurfelAdjacency<KSpace::dimension> SAdj( true );
  K3.init( image.domain().lowerBound(), image.domain().upperBound(), true );

  // Tracks the shape boundary.
  SCell intvoxel = K3.sSpel( first );
  SCell surfel = K3.sIncident( intvoxel, 0, false );
  std::set<SCell> bdry;
  Surfaces<KSpace>::trackBoundary( bdry,
           K3, SAdj, shape_set, surfel );

  trace.info() << "tracking finished, size=" << bdry.size() << endl; 

  // Display surface.
  DGtalInventor<Space> inventor;
  typedef DGtalInventor<Space>::Color Color;
  addBounds( inventor, image.domain().lowerBound(), image.domain().upperBound() );
  for ( std::set<SCell>::const_iterator it = bdry.begin(), itend = bdry.end();
  it != itend;
  ++it )
    {
      inventor.setDiffuseColor( Color( 0.7, 0.7, 1.0 ) );
      inventor.drawCell( K3.sKCoords( *it ), 
       ! K3.sDirect( *it, K3.sOrthDir( *it ) ) );
    }
  // start surfel is red.
  inventor.setDiffuseColor( Color( 1.0, 0.0, 0.0 ) );
  inventor.drawCell( K3.sKCoords( surfel ), 
         ! K3.sDirect( surfel, K3.sOrthDir( surfel ) ) );
  inventor.generate( ivv.root() );
  // Qt will get the hand.
  ivv.show();

  bool res = true;
  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
  trace.endBlock();
  return res ? 0 : 1;
}

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

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")
    ("input,i", po::value<std::string>(), "vol file (.vol) , pgm3d (.p3d or .pgm3d, pgm (with 3 dims)) file or sdp (sequence of discrete points)" )
    ("output,o", po::value<std::string>(), "Output OBJ filename" )
    ("thresholdMin,m",  po::value<int>()->default_value(0), "threshold min to define binary shape" )
    ("thresholdMax,M",  po::value<int>()->default_value(255), "threshold max to define binary shape" );

  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"
                << "Convert a  volume file into OBJ format\n"
                << general_opt << "\n";
      return 0;
    }

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

  string inputFilename = vm["input"].as<std::string>();
  string outputFilename = vm["output"].as<std::string>();
  int thresholdMin = vm["thresholdMin"].as<int>();
  int thresholdMax = vm["thresholdMax"].as<int>();

  Board3D<> board;

  typedef ImageSelector<Domain, unsigned char>::Type Image;
  string extension = inputFilename.substr(inputFilename.find_last_of(".") + 1);
  if(extension!="vol" && extension != "p3d" && extension != "pgm3D" &&
     extension != "pgm3d" && extension != "sdp" && extension != "pgm")
    {
      trace.info() << "File extension not recognized: "<< extension << std::endl;
      return 0;
    }

  if(extension=="vol" || extension=="pgm3d" || extension=="pgm3D")
    {
      Image image = GenericReader<Image>::import (inputFilename );
      trace.info() << "Image loaded: "<<image<< std::endl;
      Domain domain = image.domain();
      for(Domain::ConstIterator it = domain.begin(), itend=domain.end(); it!=itend; ++it){
        unsigned char  val= image( (*it) );
        if(val<=thresholdMax && val >=thresholdMin){
          board << *it;
        }
      }
    }
  else
    if(extension=="sdp")
      {
        vector<Z3i::Point> vectVoxels = PointListReader<Z3i::Point>::getPointsFromFile(inputFilename);
        for(unsigned int i=0;i< vectVoxels.size(); i++){
          board << vectVoxels.at(i);
        }
      }


  board.saveOBJ(outputFilename);
  return 0;
}