Fix #9 Remove recursion from Delaunator::legalize
diff --git a/.gitignore b/.gitignore index 7dde577..375f2a1 100644 --- a/.gitignore +++ b/.gitignore
@@ -14,4 +14,5 @@ mason_packages .toolchain .mason -local.env \ No newline at end of file +local.env +xcode-project \ No newline at end of file
diff --git a/Makefile b/Makefile index 5abaceb..8a0ac69 100644 --- a/Makefile +++ b/Makefile
@@ -2,6 +2,7 @@ # Whether to turn compiler warnings into errors export WERROR ?= true export BUILD_DIR ?= cmake-build +export XCODE_PROJ_DIR ?= xcode-project default: release @@ -11,6 +12,9 @@ debug: mkdir -p ./$(BUILD_DIR) && cd ./$(BUILD_DIR) && cmake ../ -DCMAKE_BUILD_TYPE=Debug -DWERROR=$(WERROR) && VERBOSE=1 cmake --build . +xcode: + mkdir -p ./$(XCODE_PROJ_DIR) && cd ./$(XCODE_PROJ_DIR) && cmake -G Xcode ../ + test: @if [ -f ./$(BUILD_DIR)/unit-tests ]; then ./$(BUILD_DIR)/unit-tests; else echo "Please run 'make release' or 'make debug' first" && exit 1; fi
diff --git a/examples/utils.hpp b/examples/utils.hpp index 02c32df..d22cb42 100644 --- a/examples/utils.hpp +++ b/examples/utils.hpp
@@ -22,7 +22,7 @@ } } -inline std::vector<double> get_geo_json_points(std::string const& json) { +inline std::vector< double> get_geo_json_points(std::string const& json) { rapidjson::Document document; if(document.Parse(json.c_str()).HasParseError()) { throw std::runtime_error("Cannot parse JSON");
diff --git a/include/delaunator.hpp b/include/delaunator.hpp index e16e0c3..11fb812 100644 --- a/include/delaunator.hpp +++ b/include/delaunator.hpp
@@ -93,47 +93,38 @@ return std::make_pair(x, y); } -inline double compare( - std::vector<double> const& coords, - std::size_t i, - std::size_t j, - double cx, - double cy) { - 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; - } else if (diff2 > 0.0 || diff2 < 0.0) { - return diff2; - } else { - return diff3; - } -} - -struct sort_to_center { +struct compare { std::vector<double> const& coords; double cx; double cy; bool operator()(std::size_t i, std::size_t j) { - return compare(coords, i, j, cx, cy) < 0; + 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( - double ax, - double ay, - double bx, - double by, - double cx, - double cy, - double px, - double py) { + 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; @@ -152,6 +143,7 @@ constexpr double EPSILON = std::numeric_limits<double>::epsilon(); constexpr std::size_t INVALID_INDEX = std::numeric_limits<std::size_t>::max(); +constexpr std::size_t LEGALIZE_STACK_SIZE = 100; inline bool check_pts_equal(double x1, double y1, double x2, double y2) { return std::fabs(x1 - x2) <= EPSILON && @@ -194,6 +186,7 @@ double m_center_x; double m_center_y; std::size_t m_hash_size; + std::size_t m_legalize_stack[LEGALIZE_STACK_SIZE]; std::size_t legalize(std::size_t a); std::size_t hash_key(double x, double y); @@ -218,7 +211,8 @@ m_hash(), m_center_x(), m_center_y(), - m_hash_size() { + m_hash_size(), + m_legalize_stack() { std::size_t n = coords.size() >> 1; double max_x = std::numeric_limits<double>::min(); @@ -305,7 +299,7 @@ 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(), sort_to_center{ coords, m_center_x, m_center_y }); + 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::size_t>(std::llround(std::ceil(std::sqrt(n)))); @@ -439,75 +433,104 @@ return sum(hull_area); } -std::size_t Delaunator::legalize(std::size_t a) { - const std::size_t b = halfedges[a]; +std::size_t Delaunator::legalize(std::size_t ia) { + std::size_t b; + std::size_t a0; + std::size_t b0; + std::size_t al; + std::size_t ar; + std::size_t bl; - /* 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 std::size_t a0 = a - a % 3; - const std::size_t b0 = b - b % 3; + size_t i = 0; + m_legalize_stack[i] = ia; + size_t size = 1; + while(i < size) { - const std::size_t al = a0 + (a + 1) % 3; - const std::size_t ar = a0 + (a + 2) % 3; - const std::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]; - - if (b == INVALID_INDEX) { - return ar; - } - - 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); + if (i >= LEGALIZE_STACK_SIZE) { + throw std::runtime_error("Legalize stack overflow"); } - link(a, hbl); - link(b, halfedges[ar]); - link(ar, bl); - std::size_t br = b0 + (b + 1) % 3; + auto const a = m_legalize_stack[i]; + i++; - legalize(a); - return legalize(br); + /* 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 + */ + + b = halfedges[a]; + + a0 = a - a % 3; + b0 = b - b % 3; + + al = a0 + (a + 1) % 3; + ar = a0 + (a + 2) % 3; + 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]; + + if (b == INVALID_INDEX) { + continue; + } + + 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 < size) { + //move elements down the stack + for(auto mi = size - 1; mi >= i; mi--) { + m_legalize_stack[mi + 2] = m_legalize_stack[mi]; + } + } + + m_legalize_stack[i] = a; + m_legalize_stack[i + 1] = br; + size+=2; + } } return ar; } @@ -537,7 +560,7 @@ return t; } -void Delaunator::link(std::size_t a, std::size_t b) { +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);