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