31 #include "DGtal/base/Common.h"
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>
47 #include "DGtal/images/ImageHelper.h"
48 #include "DGtal/topology/DigitalSurface.h"
49 #include "DGtal/graph/DepthFirstVisitor.h"
50 #include "DGtal/graph/GraphVisitorRange.h"
53 #include "DGtal/geometry/volumes/KanungoNoise.h"
56 #include "DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h"
57 #include "DGtal/geometry/surfaces/estimation/IntegralInvariantVolumeEstimator.h"
58 #include "DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h"
61 #include "DGtal/io/boards/Board3D.h"
62 #include "DGtal/io/colormaps/GradientColorMap.h"
64 #ifdef WITH_VISU3D_QGLVIEWER
65 #include "DGtal/io/viewers/Viewer3D.h"
68 using namespace DGtal;
69 using namespace functors;
146 const Color AXIS_COLOR_RED( 200, 20, 20, 255 );
147 const Color AXIS_COLOR_GREEN( 20, 200, 20, 255 );
148 const Color AXIS_COLOR_BLUE( 20, 20, 200, 255 );
149 const double AXIS_LINESIZE = 0.05;
158 void missingParam( std::string param )
160 trace.error() <<
" Parameter: " << param <<
" is required.";
161 trace.info() << std::endl;
165 int main(
int argc,
char** argv )
169 std::string inputFileName;
170 double re_convolution_kernel;
171 double noiseLevel {0.5};
172 unsigned int threshold {8};
173 int minImageThreshold {0};
174 int maxImageThreshold {255};
175 std::string mode {
"mean"};
176 std::string export_obj_filename;
177 std::string export_dat_filename;
178 bool exportOnly {
false};
179 std::vector< double> vectScale;
180 bool normalization {
false};
183 app.description(
"Visualisation of 3d curvature from .vol file using curvature from Integral Invariant\nBasic usage:\n \t3dCurvatureViewerNoise file.vol --radius 5 --mode mean --noise 0.5 \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");
187 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." )
189 ->check(CLI::ExistingFile);
191 app.add_option(
"--radius,-r", re_convolution_kernel,
"Kernel radius for IntegralInvariant" )
193 app.add_option(
"--noise,-k", noiseLevel,
"Level of Kanungo noise ]0;1[",
true);
194 app.add_option(
"--threshold,-t", threshold,
"Min size of SCell boundary of an object",
true);
195 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);
196 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);
197 app.add_option(
"--mode,-m", mode,
"type of output : mean, gaussian, k1, k2, prindir1, prindir2 or normal(default mean)",
true)
198 -> check(CLI::IsMember({
"mean",
"gaussian",
"k1",
"k2",
"prindir1",
"prindir2",
"normal" }));
199 app.add_option(
"--exportOBJ,-o", export_obj_filename,
"Export the scene to specified OBJ/MTL filename (extensions added).");
200 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." );
201 app.add_flag(
"--exportOnly", exportOnly,
"Used to only export the result without the 3d Visualisation (usefull for scripts).");
203 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.")
206 app.add_option(
"--normalization,-n",normalization,
"When exporting to OBJ, performs a normalization so that the geometry fits in [-1/2,1/2]^3") ;
208 app.get_formatter()->column_width(40);
209 CLI11_PARSE(app, argc, argv);
212 bool neededArgsGiven=
true;
214 if( noiseLevel < 0.0 || noiseLevel > 1.0 )
216 trace.error() <<
"The noise level should be in the interval: ]0, 1["<< std::endl;
220 #ifndef WITH_VISU3D_QGLVIEWER
221 bool enable_visu =
false;
223 bool enable_visu = !exportOnly;
225 bool enable_obj = export_obj_filename !=
"";
226 bool enable_dat = export_dat_filename !=
"";
228 if( !enable_visu && !enable_obj && !enable_dat )
230 #ifndef WITH_VISU3D_QGLVIEWER
231 trace.error() <<
"You should specify what you want to export with --export and/or --exportDat." << std::endl;
233 trace.error() <<
"You should specify what you want to export with --export and/or --exportDat, or remove --exportOnly." << std::endl;
235 neededArgsGiven =
false;
242 if( export_obj_filename.find(
".obj") == std::string::npos )
244 std::ostringstream oss;
245 oss << export_obj_filename <<
".obj" << std::endl;
246 export_obj_filename = oss.str();
250 std::vector< double > aGridSizeReSample;
251 if( vectScale.size() == 3)
253 aGridSizeReSample.push_back(1.0/vectScale.at(0));
254 aGridSizeReSample.push_back(1.0/vectScale.at(1));
255 aGridSizeReSample.push_back(1.0/vectScale.at(2));
260 aGridSizeReSample.push_back(1.0);
261 aGridSizeReSample.push_back(1.0);
262 aGridSizeReSample.push_back(1.0);
268 typedef Z3i::Space::RealPoint RealPoint;
269 typedef Z3i::Point Point;
270 typedef ImageSelector< Z3i::Domain, int>::Type Image;
271 typedef DGtal::functors::BasicDomainSubSampler< HyperRectDomain<SpaceND<3, int> >,
272 DGtal::int32_t,
double > ReSampler;
273 typedef DGtal::ConstImageAdapter<Image, Image::Domain, ReSampler,
274 Image::Value, DGtal::functors::Identity > SamplerImageAdapter;
275 typedef IntervalForegroundPredicate< SamplerImageAdapter > ImagePredicate;
276 typedef KanungoNoise< ImagePredicate, Z3i::Domain > KanungoPredicate;
277 typedef BinaryPointPredicate<DomainPredicate<Image::Domain>, KanungoPredicate, AndBoolFct2 > Predicate;
278 typedef Z3i::KSpace KSpace;
279 typedef KSpace::SCell SCell;
280 typedef KSpace::Cell Cell;
282 trace.beginBlock(
"Loading the file");
284 Image image = GenericReader<Image>::import( inputFileName );
286 PointVector<3,int> shiftVector3D( 0 ,0, 0 );
287 DGtal::functors::BasicDomainSubSampler< HyperRectDomain< SpaceND< 3, int > >,
288 DGtal::int32_t,
double > reSampler(image.domain(),
289 aGridSizeReSample, shiftVector3D);
290 const functors::Identity identityFunctor{};
291 SamplerImageAdapter sampledImage ( image, reSampler.getSubSampledDomain(), reSampler, identityFunctor );
292 ImagePredicate predicateIMG = ImagePredicate( sampledImage, minImageThreshold, maxImageThreshold );
293 KanungoPredicate noisifiedPredicateIMG( predicateIMG, sampledImage.domain(), noiseLevel );
294 DomainPredicate<Z3i::Domain> domainPredicate( sampledImage.domain() );
296 Predicate predicate(domainPredicate, noisifiedPredicateIMG, andF );
299 Z3i::Domain domain = sampledImage.domain();
301 bool space_ok = K.init( domain.lowerBound()-Z3i::Domain::Point::diagonal(),
302 domain.upperBound()+Z3i::Domain::Point::diagonal(),
true );
305 trace.error() <<
"Error in the Khalimsky space construction."<<std::endl;
308 CanonicSCellEmbedder< KSpace > embedder( K );
309 SurfelAdjacency< Z3i::KSpace::dimension > Sadj(
true );
315 typedef KSpace::SurfelSet SurfelSet;
316 typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
317 typedef DigitalSurface< MySetOfSurfels > MyDigitalSurface;
321 trace.beginBlock(
"Extracting surfaces");
322 std::vector< std::vector<SCell > > vectConnectedSCell;
323 Surfaces<KSpace>::extractAllConnectedSCell(vectConnectedSCell,K, Sadj, predicate,
false);
324 std::ofstream outDat;
327 trace.info() <<
"Exporting curvature as dat file: "<< export_dat_filename <<std::endl;
328 outDat.open( export_dat_filename.c_str() );
329 outDat <<
"# data exported from 3dCurvatureViewer implementing the II curvature estimator (Coeurjolly, D.; Lachaud, J.O; Levallois, J., (2013). Integral based Curvature"
330 <<
" Estimators in Digital Geometry. DGCI 2013.) " << std::endl;
331 outDat <<
"# format: surfel coordinates (in Khalimsky space) curvature: "<< mode << std::endl;
334 trace.info()<<
"Number of components= "<<vectConnectedSCell.size()<<std::endl;
337 if( vectConnectedSCell.size() == 0 )
339 trace.error()<<
"No surface component exists. Please check the vol file threshold parameter.";
340 trace.info()<<std::endl;
344 #ifdef WITH_VISU3D_QGLVIEWER
345 QApplication application( argc, argv );
346 typedef Viewer3D<Z3i::Space, Z3i::KSpace> Viewer;
348 typedef Board3D<Z3i::Space, Z3i::KSpace> Board;
350 #ifdef WITH_VISU3D_QGLVIEWER
355 #ifdef WITH_VISU3D_QGLVIEWER
363 unsigned int max_size = 0;
364 for(
unsigned int ii = 0; ii<vectConnectedSCell.size(); ++ii )
366 if( vectConnectedSCell[ii].size() <= threshold )
370 if( vectConnectedSCell[ii].size() > max_size )
372 max_size = vectConnectedSCell[ii].size();
377 MySetOfSurfels aSet(K, Sadj);
379 for( std::vector<SCell>::const_iterator it = vectConnectedSCell.at(i).begin();
380 it != vectConnectedSCell.at(i).end();
383 aSet.surfelSet().insert( *it);
386 MyDigitalSurface digSurf( aSet );
389 typedef DepthFirstVisitor<MyDigitalSurface> Visitor;
390 typedef GraphVisitorRange< Visitor > VisitorRange;
391 typedef VisitorRange::ConstIterator SurfelConstIterator;
392 VisitorRange range(
new Visitor( digSurf, *digSurf.begin() ) );
393 SurfelConstIterator abegin = range.begin();
394 SurfelConstIterator aend = range.end();
396 VisitorRange range2(
new Visitor( digSurf, *digSurf.begin() ) );
397 SurfelConstIterator abegin2 = range2.begin();
399 trace.beginBlock(
"Curvature computation on a component");
400 if( ( mode.compare(
"gaussian") == 0 ) || ( mode.compare(
"mean") == 0 )
401 || ( mode.compare(
"k1") == 0 ) || ( mode.compare(
"k2") == 0 ))
403 typedef double Quantity;
404 std::vector< Quantity > results;
405 std::back_insert_iterator< std::vector< Quantity > > resultsIterator( results );
406 if ( mode.compare(
"mean") == 0 )
408 typedef functors::IIMeanCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
409 typedef IntegralInvariantVolumeEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
411 MyIICurvatureFunctor functor;
412 functor.init( h, re_convolution_kernel );
414 MyIIEstimator estimator( functor );
415 estimator.attach( K, predicate );
416 estimator.setParams( re_convolution_kernel/h );
417 estimator.init( h, abegin, aend );
419 estimator.eval( abegin, aend, resultsIterator );
421 else if ( mode.compare(
"gaussian") == 0 )
423 typedef functors::IIGaussianCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
424 typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
426 MyIICurvatureFunctor functor;
427 functor.init( h, re_convolution_kernel );
429 MyIIEstimator estimator( functor ); estimator.attach( K,
430 predicate ); estimator.setParams( re_convolution_kernel/h );
431 estimator.init( h, abegin, aend );
433 estimator.eval( abegin, aend, resultsIterator );
435 else if ( mode.compare(
"k1") == 0 )
437 typedef functors::IIFirstPrincipalCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
438 typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
440 MyIICurvatureFunctor functor;
441 functor.init( h, re_convolution_kernel );
443 MyIIEstimator estimator( functor );
444 estimator.attach( K, predicate );
445 estimator.setParams( re_convolution_kernel/h );
446 estimator.init( h, abegin, aend );
448 estimator.eval( abegin, aend, resultsIterator );
450 else if ( mode.compare(
"k2") == 0 )
452 typedef functors::IISecondPrincipalCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
453 typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
455 MyIICurvatureFunctor functor;
456 functor.init( h, re_convolution_kernel );
458 MyIIEstimator estimator( functor );
459 estimator.attach( K, predicate );
460 estimator.setParams( re_convolution_kernel/h );
461 estimator.init( h, abegin, aend );
463 estimator.eval( abegin, aend, resultsIterator );
469 trace.beginBlock(
"Visualisation");
470 Quantity min = results[ 0 ];
471 Quantity max = results[ 0 ];
472 for (
unsigned int i = 1; i < results.size(); ++i )
474 if ( results[ i ] < min )
478 else if ( results[ i ] > max )
483 trace.info() <<
"Max value= "<<max<<
" min value= "<<min<<std::endl;
484 ASSERT( min <= max );
485 typedef GradientColorMap< Quantity > Gradient;
486 Gradient cmap_grad( min, (max==min)? max+1: max );
487 cmap_grad.addColor( Color( 50, 50, 255 ) );
488 cmap_grad.addColor( Color( 255, 0, 0 ) );
489 cmap_grad.addColor( Color( 255, 255, 10 ) );
491 #ifdef WITH_VISU3D_QGLVIEWER
494 viewer << SetMode3D((*abegin2).className(),
"Basic" );
499 board << SetMode3D((K.unsigns(*abegin2)).className(),
"Basic" );
503 for (
unsigned int i = 0; i < results.size(); ++i )
505 #ifdef WITH_VISU3D_QGLVIEWER
508 viewer << CustomColors3D( Color::Black, cmap_grad( results[ i ] ));
515 board << CustomColors3D( Color::Black, cmap_grad( results[ i ] ));
516 board << K.unsigns(*abegin2);
521 Point kCoords = K.uKCoords(K.unsigns(*abegin2));
522 outDat << kCoords[0] <<
" " << kCoords[1] <<
" " << kCoords[2] <<
" " << results[i] << std::endl;
530 typedef Z3i::Space::RealVector Quantity;
531 std::vector< Quantity > results;
532 std::back_insert_iterator< std::vector< Quantity > > resultsIterator( results );
534 if( mode.compare(
"prindir1") == 0 )
536 typedef functors::IIFirstPrincipalDirectionFunctor<Z3i::Space> MyIICurvatureFunctor;
537 typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
539 MyIICurvatureFunctor functor;
540 functor.init( h, re_convolution_kernel );
542 MyIIEstimator estimator( functor );
543 estimator.attach( K, predicate );
544 estimator.setParams( re_convolution_kernel/h );
545 estimator.init( h, abegin, aend );
547 estimator.eval( abegin, aend, resultsIterator );
549 else if( mode.compare(
"prindir2") == 0 )
551 typedef functors::IISecondPrincipalDirectionFunctor<Z3i::Space> MyIICurvatureFunctor;
552 typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
554 MyIICurvatureFunctor functor;
555 functor.init( h, re_convolution_kernel );
557 MyIIEstimator estimator( functor );
558 estimator.attach( K, predicate );
559 estimator.setParams( re_convolution_kernel/h );
560 estimator.init( h, abegin, aend );
562 estimator.eval( abegin, aend, resultsIterator );
563 }
else if( mode.compare(
"normal") == 0 )
565 typedef functors::IINormalDirectionFunctor<Z3i::Space> MyIICurvatureFunctor;
566 typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
568 MyIICurvatureFunctor functor;
569 functor.init( h, re_convolution_kernel );
571 MyIIEstimator estimator( functor );
572 estimator.attach( K, predicate );
573 estimator.setParams( re_convolution_kernel/h );
574 estimator.init( h, abegin, aend );
576 estimator.eval( abegin, aend, resultsIterator );
581 #ifdef WITH_VISU3D_QGLVIEWER
584 viewer << SetMode3D(K.uCell( K.sKCoords(*abegin2) ).className(),
"Basic" );
590 board << SetMode3D(K.uCell( K.sKCoords(*abegin2) ).className(),
"Basic" );
593 for (
unsigned int i = 0; i < results.size(); ++i )
595 DGtal::Dimension kDim = K.sOrthDir( *abegin2 );
596 SCell outer = K.sIndirectIncident( *abegin2, kDim);
597 if ( predicate(Z3i::Point(embedder(outer), functors::Round<>()) ))
599 outer = K.sDirectIncident( *abegin2, kDim);
602 Cell unsignedSurfel = K.uCell( K.sKCoords(*abegin2) );
604 #ifdef WITH_VISU3D_QGLVIEWER
607 viewer << CustomColors3D( DGtal::Color(255,255,255,255),
608 DGtal::Color(255,255,255,255))
615 board << CustomColors3D( DGtal::Color(255,255,255,255),
616 DGtal::Color(255,255,255,255))
622 Point kCoords = K.uKCoords(K.unsigns(*abegin2));
623 outDat << kCoords[0] <<
" " << kCoords[1] <<
" " << kCoords[2] <<
" "
624 << results[i][0] <<
" " << results[i][1] <<
" " << results[i][2]
628 RealPoint center = embedder( outer );
630 #ifdef WITH_VISU3D_QGLVIEWER
633 if( mode.compare(
"prindir1") == 0 )
635 viewer.setLineColor( AXIS_COLOR_BLUE );
637 else if( mode.compare(
"prindir2") == 0 )
639 viewer.setLineColor( AXIS_COLOR_RED );
641 else if( mode.compare(
"normal") == 0 )
643 viewer.setLineColor( AXIS_COLOR_GREEN );
648 center[0] - 0.5 * results[i][0],
649 center[1] - 0.5 * results[i][1],
650 center[2] - 0.5 * results[i][2]
653 center[0] + 0.5 * results[i][0],
654 center[1] + 0.5 * results[i][1],
655 center[2] + 0.5 * results[i][2]
663 if( mode.compare(
"prindir1") == 0 )
665 board.setFillColor( AXIS_COLOR_BLUE );
667 else if( mode.compare(
"prindir2") == 0 )
669 board.setFillColor( AXIS_COLOR_RED );
671 else if( mode.compare(
"normal") == 0 )
673 board.setFillColor( AXIS_COLOR_GREEN );
678 center[0] - 0.5 * results[i][0],
679 center[1] - 0.5 * results[i][1],
680 center[2] - 0.5 * results[i][2]),
682 center[0] + 0.5 * results[i][0],
683 center[1] + 0.5 * results[i][1],
684 center[2] + 0.5 * results[i][2]),
693 #ifdef WITH_VISU3D_QGLVIEWER
696 viewer << Viewer3D<>::updateDisplay;
701 trace.info()<<
"Exporting object: " << export_obj_filename <<
" ...";
702 board.saveOBJ(export_obj_filename,normalization);
703 trace.info() <<
"[done]" << std::endl;
710 #ifdef WITH_VISU3D_QGLVIEWER
713 return application.exec();