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();