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

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

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

/**
 * Example of a test. To be completed.
 *
 */
bool testLocalConvolutionNormalVectorEstimator ( int argc, char**argv )
{
    unsigned int nbok = 0;
    unsigned int nb = 0;

    trace.beginBlock ( "Testing convolution neighborhood ..." );

    QApplication application ( argc,argv );
    DGtal::Viewer3D<> viewer;

    std::string filename = testPath + "samples/cat10.vol";

    typedef ImageSelector < Z3i::Domain, int>::Type Image;
    Image image = VolReader<Image>::importVol ( filename );
    trace.info() <<image<<std::endl;
    DigitalSet set3d ( image.domain() );
    SetFromImage<DigitalSet>::append<Image> ( set3d, image,
            0,256 );

    KSpace ks;
    bool space_ok = ks.init ( image.domain().lowerBound(),
                              image.domain().upperBound(), true );
    if ( !space_ok )
    {
        trace.error() << "Error in the Khamisky space construction."<<std::endl;
        return true; //2; (return a bool !!!)
    }
    trace.endBlock();
    typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
    MySurfelAdjacency surfAdj ( true ); // interior in all directions.

    trace.beginBlock ( "Set up digital surface." );
    typedef LightImplicitDigitalSurface<KSpace, DigitalSet >
      MyDigitalSurfaceContainer;
    typedef DigitalSurface<MyDigitalSurfaceContainer> MyDigitalSurface;
    SCell bel = Surfaces<KSpace>::findABel ( ks, set3d, 100000 );
    MyDigitalSurfaceContainer* ptrSurfContainer =
        new MyDigitalSurfaceContainer ( ks, set3d, surfAdj, bel );
    MyDigitalSurface digSurf ( ptrSurfContainer ); // acquired

    MyDigitalSurface::ConstIterator it = digSurf.begin();


    //Convolution kernel
    deprecated::ConstantConvolutionWeights< MyDigitalSurface::Size > kernel;

    //Estimator definition
    typedef deprecated::LocalConvolutionNormalVectorEstimator<MyDigitalSurface,
                                                  deprecated::ConstantConvolutionWeights< MyDigitalSurface::Size > > MyEstimator;
    MyEstimator myNormalEstimator ( digSurf, kernel );

    myNormalEstimator.init ( 1.0, 5 );

    MyEstimator::Quantity res = myNormalEstimator.eval ( it );
    trace.info() << "Normal vector at begin() : "<< res << std::endl;

    viewer.show();

    DGtal::Color lineColorSave = viewer.getLineColor();
    viewer.setLineColor( DGtal::Color ( 200,20,20 ));
    for ( MyDigitalSurface::ConstIterator itbis = digSurf.begin(),itend=digSurf.end();
            itbis!=itend; ++itbis )
    {
        viewer << ks.unsigns ( *itbis );

        Point center = ks.sCoords ( *itbis );
        MyEstimator::Quantity normal = myNormalEstimator.eval ( itbis );

        viewer.addLine ( center,
                         DGtal::Z3i::RealPoint(center[0]-3*normal[0],
					       center[1]-3*normal[1],
					       center[2]-3*normal[2]) );
    }
    viewer.setLineColor( lineColorSave);
    viewer<< Viewer3D<>::updateDisplay;

    //Convolution kernel
    deprecated::GaussianConvolutionWeights< MyDigitalSurface::Size > Gkernel ( 14.0 );

    //Estimator definition
    typedef deprecated::LocalConvolutionNormalVectorEstimator<MyDigitalSurface,
                                                              deprecated::GaussianConvolutionWeights< MyDigitalSurface::Size > > MyEstimatorGaussian;
    MyEstimatorGaussian myNormalEstimatorG ( digSurf, Gkernel );

    myNormalEstimatorG.init ( 1.0, 15 );

    MyEstimatorGaussian::Quantity res2 = myNormalEstimatorG.eval ( it );
    trace.info() << "Normal vector at begin() : "<< res2 << std::endl;

    viewer<< CustomColors3D ( Color ( 200, 0, 0 ),Color ( 200, 0,0 ) );
    lineColorSave = viewer.getLineColor();
    viewer.setLineColor( DGtal::Color ( 200,20,20 ));
    for ( MyDigitalSurface::ConstIterator itbis = digSurf.begin(),itend=digSurf.end();
            itbis!=itend; ++itbis )
    {
        viewer << ks.unsigns ( *itbis );
//.........这里部分代码省略.........