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