fixed numerical bugs

This commit is contained in:
zzzzrrr 2009-11-20 12:13:19 -05:00
parent 93e9955a66
commit 2e0041456a
5 changed files with 105 additions and 85 deletions

View File

@ -1,5 +1,8 @@
100.0 100.0 474.80999000000003 555.15656999999999
-100.0 100.0 474.80999000000003 530.87086999999997
-100.0 -100.0 509.09570000000002 530.87086999999997
0.0 0.0 543.38142000000005 530.87086999999997
100.0 -100.0 543.38142000000005 555.15656999999999
543.38142000000005 579.44227000000001
509.09570000000002 579.44227000000001
474.80999000000003 579.44227000000001

View File

@ -50,6 +50,14 @@ def draw_polygon(verts, color):
for v in verts: for v in verts:
glVertex2f(v[0], v[1]) glVertex2f(v[0], v[1])
glEnd() glEnd()
def draw_line(p1, p2, color):
r, g, b = color
glColor3f(r, g, b)
glBegin(GL_LINES)
glVertex2f(p1[0], p1[1])
glVertex2f(p2[0], p2[1])
glEnd()
## ##
## Game engine / main loop / UI ## Game engine / main loop / UI

View File

@ -1,5 +1,5 @@
#!/usr/bin/env python2.6 #!/usr/bin/env python2.6
from framework import Game, draw_polygon, reset_zoom from framework import Game, draw_polygon, reset_zoom, draw_line
from seidel import Triangulator from seidel import Triangulator
@ -11,7 +11,7 @@ class Poly2Tri(Game):
super(Poly2Tri, self).__init__(*self.screen_size) super(Poly2Tri, self).__init__(*self.screen_size)
# Load point set # Load point set
file_name = "../data/nazca_monkey.dat" file_name = "../data/star.dat"
self.points = self.load_points(file_name) self.points = self.load_points(file_name)
# Triangulate # Triangulate
@ -20,8 +20,9 @@ class Poly2Tri(Game):
dt = (self.time - t1) * 1000.0 dt = (self.time - t1) * 1000.0
self.triangles = seidel.triangles() self.triangles = seidel.triangles()
self.trapezoids = seidel.trapezoids #self.trapezoids = seidel.trapezoids
#self.trapezoids = seidel.trapezoidal_map.map self.trapezoids = seidel.trapezoidal_map.map
self.edges = seidel.edge_list
print "time (ms) = %f , num triangles = %d" % (dt, len(self.triangles)) print "time (ms) = %f , num triangles = %d" % (dt, len(self.triangles))
self.main_loop() self.main_loop()
@ -30,19 +31,26 @@ class Poly2Tri(Game):
pass pass
def render(self): def render(self):
reset_zoom(7, (0, 0), self.screen_size) reset_zoom(2.1, (400, 100), self.screen_size)
red = 255, 0, 0 red = 255, 0, 0
yellow = 255, 255, 0
green = 0, 255, 0
for t in self.triangles: for t in self.triangles:
draw_polygon(t, red) draw_polygon(t, red)
green = 0, 255, 0
draw_polygon(self.points, green)
''' '''
yellow = 255, 255, 0
for t in self.trapezoids: for t in self.trapezoids:
#verts = self.trapezoids[key].vertices() verts = self.trapezoids[t].vertices()
verts = t.vertices() #verts = t.vertices()
draw_polygon(verts, yellow) draw_polygon(verts, yellow)
''' '''
for e in self.edges:
p1 = e.p.x, e.p.y
p2 = e.q.x, e.q.y
draw_line(p1, p2, green)
#draw_polygon(self.points, green)
def load_points(self, file_name): def load_points(self, file_name):
infile = open(file_name, "r") infile = open(file_name, "r")
@ -55,7 +63,8 @@ class Poly2Tri(Game):
points.append((float(s[0]), float(s[1]))) points.append((float(s[0]), float(s[1])))
return points return points
bad = (544.80998999999997, 579.86046999999996), (544.80998999999997, 450.57477), (594.09569999999997, 450.57477), (643.38142000000005, 450.57477), (643.38142000000005, 525.26486999999997), (643.38142000000005, 599.95487000000003), (603.67391999999995, 654.55056999999999), (563.96655999999996, 709.14621999999997), (554.38819999999998, 709.14621999999997), (544.80998999999997, 709.14621999999997) spam = (544.80998999999997, 579.86046999999996), (544.80998999999997, 450.57477), (594.09569999999997, 450.57477), (643.38142000000005, 450.57477), (643.38142000000005, 525.26486999999997), (643.38142000000005, 599.95487000000003), (603.67391999999995, 654.55056999999999), (563.96655999999996, 709.14621999999997), (554.38819999999998, 709.14621999999997), (544.80998999999997, 709.14621999999997)
eggs = [474.80999000000003, 555.15656999999999], [474.80999000000003, 530.87086999999997], [509.09570000000002, 530.87086999999997], [543.38142000000005, 530.87086999999997], [543.38142000000005, 555.15656999999999], [543.38142000000005, 579.44227000000001], [509.09570000000002, 579.44227000000001], [474.80999000000003, 579.44227000000001]
if __name__ == '__main__': if __name__ == '__main__':
demo = Poly2Tri() demo = Poly2Tri()

