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/volumes/KanungoNoise.h"    57 #include "DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h"    58 #include "DGtal/geometry/surfaces/estimation/IntegralInvariantVolumeEstimator.h"    59 #include "DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h"    62 #include "DGtal/io/boards/Board3D.h"    63 #include "DGtal/io/colormaps/GradientColorMap.h"    65 #ifdef WITH_VISU3D_QGLVIEWER    66 #include "DGtal/io/viewers/Viewer3D.h"    69 using namespace DGtal;
   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;
   177 void missingParam( std::string param )
   179   trace.error() << 
" Parameter: " << param << 
" is required.";
   180   trace.info() << std::endl;
   183 namespace po = boost::program_options;
   185 int main( 
int argc, 
char** argv )
   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") ;
   205   po::variables_map vm;
   208       po::store( po::parse_command_line( argc, argv, general_opt ), vm );
   210   catch( 
const std::exception & ex )
   213       trace.error() << 
" Error checking program options: " << ex.what() << std::endl;
   215   bool neededArgsGiven=
true;
   217   if (parseOK && !(vm.count(
"input"))){
   218     missingParam(
"--input");
   219     neededArgsGiven=
false;
   221   if (parseOK && !(vm.count(
"radius"))){
   222     missingParam(
"--radius");
   223     neededArgsGiven=
false;
   226   bool normalization = 
false;
   227   if (parseOK && vm.count(
"normalization"))
   228     normalization = 
true;
   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 ))
   238       trace.error() << 
" The selected mode ("<<mode << 
") is not defined."<<std::endl;
   241   double noiseLevel = vm[
"noise"].as< 
double >();
   242   if( noiseLevel < 0.0 || noiseLevel > 1.0 )
   245       trace.error() << 
"The noise level should be in the interval: ]0, 1["<< std::endl;
   248 #ifndef WITH_VISU3D_QGLVIEWER   249   bool enable_visu = 
false;
   251   bool enable_visu = !vm.count(
"exportOnly"); 
   253   bool enable_obj = vm.count(
"exportOBJ"); 
   254   bool enable_dat = vm.count(
"exportDAT"); 
   256   if( !enable_visu && !enable_obj && !enable_dat )
   258 #ifndef WITH_VISU3D_QGLVIEWER   259       trace.error() << 
"You should specify what you want to export with --export and/or --exportDat." << std::endl;
   261       trace.error() << 
"You should specify what you want to export with --export and/or --exportDat, or remove --exportOnly." << std::endl;
   263       neededArgsGiven = 
false;
   266   if(!neededArgsGiven || !parseOK || vm.count(
"help") || argc <= 1 )
   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
   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
   284   unsigned int threshold = vm[
"threshold"].as< 
unsigned int >();
   285   int minImageThreshold =  vm[
"minImageThreshold"].as<  
int >();
   286   int maxImageThreshold =  vm[
"maxImageThreshold"].as<  
int >();
   290   std::string export_obj_filename;
   291   std::string export_dat_filename;
   295       export_obj_filename = vm[
"exportOBJ"].as< std::string >();
   296       if( export_obj_filename.find(
".obj") == std::string::npos )
   298           std::ostringstream oss;
   299           oss << export_obj_filename << 
".obj" << std::endl;
   300           export_obj_filename = oss.str();
   307       export_dat_filename = vm[
"exportDAT"].as<std::string>();
   310   double re_convolution_kernel = vm[
"radius"].as< 
double >();
   313   std::vector<  double > aGridSizeReSample;
   314   if( vm.count( 
"imageScale" ))
   316       std::vector< double> vectScale = vm[
"imageScale"].as<std::vector<double > >();
   317       if( vectScale.size() != 3 )
   319           trace.error() << 
"The grid size should contains 3 elements" << std::endl;
   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));
   331       aGridSizeReSample.push_back(1.0);
   332       aGridSizeReSample.push_back(1.0);
   333       aGridSizeReSample.push_back(1.0);
   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;
   354   trace.beginBlock(
"Loading the file");
   355   std::string filename = vm[
"input"].as< std::string >();
   356   Image image = GenericReader<Image>::import( filename );
   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() );
   368   Predicate predicate(domainPredicate, noisifiedPredicateIMG, andF  );
   371   Z3i::Domain domain =  sampledImage.domain();
   373   bool space_ok = K.init( domain.lowerBound()-Z3i::Domain::Point::diagonal(),
   374                           domain.upperBound()+Z3i::Domain::Point::diagonal(), true );
   377       trace.error() << 
"Error in the Khalimsky space construction."<<std::endl;
   380   CanonicSCellEmbedder< KSpace > embedder( K );
   381   SurfelAdjacency< Z3i::KSpace::dimension > Sadj( 
true );
   387   typedef KSpace::SurfelSet SurfelSet;
   388   typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
   389   typedef DigitalSurface< MySetOfSurfels > MyDigitalSurface;
   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;
   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;
   406   trace.info()<<
