31 #include "DGtal/base/Common.h"    34 #include <boost/program_options/options_description.hpp>    35 #include <boost/program_options/parsers.hpp>    36 #include <boost/program_options/variables_map.hpp>    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>    48 #include "DGtal/images/ImageHelper.h"    49 #include "DGtal/topology/DigitalSurface.h"    50 #include "DGtal/graph/DepthFirstVisitor.h"    51 #include "DGtal/graph/GraphVisitorRange.h"    54 #include "DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h"    55 #include "DGtal/geometry/surfaces/estimation/IntegralInvariantVolumeEstimator.h"    56 #include "DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h"    59 #include "DGtal/io/boards/Board3D.h"    60 #include "DGtal/io/colormaps/GradientColorMap.h"    62 #ifdef WITH_VISU3D_QGLVIEWER    63 #include "DGtal/io/viewers/Viewer3D.h"    66 using namespace DGtal;
   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;
   172 void missingParam( std::string param )
   174   trace.error() << 
" Parameter: " << param << 
" is required.";
   175   trace.info() << std::endl;
   178 namespace po = boost::program_options;
   180 int main( 
int argc, 
char** argv )
   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") ;
   199   po::variables_map vm;
   202       po::store( po::parse_command_line( argc, argv, general_opt ), vm );
   204   catch( 
const std::exception & ex )
   207       trace.error() << 
" Error checking program options: " << ex.what() << std::endl;
   209   bool neededArgsGiven=
true;
   211   if (parseOK && !(vm.count(
"input"))){
   212     missingParam(
"--input");
   213     neededArgsGiven=
false;
   215   if (parseOK && !(vm.count(
"radius"))){
   216     missingParam(
"--radius");
   217     neededArgsGiven=
false;
   220   bool normalization = 
false;
   221   if (parseOK && vm.count(
"normalization"))
   222     normalization = 
true;
   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 ))
   232       trace.error() << 
" The selected mode ("<<mode << 
") is not defined."<<std::endl;
   235 #ifndef WITH_VISU3D_QGLVIEWER   236   bool enable_visu = 
false;
   238   bool enable_visu = !vm.count(
"exportOnly"); 
   240   bool enable_obj = vm.count(
"exportOBJ"); 
   241   bool enable_dat = vm.count(
"exportDAT"); 
   243   if( !enable_visu && !enable_obj && !enable_dat )
   245 #ifndef WITH_VISU3D_QGLVIEWER   246       trace.error() << 
"You should specify what you want to export with --export and/or --exportDat." << std::endl;
   248       trace.error() << 
"You should specify what you want to export with --export and/or --exportDat, or remove --exportOnly." << std::endl;
   250       neededArgsGiven = 
false;
   253   if(!neededArgsGiven || !parseOK || vm.count(
"help") || argc <= 1 )
   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
   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
   271   unsigned int threshold = vm[
"threshold"].as< 
unsigned int >();
   272   int minImageThreshold =  vm[
"minImageThreshold"].as<  
int >();
   273   int maxImageThreshold =  vm[
"maxImageThreshold"].as<  
int >();
   277   std::string export_obj_filename;
   278   std::string export_dat_filename;
   282       export_obj_filename = vm[
"exportOBJ"].as< std::string >();
   283       if( export_obj_filename.find(
".obj") == std::string::npos )
   285           std::ostringstream oss;
   286           oss << export_obj_filename << 
".obj" << std::endl;
   287           export_obj_filename = oss.str();
   294       export_dat_filename = vm[
"exportDAT"].as<std::string>();
   297   double re_convolution_kernel = vm[
"radius"].as< 
double >();
   300   std::vector<  double > aGridSizeReSample;
   301   if( vm.count( 
"imageScale" ))
   303       std::vector< double> vectScale = vm[
"imageScale"].as<std::vector<double > >();
   304       if( vectScale.size() != 3 )
   306           trace.error() << 
"The grid size should contains 3 elements" << std::endl;
   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));
   318       aGridSizeReSample.push_back(1.0);
   319       aGridSizeReSample.push_back(1.0);
   320       aGridSizeReSample.push_back(1.0);
   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;
   340   trace.beginBlock(
"Loading the file");
   341   std::string filename = vm[
"input"].as< std::string >();
   342   Image image = GenericReader<Image>::import( filename );
   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() );
   353   Predicate predicate(domainPredicate, predicateIMG, andF  );
   356   Z3i::Domain domain =  sampledImage.domain();
   358   bool space_ok = K.init( domain.lowerBound()-Z3i::Domain::Point::diagonal(),
   359                           domain.upperBound()+Z3i::Domain::Point::diagonal(), true );
   362       trace.error() << 
