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