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