Implement Robust orientation

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

@@ -87,7 +87,7 @@ function _delaunator!(
     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 +157,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 +180,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 +193,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)

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