DGtalTools  1.2.0
3dCurveTangentEstimator.cpp
1 
38 #include <iostream>
39 #include <iterator>
40 #include <cstdio>
41 #include <cmath>
42 #include <string>
43 #include <map>
44 #include <fstream>
45 #include <vector>
46 #include <utility>
47 
48 #include "CLI11.hpp"
49 
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"
68 
69 
70 using namespace DGtal;
71 using namespace std;
72 
73 
74 
75 
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;
145 
147 // Functions for displaying the tangential cover of a 3D curve.
149 template <typename Point, typename RealPoint, typename space, typename kspace >
150 void displayAxes( Viewer3D<space, kspace> & viewer,
151  const Point & lowerBound, const Point & upperBound,
152  const std::string & mode )
153 {
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" ) )
161  {
162  viewer.setLineColor(AXIS_COLOR);
163  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p0[ 2 ]),
164  DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p0[ 2 ]), AXIS_LINESIZE );
165  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p0[ 2 ]),
166  DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p0[ 2 ]), AXIS_LINESIZE );
167  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p0[ 2 ]),
168  DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
169  viewer.addLine( DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p0[ 2 ]),
170  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p0[ 2 ]), AXIS_LINESIZE );
171  viewer.addLine( DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p0[ 2 ]),
172  DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
173  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p0[ 2 ]),
174  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p0[ 2 ]), AXIS_LINESIZE );
175  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p0[ 2 ]),
176  DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
177  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p1[ 2 ]),
178  DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
179  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p1[ 2 ]),
180  DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
181  viewer.addLine( DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p0[ 2 ]),
182  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
183  viewer.addLine( DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p1[ 2 ]),
184  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
185  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p1[ 2 ]),
186  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
187  }
188  if ( mode == "COLORED" )
189  {
190  viewer.setFillColor(XY_COLOR);
191  viewer.addQuad(DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p1[ 2 ]),
192  DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p1[ 2 ]),
193  DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p1[ 2 ]),
194  DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p1[ 2 ]) );
195  viewer.setFillColor(XZ_COLOR);
196  viewer.addQuad(DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p1[ 2 ]),
197  DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p1[ 2 ]),
198  DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p0[ 2 ]),
199  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p0[ 2 ]));
200  viewer.setFillColor(YZ_COLOR);
201  viewer.addQuad(DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p1[ 2 ]),
202  DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p1[ 2 ]),
203  DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p0[ 2 ]),
204  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p0[ 2 ]));
205  }
206 }
207 
208 template <typename KSpace, typename Naive3DDSSComputer, typename space, typename kspace >
209 void displayDSS3d( Viewer3D<space, kspace> & viewer,
210  const KSpace & ks, const Naive3DDSSComputer & dss3d,
211  const DGtal::Color & color3d )
212 {
213  viewer << CustomColors3D( color3d, color3d ) << dss3d;
214 }
215 
216 template <typename Point1, typename Point2>
217 void assign( Point1 & p1, const Point2 & p2 )
218 {
219  p1[ 0 ] = p2[ 0 ];
220  p1[ 1 ] = p2[ 1 ];
221  p1[ 2 ] = p2[ 2 ];
222 }
223 
224 template <typename KSpace, typename Naive3DDSSComputer, typename space, typename kspace >
225 void displayDSS3dTangent( Viewer3D<space, kspace> & viewer,
226  const KSpace & ks, const Naive3DDSSComputer & dss3d,
227  const DGtal::Color & color3d )
228 {
229  typedef typename Naive3DDSSComputer::Point3d Point;
230  typedef typename Naive3DDSSComputer::Integer Integer;
231  typedef typename Naive3DDSSComputer::PointR3d PointR3d;
232  typedef typename Display3D<>::BallD3D PointD3D;
233  Point directionZ3;
234  PointR3d interceptR, thicknessR;
235  Z3i::RealPoint P1, P2, direction, intercept, thickness;
236  dss3d.getParameters( directionZ3, interceptR, thicknessR );
237 
238  intercept[0] = (double) NumberTraits<Integer>::castToDouble ( interceptR[0].first ) / (double) NumberTraits<Integer>::castToDouble ( interceptR[0].second );
239  intercept[1] = (double) NumberTraits<Integer>::castToDouble ( interceptR[1].first ) / (double) NumberTraits<Integer>::castToDouble ( interceptR[1].second );
240  intercept[2] = (double) NumberTraits<Integer>::castToDouble ( interceptR[2].first ) / (double) NumberTraits<Integer>::castToDouble ( interceptR[2].second );
241  thickness[0] = (double) NumberTraits<Integer>::castToDouble ( thicknessR[0].first ) / (double) NumberTraits<Integer>::castToDouble ( thicknessR[0].second );
242  thickness[1] = (double) NumberTraits<Integer>::castToDouble ( thicknessR[1].first ) / (double) NumberTraits<Integer>::castToDouble ( thicknessR[1].second );
243  thickness[2] = (double) NumberTraits<Integer>::castToDouble ( thicknessR[2].first ) / (double) NumberTraits<Integer>::castToDouble ( thicknessR[2].second );
244 
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 );
251 
252  Z3i::RealPoint Q1 = intercept + t1 * direction;
253  Z3i::RealPoint Q2 = intercept + t2 * direction;
254  viewer.setLineColor(color3d);
255  viewer.addLine( Z3i::RealPoint(Q1[ 0 ]-0.5, Q1[ 1 ]-0.5, Q1[ 2 ]-0.5),
256  Z3i::RealPoint(Q2[ 0 ]-0.5, Q2[ 1 ]-0.5, Q2[ 2 ]-0.5),
257  MS3D_LINESIZE );
258 }
259 
260 template <typename KSpace, typename StandardDSS6Computer, typename space, typename kspace >
261 void displayProj2d6( Viewer3D<space, kspace> & viewer,
262  const KSpace & ks, const StandardDSS6Computer & dss3d,
263  const DGtal::Color & color2d )
264 {
265  typedef typename StandardDSS6Computer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
266  typedef typename ArithmeticalDSSComputer2d::ConstIterator ConstIterator2d;
267  typedef typename ArithmeticalDSSComputer2d::Point Point2d;
268  typedef typename KSpace::Cell Cell;
269  typedef typename KSpace::Point Point3d;
270  Point3d b = ks.lowerBound();
271  for ( DGtal::Dimension i = 0; i < 3; ++i )
272  {
273  const ArithmeticalDSSComputer2d & dss2d = dss3d.arithmeticalDSS2d( i );
274  for ( ConstIterator2d itP = dss2d.begin(), itPEnd = dss2d.end(); itP != itPEnd; ++itP )
275  {
276  Point2d p = *itP;
277  Point3d q;
278  switch (i)
279  {
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;
283  }
284  Cell c = ks.uCell( q );
285  viewer << CustomColors3D( color2d, color2d ) << c;
286  }
287  }
288 }
289 
290 template <typename KSpace, typename StandardDSS6Computer, typename space, typename kspace >
291 void displayDSS2d6( Viewer3D<space, kspace> & viewer,
292  const KSpace & ks, const StandardDSS6Computer & dss3d,
293  const DGtal::Color & color2d )
294 {
295  typedef typename StandardDSS6Computer::ConstIterator ConstIterator3d;
296  typedef typename StandardDSS6Computer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
297  typedef typename ArithmeticalDSSComputer2d::ConstIterator ConstIterator2d;
298  typedef typename ArithmeticalDSSComputer2d::Point Point2d;
299  typedef typename KSpace::Cell Cell;
300  typedef typename KSpace::Point Point3d;
301  typedef DGtal::PointVector<2,double> PointD2d;
302  typedef typename Display3D<>::BallD3D PointD3D;
303  Point3d b = ks.lowerBound();
304  for ( DGtal::Dimension i = 0; i < 3; ++i )
305  {
306  const typename ArithmeticalDSSComputer2d::Primitive & dss2d
307  = dss3d.arithmeticalDSS2d( i ).primitive();
308  // draw 2D bounding boxes for each arithmetical dss 2D.
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;
315  PointD3D p3;
316  for ( unsigned int j = 0; j < pts2d.size(); ++j )
317  {
318  switch (i)
319  {
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;
323  }
324  bb.push_back( p3 );
325  }
326  for ( unsigned int j = 0; j < pts2d.size(); ++j ){
327  viewer.setLineColor(color2d);
328  viewer.addLine( DGtal::Z3i::RealPoint(bb[ j ].center[0], bb[ j ].center[1], bb[ j ].center[2]),
329  DGtal::Z3i::RealPoint(bb[ (j+1)%4 ].center[0], bb[ (j+1)%4 ].center[1], bb[ (j+1)%4 ].center[2]),
330  MS3D_LINESIZE );
331  }
332  } // for ( DGtal::Dimension i = 0; i < 3; ++i )
333 }
334 
335 
336 //Why not to just project 3D curve?
337 template <typename KSpace, typename Naive3DDSSComputer, typename space, typename kspace >
338 void displayProj2d26( Viewer3D<space, kspace> & viewer,
339  const KSpace & ks, const Naive3DDSSComputer & dss3d,
340  const DGtal::Color & color2d )
341 {
342  typedef typename Naive3DDSSComputer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
343  typedef typename ArithmeticalDSSComputer2d::ConstIterator ConstIterator2d;
344  typedef typename ArithmeticalDSSComputer2d::Point Point2d;
345  typedef typename KSpace::Cell Cell;
346  typedef typename KSpace::Point Point3d;
347  Point3d b = ks.lowerBound();
348 
349  const ArithmeticalDSSComputer2d & dssXY = dss3d.arithmeticalDSS2dXY();
350  const ArithmeticalDSSComputer2d & dssXZ = dss3d.arithmeticalDSS2dXZ();
351  const ArithmeticalDSSComputer2d & dssYZ = dss3d.arithmeticalDSS2dYZ();
352 
353  bool validXY = dss3d.validArithmeticalDSS2d ( 2 );
354  bool validXZ = dss3d.validArithmeticalDSS2d ( 1 );
355  bool validYZ = dss3d.validArithmeticalDSS2d ( 0 );
356 
357  if ( validXY && validXZ )
358  { //XY-plane, XZ-plane
359  for ( ConstIterator2d itXY = dssXY.begin(), itXZ = dssXZ.begin(), itPEnd = dssXY.end(); itXY != itPEnd; ++itXY, ++itXZ )
360  {
361  Point2d p1 = *itXY, p2 = *itXZ;
362  Point3d q1, q2, q3;
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 ] );
366  Cell c1 = ks.uCell( q1 ); Cell c2 = ks.uCell( q2 ); Cell c3 = ks.uCell( q3 );
367  viewer << CustomColors3D( color2d, color2d ) << c1;
368  viewer << CustomColors3D( color2d, color2d ) << c2;
369  viewer << CustomColors3D( color2d, color2d ) << c3;
370  }
371  }
372  else
373  {
374  if ( validYZ && validXY )
375  { //XY-plane, YZ-plane
376  for ( ConstIterator2d itYZ = dssYZ.begin(), itXY = dssXY.begin(), itPEnd = dssXY.end(); itXY != itPEnd; ++itXY, ++itYZ )
377  {
378  Point2d p1 = *itYZ, p2 = *itXY;
379  Point3d q1, q2, q3;
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 ] );
383  Cell c1 = ks.uCell( q1 ); Cell c2 = ks.uCell( q2 ); Cell c3 = ks.uCell( q3 );
384  viewer << CustomColors3D( color2d, color2d ) << c1;
385  viewer << CustomColors3D( color2d, color2d ) << c2;
386  viewer << CustomColors3D( color2d, color2d ) << c3;
387  }
388  }
389  else
390  {
391  for ( ConstIterator2d itYZ = dssYZ.begin(), itXZ = dssXZ.begin(), itPEnd = dssXZ.end(); itXZ != itPEnd; ++itXZ, ++itYZ )
392  {
393  Point2d p1 = *itYZ, p2 = *itXZ;
394  Point3d q1, q2, q3;
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 ] );
398  Cell c1 = ks.uCell( q1 ); Cell c2 = ks.uCell( q2 );Cell c3 = ks.uCell( q3 );
399  viewer << CustomColors3D( color2d, color2d ) << c1;
400  viewer << CustomColors3D( color2d, color2d ) << c2;
401  viewer << CustomColors3D( color2d, color2d ) << c3;
402  }
403  }
404  }
405 }
406 
407 
408 template <typename KSpace, typename Naive3DDSSComputer, typename space, typename kspace >
409 void displayDSS2d26( Viewer3D<space, kspace> & viewer,
410  const KSpace & ks, const Naive3DDSSComputer & dss3d,
411  const DGtal::Color & color2d )
412 {
413  typedef typename Naive3DDSSComputer::ConstIterator ConstIterator3d;
414  typedef typename Naive3DDSSComputer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
415  typedef typename ArithmeticalDSSComputer2d::ConstIterator ConstIterator2d;
416  typedef typename ArithmeticalDSSComputer2d::Point Point2d;
417  typedef typename KSpace::Cell Cell;
418  typedef typename KSpace::Point Point3d;
419  typedef DGtal::PointVector<2,double> PointD2d;
420  typedef typename Display3D<>::BallD3D PointD3D;
421  Point3d b = ks.lowerBound();
422  for ( DGtal::Dimension i = 0; i < 3; ++i )
423  {
424  if ( !dss3d.validArithmeticalDSS2d ( i ) )
425  continue;
426 
427  const typename ArithmeticalDSSComputer2d::Primitive & dss2d
428  = dss3d.arithmeticalDSS2d( i ).primitive();
429  // draw 2D bounding boxes for each arithmetical dss 2D.
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;
436  PointD3D p3;
437  for ( unsigned int j = 0; j < pts2d.size(); ++j )
438  {
439  switch (i)
440  {
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;
444  }
445  bb.push_back( p3 );
446  }
447  for ( unsigned int j = 0; j < pts2d.size(); ++j ){
448  viewer.setLineColor(color2d);
449  viewer.addLine( DGtal::Z3i::RealPoint(bb[ j ].center[0], bb[ j ].center[1], bb[ j ].center[2]),
450  DGtal::Z3i::RealPoint(bb[ (j+1)%4 ].center[0], bb[ (j+1)%4 ].center[1], bb[ (j+1)%4 ].center[2]),
451  MS3D_LINESIZE );
452  }
453  } // for ( DGtal::Dimension i = 0; i < 3; ++i )
454 }
455 
459 template <typename KSpace, typename PointIterator, typename space, typename kspace >
460 bool displayCover6( Viewer3D<space, kspace> & viewer,
461  const KSpace & ks, PointIterator b, PointIterator e,
462  bool dss3d, bool proj2d, bool dss2d, bool tangent,
463  int nbColors )
464 {
465  typedef typename PointIterator::value_type Point;
466  typedef StandardDSS6Computer<PointIterator,int,4> SegmentComputer;
467  typedef SaturatedSegmentation<SegmentComputer> Decomposition;
468  typedef typename Decomposition::SegmentComputerIterator SegmentComputerIterator;
469  typedef typename SegmentComputer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
470  SegmentComputer algo;
471  Decomposition theDecomposition(b, e, algo);
472 
473  viewer << SetMode3D( algo.className(), "BoundingBox" );
474  HueShadeColorMap<int> cmap_hue( 0, nbColors, 1 );
475 
476  unsigned int c = 0;
477  for ( SegmentComputerIterator i = theDecomposition.begin();
478  i != theDecomposition.end(); ++i)
479  {
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);
486  trace.info() << "- " << c
487  << " MS3D,"
488  << " [" << f[ 0 ] << "," << f[ 1 ] << ","<< f[ 2 ] << "]"
489  << "->[" << l[ 0 ] << "," << l[ 1 ] << ","<< l[ 2 ] << "]"
490  << ", XY("
491  << dssXY.a() << "," << dssXY.b() << "," << dssXY.mu()
492  << "), XZ("
493  << dssXZ.a() << "," << dssXZ.b() << "," << dssXZ.mu()
494  << "), YZ("
495  << dssYZ.a() << "," << dssYZ.b() << "," << dssYZ.mu()
496  << ")" << std::endl;
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 );
502  c++;
503  }
504  return true;
505 }
506 
511 template <typename KSpace, typename PointIterator, typename space, typename kspace >
512 bool displayCover26( Viewer3D<space, kspace> & viewer,
513  const KSpace & ks, PointIterator b, PointIterator e,
514  bool dss3d, bool proj2d, bool dss2d, bool tangent,
515  int nbColors )
516 {
517  typedef typename PointIterator::value_type Point;
518  typedef Naive3DDSSComputer<PointIterator,int, 8 > SegmentComputer;
519  typedef SaturatedSegmentation<SegmentComputer> Decomposition;
520  typedef typename Decomposition::SegmentComputerIterator SegmentComputerIterator;
521  typedef typename SegmentComputer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
522  SegmentComputer algo;
523  Decomposition theDecomposition(b, e, algo);
524 
525  viewer << SetMode3D( algo.className(), "BoundingBox" );
526  HueShadeColorMap<int> cmap_hue( 0, nbColors, 1 );
527 
528  unsigned int c = 0;
529  for ( SegmentComputerIterator i = theDecomposition.begin();
530  i != theDecomposition.end(); ++i)
531  {
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);
538  trace.info() << "- " << c
539  << " MS3D,"
540  << " [" << f[ 0 ] << "," << f[ 1 ] << ","<< f[ 2 ] << "]"
541  << "->[" << l[ 0 ] << "," << l[ 1 ] << ","<< l[ 2 ] << "]"
542  << ", XY("
543  << dssXY.a() << "," << dssXY.b() << "," << dssXY.mu()
544  << "), XZ("
545  << dssXZ.a() << "," << dssXZ.b() << "," << dssXZ.mu()
546  << "), YZ("
547  << dssYZ.a() << "," << dssYZ.b() << "," << dssYZ.mu()
548  << ")" << std::endl;
549  //trace.info() << ms3d << std::endl; // information
550 
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 );
556  c++;
557  }
558  return true;
559 }
560 
561 template < typename PointIterator, 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 )
564 {
565  using namespace DGtal;
566  typedef ExactPredicateLpSeparableMetric < Space, 2 > Metric; // L2-metric
567  typedef VoronoiCovarianceMeasure < Space, Metric > VCM;
568  typedef HyperRectDomain < Space > Domain;
569  typedef functors::HatPointFunction < typename Space::Point, double > KernelFunction;
570  typedef EigenDecomposition < 3, double > LinearAlgebraTool;
571  typedef LinearAlgebraTool::Matrix Matrix;
572 
573  fstream outputStream;
574  outputStream.open ( ( output + ".vcm" ).c_str(), std::ios::out );
575  outputStream << "# VCM estimation R=" << R << " r=" << r << " chi=hat" << endl;
576 
577  Metric l2;
578  VCM vcm ( R, ceil( r ), l2, true );
579  vcm.init( begin, end );
580  KernelFunction chi( 1.0, r );
581  Matrix vcm_r, evec;
582  typename Space::RealVector eval;
583  for ( PointIterator it = begin; it != end; ++it )
584  {
585  // Compute VCM and diagonalize 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;
592  }
593  outputStream.close();
594 }
595 
596 template < typename PointIterator, typename Space, typename TangentSequence >
597 void ComputeLMST6 ( const PointIterator & begin, const PointIterator & end, TangentSequence & tangents, const std::string & output )
598 {
599  typedef typename PointIterator::value_type Point;
600  typedef StandardDSS6Computer<PointIterator,int, 4 > SegmentComputer;
601  typedef SaturatedSegmentation<SegmentComputer> Decomposition;
602  typedef typename Decomposition::SegmentComputerIterator SegmentComputerIterator;
603  typedef typename SegmentComputer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
604  SegmentComputer algo;
605  Decomposition theDecomposition ( begin, end, algo );
606 
607  fstream outputStream;
608  outputStream.open ( ( output + ".lmst" ).c_str(), std::ios::out );
609  outputStream << "X Y Z X Y Z" << endl;
610 
611  LambdaMST3D < Decomposition > lmst64;
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 )
617  {
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;
623  }
624 }
625 
626 template < typename PointIterator, typename Space, typename TangentSequence >
627 void ComputeLMST26 ( const PointIterator & begin, const PointIterator & end, TangentSequence & tangents, const std::string & output )
628 {
629  typedef typename PointIterator::value_type Point;
630  typedef Naive3DDSSComputer<PointIterator,int, 8 > SegmentComputer;
631  typedef SaturatedSegmentation<SegmentComputer> Decomposition;
632  typedef typename Decomposition::SegmentComputerIterator SegmentComputerIterator;
633  typedef typename SegmentComputer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
634  SegmentComputer algo;
635  Decomposition theDecomposition ( begin, end, algo );
636 
637  fstream outputStream;
638  outputStream.open ( ( output + ".lmst" ).c_str(), std::ios::out );
639  outputStream << "X Y Z X Y Z" << endl;
640 
641  LambdaMST3D < Decomposition > lmst64;
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 )
647  {
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;
653  }
654 }
655 
656 template < typename PointIterator, typename Space, int CONNECT = 8 >
657 void find_main_axes ( PointIterator begin, PointIterator end, multimap < typename Space::Point, string > & axes )
658 {
659  typedef Naive3DDSSComputer < PointIterator, int, CONNECT > SegmentComputer;
660  typedef SaturatedSegmentation < SegmentComputer > Segmentation;
661 
662  Segmentation segmenter ( begin, end, SegmentComputer ( ) );
663 
664  for ( auto it = begin; it != end; ++it )
665  {
666  unsigned int cX = 0, cY = 0, cZ = 0;
667  for ( typename Segmentation::SegmentComputerIterator idss = segmenter.begin ( ); idss != segmenter.end ( ); ++idss )
668  {
669  if ( idss->isInDSS ( *it ) )
670  {
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() );
677 
678  if ( lenXY >= lenYZ && lenXZ >= lenYZ )
679  cX++;
680  else if ( lenXY >= lenXZ && lenYZ >= lenXZ )
681  cY++;
682  else
683  cZ++;
684  }
685  }
686  if ( cX > 0 )
687  axes.insert ( pair <typename Space::Point, string > ( *it, string ( "X" ) ) );
688  if ( cY > 0 )
689  axes.insert ( pair <typename Space::Point, string > ( *it, string ( "Y" ) ) );
690  if ( cZ > 0 )
691  axes.insert ( pair <typename Space::Point, string > ( *it, string ( "Z" ) ) );
692 
693  }
694 }
695 
696 template < typename PointIterator, typename Space >
697 void print_main_axes ( PointIterator begin, PointIterator end, multimap < typename Space::Point, string > & axes )
698 {
699  for ( PointIterator itt = begin; itt != end; ++itt )
700  {
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++ )
705  {
706  if ( distance( it, it2 ) > 1 )
707  cout << it->second << ", ";
708  else
709  cout << it->second << ")";
710  }
711  cout << endl;
712  }
713 }
714 
723 int main(int argc, char **argv)
724 {
725  using namespace std;
726  using namespace DGtal;
727  typedef SpaceND<3,int> Z3;
728  typedef KhalimskySpaceND<3,int> K3;
729  typedef Z3::Point Point;
730  typedef Z3::RealPoint RealPoint;
731  typedef Z3::RealVector RealVector;
732  typedef HyperRectDomain<Z3> Domain;
733  typedef typename vector<Point>::const_iterator PointIterator;
734 
735  // specify command line ----------------------------------------------
736  QApplication application(argc,argv); // remove Qt arguments.
737 
738  // parse command line using CLI ----------------------------------------------
739  CLI::App app;
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};
754  int box {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." )
758  ->required()
759  ->check(CLI::ExistingFile);
760 
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);
773 
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);
781 
782 
783 
784  app.get_formatter()->column_width(40);
785  CLI11_PARSE(app, argc, argv);
786  // END parse command line using CLI ----------------------------------------------
787 
788 
789 
790  // process command line ----------------------------------------------
791  vector<Point> sequence;
792  fstream inputStream;
793  inputStream.open ( inputFileName.c_str(), ios::in);
794  try
795  {
796  sequence = PointListReader<Point>::getPointsFromInputStream( inputStream );
797  if ( sequence.size() == 0) throw IOException();
798  }
799  catch (DGtal::IOException & ioe)
800  {
801  trace.error() << "Size is null." << endl;
802  }
803  inputStream.close();
804 
805  // start viewer
806  bool view = viewFlag == "ON";
807 
808  // axes
809  bool axes = axesFlag == "ON";
810  multimap < Z3i::Point, string > mainAxes;
811 
812  if ( axes && connectivity == "26" )
813  {
814  find_main_axes<PointIterator, Z3i::Space>(sequence.begin(), sequence.end(), mainAxes);
815  print_main_axes<PointIterator, Z3i::Space>(sequence.begin(), sequence.end(), mainAxes);
816  }
817  else if ( axes && connectivity == "6" )
818  {
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);
821  }
822 
823  // ----------------------------------------------------------------------
824  // Create domain and curve.
825  Point lowerBound = sequence[ 0 ];
826  Point upperBound = sequence[ 0 ];
827  for ( unsigned int j = 1; j < sequence.size(); ++j )
828  {
829  lowerBound = lowerBound.inf( sequence[ j ] );
830  upperBound = upperBound.sup( sequence[ j ] );
831  }
832  lowerBound -= Point::diagonal( box );
833  upperBound += Point::diagonal( box + 1 );
834  K3 ks; ks.init( lowerBound, upperBound, true );
835  Viewer3D<Z3,K3> viewer( ks );
836  trace.beginBlock ( "Tool 3dCurveTangentEstimator" );
837 
838  std::vector< RealVector > tangents;
839 
840  if ( method == "VCM" )
841  {
842  // input points of the curve are in sequence vector.
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;
847 
848  ComputeVCM < PointIterator, Z3, std::vector< RealVector > > ( R, r, sequence.begin(), sequence.end(), tangents, outputFileName );
849  }
850  else if ( method == "L-MST" )
851  {
852  if (connectivity == "6")
853  ComputeLMST6 < PointIterator, Z3, std::vector< RealVector > > ( sequence.begin(), sequence.end(), tangents, outputFileName );
854  else
855  ComputeLMST26 < PointIterator, Z3, std::vector< RealVector > > ( sequence.begin(), sequence.end(), tangents, outputFileName );
856  }
857  else
858  {
859  trace.info() << "Wrong method! Try: VCM or L-MST" << endl;
860  abort();
861  }
862 
863  if ( view )
864  {
865  for ( unsigned int i = 0; i < tangents.size(); i++ )
866  {
867  // Display normal
868  RealPoint p = sequence[i];
869  RealVector tangent = tangents[i];
870  viewer.setFillColor( Color ( 255, 0, 0, 255 ) );
871  viewer.setLineColor( Color ( 255, 0, 0, 255 ) );
872  viewer.addLine( p + 2.0 * tangent, p - 2.0 * tangent, 5.0 );
873  viewer.setFillColor( Color ( 100, 100, 140, 255 ) );
874  viewer.setLineColor( Color ( 100, 100, 140, 255 ) );
875  viewer.addBall( p, 0.125, 8 );
876  }
877  }
878 
879  GridCurve<K3> gc( ks );
880  try
881  {
882  gc.initFromPointsVector( sequence );
883  }
884  catch (DGtal::ConnectivityException& /*ce*/)
885  {
886  trace.warning() << "[ConnectivityException] GridCurve only accepts a sequence of face adjacent points. Try connectivity=6 instead." << endl;
887  }
888 
889  // ----------------------------------------------------------------------
890  // Displays everything.
891  bool res = true;
892  if ( view )
893  {
894  viewer.show();
895  // Display axes.
896  if ( viewBox != "" )
897  displayAxes<Point,RealPoint, Z3i::Space, Z3i::KSpace>( viewer, lowerBound, upperBound, viewBox );
898  // Display 3D tangential cover.
899  res = connectivity == "6"
900  ? displayCover6( viewer, ks, sequence.begin(), sequence.end(),
901  cover3d,
902  curve2d,
903  cover2d,
904  displayTangent,
905  nbColors )
906  : displayCover26( viewer, ks, sequence.begin(), sequence.end(),
907  cover3d,
908  curve2d,
909  cover2d,
910  displayTangent,
911  nbColors );
912  // Display 3D curve points.
913  if ( curve3d && ! axes )
914  {
915  viewer << CustomColors3D( CURVE3D_COLOR, CURVE3D_COLOR );
916  for ( vector<Point>::const_iterator it = sequence.begin(); it != sequence.end(); ++it )
917  viewer << *it;
918  }
919  else if ( curve3d && axes )
920  {
921  for ( vector<Point>::const_iterator itt = sequence.begin(); itt != sequence.end(); ++itt ) {
922 
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);
925  Z3i::Point pColor;
926  for (; it != it2; it++ )
927  {
928  if (it->second == "X")
929  pColor[0] = 255;
930  if (it->second == "Y")
931  pColor[1] = 255;
932  if (it->second == "Z")
933  pColor[2] = 255;
934 
935  }
936  Color c ( pColor[0], pColor[1], pColor[2], 255 );
937  viewer.setFillColor(c);
938  viewer << *itt;
939  }
940  }
941  // ----------------------------------------------------------------------
942  // User "interaction".
943  viewer << Viewer3D<Z3,K3>::updateDisplay;
944  application.exec();
945  }
946  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
947  trace.endBlock();
948 
949  return res ? 0 : 1;
950 }
DGtal
Definition: ATu0v1.h:57
Companion project associated to DGtal ; Generated on Wed Jun 2 2021 15:34:51 for DGtalTools by   doxygen 1.9.1