hull_almost_exact.cpp

#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

//generic vector
template<typename number>
struct generic_vec {
  generic_vec() {} ;

  generic_vec(number x, number y) : x(x), y(y) {}
  number x ;
  number y ;

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

  template<typename other_number>
  generic_vec& operator=(const generic_vec<other_number>& v) {
    x = v.x ;
    y = v.y ;
    return *this ;
  }

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

  generic_vec operator*(number a) const {
    return generic_vec(a*x, a*y) ;
  }

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

  bool lexicoinf(const generic_vec& v) const {
    if(x == v.x) return y < v.y ;
    return x < v.x ;
  }
} ;

//redefine the previous vector with floats
using vec = generic_vec<float> ;

//exact vector on rationals
using rational = boost::multiprecision::cpp_rational ;
using exact_vec = generic_vec<rational> ;

//}}}

template<typename number>
int sign(number x) {
  if(x > 0) return 1 ;
  if(x < 0) return -1 ;
  return 0 ;
}

int orient(const vec& v0_in, const vec& v1_in) {
  //sort the inputs to resist permutation and use exactness
  exact_vec v0, v1 ;
  int signature = 1 ;
  if(v0_in.lexicoinf(v1_in)) {
    v0 = v0_in ;
    v1 = v1_in ;
  } else {
    v1 = v0_in ;
    v0 = v1_in ;
    signature = -1 ;
  }
  //determinant
  rational d = v0.x*v1.y - v0.y*v1.x ;
  if(d != 0) return signature*sign(d) ;

  //enter perturbation with perturbation matrix
  //  [ e   e^2 ]
  //  [ e^2 e^2 ]
  if(v1.y != 0) return signature*sign(v1.y) ;
  rational a = v0.x - v0.y - v1.x ;
  if(a != 0) return signature*sign(a) ;
  return signature ;

  //other perturbations work, for instance
  //  [ e  0   ]
  //  [ 0  e^2 ]
  if(v1.y != 0) return signature*sign(v1.y) ;
  if(v0.x != 0) return signature*sign(v0.x) ;
  return signature ;
}

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 ;
}