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);