From 260e7935b9651406b01f26ea4537bc41aa34460c Mon Sep 17 00:00:00 2001 From: zzzzrrr Date: Tue, 11 Aug 2009 11:29:33 -0400 Subject: [PATCH] added more robust orient test --- src/org/poly2tri/Poly2Tri.scala | 8 +-- src/org/poly2tri/cdt/CDT.scala | 14 ++++- src/org/poly2tri/shapes/Triangle.scala | 34 ++++------- src/org/poly2tri/utils/Util.scala | 82 ++++++++++++++++++++++++-- 4 files changed, 105 insertions(+), 33 deletions(-) diff --git a/src/org/poly2tri/Poly2Tri.scala b/src/org/poly2tri/Poly2Tri.scala index da29077..4e91765 100644 --- a/src/org/poly2tri/Poly2Tri.scala +++ b/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 - doCDT = false; drawCDT = false; drawcdtMesh = false + CDT.clearPoint = 50 + //doCDT = false; drawCDT = false; drawcdtMesh = false loadModel(nazcaMonkey, 4.5f, Point(400, 300), 1500) case "data/bird.dat" => - doCDT = false; drawCDT = false; drawcdtMesh = false - //CDT.clearPoint = 80 + //doCDT = false; drawCDT = false; drawcdtMesh = false + CDT.clearPoint = 80 loadModel(bird, 25f, Point(400, 300), 350) case "data/i.snake" => doCDT = true; drawCDT = true diff --git a/src/org/poly2tri/cdt/CDT.scala b/src/org/poly2tri/cdt/CDT.scala index 488d8a0..830aea9 100644 --- a/src/org/poly2tri/cdt/CDT.scala +++ b/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) { diff --git a/src/org/poly2tri/shapes/Triangle.scala b/src/org/poly2tri/shapes/Triangle.scala index ddc23f5..b78dadf 100644 --- a/src/org/poly2tri/shapes/Triangle.scala +++ b/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 = { diff --git a/src/org/poly2tri/utils/Util.scala b/src/org/poly2tri/utils/Util.scala index a24ba27..bbb72c9 100644 --- a/src/org/poly2tri/utils/Util.scala +++ b/src/org/poly2tri/utils/Util.scala @@ -8,8 +8,11 @@ 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] = { def merge(xs1: List[A], xs2: List[A]): List[A] = @@ -49,9 +52,37 @@ object Util { } - // Return: positive if point p is left of ab - // negative if point p is right of ab - // zero if points are colinear + /* From Jonathan Shewchuk's "Adaptive Precision Floating-Point Arithmetic + * and Fast Robust Predicates for Computational Geometry" + * 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 @@ -60,7 +91,50 @@ object Util { val bcy = b.y - p.y 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 Random offers a default implementation