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