DGtalTools  0.9.2
3dCurveTangentEstimator.cpp
1 
38 #include <iostream>
39 #include <iterator>
40 #include <cstdio>
41 #include <cmath>
42 #include <fstream>
43 #include <vector>
44 
45 #include <boost/program_options/options_description.hpp>
46 #include <boost/program_options/parsers.hpp>
47 #include <boost/program_options/variables_map.hpp>
48 
49 #include "DGtal/base/Common.h"
50 #include "DGtal/base/Exceptions.h"
51 #include "DGtal/kernel/SpaceND.h"
52 #include "DGtal/kernel/domains/HyperRectDomain.h"
53 #include "DGtal/kernel/BasicPointPredicates.h"
54 #include "DGtal/math/linalg/EigenDecomposition.h"
55 #include "DGtal/topology/KhalimskySpaceND.h"
56 #include "DGtal/geometry/curves/GridCurve.h"
57 #include "DGtal/geometry/curves/Naive3DDSSComputer.h"
58 #include "DGtal/geometry/curves/StandardDSS6Computer.h"
59 #include "DGtal/geometry/curves/SaturatedSegmentation.h"
60 #include "DGtal/geometry/volumes/distance/ExactPredicateLpSeparableMetric.h"
61 #include "DGtal/geometry/volumes/estimation/VoronoiCovarianceMeasure.h"
62 #include "DGtal/geometry/curves/estimation/LambdaMST3D.h"
63 #include "DGtal/geometry/curves/estimation/FunctorsLambdaMST.h"
64 #include "DGtal/io/viewers/Viewer3D.h"
65 #include "DGtal/io/readers/PointListReader.h"
66 #include "DGtal/io/DrawWithDisplay3DModifier.h"
67 
68 
69 using namespace DGtal;
70 using namespace std;
71 
72 
73 
74 
150 const Color AXIS_COLOR( 0, 0, 0, 255 );
151 const double AXIS_LINESIZE = 0.1;
152 const Color XY_COLOR( 0, 0, 255, 50 );
153 const Color XZ_COLOR( 0, 255, 0, 50 );
154 const Color YZ_COLOR( 255, 0, 0, 50 );
155 const Color CURVE3D_COLOR( 100, 100, 140, 128 );
156 const Color CURVE2D_COLOR( 200, 200, 200, 100 );
157 const double MS3D_LINESIZE = 0.05;
158 
160 // Functions for displaying the tangential cover of a 3D curve.
162 template <typename Point, typename RealPoint, typename space, typename kspace >
163 void displayAxes( Viewer3D<space, kspace> & viewer,
164  const Point & lowerBound, const Point & upperBound,
165  const std::string & mode )
166 {
167  RealPoint p0( (double)lowerBound[ 0 ]-0.5,
168  (double)lowerBound[ 1 ]-0.5,
169  (double)lowerBound[ 2 ]-0.5 );
170  RealPoint p1( (double)upperBound[ 0 ]-0.5,
171  (double)upperBound[ 1 ]-0.5,
172  (double)upperBound[ 2 ]-0.5 );
173  if ( ( mode == "WIRED" ) || ( mode == "COLORED" ) )
174  {
175  viewer.setLineColor(AXIS_COLOR);
176  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p0[ 2 ]),
177  DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p0[ 2 ]), AXIS_LINESIZE );
178  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p0[ 2 ]),
179  DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p0[ 2 ]), AXIS_LINESIZE );
180  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p0[ 2 ]),
181  DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
182  viewer.addLine( DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p0[ 2 ]),
183  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p0[ 2 ]), AXIS_LINESIZE );
184  viewer.addLine( DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p0[ 2 ]),
185  DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
186  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p0[ 2 ]),
187  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p0[ 2 ]), AXIS_LINESIZE );
188  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p0[ 2 ]),
189  DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
190  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p1[ 2 ]),
191  DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
192  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p1[ 2 ]),
193  DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
194  viewer.addLine( DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p0[ 2 ]),
195  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
196  viewer.addLine( DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p1[ 2 ]),
197  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
198  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p1[ 2 ]),
199  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
200  }
201  if ( mode == "COLORED" )
202  {
203  viewer.setFillColor(XY_COLOR);
204  viewer.addQuad(DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p1[ 2 ]),
205  DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p1[ 2 ]),
206  DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p1[ 2 ]),
207  DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p1[ 2 ]) );
208  viewer.setFillColor(XZ_COLOR);
209  viewer.addQuad(DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p1[ 2 ]),
210  DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p1[ 2 ]),
211  DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p0[ 2 ]),
212  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p0[ 2 ]));
213  viewer.setFillColor(YZ_COLOR);
214  viewer.addQuad(DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p1[ 2 ]),
215  DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p1[ 2 ]),
216  DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p0[ 2 ]),
217  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p0[ 2 ]));
218  }
219 }
220 
221 template <typename KSpace, typename Naive3DDSSComputer, typename space, typename kspace >
222 void displayDSS3d( Viewer3D<space, kspace> & viewer,
223  const KSpace & ks, const Naive3DDSSComputer & dss3d,
224  const DGtal::Color & color3d )
225 {
226  viewer << CustomColors3D( color3d, color3d ) << dss3d;
227 }
228 
229 template <typename Point1, typename Point2>
230 void assign( Point1 & p1, const Point2 & p2 )
231 {
232  p1[ 0 ] = p2[ 0 ];
233  p1[ 1 ] = p2[ 1 ];
234  p1[ 2 ] = p2[ 2 ];
235 }
236 
237 template <typename KSpace, typename Naive3DDSSComputer, typename space, typename kspace >
238 void displayDSS3dTangent( Viewer3D<space, kspace> & viewer,
239  const KSpace & ks, const Naive3DDSSComputer & dss3d,
240  const DGtal::Color & color3d )
241 {
242  typedef typename Naive3DDSSComputer::Point3d Point;
243  typedef typename Naive3DDSSComputer::Integer Integer;
244  typedef typename Naive3DDSSComputer::PointR3d PointR3d;
245  typedef typename Display3D<>::BallD3D PointD3D;
246  Point directionZ3;
247  PointR3d interceptR, thicknessR;
248  Z3i::RealPoint P1, P2, direction, intercept, thickness;
249  dss3d.getParameters( directionZ3, interceptR, thicknessR );
250 
251  intercept[0] = (double) NumberTraits<Integer>::castToDouble ( interceptR[0].first ) / (double) NumberTraits<Integer>::castToDouble ( interceptR[0].second );
252  intercept[1] = (double) NumberTraits<Integer>::castToDouble ( interceptR[1].first ) / (double) NumberTraits<Integer>::castToDouble ( interceptR[1].second );
253  intercept[2] = (double) NumberTraits<Integer>::castToDouble ( interceptR[2].first ) / (double) NumberTraits<Integer>::castToDouble ( interceptR[2].second );
254  thickness[0] = (double) NumberTraits<Integer>::castToDouble ( thicknessR[0].first ) / (double) NumberTraits<Integer>::castToDouble ( thicknessR[0].second );
255  thickness[1] = (double) NumberTraits<Integer>::castToDouble ( thicknessR[1].first ) / (double) NumberTraits<Integer>::castToDouble ( thicknessR[1].second );
256  thickness[2] = (double) NumberTraits<Integer>::castToDouble ( thicknessR[2].first ) / (double) NumberTraits<Integer>::castToDouble ( thicknessR[2].second );
257 
258  assign( direction, directionZ3 );
259  direction /= direction.norm();
260  assign( P1, *dss3d.begin() );
261  assign( P2, *(dss3d.end()-1) );
262  double t1 = (P1 - intercept).dot( direction );
263  double t2 = (P2 - intercept).dot( direction );
264 
265  Z3i::RealPoint Q1 = intercept + t1 * direction;
266  Z3i::RealPoint Q2 = intercept + t2 * direction;
267  viewer.setLineColor(color3d);
268  viewer.addLine( Z3i::RealPoint(Q1[ 0 ]-0.5, Q1[ 1 ]-0.5, Q1[ 2 ]-0.5),
269  Z3i::RealPoint(Q2[ 0 ]-0.5, Q2[ 1 ]-0.5, Q2[ 2 ]-0.5),
270  MS3D_LINESIZE );
271 }
272 
273 template <typename KSpace, typename StandardDSS6Computer, typename space, typename kspace >
274 void displayProj2d6( Viewer3D<space, kspace> & viewer,
275  const KSpace & ks, const StandardDSS6Computer & dss3d,
276  const DGtal::Color & color2d )
277 {
278  typedef typename StandardDSS6Computer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
279  typedef typename ArithmeticalDSSComputer2d::ConstIterator ConstIterator2d;
280  typedef typename ArithmeticalDSSComputer2d::Point Point2d;
281  typedef typename KSpace::Cell Cell;
282  typedef typename KSpace::Point Point3d;
283  Point3d b = ks.lowerBound();
284  for ( DGtal::Dimension i = 0; i < 3; ++i )
285  {
286  const ArithmeticalDSSComputer2d & dss2d = dss3d.arithmeticalDSS2d( i );
287  for ( ConstIterator2d itP = dss2d.begin(), itPEnd = dss2d.end(); itP != itPEnd; ++itP )
288  {
289  Point2d p = *itP;
290  Point3d q;
291  switch (i)
292  {
293  case 0: q = Point3d( 2*b[ i ] , 2*p[ 0 ]+1, 2*p[ 1 ]+1 ); break;
294  case 1: q = Point3d( 2*p[ 0 ]+1, 2*b[ i ] , 2*p[ 1 ]+1 ); break;
295  case 2: q = Point3d( 2*p[ 0 ]+1, 2*p[ 1 ]+1, 2*b[ i ] ); break;
296  }
297  Cell c = ks.uCell( q );
298  viewer << CustomColors3D( color2d, color2d ) << c;
299  }
300  }
301 }
302 
303 template <typename KSpace, typename StandardDSS6Computer, typename space, typename kspace >
304 void displayDSS2d6( Viewer3D<space, kspace> & viewer,
305  const KSpace & ks, const StandardDSS6Computer & dss3d,
306  const DGtal::Color & color2d )
307 {
308  typedef typename StandardDSS6Computer::ConstIterator ConstIterator3d;
309  typedef typename StandardDSS6Computer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
310  typedef typename ArithmeticalDSSComputer2d::ConstIterator ConstIterator2d;
311  typedef typename ArithmeticalDSSComputer2d::Point Point2d;
312  typedef typename KSpace::Cell Cell;
313  typedef typename KSpace::Point Point3d;
314  typedef DGtal::PointVector<2,double> PointD2d;
315  typedef typename Display3D<>::BallD3D PointD3D;
316  Point3d b = ks.lowerBound();
317  for ( DGtal::Dimension i = 0; i < 3; ++i )
318  {
319  const typename ArithmeticalDSSComputer2d::Primitive & dss2d
320  = dss3d.arithmeticalDSS2d( i ).primitive();
321  // draw 2D bounding boxes for each arithmetical dss 2D.
322  std::vector<PointD2d> pts2d;
323  pts2d.push_back( dss2d.project(dss2d.back(), dss2d.Uf()) );
324  pts2d.push_back( dss2d.project(dss2d.back(), dss2d.Lf()) );
325  pts2d.push_back( dss2d.project(dss2d.front(), dss2d.Lf()) );
326  pts2d.push_back( dss2d.project(dss2d.front(), dss2d.Uf()) );
327  std::vector<PointD3D> bb;
328  PointD3D p3;
329  for ( unsigned int j = 0; j < pts2d.size(); ++j )
330  {
331  switch (i)
332  {
333  case 0: p3.center[0] = (double) b[ i ]-0.5; p3.center[1] = pts2d[ j ][ 0 ]; p3.center[2] = pts2d[ j ][ 1 ]; break;
334  case 1: p3.center[0] = pts2d[ j ][ 0 ]; p3.center[1] = (double) b[ i ]-0.5; p3.center[2] = pts2d[ j ][ 1 ]; break;
335  case 2: p3.center[0] = pts2d[ j ][ 0 ]; p3.center[1] = pts2d[ j ][ 1 ]; p3.center[2] = (double) b[ i ]-0.5; break;
336  }
337  bb.push_back( p3 );
338  }
339  for ( unsigned int j = 0; j < pts2d.size(); ++j ){
340  viewer.setLineColor(color2d);
341  viewer.addLine( DGtal::Z3i::RealPoint(bb[ j ].center[0], bb[ j ].center[1], bb[ j ].center[2]),
342  DGtal::Z3i::RealPoint(bb[ (j+1)%4 ].center[0], bb[ (j+1)%4 ].center[1], bb[ (j+1)%4 ].center[2]),
343  MS3D_LINESIZE );
344  }
345  } // for ( DGtal::Dimension i = 0; i < 3; ++i )
346 }
347 
348 
349 //Why not to just project 3D curve?
350 template <typename KSpace, typename Naive3DDSSComputer, typename space, typename kspace >
351 void displayProj2d26( Viewer3D<space, kspace> & viewer,
352  const KSpace & ks, const Naive3DDSSComputer & dss3d,
353  const DGtal::Color & color2d )
354 {
355  typedef typename Naive3DDSSComputer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
356  typedef typename ArithmeticalDSSComputer2d::ConstIterator ConstIterator2d;
357  typedef typename ArithmeticalDSSComputer2d::Point Point2d;
358  typedef typename KSpace::Cell Cell;
359  typedef typename KSpace::Point Point3d;
360  Point3d b = ks.lowerBound();
361 
362  const ArithmeticalDSSComputer2d & dssXY = dss3d.arithmeticalDSS2dXY();
363  const ArithmeticalDSSComputer2d & dssXZ = dss3d.arithmeticalDSS2dXZ();
364  const ArithmeticalDSSComputer2d & dssYZ = dss3d.arithmeticalDSS2dYZ();
365 
366  bool validXY = dss3d.validArithmeticalDSS2d ( 2 );
367  bool validXZ = dss3d.validArithmeticalDSS2d ( 1 );
368  bool validYZ = dss3d.validArithmeticalDSS2d ( 0 );
369 
370  if ( validXY && validXZ )
371  { //XY-plane, XZ-plane
372  for ( ConstIterator2d itXY = dssXY.begin(), itXZ = dssXZ.begin(), itPEnd = dssXY.end(); itXY != itPEnd; ++itXY, ++itXZ )
373  {
374  Point2d p1 = *itXY, p2 = *itXZ;
375  Point3d q1, q2, q3;
376  q1 = Point3d ( 2*b[ 0 ], 2*p1[ 1 ]+1, 2*p2[ 1 ]+1 );
377  q2 = Point3d ( 2*p2[ 0 ]+1, 2*b[ 1 ], 2*p2[ 1 ]+1 );
378  q3 = Point3d ( 2*p1[ 0 ]+1, 2*p1[ 1 ]+1, 2*b[ 2 ] );
379  Cell c1 = ks.uCell( q1 ); Cell c2 = ks.uCell( q2 ); Cell c3 = ks.uCell( q3 );
380  viewer << CustomColors3D( color2d, color2d ) << c1;
381  viewer << CustomColors3D( color2d, color2d ) << c2;
382  viewer << CustomColors3D( color2d, color2d ) << c3;
383  }
384  }
385  else
386  {
387  if ( validYZ && validXY )
388  { //XY-plane, YZ-plane
389  for ( ConstIterator2d itYZ = dssYZ.begin(), itXY = dssXY.begin(), itPEnd = dssXY.end(); itXY != itPEnd; ++itXY, ++itYZ )
390  {
391  Point2d p1 = *itYZ, p2 = *itXY;
392  Point3d q1, q2, q3;
393  q1 = Point3d ( 2*b[ 0 ], 2*p1[ 0 ]+1, 2*p1[ 1 ]+1 );
394  q2 = Point3d ( 2*p2[ 0 ]+1, 2*b[ 1 ], 2*p1[ 1 ]+1 );
395  q3 = Point3d ( 2*p2[ 0 ]+1, 2*p2[ 1 ]+1, 2*b[ 2 ] );
396  Cell c1 = ks.uCell( q1 ); Cell c2 = ks.uCell( q2 ); Cell c3 = ks.uCell( q3 );
397  viewer << CustomColors3D( color2d, color2d ) << c1;
398  viewer << CustomColors3D( color2d, color2d ) << c2;
399  viewer << CustomColors3D( color2d, color2d ) << c3;
400  }
401  }
402  else
403  {
404  for ( ConstIterator2d itYZ = dssYZ.begin(), itXZ = dssXZ.begin(), itPEnd = dssXZ.end(); itXZ != itPEnd; ++itXZ, ++itYZ )
405  {
406  Point2d p1 = *itYZ, p2 = *itXZ;
407  Point3d q1, q2, q3;
408  q1 = Point3d ( 2*b[ 0 ], 2*p1[ 0 ]+1, 2*p1[ 1 ]+1 );
409  q2 = Point3d ( 2*p2[ 0 ]+1, 2*b[ 1 ], 2*p2[ 1 ]+1 );
410  q3 = Point3d ( 2*p2[ 0 ]+1, 2*p1[ 0 ]+1, 2*b[ 2 ] );
411  Cell c1 = ks.uCell( q1 ); Cell c2 = ks.uCell( q2 );Cell c3 = ks.uCell( q3 );
412  viewer << CustomColors3D( color2d, color2d ) << c1;
413  viewer << CustomColors3D( color2d, color2d ) << c2;
414  viewer << CustomColors3D( color2d, color2d ) << c3;
415  }
416  }
417  }
418 }
419 
420 
421 template <typename KSpace, typename Naive3DDSSComputer, typename space, typename kspace >
422 void displayDSS2d26( Viewer3D<space, kspace> & viewer,
423  const KSpace & ks, const Naive3DDSSComputer & dss3d,
424  const DGtal::Color & color2d )
425 {
426  typedef typename Naive3DDSSComputer::ConstIterator ConstIterator3d;
427  typedef typename Naive3DDSSComputer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
428  typedef typename ArithmeticalDSSComputer2d::ConstIterator ConstIterator2d;
429  typedef typename ArithmeticalDSSComputer2d::Point Point2d;
430  typedef typename KSpace::Cell Cell;
431  typedef typename KSpace::Point Point3d;
432  typedef DGtal::PointVector<2,double> PointD2d;
433  typedef typename Display3D<>::BallD3D PointD3D;
434  Point3d b = ks.lowerBound();
435  for ( DGtal::Dimension i = 0; i < 3; ++i )
436  {
437  if ( !dss3d.validArithmeticalDSS2d ( i ) )
438  continue;
439 
440  const typename ArithmeticalDSSComputer2d::Primitive & dss2d
441  = dss3d.arithmeticalDSS2d( i ).primitive();
442  // draw 2D bounding boxes for each arithmetical dss 2D.
443  std::vector<PointD2d> pts2d;
444  pts2d.push_back( dss2d.project(dss2d.back(), dss2d.Uf()) );
445  pts2d.push_back( dss2d.project(dss2d.back(), dss2d.Lf()) );
446  pts2d.push_back( dss2d.project(dss2d.front(), dss2d.Lf()) );
447  pts2d.push_back( dss2d.project(dss2d.front(), dss2d.Uf()) );
448  std::vector<PointD3D> bb;
449  PointD3D p3;
450  for ( unsigned int j = 0; j < pts2d.size(); ++j )
451  {
452  switch (i)
453  {
454  case 0: p3.center[0] = (double) b[ i ]-0.5; p3.center[1] = pts2d[ j ][ 0 ]; p3.center[2] = pts2d[ j ][ 1 ]; break;
455  case 1: p3.center[0] = pts2d[ j ][ 0 ]; p3.center[1] = (double) b[ i ]-0.5; p3.center[2] = pts2d[ j ][ 1 ]; break;
456  case 2: p3.center[0] = pts2d[ j ][ 0 ]; p3.center[1] = pts2d[ j ][ 1 ]; p3.center[2] = (double) b[ i ]-0.5; break;
457  }
458  bb.push_back( p3 );
459  }
460  for ( unsigned int j = 0; j < pts2d.size(); ++j ){
461  viewer.setLineColor(color2d);
462  viewer.addLine( DGtal::Z3i::RealPoint(bb[ j ].center[0], bb[ j ].center[1], bb[ j ].center[2]),
463  DGtal::Z3i::RealPoint(bb[ (j+1)%4 ].center[0], bb[ (j+1)%4 ].center[1], bb[ (j+1)%4 ].center[2]),
464  MS3D_LINESIZE );
465  }
466  } // for ( DGtal::Dimension i = 0; i < 3; ++i )
467 }
468 
472 template <typename KSpace, typename PointIterator, typename space, typename kspace >
473 bool displayCover6( Viewer3D<space, kspace> & viewer,
474  const KSpace & ks, PointIterator b, PointIterator e,
475  bool dss3d, bool proj2d, bool dss2d, bool tangent,
476  int nbColors )
477 {
478  typedef typename PointIterator::value_type Point;
479  typedef StandardDSS6Computer<PointIterator,int,4> SegmentComputer;
480  typedef SaturatedSegmentation<SegmentComputer> Decomposition;
481  typedef typename Decomposition::SegmentComputerIterator SegmentComputerIterator;
482  typedef typename SegmentComputer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
483  SegmentComputer algo;
484  Decomposition theDecomposition(b, e, algo);
485 
486  viewer << SetMode3D( algo.className(), "BoundingBox" );
487  HueShadeColorMap<int> cmap_hue( 0, nbColors, 1 );
488 
489  unsigned int c = 0;
490  for ( SegmentComputerIterator i = theDecomposition.begin();
491  i != theDecomposition.end(); ++i)
492  {
493  SegmentComputer ms3d(*i);
494  const ArithmeticalDSSComputer2d & dssXY = ms3d.arithmeticalDSS2dXY();
495  const ArithmeticalDSSComputer2d & dssXZ = ms3d.arithmeticalDSS2dXZ();
496  const ArithmeticalDSSComputer2d & dssYZ = ms3d.arithmeticalDSS2dYZ();
497  Point f = *ms3d.begin();
498  Point l = *(ms3d.end() - 1);
499  trace.info() << "- " << c
500  << " MS3D,"
501  << " [" << f[ 0 ] << "," << f[ 1 ] << ","<< f[ 2 ] << "]"
502  << "->[" << l[ 0 ] << "," << l[ 1 ] << ","<< l[ 2 ] << "]"
503  << ", XY("
504  << dssXY.a() << "," << dssXY.b() << "," << dssXY.mu()
505  << "), XZ("
506  << dssXZ.a() << "," << dssXZ.b() << "," << dssXZ.mu()
507  << "), YZ("
508  << dssYZ.a() << "," << dssYZ.b() << "," << dssYZ.mu()
509  << ")" << std::endl;
510  Color color = cmap_hue( c );
511  if ( tangent ) displayDSS3dTangent( viewer, ks, ms3d, color );
512  if ( dss3d ) displayDSS3d( viewer, ks, ms3d, color );
513  if ( dss2d ) displayDSS2d6( viewer, ks, ms3d, color );
514  if ( proj2d ) displayProj2d6( viewer, ks, ms3d, CURVE2D_COLOR );
515  c++;
516  }
517  return true;
518 }
519 
524 template <typename KSpace, typename PointIterator, typename space, typename kspace >
525 bool displayCover26( Viewer3D<space, kspace> & viewer,
526  const KSpace & ks, PointIterator b, PointIterator e,
527  bool dss3d, bool proj2d, bool dss2d, bool tangent,
528  int nbColors )
529 {
530  typedef typename PointIterator::value_type Point;
531  typedef Naive3DDSSComputer<PointIterator,int, 8 > SegmentComputer;
532  typedef SaturatedSegmentation<SegmentComputer> Decomposition;
533  typedef typename Decomposition::SegmentComputerIterator SegmentComputerIterator;
534  typedef typename SegmentComputer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
535  SegmentComputer algo;
536  Decomposition theDecomposition(b, e, algo);
537 
538  viewer << SetMode3D( algo.className(), "BoundingBox" );
539  HueShadeColorMap<int> cmap_hue( 0, nbColors, 1 );
540 
541  unsigned int c = 0;
542  for ( SegmentComputerIterator i = theDecomposition.begin();
543  i != theDecomposition.end(); ++i)
544  {
545  SegmentComputer ms3d(*i);
546  const ArithmeticalDSSComputer2d & dssXY = ms3d.arithmeticalDSS2dXY();
547  const ArithmeticalDSSComputer2d & dssXZ = ms3d.arithmeticalDSS2dXZ();
548  const ArithmeticalDSSComputer2d & dssYZ = ms3d.arithmeticalDSS2dYZ();
549  Point f = *ms3d.begin();
550  Point l = *(ms3d.end() - 1);
551  trace.info() << "- " << c
552  << " MS3D,"
553  << " [" << f[ 0 ] << "," << f[ 1 ] << ","<< f[ 2 ] << "]"
554  << "->[" << l[ 0 ] << "," << l[ 1 ] << ","<< l[ 2 ] << "]"
555  << ", XY("
556  << dssXY.a() << "," << dssXY.b() << "," << dssXY.mu()
557  << "), XZ("
558  << dssXZ.a() << "," << dssXZ.b() << "," << dssXZ.mu()
559  << "), YZ("
560  << dssYZ.a() << "," << dssYZ.b() << "," << dssYZ.mu()
561  << ")" << std::endl;
562  //trace.info() << ms3d << std::endl; // information
563 
564  Color color = cmap_hue( c );
565  if ( tangent ) displayDSS3dTangent( viewer, ks, ms3d, color );
566  if ( dss3d ) displayDSS3d( viewer, ks, ms3d, color );
567  if ( dss2d ) displayDSS2d26( viewer, ks, ms3d, color );
568  if ( proj2d ) displayProj2d26( viewer, ks, ms3d, CURVE2D_COLOR );
569  c++;
570  }
571  return true;
572 }
573 
574 template < typename PointIterator, typename Space, typename TangentSequence >
575 void ComputeVCM ( const double & R, const double & r,
576  const PointIterator & begin, const PointIterator & end, TangentSequence & tangents, const std::string & output )
577 {
578  using namespace DGtal;
579  typedef ExactPredicateLpSeparableMetric < Space, 2 > Metric; // L2-metric
580  typedef VoronoiCovarianceMeasure < Space, Metric > VCM;
581  typedef HyperRectDomain < Space > Domain;
582  typedef functors::HatPointFunction < typename Space::Point, double > KernelFunction;
583  typedef EigenDecomposition < 3, double > LinearAlgebraTool;
584  typedef LinearAlgebraTool::Matrix Matrix;
585 
586  fstream outputStream;
587  outputStream.open ( ( output + ".vcm" ).c_str(), std::ios::out );
588  outputStream << "# VCM estimation R=" << R << " r=" << r << " chi=hat" << endl;
589 
590  Metric l2;
591  VCM vcm ( R, ceil( r ), l2, true );
592  vcm.init( begin, end );
593  KernelFunction chi( 1.0, r );
594  Matrix vcm_r, evec;
595  typename Space::RealVector eval;
596  for ( PointIterator it = begin; it != end; ++it )
597  {
598  // Compute VCM and diagonalize it.
599  vcm_r = vcm.measure( chi, *it );
600  LinearAlgebraTool::getEigenDecomposition ( vcm_r, evec, eval );
601  typename Space::RealVector tangent = evec.column( 0 );
602  tangents.push_back ( tangent );
603  outputStream << (*it)[0] << " " << (*it)[1] << " " << (*it)[2] << " "
604  << tangent[0] << " " << tangent[1] << " " << tangent[2] << endl;
605  }
606  outputStream.close();
607 }
608 
609 template < typename PointIterator, typename Space, typename TangentSequence >
610 void ComputeLMST6 ( const PointIterator & begin, const PointIterator & end, TangentSequence & tangents, const std::string & output )
611 {
612  typedef typename PointIterator::value_type Point;
613  typedef StandardDSS6Computer<PointIterator,int, 4 > SegmentComputer;
614  typedef SaturatedSegmentation<SegmentComputer> Decomposition;
615  typedef typename Decomposition::SegmentComputerIterator SegmentComputerIterator;
616  typedef typename SegmentComputer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
617  SegmentComputer algo;
618  Decomposition theDecomposition ( begin, end, algo );
619 
620  fstream outputStream;
621  outputStream.open ( ( output + ".lmst" ).c_str(), std::ios::out );
622  outputStream << "X Y Z X Y Z" << endl;
623 
624  LambdaMST3D < Decomposition > lmst64;
625  lmst64.attach ( theDecomposition );
626  lmst64.init ( begin, end );
627  lmst64.eval ( begin, end, std::back_inserter ( tangents ) );
628  typename TangentSequence::iterator itt = tangents.begin();
629  for ( PointIterator it = begin; it != end; ++it, ++itt )
630  {
631  typename Space::RealVector & tangent = (*itt);
632  if ( tangent.norm() != 0. )
633  tangent = tangent.getNormalized();
634  outputStream << (*it)[0] << " " << (*it)[1] << " " << (*it)[2] << " "
635  << tangent[0] << " " << tangent[1] << " " << tangent[2] << endl;
636  }
637 }
638 
639 template < typename PointIterator, typename Space, typename TangentSequence >
640 void ComputeLMST26 ( const PointIterator & begin, const PointIterator & end, TangentSequence & tangents, const std::string & output )
641 {
642  typedef typename PointIterator::value_type Point;
643  typedef Naive3DDSSComputer<PointIterator,int, 8 > SegmentComputer;
644  typedef SaturatedSegmentation<SegmentComputer> Decomposition;
645  typedef typename Decomposition::SegmentComputerIterator SegmentComputerIterator;
646  typedef typename SegmentComputer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
647  SegmentComputer algo;
648  Decomposition theDecomposition ( begin, end, algo );
649 
650  fstream outputStream;
651  outputStream.open ( ( output + ".lmst" ).c_str(), std::ios::out );
652  outputStream << "X Y Z X Y Z" << endl;
653 
654  LambdaMST3D < Decomposition > lmst64;
655  lmst64.attach ( theDecomposition );
656  lmst64.init ( begin, end );
657  lmst64.eval ( begin, end, std::back_inserter ( tangents ) );
658  typename TangentSequence::iterator itt = tangents.begin();
659  for ( PointIterator it = begin; it != end; ++it, ++itt )
660  {
661  typename Space::RealVector & tangent = (*itt);
662  if ( tangent.norm() != 0. )
663  tangent = tangent.getNormalized();
664  outputStream << (*it)[0] << " " << (*it)[1] << " " << (*it)[2] << " "
665  << tangent[0] << " " << tangent[1] << " " << tangent[2] << endl;
666  }
667 }
668 
670 namespace po = boost::program_options;
671 
680 int main(int argc, char **argv)
681 {
682  using namespace std;
683  using namespace DGtal;
684  typedef SpaceND<3,int> Z3;
685  typedef KhalimskySpaceND<3,int> K3;
686  typedef Z3::Point Point;
687  typedef Z3::RealPoint RealPoint;
688  typedef Z3::RealVector RealVector;
689  typedef HyperRectDomain<Z3> Domain;
690  typedef typename vector<Point>::const_iterator PointIterator;
691 
692  // specify command line ----------------------------------------------
693  QApplication application(argc,argv); // remove Qt arguments.
694  po::options_description general_opt("Specific allowed options are ");
695  general_opt.add_options()
696  ("help,h", "display this message")
697  ("input,i", po::value<string>(), "the name of the text file containing the list of 3D points: (x y z) per line." )
698  ("view,V", po::value<string>()->default_value( "OFF" ), "toggles display ON/OFF" )
699  ("box,b", po::value<int>()->default_value( 0 ), "specifies the the tightness of the bounding box around the curve with a given integer displacement <arg> to enlarge it (0 is tight)" )
700  ("viewBox,v", po::value<string>()->default_value( "WIRED" ), "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)." )
701  ("connectivity,T", po::value<string>()->default_value( "6" ), "specifies whether it is a 6-connected curve or a 26-connected curve: arg=6 | 26.")
702  ("curve3d,C", "displays the 3D curve")
703  ("curve2d,c", "displays the 2D projections of the 3D curve on the bounding box")
704  ("cover3d,3", "displays the 3D tangential cover of the curve" )
705  ("cover2d,2", "displays the 2D projections of the 3D tangential cover of the curve" )
706  ("nbColors,n", po::value<int>()->default_value( 3 ), "sets the number of successive colors used for displaying 2d and 3d maximal segments (default is 3: red, green, blue)" )
707  ("tangent,t", "displays the tangents to the curve" )
708  ("big-radius,R", po::value<double>()->default_value( 10.0 ), "the radius parameter R in the VCM estimator." )
709  ("small-radius,r", po::value<double>()->default_value( 3.0 ), "the radius parameter r in the VCM estimator." )
710  ("method,m", po::value<string>()->default_value( "VCM" ), "the method of tangent computation: VCM (default), L-MST." )
711  ("output,o", po::value<string>()->default_value( "3d-curve-tangent-estimations"), "the basename of the output text file which will contain points and tangent vectors: (x y z tx ty tz) per line" )
712  ;
713  po::positional_options_description pos_opt;
714  pos_opt.add("input", 1);
715 
716  // parse command line ----------------------------------------------
717  bool parseOK=true;
718  po::variables_map vm;
719  try
720  {
721  po::command_line_parser clp( argc, argv );
722  clp.options( general_opt ).positional( pos_opt );
723  po::store( clp.run(), vm );
724  }
725  catch( const exception& ex )
726  {
727  parseOK = false;
728  trace.info() << "Error checking program options: "<< ex.what() << endl;
729  }
730  po::notify( vm );
731  if( !parseOK || vm.count("help")||argc<=1)
732  {
733  cout << "Usage: " << argv[0] << " [options] --input <filename>\n"
734  << "This program estimates the tangent vector to a set of 3D integer points, which are supposed to approximate a 3D curve. "
735  << "This set of points is given as a list of points in file <input>."
736  << endl
737  << "The tangent estimator uses either the digital Voronoi Covariance Measure (VCM) or the 3D lambda-Maximal Segment Tangent (L-MST)."
738  << endl
739  << "This program can also displays the curve and tangent estimations, "
740  << "and it can also extract maximal digital straight segments (2D and 3D)."
741  << endl
742  << "Note: It is not compulsory for the points to be ordered in sequence, except if you wish to compute maximal digital straight segments."
743  << endl
744  << "In this case, you can select the connectivity of your curve between 6 (standard) or 26 (naive).\n"
745  << general_opt << "\n\n";
746  cout << "Example:\n"
747  << "3dCurveTangentEstimator -i ${DGtal}/examples/samples/sinus.dat -V ON -c -R 20 -r 3 -T 6\n";
748  return 0;
749  }
750 
751  // process command line ----------------------------------------------
752  string input = vm["input"].as<string>();
753  int b = vm["box"].as<int>();
754  // Create curve 3D.
755  vector<Point> sequence;
756  fstream inputStream;
757  inputStream.open ( input.c_str(), ios::in);
758  try
759  {
760  sequence = PointListReader<Point>::getPointsFromInputStream( inputStream );
761  if ( sequence.size() == 0) throw IOException();
762  }
763  catch (DGtal::IOException & ioe)
764  {
765  trace.error() << "Size is null." << endl;
766  }
767  inputStream.close();
768 
769  // start viewer
770  bool view = vm[ "view" ].as<string>() == "ON";
771 
772  // ----------------------------------------------------------------------
773  // Create domain and curve.
774  Point lowerBound = sequence[ 0 ];
775  Point upperBound = sequence[ 0 ];
776  for ( unsigned int j = 1; j < sequence.size(); ++j )
777  {
778  lowerBound = lowerBound.inf( sequence[ j ] );
779  upperBound = upperBound.sup( sequence[ j ] );
780  }
781  lowerBound -= Point::diagonal( b );
782  upperBound += Point::diagonal( b + 1 );
783  K3 ks; ks.init( lowerBound, upperBound, true );
784  Viewer3D<Z3,K3> viewer( ks );
785  trace.beginBlock ( "Tool 3dCurveTangentEstimator" );
786 
787  std::vector< RealVector > tangents;
788  string output = vm["output"].as<string>();
789  string method = vm["method"].as<string>();
790 
791  if ( method == "VCM" )
792  {
793  // input points of the curve are in sequence vector.
794  const double R = vm["big-radius"].as<double>();
795  trace.info() << "Big radius R = " << R << endl;
796  const double r = vm["small-radius"].as<double>();
797  trace.info() << "Small radius r = " << r << endl;
798 
799  ComputeVCM < PointIterator, Z3, std::vector< RealVector > > ( R, r, sequence.begin(), sequence.end(), tangents, output );
800  }
801  else if ( method == "L-MST" )
802  {
803  if (vm[ "connectivity" ].as<string>() == "6")
804  ComputeLMST6 < PointIterator, Z3, std::vector< RealVector > > ( sequence.begin(), sequence.end(), tangents, output );
805  else
806  ComputeLMST26 < PointIterator, Z3, std::vector< RealVector > > ( sequence.begin(), sequence.end(), tangents, output );
807  }
808  else
809  {
810  trace.info() << "Wrong method! Try: VCM or L-MST" << endl;
811  abort();
812  }
813 
814  if ( view )
815  {
816  for ( unsigned int i = 0; i < tangents.size(); i++ )
817  {
818  // Display normal
819  RealPoint p = sequence[i];
820  RealVector tangent = tangents[i];
821  viewer.setFillColor( Color ( 255, 0, 0, 255 ) );
822  viewer.setLineColor( Color ( 255, 0, 0, 255 ) );
823  viewer.addLine( p + 2.0 * tangent, p - 2.0 * tangent, 5.0 );
824  viewer.setFillColor( Color ( 100, 100, 140, 255 ) );
825  viewer.setLineColor( Color ( 100, 100, 140, 255 ) );
826  viewer.addBall( p, 0.125, 8 );
827  }
828  }
829 
830  GridCurve<K3> gc( ks );
831  try
832  {
833  gc.initFromPointsVector( sequence );
834  }
835  catch (DGtal::ConnectivityException& /*ce*/)
836  {
837  trace.warning() << "[ConnectivityException] GridCurve only accepts a sequence of face adjacent points. Try connectivity=6 instead." << endl;
838  }
839 
840  // ----------------------------------------------------------------------
841  // Displays everything.
842  bool res = true;
843  if ( view )
844  {
845  viewer.show();
846  // Display axes.
847  if ( vm.count( "viewBox" ) )
848  displayAxes<Point,RealPoint, Z3i::Space, Z3i::KSpace>( viewer, lowerBound, upperBound, vm[ "viewBox" ].as<std::string>() );
849  // Display 3D tangential cover.
850  res = vm[ "connectivity" ].as<string>() == "6"
851  ? displayCover6( viewer, ks, sequence.begin(), sequence.end(),
852  vm.count( "cover3d" ),
853  vm.count( "curve2d" ),
854  vm.count( "cover2d" ),
855  vm.count( "tangent" ),
856  vm["nbColors"].as<int>() )
857  : displayCover26( viewer, ks, sequence.begin(), sequence.end(),
858  vm.count( "cover3d" ),
859  vm.count( "curve2d" ),
860  vm.count( "cover2d" ),
861  vm.count( "tangent" ),
862  vm["nbColors"].as<int>() );
863  // Display 3D curve points.
864  if ( vm.count( "curve3d" ) )
865  {
866  viewer << CustomColors3D( CURVE3D_COLOR, CURVE3D_COLOR );
867  for ( vector<Point>::const_iterator it = sequence.begin(), itE = sequence.end();
868  it != itE; ++it )
869  viewer << *it;
870  }
871  // ----------------------------------------------------------------------
872  // User "interaction".
873  viewer << Viewer3D<Z3,K3>::updateDisplay;
874  application.exec();
875  }
876  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
877  trace.endBlock();
878 
879  return res ? 0 : 1;
880 }
std
STL namespace.
Companion project associated to DGtal ; Generated on Mon Jun 27 2016 12:21:49 for DGtalTools by   doxygen 1.8.11