patternTriangulation.cpp

#include <iostream>
#include <iterator>
#include <cstdio>
#include <cmath>
#include <fstream>
#include <vector>

#include <boost/program_options/options_description.hpp>
#include <boost/program_options/parsers.hpp>
#include <boost/program_options/variables_map.hpp>

#include "DGtal/base/Common.h"
#include "DGtal/helpers/StdDefs.h"

#include "DGtal/arithmetic/LighterSternBrocot.h"
#include "DGtal/arithmetic/Pattern.h"

#include "DGtal/io/boards/Board2D.h"
#include "DGtal/io/Color.h"

using namespace DGtal;
using namespace std;

namespace po = boost::program_options;

template <typename Point>
struct OddConvexHullMap {

public: 
  typedef Point Vector; 

public: 
  OddConvexHullMap(const Vector& aShift = Vector(1,-1))
    : myS(aShift) {}
  
  Point operator()(const Point& aP) const
  {
    return aP + myS; 
  }

private: 
  Vector myS; 

}; 

template <typename Point>
struct EvenConvexHullMap {

public: 
  typedef Point Vector; 

public: 
  EvenConvexHullMap(const Point& aPoint, const Vector& aShift = Vector(1,-1))
    : myZn(aPoint), myS(aShift) {}
  
  Point operator()(const Point& aP) const
  {
    return myZn - aP + myS; 
  }

private: 
  Point myZn; 
  Vector myS; 

}; 


template <typename Board, typename Point>
void drawSegment(Board& aBoard, const Point& aP, const Point& aQ)
{
  aBoard.drawLine(aP[0], aP[1], aQ[0], aQ[1]);   
  aBoard << aP << aQ; 
}

template <typename Board, typename Iterator, typename Functor>
void drawSegments(Board& aBoard, 
                  const Iterator& itb, const Iterator& ite, 
                  const Functor& aF)
{
  Iterator it = itb; 
  if (it != ite) 
    {
      Iterator prevIt = it; 
      ++it; 
      for ( ; it != ite; ++prevIt, ++it) 
        {
          drawSegment( aBoard, aF(*prevIt), aF(*it) ); 
        }
    }
}



template<typename Fraction, typename Container>
void fillConvergents(const Fraction& aZn, 
                Container& odd, Container& even)
{
  typedef typename Container::value_type Vector; 

  for (typename Fraction::Quotient i = 1; i <= aZn.k(); ++i)
    {
      Fraction zk = aZn.reduced(i); 
      if (((aZn.k() - i)%2) == 1 )
        { //odd
          odd.push_back(Vector(zk.q(),zk.p())); 
        }
      else 
        { //even
          even.push_back(Vector(zk.q(),zk.p())); 
        }
    }
}

template <typename Pattern, typename Integer>
void displayConvexHull(const Pattern& aP, const Integer& aD, bool withFDT = false)
{
  typedef typename Pattern::Fraction Fraction; 
  typedef typename Pattern::Quotient Quotient; 
  typedef typename Pattern::Vector2I Vector;  
  typedef std::vector<Vector> Convergents; 
  typedef typename Pattern::Point2I Point;  
  
  Board2D aBoard;

  //list of convergents
  Convergents oddConvergents, evenConvergents;

  //fill the lists
  Fraction zn = aP.slope();
  // aD > 1
  if (aD > 1)
    {
      Fraction znm1 = zn.father(); 
      if (zn.odd())
        oddConvergents.push_back(Vector(znm1.q(),znm1.p())); 
      else 
        evenConvergents.push_back(Vector(znm1.q(),znm1.p())); 
    }
  // aD >= 1
  fillConvergents( zn, oddConvergents, evenConvergents ); 

  //display the domain
  Point Zn = Point(aD*zn.q(),aD*zn.p()); 
  aBoard << Z2i::Domain( Point(0,0), Zn ); 
  aBoard << SetMode( Zn.className(), "Grid" );
  aBoard.setLineStyle(Board2D::Shape::SolidStyle); 
  aBoard.setPenColor(DGtal::Color::Black);

  //display the two main lists
  OddConvexHullMap<Point> oddH;
  drawSegments( aBoard, oddConvergents.begin(), oddConvergents.end(), oddH ); 

  EvenConvexHullMap<Point> evenH(Zn);
  drawSegments( aBoard, evenConvergents.begin(), evenConvergents.end(), evenH ); 

  //display the four last segments
  drawSegment( aBoard, Point(0,0), Zn ); 
  if (evenConvergents.size() != 0)
    {
      drawSegment( aBoard, evenH(*evenConvergents.rbegin()), Zn ); 
      if (oddConvergents.size() != 0)
        {
          drawSegment( aBoard, Point(0,0), oddH(*oddConvergents.rbegin()) );
          drawSegment( aBoard, 
                       oddH(*oddConvergents.begin()), 
                       evenH(*evenConvergents.begin()) ); 
        }
      else 
        {
          drawSegment( aBoard, Point(0,0), evenH(*evenConvergents.rbegin()) );
        }
    }

  if (withFDT)
    {//display internal edges

      //first internal edge
      if (evenConvergents.size() != 0)
        drawSegment( aBoard, 
                     Point(0,0), 
                     evenH(*evenConvergents.rbegin()) ); 
      //other internal edges
      typedef typename Convergents::const_reverse_iterator RI; 
      RI oddRit = oddConvergents.rbegin(); 
      RI evenRit = evenConvergents.rbegin(); 
      bool hasToContinue = true; 
      while (hasToContinue) 
        {
          if (oddRit != oddConvergents.rend())
            {         
              drawSegment( aBoard, 
                           oddH(*oddRit), 
                           evenH(*evenRit) ); 
            }
          else 
            hasToContinue = false; 
          ++evenRit; 

          if ( (hasToContinue) && (evenRit != evenConvergents.rend()) )
            {         
              drawSegment( aBoard, 
                           oddH(*oddRit), 
                           evenH(*evenRit) ); 
            }
          else 
            hasToContinue = false; 
          ++oddRit; 
        }

      aBoard.saveEPS("FDT.eps");
    }
  else
    {
      aBoard.saveEPS("CH.eps");
    }  
}


