C++ KSpace类代码示例

void
ballGenerator(const int& size, double aCx, double aCy, double aR, GridCurve<TKSpace>& gc)
{

  ASSERT( aR < (double) size ); 
 
  // Types
  typedef TKSpace KSpace;  
  typedef typename KSpace::SCell SCell;
  typedef typename KSpace::Space Space;  
  typedef typename Space::Point Point;

  KSpace K;
  bool ok = K.init( Point(-size,-size), Point(size,size), true );
  if ( ! ok )
    {
      std::cerr << " error in creating KSpace." << std::endl;
    }
  try 
    {
      BallPredicate<Point> dig(aCx, aCy, aR); 
      // Extracts shape boundary
      SurfelAdjacency<KSpace::dimension> SAdj( true );
      SCell bel = Surfaces<KSpace>::findABel( K, dig, 10000 );
      // Getting the consecutive surfels of the 2D boundary
      std::vector<Point> points;
      Surfaces<KSpace>::track2DBoundaryPoints( points, K, SAdj, dig, bel );
      gc.initFromVector(points); 
    }
  catch ( InputException e )
    {
      std::cerr << " error in finding a bel." << std::endl;
    }
}

int main( int argc, char** argv )
{
  if ( argc < 4 )
    {
      usage( argc, argv );
      return 1;
    }
  std::string inputFilename = argv[ 1 ];
  unsigned int minThreshold = atoi( argv[ 2 ] );
  unsigned int maxThreshold = atoi( argv[ 3 ] );

  //! [volScanBoundary-readVol]
  trace.beginBlock( "Reading vol file into an image." );
  typedef ImageSelector < Domain, int>::Type Image;
  Image image = VolReader<Image>::importVol(inputFilename);
  DigitalSet set3d (image.domain());
  SetFromImage<DigitalSet>::append<Image>(set3d, image, 
                                          minThreshold, maxThreshold);
  trace.endBlock();
  //! [volScanBoundary-readVol]
  
  
  //! [volScanBoundary-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."<<std::endl;
      return 2;
    }
  trace.endBlock();
  //! [volScanBoundary-KSpace]

  //! [volScanBoundary-ExtractingSurface]
  trace.beginBlock( "Extracting boundary by scanning the space. " );
  KSpace::SCellSet boundary;
  Surfaces<KSpace>::sMakeBoundary( boundary, ks, set3d,
                                   image.domain().lowerBound(), 
                                   image.domain().upperBound() );
  trace.endBlock();
  //! [volScanBoundary-ExtractingSurface]

  //! [volScanBoundary-DisplayingSurface]
  trace.beginBlock( "Displaying surface in Viewer3D." );
  QApplication application(argc,argv);
  Viewer3D viewer;
  viewer.show(); 
  viewer << CustomColors3D(Color(250, 0, 0 ), Color( 128, 128, 128 ) );
  unsigned long nbSurfels = 0;
  for ( KSpace::SCellSet::const_iterator it = boundary.begin(),
          it_end = boundary.end(); it != it_end; ++it, ++nbSurfels )
    viewer << *it;
  viewer << Viewer3D::updateDisplay;
  trace.info() << "nb surfels = " << nbSurfels << std::endl;
  trace.endBlock();
  return application.exec();
  //! [volScanBoundary-DisplayingSurface]
}

