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