template <typename Board, typename Point>
void drawTriangle(Board& aBoard, 
                  const Point& aP, const Point& aQ, const Point& aR)
{
  aBoard.drawTriangle( aP[0], aP[1], aQ[0], aQ[1], aR[0], aR[1] ); 
}

struct InDirectBezoutComputer
{
  template<typename Pattern>
  typename Pattern::Vector2I
  getPositiveBezout(const Pattern& aP) const
  {
    return aP.bezout(); 
  }

  template<typename Pattern>
  typename Pattern::Vector2I
  getNegativeBezout(const Pattern& aP) const
  {
    return aP.slope().even() 
      ? aP.U( 1 ) - aP.previousPattern().U( 1 ) 
      : aP.previousPattern().U( 1 );
  }

}; 
 
struct DirectBezoutComputer
{
  template<typename Pattern>
  typename Pattern::Vector2I
  getPositiveBezout(const Pattern& aP) const
  {
    return aP.slope().even() 
      ? aP.U( 1 ) - aP.previousPattern().U( 1 ) 
      : aP.previousPattern().U( 1 );
  }

  template<typename Pattern>
  typename Pattern::Vector2I
  getNegativeBezout(const Pattern& aP) const
  {
    return aP.bezout(); 
  }

}; 

template <typename Board, typename Pattern, typename BezoutComputer>
void displayPartialCDT(Board& aBoard, 
                       const Pattern& aP, 
                       const typename Pattern::Point2I& aStartingPoint,
                       const BezoutComputer& aBC )
{
  typedef typename Pattern::Fraction Fraction; 
  typedef typename Pattern::Vector2I Vector;  
  typedef typename Pattern::Point2I Point;  
  
  Fraction f = aP.slope();
  if ( (f.p() >= 1)&&(f.q() >= 1) )
    {
      Point O = aStartingPoint; 
      Point Zn = aStartingPoint + aP.v();
      Vector v1 = aBC.getNegativeBezout(aP); 
      Point B = aStartingPoint + v1; 
      drawTriangle(aBoard, O, Zn, B); 

      if ( (f.p() > 1)||(f.q() > 1) )
        {
          //recursive calls 
          // - first pattern
          displayPartialCDT(aBoard, Pattern(Fraction(v1[1],v1[0])), O, aBC); 
          // - second pattern 
          Vector v2 = aBC.getPositiveBezout(aP); 
          displayPartialCDT(aBoard, Pattern(Fraction(v2[1],v2[0])), B, aBC); 
        }
    }
}