bool testFindABel()
{
  typedef typename KSpace::Point Point;
  typedef SpaceND< KSpace::dimension, typename KSpace::Integer > Space;
  typedef HyperRectDomain<Space> Domain;
  typedef typename DigitalSetSelector< Domain, BIG_DS+HIGH_BEL_DS >::Type DigitalSet;
  typedef typename KSpace::SCell SCell;
  typedef SetPredicate<DigitalSet> PointPredicate;

  trace.beginBlock("Test FindABel");
  Point low(-3,-3,-3), high(3,3,3);
  Domain domain( low, high );
  DigitalSet shape_set( domain );
  PointPredicate pp( shape_set );
  KSpace K;
  K.init( low, high, true );

  Point p000(0,0,0), p001(0,0,1), p010(0,1,0), p011(0,1,1),
        p100(1,0,0), p101(1,0,1), p110(1,1,0), p111(1,1,1);

  shape_set.insert( p000 );
  shape_set.insert( p100 );

  Surfaces<KSpace>::findABel( K, pp , p000 , p011 );
  Surfaces<KSpace>::findABel( K, pp , p000 , p110 );
  Surfaces<KSpace>::findABel( K, pp , p000 , p111 );
  Surfaces<KSpace>::findABel( K, pp , p000 , p101 );
  SCell s010 = Surfaces<KSpace>::findABel( K, pp , p000 , p010 );
  SCell s001 = Surfaces<KSpace>::findABel( K, pp , p000 , p001 );

  trace.endBlock();
  return ( (s010 == SCell( Point(1,2,1), true  ) ) &&
           (s001 == SCell( Point(1,1,2), false ) ) );
}

bool testNormaliation()
{
  trace.beginBlock ( "Testing normalization ..." );

  Point p1( 0, 0, 0 );
  Point p2( 0, 10 , 0);
  Point p3( 10, 10, 0);
  Point p4(10, 0, 100 );
  Point p5( 20, 0 , 0);
  Point p6( 20, 10, 0);
  
  KSpace k;

  k.init(Point(2,2,2), Point(4,4,4), true);

  Board3D<Space,KSpace> board(k);
  
  board << p1<<p2<<p3<<p4; 

  board.saveOBJ("dgtalBoard3D-norm.obj", true);
  board.saveOBJ("dgtalBoard3D-wonorm.obj");

  trace.endBlock();

  return true;
}

void displayProj2d( Viewer3D<space, kspace> & viewer,
		    const KSpace & ks, const StandardDSS6Computer & dss3d,
		    const DGtal::Color & color2d )
{
  typedef typename StandardDSS6Computer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
  typedef typename ArithmeticalDSSComputer2d::ConstIterator ConstIterator2d;
  typedef typename ArithmeticalDSSComputer2d::Point Point2d;
  typedef typename KSpace::Cell Cell;
  typedef typename KSpace::Point Point3d;
  Point3d b = ks.lowerBound();
  for ( DGtal::Dimension i = 0; i < 3; ++i )
    {
      const ArithmeticalDSSComputer2d & dss2d = dss3d.arithmeticalDSS2d( i );
      for ( ConstIterator2d itP = dss2d.begin(), itPEnd = dss2d.end(); itP != itPEnd; ++itP )
	{
	  Point2d p = *itP;
	  Point3d q;
	  switch (i) {
	  case 0: q = Point3d( 2*b[ i ]  , 2*p[ 0 ]+1, 2*p[ 1 ]+1 ); break;
	  case 1: q = Point3d( 2*p[ 0 ]+1, 2*b[ i ]  , 2*p[ 1 ]+1 ); break;
	  case 2: q = Point3d( 2*p[ 0 ]+1, 2*p[ 1 ]+1, 2*b[ i ]   ); break;
	  }
	  Cell c = ks.uCell( q );
	  viewer << CustomColors3D( color2d, color2d ) << c;
	}
    }
}

