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