Implement Robust orientation (#18)

* Update for Robust Orientation
* Add new robust incircle test
* Add tests for robust codes
* Add check for julia version for supposition tests.
* Fix for tests on 1.6
* Only run supposition locally...
* Add performance test to look for regressions.

.github/workflows/CI.yml

@@ -1,4 +1,10 @@
 name: CI
+
+env:
+  OPENBLAS_NUM_THREADS: 1
+  OMP_NUM_THREADS: 1
+  JULIA_ACTIONS_RUNTEST_ARGS: CI
+
 on:
   pull_request:
     branches:
@@ -17,6 +23,7 @@ jobs:
         version:
           - '1.6'
           - '1.8'
+          - '1.10'
           #- 'nightly'
         os:
           - ubuntu-latest

Project.toml

@@ -1,26 +1,8 @@
 name = "Delaunator"
 uuid = "466f8f70-d5e3-4806-ac0b-a54b75a91218"
 authors = ["David Gleich", "Steve Kelly"]
-version = "0.1.1"
+version = "0.1.2"
 
 [compat]
-Aqua = "0.8"
-GeometryBasics = "0.4"
-JET = "0.4, 0.7, 0.8"
-JSON = "0.21"
-StableRNGs = "1"
-Test = "1"
-Pkg = "1"
-julia = "1"
+julia = "1.6"
 
-[extras]
-Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
-GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
-JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
-JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
-Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[targets]
-test = ["Aqua", "GeometryBasics", "JET", "JSON", "Pkg", "Test", "StableRNGs"]

README.jmd

@@ -86,7 +86,7 @@ rval = triangulate(pts)
 f = scatter(pts); hidespines!(f.axis); hidedecorations!(f.axis) 
 text!(f.axis, pts, text=map(i->"$i", 1:length(pts)))
 for i in eachindex(rval) 
-  p = Delaunator.dualcell(rval, i)
+  local p = Delaunator.dualcell(rval, i)
   # use clipped poly to get closed polygons
   # for the dualcells... 
   linesegments!(f.axis, collect(segments(p)),
@@ -159,6 +159,13 @@ might be suitable for those with computational geometry applications that need b
 accurate primitives, although these are relaxed slightly in the d3-delaunay usage -- especially in terms of the
 Voronoi cells. 
 
+Robust orientation
+==================
+The Delaunator.jl code direct includes the `robust_orient` function from 
+[`AdaptivePredicates.jl`](https://github.com/vchuravy/AdaptivePredicates.jl), which 
+ports the  
+[robust orientation routines from Jonathan Richard Shewchuk](https://www.cs.cmu.edu/~quake/robust.html)
+
 > This readme is auto-generated by weave from `README.jmd`
 ```julia; eval=false, echo=false 
 weave("README.jmd", doctype = "github", fig_path = "docs")

README.md

@@ -128,7 +128,7 @@ rval = triangulate(pts)
 f = scatter(pts); hidespines!(f.axis); hidedecorations!(f.axis) 
 text!(f.axis, pts, text=map(i->"$i", 1:length(pts)))
 for i in eachindex(rval) 
-  p = Delaunator.dualcell(rval, i)
+  local p = Delaunator.dualcell(rval, i)
   # use clipped poly to get closed polygons
   # for the dualcells... 
   linesegments!(f.axis, collect(segments(p)),
@@ -207,4 +207,11 @@ might be suitable for those with computational geometry applications that need b
 accurate primitives, although these are relaxed slightly in the d3-delaunay usage -- especially in terms of the
 Voronoi cells. 
 
+Robust orientation
+==================
+The Delaunator.jl code direct includes the `robust_orient` function from 
+[`AdaptivePredicates.jl`](https://github.com/vchuravy/AdaptivePredicates.jl), which 
+ports the  
+[robust orientation routines from Jonathan Richard Shewchuk](https://www.cs.cmu.edu/~quake/robust.html)
+
 > This readme is auto-generated by weave from `README.jmd`

src/Delaunator.jl

@@ -5,6 +5,7 @@ export triangulate, basictriangulation, update!, triangles, points, inhull
 
 include("quicksort.jl")
 include("geometry.jl")
+include("robust.jl")
 export isinfinite, dualcell, firstpoint, lastpoint, segments, cellarea
 
 include("clipping.jl")

src/algorithm.jl

@@ -79,15 +79,16 @@ function _delaunator!(
 ) where {FloatType}
     #FloatType = eltype(_dists)
     n = length(coords)
-    legalize = (t,_hullStart)->_legalize(FloatType, t, _triangles, _halfedges, coords, edgeStack, hullPrev, hullTri, _hullStart)
+    incircle_mempool = incircle_cache() 
+    legalize = (t,_hullStart)->_legalize(FloatType, t, _triangles, _halfedges, coords, edgeStack, hullPrev, hullTri, _hullStart, incircle_mempool)
 
     i0, i1, i2 = seeds
     i0x,i0y = point(FloatType,coords,i0)
     i1x,i1y = point(FloatType,coords,i1)
     i2x,i2y = point(FloatType,coords,i2)
 
     # swap the order of the seed points for counter-clockwise orientation
-    if orient(i0x, i0y, i1x, i1y, i2x, i2y)
+    if robust_orient(i0x, i0y, i1x, i1y, i2x, i2y)
         i = i1
         x = i1x
         y = i1y
@@ -157,7 +158,7 @@ function _delaunator!(
         start = hullPrev[start]
         e = start
         q = hullNext[e]
-        while !orient(x, y, point(FloatType,coords,e)..., point(FloatType,coords,q)...)
+        while !robust_orient(x, y, point(FloatType,coords,e)..., point(FloatType,coords,q)...)
             e = q
             if e == start
                 e = -1
@@ -180,7 +181,7 @@ function _delaunator!(
         # walk forward through the hull, adding more triangles and flipping recursively
         n = hullNext[e]
         q = hullNext[n]
-        while orient(x, y, point(FloatType,coords,n)..., point(FloatType,coords,q)...)
+        while robust_orient(x, y, point(FloatType,coords,n)..., point(FloatType,coords,q)...)
             trianglesLen = _addTriangle(_triangles, _halfedges, trianglesLen, n, i, q, hullTri[i], -1, hullTri[n])
             t = trianglesLen-3
             hullTri[i] = legalize(t + 2, _hullStart)
@@ -193,7 +194,7 @@ function _delaunator!(
         # walk backward from the other side, adding more triangles and flipping
         if e == start
             q = hullPrev[e]
-            while orient(x, y, point(FloatType,coords,q)...,  point(FloatType,coords,e)...)
+            while robust_orient(x, y, point(FloatType,coords,q)...,  point(FloatType,coords,e)...)
                 trianglesLen = _addTriangle(_triangles, _halfedges, trianglesLen, q, i, e, -1, hullTri[e], hullTri[q])
                 t = trianglesLen-3
                 legalize(t + 2, _hullStart)
@@ -252,7 +253,7 @@ function _hashKey(x, y, cx, cy, _hashSize)
     return (floor(Int, pseudoAngle(x - cx, y - cy) * _hashSize) % _hashSize)+1
 end
 
-function _legalize(::Type{FloatType}, a, triangles, halfedges, coords, EDGE_STACK, hullPrev, hullTri, _hullStart) where {FloatType}
+function _legalize(::Type{FloatType}, a, triangles, halfedges, coords, EDGE_STACK, hullPrev, hullTri, _hullStart, mempool) where {FloatType}
 
     i = 1
     ar = 0
@@ -296,10 +297,11 @@ function _legalize(::Type{FloatType}, a, triangles, halfedges, coords, EDGE_STAC
         pl = triangles[al]
         p1 = triangles[bl]
 
-        illegal = inCircle(point(FloatType,coords,p0)..., 
-                           point(FloatType,coords,pr)..., 
-                           point(FloatType,coords,pl)..., 
-                           point(FloatType,coords,p1)...)
+        illegal = incircle(point(FloatType,coords,p0), 
+                           point(FloatType,coords,pr),
+                           point(FloatType,coords,pl),
+                           point(FloatType,coords,p1), 
+                           mempool) < 0 
 
         if illegal
             triangles[a] = p1

src/geometry.jl

@@ -84,7 +84,7 @@ function _isright(pa, pb, pt)
     return (x2 - x1)*(y - y1) < (y2 - y1)*(x - x1)
 end 
 function _isleft(pa, pb, pt)
-    return orient(pa...,pb...,pt...)
+    return robust_orient(pa...,pb...,pt...)
 end 
 
 
@@ -302,6 +302,7 @@ function orient(rx, ry, qx, qy, px, py)
            orientIfSure(qx, qy, px, py, rx, ry) < 0
 end
 
+
 function inCircle(ax, ay, bx, by, cx, cy, px, py)
     dx = ax - px
     dy = ay - py
@@ -319,6 +320,7 @@ function inCircle(ax, ay, bx, by, cx, cy, px, py)
            ap * (ex * fy - ey * fx) < 0
 end
 
+
 function circumradius(ax, ay, bx, by, cx, cy)
     dx = bx - ax
     dy = by - ay