From 3425064ebe231d03c6e1fa1ac01b7fa5d00a9674 Mon Sep 17 00:00:00 2001
From: Vincent Nivoliers <vincent.nivoliers@univ-lyon1.fr>
Date: Mon, 18 Nov 2019 19:08:06 +0100
Subject: [PATCH] multiple hull versions
---
Code/Makefile | 11 +-
Code/hull_almost_exact.cpp | 310 +++++++++++++++++++++++++++++++++++++
Code/hull_permutation.cpp | 293 +++++++++++++++++++++++++++++++++++
Code/hull_perturbation.cpp | 277 +++++++++++++++++++++++++++++++++
4 files changed, 890 insertions(+), 1 deletion(-)
create mode 100755 Code/hull_almost_exact.cpp
create mode 100755 Code/hull_permutation.cpp
create mode 100755 Code/hull_perturbation.cpp
diff --git a/Code/Makefile b/Code/Makefile
index cc59505..1d4ab52 100644
--- a/Code/Makefile
+++ b/Code/Makefile
@@ -1,4 +1,4 @@
-all: precision hull
+all: precision hull hull_perturbation hull_permutation hull_almost_exact
precision: precision.cpp
g++ -g -Wall precision.cpp -o precision
@@ -6,5 +6,14 @@ precision: precision.cpp
hull: hull.cpp svg.cpp
g++ -g -Wall hull.cpp svg.cpp -o hull
+hull_perturbation: hull_perturbation.cpp svg.cpp
+ g++ -g -Wall hull_perturbation.cpp svg.cpp -o hull_perturbation
+
+hull_permutation: hull_permutation.cpp svg.cpp
+ g++ -g -Wall hull_permutation.cpp svg.cpp -o hull_permutation
+
+hull_almost_exact: hull_almost_exact.cpp svg.cpp
+ g++ -g -Wall hull_almost_exact.cpp svg.cpp -o hull_almost_exact
+
clean:
rm -f precision hull
diff --git a/Code/hull_almost_exact.cpp b/Code/hull_almost_exact.cpp
new file mode 100755
index 0000000..a561335
--- /dev/null
+++ b/Code/hull_almost_exact.cpp
@@ -0,0 +1,310 @@
+#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 ;
+}
diff --git a/Code/hull_permutation.cpp b/Code/hull_permutation.cpp
new file mode 100755
index 0000000..83b0fe5
--- /dev/null
+++ b/Code/hull_permutation.cpp
@@ -0,0 +1,293 @@
+#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()) ;
+ }
+
+ bool lexicoinf(const vec& v) const {
+ if(x == v.x) return y < v.y ;
+ return x < v.x ;
+ }
+} ;
+
+//}}}
+
+int sign(float 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
+ 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
+ float 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) ;
+ float 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 ;
+}
diff --git a/Code/hull_perturbation.cpp b/Code/hull_perturbation.cpp
new file mode 100755
index 0000000..54ffc51
--- /dev/null
+++ b/Code/hull_perturbation.cpp
@@ -0,0 +1,277 @@
+#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 ;
+ if(d != 0) return sign(d) ;
+
+ //enter perturbation with perturbation matrix
+ // [ e e^2 ]
+ // [ e^2 e^2 ]
+ if(v1.y != 0) return sign(v1.y) ;
+ float a = v0.x - v0.y - v1.x ;
+ if(a != 0) return sign(a) ;
+ return 1 ;
+
+ //other perturbations work, for instance
+ // [ e 0 ]
+ // [ 0 e^2 ]
+ if(v1.y != 0) return sign(v1.y) ;
+ if(v0.x != 0) return sign(v0.x) ;
+ return 1 ;
+}
+
+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 ;
+}
--
GitLab