DGtalTools  0.9.2
3dCurvatureViewer.cpp
1 
29 #include <iostream>
31 #include "DGtal/base/Common.h"
32 #include <cstring>
33 
34 #include <boost/program_options/options_description.hpp>
35 #include <boost/program_options/parsers.hpp>
36 #include <boost/program_options/variables_map.hpp>
37 
38 // Shape constructors
39 #include "DGtal/io/readers/GenericReader.h"
40 #include "DGtal/images/ImageSelector.h"
41 #include "DGtal/images/imagesSetsUtils/SetFromImage.h"
42 #include "DGtal/images/IntervalForegroundPredicate.h"
43 #include "DGtal/topology/SurfelAdjacency.h"
44 #include "DGtal/topology/helpers/Surfaces.h"
45 #include "DGtal/topology/LightImplicitDigitalSurface.h"
46 #include <DGtal/topology/SetOfSurfels.h>
47 
48 #include "DGtal/images/ImageHelper.h"
49 #include "DGtal/topology/DigitalSurface.h"
50 #include "DGtal/graph/DepthFirstVisitor.h"
51 #include "DGtal/graph/GraphVisitorRange.h"
52 
53 // Integral Invariant includes
54 #include "DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h"
55 #include "DGtal/geometry/surfaces/estimation/IntegralInvariantVolumeEstimator.h"
56 #include "DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h"
57 
58 // Drawing
59 #include "DGtal/io/boards/Board3D.h"
60 #include "DGtal/io/colormaps/GradientColorMap.h"
61 
62 #ifdef WITH_VISU3D_QGLVIEWER
63 #include "DGtal/io/viewers/Viewer3D.h"
64 #endif
65 
66 using namespace DGtal;
67 using namespace functors;
68 
160 const Color AXIS_COLOR_RED( 200, 20, 20, 255 );
161 const Color AXIS_COLOR_GREEN( 20, 200, 20, 255 );
162 const Color AXIS_COLOR_BLUE( 20, 20, 200, 255 );
163 const double AXIS_LINESIZE = 0.05;
164 
166 
172 void missingParam( std::string param )
173 {
174  trace.error() << " Parameter: " << param << " is required.";
175  trace.info() << std::endl;
176 }
177 
178 namespace po = boost::program_options;
179 
180 int main( int argc, char** argv )
181 {
182  // parse command line ----------------------------------------------
183  po::options_description general_opt("Allowed options are");
184  general_opt.add_options()
185  ("help,h", "display this message")
186  ("input,i", po::value< std::string >(), ".vol file")
187  ("radius,r", po::value< double >(), "Kernel radius for IntegralInvariant" )
188  ("threshold,t", po::value< unsigned int >()->default_value(8), "Min size of SCell boundary of an object" )
189  ("minImageThreshold,l", po::value< int >()->default_value(0), "set the minimal image threshold to define the image object (object defined by the voxel with intensity belonging to ]minImageThreshold, maxImageThreshold ] )." )
190  ("maxImageThreshold,u", po::value< int >()->default_value(255), "set the minimal image threshold to define the image object (object defined by the voxel with intensity belonging to ]minImageThreshold, maxImageThreshold] )." )
191  ("mode,m", po::value< std::string >()->default_value("mean"), "type of output : mean, gaussian, k1, k2, prindir1, prindir2 or normal (default mean)")
192  ("exportOBJ,o", po::value< std::string >(), "Export the scene to specified OBJ/MTL filename (extensions added)." )
193  ("exportDAT,d", po::value<std::string>(), "Export resulting curvature (for mean, gaussian, k1 or k2 mode) in a simple data file each line representing a surfel. ")
194  ("exportOnly", "Used to only export the result without the 3d Visualisation (usefull for scripts)." )
195  ("imageScale,s", po::value<std::vector<double> >()->multitoken(), "scaleX, scaleY, scaleZ: re sample the source image according with a grid of size 1.0/scale (usefull to compute curvature on image defined on anisotropic grid). Set by default to 1.0 for the three axis. ")
196  ("normalization,n", "When exporting to OBJ, performs a normalization so that the geometry fits in [-1/2,1/2]^3") ;
197 
198  bool parseOK = true;
199  po::variables_map vm;
200  try
201  {
202  po::store( po::parse_command_line( argc, argv, general_opt ), vm );
203  }
204  catch( const std::exception & ex )
205  {
206  parseOK = false;
207  trace.error() << " Error checking program options: " << ex.what() << std::endl;
208  }
209  bool neededArgsGiven=true;
210 
211  if (parseOK && !(vm.count("input"))){
212  missingParam("--input");
213  neededArgsGiven=false;
214  }
215  if (parseOK && !(vm.count("radius"))){
216  missingParam("--radius");
217  neededArgsGiven=false;
218  }
219 
220  bool normalization = false;
221  if (parseOK && vm.count("normalization"))
222  normalization = true;
223 
224  std::string mode;
225  if( parseOK )
226  mode = vm["mode"].as< std::string >();
227  if ( parseOK && ( mode.compare("gaussian") != 0 ) && ( mode.compare("mean") != 0 ) &&
228  ( mode.compare("k1") != 0 ) && ( mode.compare("k2") != 0 ) &&
229  ( mode.compare("prindir1") != 0 ) && ( mode.compare("prindir2") != 0 ) && ( mode.compare("normal") != 0 ))
230  {
231  parseOK = false;
232  trace.error() << " The selected mode ("<<mode << ") is not defined."<<std::endl;
233  }
234 
235 #ifndef WITH_VISU3D_QGLVIEWER
236  bool enable_visu = false;
237 #else
238  bool enable_visu = !vm.count("exportOnly");
239 #endif
240  bool enable_obj = vm.count("exportOBJ");
241  bool enable_dat = vm.count("exportDAT");
242 
243  if( !enable_visu && !enable_obj && !enable_dat )
244  {
245 #ifndef WITH_VISU3D_QGLVIEWER
246  trace.error() << "You should specify what you want to export with --export and/or --exportDat." << std::endl;
247 #else
248  trace.error() << "You should specify what you want to export with --export and/or --exportDat, or remove --exportOnly." << std::endl;
249 #endif
250  neededArgsGiven = false;
251  }
252 
253  if(!neededArgsGiven || !parseOK || vm.count("help") || argc <= 1 )
254  {
255  trace.info()<< "Visualisation of 3d curvature from .vol file using curvature from Integral Invariant" <<std::endl
256  << general_opt << "\n"
257  << "Basic usage: "<<std::endl
258  << "\t3dCurvatureViewer -i file.vol --radius 5 --mode mean"<<std::endl
259  << std::endl
260  << "Below are the different available modes: " << std::endl
261  << "\t - \"mean\" for the mean curvature" << std::endl
262  << "\t - \"gaussian\" for the Gaussian curvature" << std::endl
263  << "\t - \"k1\" for the first principal curvature" << std::endl
264  << "\t - \"k2\" for the second principal curvature" << std::endl
265  << "\t - \"prindir1\" for the first principal curvature direction" << std::endl
266  << "\t - \"prindir2\" for the second principal curvature direction" << std::endl
267  << "\t - \"normal\" for the normal vector" << std::endl
268  << std::endl;
269  return 0;
270  }
271  unsigned int threshold = vm["threshold"].as< unsigned int >();
272  int minImageThreshold = vm["minImageThreshold"].as< int >();
273  int maxImageThreshold = vm["maxImageThreshold"].as< int >();
274 
275  double h = 1.0;
276 
277  std::string export_obj_filename;
278  std::string export_dat_filename;
279 
280  if( enable_obj )
281  {
282  export_obj_filename = vm["exportOBJ"].as< std::string >();
283  if( export_obj_filename.find(".obj") == std::string::npos )
284  {
285  std::ostringstream oss;
286  oss << export_obj_filename << ".obj" << std::endl;
287  export_obj_filename = oss.str();
288  }
289  }
290 
291 
292  if( enable_dat )
293  {
294  export_dat_filename = vm["exportDAT"].as<std::string>();
295  }
296 
297  double re_convolution_kernel = vm["radius"].as< double >();
298 
299 
300  std::vector< double > aGridSizeReSample;
301  if( vm.count( "imageScale" ))
302  {
303  std::vector< double> vectScale = vm["imageScale"].as<std::vector<double > >();
304  if( vectScale.size() != 3 )
305  {
306  trace.error() << "The grid size should contains 3 elements" << std::endl;
307  return 0;
308  }
309  else
310  {
311  aGridSizeReSample.push_back(1.0/vectScale.at(0));
312  aGridSizeReSample.push_back(1.0/vectScale.at(1));
313  aGridSizeReSample.push_back(1.0/vectScale.at(2));
314  }
315  }
316  else
317  {
318  aGridSizeReSample.push_back(1.0);
319  aGridSizeReSample.push_back(1.0);
320  aGridSizeReSample.push_back(1.0);
321  }
322 
323 
324 
325  // Construction of the shape from vol file
326  typedef Z3i::Space::RealPoint RealPoint;
327  typedef Z3i::Point Point;
328  typedef ImageSelector< Z3i::Domain, int>::Type Image;
329  typedef DGtal::functors::BasicDomainSubSampler< HyperRectDomain<SpaceND<3, int> >,
330  DGtal::int32_t, double > ReSampler;
331  typedef DGtal::ConstImageAdapter<Image, Image::Domain, ReSampler,
332  Image::Value, DGtal::functors::Identity > SamplerImageAdapter;
333  typedef IntervalForegroundPredicate< SamplerImageAdapter > ImagePredicate;
334  typedef BinaryPointPredicate<DomainPredicate<Image::Domain>, ImagePredicate, AndBoolFct2 > Predicate;
335  typedef Z3i::KSpace KSpace;
336  typedef KSpace::SCell SCell;
337  typedef KSpace::Cell Cell;
338  typedef KSpace::Surfel Surfel;
339 
340  trace.beginBlock("Loading the file");
341  std::string filename = vm["input"].as< std::string >();
342  Image image = GenericReader<Image>::import( filename );
343 
344  PointVector<3,int> shiftVector3D( 0 ,0, 0 );
345  DGtal::functors::BasicDomainSubSampler< HyperRectDomain< SpaceND< 3, int > >,
346  DGtal::int32_t, double > reSampler(image.domain(),
347  aGridSizeReSample, shiftVector3D);
348  const functors::Identity identityFunctor{};
349  SamplerImageAdapter sampledImage ( image, reSampler.getSubSampledDomain(), reSampler, identityFunctor );
350  ImagePredicate predicateIMG = ImagePredicate( sampledImage, minImageThreshold, maxImageThreshold );
351  DomainPredicate<Z3i::Domain> domainPredicate( sampledImage.domain() );
352  AndBoolFct2 andF;
353  Predicate predicate(domainPredicate, predicateIMG, andF );
354 
355 
356  Z3i::Domain domain = sampledImage.domain();
357  Z3i::KSpace K;
358  bool space_ok = K.init( domain.lowerBound()-Z3i::Domain::Point::diagonal(),
359  domain.upperBound()+Z3i::Domain::Point::diagonal(), true );
360  if (!space_ok)
361  {
362  trace.error() << "Error in the Khalimsky space construction."<<std::endl;
363  return 2;
364  }
365  CanonicSCellEmbedder< KSpace > embedder( K );
366  SurfelAdjacency< Z3i::KSpace::dimension > Sadj( true );
367  trace.endBlock();
368  // Viewer settings
369 
370 
371  // Extraction of components
372  typedef KSpace::SurfelSet SurfelSet;
373  typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
374  typedef DigitalSurface< MySetOfSurfels > MyDigitalSurface;
375 
376  trace.beginBlock("Extracting surfaces");
377  std::vector< std::vector<SCell > > vectConnectedSCell;
378  Surfaces<KSpace>::extractAllConnectedSCell(vectConnectedSCell,K, Sadj, predicate, false);
379  std::ofstream outDat;
380  if( enable_dat )
381  {
382  trace.info() << "Exporting curvature as dat file: "<< export_dat_filename <<std::endl;
383  outDat.open( export_dat_filename.c_str() );
384  outDat << "# data exported from 3dCurvatureViewer implementing the II curvature estimator (Coeurjolly, D.; Lachaud, J.O; Levallois, J., (2013). Integral based Curvature"
385  << " Estimators in Digital Geometry. DGCI 2013.) " << std::endl;
386  outDat << "# format: surfel coordinates (in Khalimsky space) curvature: "<< mode << std::endl;
387  }
388 
389  trace.info()<<"Number of components= "<<vectConnectedSCell.size()<<std::endl;
390  trace.endBlock();
391 
392  if( vectConnectedSCell.size() == 0 )
393  {
394  trace.error()<< "No surface component exists. Please check the vol file threshold parameter.";
395  trace.info()<<std::endl;
396  exit(2);
397  }
398 
399 #ifdef WITH_VISU3D_QGLVIEWER
400  QApplication application( argc, argv );
401  typedef Viewer3D<Z3i::Space, Z3i::KSpace> Viewer;
402 #endif
403  typedef Board3D<Z3i::Space, Z3i::KSpace> Board;
404 
405 #ifdef WITH_VISU3D_QGLVIEWER
406  Viewer viewer( K );
407 #endif
408  Board board( K );
409 
410 #ifdef WITH_VISU3D_QGLVIEWER
411  if( enable_visu )
412  {
413  viewer.show();
414  }
415 #endif
416 
417  for( unsigned int i = 0; i<vectConnectedSCell.size(); ++i )
418  {
419  if( vectConnectedSCell[i].size() <= threshold )
420  {
421  continue;
422  }
423 
424  MySetOfSurfels aSet(K, Sadj);
425 
426  for( std::vector<SCell>::const_iterator it = vectConnectedSCell.at(i).begin();
427  it != vectConnectedSCell.at(i).end();
428  ++it )
429  {
430  aSet.surfelSet().insert( *it);
431  }
432 
433  MyDigitalSurface digSurf( aSet );
434 
435 
436  typedef DepthFirstVisitor<MyDigitalSurface> Visitor;
437  typedef GraphVisitorRange< Visitor > VisitorRange;
438  typedef VisitorRange::ConstIterator SurfelConstIterator;
439  VisitorRange range( new Visitor( digSurf, *digSurf.begin() ) );
440  SurfelConstIterator abegin = range.begin();
441  SurfelConstIterator aend = range.end();
442 
443  VisitorRange range2( new Visitor( digSurf, *digSurf.begin() ) );
444  SurfelConstIterator abegin2 = range2.begin();
445 
446  trace.beginBlock("Curvature computation on a component");
447  if( ( mode.compare("gaussian") == 0 ) || ( mode.compare("mean") == 0 )
448  || ( mode.compare("k1") == 0 ) || ( mode.compare("k2") == 0 ))
449  {
450  typedef double Quantity;
451  std::vector< Quantity > results;
452  std::back_insert_iterator< std::vector< Quantity > > resultsIterator( results );
453  if ( mode.compare("mean") == 0 )
454  {
455  typedef functors::IIMeanCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
456  typedef IntegralInvariantVolumeEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
457 
458  MyIICurvatureFunctor functor;
459  functor.init( h, re_convolution_kernel );
460 
461  MyIIEstimator estimator( functor );
462  estimator.attach( K, predicate );
463  estimator.setParams( re_convolution_kernel/h );
464  estimator.init( h, abegin, aend );
465 
466  estimator.eval( abegin, aend, resultsIterator );
467  }
468  else if ( mode.compare("gaussian") == 0 )
469  {
470  typedef functors::IIGaussianCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
471  typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
472 
473  MyIICurvatureFunctor functor;
474  functor.init( h, re_convolution_kernel );
475 
476  MyIIEstimator estimator( functor ); estimator.attach( K,
477  predicate ); estimator.setParams( re_convolution_kernel/h );
478  estimator.init( h, abegin, aend );
479 
480  estimator.eval( abegin, aend, resultsIterator );
481  }
482  else if ( mode.compare("k1") == 0 )
483  {
484  typedef functors::IIFirstPrincipalCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
485  typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
486 
487  MyIICurvatureFunctor functor;
488  functor.init( h, re_convolution_kernel );
489 
490  MyIIEstimator estimator( functor );
491  estimator.attach( K, predicate );
492  estimator.setParams( re_convolution_kernel/h );
493  estimator.init( h, abegin, aend );
494 
495  estimator.eval( abegin, aend, resultsIterator );
496  }
497  else if ( mode.compare("k2") == 0 )
498  {
499  typedef functors::IISecondPrincipalCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
500  typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
501 
502  MyIICurvatureFunctor functor;
503  functor.init( h, re_convolution_kernel );
504 
505  MyIIEstimator estimator( functor );
506  estimator.attach( K, predicate );
507  estimator.setParams( re_convolution_kernel/h );
508  estimator.init( h, abegin, aend );
509 
510  estimator.eval( abegin, aend, resultsIterator );
511  }
512  trace.endBlock();
513 
514 
515  // Drawing results
516  trace.beginBlock("Visualisation");
517  Quantity min = results[ 0 ];
518  Quantity max = results[ 0 ];
519  for ( unsigned int i = 1; i < results.size(); ++i )
520  {
521  if ( results[ i ] < min )
522  {
523  min = results[ i ];
524  }
525  else if ( results[ i ] > max )
526  {
527  max = results[ i ];
528  }
529  }
530  trace.info() << "Max value= "<<max<<" min value= "<<min<<std::endl;
531  ASSERT( min <= max );
532  typedef GradientColorMap< Quantity > Gradient;
533  Gradient cmap_grad( min, (max==min)? max+1: max );
534  cmap_grad.addColor( Color( 50, 50, 255 ) );
535  cmap_grad.addColor( Color( 255, 0, 0 ) );
536  cmap_grad.addColor( Color( 255, 255, 10 ) );
537 
538 #ifdef WITH_VISU3D_QGLVIEWER
539  if( enable_visu )
540  {
541  viewer << SetMode3D((*abegin2).className(), "Basic" );
542  }
543 #endif
544  if( enable_obj )
545  {
546  board << SetMode3D((K.unsigns(*abegin2)).className(), "Basic" );
547  }
548 
549 
550  for ( unsigned int i = 0; i < results.size(); ++i )
551  {
552 #ifdef WITH_VISU3D_QGLVIEWER
553  if( enable_visu )
554  {
555  viewer << CustomColors3D( Color::Black, cmap_grad( results[ i ] ));
556  viewer << *abegin2;
557  }
558 #endif
559 
560  if( enable_obj )
561  {
562  board << CustomColors3D( Color::Black, cmap_grad( results[ i ] ));
563  board << K.unsigns(*abegin2);
564  }
565 
566  if( enable_dat )
567  {
568  Point kCoords = K.uKCoords(K.unsigns(*abegin2));
569  outDat << kCoords[0] << " " << kCoords[1] << " " << kCoords[2] << " " << results[i] << std::endl;
570  }
571 
572  ++abegin2;
573  }
574  }
575  else
576  {
577  typedef Z3i::Space::RealVector Quantity;
578  std::vector< Quantity > results;
579  std::back_insert_iterator< std::vector< Quantity > > resultsIterator( results );
580 
581  if( mode.compare("prindir1") == 0 )
582  {
583  typedef functors::IIFirstPrincipalDirectionFunctor<Z3i::Space> MyIICurvatureFunctor;
584  typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
585 
586  MyIICurvatureFunctor functor;
587  functor.init( h, re_convolution_kernel );
588 
589  MyIIEstimator estimator( functor );
590  estimator.attach( K, predicate );
591  estimator.setParams( re_convolution_kernel/h );
592  estimator.init( h, abegin, aend );
593 
594  estimator.eval( abegin, aend, resultsIterator );
595  }
596  else if( mode.compare("prindir2") == 0 )
597  {
598  typedef functors::IISecondPrincipalDirectionFunctor<Z3i::Space> MyIICurvatureFunctor;
599  typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
600 
601  MyIICurvatureFunctor functor;
602  functor.init( h, re_convolution_kernel );
603 
604  MyIIEstimator estimator( functor );
605  estimator.attach( K, predicate );
606  estimator.setParams( re_convolution_kernel/h );
607  estimator.init( h, abegin, aend );
608 
609  estimator.eval( abegin, aend, resultsIterator );
610  }
611  else
612  if( mode.compare("normal") == 0 )
613  {
614  typedef functors::IINormalDirectionFunctor<Z3i::Space> MyIICurvatureFunctor;
615  typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
616 
617  MyIICurvatureFunctor functor;
618  functor.init( h, re_convolution_kernel );
619 
620  MyIIEstimator estimator( functor );
621  estimator.attach( K, predicate );
622  estimator.setParams( re_convolution_kernel/h );
623  estimator.init( h, abegin, aend );
624 
625  estimator.eval( abegin, aend, resultsIterator );
626  }
627 
628 
629 
631 
632 #ifdef WITH_VISU3D_QGLVIEWER
633  if( enable_visu )
634  {
635  viewer << SetMode3D(K.uCell( K.sKCoords(*abegin2) ).className(), "Basic" );
636  }
637 #endif
638 
639  if( enable_obj )
640  {
641  board << SetMode3D(K.uCell( K.sKCoords(*abegin2) ).className(), "Basic" );
642  }
643 
644  for ( unsigned int i = 0; i < results.size(); ++i )
645  {
646  DGtal::Dimension kDim = K.sOrthDir( *abegin2 );
647  SCell outer = K.sIndirectIncident( *abegin2, kDim);
648  if ( predicate(embedder(outer)) )
649  {
650  outer = K.sDirectIncident( *abegin2, kDim);
651  }
652 
653  Cell unsignedSurfel = K.uCell( K.sKCoords(*abegin2) );
654 
655 #ifdef WITH_VISU3D_QGLVIEWER
656  if( enable_visu )
657  {
658  viewer << CustomColors3D( DGtal::Color(255,255,255,255),
659  DGtal::Color(255,255,255,255))
660  << unsignedSurfel;
661  }
662 #endif
663 
664  if( enable_obj )
665  {
666  board << CustomColors3D( DGtal::Color(255,255,255,255),
667  DGtal::Color(255,255,255,255))
668  << unsignedSurfel;
669  }
670 
671  if( enable_dat )
672  {
673  Point kCoords = K.uKCoords(K.unsigns(*abegin2));
674  outDat << kCoords[0] << " " << kCoords[1] << " " << kCoords[2] << " "
675  << results[i][0] << " " << results[i][1] << " " << results[i][2]
676  << std::endl;
677  }
678 
679  RealPoint center = embedder( outer );
680 
681 #ifdef WITH_VISU3D_QGLVIEWER
682  if( enable_visu )
683  {
684  if( mode.compare("prindir1") == 0 )
685  {
686  viewer.setLineColor( AXIS_COLOR_BLUE );
687  }
688  else if( mode.compare("prindir2") == 0 )
689  {
690  viewer.setLineColor( AXIS_COLOR_RED );
691  }
692  else if( mode.compare("normal") == 0 )
693  {
694  viewer.setLineColor( AXIS_COLOR_GREEN );
695  }
696 
697 
698  viewer.addLine (
699  RealPoint(
700  center[0] - 0.5 * results[i][0],
701  center[1] - 0.5 * results[i][1],
702  center[2] - 0.5 * results[i][2]
703  ),
704  RealPoint(
705  center[0] + 0.5 * results[i][0],
706  center[1] + 0.5 * results[i][1],
707  center[2] + 0.5 * results[i][2]
708  ),
709  AXIS_LINESIZE );
710  }
711 #endif
712 
713  if( enable_obj )
714  {
715  if( mode.compare("prindir1") == 0 )
716  {
717  board.setFillColor( AXIS_COLOR_BLUE );
718  }
719  else if( mode.compare("prindir2") == 0 )
720  {
721  board.setFillColor( AXIS_COLOR_RED );
722  }
723  else if( mode.compare("normal") == 0 )
724  {
725  board.setFillColor( AXIS_COLOR_GREEN );
726  }
727 
728  board.addCylinder (
729  RealPoint(
730  center[0] - 0.5 * results[i][0],
731  center[1] - 0.5 * results[i][1],
732  center[2] - 0.5 * results[i][2]),
733  RealPoint(
734  center[0] + 0.5 * results[i][0],
735  center[1] + 0.5 * results[i][1],
736  center[2] + 0.5 * results[i][2]),
737  0.2 );
738  }
739 
740  ++abegin2;
741  }
742  }
743  trace.endBlock();
744  }
745 
746 #ifdef WITH_VISU3D_QGLVIEWER
747  if( enable_visu )
748  {
749  viewer << Viewer3D<>::updateDisplay;
750  }
751 #endif
752  if( enable_obj )
753  {
754  trace.info()<< "Exporting object: " << export_obj_filename << " ...";
755  board.saveOBJ(export_obj_filename,normalization);
756  trace.info() << "[done]" << std::endl;
757  }
758  if( enable_dat )
759  {
760  outDat.close();
761  }
762 
763 #ifdef WITH_VISU3D_QGLVIEWER
764  if( enable_visu )
765  {
766  return application.exec();
767  }
768 #endif
769 
770  return 0;
771 }
772