//-----------------------------------------------------------------------------
// Testing LightImplicitDigitalSurface
//-----------------------------------------------------------------------------
bool testLightImplicitDigitalSurface()
{
  using namespace Z3i;
  typedef ImplicitDigitalEllipse3<Point> ImplicitDigitalEllipse;
  typedef LightImplicitDigitalSurface<KSpace,ImplicitDigitalEllipse> Boundary;
  typedef Boundary::SurfelConstIterator ConstIterator;
  typedef Boundary::Tracker Tracker;
  typedef Boundary::Surfel Surfel;

  unsigned int nbok = 0;
  unsigned int nb = 0;
  trace.beginBlock ( "Testing block ... LightImplicitDigitalSurface" );
  Point p1( -10, -10, -10 );
  Point p2( 10, 10, 10 );
  KSpace K;
  nbok += K.init( p1, p2, true ) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
	       << "K.init() is ok" << std::endl;
  ImplicitDigitalEllipse ellipse( 6.0, 4.5, 3.4 );
  Surfel bel = Surfaces<KSpace>::findABel( K, ellipse, 10000 );
  Boundary boundary( K, ellipse, 
                     SurfelAdjacency<KSpace::dimension>( true ), bel );
  unsigned int nbsurfels = 0;
  for ( ConstIterator it = boundary.begin(), it_end = boundary.end();
        it != it_end; ++it )
    {
      ++nbsurfels;
    }
  trace.info() << nbsurfels << " surfels found." << std::endl;
  trace.beginBlock ( "Checks if adjacent surfels are part of the surface." );

  for ( ConstIterator it = boundary.begin(), it_end = boundary.end();
        it != it_end; ++it )
    {
      Tracker* ptrTracker = boundary.newTracker( *it );
      Surfel s = ptrTracker->current();
      Dimension trackDir = * K.sDirs( s );
      Surfel s1, s2;
      // unsigned int m1 = 
      ptrTracker->adjacent( s1, trackDir, true ); 
      // unsigned int m2 = 
      ptrTracker->adjacent( s2, trackDir, false ); 
      // trace.info() << "s = " << s << std::endl;
      // trace.info() << "s1 = " << s1 << " m1 = " << m1 << std::endl;
      // trace.info() << "s2 = " << s2 << " m2 = " << m2 << std::endl;
      nb++, nbok += boundary.isInside( s1 ) ? 1 : 0;
      // trace.info() << "(" << nbok << "/" << nb << ") "
      //              << "boundary.isInside( s1 )" << std::endl;
      nb++, nbok += boundary.isInside( s2 ) ? 1 : 0;
      // trace.info() << "(" << nbok << "/" << nb << ") "
      //              << "boundary.isInside( s2 )" << std::endl;
      delete ptrTracker;
    }
  trace.info() << "(" << nbok << "/" << nb << ") isInside tests." << std::endl;
  trace.endBlock();
  trace.endBlock();
  return nbok == nb;
}

/**
 * Example of a test. To be completed.
 *
 */
bool testDigitalSetBoundary()
{
  unsigned int nbok = 0;
  unsigned int nb = 0;
  
  trace.beginBlock ( "Testing block ... DigitalSetBoundary" );
  using namespace Z2i;
  typedef DigitalSetBoundary<KSpace,DigitalSet> Boundary;
  typedef Boundary::SurfelConstIterator ConstIterator;
  typedef Boundary::Tracker Tracker;
  typedef Boundary::Surfel Surfel;
  Point p1( -10, -10 );
  Point p2( 10, 10 );
  Domain domain( p1, p2 );
  DigitalSet dig_set( domain );
  Shapes<Domain>::addNorm2Ball( dig_set, Point( 0, 0 ), 5 );
  Shapes<Domain>::removeNorm2Ball( dig_set, Point( 0, 0 ), 1 );
  KSpace K;
  nbok += K.init( domain.lowerBound(), domain.upperBound(), true ) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
	       << "K.init() is ok" << std::endl;
  Boundary boundary( K, dig_set );
  unsigned int nbsurfels = 0;
  for ( ConstIterator it = boundary.begin(), it_end = boundary.end();
        it != it_end; ++it )
    {
      ++nbsurfels;
    }
  trace.info() << nbsurfels << " surfels found." << std::endl;
  nb++, nbok += nbsurfels == ( 12 + 44 ) ? 1 : 0;
  trace.info() << "(" << nbok << "/" << nb << ") "
	       << "nbsurfels == (12 + 44 )" << std::endl;
  for ( ConstIterator it = boundary.begin(), it_end = boundary.end();
        it != it_end; ++it )
    {
      Tracker* ptrTracker = boundary.newTracker( *it );
      Surfel s = ptrTracker->current();
      Dimension trackDir = * K.sDirs( s );
      Surfel s1, s2;
      unsigned int m1 = ptrTracker->adjacent( s1, trackDir, true ); 
      unsigned int m2 = ptrTracker->adjacent( s2, trackDir, false ); 
      trace.info() << "s = " << s << std::endl;
      trace.info() << "s1 = " << s1 << " m1 = " << m1 << std::endl;
      trace.info() << "s2 = " << s2 << " m2 = " << m2 << std::endl;
      nb++, nbok += boundary.isInside( s1 ) ? 1 : 0;
      trace.info() << "(" << nbok << "/" << nb << ") "
                   << "boundary.isInside( s1 )" << std::endl;
      nb++, nbok += boundary.isInside( s2 ) ? 1 : 0;
      trace.info() << "(" << nbok << "/" << nb << ") "
                   << "boundary.isInside( s2 )" << std::endl;
      delete ptrTracker;
    }
  trace.endBlock();
  return nbok == nb;
}

