#include "svg.hpp"

#include <boost/numeric/interval.hpp>
#include <boost/multiprecision/cpp_int.hpp>

#include <iostream>
#include <algorithm>
#include <random>
#include <vector>
#include <sstream>
#include <cmath>

//{{{ struct vec

struct vec {
  vec() {} ;

  vec(float x, float y) : x(x), y(y) {}
  float x ;
  float y ;

  vec operator-(const vec& v) const {
    return vec(x - v.x, y - v.y) ;
  }

  vec operator+(const vec& v) const {
    return vec(x + v.x, y + v.y) ;
  }

  vec operator*(float a) const {
    return vec(a*x, a*y) ;
  }

  float length2() const {
    return x*x + y*y ;
  }
  
  float length() const {
    return sqrt(length2()) ;
  }
} ;

//}}}

int sign(float x) {
  if(x > 0) return 1 ;
  if(x < 0) return -1 ;
  return 0 ;
}

int orient(const vec& v0, const vec& v1) {
  //determinant
  float d = v0.x*v1.y - v0.y*v1.x ;
  return sign(d) ;
}

struct vec_compare {
  //functor given leftmost point for sort
  vec_compare(const vec& origin) : origin(origin) {}

  //comparison
  bool operator()(const vec& p0, const vec& p1) {
    int o = orient(p0-origin, p1-origin) ;
    return o < 0 ;
  }
  const vec& origin ;
} ;

//{{{leftmost

size_t leftmost(const std::vector<vec>& points) {
  size_t min =  0 ;
  //iterato on all the points
  for(size_t i = 1; i < points.size(); ++i) {
    if(points[i].x < points[min].x) {
      //update the leftmost
      min = i ;
    }
  }
  return min ;
}

//}}}

void compute_hull(std::vector<vec>& points, std::vector<vec>& hull) {
  //leftmost point ahead
  size_t l = leftmost(points) ;
  std::swap(*points.begin(), *(points.begin() + l)) ;
  //sort the rest by angle
  std::sort(points.begin() + 1, points.end(), vec_compare(points[0])) ;

  //add the first two points to the hull
  hull.push_back(points[0]) ;
  hull.push_back(points[1]) ;
  //sweep the remaining points
  for(size_t i = 2; i < points.size(); ++i) {
    //check whether the last segment creates a wrong turn
    const vec& p = points[i] ;
    vec v1 = p - *(hull.end() - 2) ;
    vec v2 = *(hull.end() - 1) - *(hull.end() - 2) ;
    while(hull.size() > 1 && orient(v1, v2) < 0) {
      //in case of wrong turn, rease the last point of the hull and iterate
      hull.erase(hull.end() - 1) ;
      v1 = p - *(hull.end() - 2) ;
      v2 = *(hull.end() - 1) - *(hull.end() - 2) ;
    }
    hull.push_back(p) ;
  }
}

//{{{ generate svg

void svg_hull(
    const std::string& filename, 
    const std::vector<vec>& points, 
    const std::vector<vec>& hull,
    float salt = 0
    ) 
{
  //open file for export
  std::ofstream file(filename) ;
  SvgPad svg(file, 1024, 1024) ;
  svg.open() ;

  //plot data points
  for(const vec& p : points) {
    svg.solid_point(p.x, p.y) ;
  }

  //randomness in case salt is added for hull debugging
  std::mt19937 mt ;
  std::uniform_real_distribution<float> rf(-salt, salt) ;

  //plot hull
  size_t hs = hull.size() ;
  for(size_t i = 0; i < hs; ++i) {

    //salt to perturb the line endpoints and avoid overlaps
    vec salti(0,0) ;
    vec saltip(0,0) ;
    if(salt > 0) {
      salti.x = rf(mt) ;
      salti.y = rf(mt) ;
      saltip.x = rf(mt) ;
      saltip.y = rf(mt) ;
    }

    //hull points
    svg.contour_point(hull[i].x, hull[i].y) ;

    //hull sides with salt
    svg.line(
        hull[i].x + salti.x, 
        hull[i].y + salti.y, 
        hull[(i+1)%hs].x + saltip.x,
        hull[(i+1)%hs].y + saltip.y
    ) ;
  }

  //close file
  svg.close() ;
}

//utility to illustrate angle sorting
void svg_sort(
    const std::string& filename, 
    const std::vector<vec>& points) 
{
  //open svg file for export
  std::ofstream file(filename) ;
  SvgPad svg(file, 1024, 1024) ;
  svg.open() ;

  //colors for the linear gradient of lines
  float col0[3] = { 1., 0.9, 0.3 } ;
  float col1[3] = { 0.35, 0.2, 0.4 } ;

  //iterate on every point except the leftmost at 0
  for(size_t i = 1; i < points.size(); ++i) {
    //hue
    float b = ((float) i-1) / (points.size() - 1) ;
    float a = 1-b ;
    svg.stroke(a*col0[0] + b*col1[0], a*col0[1] + b*col1[1], a*col0[2] + b*col1[2]) ;
    //line defining the sorting angle
    svg.line(points[0].x, points[0].y, points[i].x, points[i].y) ;
  }

  svg.close() ;
}

//}}}

int main() {

  unsigned int size = 100 ;

  //{{{ random test

  {
    //random numbers
    std::random_device rd ;
    std::mt19937 alea(rd()) ;
    std::uniform_real_distribution<float> uniform_float(0, 1) ;
    std::normal_distribution<float> normal_float(0, 1) ;

    std::vector<vec> points ;
    for(unsigned int i = 0; i < size; ++i) {
      //generate a uniform point in a ball
      vec p(normal_float(alea), normal_float(alea)) ;
      float n = sqrt(uniform_float(alea)) / 2 ;
      p = p * (n / p.length()) + vec(0.5,0.5);
      points.push_back(p) ;
    }

    //hull
    std::vector<vec> hull ;
    compute_hull(points, hull) ;
    svg_hull("/tmp/hull_test.svg", points, hull) ;
    svg_sort("/tmp/hull_lines.svg", points) ;
  }

  //}}}

  //{{{ wicked test

  {
    std::vector<vec> aligned_points ;
    //aligned points
    aligned_points.emplace_back(0, 1) ;
    aligned_points.emplace_back(0, 0) ;
    aligned_points.emplace_back(0, 0.8) ;
    aligned_points.emplace_back(0, 0.5) ;

    //test every permutation of the input
    int indices[4] = {0, 1, 2, 3} ;
    int permutation_index = 0 ;
    do {
      std::vector<vec> points ;
      for(unsigned int i = 0; i < 4; ++i) {
        points.push_back(aligned_points[indices[i]]) ;
      }
      //make the hull non flat
      points.emplace_back(1, 0.5) ;

      //hull
      std::vector<vec> hull ;
      compute_hull(points, hull) ;
      
      //export with salt for readability
      std::stringstream ss ;
      ss << "/tmp/hull_aligned_" << permutation_index << ".svg" ;
      svg_hull(ss.str(), points, hull, 0.02) ;

      //move to next permutation
      ++ permutation_index ;
    }while(std::next_permutation(indices, indices + 4)) ;
  }

  //}}}

  return 0 ;
}