DGtalTools 2.0.0
Loading...
Searching...
No Matches
patternTriangulation.cpp
1
39#include <iostream>
40#include <iterator>
41#include <cstdio>
42#include <cmath>
43#include <fstream>
44#include <vector>
45
46#include "DGtal/base/Common.h"
47#include "DGtal/helpers/StdDefs.h"
48
49#include "DGtal/arithmetic/LighterSternBrocot.h"
50#include "DGtal/arithmetic/Pattern.h"
51
52#include "DGtal/io/boards/Board2D.h"
53#include "DGtal/io/Color.h"
54
55#include "CLI11.hpp"
56
57
58using namespace DGtal;
59using namespace std;
60
61
62
111template <typename Point>
112struct OddConvexHullMap {
113
114public:
115 typedef Point Vector;
116
117public:
124 OddConvexHullMap(const Vector& aShift = Vector(1,-1))
125 : myS(aShift) {}
126
133 Point operator()(const Point& aP) const
134 {
135 return aP + myS;
136 }
137
138private:
139 Vector myS;
140
141};
142
150template <typename Point>
151struct EvenConvexHullMap {
152
153public:
154 typedef Point Vector;
155
156public:
164 EvenConvexHullMap(const Point& aPoint, const Vector& aShift = Vector(1,-1))
165 : myZn(aPoint), myS(aShift) {}
166
173 Point operator()(const Point& aP) const
174 {
175 return myZn - aP + myS;
176 }
177
178private:
179 Point myZn;
180 Vector myS;
181
182};
183
185
196template <typename Board, typename Point>
197void drawSegment(Board& aBoard, const Point& aP, const Point& aQ)
198{
199 aBoard.drawLine(aP[0], aP[1], aQ[0], aQ[1]);
200 aBoard << aP << aQ;
201}
202
217template <typename Board, typename Iterator, typename Functor>
218void drawSegments(Board& aBoard,
219 const Iterator& itb, const Iterator& ite,
220 const Functor& aF)
221{
222 Iterator it = itb;
223 if (it != ite)
224 {
225 Iterator prevIt = it;
226 ++it;
227 for ( ; it != ite; ++prevIt, ++it)
228 {
229 drawSegment( aBoard, aF(*prevIt), aF(*it) );
230 }
231 }
232}
233
234
236
251template<typename Fraction, typename Container>
252void fillConvergents(const Fraction& aZn,
253 Container& odd, Container& even)
254{
255 typedef typename Container::value_type Vector;
256
257 for (typename Fraction::Quotient i = 1; i <= aZn.k(); ++i)
258 {
259 Fraction zk = aZn.reduced(i);
260 if (((aZn.k() - i)%2) == 1 )
261 { //odd
262 odd.push_back(Vector(zk.q(),zk.p()));
263 }
264 else
265 { //even
266 even.push_back(Vector(zk.q(),zk.p()));
267 }
268 }
269}
270
285template <typename Pattern, typename Integer>
286void displayConvexHull(const Pattern& aP, const Integer& aD, bool withFDT = false)
287{
288 typedef typename Pattern::Fraction Fraction;
289 typedef typename Pattern::Vector2I Vector;
290 typedef std::vector<Vector> Convergents;
291 typedef typename Pattern::Point2I Point;
292
293 Board2D aBoard;
294
295 //list of convergents
296 Convergents oddConvergents, evenConvergents;
297
298 //fill the lists
299 Fraction zn = aP.slope();
300 // aD > 1
301 if (aD > 1)
302 {
303 Fraction znm1 = zn.father();
304 if (zn.odd())
305 oddConvergents.push_back(Vector(znm1.q(),znm1.p()));
306 else
307 evenConvergents.push_back(Vector(znm1.q(),znm1.p()));
308 }
309 // aD >= 1
310 fillConvergents( zn, oddConvergents, evenConvergents );
311
312 //display the domain
313 Point Zn = Point(aD*zn.q(),aD*zn.p());
314 aBoard << Z2i::Domain( Point(0,0), Zn );
315 aBoard << SetMode( Zn.className(), "Grid" );
316 aBoard.setLineStyle(Board2D::Shape::SolidStyle);
317 aBoard.setPenColor(DGtal::Color::Black);
318
319 //display the two main lists
320 OddConvexHullMap<Point> oddH;
321 drawSegments( aBoard, oddConvergents.begin(), oddConvergents.end(), oddH );
322
323 EvenConvexHullMap<Point> evenH(Zn);
324 drawSegments( aBoard, evenConvergents.begin(), evenConvergents.end(), evenH );
325
326 //display the four last segments
327 drawSegment( aBoard, Point(0,0), Zn );
328 if (evenConvergents.size() != 0)
329 {
330 drawSegment( aBoard, evenH(*evenConvergents.rbegin()), Zn );
331 if (oddConvergents.size() != 0)
332 {
333 drawSegment( aBoard, Point(0,0), oddH(*oddConvergents.rbegin()) );
334 drawSegment( aBoard,
335 oddH(*oddConvergents.begin()),
336 evenH(*evenConvergents.begin()) );
337 }
338 else
339 {
340 drawSegment( aBoard, Point(0,0), evenH(*evenConvergents.rbegin()) );
341 }
342 }
343
344 if (withFDT)
345 {//display internal edges
346
347 //first internal edge
348 if (evenConvergents.size() != 0)
349 drawSegment( aBoard,
350 Point(0,0),
351 evenH(*evenConvergents.rbegin()) );
352 //other internal edges
353 typedef typename Convergents::const_reverse_iterator RI;
354 RI oddRit = oddConvergents.rbegin();
355 RI evenRit = evenConvergents.rbegin();
356 bool hasToContinue = true;
357 while (hasToContinue)
358 {
359 if (oddRit != oddConvergents.rend())
360 {
361 drawSegment( aBoard,
362 oddH(*oddRit),
363 evenH(*evenRit) );
364 }
365 else
366 hasToContinue = false;
367 ++evenRit;
368
369 if ( (hasToContinue) && (evenRit != evenConvergents.rend()) )
370 {
371 drawSegment( aBoard,
372 oddH(*oddRit),
373 evenH(*evenRit) );
374 }
375 else
376 hasToContinue = false;
377 ++oddRit;
378 }
379
380 aBoard.saveEPS("FDT.eps");
381 }
382 else
383 {
384 aBoard.saveEPS("CH.eps");
385 }
386}
387
389
400template <typename Board, typename Point>
401void drawTriangle(Board& aBoard,
402 const Point& aP, const Point& aQ, const Point& aR)
403{
404 aBoard.drawTriangle( aP[0], aP[1], aQ[0], aQ[1], aR[0], aR[1] );
405}
406
416struct InDirectBezoutComputer
417{
425 template<typename Pattern>
426 typename Pattern::Vector2I
427 getPositiveBezout(const Pattern& aP) const
428 {
429 return aP.bezout();
430 }
431
439 template<typename Pattern>
440 typename Pattern::Vector2I
441 getNegativeBezout(const Pattern& aP) const
442 {
443 return aP.slope().even()
444 ? aP.U( 1 ) - aP.previousPattern().U( 1 )
445 : aP.previousPattern().U( 1 );
446 }
447
448};
449
459struct DirectBezoutComputer
460{
468 template<typename Pattern>
469 typename Pattern::Vector2I
470 getPositiveBezout(const Pattern& aP) const
471 {
472 return aP.slope().even()
473 ? aP.U( 1 ) - aP.previousPattern().U( 1 )
474 : aP.previousPattern().U( 1 );
475 }
476
484 template<typename Pattern>
485 typename Pattern::Vector2I
486 getNegativeBezout(const Pattern& aP) const
487 {
488 return aP.bezout();
489 }
490
491};
492
511template <typename Board, typename Pattern, typename BezoutComputer>
512void displayPartialCDT(Board& aBoard,
513 const Pattern& aP,
514 const typename Pattern::Point2I& aStartingPoint,
515 const BezoutComputer& aBC )
516{
517 typedef typename Pattern::Fraction Fraction;
518 typedef typename Pattern::Vector2I Vector;
519 typedef typename Pattern::Point2I Point;
520
521 Fraction f = aP.slope();
522 if ( (f.p() >= 1)&&(f.q() >= 1) )
523 {
524 Point O = aStartingPoint;
525 Point Zn = aStartingPoint + aP.v();
526 Vector v1 = aBC.getNegativeBezout(aP);
527 Point B = aStartingPoint + v1;
528 drawTriangle(aBoard, O, Zn, B);
529
530 if ( (f.p() > 1)||(f.q() > 1) )
531 {
532 //recursive calls
533 // - first pattern
534 displayPartialCDT(aBoard, Pattern(Fraction(v1[1],v1[0])), O, aBC);
535 // - second pattern
536 Vector v2 = aBC.getPositiveBezout(aP);
537 displayPartialCDT(aBoard, Pattern(Fraction(v2[1],v2[0])), B, aBC);
538 }
539 }
540}
541
553template <typename Pattern, typename Integer>
554void displayCDT(const Pattern& aP, const Integer& aD)
555{
556 Board2D aBoard;
557
558 typedef typename Pattern::Point2I Point;
559 typedef typename Pattern::Fraction Fraction;
560 Fraction zn = aP.slope();
561
562 //display the domain
563 Point Zn = Point(aD*zn.q(),aD*zn.p());
564 aBoard << Z2i::Domain( Point(0,0), Zn );
565 aBoard << SetMode( Zn.className(), "Grid" );
566 aBoard.setLineStyle(Board2D::Shape::SolidStyle);
567 aBoard.setPenColor(DGtal::Color::Black);
568
569 //display the upper part of the triangulation
570 for (Integer i = 0; i < aD; ++i)
571 {
572 displayPartialCDT( aBoard, aP, Point(i*zn.q(), i*zn.p()), InDirectBezoutComputer() );
573 }
574
575 //display the lower parts of the triangulation
576 // - getting the reversed patterns
577 typedef typename Pattern::Vector2I Vector;
578 typedef std::vector<Vector> Convergents;
579 typedef typename std::vector<Vector>::const_reverse_iterator ReverseIterator;
580 Convergents oddConvergents, evenConvergents;
581 fillConvergents( zn, oddConvergents, evenConvergents );
582 // - displaying the reversed patterns...
583 Point rPStartingPoint;
584 // - of odd slope:
585 OddConvexHullMap<Point> oddH;
586 {
587 ReverseIterator itb = evenConvergents.rbegin();
588 ReverseIterator ite = evenConvergents.rend();
589 Point aStartingPoint = oddH( *oddConvergents.rbegin() );
590 //
591 Point p = aStartingPoint;
592 ReverseIterator it = itb;
593 if (it != ite)
594 {
595 for (++it; it != ite; ++it)
596 {
597 displayPartialCDT( aBoard, Pattern( (*it)[1], (*it)[0] ), p, DirectBezoutComputer() );
598 p += *it;
599 }
600 }
601 rPStartingPoint = p;
602 }
603 // - of even slope:
604 EvenConvexHullMap<Point> evenH(Zn);
605 {
606 ReverseIterator itb = oddConvergents.rbegin();
607 ReverseIterator ite = oddConvergents.rend();
608 Point aStartingPoint = evenH( *evenConvergents.rbegin() );
609 //
610 Point p = aStartingPoint;
611 ReverseIterator it = itb;
612 for ( ; it != ite; ++it)
613 {
614 p -= *it;
615 displayPartialCDT( aBoard, Pattern( (*it)[1], (*it)[0] ), p, DirectBezoutComputer() );
616 }
617 }
618 // - of same slope:
619 {
620 for (Integer i = 0; i < (aD-1); ++i)
621 {
622 Point p = rPStartingPoint + i*aP.v();
623 displayPartialCDT( aBoard, aP, p, DirectBezoutComputer() );
624 }
625 }
626
627
628 aBoard.saveEPS("CDT.eps");
629}
630
631
633
638void missingParam(std::string param)
639{
640 trace.error() << " Parameter: " << param << " is required...";
641 trace.info() << std::endl;
642 exit(1);
643}
644
645
647
655int main(int argc, char **argv)
656{
657 typedef DGtal::int32_t Integer;
658
659 Integer a;
660 Integer b;
661 Integer d {1};
662 string type {"CH"};
663
664
665 // parse command line using CLI ----------------------------------------------
666 CLI::App app;
667 app.description("Draw the Delaunay triangulation of a pattern using DGtal library. Example: patternTriangulation -a 5 -b 8 ");
668 app.add_option("-a,--aparam",a,"pattern a parameter")
669 ->required();
670 app.add_option("-b,--bparam",b,"pattern b parameter")
671 ->required();
672 app.add_option("-d,--delta",d,"number of repetitions");
673 app.add_option("--triangulation,-t",type,"output:\n\tClosest-point Delaunay triangulation {CDT}\n\tFarthest-point Delaunay triangulation {FDT}\n\tConvex hull {CH}");
674
675
676 app.get_formatter()->column_width(40);
677 CLI11_PARSE(app, argc, argv);
678 // END parse command line using CLI ----------------------------------------------
679
680
681
682
683 if ( (a < 0) || (b <= 0) )
684 {
685 trace.error() << " parameters a and b should be strictly positive.";
686 trace.info() << std::endl;
687 return 0;
688 }
689 if (a >= b)
690 {
691 trace.error() << " parameter a should be strictly less than b.";
692 trace.info() << std::endl;
693 return 0;
694 }
695
696 if (d <= 0)
697 {
698 trace.error() << " parameter d should be strictly positive";
699 trace.info() << std::endl;
700 return 0;
701 }
702 trace.info() << "a=" << a << ", b=" << b << ", d=" << d << std::endl;
703
704 typedef DGtal::int32_t Quotient;
705 typedef LighterSternBrocot<Integer, Quotient, StdMapRebinder> SB;
706 typedef SB::Fraction Fraction; // the type for fractions
707 typedef Pattern<Fraction> Pattern; // the type for patterns
708
709 Pattern pattern( a, b );
710
711 if (type == "CDT")
712 {
713 displayCDT(pattern, d);
714 }
715 else if (type == "FDT")
716 {
717 displayConvexHull(pattern, d, true);
718 }
719 else if (type == "CH")
720 {
721 displayConvexHull(pattern, d);
722 }
723 else
724 {
725 trace.error() << " unknown output type. Try -h option. ";
726 trace.info() << std::endl;
727 return EXIT_FAILURE;
728 }
729
730 return EXIT_SUCCESS;
731}
DGtal
Definition ATu0v1.h:57