"Number of components= "<<vectConnectedSCell.size()<<std::endl;
   409   if( vectConnectedSCell.size() == 0 )
   411       trace.error()<< 
"No surface component exists. Please check the vol file threshold parameter.";
   412       trace.info()<<std::endl;
   416 #ifdef WITH_VISU3D_QGLVIEWER   417   QApplication application( argc, argv );
   418   typedef Viewer3D<Z3i::Space, Z3i::KSpace> Viewer;
   420   typedef Board3D<Z3i::Space, Z3i::KSpace> Board;
   422 #ifdef WITH_VISU3D_QGLVIEWER   427 #ifdef WITH_VISU3D_QGLVIEWER   435   unsigned int max_size = 0;
   436   for( 
unsigned int ii = 0; ii<vectConnectedSCell.size(); ++ii )
   438       if( vectConnectedSCell[ii].size() <= threshold )
   442       if( vectConnectedSCell[ii].size() > max_size )
   444           max_size = vectConnectedSCell[ii].size();
   449   MySetOfSurfels  aSet(K, Sadj);
   451   for( std::vector<SCell>::const_iterator it = vectConnectedSCell.at(i).begin();
   452        it != vectConnectedSCell.at(i).end();
   455       aSet.surfelSet().insert( *it);
   458   MyDigitalSurface digSurf( aSet );
   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();
   468   VisitorRange range2( 
new Visitor( digSurf, *digSurf.begin() ) );
   469   SurfelConstIterator abegin2 = range2.begin();
   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 ))
   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 )
   480           typedef functors::IIMeanCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
   481           typedef IntegralInvariantVolumeEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
   483           MyIICurvatureFunctor functor;
   484           functor.init( h, re_convolution_kernel );
   486           MyIIEstimator estimator( functor );
   487           estimator.attach( K, predicate );
   488           estimator.setParams( re_convolution_kernel/h );
   489           estimator.init( h, abegin, aend );
   491           estimator.eval( abegin, aend, resultsIterator );
   493       else if ( mode.compare(
"gaussian") == 0 )
   495           typedef functors::IIGaussianCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
   496           typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
   498           MyIICurvatureFunctor functor;
   499           functor.init( h, re_convolution_kernel );
   501           MyIIEstimator estimator( functor ); estimator.attach( K,
   502                                                                 predicate ); estimator.setParams( re_convolution_kernel/h );
   503           estimator.init( h, abegin, aend );
   505           estimator.eval( abegin, aend, resultsIterator );
   507       else if ( mode.compare(
"k1") == 0 )
   509           typedef functors::IIFirstPrincipalCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
   510           typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
   512           MyIICurvatureFunctor functor;
   513           functor.init( h, re_convolution_kernel );
   515           MyIIEstimator estimator( functor );
   516           estimator.attach( K, predicate );
   517           estimator.setParams( re_convolution_kernel/h );
   518           estimator.init( h, abegin, aend );
   520           estimator.eval( abegin, aend, resultsIterator );
   522       else if ( mode.compare(
"k2") == 0 )
   524           typedef functors::IISecondPrincipalCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
   525           typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
   527           MyIICurvatureFunctor functor;
   528           functor.init( h, re_convolution_kernel );
   530           MyIIEstimator estimator( functor );
   531           estimator.attach( K, predicate );
   532           estimator.setParams( re_convolution_kernel/h );
   533           estimator.init( h, abegin, aend );
   535           estimator.eval( abegin, aend, resultsIterator );
   541       trace.beginBlock(
"Visualisation");
   542       Quantity min = results[ 0 ];
   543       Quantity max = results[ 0 ];
   544       for ( 
unsigned int i = 1; i < results.size(); ++i )
   546           if ( results[ i ] < min )
   550           else if ( results[ i ] > max )
   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 ) );
   563 #ifdef WITH_VISU3D_QGLVIEWER   566           viewer << SetMode3D((*abegin2).className(), 
"Basic" );
   571           board << SetMode3D((K.unsigns(*abegin2)).className(), 
"Basic" );
   575       for ( 
unsigned int i = 0; i < results.size(); ++i )
   577 #ifdef WITH_VISU3D_QGLVIEWER   580               viewer << CustomColors3D( Color::Black, cmap_grad( results[ i ] ));
   587               board << CustomColors3D( Color::Black, cmap_grad( results[ i ] ));
   588               board      << K.unsigns(*abegin2);
   593               Point kCoords = K.uKCoords(K.unsigns(*abegin2));
   594               outDat << kCoords[0] << 
