diff --git a/tests/d2g b/tests/d2g deleted file mode 100755 index 9cad7b7..0000000 Binary files a/tests/d2g and /dev/null differ diff --git a/tests/delaunator.hpp b/tests/delaunator.hpp deleted file mode 100644 index 6bc5221..0000000 --- a/tests/delaunator.hpp +++ /dev/null @@ -1,585 +0,0 @@ -#pragma once - -#include -#include -#include -#include -#include -#include -#include -#include - -namespace delaunator { - -//@see https://stackoverflow.com/questions/33333363/built-in-mod-vs-custom-mod-function-improve-the-performance-of-modulus-op/33333636#33333636 -inline size_t fast_mod(const size_t i, const size_t c) { - return i >= c ? i % c : i; -} - -// Kahan and Babuska summation, Neumaier variant; accumulates less FP error -inline double sum(const std::vector& x) { - double sum = x[0]; - double err = 0.0; - - for (size_t i = 1; i < x.size(); i++) { - const double k = x[i]; - const double m = sum + k; - err += std::fabs(sum) >= std::fabs(k) ? sum - m + k : k - m + sum; - sum = m; - } - return sum + err; -} - -inline double dist( - const double ax, - const double ay, - const double bx, - const double by) { - const double dx = ax - bx; - const double dy = ay - by; - return dx * dx + dy * dy; -} - -inline double circumradius( - const double ax, - const double ay, - const double bx, - const double by, - const double cx, - const double cy) { - const double dx = bx - ax; - const double dy = by - ay; - const double ex = cx - ax; - const double ey = cy - ay; - - const double bl = dx * dx + dy * dy; - const double cl = ex * ex + ey * ey; - const double d = dx * ey - dy * ex; - - const double x = (ey * bl - dy * cl) * 0.5 / d; - const double y = (dx * cl - ex * bl) * 0.5 / d; - - if ((bl > 0.0 || bl < 0.0) && (cl > 0.0 || cl < 0.0) && (d > 0.0 || d < 0.0)) { - return x * x + y * y; - } else { - return std::numeric_limits::max(); - } -} - -inline bool orient( - const double px, - const double py, - const double qx, - const double qy, - const double rx, - const double ry) { - return (qy - py) * (rx - qx) - (qx - px) * (ry - qy) < 0.0; -} - -inline std::pair circumcenter( - const double ax, - const double ay, - const double bx, - const double by, - const double cx, - const double cy) { - const double dx = bx - ax; - const double dy = by - ay; - const double ex = cx - ax; - const double ey = cy - ay; - - const double bl = dx * dx + dy * dy; - const double cl = ex * ex + ey * ey; - const double d = dx * ey - dy * ex; - - const double x = ax + (ey * bl - dy * cl) * 0.5 / d; - const double y = ay + (dx * cl - ex * bl) * 0.5 / d; - - return std::make_pair(x, y); -} - -struct compare { - - std::vector const& coords; - double cx; - double cy; - - bool operator()(std::size_t i, std::size_t j) { - const double d1 = dist(coords[2 * i], coords[2 * i + 1], cx, cy); - const double d2 = dist(coords[2 * j], coords[2 * j + 1], cx, cy); - const double diff1 = d1 - d2; - const double diff2 = coords[2 * i] - coords[2 * j]; - const double diff3 = coords[2 * i + 1] - coords[2 * j + 1]; - - if (diff1 > 0.0 || diff1 < 0.0) { - return diff1 < 0; - } else if (diff2 > 0.0 || diff2 < 0.0) { - return diff2 < 0; - } else { - return diff3 < 0; - } - } -}; - -inline bool in_circle( - const double ax, - const double ay, - const double bx, - const double by, - const double cx, - const double cy, - const double px, - const double py) { - const double dx = ax - px; - const double dy = ay - py; - const double ex = bx - px; - const double ey = by - py; - const double fx = cx - px; - const double fy = cy - py; - - const double ap = dx * dx + dy * dy; - const double bp = ex * ex + ey * ey; - const double cp = fx * fx + fy * fy; - - return (dx * (ey * cp - bp * fy) - - dy * (ex * cp - bp * fx) + - ap * (ex * fy - ey * fx)) < 0.0; -} - -constexpr double EPSILON = std::numeric_limits::epsilon(); -constexpr std::size_t INVALID_INDEX = std::numeric_limits::max(); - -inline bool check_pts_equal(double x1, double y1, double x2, double y2) { - return std::fabs(x1 - x2) <= EPSILON && - std::fabs(y1 - y2) <= EPSILON; -} - -// monotonically increases with real angle, but doesn't need expensive trigonometry -inline double pseudo_angle(const double dx, const double dy) { - const double p = dx / (std::abs(dx) + std::abs(dy)); - return (dy > 0.0 ? 3.0 - p : 1.0 + p) / 4.0; // [0..1) -} - -struct DelaunatorPoint { - std::size_t i; - double x; - double y; - std::size_t t; - std::size_t prev; - std::size_t next; - bool removed; -}; - -class Delaunator { - -public: - std::vector const& coords; - std::vector triangles; - std::vector halfedges; - std::vector hull_prev; - std::vector hull_next; - std::vector hull_tri; - std::size_t hull_start; - - Delaunator(std::vector const& in_coords); - - double get_hull_area(); - -private: - std::vector m_hash; - double m_center_x; - double m_center_y; - std::size_t m_hash_size; - std::vector m_edge_stack; - - std::size_t legalize(std::size_t a); - std::size_t hash_key(double x, double y) const; - std::size_t add_triangle( - std::size_t i0, - std::size_t i1, - std::size_t i2, - std::size_t a, - std::size_t b, - std::size_t c); - void link(std::size_t a, std::size_t b); -}; - -Delaunator::Delaunator(std::vector const& in_coords) - : coords(in_coords), - triangles(), - halfedges(), - hull_prev(), - hull_next(), - hull_tri(), - hull_start(), - m_hash(), - m_center_x(), - m_center_y(), - m_hash_size(), - m_edge_stack() { - std::size_t n = coords.size() >> 1; - - double max_x = std::numeric_limits::min(); - double max_y = std::numeric_limits::min(); - double min_x = std::numeric_limits::max(); - double min_y = std::numeric_limits::max(); - std::vector ids; - ids.reserve(n); - - for (std::size_t i = 0; i < n; i++) { - const double x = coords[2 * i]; - const double y = coords[2 * i + 1]; - - if (x < min_x) min_x = x; - if (y < min_y) min_y = y; - if (x > max_x) max_x = x; - if (y > max_y) max_y = y; - - ids.push_back(i); - } - const double cx = (min_x + max_x) / 2; - const double cy = (min_y + max_y) / 2; - double min_dist = std::numeric_limits::max(); - - std::size_t i0 = INVALID_INDEX; - std::size_t i1 = INVALID_INDEX; - std::size_t i2 = INVALID_INDEX; - - // pick a seed point close to the centroid - for (std::size_t i = 0; i < n; i++) { - const double d = dist(cx, cy, coords[2 * i], coords[2 * i + 1]); - if (d < min_dist) { - i0 = i; - min_dist = d; - } - } - - const double i0x = coords[2 * i0]; - const double i0y = coords[2 * i0 + 1]; - - min_dist = std::numeric_limits::max(); - - // find the point closest to the seed - for (std::size_t i = 0; i < n; i++) { - if (i == i0) continue; - const double d = dist(i0x, i0y, coords[2 * i], coords[2 * i + 1]); - if (d < min_dist && d > 0.0) { - i1 = i; - min_dist = d; - } - } - - double i1x = coords[2 * i1]; - double i1y = coords[2 * i1 + 1]; - - double min_radius = std::numeric_limits::max(); - - // find the third point which forms the smallest circumcircle with the first two - for (std::size_t i = 0; i < n; i++) { - if (i == i0 || i == i1) continue; - - const double r = circumradius( - i0x, i0y, i1x, i1y, coords[2 * i], coords[2 * i + 1]); - - if (r < min_radius) { - i2 = i; - min_radius = r; - } - } - - if (!(min_radius < std::numeric_limits::max())) { - throw std::runtime_error("not triangulation"); - } - - double i2x = coords[2 * i2]; - double i2y = coords[2 * i2 + 1]; - - if (orient(i0x, i0y, i1x, i1y, i2x, i2y)) { - std::swap(i1, i2); - std::swap(i1x, i2x); - std::swap(i1y, i2y); - } - - std::tie(m_center_x, m_center_y) = circumcenter(i0x, i0y, i1x, i1y, i2x, i2y); - - // sort the points by distance from the seed triangle circumcenter - std::sort(ids.begin(), ids.end(), compare{ coords, m_center_x, m_center_y }); - - // initialize a hash table for storing edges of the advancing convex hull - m_hash_size = static_cast(std::llround(std::ceil(std::sqrt(n)))); - m_hash.resize(m_hash_size); - std::fill(m_hash.begin(), m_hash.end(), INVALID_INDEX); - - // initialize arrays for tracking the edges of the advancing convex hull - hull_prev.resize(n); - hull_next.resize(n); - hull_tri.resize(n); - - hull_start = i0; - - size_t hull_size = 3; - - hull_next[i0] = hull_prev[i2] = i1; - hull_next[i1] = hull_prev[i0] = i2; - hull_next[i2] = hull_prev[i1] = i0; - - hull_tri[i0] = 0; - hull_tri[i1] = 1; - hull_tri[i2] = 2; - - m_hash[hash_key(i0x, i0y)] = i0; - m_hash[hash_key(i1x, i1y)] = i1; - m_hash[hash_key(i2x, i2y)] = i2; - - std::size_t max_triangles = n < 3 ? 1 : 2 * n - 5; - triangles.reserve(max_triangles * 3); - halfedges.reserve(max_triangles * 3); - add_triangle(i0, i1, i2, INVALID_INDEX, INVALID_INDEX, INVALID_INDEX); - double xp = std::numeric_limits::quiet_NaN(); - double yp = std::numeric_limits::quiet_NaN(); - for (std::size_t k = 0; k < n; k++) { - const std::size_t i = ids[k]; - const double x = coords[2 * i]; - const double y = coords[2 * i + 1]; - - // skip near-duplicate points - if (k > 0 && check_pts_equal(x, y, xp, yp)) continue; - xp = x; - yp = y; - - // skip seed triangle points - if ( - check_pts_equal(x, y, i0x, i0y) || - check_pts_equal(x, y, i1x, i1y) || - check_pts_equal(x, y, i2x, i2y)) continue; - - // find a visible edge on the convex hull using edge hash - std::size_t start = 0; - - size_t key = hash_key(x, y); - for (size_t j = 0; j < m_hash_size; j++) { - start = m_hash[fast_mod(key + j, m_hash_size)]; - if (start != INVALID_INDEX && start != hull_next[start]) break; - } - - start = hull_prev[start]; - size_t e = start; - size_t q; - - while (q = hull_next[e], !orient(x, y, coords[2 * e], coords[2 * e + 1], coords[2 * q], coords[2 * q + 1])) { //TODO: does it works in a same way as in JS - e = q; - if (e == start) { - e = INVALID_INDEX; - break; - } - } - - if (e == INVALID_INDEX) continue; // likely a near-duplicate point; skip it - - // add the first triangle from the point - std::size_t t = add_triangle( - e, - i, - hull_next[e], - INVALID_INDEX, - INVALID_INDEX, - hull_tri[e]); - - hull_tri[i] = legalize(t + 2); - hull_tri[e] = t; - hull_size++; - - // walk forward through the hull, adding more triangles and flipping recursively - std::size_t next = hull_next[e]; - while ( - q = hull_next[next], - orient(x, y, coords[2 * next], coords[2 * next + 1], coords[2 * q], coords[2 * q + 1])) { - t = add_triangle(next, i, q, hull_tri[i], INVALID_INDEX, hull_tri[next]); - hull_tri[i] = legalize(t + 2); - hull_next[next] = next; // mark as removed - hull_size--; - next = q; - } - - // walk backward from the other side, adding more triangles and flipping - if (e == start) { - while ( - q = hull_prev[e], - orient(x, y, coords[2 * q], coords[2 * q + 1], coords[2 * e], coords[2 * e + 1])) { - t = add_triangle(q, i, e, INVALID_INDEX, hull_tri[e], hull_tri[q]); - legalize(t + 2); - hull_tri[q] = t; - hull_next[e] = e; // mark as removed - hull_size--; - e = q; - } - } - - // update the hull indices - hull_prev[i] = e; - hull_start = e; - hull_prev[next] = i; - hull_next[e] = i; - hull_next[i] = next; - - m_hash[hash_key(x, y)] = i; - m_hash[hash_key(coords[2 * e], coords[2 * e + 1])] = e; - } -} - -double Delaunator::get_hull_area() { - std::vector hull_area; - size_t e = hull_start; - do { - hull_area.push_back((coords[2 * e] - coords[2 * hull_prev[e]]) * (coords[2 * e + 1] + coords[2 * hull_prev[e] + 1])); - e = hull_next[e]; - } while (e != hull_start); - return sum(hull_area); -} - -std::size_t Delaunator::legalize(std::size_t a) { - std::size_t i = 0; - std::size_t ar = 0; - m_edge_stack.clear(); - - // recursion eliminated with a fixed-size stack - while (true) { - const size_t b = halfedges[a]; - - /* if the pair of triangles doesn't satisfy the Delaunay condition - * (p1 is inside the circumcircle of [p0, pl, pr]), flip them, - * then do the same check/flip recursively for the new pair of triangles - * - * pl pl - * /||\ / \ - * al/ || \bl al/ \a - * / || \ / \ - * / a||b \ flip /___ar___\ - * p0\ || /p1 => p0\---bl---/p1 - * \ || / \ / - * ar\ || /br b\ /br - * \||/ \ / - * pr pr - */ - const size_t a0 = 3 * (a / 3); - ar = a0 + (a + 2) % 3; - - if (b == INVALID_INDEX) { - if (i > 0) { - i--; - a = m_edge_stack[i]; - continue; - } else { - //i = INVALID_INDEX; - break; - } - } - - const size_t b0 = 3 * (b / 3); - const size_t al = a0 + (a + 1) % 3; - const size_t bl = b0 + (b + 2) % 3; - - const std::size_t p0 = triangles[ar]; - const std::size_t pr = triangles[a]; - const std::size_t pl = triangles[al]; - const std::size_t p1 = triangles[bl]; - - const bool illegal = in_circle( - coords[2 * p0], - coords[2 * p0 + 1], - coords[2 * pr], - coords[2 * pr + 1], - coords[2 * pl], - coords[2 * pl + 1], - coords[2 * p1], - coords[2 * p1 + 1]); - - if (illegal) { - triangles[a] = p1; - triangles[b] = p0; - - auto hbl = halfedges[bl]; - - // edge swapped on the other side of the hull (rare); fix the halfedge reference - if (hbl == INVALID_INDEX) { - std::size_t e = hull_start; - do { - if (hull_tri[e] == bl) { - hull_tri[e] = a; - break; - } - e = hull_next[e]; - } while (e != hull_start); - } - link(a, hbl); - link(b, halfedges[ar]); - link(ar, bl); - std::size_t br = b0 + (b + 1) % 3; - - if (i < m_edge_stack.size()) { - m_edge_stack[i] = br; - } else { - m_edge_stack.push_back(br); - } - i++; - - } else { - if (i > 0) { - i--; - a = m_edge_stack[i]; - continue; - } else { - break; - } - } - } - return ar; -} - -inline std::size_t Delaunator::hash_key(const double x, const double y) const { - const double dx = x - m_center_x; - const double dy = y - m_center_y; - return fast_mod( - static_cast(std::llround(std::floor(pseudo_angle(dx, dy) * static_cast(m_hash_size)))), - m_hash_size); -} - -std::size_t Delaunator::add_triangle( - std::size_t i0, - std::size_t i1, - std::size_t i2, - std::size_t a, - std::size_t b, - std::size_t c) { - std::size_t t = triangles.size(); - triangles.push_back(i0); - triangles.push_back(i1); - triangles.push_back(i2); - link(t, a); - link(t + 1, b); - link(t + 2, c); - return t; -} - -void Delaunator::link(const std::size_t a, const std::size_t b) { - std::size_t s = halfedges.size(); - if (a == s) { - halfedges.push_back(b); - } else if (a < s) { - halfedges[a] = b; - } else { - throw std::runtime_error("Cannot link edge"); - } - if (b != INVALID_INDEX) { - std::size_t s2 = halfedges.size(); - if (b == s2) { - halfedges.push_back(a); - } else if (b < s2) { - halfedges[b] = a; - } else { - throw std::runtime_error("Cannot link edge"); - } - } -} - -} //namespace delaunator diff --git a/tests/delaunay2graph.cpp b/tests/delaunay2graph.cpp deleted file mode 100644 index 2cfb4e2..0000000 --- a/tests/delaunay2graph.cpp +++ /dev/null @@ -1,143 +0,0 @@ -/***************************************************************** - * IST: ASA 2022/23 - P2 - Pedro T. Monteiro | graph generator * - ***************************************************************** - * Random points between (0,0)-(M,M) -> Delaunay -> Graph * - * Using: https://github.com/delfrrr/delaunator-cpp * - *****************************************************************/ -#include "delaunator.hpp" -#include -#include -#include - -using namespace std; - -unsigned int _V, _M, _currV; -int randLimit; -map edgePool; - -//------------------------------------------------------------------- -void printUsage(char *progname) { - cout << "Usage: " << progname << " #V M p seed" << endl; - cout << "\t#V: total number of vertices" << endl; - cout << "\tM: max (M,M) coordinate where M >= V" << endl; - cout << "\tp: probability to create new graph" << endl; - cout << "\tseed: random seed number (optional)" << endl; - exit(0); -} - -void parseArgs(int argc, char **argv) { - int seed = 0; - - if (argc < 4 || argc > 5) { - cout << "ERROR: invalid # arguments" << endl; - printUsage(argv[0]); - } - - sscanf(argv[1], "%d", &_V); - if (_V < 3) { - cout << "ERROR: #V vertices must be >= 3" << endl; - printUsage(argv[0]); - } - - sscanf(argv[2], "%d", &_M); - if (_M < _V) { - cout << "ERROR: (M,M) coordinates must be >= V" << endl; - printUsage(argv[0]); - } - - float _p; - sscanf(argv[3], "%f", &_p); - if (_p < 0 || _p >= 1) { - cout << "ERROR: probability p must be between [0, 1[" << endl; - printUsage(argv[0]); - } - randLimit = RAND_MAX * _p; - - if (argc == 4) { // Init rand seed - sscanf(argv[3], "%d", &seed); - srand(seed); - } else { // pseudo-rand seed - srand((unsigned int)time(NULL)); - } -} - -//------------------------------------------------------------------- -inline int randomValue(int max) { - return rand() % max; // [0, max - 1] -} - -inline void addEdge(map, unsigned int> &vertexPool, - int x0, int y0, int x1, int y1) { - map, unsigned int>::iterator it0 = - vertexPool.find({x0, y0}); - map, unsigned int>::iterator it1 = - vertexPool.find({x1, y1}); - if (it1->second < it0->second) - swap(it0, it1); // avoid repeated edges - - string key = to_string(it0->second) + " " + to_string(it1->second); - if (edgePool.find(key) != edgePool.end()) - return; // (Reversed?) edge already inserted - - // Compute distance/weight - unsigned int d = sqrt(pow((x0-x1),2) + pow((y0-y1),2)); - edgePool.insert({key, d}); -} - - -//------------------------------------------------------------------- -void generateGraph() { - // actual coordinates for Delaunay triangulation - vector coords; - - // to avoid repeated coordinates - map, unsigned int> vertexPool; - - // Create a pool of points on a (0,0) x (_M, _M) square - for (int i = 0; _currV < _V; i++) { - pair point = {randomValue(_M - 1), randomValue(_M - 1)}; - int id = ++_currV; - vertexPool.insert({point, id}); - coords.push_back(point.first); - coords.push_back(point.second); - if (i > 3 && (_V - _currV > 3) && rand() < randLimit) - break; // ensures enough vertices for a triangle - } - - // triangulation happens here - delaunator::Delaunator d(coords); - - // Extract edges from triangles (may be repeated) - for(size_t i = 0; i < d.triangles.size(); i += 3) { - int px0 = d.coords[2 * d.triangles[i]]; - int py0 = d.coords[2 * d.triangles[i]+1]; - //-- - int px1 = d.coords[2 * d.triangles[i+1]]; - int py1 = d.coords[2 * d.triangles[i+1]+1]; - //-- - int px2 = d.coords[2 * d.triangles[i+2]]; - int py2 = d.coords[2 * d.triangles[i+2]+1]; - //-- - addEdge(vertexPool, px0, py0, px1, py1); - addEdge(vertexPool, px1, py1, px2, py2); - addEdge(vertexPool, px2, py2, px0, py0); - } -} - -//------------------------------------------------------------------- -int main(int argc, char *argv[]) { - parseArgs(argc, argv); - _currV = 0; - - // Possibility to generate disconnected graph - while (_currV < _V) { - generateGraph(); - } - - // print (set of) graph(s) - cout << _V << endl << edgePool.size() << endl; - for (auto const& elem : edgePool) { - cout << elem.first << " " << elem.second << endl; - } - return 0; -} diff --git a/tests/v10.in b/tests/v10.in deleted file mode 100644 index 74c8157..0000000 --- a/tests/v10.in +++ /dev/null @@ -1,17 +0,0 @@ -10 -15 -1 3 34 -1 4 58 -1 5 28 -2 3 60 -2 4 7 -3 4 62 -3 5 24 -6 10 33 -6 7 50 -6 8 63 -7 10 57 -7 8 13 -7 9 37 -8 9 24 -9 10 90 diff --git a/tests/v15.in b/tests/v15.in deleted file mode 100644 index afdbea2..0000000 --- a/tests/v15.in +++ /dev/null @@ -1,26 +0,0 @@ -15 -24 -1 2 1472 -1 4 2148 -1 5 1108 -10 11 2288 -12 13 2422 -12 14 2715 -12 15 2004 -13 14 2483 -13 15 1996 -14 15 738 -2 3 2128 -2 5 1495 -3 4 2282 -3 5 942 -3 6 1524 -4 5 2035 -4 6 2162 -7 10 1974 -7 11 2273 -7 9 819 -8 10 2390 -8 11 1505 -9 10 1166 -9 11 1956 diff --git a/tests/v6.in b/tests/v6.in deleted file mode 100644 index 37cf81e..0000000 --- a/tests/v6.in +++ /dev/null @@ -1,12 +0,0 @@ -6 -10 -1 3 720 -1 5 498 -1 6 527 -2 3 1904 -2 4 1504 -2 6 1951 -3 4 2925 -3 5 1150 -3 6 833 -5 6 954