int main()
{
  KSpace K;
  Point plow(-3,-2);
  Point pup(5,3);
  Domain domain( plow, pup );
  Board2D board; // for 2D display
  K.init( plow, pup, true );
  board << SetMode( domain.styleName(), "Paving" )
	<< domain;
  Cell pixlow = K.uSpel( plow ); // pixel (-3*2+1,-2*2+1)
  Cell ptlow = K.uPointel( plow ); // pointel (-3*2,-2*2)
  Cell pixup = K.uSpel( pup ); // pixel (5*2+1,3*2+1)
  Cell ptup1 = K.uPointel( pup );   // pointel (5*2,3*2)
  Cell ptup2 = K.uTranslation( ptup1, Point::diagonal() ); // pointel (6*2,4*2)
  Cell linelb = K.uCell( Point( 1, 0 ) ); // linel (1,0) bottom
  Cell linelt = K.uCell( Point( 1, 2 ) ); // linel (1,2) top
  Cell linell = K.uCell( Point( 0, 1 ) ); // linel (0,1) left
  Cell linelr = K.uCell( Point( 2, 1 ) ); // linel (2,1) right
  board << CustomStyle( ptlow.styleName(), 
			new CustomColors( Color( 0, 0, 200 ), 
					  Color( 100, 100, 255 ) ) )
	<< ptlow << ptup2;
  board << CustomStyle( pixlow.styleName(), 
			new CustomColors( Color( 200, 0, 0 ), 
					  Color( 255, 100, 100 ) ) )
	<< pixlow << pixup;
  board << CustomStyle( linelb.styleName(), 
			new CustomColors( Color( 0, 200, 0 ), 
					  Color( 100, 255, 100 ) ) )
	<< linelb << linelt << linell << linelr;
  board.saveSVG("ctopo-1.svg");
  board.saveEPS("ctopo-1.eps");
  return 0;
}

bool testBallQuad()
{
  unsigned int nbok = 0;
  unsigned int nb = 0;

  trace.beginBlock ( "Testing... Ball with quadnormal");
  using namespace Z3i;
  typedef ImplicitDigitalBall3<Point> ImplicitDigitalBall;
  typedef ImplicitDigitalSurface<KSpace,ImplicitDigitalBall> Boundary;
  typedef Boundary::SurfelConstIterator ConstIterator;
  typedef Boundary::Surfel Surfel;
  Point p1( -50, -50, -50 );
  Point p2( 50, 50, 50 );
  KSpace K;
  nbok += K.init( p1, p2, true ) ? 1 : 0;
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << "K.init() is ok" << std::endl;
  ImplicitDigitalBall ball( 30.0 );
  Surfel bel = Surfaces<KSpace>::findABel( K, ball, 10000 );
  Boundary boundary( K, ball,
                     SurfelAdjacency<KSpace::dimension>( true ), bel );
  unsigned int nbsurfels = 0;

  Board3D<Space,KSpace> board(K);

  for ( ConstIterator it = boundary.begin(), it_end = boundary.end();
        it != it_end; ++it )
    {
      ++nbsurfels;

      Display3DFactory<>::drawOrientedSurfelWithNormal(board,
                                                       *it,
                                                       board.embedKS(*it).getNormalized());
    }

  trace.info() << nbsurfels << " surfels found." << std::endl;

  board.saveOBJ("testball.obj");


  nbok += true ? 1 : 0;
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << "true == true" << std::endl;
  trace.endBlock();

  return nbok == nb;
}

