Skip to content
Snippets Groups Projects
Commit 3425064e authored by Vincent Nivoliers's avatar Vincent Nivoliers
Browse files

multiple hull versions

parent 5e62402e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Code/Makefile
+ 10
− 1
View file @ 3425064e
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
Code/hull_almost_exact.cpp 0 → 100755
+ 310
− 0
View file @ 3425064e
#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 ;
}
Code/hull_permutation.cpp 0 → 100755
+ 293
− 0
View file @ 3425064e
#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 ;
}
Code/hull_perturbation.cpp 0 → 100755
+ 277
− 0
View file @ 3425064e
#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 ;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment