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