bool testBallQuad(int argc, char **argv)
{
  unsigned int nbok = 0;
  unsigned int nb = 0;

  QApplication application(argc, argv);

  trace.beginBlock ( "Testing... Ball with quadnormal");
  using namespace Z3i;
  typedef ImplicitDigitalBall3<Point> ImplicitDigitalBall;
  typedef ImplicitDigitalSurface<KSpace,ImplicitDigitalBall> Boundary;
  typedef Boundary::SurfelConstIterator ConstIterator;
  typedef Boundary::Surfel Surfel;
  Point p1( -100, -100, -100 );
  Point p2( 100, 100, 100 );
  KSpace K;
  nbok += K.init( p1, p2, true ) ? 1 : 0;
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << "K.init() is ok" << std::endl;
  ImplicitDigitalBall ball( 60.0 );
  Surfel bel = Surfaces<KSpace>::findABel( K, ball, 10000 );
  Boundary boundary( K, ball,
                     SurfelAdjacency<KSpace::dimension>( true ), bel );
  unsigned int nbsurfels = 0;

  Viewer3D<Space,KSpace> viewer(K);
  viewer.setWindowTitle("simpleViewer");
  viewer.show();


  for ( ConstIterator it = boundary.begin(), it_end = boundary.end();
        it != it_end; ++it )
    {
      ++nbsurfels;

      Display3DFactory<>::drawOrientedSurfelWithNormal(viewer,
                                                       *it,
                                                       viewer.embedKS(*it).getNormalized());
    }

  trace.info() << nbsurfels << " surfels found." << std::endl;
  viewer  << Display3D<Space, KSpace>::updateDisplay;

  bool res = application.exec();

  return res;
}

bool testQuadNorm()
{
  unsigned int nbok = 0;
  unsigned int nb = 0;
  trace.beginBlock ( "Testing Board3D Quads ..." );



  Point p1( 0, 0, 0 );
  Point p2( 0, 1 , 0);
  Point p3( 1, 1, 0);
  Point p4(1, 0, 0 );
  Point p5( 2, 0 , 0);
  Point p6( 2, 1, 0);
  RealVector n(1,1,1);
  RealVector n2(0,1,1);

  KSpace k;

  k.init(Point(2,2,2), Point(4,4,4), true);

  Board3D<Space,KSpace> board(k);
  board << CustomColors3D(Color(0, 255,0),Color(0, 255, 0));
  board.addQuadWithNormal(p1,p2,p3,p4, n.getNormalized(), true);
  board << CustomColors3D(Color(0, 0, 255),Color(0, 0, 255));
  board.addQuadWithNormal(p4,p5,p6,p3, n2.getNormalized(), true);

  Cell surfel = k.uCell( Point( 4,5,5) );
  Display3DFactory<Space,KSpace>::drawUnorientedSurfelWithNormal( board, surfel, n2.getNormalized());

  board.saveOBJ("dgtalBoard3D.quad.obj");


  nbok += true ? 1 : 0;
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
	       << "true == true" << std::endl;
  trace.endBlock();

  return nbok == nb;
}

bool testRaySurfelIntersection()
{
  unsigned int nbok = 0;
  unsigned int nb = 0;
  
  trace.beginBlock ( "Testing RaySurfel ..." );

  using namespace Z3i;
  
  KSpace k;

  k.init(Point(0,0,0), Point(10,10,10), true);
  
  typedef RayIntersectionPredicate<KSpace::Cell::Point> Ray;
  Ray  ray(KSpace::Cell::Point(0,0,0), 
           KSpace::Cell::Point(2,1,1));

  KSpace::Surfel surf =  k.sCell( Point( 2,1,1) );
  KSpace::Surfel surf2 =  k.sCell( Point( 2,7,7) );
  
  trace.info() << "Ray intersection with surf   "<<std::endl;
  nbok += ray(surf)  ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
               << "true " << std::endl;
  trace.info()<<std::endl;
  
  trace.info() << "Ray intersection with surf2  "<<std::endl;
  nbok += !ray(surf2 ) ? 1 : 0; 
  nb++;
  trace.info() << "(" << nbok << "/" << nb << ") "
	       << "false " << std::endl;
  trace.info()<<std::endl;

  trace.endBlock();
  
  return nbok == nb;
}

