地形篇-构网方法
地形里常见的构网方法就是delaunay方法了，而且该方法也较为重要，譬如地形成果（dem\等高线）、土方计算、以及断面线等重要成果都依赖于此，地形生产相关软件的构网方法均采用该方法，优点:算法较为成熟，网上资源很多，效率较高。其原理网上的教程很多，感兴趣的可以自行百度搜索。博主本人也查找过一些关于delaunay构网的方法，目前运用较多的是cgal,但cgal作为一个庞大的算法库配置较为麻烦、而且还涉及到商业协议所以该库的使用面临着诸多难题，所以也百度搜索过很多技术大牛提供的源代码，但是该类方法效率一般较低且不稳定。寻找一份高效、且方便移植的源代码变得很有必要。
有人可能会好奇：博主本人为啥不自己实现一个高效且稳定的delaunay方法，一方面水平所限吧；另一方面我所需要解决的核心问题与此无关（并不是我实现了该方法我所遇到的难题就迎刃而解了），所以没有花大量的精力来研究该问题。最近通过对一系列的资源搜索筛选，终于查找到一段高效的delaunay构网源码；所以喜出望外进行博客的整理。
实测1700w数据点，耗时15s;效率与cgal相当，可能更高；
但是对于一些较为复杂的构网方法（例如约束构网、金字塔模型）、甚至是三角网的处理还是得依赖一些成熟的三方库，不过对于这方面的资源自己日常也是在努力整合当中，毕竟这是一个系统性的问题，不是一朝一夕能轻松解决的。
源码如下，经本人测试，稳定可用。
#pragma once #ifdef DELAUNATOR_HEADER_ONLY #define INLINE inline #else #define INLINE #endif #if WINDOWS #undef max #endif // WINDOWS #include <limits> #include <vector> #include <stdint.h> namespace delaunator { // using shorter value type for indices rather than size_t // to consume less memory (it is unlikely that the triangulation // would be done on more than 4 billion points) typedef uint32_t index_t; constexpr index_t INVALID_INDEX = (std::numeric_limits<index_t>::max)(); class Point { public: Point(double x, double y) : m_x(x), m_y(y) {} Point() : m_x(0), m_y(0) {} double x() const { return m_x; } double y() const { return m_y; } private: double m_x; double m_y; }; class Delaunator { public: std::vector<double> const& coords; std::vector<index_t> triangles; std::vector<index_t> halfedges; std::vector<index_t> hull_prev; std::vector<index_t> hull_next; std::vector<index_t> hull_tri; index_t hull_start; INLINE Delaunator(std::vector<double> const& in_coords); INLINE double get_hull_area(); private: std::vector<index_t> m_hash; Point m_center; index_t m_hash_size; std::vector<index_t> m_edge_stack; INLINE index_t legalize(index_t a); INLINE index_t hash_key(double x, double y) const; INLINE index_t add_triangle( index_t i0, index_t i1, index_t i2, index_t a, index_t b, index_t c); INLINE void link(index_t a, index_t b); }; } //namespace delaunator #undef INLINE
#include "delaunator.hpp" #include <algorithm> #include <cmath> #include <limits> #include <stdexcept> #include <tuple> namespace delaunator { //@see https://stackoverflow.com/questions/33333363/built-in-mod-vs-custom-mod-function-improve-the-performance-of-modulus-op/33333636#33333636 inline index_t fast_mod(const index_t i, const index_t c) { return i >= c ? i % c : i; } // Kahan and Babuska summation, Neumaier variant; accumulates less FP error inline double sum(const std::vector<double>& x) { double sum = x[0]; double err = 0.0; for (index_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<double>::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 Point 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 Point(x, y); } 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<double>::epsilon(); 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 { index_t i; double x; double y; index_t t; index_t prev; index_t next; bool removed; }; Delaunator::Delaunator(std::vector<double> const& in_coords) : coords(in_coords), triangles(), halfedges(), hull_prev(), hull_next(), hull_tri(), hull_start(), m_hash(), m_hash_size(), m_edge_stack() { index_t n = coords.size() >> 1; double max_x = (std::numeric_limits<double>::min)(); double max_y = (std::numeric_limits<double>::min)(); double min_x = (std::numeric_limits<double>::max)(); double min_y = (std::numeric_limits<double>::max)(); std::vector<index_t> ids; ids.reserve(n); for (index_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<double>::max)(); index_t i0 = INVALID_INDEX; index_t i1 = INVALID_INDEX; index_t i2 = INVALID_INDEX; // pick a seed point close to the centroid for (index_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<double>::max)(); // find the point closest to the seed for (index_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<double>::max)(); // find the third point which forms the smallest circumcircle with the first two for (index_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<double>::max()))) { throw std::runtime_error("All points collinear"); } 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); } m_center = circumcenter(i0x, i0y, i1x, i1y, i2x, i2y); // Calculate the distances from the center once to avoid having to // calculate for each compare. This used to be done in the comparator, // but GCC 7.5+ would copy the comparator to iterators used in the // sort, and this was excruciatingly slow when there were many points // because you had to copy the vector of distances. std::vector<double> dists; dists.reserve(n); for (index_t i = 0; i < n; i++) { const double& x = coords[2 * i]; const double& y = coords[2 * i + 1]; dists.push_back(dist(x, y, m_center.x(), m_center.y())); } // sort the points by distance from the seed triangle circumcenter std::sort(ids.begin(), ids.end(), [&dists](index_t i, index_t j) { return dists[i] < dists[j]; }); // initialize a hash table for storing edges of the advancing convex hull m_hash_size = static_cast<index_t>(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; 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; index_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<double>::quiet_NaN(); double yp = std::numeric_limits<double>::quiet_NaN(); for (index_t k = 0; k < n; k++) { const index_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 index_t start = 0; index_t key = hash_key(x, y); for (index_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]; index_t e = start; index_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 index_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; // walk forward through the hull, adding more triangles and flipping recursively index_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 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 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<double> hull_area; index_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); } index_t Delaunator::legalize(index_t a) { index_t i = 0; index_t ar = 0; m_edge_stack.clear(); // recursion eliminated with a fixed-size stack while (true) { const index_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 index_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 index_t b0 = 3 * (b / 3); const index_t al = a0 + (a + 1) % 3; const index_t bl = b0 + (b + 2) % 3; const index_t p0 = triangles[ar]; const index_t pr = triangles[a]; const index_t pl = triangles[al]; const index_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) { index_t e = hull_start; do { if (hull_tri[e] == bl) { hull_tri[e] = a; break; } e = hull_prev[e]; } while (e != hull_start); } link(a, hbl); link(b, halfedges[ar]); link(ar, bl); index_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; } index_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<index_t>(std::llround(std::floor(pseudo_angle(dx, dy) * static_cast<double>(m_hash_size)))), m_hash_size); } index_t Delaunator::add_triangle( index_t i0, index_t i1, index_t i2, index_t a, index_t b, index_t c) { index_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 index_t a, const index_t b) { index_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) { index_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
