Skip to content

Commit 36df20a

Browse files
delfrrrdelfrrr
committed
Fix #9 Remove recursion from Delaunator::legalize
1 parent 393cfc3 commit 36df20a

File tree

4 files changed

+129
-101
lines changed

4 files changed

+129
-101
lines changed

‎.gitignore‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -14,4 +14,5 @@ node_modules
1414
mason_packages
1515
.toolchain
1616
.mason
17-
local.env
17+
local.env
18+
xcode-project

‎Makefile‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22
# Whether to turn compiler warnings into errors
33
export WERROR ?= true
44
export BUILD_DIR ?= cmake-build
5+
export XCODE_PROJ_DIR ?= xcode-project
56

67
default: release
78

@@ -11,6 +12,9 @@ release:
1112
debug:
1213
mkdir -p ./$(BUILD_DIR) && cd ./$(BUILD_DIR) && cmake ../ -DCMAKE_BUILD_TYPE=Debug -DWERROR=$(WERROR) && VERBOSE=1 cmake --build .
1314

15+
xcode:
16+
mkdir -p ./$(XCODE_PROJ_DIR) && cd ./$(XCODE_PROJ_DIR) && cmake -G Xcode ../
17+
1418
test:
1519
@if [ -f ./$(BUILD_DIR)/unit-tests ]; then ./$(BUILD_DIR)/unit-tests; else echo "Please run 'make release' or 'make debug' first" && exit 1; fi
1620

‎examples/utils.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ inline std::string read_file(const char* filename) {
2222
}
2323
}
2424