int main( int, char** )
{
  using namespace Z3i;
  typedef DGtal::ImplicitDigitalEllipse3<Point> ImplicitDigitalEllipse;
  typedef KSpace::SCell Surfel;
  bool res;
  trace.beginBlock ( "Testing class Object" );
  Point p1( -200, -200, -200 );
  Point p2( 200, 200, 200 );
  KSpace K;
  if ( K.init( p1, p2, true ) )
    {
      ImplicitDigitalEllipse ellipse( 180.0, 135.0, 102.0 );
      Surfel bel = Surfaces<KSpace>::findABel( K, ellipse, 10000 );
      res = testLightImplicitDigitalSurface<KSpace, ImplicitDigitalEllipse>
        ( K, ellipse, bel );
    }
  else
    res = false;
  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
  trace.endBlock();
  return res ? 0 : 1;
}

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

void displayDSS2d( Viewer3D<space, kspace> & viewer,
		   const KSpace & ks, const StandardDSS6Computer & dss3d,
		   const DGtal::Color & color2d )
{
  typedef typename StandardDSS6Computer::ConstIterator ConstIterator3d;
  typedef typename StandardDSS6Computer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
  typedef typename ArithmeticalDSSComputer2d::ConstIterator ConstIterator2d;
  typedef typename ArithmeticalDSSComputer2d::Point Point2d;
  typedef typename KSpace::Cell Cell;
  typedef typename KSpace::Point Point3d;
  typedef DGtal::PointVector<2,double> PointD2d;
  typedef typename Display3D<>::BallD3D PointD3D;
  Point3d b = ks.lowerBound();
  for ( DGtal::Dimension i = 0; i < 3; ++i )
    {
      const typename ArithmeticalDSSComputer2d::Primitive & dss2d 
	= dss3d.arithmeticalDSS2d( i ).primitive();
      // draw 2D bounding boxes for each arithmetical dss 2D.
      std::vector<PointD2d> pts2d;
      pts2d.push_back( dss2d.project(dss2d.back(), dss2d.Uf()) );
      pts2d.push_back( dss2d.project(dss2d.back(), dss2d.Lf()) );
      pts2d.push_back( dss2d.project(dss2d.front(), dss2d.Lf()) );
      pts2d.push_back( dss2d.project(dss2d.front(), dss2d.Uf()) );
      std::vector<PointD3D> bb;
      PointD3D p3;
      for ( unsigned int j = 0; j < pts2d.size(); ++j )
	{
	  switch (i) {
	  case 0: p3.center[0] = (double) b[ i ]-0.5; p3.center[1] = pts2d[ j ][ 0 ];  p3.center[2] = pts2d[ j ][ 1 ]; break;
	  case 1: p3.center[0] = pts2d[ j ][ 0 ];  p3.center[1] = (double) b[ i ]-0.5; p3.center[2] = pts2d[ j ][ 1 ];     break;
	  case 2: p3.center[0] = pts2d[ j ][ 0 ];  p3.center[1] = pts2d[ j ][ 1 ];     p3.center[2] = (double) b[ i ]-0.5; break;
	  }
	  bb.push_back( p3 );
	}
      for ( unsigned int j = 0; j < pts2d.size(); ++j ){
	viewer.setLineColor(color2d);
	viewer.addLine( DGtal::Z3i::RealPoint(bb[ j ].center[0], bb[ j ].center[1], bb[ j ].center[2]),
                        DGtal::Z3i::RealPoint(bb[ (j+1)%4 ].center[0], bb[ (j+1)%4 ].center[1], bb[ (j+1)%4 ].center[2]),
			MS3D_LINESIZE );
      }
    } // for ( DGtal::Dimension i = 0; i < 3; ++i )
}