DGtalTools 2.0.0
Loading...
Searching...
No Matches
3dCurvatureViewerNoise.cpp
1
30#include <iostream>
31#include "DGtal/base/Common.h"
32#include <cstring>
33
34#include "CLI11.hpp"
35
36
37// Shape constructors
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>
46
47#include "DGtal/images/ImageHelper.h"
48#include "DGtal/topology/DigitalSurface.h"
49#include "DGtal/graph/DepthFirstVisitor.h"
50#include "DGtal/graph/GraphVisitorRange.h"
51
52// Noise
53#include "DGtal/geometry/volumes/KanungoNoise.h"
54
55// Integral Invariant includes
56#include "DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h"
57#include "DGtal/geometry/surfaces/estimation/IntegralInvariantVolumeEstimator.h"
58#include "DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h"
59
60// Drawing
61#include "DGtal/io/viewers/PolyscopeViewer.h"
62
63using namespace DGtal;
64using namespace functors;
65
142const Color AXIS_COLOR_RED( 200, 20, 20, 255 );
143const Color AXIS_COLOR_GREEN( 20, 200, 20, 255 );
144const Color AXIS_COLOR_BLUE( 20, 20, 200, 255 );
145const double AXIS_LINESIZE = 0.05;
146
148
154void missingParam( std::string param )
155{
156 trace.error() << " Parameter: " << param << " is required.";
157 trace.info() << std::endl;
158}
159
160
161int main( int argc, char** argv )
162{
163 // parse command line using CLI ----------------------------------------------
164 CLI::App app;
165 std::string inputFileName;
166 double re_convolution_kernel;
167 double noiseLevel {0.5};
168 unsigned int threshold {8};
169 int minImageThreshold {0};
170 int maxImageThreshold {255};
171 std::string mode {"mean"};
172 std::string export_dat_filename;
173 bool exportOnly {false};
174 std::vector< double> vectScale;
175 bool normalization {false};
176
177
178 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");
179
180
181
182 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." )
183 ->required()
184 ->check(CLI::ExistingFile);
185
186 app.add_option("--radius,-r", re_convolution_kernel, "Kernel radius for IntegralInvariant" )
187 ->required();
188 app.add_option("--noise,-k", noiseLevel, "Level of Kanungo noise ]0;1[");
189 app.add_option("--threshold,-t", threshold, "Min size of SCell boundary of an object");
190 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 ] ).");
191 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] ).");
192 app.add_option("--mode,-m", mode, "type of output : mean, gaussian, k1, k2, prindir1, prindir2 or normal(default mean)")
193 -> check(CLI::IsMember({"mean","gaussian", "k1", "k2", "prindir1","prindir2", "normal" }));
194 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." );
195 app.add_flag("--exportOnly", exportOnly, "Used to only export the result without the 3d Visualisation (usefull for scripts).");
196
197 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.")
198 ->expected(3);
199
200 app.add_option("--normalization,-n",normalization, "When exporting to OBJ, performs a normalization so that the geometry fits in [-1/2,1/2]^3") ;
201
202 app.get_formatter()->column_width(40);
203 CLI11_PARSE(app, argc, argv);
204 // END parse command line using CLI ----------------------------------------------
205
206 bool neededArgsGiven=true;
207
208 if( noiseLevel < 0.0 || noiseLevel > 1.0 )
209 {
210 trace.error() << "The noise level should be in the interval: ]0, 1["<< std::endl;
211 exit(EXIT_FAILURE);
212 }
213 bool enable_visu = !exportOnly;
214 bool enable_dat = export_dat_filename != "";
215
216 if( !enable_visu && !enable_dat )
217 {
218 trace.error() << "You should specify what you want to export with --export and/or --exportDat, or remove --exportOnly." << std::endl;
219 neededArgsGiven = false;
220 }
221
222 double h = 1.0;
223
224 std::vector< double > aGridSizeReSample;
225 if( vectScale.size() == 3)
226 {
227 aGridSizeReSample.push_back(1.0/vectScale.at(0));
228 aGridSizeReSample.push_back(1.0/vectScale.at(1));
229 aGridSizeReSample.push_back(1.0/vectScale.at(2));
230
231 }
232 else
233 {
234 aGridSizeReSample.push_back(1.0);
235 aGridSizeReSample.push_back(1.0);
236 aGridSizeReSample.push_back(1.0);
237 }
238
239
240
241 // Construction of the shape from vol file
242 typedef Z3i::Space::RealPoint RealPoint;
243 typedef Z3i::Point Point;
244 typedef ImageSelector< Z3i::Domain, int>::Type Image;
245 typedef DGtal::functors::BasicDomainSubSampler< HyperRectDomain<SpaceND<3, int> >,
246 DGtal::int32_t, double > ReSampler;
247 typedef DGtal::ConstImageAdapter<Image, Image::Domain, ReSampler,
248 Image::Value, DGtal::functors::Identity > SamplerImageAdapter;
249 typedef IntervalForegroundPredicate< SamplerImageAdapter > ImagePredicate;
250 typedef KanungoNoise< ImagePredicate, Z3i::Domain > KanungoPredicate;
251 typedef BinaryPointPredicate<DomainPredicate<Image::Domain>, KanungoPredicate, AndBoolFct2 > Predicate;
252 typedef Z3i::KSpace KSpace;
253 typedef KSpace::SCell SCell;
254 typedef KSpace::Cell Cell;
255
256 trace.beginBlock("Loading the file");
257
258 Image image = GenericReader<Image>::import( inputFileName );
259
260 PointVector<3,int> shiftVector3D( 0 ,0, 0 );
261 DGtal::functors::BasicDomainSubSampler< HyperRectDomain< SpaceND< 3, int > >,
262 DGtal::int32_t, double > reSampler(image.domain(),
263 aGridSizeReSample, shiftVector3D);
264 const functors::Identity identityFunctor{};
265 SamplerImageAdapter sampledImage ( image, reSampler.getSubSampledDomain(), reSampler, identityFunctor );
266 ImagePredicate predicateIMG = ImagePredicate( sampledImage, minImageThreshold, maxImageThreshold );
267 KanungoPredicate noisifiedPredicateIMG( predicateIMG, sampledImage.domain(), noiseLevel );
268 DomainPredicate<Z3i::Domain> domainPredicate( sampledImage.domain() );
269 AndBoolFct2 andF;
270 Predicate predicate(domainPredicate, noisifiedPredicateIMG, andF );
271
272
273 Z3i::Domain domain = sampledImage.domain();
274 Z3i::KSpace K;
275 bool space_ok = K.init( domain.lowerBound()-Z3i::Domain::Point::diagonal(),
276 domain.upperBound()+Z3i::Domain::Point::diagonal(), true );
277 if (!space_ok)
278 {
279 trace.error() << "Error in the Khalimsky space construction."<<std::endl;
280 return 2;
281 }
282 CanonicSCellEmbedder< KSpace > embedder( K );
283 SurfelAdjacency< Z3i::KSpace::dimension > Sadj( true );
284 trace.endBlock();
285 // Viewer settings
286
287
288 // Extraction of components
289 typedef KSpace::SurfelSet SurfelSet;
290 typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
291 typedef DigitalSurface< MySetOfSurfels > MyDigitalSurface;
292
293
294
295 trace.beginBlock("Extracting surfaces");
296 std::vector< std::vector<SCell > > vectConnectedSCell;
297 Surfaces<KSpace>::extractAllConnectedSCell(vectConnectedSCell,K, Sadj, predicate, false);
298 std::ofstream outDat;
299 if( enable_dat )
300 {
301 trace.info() << "Exporting curvature as dat file: "<< export_dat_filename <<std::endl;
302 outDat.open( export_dat_filename.c_str() );
303 outDat << "# data exported from 3dCurvatureViewer implementing the II curvature estimator (Coeurjolly, D.; Lachaud, J.O; Levallois, J., (2013). Integral based Curvature"
304 << " Estimators in Digital Geometry. DGCI 2013.) " << std::endl;
305 outDat << "# format: surfel coordinates (in Khalimsky space) curvature: "<< mode << std::endl;
306 }
307
308 trace.info()<<"Number of components= "<<vectConnectedSCell.size()<<std::endl;
309 trace.endBlock();
310
311 if( vectConnectedSCell.size() == 0 )
312 {
313 trace.error()<< "No surface component exists. Please check the vol file threshold parameter.";
314 trace.info()<<std::endl;
315 exit(2);
316 }
317
318 typedef PolyscopeViewer<Z3i::Space, Z3i::KSpace> Viewer;
319 Viewer viewer( K );
320 viewer.allowReuseList = true;
321
322 unsigned int i = 0;
323 unsigned int max_size = 0;
324 for( unsigned int ii = 0; ii<vectConnectedSCell.size(); ++ii )
325 {
326 if( vectConnectedSCell[ii].size() <= threshold )
327 {
328 continue;
329 }
330 if( vectConnectedSCell[ii].size() > max_size )
331 {
332 max_size = vectConnectedSCell[ii].size();
333 i = ii;
334 }
335 }
336
337 MySetOfSurfels aSet(K, Sadj);
338
339 for( std::vector<SCell>::const_iterator it = vectConnectedSCell.at(i).begin();
340 it != vectConnectedSCell.at(i).end();
341 ++it )
342 {
343 aSet.surfelSet().insert( *it);
344 }
345
346 MyDigitalSurface digSurf( aSet );
347
348
349 typedef DepthFirstVisitor<MyDigitalSurface> Visitor;
350 typedef GraphVisitorRange< Visitor > VisitorRange;
351 typedef VisitorRange::ConstIterator SurfelConstIterator;
352 VisitorRange range( new Visitor( digSurf, *digSurf.begin() ) );
353 SurfelConstIterator abegin = range.begin();
354 SurfelConstIterator aend = range.end();
355
356 VisitorRange range2( new Visitor( digSurf, *digSurf.begin() ) );
357 SurfelConstIterator abegin2 = range2.begin();
358
359 trace.beginBlock("Curvature computation on a component");
360 if( ( mode.compare("gaussian") == 0 ) || ( mode.compare("mean") == 0 )
361 || ( mode.compare("k1") == 0 ) || ( mode.compare("k2") == 0 ))
362 {
363 typedef double Quantity;
364 std::vector< Quantity > results;
365 std::back_insert_iterator< std::vector< Quantity > > resultsIterator( results );
366 if ( mode.compare("mean") == 0 )
367 {
368 typedef functors::IIMeanCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
369 typedef IntegralInvariantVolumeEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
370
371 MyIICurvatureFunctor functor;
372 functor.init( h, re_convolution_kernel );
373
374 MyIIEstimator estimator( functor );
375 estimator.attach( K, predicate );
376 estimator.setParams( re_convolution_kernel/h );
377 estimator.init( h, abegin, aend );
378
379 estimator.eval( abegin, aend, resultsIterator );
380 }
381 else if ( mode.compare("gaussian") == 0 )
382 {
383 typedef functors::IIGaussianCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
384 typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
385
386 MyIICurvatureFunctor functor;
387 functor.init( h, re_convolution_kernel );
388
389 MyIIEstimator estimator( functor ); estimator.attach( K,
390 predicate ); estimator.setParams( re_convolution_kernel/h );
391 estimator.init( h, abegin, aend );
392
393 estimator.eval( abegin, aend, resultsIterator );
394 }
395 else if ( mode.compare("k1") == 0 )
396 {
397 typedef functors::IIFirstPrincipalCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
398 typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
399
400 MyIICurvatureFunctor functor;
401 functor.init( h, re_convolution_kernel );
402
403 MyIIEstimator estimator( functor );
404 estimator.attach( K, predicate );
405 estimator.setParams( re_convolution_kernel/h );
406 estimator.init( h, abegin, aend );
407
408 estimator.eval( abegin, aend, resultsIterator );
409 }
410 else if ( mode.compare("k2") == 0 )
411 {
412 typedef functors::IISecondPrincipalCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
413 typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
414
415 MyIICurvatureFunctor functor;
416 functor.init( h, re_convolution_kernel );
417
418 MyIIEstimator estimator( functor );
419 estimator.attach( K, predicate );
420 estimator.setParams( re_convolution_kernel/h );
421 estimator.init( h, abegin, aend );
422
423 estimator.eval( abegin, aend, resultsIterator );
424 }
425 trace.endBlock();
426
427
428 // Drawing results
429 trace.beginBlock("Visualisation");
430 Quantity min = results[ 0 ];
431 Quantity max = results[ 0 ];
432 for ( unsigned int i = 1; i < results.size(); ++i )
433 {
434 if ( results[ i ] < min )
435 {
436 min = results[ i ];
437 }
438 else if ( results[ i ] > max )
439 {
440 max = results[ i ];
441 }
442 }
443 trace.info() << "Max value= "<<max<<" min value= "<<min<<std::endl;
444 ASSERT( min <= max );
445
446 for ( unsigned int i = 0; i < results.size(); ++i )
447 {
448 if( enable_visu )
449 {
450 viewer << WithQuantity(*abegin2, "curvature", results[i]);
451 }
452
453 if( enable_dat )
454 {
455 Point kCoords = K.uKCoords(K.unsigns(*abegin2));
456 outDat << kCoords[0] << " " << kCoords[1] << " " << kCoords[2] << " " << results[i] << std::endl;
457 }
458
459 ++abegin2;
460 }
461 }
462 else
463 {
464 typedef Z3i::Space::RealVector Quantity;
465 std::vector< Quantity > results;
466 std::back_insert_iterator< std::vector< Quantity > > resultsIterator( results );
467
468 if( mode.compare("prindir1") == 0 )
469 {
470 typedef functors::IIFirstPrincipalDirectionFunctor<Z3i::Space> MyIICurvatureFunctor;
471 typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
472
473 MyIICurvatureFunctor functor;
474 functor.init( h, re_convolution_kernel );
475
476 MyIIEstimator estimator( functor );
477 estimator.attach( K, predicate );
478 estimator.setParams( re_convolution_kernel/h );
479 estimator.init( h, abegin, aend );
480
481 estimator.eval( abegin, aend, resultsIterator );
482 }
483 else if( mode.compare("prindir2") == 0 )
484 {
485 typedef functors::IISecondPrincipalDirectionFunctor<Z3i::Space> MyIICurvatureFunctor;
486 typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
487
488 MyIICurvatureFunctor functor;
489 functor.init( h, re_convolution_kernel );
490
491 MyIIEstimator estimator( functor );
492 estimator.attach( K, predicate );
493 estimator.setParams( re_convolution_kernel/h );
494 estimator.init( h, abegin, aend );
495
496 estimator.eval( abegin, aend, resultsIterator );
497 }else if( mode.compare("normal") == 0 )
498 {
499 typedef functors::IINormalDirectionFunctor<Z3i::Space> MyIICurvatureFunctor;
500 typedef IntegralInvariantCovarianceEstimator<Z3i::KSpace, Predicate, MyIICurvatureFunctor> MyIIEstimator;
501
502 MyIICurvatureFunctor functor;
503 functor.init( h, re_convolution_kernel );
504
505 MyIIEstimator estimator( functor );
506 estimator.attach( K, predicate );
507 estimator.setParams( re_convolution_kernel/h );
508 estimator.init( h, abegin, aend );
509
510 estimator.eval( abegin, aend, resultsIterator );
511 }
512
513 //Visualization + export
514 for ( unsigned int i = 0; i < results.size(); ++i )
515 {
516 DGtal::Dimension kDim = K.sOrthDir( *abegin2 );
517 SCell outer = K.sIndirectIncident( *abegin2, kDim);
518 if ( predicate(Z3i::Point(embedder(outer), functors::Round<>()) ))
519 {
520 outer = K.sDirectIncident( *abegin2, kDim);
521 }
522
523 Cell unsignedSurfel = K.uCell( K.sKCoords(*abegin2) );
524
525 if( enable_visu )
526 {
527 viewer << DGtal::Color(255,255,255,255)
528 << unsignedSurfel;
529 }
530
531 if( enable_dat )
532 {
533 Point kCoords = K.uKCoords(K.unsigns(*abegin2));
534 outDat << kCoords[0] << " " << kCoords[1] << " " << kCoords[2] << " "
535 << results[i][0] << " " << results[i][1] << " " << results[i][2]
536 << std::endl;
537 }
538
539 RealPoint center = embedder( outer );
540
541 if( enable_visu )
542 {
543 if( mode.compare("prindir1") == 0 )
544 {
545 viewer.drawColor( AXIS_COLOR_BLUE );
546 }
547 else if( mode.compare("prindir2") == 0 )
548 {
549 viewer.drawColor( AXIS_COLOR_RED );
550 }
551 else if( mode.compare("normal") == 0 )
552 {
553 viewer.drawColor( AXIS_COLOR_GREEN );
554 }
555
556 viewer.drawLine (
557 RealPoint(
558 center[0] - 0.5 * results[i][0],
559 center[1] - 0.5 * results[i][1],
560 center[2] - 0.5 * results[i][2]
561 ),
562 RealPoint(
563 center[0] + 0.5 * results[i][0],
564 center[1] + 0.5 * results[i][1],
565 center[2] + 0.5 * results[i][2]
566 ));
567 }
568
569 ++abegin2;
570 }
571 }
572 trace.endBlock();
573
574 if( enable_dat )
575 {
576 outDat.close();
577 }
578
579 if( enable_visu )
580 {
581 viewer.show();
582 }
583
584 return 0;
585}
586
DGtal
Definition ATu0v1.h:57