"Error in the Khalimsky space construction."<<std::endl;
   365   CanonicSCellEmbedder< KSpace > embedder( K );
   366   SurfelAdjacency< Z3i::KSpace::dimension > Sadj( 
true );
   372   typedef KSpace::SurfelSet SurfelSet;
   373   typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
   374   typedef DigitalSurface< MySetOfSurfels > MyDigitalSurface;
   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;
   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;
   389   trace.info()<<
"Number of components= "<<vectConnectedSCell.size()<<std::endl;
   392   if( vectConnectedSCell.size() == 0 )
   394       trace.error()<< 
"No surface component exists. Please check the vol file threshold parameter.";
   395       trace.info()<<std::endl;
   399 #ifdef WITH_VISU3D_QGLVIEWER   400   QApplication application( argc, argv );
   401   typedef Viewer3D<Z3i::Space, Z3i::KSpace> Viewer;
   403   typedef Board3D<Z3i::Space, Z3i::KSpace> Board;
   405 #ifdef WITH_VISU3D_QGLVIEWER   410 #ifdef WITH_VISU3D_QGLVIEWER   417   for( 
unsigned int i = 0; i<vectConnectedSCell.size(); ++i )
   419       if( vectConnectedSCell[i].size() <= threshold )
   424       MySetOfSurfels  aSet(K, Sadj);
   426       for( std::vector<SCell>::const_iterator it = vectConnectedSCell.at(i).begin();
   427            it != vectConnectedSCell.at(i).end();
   430           aSet.surfelSet().insert( *it);
   433       MyDigitalSurface digSurf( aSet );
   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();
   443       VisitorRange range2( 
new Visitor( digSurf, *digSurf.begin() ) );
   444       SurfelConstIterator abegin2 = range2.begin();
   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 ))
   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 )
   455               typedef functors::IIMeanCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
   456               typedef IntegralInvariantVolumeEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
   458               MyIICurvatureFunctor functor;
   459               functor.init( h, re_convolution_kernel );
   461               MyIIEstimator estimator( functor );
   462               estimator.attach( K, predicate );
   463               estimator.setParams( re_convolution_kernel/h );
   464               estimator.init( h, abegin, aend );
   466               estimator.eval( abegin, aend, resultsIterator );
   468           else if ( mode.compare(
"gaussian") == 0 )
   470               typedef functors::IIGaussianCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
   471               typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
   473               MyIICurvatureFunctor functor;
   474               functor.init( h, re_convolution_kernel );
   476               MyIIEstimator estimator( functor ); estimator.attach( K,
   477                                                                     predicate ); estimator.setParams( re_convolution_kernel/h );
   478               estimator.init( h, abegin, aend );
   480               estimator.eval( abegin, aend, resultsIterator );
   482           else if ( mode.compare(
"k1") == 0 )
   484               typedef functors::IIFirstPrincipalCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
   485               typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
   487               MyIICurvatureFunctor functor;
   488               functor.init( h, re_convolution_kernel );
   490               MyIIEstimator estimator( functor );
   491               estimator.attach( K, predicate );
   492               estimator.setParams( re_convolution_kernel/h );
   493               estimator.init( h, abegin, aend );
   495               estimator.eval( abegin, aend, resultsIterator );
   497           else if ( mode.compare(
"k2") == 0 )
   499               typedef functors::IISecondPrincipalCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
   500               typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
   502               MyIICurvatureFunctor functor;
   503               functor.init( h, re_convolution_kernel );
   505               MyIIEstimator estimator( functor );
   506               estimator.attach( K, predicate );
   507               estimator.setParams( re_convolution_kernel/h );
   508               estimator.init( h, abegin, aend );
   510               estimator.eval( abegin, aend, resultsIterator );
   516           trace.beginBlock(
"Visualisation");
   517           Quantity min = results[ 0 ];
   518           Quantity max = results[ 0 ];
   519           for ( 
unsigned int i = 1; i < results.size(); ++i )
   521               if ( results[ i ] < min )
   525               else if ( results[ i ] > max )
   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 ) );
   538 #ifdef WITH_VISU3D_QGLVIEWER   541               viewer << SetMode3D((*abegin2).className(), 
"Basic" );
   546               board << SetMode3D((K.unsigns(*abegin2)).className(), 
"Basic" );
   550           for ( 
unsigned int i = 0; i < results.size(); ++i )
   552 #ifdef WITH_VISU3D_QGLVIEWER   555                   viewer << CustomColors3D( Color::Black, cmap_grad( results[ i ] ));
   562                   board << CustomColors3D( Color::Black, cmap_grad( results[ i ] ));
   563                   board      << K.unsigns(*abegin2);
   568                   Point kCoords = K.uKCoords(K.unsigns(*abegin2));
   569                   outDat << kCoords[0] << 
