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