DGtalTools 2.1.0
Loading...
Searching...
No Matches
3dCurvatureViewerNoise.cpp
1
30#include <iostream>
31#include <sstream>
32#include "DGtal/base/Common.h"
33#include <cstring>
34
35#include "CLI11.hpp"
36
37
38// Shape constructors
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>
47
48#include "DGtal/images/ImageHelper.h"
49#include "DGtal/topology/DigitalSurface.h"
50#include "DGtal/graph/DepthFirstVisitor.h"
51#include "DGtal/graph/GraphVisitorRange.h"
52
53// Noise
54#include "DGtal/geometry/volumes/KanungoNoise.h"
55
56// Integral Invariant includes
57#include "DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h"
58#include "DGtal/geometry/surfaces/estimation/IntegralInvariantVolumeEstimator.h"
59#include "DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h"
60
61// Drawing
62#include "DGtal/io/viewers/PolyscopeViewer.h"
63
64using namespace DGtal;
65using namespace functors;
66
141const Color AXIS_COLOR_RED( 200, 20, 20, 255 );
142const Color AXIS_COLOR_GREEN( 20, 200, 20, 255 );
143const Color AXIS_COLOR_BLUE( 20, 20, 200, 255 );
144const double AXIS_LINESIZE = 0.05;
145
147
148// used to enable/disable the polyscope UI
149bool show_ui = false;
150void myCallback() {
151 ImGuiIO& io = ImGui::GetIO();
152 if (ImGui::IsKeyPressed(ImGuiKey_W)) {
153 show_ui = !show_ui;
154 polyscope::options::buildGui = show_ui;
155 }
156}
157
158
164void missingParam( std::string param )
165{
166 trace.error() << " Parameter: " << param << " is required.";
167 trace.info() << std::endl;
168}
169
170
171int main( int argc, char** argv )
172{
173 // parse command line using CLI ----------------------------------------------
174 CLI::App app;
175 std::string inputFileName;
176 double re_convolution_kernel;
177 double noiseLevel {0.5};
178 unsigned int threshold {8};
179 int minImageThreshold {0};
180 int maxImageThreshold {255};
181 std::string mode {"mean"};
182 std::string export_dat_filename;
183 bool exportOnly {false};
184 std::vector< double> vectScale;
185
186
187
188 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\n Example: 3dCurvatureViewerNoise $DGtal/examples/samples/Al.100.vol -r 15 -m mean -k 0.5");
189
190
191
192 app.add_option("-i,--input,1", inputFileName, "vol file (.vol, .longvol .p3d, .pgm3d and if DGTAL_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." )
193 ->required()
194 ->check(CLI::ExistingFile);
195
196 app.add_option("--radius,-r", re_convolution_kernel, "Kernel radius for IntegralInvariant" )
197 ->required();
198 app.add_option("--noise,-k", noiseLevel, "Level of Kanungo noise ]0;1[");
199 app.add_option("--threshold,-t", threshold, "Min size of SCell boundary of an object");
200 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 ] ).");
201 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] ).");
202 app.add_option("--mode,-m", mode, "type of output : mean, gaussian, k1, k2, prindir1, prindir2 or normal(default mean)")
203 -> check(CLI::IsMember({"mean","gaussian", "k1", "k2", "prindir1","prindir2", "normal" }));
204 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." );
205 app.add_flag("--exportOnly", exportOnly, "Used to only export the result without the 3d Visualisation (usefull for scripts).");
206
207 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.")
208 ->expected(3);
209
210 app.get_formatter()->column_width(40);
211 CLI11_PARSE(app, argc, argv);
212 // END parse command line using CLI ----------------------------------------------
213
214 bool neededArgsGiven=true;
215
216 if( noiseLevel < 0.0 || noiseLevel > 1.0 )
217 {
218 trace.error() << "The noise level should be in the interval: ]0, 1["<< std::endl;
219 exit(EXIT_FAILURE);
220 }
221 bool enable_visu = !exportOnly;
222 bool enable_dat = export_dat_filename != "";
223
224 if( !enable_visu && !enable_dat )
225 {
226 trace.error() << "You should specify what you want to export with --export and/or --exportDat, or remove --exportOnly." << std::endl;
227 neededArgsGiven = false;
228 }
229
230 double h = 1.0;
231
232 std::vector< double > aGridSizeReSample;
233 if( vectScale.size() == 3)
234 {
235 aGridSizeReSample.push_back(1.0/vectScale.at(0));
236 aGridSizeReSample.push_back(1.0/vectScale.at(1));
237 aGridSizeReSample.push_back(1.0/vectScale.at(2));
238
239 }
240 else
241 {
242 aGridSizeReSample.push_back(1.0);
243 aGridSizeReSample.push_back(1.0);
244 aGridSizeReSample.push_back(1.0);
245 }
246
247
248
249 // Construction of the shape from vol file
250 typedef Z3i::Space::RealPoint RealPoint;
251 typedef Z3i::Point Point;
252 typedef ImageSelector< Z3i::Domain, int>::Type Image;
253 typedef DGtal::functors::BasicDomainSubSampler< HyperRectDomain<SpaceND<3, int> >,
254 DGtal::int32_t, double > ReSampler;
255 typedef DGtal::ConstImageAdapter<Image, Image::Domain, ReSampler,
256 Image::Value, DGtal::functors::Identity > SamplerImageAdapter;
257 typedef IntervalForegroundPredicate< SamplerImageAdapter > ImagePredicate;
258 typedef KanungoNoise< ImagePredicate, Z3i::Domain > KanungoPredicate;
259 typedef BinaryPointPredicate<DomainPredicate<Image::Domain>, KanungoPredicate, AndBoolFct2 > Predicate;
260 typedef Z3i::KSpace KSpace;
261 typedef KSpace::SCell SCell;
262 typedef KSpace::Cell Cell;
263
264 trace.beginBlock("Loading the file");
265
266 Image image = GenericReader<Image>::import( inputFileName );
267
268 PointVector<3,int> shiftVector3D( 0 ,0, 0 );
269 DGtal::functors::BasicDomainSubSampler< HyperRectDomain< SpaceND< 3, int > >,
270 DGtal::int32_t, double > reSampler(image.domain(),
271 aGridSizeReSample, shiftVector3D);
272 const functors::Identity identityFunctor{};
273 SamplerImageAdapter sampledImage ( image, reSampler.getSubSampledDomain(), reSampler, identityFunctor );
274 ImagePredicate predicateIMG = ImagePredicate( sampledImage, minImageThreshold, maxImageThreshold );
275 KanungoPredicate noisifiedPredicateIMG( predicateIMG, sampledImage.domain(), noiseLevel );
276 DomainPredicate<Z3i::Domain> domainPredicate( sampledImage.domain() );
277 AndBoolFct2 andF;
278 Predicate predicate(domainPredicate, noisifiedPredicateIMG, andF );
279
280
281 Z3i::Domain domain = sampledImage.domain();
282 Z3i::KSpace K;
283 bool space_ok = K.init( domain.lowerBound()-Z3i::Domain::Point::diagonal(),
284 domain.upperBound()+Z3i::Domain::Point::diagonal(), true );
285 if (!space_ok)
286 {
287 trace.error() << "Error in the Khalimsky space construction."<<std::endl;
288 return 2;
289 }
290 CanonicSCellEmbedder< KSpace > embedder( K );
291 SurfelAdjacency< Z3i::KSpace::dimension > Sadj( true );
292 trace.endBlock();
293 // Viewer settings
294
295
296 // Extraction of components
297 typedef KSpace::SurfelSet SurfelSet;
298 typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
299 typedef DigitalSurface< MySetOfSurfels > MyDigitalSurface;
300
301
302
303 trace.beginBlock("Extracting surfaces");
304 std::vector< std::vector<SCell > > vectConnectedSCell;
305 Surfaces<KSpace>::extractAllConnectedSCell(vectConnectedSCell,K, Sadj, predicate, false);
306 std::ofstream outDat;
307 if( enable_dat )
308 {
309 trace.info() << "Exporting curvature as dat file: "<< export_dat_filename <<std::endl;
310 outDat.open( export_dat_filename.c_str() );
311 outDat << "# data exported from 3dCurvatureViewer implementing the II curvature estimator (Coeurjolly, D.; Lachaud, J.O; Levallois, J., (2013). Integral based Curvature"
312 << " Estimators in Digital Geometry. DGCI 2013.) " << std::endl;
313 outDat << "# format: surfel coordinates (in Khalimsky space) curvature: "<< mode << std::endl;
314 }
315
316 trace.info()<<"Number of components= "<<vectConnectedSCell.size()<<std::endl;
317 trace.endBlock();
318
319 if( vectConnectedSCell.size() == 0 )
320 {
321 trace.error()<< "No surface component exists. Please check the vol file threshold parameter.";
322 trace.info()<<std::endl;
323 exit(2);
324 }
325 std::stringstream s;
326 s << "3dCurvatureViewerNoise - DGtalTools: ";
327 std::string name = inputFileName.substr(inputFileName.find_last_of("/")+1,inputFileName.size()) ;
328 s << " " << name << " (W key to display settings)" ;
329 polyscope::options::programName = s.str();
330 polyscope::options::buildGui=false;
331 polyscope::options::groundPlaneMode = polyscope::GroundPlaneMode::None;
332
333 typedef PolyscopeViewer<Z3i::Space, Z3i::KSpace> Viewer;
334 Viewer viewer( K );
335 viewer.allowReuseList = true;
336 polyscope::state::userCallback = myCallback;
337
338 unsigned int i = 0;
339 unsigned int max_size = 0;
340 for( unsigned int ii = 0; ii<vectConnectedSCell.size(); ++ii )
341 {
342 if( vectConnectedSCell[ii].size() <= threshold )
343 {
344 continue;
345 }
346 if( vectConnectedSCell[ii].size() > max_size )
347 {
348 max_size = vectConnectedSCell[ii].size();
349 i = ii;
350 }
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
462 for ( unsigned int i = 0; i < results.size(); ++i )
463 {
464 if( enable_visu )
465 {
466 viewer << WithQuantity(*abegin2, "curvature", results[i]);
467 }
468
469 if( enable_dat )
470 {
471 Point kCoords = K.uKCoords(K.unsigns(*abegin2));
472 outDat << kCoords[0] << " " << kCoords[1] << " " << kCoords[2] << " " << results[i] << std::endl;
473 }
474
475 ++abegin2;
476 }
477 }
478 else
479 {
480 typedef Z3i::Space::RealVector Quantity;
481 std::vector< Quantity > results;
482 std::back_insert_iterator< std::vector< Quantity > > resultsIterator( results );
483
484 if( mode.compare("prindir1") == 0 )
485 {
486 typedef functors::IIFirstPrincipalDirectionFunctor<Z3i::Space> MyIICurvatureFunctor;
487 typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
488
489 MyIICurvatureFunctor functor;
490 functor.init( h, re_convolution_kernel );
491
492 MyIIEstimator estimator( functor );
493 estimator.attach( K, predicate );
494 estimator.setParams( re_convolution_kernel/h );
495 estimator.init( h, abegin, aend );
496
497 estimator.eval( abegin, aend, resultsIterator );
498 }
499 else if( mode.compare("prindir2") == 0 )
500 {
501 typedef functors::IISecondPrincipalDirectionFunctor<Z3i::Space> MyIICurvatureFunctor;
502 typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
503
504 MyIICurvatureFunctor functor;
505 functor.init( h, re_convolution_kernel );
506
507 MyIIEstimator estimator( functor );
508 estimator.attach( K, predicate );
509 estimator.setParams( re_convolution_kernel/h );
510 estimator.init( h, abegin, aend );
511
512 estimator.eval( abegin, aend, resultsIterator );
513 }else if( mode.compare("normal") == 0 )
514 {
515 typedef functors::IINormalDirectionFunctor<Z3i::Space> MyIICurvatureFunctor;
516 typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
517
518 MyIICurvatureFunctor functor;
519 functor.init( h, re_convolution_kernel );
520
521 MyIIEstimator estimator( functor );
522 estimator.attach( K, predicate );
523 estimator.setParams( re_convolution_kernel/h );
524 estimator.init( h, abegin, aend );
525
526 estimator.eval( abegin, aend, resultsIterator );
527 }
528
529 //Visualization + export
530 for ( unsigned int i = 0; i < results.size(); ++i )
531 {
532 DGtal::Dimension kDim = K.sOrthDir( *abegin2 );
533 SCell outer = K.sIndirectIncident( *abegin2, kDim);
534 if ( predicate(Z3i::Point(embedder(outer), functors::Round<>()) ))
535 {
536 outer = K.sDirectIncident( *abegin2, kDim);
537 }
538
539 Cell unsignedSurfel = K.uCell( K.sKCoords(*abegin2) );
540
541 if( enable_visu )
542 {
543 viewer << DGtal::Color(255,255,255,255)
544 << unsignedSurfel;
545 }
546
547 if( enable_dat )
548 {
549 Point kCoords = K.uKCoords(K.unsigns(*abegin2));
550 outDat << kCoords[0] << " " << kCoords[1] << " " << kCoords[2] << " "
551 << results[i][0] << " " << results[i][1] << " " << results[i][2]
552 << std::endl;
553 }
554
555 RealPoint center = embedder( outer );
556
557 if( enable_visu )
558 {
559 if( mode.compare("prindir1") == 0 )
560 {
561 viewer.drawColor( AXIS_COLOR_BLUE );
562 }
563 else if( mode.compare("prindir2") == 0 )
564 {
565 viewer.drawColor( AXIS_COLOR_RED );
566 }
567 else if( mode.compare("normal") == 0 )
568 {
569 viewer.drawColor( AXIS_COLOR_GREEN );
570 }
571
572 viewer.drawLine (
573 RealPoint(
574 center[0] - 0.5 * results[i][0],
575 center[1] - 0.5 * results[i][1],
576 center[2] - 0.5 * results[i][2]
577 ),
578 RealPoint(
579 center[0] + 0.5 * results[i][0],
580 center[1] + 0.5 * results[i][1],
581 center[2] + 0.5 * results[i][2]
582 ));
583 }
584
585 ++abegin2;
586 }
587 }
588 trace.endBlock();
589
590 if( enable_dat )
591 {
592 outDat.close();
593 }
594
595 if( enable_visu )
596 {
597 viewer.show();
598 }
599
600 return 0;
601}
602
DGtal
Definition ATu0v1.h:57