25-
inline std::vector<double> get_geo_json_points(std::string const& json) {
25+
inline std::vector< double> get_geo_json_points(std::string const& json) {
2626
rapidjson::Document document;
2727
if(document.Parse(json.c_str()).HasParseError()) {
2828
throw std::runtime_error("Cannot parse JSON");

‎include/delaunator.hpp‎

Lines changed: 122 additions & 99 deletions
Original file line numberDiff line numberDiff line change
@@ -93,47 +93,38 @@ inline std::pair<double, double> circumcenter(
9393
return std::make_pair(x, y);
9494
}
9595

96-
inline double compare(
97-
std::vector<double> const& coords,
98-
std::size_t i,
99-
std::size_t j,
100-
double cx,
101-
double cy) {
102-
const double d1 = dist(coords[2 * i], coords[2 * i + 1], cx, cy);
103-
const double d2 = dist(coords[2 * j], coords[2 * j + 1], cx, cy);
104-
const double diff1 = d1 - d2;
105-
const double diff2 = coords[2 * i] - coords[2 * j];
106-
const double diff3 = coords[2 * i + 1] - coords[2 * j + 1];
107-
108-
if (diff1 > 0.0 || diff1 < 0.0) {
109-
return diff1;
110-
} else if (diff2 > 0.0 || diff2 < 0.0) {
111-
return diff2;
112-
} else {
113-
return diff3;
114-
}
115-
}
116-
117-
struct sort_to_center {
96+
struct compare {
11897

11998
std::vector<double> const& coords;
12099
double cx;
121100
double cy;
122101

123102
bool operator()(std::size_t i, std::size_t j) {
124-
return compare(coords, i, j, cx, cy) < 0;
103+
const double d1 = dist(coords[2 * i], coords[2 * i + 1], cx, cy);
104+
const double d2 = dist(coords[2 * j], coords[2 * j + 1], cx, cy);
105+
const double diff1 = d1 - d2;
106+
const double diff2 = coords[2 * i] - coords[2 * j];
107+
const double diff3 = coords[2 * i + 1] - coords[2 * j + 1];
108+
109+
if (diff1 > 0.0 || diff1 < 0.0) {
110+
return diff1 < 0;
111+
} else if (diff2 > 0.0 || diff2 < 0.0) {
112+
return diff2 < 0;
113+
} else {
114+
return diff3 < 0;
115+
}
125116
}
126117
};
127118

128119
inline bool in_circle(
129-
double ax,
130-
double ay,
131-
double bx,
132-
double by,
133-
double cx,
134-
double cy,
135-
double px,
136-
double py) {
120+
const double ax,
121+
const double ay,
122+
const double bx,
123+
const double by,
124+
const double cx,
125+
const double cy,
126+
const double px,
127+
const double py) {
137128
const double dx = ax - px;
138129
const double dy = ay - py;
139130
const double ex = bx - px;
@@ -152,6 +143,7 @@ inline bool in_circle(
152143

153144
constexpr double EPSILON = std::numeric_limits<double>::epsilon();
154145
constexpr std::size_t INVALID_INDEX = std::numeric_limits<std::size_t>::max();
146+
constexpr std::size_t LEGALIZE_STACK_SIZE = 100;
155147

156148
inline bool check_pts_equal(double x1, double y1, double x2, double y2) {
157149
return std::fabs(x1 - x2) <= EPSILON &&
@@ -194,6 +186,7 @@ class Delaunator {
194186
double m_center_x;
195187
double m_center_y;
196188
std::size_t m_hash_size;
189+
std::size_t m_legalize_stack[LEGALIZE_STACK_SIZE];
197190

198191
std::size_t legalize(std::size_t a);
199192
std::size_t hash_key(double x, double y);
@@ -218,7 +211,8 @@ Delaunator::Delaunator(std::vector<double> const& in_coords)
218211
m_hash(),
219212
m_center_x(),
220213
m_center_y(),
221-
m_hash_size() {
214+
m_hash_size(),
215+
m_legalize_stack() {
222216
std::size_t n = coords.size() >> 1;
223217

224218
double max_x = std::numeric_limits<double>::min();
@@ -305,7 +299,7 @@ Delaunator::Delaunator(std::vector<double> const& in_coords)
305299
std::tie(m_center_x, m_center_y) = circumcenter(i0x, i0y, i1x, i1y, i2x, i2y);
306300

307301
// sort the points by distance from the seed triangle circumcenter
308-
std::sort(ids.begin(), ids.end(), sort_to_center{ coords, m_center_x, m_center_y });
302+
std::sort(ids.begin(), ids.end(), compare{ coords, m_center_x, m_center_y });
309303

310304
// initialize a hash table for storing edges of the advancing convex hull
311305
m_hash_size = static_cast<std::size_t>(std::llround(std::ceil(std::sqrt(n))));
@@ -439,75 +433,104 @@ double Delaunator::get_hull_area() {
439433
return sum(hull_area);
440434
}
441435

442-
std::size_t Delaunator::legalize(std::size_t a) {
443-
const std::size_t b = halfedges[a];
444-
445-
/* if the pair of triangles doesn't satisfy the Delaunay condition
446-
* (p1 is inside the circumcircle of [p0, pl, pr]), flip them,
447-
* then do the same check/flip recursively for the new pair of triangles
448-
*
449-
* pl pl
450-
* /||\ / \
451-
* al/ || \bl al/ \a
452-
* / || \ / \
453-
* / a||b \ flip /___ar___\
454-
* p0\ || /p1 => p0\---bl---/p1
455-
* \ || / \ /
456-
* ar\ || /br b\ /br
457-
* \||/ \ /
458-
* pr pr
459-
*/
460-
const std::size_t a0 = a - a % 3;
461-
const std::size_t b0 = b - b % 3;
462-
463-
const std::size_t al = a0 + (a + 1) % 3;
464-
const std::size_t ar = a0 + (a + 2) % 3;
465-
const std::size_t bl = b0 + (b + 2) % 3;
466-
467-
const std::size_t p0 = triangles[ar];
468-
const std::size_t pr = triangles[a];
469-
const std::size_t pl = triangles[al];
470-
const std::size_t p1 = triangles[bl];
471-
472-
if (b == INVALID_INDEX) {
473-
return ar;
474-
}
436+
std::size_t Delaunator::legalize(std::size_t ia) {
437+
std::size_t b;
438+
std::size_t a0;
439+
std::size_t b0;
440+
std::size_t al;
441+
std::size_t ar;
442+
std::size_t bl;
443+
444+
size_t i = 0;
445+
m_legalize_stack[i] = ia;
446+
size_t size = 1;
447+
while(i < size) {
448+
449+
if (i >= LEGALIZE_STACK_SIZE) {
450+
throw std::runtime_error("Legalize stack overflow");
451+
}
475452

476-
const bool illegal = in_circle(
477-
coords[2 * p0],
478-
coords[2 * p0 + 1],
479-
coords[2 * pr],
480-
coords[2 * pr + 1],
481-
coords[2 * pl],
482-
coords[2 * pl + 1],
483-
coords[2 * p1],
484-
coords[2 * p1 + 1]);
485-
486-
if (illegal) {
487-
triangles[a] = p1;
488-
triangles[b] = p0;
489-
490-
auto hbl = halfedges[bl];
491-
492-
// edge swapped on the other side of the hull (rare); fix the halfedge reference
493-
if (hbl == INVALID_INDEX) {
494-
std::size_t e = hull_start;
495-
do {
496-
if (hull_tri[e] == bl) {
497-
hull_tri[e] = a;
498-
break;
499-
}
500-
e = hull_next[e];
501-
} while (e != hull_start);
453+
auto const a = m_legalize_stack[i];
454+
i++;
455+
456+
/* if the pair of triangles doesn't satisfy the Delaunay condition
457+
* (p1 is inside the circumcircle of [p0, pl, pr]), flip them,
458+
* then do the same check/flip recursively for the new pair of triangles
459+
*
460+
* pl pl
461+
* /||\ / \
462+
* al/ || \bl al/ \a
463+
* / || \ / \
464+
* / a||b \ flip /___ar___\
465+
* p0\ || /p1 => p0\---bl---/p1
466+
* \ || / \ /
467+
* ar\ || /br b\ /br
468+
* \||/ \ /
469+
* pr pr
470+
*/
471+
472+
b = halfedges[a];
473+
474+
a0 = a - a % 3;
475+
b0 = b - b % 3;
476+
477+
al = a0 + (a + 1) % 3;
478+
ar = a0 + (a + 2) % 3;
479+
bl = b0 + (b + 2) % 3;
480+
481+
const std::size_t p0 = triangles[ar];
482+
const std::size_t pr = triangles[a];
483+
const std::size_t pl = triangles[al];
484+
const std::size_t p1 = triangles[bl];
485+
486+
if (b == INVALID_INDEX) {
487+
continue;
502488
}
503-
link(a, hbl);
504-
link(b, halfedges[ar]);
505-
link(ar, bl);
506489

507-
std::size_t br = b0 + (b + 1) % 3;
490+
const bool illegal = in_circle(
491+
coords[2 * p0],
492+
coords[2 * p0 + 1],
493+
coords[2 * pr],
494+
coords[2 * pr + 1],
495+
coords[2 * pl],
496+
coords[2 * pl + 1],
497+
coords[2 * p1],
498+
coords[2 * p1 + 1]);
499+
500+
if (illegal) {
501+
triangles[a] = p1;
502+
triangles[b] = p0;
503+
504+
auto hbl = halfedges[bl];
505+
506+
// edge swapped on the other side of the hull (rare); fix the halfedge reference
507+
if (hbl == INVALID_INDEX) {
508+
std::size_t e = hull_start;
509+
do {
510+
if (hull_tri[e] == bl) {
511+
hull_tri[e] = a;
512+
break;
513+
}
514+
e = hull_next[e];
515+
} while (e != hull_start);
516+
}
517+
link(a, hbl);
518+
link(b, halfedges[ar]);
519+
link(ar, bl);
520+
521+
std::size_t br = b0 + (b + 1) % 3;
508522

509-
legalize(a);
510-
return legalize(br);
523+
if (i < size) {
524+
//move elements down the stack
525+
for(auto mi = size - 1; mi >= i; mi--) {
526+
m_legalize_stack[mi + 2] = m_legalize_stack[mi];
527+
}
528+
}
529+
530+
m_legalize_stack[i] = a;
531+
m_legalize_stack[i + 1] = br;
532+
size+=2;
533+
}
511534
}
512535
return ar;
513536
}
@@ -537,7 +560,7 @@ std::size_t Delaunator::add_triangle(
537560
return t;
538561
}
539562

540-
void Delaunator::link(std::size_t a, std::size_t b) {
563+
void Delaunator::link(const std::size_t a, const std::size_t b) {
541564
std::size_t s = halfedges.size();
542565
if (a == s) {
543566
halfedges.push_back(b);

0 commit comments

Comments
 (0)