" " << kCoords[1] << 
" " << kCoords[2] <<  
" " <<  results[i] << std::endl;
   602       typedef Z3i::Space::RealVector Quantity;
   603       std::vector< Quantity > results;
   604       std::back_insert_iterator< std::vector< Quantity > > resultsIterator( results );
   606       if( mode.compare(
"prindir1") == 0 )
   608           typedef functors::IIFirstPrincipalDirectionFunctor<Z3i::Space> MyIICurvatureFunctor;
   609           typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
   611           MyIICurvatureFunctor functor;
   612           functor.init( h, re_convolution_kernel );
   614           MyIIEstimator estimator( functor );
   615           estimator.attach( K, predicate );
   616           estimator.setParams( re_convolution_kernel/h );
   617           estimator.init( h, abegin, aend );
   619           estimator.eval( abegin, aend, resultsIterator );
   621       else if( mode.compare(
"prindir2") == 0 )
   623           typedef functors::IISecondPrincipalDirectionFunctor<Z3i::Space> MyIICurvatureFunctor;
   624           typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
   626           MyIICurvatureFunctor functor;
   627           functor.init( h, re_convolution_kernel );
   629           MyIIEstimator estimator( functor );
   630           estimator.attach( K, predicate );
   631           estimator.setParams( re_convolution_kernel/h );
   632           estimator.init( h, abegin, aend );
   634           estimator.eval( abegin, aend, resultsIterator );
   635         }
else if( mode.compare(
"normal") == 0 )
   637           typedef functors::IINormalDirectionFunctor<Z3i::Space> MyIICurvatureFunctor;
   638           typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
   640           MyIICurvatureFunctor functor;
   641           functor.init( h, re_convolution_kernel );
   643           MyIIEstimator estimator( functor );
   644           estimator.attach( K, predicate );
   645           estimator.setParams( re_convolution_kernel/h );
   646           estimator.init( h, abegin, aend );
   648           estimator.eval( abegin, aend, resultsIterator );
   653 #ifdef WITH_VISU3D_QGLVIEWER   656           viewer << SetMode3D(K.uCell( K.sKCoords(*abegin2) ).className(), 
"Basic" );
   662           board << SetMode3D(K.uCell( K.sKCoords(*abegin2) ).className(), 
"Basic" );
   665       for ( 
unsigned int i = 0; i < results.size(); ++i )
   667           DGtal::Dimension kDim = K.sOrthDir( *abegin2 );
   668           SCell outer = K.sIndirectIncident( *abegin2, kDim);
   669           if ( predicate(embedder(outer)) )
   671               outer = K.sDirectIncident( *abegin2, kDim);
   674           Cell unsignedSurfel = K.uCell( K.sKCoords(*abegin2) );
   676 #ifdef WITH_VISU3D_QGLVIEWER   679               viewer << CustomColors3D( DGtal::Color(255,255,255,255),
   680                                         DGtal::Color(255,255,255,255))
   687               board << CustomColors3D( DGtal::Color(255,255,255,255),
   688                                        DGtal::Color(255,255,255,255))
   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]
   700           RealPoint center = embedder( outer );
   702 #ifdef WITH_VISU3D_QGLVIEWER   705               if( mode.compare(
"prindir1") == 0 )
   707                   viewer.setLineColor( AXIS_COLOR_BLUE );
   709               else if( mode.compare(
"prindir2") == 0 )
   711                   viewer.setLineColor( AXIS_COLOR_RED );
   713               else if( mode.compare(
"normal") == 0 )
   715                   viewer.setLineColor( AXIS_COLOR_GREEN );
   720                                         center[0] -  0.5 * results[i][0],
   721                                         center[1] -  0.5 * results[i][1],
   722                                         center[2] -  0.5 * results[i][2]
   725                                         center[0] +  0.5 * results[i][0],
   726                                         center[1] +  0.5 * results[i][1],
   727                                         center[2] +  0.5 * results[i][2]
   735               if( mode.compare(
"prindir1") == 0 )
   737                   board.setFillColor( AXIS_COLOR_BLUE );
   739               else if( mode.compare(
"prindir2") == 0 )
   741                   board.setFillColor( AXIS_COLOR_RED );
   743               else if( mode.compare(
"normal") == 0 )
   745                   board.setFillColor( AXIS_COLOR_GREEN );
   750                                            center[0] -  0.5 * results[i][0],
   751                                            center[1] -  0.5 * results[i][1],
   752                                            center[2] -  0.5 * results[i][2]),
   754                                            center[0] +  0.5 * results[i][0],
   755                                            center[1] +  0.5 * results[i][1],
   756                                            center[2] +  0.5 * results[i][2]),
   765 #ifdef WITH_VISU3D_QGLVIEWER   768       viewer << Viewer3D<>::updateDisplay;
   773       trace.info()<< 
"Exporting object: " << export_obj_filename << 
" ...";
   774       board.saveOBJ(export_obj_filename,normalization);
   775       trace.info() << 
"[done]" << std::endl;
   782 #ifdef WITH_VISU3D_QGLVIEWER   785       return application.exec();