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