template <typename Pattern, typename Integer>
void displayCDT(const Pattern& aP, const Integer& aD)
{
  Board2D aBoard; 

  typedef typename Pattern::Point2I Point;  
  typedef typename Pattern::Quotient Quotient;  
  typedef typename Pattern::Fraction Fraction;  
  Fraction zn = aP.slope(); 

  //display the domain
  Point Zn = Point(aD*zn.q(),aD*zn.p()); 
  aBoard << Z2i::Domain( Point(0,0), Zn ); 
  aBoard << SetMode( Zn.className(), "Grid" );
  aBoard.setLineStyle(Board2D::Shape::SolidStyle); 
  aBoard.setPenColor(DGtal::Color::Black);

  //display the upper part of the triangulation
  for (Integer i = 0; i < aD; ++i)
    {
      displayPartialCDT( aBoard, aP, Point(i*zn.q(), i*zn.p()), InDirectBezoutComputer() );
    }

  //display the lower parts of the triangulation
  // - getting the reversed patterns
  typedef typename Pattern::Vector2I Vector;  
  typedef std::vector<Vector> Convergents; 
  typedef typename std::vector<Vector>::const_reverse_iterator ReverseIterator; 
  Convergents oddConvergents, evenConvergents;
  fillConvergents( zn, oddConvergents, evenConvergents ); 
  // - displaying the reversed patterns...
  Point rPStartingPoint; 
  //  - of odd slope: 
  OddConvexHullMap<Point> oddH;
  {
    ReverseIterator itb = evenConvergents.rbegin(); 
    ReverseIterator ite = evenConvergents.rend(); 
    Point aStartingPoint = oddH( *oddConvergents.rbegin() ); 
    //
    Point p = aStartingPoint; 
    ReverseIterator it = itb; 
    if (it != ite)
      {
        for (++it; it != ite; ++it) 
          {
            displayPartialCDT( aBoard, Pattern( (*it)[1], (*it)[0] ), p, DirectBezoutComputer() ); 
            p += *it; 
          }
      }
    rPStartingPoint = p; 
  }
  //  - of even slope: 
  EvenConvexHullMap<Point> evenH(Zn);
  {
    ReverseIterator itb = oddConvergents.rbegin(); 
    ReverseIterator ite = oddConvergents.rend(); 
    Point aStartingPoint = evenH( *evenConvergents.rbegin() ); 
    //
    Point p = aStartingPoint;
    ReverseIterator it = itb; 
    for ( ; it != ite; ++it) 
      {
        p -= *it; 
        displayPartialCDT( aBoard, Pattern( (*it)[1], (*it)[0] ), p, DirectBezoutComputer() ); 
      }
  }
  //  - of same slope: 
  {
    for (Integer i = 0; i < (aD-1); ++i)
      {
        Point p = rPStartingPoint + i*aP.v(); 
        displayPartialCDT( aBoard, aP, p, DirectBezoutComputer() );
      }
  }

 
  aBoard.saveEPS("CDT.eps");
}



void missingParam(std::string param)
{
  trace.error() << " Parameter: " << param << " is required...";
  trace.info() << std::endl;
  exit(1);
}



int main(int argc, char **argv)
{
  // parse command line ----------------------------------------------
  po::options_description general_opt("Allowed options are");
  general_opt.add_options()
    ("help,h", "display this message")
    ("aparam,a",   po::value<int>(), "pattern a parameter" )
    ("bparam,b",   po::value<int>(), "pattern b parameter" )
    ("delta,d",   po::value<int>()->default_value(1), "number of repetitions" )
    ("triangulation,t",   po::value<string>()->default_value("CH"), "output:\n\tClosest-point Delaunay triangulation {CDT}\n\tFarthest-point Delaunay triangulation {FDT}\n\tConvex hull {CH}" );

  bool parseOK=true;
  po::variables_map vm;
  try{
    po::store(po::parse_command_line(argc, argv, general_opt), vm);  
  }catch(const std::exception& ex){
    parseOK=false;
    trace.info()<< "Error checking program options: "<< ex.what()<< endl;
  }
 
  po::notify(vm);    
  if(!parseOK || vm.count("help")||argc<=1)
    {
      trace.info()<< "Draw the Delaunay triangulation of a pattern using DGtal library" <<std::endl 
                  << "Basic usage: "<<std::endl
                  << "\t" << argv[0] << " -a 5 -b 8 "<<std::endl
                  << general_opt << "\n";
      return 0;
    }


  if (!(vm.count("aparam"))) missingParam("-a");
  if (!(vm.count("bparam"))) missingParam("-b");

  typedef DGtal::int32_t Integer;
  Integer a = vm["aparam"].as<int>();
  Integer b = vm["bparam"].as<int>();
  if ( (a < 0) || (b <= 0) )
    {
      trace.error() << " parameters a and b should be strictly positive.";
      trace.info() << std::endl;
      return 0;
    }
  if (a >= b)
    {
      trace.error() << " parameter a should be strictly less than b.";
      trace.info() << std::endl;
      return 0;
    }
  Integer d = vm["delta"].as<int>(); 
  if (d <= 0)
    {
      trace.error() << " parameter d should be strictly positive";
      trace.info() << std::endl;
      return 0;
    }
  trace.info() << "a=" << a << ", b=" << b << ", d=" << d << std::endl; 

  typedef DGtal::int32_t Quotient;
  typedef LighterSternBrocot<Integer, Quotient, StdMapRebinder> SB;
  typedef SB::Fraction Fraction; // the type for fractions
  typedef Pattern<Fraction> Pattern; // the type for patterns
 
  Pattern pattern( a, b );

  string type = vm["triangulation"].as<string>();
  if (type == "CDT")
    {
      displayCDT(pattern, d); 
    }
  else if (type == "FDT")
    {
      displayConvexHull(pattern, d, true); 
    }
  else if (type == "CH")
    {
      displayConvexHull(pattern, d); 
    }
  else 
    {
      trace.error() << " unknown output type. Try -h option. ";
      trace.info() << std::endl;
      return 0;
    }

  return 1;
}