DGtalTools  0.9.2
3dCurveViewer.cpp
1 
84 #include <iostream>
85 #include <iterator>
86 #include <cstdio>
87 #include <cmath>
88 #include <fstream>
89 #include <vector>
90 
91 #include <boost/program_options/options_description.hpp>
92 #include <boost/program_options/parsers.hpp>
93 #include <boost/program_options/variables_map.hpp>
94 
95 #include "DGtal/base/Common.h"
96 #include "DGtal/base/Exceptions.h"
97 #include "DGtal/kernel/SpaceND.h"
98 #include "DGtal/kernel/domains/HyperRectDomain.h"
99 #include "DGtal/topology/KhalimskySpaceND.h"
100 #include "DGtal/geometry/curves/GridCurve.h"
101 #include "DGtal/geometry/curves/StandardDSS6Computer.h"
102 #include "DGtal/geometry/curves/SaturatedSegmentation.h"
103 #include "DGtal/io/viewers/Viewer3D.h"
104 #include "DGtal/io/readers/PointListReader.h"
105 #include "DGtal/io/DrawWithDisplay3DModifier.h"
106 
107 
108 using namespace DGtal;
109 using namespace std;
110 
111 const Color AXIS_COLOR( 0, 0, 0, 255 );
112 const double AXIS_LINESIZE = 0.1;
113 const Color XY_COLOR( 0, 0, 255, 50 );
114 const Color XZ_COLOR( 0, 255, 0, 50 );
115 const Color YZ_COLOR( 255, 0, 0, 50 );
116 const Color CURVE3D_COLOR( 100, 100, 140, 128 );
117 const Color CURVE2D_COLOR( 200, 200, 200, 100 );
118 const double MS3D_LINESIZE = 0.05;
119 
121 // Functions for displaying the tangential cover of a 3D curve.
123 template <typename Point, typename RealPoint, typename space, typename kspace >
124 void displayAxes( Viewer3D<space, kspace> & viewer,
125  const Point & lowerBound, const Point & upperBound,
126  const std::string & mode )
127 {
128  RealPoint p0( (double)lowerBound[ 0 ]-0.5,
129  (double)lowerBound[ 1 ]-0.5,
130  (double)lowerBound[ 2 ]-0.5 );
131  RealPoint p1( (double)upperBound[ 0 ]-0.5,
132  (double)upperBound[ 1 ]-0.5,
133  (double)upperBound[ 2 ]-0.5 );
134  if ( ( mode == "WIRED" ) || ( mode == "COLORED" ) )
135  {
136  viewer.setLineColor(AXIS_COLOR);
137  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p0[ 2 ]),
138  DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p0[ 2 ]), AXIS_LINESIZE );
139  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p0[ 2 ]),
140  DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p0[ 2 ]), AXIS_LINESIZE );
141  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p0[ 2 ]),
142  DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
143  viewer.addLine( DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p0[ 2 ]),
144  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p0[ 2 ]), AXIS_LINESIZE );
145  viewer.addLine( DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p0[ 2 ]),
146  DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
147  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p0[ 2 ]),
148  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p0[ 2 ]), AXIS_LINESIZE );
149  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p0[ 2 ]),
150  DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
151  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p1[ 2 ]),
152  DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
153  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p1[ 2 ]),
154  DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
155  viewer.addLine( DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p0[ 2 ]),
156  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
157  viewer.addLine( DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p1[ 2 ]),
158  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
159  viewer.addLine( DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p1[ 2 ]),
160  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p1[ 2 ]), AXIS_LINESIZE );
161  }
162  if ( mode == "COLORED" )
163  {
164  viewer.setFillColor(XY_COLOR);
165  viewer.addQuad(DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p1[ 2 ]),
166  DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p1[ 2 ]),
167  DGtal::Z3i::RealPoint(p0[ 0 ], p0[ 1 ], p1[ 2 ]),
168  DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p1[ 2 ]) );
169  viewer.setFillColor(XZ_COLOR);
170  viewer.addQuad(DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p1[ 2 ]),
171  DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p1[ 2 ]),
172  DGtal::Z3i::RealPoint(p0[ 0 ], p1[ 1 ], p0[ 2 ]),
173  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p0[ 2 ]));
174  viewer.setFillColor(YZ_COLOR);
175  viewer.addQuad(DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p1[ 2 ]),
176  DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p1[ 2 ]),
177  DGtal::Z3i::RealPoint(p1[ 0 ], p0[ 1 ], p0[ 2 ]),
178  DGtal::Z3i::RealPoint(p1[ 0 ], p1[ 1 ], p0[ 2 ]));
179  }
180 }
181 
182 template <typename KSpace, typename StandardDSS6Computer, typename space, typename kspace >
183 void displayDSS3d( Viewer3D<space, kspace> & viewer,
184  const KSpace & ks, const StandardDSS6Computer & dss3d,
185  const DGtal::Color & color3d )
186 {
187  viewer << CustomColors3D( color3d, color3d ) << dss3d;
188 }
189 
190 template <typename Point1, typename Point2>
191 void assign( Point1 & p1, const Point2 & p2 )
192 {
193  p1[ 0 ] = p2[ 0 ];
194  p1[ 1 ] = p2[ 1 ];
195  p1[ 2 ] = p2[ 2 ];
196 }
197 
198 template <typename KSpace, typename StandardDSS6Computer, typename space, typename kspace >
199 void displayDSS3dTangent( Viewer3D<space, kspace> & viewer,
200  const KSpace & ks, const StandardDSS6Computer & dss3d,
201  const DGtal::Color & color3d )
202 {
203  typedef typename StandardDSS6Computer::Point3d Point;
204  typedef typename StandardDSS6Computer::PointR3d PointR3d;
205  typedef DGtal::PointVector<3,double> PointD3d;
206  typedef typename Display3D<>::BallD3D PointD3D;
207  Point directionZ3;
208  PointR3d interceptR, thicknessR;
209  PointD3d P1, P2, direction;
210  dss3d.getParameters( directionZ3, interceptR, thicknessR );
211 
212  PointD3d intercept;
213  intercept[0] = (double) NumberTraits<int>::castToInt64_t ( interceptR[0].first ) / (double) NumberTraits<int>::castToInt64_t ( interceptR[0].second );
214  intercept[1] = (double) NumberTraits<int>::castToInt64_t ( interceptR[1].first ) / (double) NumberTraits<int>::castToInt64_t ( interceptR[1].second );
215  intercept[2] = (double) NumberTraits<int>::castToInt64_t ( interceptR[2].first ) / (double) NumberTraits<int>::castToInt64_t ( interceptR[2].second );
216 
217  PointD3d thickness;
218  thickness[0] = (double) NumberTraits<int>::castToInt64_t ( thicknessR[0].first ) / (double) NumberTraits<int>::castToInt64_t ( thicknessR[0].second );
219  thickness[1] = (double) NumberTraits<int>::castToInt64_t ( thicknessR[1].first ) / (double) NumberTraits<int>::castToInt64_t ( thicknessR[1].second );
220  thickness[2] = (double) NumberTraits<int>::castToInt64_t ( thicknessR[2].first ) / (double) NumberTraits<int>::castToInt64_t ( thicknessR[2].second );
221 
222  assign( direction, directionZ3 );
223  direction /= direction.norm();
224  assign( P1, *dss3d.begin() );
225  assign( P2, *(dss3d.end()-1) );
226  double t1 = (P1 - intercept).dot( direction );
227  double t2 = (P2 - intercept).dot( direction );
228 
229  PointD3d Q1 = intercept + t1 * direction;
230  PointD3d Q2 = intercept + t2 * direction;
231  viewer.setLineColor(color3d);
232  viewer.addLine( DGtal::Z3i::RealPoint(Q1[ 0 ]-0.5, Q1[ 1 ]-0.5, Q1[ 2 ]-0.5),
233  DGtal::Z3i::RealPoint(Q2[ 0 ]-0.5, Q2[ 1 ]-0.5, Q2[ 2 ]-0.5),
234  MS3D_LINESIZE );
235 }
236 
237 template <typename KSpace, typename StandardDSS6Computer, typename space, typename kspace >
238 void displayProj2d( Viewer3D<space, kspace> & viewer,
239  const KSpace & ks, const StandardDSS6Computer & dss3d,
240  const DGtal::Color & color2d )
241 {
242  typedef typename StandardDSS6Computer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
243  typedef typename ArithmeticalDSSComputer2d::ConstIterator ConstIterator2d;
244  typedef typename ArithmeticalDSSComputer2d::Point Point2d;
245  typedef typename KSpace::Cell Cell;
246  typedef typename KSpace::Point Point3d;
247  Point3d b = ks.lowerBound();
248  for ( DGtal::Dimension i = 0; i < 3; ++i )
249  {
250  const ArithmeticalDSSComputer2d & dss2d = dss3d.arithmeticalDSS2d( i );
251  for ( ConstIterator2d itP = dss2d.begin(), itPEnd = dss2d.end(); itP != itPEnd; ++itP )
252  {
253  Point2d p = *itP;
254  Point3d q;
255  switch (i) {
256  case 0: q = Point3d( 2*b[ i ] , 2*p[ 0 ]+1, 2*p[ 1 ]+1 ); break;
257  case 1: q = Point3d( 2*p[ 0 ]+1, 2*b[ i ] , 2*p[ 1 ]+1 ); break;
258  case 2: q = Point3d( 2*p[ 0 ]+1, 2*p[ 1 ]+1, 2*b[ i ] ); break;
259  }
260  Cell c = ks.uCell( q );
261  viewer << CustomColors3D( color2d, color2d ) << c;
262  }
263  }
264 }
265 
266 template <typename KSpace, typename StandardDSS6Computer, typename space, typename kspace >
267 void displayDSS2d( Viewer3D<space, kspace> & viewer,
268  const KSpace & ks, const StandardDSS6Computer & dss3d,
269  const DGtal::Color & color2d )
270 {
271  typedef typename StandardDSS6Computer::ConstIterator ConstIterator3d;
272  typedef typename StandardDSS6Computer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
273  typedef typename ArithmeticalDSSComputer2d::ConstIterator ConstIterator2d;
274  typedef typename ArithmeticalDSSComputer2d::Point Point2d;
275  typedef typename KSpace::Cell Cell;
276  typedef typename KSpace::Point Point3d;
277  typedef DGtal::PointVector<2,double> PointD2d;
278  typedef typename Display3D<>::BallD3D PointD3D;
279  Point3d b = ks.lowerBound();
280  for ( DGtal::Dimension i = 0; i < 3; ++i )
281  {
282  const typename ArithmeticalDSSComputer2d::Primitive & dss2d
283  = dss3d.arithmeticalDSS2d( i ).primitive();
284  // draw 2D bounding boxes for each arithmetical dss 2D.
285  std::vector<PointD2d> pts2d;
286  pts2d.push_back( dss2d.project(dss2d.back(), dss2d.Uf()) );
287  pts2d.push_back( dss2d.project(dss2d.back(), dss2d.Lf()) );
288  pts2d.push_back( dss2d.project(dss2d.front(), dss2d.Lf()) );
289  pts2d.push_back( dss2d.project(dss2d.front(), dss2d.Uf()) );
290  std::vector<PointD3D> bb;
291  PointD3D p3;
292  for ( unsigned int j = 0; j < pts2d.size(); ++j )
293  {
294  switch (i) {
295  case 0: p3.center[0] = (double) b[ i ]-0.5; p3.center[1] = pts2d[ j ][ 0 ]; p3.center[2] = pts2d[ j ][ 1 ]; break;
296  case 1: p3.center[0] = pts2d[ j ][ 0 ]; p3.center[1] = (double) b[ i ]-0.5; p3.center[2] = pts2d[ j ][ 1 ]; break;
297  case 2: p3.center[0] = pts2d[ j ][ 0 ]; p3.center[1] = pts2d[ j ][ 1 ]; p3.center[2] = (double) b[ i ]-0.5; break;
298  }
299  bb.push_back( p3 );
300  }
301  for ( unsigned int j = 0; j < pts2d.size(); ++j ){
302  viewer.setLineColor(color2d);
303  viewer.addLine( DGtal::Z3i::RealPoint(bb[ j ].center[0], bb[ j ].center[1], bb[ j ].center[2]),
304  DGtal::Z3i::RealPoint(bb[ (j+1)%4 ].center[0], bb[ (j+1)%4 ].center[1], bb[ (j+1)%4 ].center[2]),
305  MS3D_LINESIZE );
306  }
307  } // for ( DGtal::Dimension i = 0; i < 3; ++i )
308 }
309 
314 template <typename KSpace, typename PointIterator, typename space, typename kspace >
315 bool displayCover( Viewer3D<space, kspace> & viewer,
316  const KSpace & ks, PointIterator b, PointIterator e,
317  bool dss3d, bool proj2d, bool dss2d, bool tangent,
318  int nbColors )
319 {
320  typedef typename PointIterator::value_type Point;
321  typedef StandardDSS6Computer<PointIterator,int,4> SegmentComputer;
322  typedef SaturatedSegmentation<SegmentComputer> Decomposition;
323  typedef typename Decomposition::SegmentComputerIterator SegmentComputerIterator;
324  typedef typename SegmentComputer::ArithmeticalDSSComputer2d ArithmeticalDSSComputer2d;
325  SegmentComputer algo;
326  Decomposition theDecomposition(b, e, algo);
327 
328  viewer << SetMode3D( algo.className(), "BoundingBox" );
329  HueShadeColorMap<int> cmap_hue( 0, nbColors, 1 );
330 
331  unsigned int c = 0;
332  for ( SegmentComputerIterator i = theDecomposition.begin();
333  i != theDecomposition.end(); ++i)
334  {
335  SegmentComputer ms3d(*i);
336  const ArithmeticalDSSComputer2d & dssXY = ms3d.arithmeticalDSS2dXY();
337  const ArithmeticalDSSComputer2d & dssXZ = ms3d.arithmeticalDSS2dXZ();
338  const ArithmeticalDSSComputer2d & dssYZ = ms3d.arithmeticalDSS2dYZ();
339  Point f = *ms3d.begin();
340  Point l = *(ms3d.end() - 1);
341  trace.info() << "- " << c
342  << " MS3D,"
343  << " [" << f[ 0 ] << "," << f[ 1 ] << ","<< f[ 2 ] << "]"
344  << "->[" << l[ 0 ] << "," << l[ 1 ] << ","<< l[ 2 ] << "]"
345  << ", XY("
346  << dssXY.a() << "," << dssXY.b() << "," << dssXY.mu()
347  << "), XZ("
348  << dssXZ.a() << "," << dssXZ.b() << "," << dssXZ.mu()
349  << "), YZ("
350  << dssYZ.a() << "," << dssYZ.b() << "," << dssYZ.mu()
351  << ")" << std::endl;
352  //trace.info() << ms3d << std::endl; // information
353 
354  Color color = cmap_hue( c );
355  if ( tangent ) displayDSS3dTangent( viewer, ks, ms3d, color );
356  if ( dss3d ) displayDSS3d( viewer, ks, ms3d, color );
357  if ( dss2d ) displayDSS2d( viewer, ks, ms3d, color );
358  if ( proj2d ) displayProj2d( viewer, ks, ms3d, CURVE2D_COLOR );
359  c++;
360  }
361  return true;
362 }
363 
364 
366 namespace po = boost::program_options;
367 
376 int main(int argc, char **argv)
377 {
378  typedef SpaceND<3,int> Z3;
379  typedef KhalimskySpaceND<3,int> K3;
380  typedef Z3::Point Point;
381  typedef Z3::RealPoint RealPoint;
382 
383  // specify command line ----------------------------------------------
384  QApplication application(argc,argv); // remove Qt arguments.
385  po::options_description general_opt("Specific allowed options (for Qt options, see Qt official site) are");
386  general_opt.add_options()
387  ("help,h", "display this message")
388  ("input,i", po::value<std::string>(), "the name of the text file containing the list of 3D points (x y z per line)" )
389  ("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)" )
390  ("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)." )
391  ("curve3d,C", "displays the 3D curve")
392  ("curve2d,c", "displays the 2D projections of the 3D curve on the bounding box")
393  ("cover3d,3", "displays the 3D tangential cover of the curve" )
394  ("cover2d,2", "displays the 2D projections of the 3D tangential cover of the curve" )
395  ("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)" )
396  ("tangent,t", "displays the tangents to the curve" )
397  ;
398  po::positional_options_description pos_opt;
399  pos_opt.add("input", 1);
400 
401  // parse command line ----------------------------------------------
402  bool parseOK=true;
403  po::variables_map vm;
404  try {
405  po::command_line_parser clp( argc, argv );
406  clp.options( general_opt ).positional( pos_opt );
407  po::store( clp.run(), vm );
408  } catch( const std::exception& ex ) {
409  parseOK = false;
410  trace.info() << "Error checking program options: "<< ex.what() << endl;
411  }
412  po::notify( vm );
413  if( !parseOK || vm.count("help")||argc<=1)
414  {
415  std::cout << "Usage: " << argv[0] << " [options] input\n"
416  << "Display a 3D curve given as the <input> filename (with possibly projections and/or tangent information) by using QGLviewer.\n"
417  << general_opt << "\n\n";
418  std::cout << "Example:\n"
419  << "3dCurveViewer -C -b 1 -3 -2 -c ${DGtal}/examples/samples/sinus.dat\n";
420  return 0;
421  }
422 
423  // process command line ----------------------------------------------
424  string input = vm["input"].as<std::string>();
425  int b = vm["box"].as<int>();
426  // Create curve 3D.
427  vector<Point> sequence;
428  fstream inputStream;
429  inputStream.open ( input.c_str(), ios::in);
430  try {
431  sequence = PointListReader<Point>::getPointsFromInputStream( inputStream );
432  if ( sequence.size() == 0) throw IOException();
433  }
434  catch (DGtal::IOException & ioe) {
435  trace.error() << "Size is null." << std::endl;
436  }
437  inputStream.close();
438 
439  // start viewer
440  Viewer3D<> viewer;
441  trace.beginBlock ( "Tool 3dCurveViewer" );
442 
443  // ----------------------------------------------------------------------
444  // Create domain and curve.
445  Point lowerBound = sequence[ 0 ];
446  Point upperBound = sequence[ 0 ];
447  for ( unsigned int j = 1; j < sequence.size(); ++j )
448  {
449  lowerBound = lowerBound.inf( sequence[ j ] );
450  upperBound = upperBound.sup( sequence[ j ] );
451  }
452  lowerBound -= Point::diagonal( b );
453  upperBound += Point::diagonal( b+1 );
454  K3 ks; ks.init( lowerBound, upperBound, true );
455  GridCurve<K3> gc( ks );
456  try {
457  gc.initFromPointsVector( sequence );
458  } catch (DGtal::ConnectivityException& /*ce*/) {
459  throw ConnectivityException();
460  return false;
461  }
462 
463  // ----------------------------------------------------------------------
464  // Displays everything.
465  viewer.show();
466  // Display axes.
467  if ( vm.count( "viewBox" ) )
468  displayAxes<Point,RealPoint, Z3i::Space, Z3i::KSpace>( viewer, lowerBound, upperBound, vm[ "viewBox" ].as<std::string>() );
469  // Display 3D tangential cover.
470  bool res = displayCover( viewer, ks, sequence.begin(), sequence.end(),
471  vm.count( "cover3d" ),
472  vm.count( "curve2d" ),
473  vm.count( "cover2d" ),
474  vm.count( "tangent" ),
475  vm["nbColors"].as<int>() );
476  // Display 3D curve points.
477  if ( vm.count( "curve3d" ) )
478  viewer << CustomColors3D( CURVE3D_COLOR, CURVE3D_COLOR )
479  << gc.getPointsRange()
480  << sequence.back(); // curiously, last point is not displayed.
481 
482  // ----------------------------------------------------------------------
483  // User "interaction".
484  viewer << Viewer3D<>::updateDisplay;
485  application.exec();
486  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
487  trace.endBlock();
488 
489  return res ? 0 : 1;
490 }
std
STL namespace.
Companion project associated to DGtal ; Generated on Mon Jun 27 2016 12:21:49 for DGtalTools by   doxygen 1.8.11