" " << kCoords[1] << 
" " << kCoords[2] <<  
" " <<  results[i] << std::endl;
   577           typedef Z3i::Space::RealVector Quantity;
   578           std::vector< Quantity > results;
   579           std::back_insert_iterator< std::vector< Quantity > > resultsIterator( results );
   581           if( mode.compare(
"prindir1") == 0 )
   583               typedef functors::IIFirstPrincipalDirectionFunctor<Z3i::Space> MyIICurvatureFunctor;
   584               typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
   586               MyIICurvatureFunctor functor;
   587               functor.init( h, re_convolution_kernel );
   589               MyIIEstimator estimator( functor );
   590               estimator.attach( K, predicate );
   591               estimator.setParams( re_convolution_kernel/h );
   592               estimator.init( h, abegin, aend );
   594               estimator.eval( abegin, aend, resultsIterator );
   596           else if( mode.compare(
"prindir2") == 0 )
   598               typedef functors::IISecondPrincipalDirectionFunctor<Z3i::Space> MyIICurvatureFunctor;
   599               typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
   601               MyIICurvatureFunctor functor;
   602               functor.init( h, re_convolution_kernel );
   604               MyIIEstimator estimator( functor );
   605               estimator.attach( K, predicate );
   606               estimator.setParams( re_convolution_kernel/h );
   607               estimator.init( h, abegin, aend );
   609               estimator.eval( abegin, aend, resultsIterator );
   612             if( mode.compare(
"normal") == 0 )
   614                 typedef functors::IINormalDirectionFunctor<Z3i::Space> MyIICurvatureFunctor;
   615                 typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
   617                 MyIICurvatureFunctor functor;
   618                 functor.init( h, re_convolution_kernel );
   620                 MyIIEstimator estimator( functor );
   621                 estimator.attach( K, predicate );
   622                 estimator.setParams( re_convolution_kernel/h );
   623                 estimator.init( h, abegin, aend );
   625                 estimator.eval( abegin, aend, resultsIterator );
   632 #ifdef WITH_VISU3D_QGLVIEWER   635               viewer << SetMode3D(K.uCell( K.sKCoords(*abegin2) ).className(), 
"Basic" );
   641               board << SetMode3D(K.uCell( K.sKCoords(*abegin2) ).className(), 
"Basic" );
   644           for ( 
unsigned int i = 0; i < results.size(); ++i )
   646               DGtal::Dimension kDim = K.sOrthDir( *abegin2 );
   647               SCell outer = K.sIndirectIncident( *abegin2, kDim);
   648               if ( predicate(embedder(outer)) )
   650                   outer = K.sDirectIncident( *abegin2, kDim);
   653               Cell unsignedSurfel = K.uCell( K.sKCoords(*abegin2) );
   655 #ifdef WITH_VISU3D_QGLVIEWER   658                   viewer << CustomColors3D( DGtal::Color(255,255,255,255),
   659                                             DGtal::Color(255,255,255,255))
   666                   board << CustomColors3D( DGtal::Color(255,255,255,255),
   667                                            DGtal::Color(255,255,255,255))
   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]
   679               RealPoint center = embedder( outer );
   681 #ifdef WITH_VISU3D_QGLVIEWER   684                   if( mode.compare(
"prindir1") == 0 )
   686                       viewer.setLineColor( AXIS_COLOR_BLUE );
   688                   else if( mode.compare(
"prindir2") == 0 )
   690                       viewer.setLineColor( AXIS_COLOR_RED );
   692                   else if( mode.compare(
"normal") == 0 )
   694                       viewer.setLineColor( AXIS_COLOR_GREEN );
   700                                             center[0] -  0.5 * results[i][0],
   701                                             center[1] -  0.5 * results[i][1],
   702                                             center[2] -  0.5 * results[i][2]
   705                                             center[0] +  0.5 * results[i][0],
   706                                             center[1] +  0.5 * results[i][1],
   707                                             center[2] +  0.5 * results[i][2]
   715                   if( mode.compare(
"prindir1") == 0 )
   717                       board.setFillColor( AXIS_COLOR_BLUE );
   719                   else if( mode.compare(
"prindir2") == 0 )
   721                       board.setFillColor( AXIS_COLOR_RED );
   723                   else if( mode.compare(
"normal") == 0 )
   725                       board.setFillColor( AXIS_COLOR_GREEN );
   730                                                center[0] -  0.5 * results[i][0],
   731                                                center[1] -  0.5 * results[i][1],
   732                                                center[2] -  0.5 * results[i][2]),
   734                                                center[0] +  0.5 * results[i][0],
   735                                                center[1] +  0.5 * results[i][1],
   736                                                center[2] +  0.5 * results[i][2]),
   746 #ifdef WITH_VISU3D_QGLVIEWER   749       viewer << Viewer3D<>::updateDisplay;
   754       trace.info()<< 
"Exporting object: " << export_obj_filename << 
" ...";
   755       board.saveOBJ(export_obj_filename,normalization);
   756       trace.info() << 
"[done]" << std::endl;
   763 #ifdef WITH_VISU3D_QGLVIEWER   766       return application.exec();