added more robust orient test

src/org/poly2tri/Poly2Tri.scala

@@ -279,12 +279,12 @@ class Poly2TriDemo extends BasicGame("Poly2Tri") {
   def selectModel(model: String) {
     model match {
       case "data/nazca_monkey.dat" => 
-        //CDT.clearPoint = 50
+        CDT.clearPoint = 50
-        doCDT = false; drawCDT = false; drawcdtMesh = false
+        //doCDT = false; drawCDT = false; drawcdtMesh = false
         loadModel(nazcaMonkey, 4.5f, Point(400, 300), 1500)
       case "data/bird.dat" => 
-        doCDT = false; drawCDT = false; drawcdtMesh = false
+        //doCDT = false; drawCDT = false; drawcdtMesh = false
-        //CDT.clearPoint = 80
+        CDT.clearPoint = 80
         loadModel(bird, 25f, Point(400, 300), 350)
       case "data/i.snake" => 
         doCDT = true; drawCDT = true

src/org/poly2tri/cdt/CDT.scala

@@ -133,6 +133,8 @@ class CDT(val points: List[Point], val segments: List[Segment], iTriangle: Trian
   private def sweep {
     
     for(i <- 1 until points.size) {
+      
+      try {
       val point = points(i)
       // Process Point event
       var triangle = pointEvent(point)
@@ -140,8 +142,13 @@ class CDT(val points: List[Point], val segments: List[Segment], iTriangle: Trian
       if(triangle != null)
         point.edges.foreach(e => triangle = edgeEvent(e, triangle))
       if(i == CDT.clearPoint) {cleanTri = triangle; mesh.debug += cleanTri}
-    }
       
+      } catch {
+        case e: Exception =>
+          throw new Exception("Systect point = " + i)
+      }
+      
+    }
   }  
   
   // Final step in the sweep-line CDT algo
@@ -208,7 +215,7 @@ class CDT(val points: List[Point], val segments: List[Segment], iTriangle: Trian
 		tList.foreach(t => {
 		  t.points.foreach(p => {
 		    if(p != edge.q && p != edge.p) {
-		      if(t.orient(edge.q, edge.p, p) > 0 ) {
+		      if(Util.orient2d(edge.q, edge.p, p) > 0 ) {
 		        // Keep duplicate points out
                 if(!lPoints.contains(p)) {
 		          lPoints += p
@@ -347,7 +354,7 @@ class CDT(val points: List[Point], val segments: List[Segment], iTriangle: Trian
     } 
     
     if(!P.isEmpty) {
-      val left = Util.orient(a, b, P(i)) < 0
+      val left = Util.orient2d(a, b, P(i)) > 0
       val pB = if(left) P(i) else b
       val pC = if(left) b else P(i)
       val points = Array(a, pB, pC)
@@ -434,6 +441,7 @@ class CDT(val points: List[Point], val segments: List[Segment], iTriangle: Trian
     val point = t1.points(0)
     val oPoint = t2 oppositePoint t1
     
+    // Try to avoid creating degenerate triangles
     val collinear = t1.collinear(oPoint) || t2.collinear(oPoint, point)
     
     if(illegal(t1.points(1), oPoint, t1.points(2), t1.points(0)) && !collinear) {

src/org/poly2tri/shapes/Triangle.scala

@@ -126,23 +126,25 @@ class Triangle(val points: Array[Point]) {
     val p = edge.p
     val q = edge.q
     val e = p - q
+    
+    val ik = points(2) - points(0)
+    val ij = points(1) - points(0)
+    val jk = points(2) - points(1)
+    val ji = points(0) - points(1)
+    val kj = points(1) - points(2)
+    val ki = points(0) - points(2)
+    
     if(q == points(0)) {
-      val ik = points(2) - points(0)
-      val ij = points(1) - points(0)
       val sameSign = Math.signum(ik cross e) == Math.signum(ij cross e)
       if(!sameSign) return this
       if(neighbors(2) == null) return null
       return neighbors(2).locateFirst(edge)
     } else if(q == points(1)) {
-      val jk = points(2) - points(1)
-      val ji = points(0) - points(1)
       val sameSign = Math.signum(jk cross e) == Math.signum(ji cross e)
       if(!sameSign) return this
       if(neighbors(0) == null) return null
       return neighbors(0).locateFirst(edge)
     } else if(q == points(2)) {
-      val kj = points(1) - points(2)
-      val ki = points(0) - points(2)
       val sameSign = Math.signum(kj cross e) == Math.signum(ki cross e)
       if(!sameSign) return this
       if(neighbors(1) == null) return null
@@ -153,29 +155,17 @@ class Triangle(val points: Array[Point]) {
   
   // Locate next triangle crossed by edge
   def findNeighbor(e: Point): Triangle = {
-    if(orient(points(0), points(1), e) >= 0)
+    if(Util.orient2d(points(0), points(1), e) > 0)
       return neighbors(2)
-    else if(orient(points(1), points(2), e) >= 0)
+    else if(Util.orient2d(points(1), points(2), e) > 0)
       return neighbors(0)
-    else if(orient(points(2), points(0), e) >= 0)
+    else if(Util.orient2d(points(2), points(0), e) > 0)
       return neighbors(1)
     else
       // Point must reside inside this triangle
       this
   }
   
-  // Return: positive if point p is left of ab
-  //         negative if point p is right of ab
-  //         zero if points are colinear
-  // See: http://www-2.cs.cmu.edu/~quake/robust.html
-  def orient(b: Point, a: Point, p: Point): Float = {
-    val acx = a.x - p.x
-    val bcx = b.x - p.x
-    val acy = a.y - p.y
-    val bcy = b.y - p.y
-    acx * bcy - acy * bcx  
-  }
-  
   // The neighbor clockwise to given point
   def neighborCW(point: Point): Triangle = {
     if(point == points(0)) {
@@ -260,7 +250,7 @@ class Triangle(val points: Array[Point]) {
   }
   
   // Make legalized triangle will not be collinear
-  def collinear(oPoint: Point): Boolean = Util.collinear(points(0), points(1), oPoint)
+  def collinear(oPoint: Point): Boolean = Util.collinear(points(1), points(0), oPoint)
   
   // Make sure legalized triangle will not be collinear
   def collinear(oPoint: Point, nPoint: Point): Boolean = {

src/org/poly2tri/utils/Util.scala

@@ -8,7 +8,10 @@ import shapes.Point
 object Util {
   
   // Almost zero
-  val COLLINEAR_SLOP = 0.1f
+  val COLLINEAR_SLOP = 0.25f
+  
+  val epsilon = exactinit
+  val ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon
 
   // From "Scala By Example," by Martin Odersky
   def msort[A](less: (A, A) => Boolean)(xs: List[A]): List[A] = {
@@ -49,9 +52,37 @@ object Util {
     
   }
   
-  // Return: positive if point p is left of ab
+  /* From Jonathan Shewchuk's "Adaptive Precision Floating-Point Arithmetic 
-  //         negative if point p is right of ab
+   * and Fast Robust Predicates for Computational Geometry"
-  //         zero if points are colinear
+   * See: http://www.cs.cmu.edu/~quake/robust.html
+   */
+  def exactinit = {
+  
+	  var every_other = true
+	  var half = 0.5
+      var splitter = 1.0
+	  var epsilon = 1.0
+	  var check = 1.0
+	  var lastcheck = 0.0
+   
+	  do {
+	    lastcheck = check
+	    epsilon *= half
+	    if (every_other) {
+	      splitter *= 2.0
+	    }
+	    every_other = !every_other
+	    check = 1.0 + epsilon
+	  } while ((check != 1.0) && (check != lastcheck))
+
+	  epsilon
+  
+  }
+  
+  // Approximate 2D orientation test.  Nonrobust. 
+  // Return: positive if point  a, b, and c are counterclockwise
+  //         negative if point a, b, and c are clockwise
+  //         zero if points are collinear
   // See: http://www-2.cs.cmu.edu/~quake/robust.html
   def orient(b: Point, a: Point, p: Point): Float = {
     val acx = a.x - p.x
@@ -61,6 +92,49 @@ object Util {
     acx * bcy - acy * bcx  
   }
    
+   // Adaptive exact 2D orientation test.  Robust. By Jonathan Shewchuk
+   // Return: positive if point  a, b, and c are counterclockwise
+   //         negative if point a, b, and c are clockwise
+   //         zero if points are collinear
+   // See: http://www-2.cs.cmu.edu/~quake/robust.html
+   def orient2d(pa: Point, pb: Point, pc: Point): Double = {
+  
+	   val detleft: Double = (pa.x - pc.x) * (pb.y - pc.y)
+	   val detright: Double = (pa.y - pc.y) * (pb.x - pc.x)
+	   val det = detleft - detright
+	   var detsum = 0.0
+	  
+	   if (detleft > 0.0) {
+	     if (detright <= 0.0) {
+	       return det;
+	     } else {
+	       detsum = detleft + detright
+	     }
+	   } else if (detleft < 0.0) {
+	     if (detright >= 0.0) {
+	       return det
+	     } else {
+	       detsum = -detleft - detright
+	     }
+	   } else {
+	     return det
+	   }
+	
+	   val errbound = Util.ccwerrboundA * detsum
+	   if ((det >= errbound) || (-det >= errbound)) {
+	     return det
+	   } else {
+	     println(pa + "," + pb + "," + pc)
+	     println(detleft + "," + detright)
+	     println("Det = " + det + " , errbound = " + errbound)
+	     throw new Exception("Degenerate triangle")
+	   }
+	   
+	  // Not implemented; 
+	  // http://www.cs.cmu.edu/afs/cs/project/quake/public/code/predicates.c
+	  // return orient2dadapt(pa, pb, pc, detsum)
+ }
+ 
 }
 
 /** The object <code>Random</code> offers a default implementation