50 #include "DGtal/base/Common.h"
51 #include "DGtal/base/Exceptions.h"
52 #include "DGtal/kernel/SpaceND.h"
53 #include "DGtal/kernel/domains/HyperRectDomain.h"
54 #include "DGtal/kernel/BasicPointPredicates.h"
55 #include "DGtal/math/linalg/EigenDecomposition.h"
56 #include "DGtal/topology/KhalimskySpaceND.h"
57 #include "DGtal/geometry/curves/GridCurve.h"
58 #include "DGtal/geometry/curves/Naive3DDSSComputer.h"
59 #include "DGtal/geometry/curves/StandardDSS6Computer.h"
60 #include "DGtal/geometry/curves/SaturatedSegmentation.h"
61 #include "DGtal/geometry/volumes/distance/ExactPredicateLpSeparableMetric.h"
62 #include "DGtal/geometry/volumes/estimation/VoronoiCovarianceMeasure.h"
63 #include "DGtal/geometry/curves/estimation/LambdaMST3D.h"
64 #include "DGtal/geometry/curves/estimation/FunctorsLambdaMST.h"
65 #include "DGtal/io/viewers/Viewer3D.h"
66 #include "DGtal/io/readers/PointListReader.h"
67 #include "DGtal/io/DrawWithDisplay3DModifier.h"
70 using namespace DGtal;
137 const Color AXIS_COLOR( 0, 0, 0, 255 );
138 const double AXIS_LINESIZE = 0.1;
139 const Color XY_COLOR( 0, 0, 255, 50 );
140 const Color XZ_COLOR( 0, 255, 0, 50 );
141 const Color YZ_COLOR( 255, 0, 0, 50 );
142 const Color CURVE3D_COLOR( 100, 100, 140, 128 );
143 const Color CURVE2D_COLOR( 200, 200, 200, 100 );
144 const double MS3D_LINESIZE = 0.05;
149 template <
typename Po
int,
typename RealPo
int,
typename space,
typename kspace >
151 const Point & lowerBound,
const Point & upperBound,
152 const std::string & mode )
154 RealPoint p0( (
double)lowerBound[ 0 ]-0.5,
155 (
double)lowerBound[ 1 ]-0.5,
156 (
double)lowerBound[ 2 ]-0.5 );
157 RealPoint p1( (
double)upperBound[ 0 ]-0.5,
158 (
double)upperBound[ 1 ]-0.5,
159 (
double)upperBound[ 2 ]-0.5 );
160 if ( ( mode ==
"WIRED" ) || ( mode ==
"COLORED" ) )
188 if ( mode ==
"COLORED" )
208 template <
typename KSpace,
typename Naive3DDSSComputer,
typename space,
typename kspace >
216 template <
typename Po
int1,
typename Po
int2>
217 void assign( Point1 & p1,
const Point2 & p2 )
224 template <
typename KSpace,
typename Naive3DDSSComputer,
typename space,
typename kspace >
234 PointR3d interceptR, thicknessR;
245 assign( direction, directionZ3 );
246 direction /= direction.
norm();
247 assign( P1, *dss3d.
begin() );
248 assign( P2, *(dss3d.
end()-1) );
249 double t1 = (P1 - intercept).dot( direction );
250 double t2 = (P2 - intercept).dot( direction );
260 template <
typename KSpace,
typename StandardDSS6Computer,
typename space,
typename kspace >
266 typedef typename ArithmeticalDSSComputer2d::ConstIterator ConstIterator2d;
267 typedef typename ArithmeticalDSSComputer2d::Point Point2d;
274 for ( ConstIterator2d itP = dss2d.begin(), itPEnd = dss2d.end(); itP != itPEnd; ++itP )
280 case 0: q = Point3d( 2*b[ i ] , 2*p[ 0 ]+1, 2*p[ 1 ]+1 );
break;
281 case 1: q = Point3d( 2*p[ 0 ]+1, 2*b[ i ] , 2*p[ 1 ]+1 );
break;
282 case 2: q = Point3d( 2*p[ 0 ]+1, 2*p[ 1 ]+1, 2*b[ i ] );
break;
290 template <
typename KSpace,
typename StandardDSS6Computer,
typename space,
typename kspace >
297 typedef typename ArithmeticalDSSComputer2d::ConstIterator ConstIterator2d;
298 typedef typename ArithmeticalDSSComputer2d::Point Point2d;
306 const typename ArithmeticalDSSComputer2d::Primitive & dss2d
309 std::vector<PointD2d> pts2d;
310 pts2d.push_back( dss2d.project(dss2d.back(), dss2d.Uf()) );
311 pts2d.push_back( dss2d.project(dss2d.back(), dss2d.Lf()) );
312 pts2d.push_back( dss2d.project(dss2d.front(), dss2d.Lf()) );
313 pts2d.push_back( dss2d.project(dss2d.front(), dss2d.Uf()) );
314 std::vector<PointD3D> bb;
316 for (
unsigned int j = 0; j < pts2d.size(); ++j )
320 case 0: p3.center[0] = (double) b[ i ]-0.5; p3.center[1] = pts2d[ j ][ 0 ]; p3.center[2] = pts2d[ j ][ 1 ];
break;
321 case 1: p3.center[0] = pts2d[ j ][ 0 ]; p3.center[1] = (double) b[ i ]-0.5; p3.center[2] = pts2d[ j ][ 1 ];
break;
322 case 2: p3.center[0] = pts2d[ j ][ 0 ]; p3.center[1] = pts2d[ j ][ 1 ]; p3.center[2] = (double) b[ i ]-0.5;
break;
326 for (
unsigned int j = 0; j < pts2d.size(); ++j ){
337 template <
typename KSpace,
typename Naive3DDSSComputer,
typename space,
typename kspace >
343 typedef typename ArithmeticalDSSComputer2d::ConstIterator ConstIterator2d;
344 typedef typename ArithmeticalDSSComputer2d::Point Point2d;
357 if ( validXY && validXZ )
359 for ( ConstIterator2d itXY = dssXY.begin(), itXZ = dssXZ.begin(), itPEnd = dssXY.end(); itXY != itPEnd; ++itXY, ++itXZ )
361 Point2d p1 = *itXY, p2 = *itXZ;
363 q1 = Point3d ( 2*b[ 0 ], 2*p1[ 1 ]+1, 2*p2[ 1 ]+1 );
364 q2 = Point3d ( 2*p2[ 0 ]+1, 2*b[ 1 ], 2*p2[ 1 ]+1 );
365 q3 = Point3d ( 2*p1[ 0 ]+1, 2*p1[ 1 ]+1, 2*b[ 2 ] );
374 if ( validYZ && validXY )
376 for ( ConstIterator2d itYZ = dssYZ.begin(), itXY = dssXY.begin(), itPEnd = dssXY.end(); itXY != itPEnd; ++itXY, ++itYZ )
378 Point2d p1 = *itYZ, p2 = *itXY;
380 q1 = Point3d ( 2*b[ 0 ], 2*p1[ 0 ]+1, 2*p1[ 1 ]+1 );
381 q2 = Point3d ( 2*p2[ 0 ]+1, 2*b[ 1 ], 2*p1[ 1 ]+1 );
382 q3 = Point3d ( 2*p2[ 0 ]+1, 2*p2[ 1 ]+1, 2*b[ 2 ] );
391 for ( ConstIterator2d itYZ = dssYZ.begin(), itXZ = dssXZ.begin(), itPEnd = dssXZ.end(); itXZ != itPEnd; ++itXZ, ++itYZ )
393 Point2d p1 = *itYZ, p2 = *itXZ;
395 q1 = Point3d ( 2*b[ 0 ], 2*p1[ 0 ]+1, 2*p1[ 1 ]+1 );
396 q2 = Point3d ( 2*p2[ 0 ]+1, 2*b[ 1 ], 2*p2[ 1 ]+1 );
397 q3 = Point3d ( 2*p2[ 0 ]+1, 2*p1[ 0 ]+1, 2*b[ 2 ] );
408 template <
typename KSpace,
typename Naive3DDSSComputer,
typename space,
typename kspace >
415 typedef typename ArithmeticalDSSComputer2d::ConstIterator ConstIterator2d;
416 typedef typename ArithmeticalDSSComputer2d::Point Point2d;
427 const typename ArithmeticalDSSComputer2d::Primitive & dss2d
430 std::vector<PointD2d> pts2d;
431 pts2d.push_back( dss2d.project(dss2d.back(), dss2d.Uf()) );
432 pts2d.push_back( dss2d.project(dss2d.back(), dss2d.Lf()) );
433 pts2d.push_back( dss2d.project(dss2d.front(), dss2d.Lf()) );
434 pts2d.push_back( dss2d.project(dss2d.front(), dss2d.Uf()) );
435 std::vector<PointD3D> bb;
437 for (
unsigned int j = 0; j < pts2d.size(); ++j )
441 case 0: p3.center[0] = (double) b[ i ]-0.5; p3.center[1] = pts2d[ j ][ 0 ]; p3.center[2] = pts2d[ j ][ 1 ];
break;
442 case 1: p3.center[0] = pts2d[ j ][ 0 ]; p3.center[1] = (double) b[ i ]-0.5; p3.center[2] = pts2d[ j ][ 1 ];
break;
443 case 2: p3.center[0] = pts2d[ j ][ 0 ]; p3.center[1] = pts2d[ j ][ 1 ]; p3.center[2] = (double) b[ i ]-0.5;
break;
447 for (
unsigned int j = 0; j < pts2d.size(); ++j ){
459 template <
typename KSpace,
typename Po
intIterator,
typename space,
typename kspace >
461 const KSpace & ks, PointIterator b, PointIterator e,
462 bool dss3d,
bool proj2d,
bool dss2d,
bool tangent,
465 typedef typename PointIterator::value_type
Point;
468 typedef typename Decomposition::SegmentComputerIterator SegmentComputerIterator;
469 typedef typename SegmentComputer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
470 SegmentComputer algo;
471 Decomposition theDecomposition(b, e, algo);
473 viewer <<
SetMode3D( algo.className(),
"BoundingBox" );
477 for ( SegmentComputerIterator i = theDecomposition.begin();
478 i != theDecomposition.end(); ++i)
480 SegmentComputer ms3d(*i);
481 const ArithmeticalDSSComputer2d & dssXY = ms3d.arithmeticalDSS2dXY();
482 const ArithmeticalDSSComputer2d & dssXZ = ms3d.arithmeticalDSS2dXZ();
483 const ArithmeticalDSSComputer2d & dssYZ = ms3d.arithmeticalDSS2dYZ();
484 Point f = *ms3d.begin();
485 Point l = *(ms3d.end() - 1);
488 <<
" [" << f[ 0 ] <<
"," << f[ 1 ] <<
","<< f[ 2 ] <<
"]"
489 <<
"->[" << l[ 0 ] <<
"," << l[ 1 ] <<
","<< l[ 2 ] <<
"]"
491 << dssXY.a() <<
"," << dssXY.b() <<
"," << dssXY.mu()
493 << dssXZ.a() <<
"," << dssXZ.b() <<
"," << dssXZ.mu()
495 << dssYZ.a() <<
"," << dssYZ.b() <<
"," << dssYZ.mu()
497 Color color = cmap_hue( c );
498 if ( tangent ) displayDSS3dTangent( viewer, ks, ms3d, color );
499 if ( dss3d ) displayDSS3d( viewer, ks, ms3d, color );
500 if ( dss2d ) displayDSS2d6( viewer, ks, ms3d, color );
501 if ( proj2d ) displayProj2d6( viewer, ks, ms3d, CURVE2D_COLOR );
511 template <
typename KSpace,
typename Po
intIterator,
typename space,
typename kspace >
513 const KSpace & ks, PointIterator b, PointIterator e,
514 bool dss3d,
bool proj2d,
bool dss2d,
bool tangent,
517 typedef typename PointIterator::value_type
Point;
520 typedef typename Decomposition::SegmentComputerIterator SegmentComputerIterator;
521 typedef typename SegmentComputer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
522 SegmentComputer algo;
523 Decomposition theDecomposition(b, e, algo);
525 viewer <<
SetMode3D( algo.className(),
"BoundingBox" );
529 for ( SegmentComputerIterator i = theDecomposition.begin();
530 i != theDecomposition.end(); ++i)
532 SegmentComputer ms3d(*i);
533 const ArithmeticalDSSComputer2d & dssXY = ms3d.arithmeticalDSS2dXY();
534 const ArithmeticalDSSComputer2d & dssXZ = ms3d.arithmeticalDSS2dXZ();
535 const ArithmeticalDSSComputer2d & dssYZ = ms3d.arithmeticalDSS2dYZ();
536 Point f = *ms3d.begin();
537 Point l = *(ms3d.end() - 1);
540 <<
" [" << f[ 0 ] <<
"," << f[ 1 ] <<
","<< f[ 2 ] <<
"]"
541 <<
"->[" << l[ 0 ] <<
"," << l[ 1 ] <<
","<< l[ 2 ] <<
"]"
543 << dssXY.a() <<
"," << dssXY.b() <<
"," << dssXY.mu()
545 << dssXZ.a() <<
"," << dssXZ.b() <<
"," << dssXZ.mu()
547 << dssYZ.a() <<
"," << dssYZ.b() <<
"," << dssYZ.mu()
551 Color color = cmap_hue( c );
552 if ( tangent ) displayDSS3dTangent( viewer, ks, ms3d, color );
553 if ( dss3d ) displayDSS3d( viewer, ks, ms3d, color );
554 if ( dss2d ) displayDSS2d26( viewer, ks, ms3d, color );
555 if ( proj2d ) displayProj2d26( viewer, ks, ms3d, CURVE2D_COLOR );
561 template <
typename Po
intIterator,
typename Space,
typename TangentSequence >
562 void ComputeVCM (
const double & R,
const double & r,
563 const PointIterator & begin,
const PointIterator & end, TangentSequence & tangents,
const std::string & output )
565 using namespace DGtal;
571 typedef LinearAlgebraTool::Matrix Matrix;
573 fstream outputStream;
574 outputStream.open ( ( output +
".vcm" ).c_str(), std::ios::out );
575 outputStream <<
"# VCM estimation R=" << R <<
" r=" << r <<
" chi=hat" << endl;
578 VCM vcm ( R, ceil( r ), l2,
true );
579 vcm.init( begin, end );
580 KernelFunction chi( 1.0, r );
582 typename Space::RealVector eval;
583 for ( PointIterator it = begin; it != end; ++it )
586 vcm_r = vcm.measure( chi, *it );
587 LinearAlgebraTool::getEigenDecomposition ( vcm_r, evec, eval );
588 typename Space::RealVector tangent = evec.column( 0 );
589 tangents.push_back ( tangent );
590 outputStream << (*it)[0] <<
" " << (*it)[1] <<
" " << (*it)[2] <<
" "
591 << tangent[0] <<
" " << tangent[1] <<
" " << tangent[2] << endl;
593 outputStream.close();
596 template <
typename Po
intIterator,
typename Space,
typename TangentSequence >
597 void ComputeLMST6 (
const PointIterator & begin,
const PointIterator & end, TangentSequence & tangents,
const std::string & output )
599 typedef typename PointIterator::value_type
Point;
602 typedef typename Decomposition::SegmentComputerIterator SegmentComputerIterator;
603 typedef typename SegmentComputer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
604 SegmentComputer algo;
605 Decomposition theDecomposition ( begin, end, algo );
607 fstream outputStream;
608 outputStream.open ( ( output +
".lmst" ).c_str(), std::ios::out );
609 outputStream <<
"X Y Z X Y Z" << endl;
612 lmst64.
attach ( theDecomposition );
613 lmst64.
init ( begin, end );
614 lmst64.
eval ( begin, end, std::back_inserter ( tangents ) );
615 typename TangentSequence::iterator itt = tangents.
begin();
616 for ( PointIterator it = begin; it != end; ++it, ++itt )
618 typename Space::RealVector & tangent = (*itt);
619 if ( tangent.norm() != 0. )
620 tangent = tangent.getNormalized();
621 outputStream << (*it)[0] <<
" " << (*it)[1] <<
" " << (*it)[2] <<
" "
622 << tangent[0] <<
" " << tangent[1] <<
" " << tangent[2] << endl;
626 template <
typename Po
intIterator,
typename Space,
typename TangentSequence >
627 void ComputeLMST26 (
const PointIterator & begin,
const PointIterator & end, TangentSequence & tangents,
const std::string & output )
629 typedef typename PointIterator::value_type
Point;
632 typedef typename Decomposition::SegmentComputerIterator SegmentComputerIterator;
633 typedef typename SegmentComputer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
634 SegmentComputer algo;
635 Decomposition theDecomposition ( begin, end, algo );
637 fstream outputStream;
638 outputStream.open ( ( output +
".lmst" ).c_str(), std::ios::out );
639 outputStream <<
"X Y Z X Y Z" << endl;
642 lmst64.
attach ( theDecomposition );
643 lmst64.
init ( begin, end );
644 lmst64.
eval ( begin, end, std::back_inserter ( tangents ) );
645 typename TangentSequence::iterator itt = tangents.
begin();
646 for ( PointIterator it = begin; it != end; ++it, ++itt )
648 typename Space::RealVector & tangent = (*itt);
649 if ( tangent.norm() != 0. )
650 tangent = tangent.getNormalized();
651 outputStream << (*it)[0] <<
" " << (*it)[1] <<
" " << (*it)[2] <<
" "
652 << tangent[0] <<
" " << tangent[1] <<
" " << tangent[2] << endl;
656 template <
typename Po
intIterator,
typename Space,
int CONNECT = 8 >
657 void find_main_axes ( PointIterator begin, PointIterator end, multimap < typename Space::Point, string > & axes )
662 Segmentation segmenter ( begin, end, SegmentComputer ( ) );
664 for (
auto it = begin; it != end; ++it )
666 unsigned int cX = 0, cY = 0, cZ = 0;
667 for (
typename Segmentation::SegmentComputerIterator idss = segmenter.begin ( ); idss != segmenter.end ( ); ++idss )
669 if ( idss->isInDSS ( *it ) )
671 const typename SegmentComputer::ArithmeticalDSSComputer2d & dssXY = (*idss).arithmeticalDSS2dXY();
672 const typename SegmentComputer::ArithmeticalDSSComputer2d & dssXZ = (*idss).arithmeticalDSS2dXZ();
673 const typename SegmentComputer::ArithmeticalDSSComputer2d & dssYZ = (*idss).arithmeticalDSS2dYZ();
674 unsigned int lenXY = distance ( dssXY.begin(), dssXY.end() );
675 unsigned int lenXZ = distance ( dssXZ.begin(), dssXZ.end() );
676 unsigned int lenYZ = distance ( dssYZ.begin(), dssYZ.end() );
678 if ( lenXY >= lenYZ && lenXZ >= lenYZ )
680 else if ( lenXY >= lenXZ && lenYZ >= lenXZ )
687 axes.insert ( pair <typename Space::Point, string > ( *it,
string (
"X" ) ) );
689 axes.insert ( pair <typename Space::Point, string > ( *it,
string (
"Y" ) ) );
691 axes.insert ( pair <typename Space::Point, string > ( *it,
string (
"Z" ) ) );
696 template <
typename Po
intIterator,
typename Space >
697 void print_main_axes ( PointIterator begin, PointIterator end, multimap < typename Space::Point, string > & axes )
699 for ( PointIterator itt = begin; itt != end; ++itt )
701 typename std::multimap<typename Space::Point, string>::const_iterator it = axes.lower_bound(*itt);
702 typename std::multimap<typename Space::Point, string>::const_iterator it2 = axes.upper_bound(*itt);
703 cout <<
"(" << it->first[0] <<
", " << it->first[1] <<
", " << it->first[2] <<
"); MAIN_AXES = (";
704 for (; it != it2; it++ )
706 if ( distance( it, it2 ) > 1 )
707 cout << it->second <<
", ";
709 cout << it->second <<
")";
723 int main(
int argc,
char **argv)
726 using namespace DGtal;
729 typedef Z3::Point
Point;
733 typedef typename vector<Point>::const_iterator PointIterator;
736 QApplication application(argc,argv);
740 std::string inputFileName;
741 std::string viewFlag {
"OFF"};
742 std::string viewBox {
"WIRED"};
743 std::string connectivity {
"6"};
744 std::string method {
"VCM"};
745 std::string axesFlag {
"OFF"};
746 std::string outputFileName {
"3d-curve-tangent-estimations"};
747 bool curve3d {
false};
748 bool curve2d {
false};
749 bool cover3d {
false};
750 bool cover2d {
false};
751 bool displayTangent {
false};
752 double bigRad {10.0};
753 double smallRad {3.0};
755 unsigned int nbColors {3};
756 app.description(
"This program estimates the tangent vector to a set of 3D integer points, which are supposed to approximate a 3D curve. This set of points is given as a list of points in file <input>.\n The tangent estimator uses either the digital Voronoi Covariance Measure (VCM) or the 3D lambda-Maximal Segment Tangent (L-MST).\n This program can also displays the curve and tangent estimations, and it can also extract maximal digital straight segments (2D and 3D).\n Note: It is not compulsory for the points to be ordered in sequence, except if you wish to compute maximal digital straight segments. In this case, you can select the connectivity of your curve between 6 (standard) or 26 (naive).\n Example:\n 3dCurveTangentEstimator -i ${DGtal}/examples/samples/sinus.dat -V ON -c -R 20 -r 3 -T 6\n" );
757 app.add_option(
"-i,--input,1", inputFileName,
"the name of the text file containing the list of 3D points: (x y z) per line." )
759 ->check(CLI::ExistingFile);
761 app.add_option(
"--view,-V", viewFlag,
"toggles display ON/OFF",
true );
762 app.add_option(
"--box,-b", box,
"specifies the tightness of the bounding box around the curve with a given integer displacement <arg> to enlarge it (0 is tight)",
true);
763 app.add_option(
"--viewBox,-v", viewBox,
"displays the bounding box, <arg>=WIRED means that only edges are displayed, <arg>=COLORED adds colors for planes (XY is red, XZ green, YZ, blue).",
true)
764 -> check(CLI::IsMember({
"WIRED" ,
"COLORED"}));
765 app.add_option(
"--connectivity,-T",connectivity,
"specifies whether it is a 6-connected curve or a 26-connected curve: arg=6 | 26.",
true )
766 -> check(CLI::IsMember({
"6",
"26"}));
767 app.add_flag(
"--curve3d,-C", curve3d,
"displays the 3D curve" );
768 app.add_flag(
"--curve2d,-c", curve2d,
"displays the 2D projections of the 3D curve on the bounding box" );
769 app.add_flag(
"--cover3d,-3", cover3d,
"displays the 3D tangential cover of the curve");
770 app.add_flag(
"--cover2d,-2", cover2d,
"displays the 2D projections of the 3D tangential cover of the curve");
771 app.add_flag(
"--tangent,-t", displayTangent,
"displays the tangents to the curve.");
772 app.add_option(
"--nbColors,-n", nbColors,
"sets the number of successive colors used for displaying 2d and 3d maximal segments (default is 3: red, green, blue)",
true);
774 app.add_option(
"--big-radius,-R",bigRad,
"the radius parameter R in the VCM estimator.",
true);
775 app.add_option(
"--small-radius,-r",smallRad,
"the radius parameter r in the VCM estimator.",
true);
776 app.add_option(
"--method,-m", method,
"the method of tangent computation: VCM (default), L-MST.",
true)
777 -> check(CLI::IsMember({
"VCM",
"L-MST"}));
778 app.add_option(
"--axes,-a", axesFlag,
"show main axes - prints list of axes for each point and color points color = (if X => 255, if Y => 255, if Z => 255)",
true)
779 -> check(CLI::IsMember({
"ON",
"OFF"}));
780 app.add_option(
"--output,-o",outputFileName,
"the basename of the output text file which will contain points and tangent vectors: (x y z tx ty tz) per line",
true);
784 app.get_formatter()->column_width(40);
785 CLI11_PARSE(app, argc, argv);
791 vector<Point> sequence;
793 inputStream.open ( inputFileName.c_str(), ios::in);
806 bool view = viewFlag ==
"ON";
809 bool axes = axesFlag ==
"ON";
810 multimap < Z3i::Point, string > mainAxes;
812 if ( axes && connectivity ==
"26" )
814 find_main_axes<PointIterator, Z3i::Space>(sequence.begin(), sequence.end(), mainAxes);
815 print_main_axes<PointIterator, Z3i::Space>(sequence.begin(), sequence.end(), mainAxes);
817 else if ( axes && connectivity ==
"6" )
819 find_main_axes < PointIterator, Z3i::Space, 4 > ( sequence.begin ( ), sequence.end ( ), mainAxes );
820 print_main_axes<PointIterator, Z3i::Space>(sequence.begin(), sequence.end(), mainAxes);
825 Point lowerBound = sequence[ 0 ];
826 Point upperBound = sequence[ 0 ];
827 for (
unsigned int j = 1; j < sequence.size(); ++j )
829 lowerBound = lowerBound.inf( sequence[ j ] );
830 upperBound = upperBound.sup( sequence[ j ] );
832 lowerBound -= Point::diagonal( box );
833 upperBound += Point::diagonal( box + 1 );
834 K3 ks; ks.
init( lowerBound, upperBound,
true );
838 std::vector< RealVector > tangents;
840 if ( method ==
"VCM" )
843 const double R = bigRad;
844 trace.
info() <<
"Big radius R = " << R << endl;
845 const double r = smallRad;
846 trace.
info() <<
"Small radius r = " << r << endl;
848 ComputeVCM < PointIterator, Z3, std::vector< RealVector > > ( R, r, sequence.begin(), sequence.end(), tangents, outputFileName );
850 else if ( method ==
"L-MST" )
852 if (connectivity ==
"6")
853 ComputeLMST6 < PointIterator, Z3, std::vector< RealVector > > ( sequence.begin(), sequence.end(), tangents, outputFileName );
855 ComputeLMST26 < PointIterator, Z3, std::vector< RealVector > > ( sequence.begin(), sequence.end(), tangents, outputFileName );
859 trace.
info() <<
"Wrong method! Try: VCM or L-MST" << endl;
865 for (
unsigned int i = 0; i < tangents.size(); i++ )
872 viewer.
addLine( p + 2.0 * tangent, p - 2.0 * tangent, 5.0 );
882 gc.initFromPointsVector( sequence );
886 trace.
warning() <<
"[ConnectivityException] GridCurve only accepts a sequence of face adjacent points. Try connectivity=6 instead." << endl;
897 displayAxes<Point,RealPoint, Z3i::Space, Z3i::KSpace>( viewer, lowerBound, upperBound, viewBox );
899 res = connectivity ==
"6"
900 ? displayCover6( viewer, ks, sequence.begin(), sequence.end(),
906 : displayCover26( viewer, ks, sequence.begin(), sequence.end(),
913 if ( curve3d && ! axes )
916 for ( vector<Point>::const_iterator it = sequence.begin(); it != sequence.end(); ++it )
919 else if ( curve3d && axes )
921 for ( vector<Point>::const_iterator itt = sequence.begin(); itt != sequence.end(); ++itt ) {
923 typename std::multimap<typename Z3i::Point, string>::const_iterator it = mainAxes.lower_bound(*itt);
924 typename std::multimap<typename Z3i::Point, string>::const_iterator it2 = mainAxes.upper_bound(*itt);
926 for (; it != it2; it++ )
928 if (it->second ==
"X")
930 if (it->second ==
"Y")
932 if (it->second ==
"Z")
936 Color c ( pColor[0], pColor[1], pColor[2], 255 );
943 viewer << Viewer3D<Z3,K3>::updateDisplay;
946 trace.
emphase() << ( res ?
"Passed." :
"Error." ) << endl;
int main(int argc, char **argv)
const Primitive & primitive() const
virtual void setLineColor(DGtal::Color aColor)
virtual void setFillColor(DGtal::Color aColor)
void addLine(const RealPoint &p1, const RealPoint &p2, const double width=0.03)
void addBall(const RealPoint ¢er, const double radius=0.5, const unsigned int resolution=30)
void addQuad(const RealPoint &p1, const RealPoint &p2, const RealPoint &p3, const RealPoint &p4)
typename Self::Domain Domain
typename Self::Point Point
bool init(const Point &lower, const Point &upper, bool isClosed)
const Point & lowerBound() const
Cell uCell(const PreCell &c) const
RealVector eval(const ConstIterator &it)
void init(const ConstIterator &itb, const ConstIterator &ite)
void attach(Alias< DSSSegmentationComputer > segmentComputer)
IteratorCirculatorTraits< ConstIterator >::Value Point3d
void getParameters(Vector3d &direction, PointR3d &intercept, PointR3d &thickness) const
ConstIterator begin() const
bool validArithmeticalDSS2d(Dimension i) const
ConstIterator end() const
const ArithmeticalDSSComputer2d & arithmeticalDSS2dXY() const
const ArithmeticalDSSComputer2d & arithmeticalDSS2d(Dimension i) const
std::array< Quotient, 3 > PointR3d
const ArithmeticalDSSComputer2d & arithmeticalDSS2dXZ() const
const ArithmeticalDSSComputer2d & arithmeticalDSS2dYZ() const
double norm(const NormType type=L_2) const
const ArithmeticalDSSComputer2d & arithmeticalDSS2d(Dimension i) const
void beginBlock(const std::string &keyword="")
Space::RealVector RealVector
Space::RealPoint RealPoint
Trace trace(traceWriterTerm)
DGtal::uint32_t Dimension
static std::vector< TPoint > getPointsFromInputStream(std::istream &in, std::vector< unsigned int > aVectPosition=std::vector< unsigned int >())