View File

@ -38,7 +38,7 @@ from math import atan2
## (Ported from poly2tri) ## (Ported from poly2tri)
## ##
# Shear transform. Can effect numerical robustness # Shear transform. May effect numerical robustness
SHEAR = 1e-6 SHEAR = 1e-6
class Point(object): class Point(object):
@ -80,6 +80,9 @@ class Point(object):
def less(self, p): def less(self, p):
return self.x < p.x return self.x < p.x
def neq(self, other):
return other.x != self.x or other.y != self.y
def clone(self): def clone(self):
return Point(self.x, self.y) return Point(self.x, self.y)
@ -89,20 +92,22 @@ class Edge(object):
def __init__(self, p, q): def __init__(self, p, q):
self.p = p self.p = p
self.q = q self.q = q
self.slope = 0.0 if (q.x - p.x) == 0.0 else (q.y - p.y)/(q.x - p.x) self.slope = (q.y - p.y) / (q.x - p.x)
self.b = p.y - (p.x * self.slope) self.b = p.y - (p.x * self.slope)
self.above, self.below = None, None self.above, self.below = None, None
self.mpoints = [] self.mpoints = []
self.mpoints.append(p) self.mpoints.append(p)
self.mpoints.append(q) self.mpoints.append(q)
##
## NOTE Rounding accuracy significantly effects numerical robustness!!!
##
def is_above(self, point): def is_above(self, point):
# NOTE rounding can effect numerical robustness return (round(point.y, 2) < round(self.slope * point.x + self.b, 2))
return (round(point.y, 10) < round(self.slope * point.x + self.b, 10))
def is_below(self, point): def is_below(self, point):
# NOTE rounding can effect numerical robustness return (round(point.y, 2) > round(self.slope * point.x + self.b, 2))
return (round(point.y, 10) > round(self.slope * point.x + self.b, 10))
def intersect(self, c, d): def intersect(self, c, d):
a = self.p a = self.p
@ -127,13 +132,13 @@ class Trapezoid(object):
self.right_point = right_point self.right_point = right_point
self.top = top self.top = top
self.bottom = bottom self.bottom = bottom
self.hash = hash(self)
self.upper_left = None self.upper_left = None
self.upper_right = None self.upper_right = None
self.lower_left = None self.lower_left = None
self.lower_right = None self.lower_right = None
self.inside = True self.inside = True
self.sink = None self.sink = None
self.key = hash(self)
def update_left(self, ul, ll): def update_left(self, ul, ll):
self.upper_left = ul self.upper_left = ul
@ -193,6 +198,7 @@ def line_intersect(edge, x):
class Triangulator(object): class Triangulator(object):
def __init__(self, poly_line): def __init__(self, poly_line):
assert len(poly_line) > 3, "Number of points must be > 3"
self.polygons = [] self.polygons = []
self.trapezoids = [] self.trapezoids = []
self.xmono_poly = [] self.xmono_poly = []
@ -200,7 +206,7 @@ class Triangulator(object):
self.trapezoidal_map = TrapezoidalMap() self.trapezoidal_map = TrapezoidalMap()
self.bounding_box = self.trapezoidal_map.bounding_box(self.edge_list) self.bounding_box = self.trapezoidal_map.bounding_box(self.edge_list)
self.query_graph = QueryGraph(isink(self.bounding_box)) self.query_graph = QueryGraph(isink(self.bounding_box))
self.process() self.process()
def triangles(self): def triangles(self):
@ -221,7 +227,7 @@ class Triangulator(object):
traps = self.query_graph.follow_edge(edge) traps = self.query_graph.follow_edge(edge)
for t in traps: for t in traps:
# Remove old trapezods # Remove old trapezods
del self.trapezoidal_map.map[t.hash] del self.trapezoidal_map.map[t.key]
# Bisect old trapezoids and create new # Bisect old trapezoids and create new
cp = t.contains(edge.p) cp = t.contains(edge.p)
cq = t.contains(edge.q) cq = t.contains(edge.q)
@ -239,12 +245,9 @@ class Triangulator(object):
self.query_graph.case4(t.sink, edge, tlist) self.query_graph.case4(t.sink, edge, tlist)
# Add new trapezoids to map # Add new trapezoids to map
for t in tlist: for t in tlist:
self.trapezoidal_map.map[t.hash] = t self.trapezoidal_map.map[t.key] = t
self.trapezoidal_map.clear() self.trapezoidal_map.clear()
#TODO remove invalid/extra trapezoids
#print len(self.trapezoidal_map.map)
# Mark outside trapezoids w/ depth-first search # Mark outside trapezoids w/ depth-first search
for k, t in self.trapezoidal_map.map.items(): for k, t in self.trapezoidal_map.map.items():
self.mark_outside(t) self.mark_outside(t)
@ -281,30 +284,30 @@ class Triangulator(object):
t.trim_neighbors() t.trim_neighbors()
def init_edges(self, points): def init_edges(self, points):
edges = [] edge_list = []
for i in range(len(points)-1): size = len(points)
p = Point(points[i][0], points[i][1]) for i in range(size):
q = Point(points[i+1][0], points[i+1][1]) j = i + 1 if i < size-1 else 0
edges.append(Edge(p, q)) p = points[i][0], points[i][1]
p = Point(points[0][0], points[0][1]) q = points[j][0], points[j][1]
q = Point(points[-1][0], points[-1][1]) edge_list.append((p, q))
edges.append(Edge(p, q)) return self.order_edges(edge_list)
return self.order_edges(edges)
def order_edges(self, edge_list): def order_edges(self, edge_list):
edges = [] edges = []
for e in edge_list: for e in edge_list:
p = self.shear_transform(e.p) p = shear_transform(e[0])
q = self.shear_transform(e.q) q = shear_transform(e[1])
if p.x > q.x: if p.x > q.x:
edges.append(Edge(q, p)) edges.append(Edge(q, p))
elif p.x < q.x: else:
edges.append(Edge(p, q)) edges.append(Edge(p, q))
# Randomized incremental algorithm
shuffle(edges) shuffle(edges)
return edges return edges
def shear_transform(self, point): def shear_transform(point):
return Point(point.x + SHEAR * point.y, point.y) return Point(point[0] + SHEAR * point[1], point[1])
def merge_sort(l): def merge_sort(l):
if len(l)>1 : if len(l)>1 :
@ -422,29 +425,29 @@ class TrapezoidalMap(object):
if e.q.y < min.y: min = Point(min.x, e.q.y - margin) if e.q.y < min.y: min = Point(min.x, e.q.y - margin)
top = Edge(Point(min.x, max.y), Point(max.x, max.y)) top = Edge(Point(min.x, max.y), Point(max.x, max.y))
bottom = Edge(Point(min.x, min.y), Point(max.x, min.y)) bottom = Edge(Point(min.x, min.y), Point(max.x, min.y))
left = bottom.p left = top.p
right = top.q right = top.q
trap = Trapezoid(left, right, top, bottom) trap = Trapezoid(left, right, top, bottom)
self.map[hash(trap)] = trap self.map[trap.key] = trap
return trap return trap
class Node(object): class Node(object):
def __init__(self, left, right): def __init__(self, lchild, rchild):
self.parent_list = [] self.parent_list = []
self.left = left self.lchild = lchild
self.right = right self.rchild = rchild
if left != None: if lchild != None:
left.parent_list.append(self) lchild.parent_list.append(self)
if right != None: if rchild != None:
right.parent_list.append(self) rchild.parent_list.append(self)
def replace(self, node): def replace(self, node):
for parent in node.parent_list: for parent in node.parent_list:
if parent.left is node: if parent.lchild is node:
parent.left = self parent.lchild = self
else: else:
parent.right = self parent.rchild = self
self.parent_list += node.parent_list self.parent_list += node.parent_list
class Sink(Node): class Sink(Node):
@ -458,9 +461,9 @@ class Sink(Node):
return self return self
def isink(trapezoid): def isink(trapezoid):
if trapezoid.sink != None: if trapezoid.sink is None:
return trapezoid.sink return Sink(trapezoid)
return Sink(trapezoid) return trapezoid.sink
class XNode(Node): class XNode(Node):
@ -470,8 +473,8 @@ class XNode(Node):
def locate(self, edge): def locate(self, edge):
if edge.p.x >= self.point.x: if edge.p.x >= self.point.x:
return self.right.locate(edge) return self.rchild.locate(edge)
return self.left.locate(edge) return self.lchild.locate(edge)
class YNode(Node): class YNode(Node):
@ -481,14 +484,12 @@ class YNode(Node):
def locate(self, edge): def locate(self, edge):
if self.edge.is_above(edge.p): if self.edge.is_above(edge.p):
return self.right.locate(edge) return self.rchild.locate(edge)
elif self.edge.is_below(edge.p): if self.edge.is_below(edge.p):
return self.left.locate(edge) return self.lchild.locate(edge)
else: if edge.slope < self.edge.slope:
if edge.slope < self.edge.slope: return self.rchild.locate(edge)
return self.right.locate(edge) return self.lchild.locate(edge)
else:
return self.left.locate(edge)
class QueryGraph: class QueryGraph:
@ -500,20 +501,18 @@ class QueryGraph:
def follow_edge(self, edge): def follow_edge(self, edge):
trapezoids = [self.locate(edge)] trapezoids = [self.locate(edge)]
j = 0 while(edge.q.x > trapezoids[-1].right_point.x):
while(edge.q.x > trapezoids[j].right_point.x): if edge.is_above(trapezoids[-1].right_point):
if edge.is_above(trapezoids[j].right_point): trapezoids.append(trapezoids[-1].upper_right)
trapezoids.append(trapezoids[j].upper_right)
else: else:
trapezoids.append(trapezoids[j].lower_right) trapezoids.append(trapezoids[-1].lower_right)
j += 1
return trapezoids return trapezoids
def replace(self, sink, node): def replace(self, sink, node):
if not sink.parent_list: if sink.parent_list:
self.head = node node.replace(sink)
else: else:
node.replace(sink) self.head = node
def case1(self, sink, edge, tlist): def case1(self, sink, edge, tlist):
yNode = YNode(edge, isink(tlist[1]), isink(tlist[2])) yNode = YNode(edge, isink(tlist[1]), isink(tlist[2]))
@ -554,13 +553,13 @@ class MonotoneMountain:
self.head = point self.head = point
self.size = 1 self.size = 1
elif self.size is 1: elif self.size is 1:
if point != self.head: if point.neq(self.head):
self.tail = point self.tail = point
self.tail.prev = self.head self.tail.prev = self.head
self.head.next = self.tail self.head.next = self.tail
self.size = 2 self.size = 2
else: else:
if point != self.tail: if point.neq(self.tail):
self.tail.next = point self.tail.next = point
point.prev = self.tail point.prev = self.tail
self.tail = point self.tail = point
@ -599,7 +598,7 @@ class MonotoneMountain:
self.convex_points.append(a) self.convex_points.append(a)
if self.valid(c): if self.valid(c):
self.convex_points.append(c) self.convex_points.append(c)
assert(self.size <= 3, "Triangulation bug, please report") #assert self.size <= 3, "Triangulation bug, please report"
def valid(self, p): def valid(self, p):
return p != self.head and p != self.tail and self.is_convex(p) return p != self.head and p != self.tail and self.is_convex(p)

View File

@ -237,6 +237,7 @@ class Poly2TriDemo extends BasicGame("Poly2Tri") {
// Correctly adjust for pan and zoom // Correctly adjust for pan and zoom
if(mouseButton == 1) { if(mouseButton == 1) {
val point = mousePos/scaleFactor + Point(deltaX, deltaY) val point = mousePos/scaleFactor + Point(deltaX, deltaY)
println(point)
slCDT.addPoint(point) slCDT.addPoint(point)
slCDT.triangulate slCDT.triangulate
} }
@ -333,8 +334,8 @@ class Poly2TriDemo extends BasicGame("Poly2Tri") {
val clearPoint = Point(365, 427) val clearPoint = Point(365, 427)
loadModel(dude, -1f, Point(100f, -200f), 10, clearPoint) loadModel(dude, -1f, Point(100f, -200f), 10, clearPoint)
case "../data/basic.dat" => case "../data/basic.dat" =>
val clearPoint = Point(365, 427) val clearPoint = Point(500, 450)
loadModel(basic, -1f, Point(400f, 300f), 10, clearPoint) loadModel(basic, -5f, Point(500, 450), 10, clearPoint)
case _ => case _ =>
} }
currentModel = model currentModel = model