Other fixes to the Newell's algorithm implementation
[vrm.git] / vrm.py
diff --git a/vrm.py b/vrm.py
index 9095307..7cf79bc 100755 (executable)
--- a/vrm.py
+++ b/vrm.py
@@ -92,7 +92,7 @@ class config:
     polygons = dict()
     polygons['SHOW'] = True
     polygons['SHADING'] = 'FLAT'
-    polygons['HSR'] = 'PAINTER' # 'PAINTER' or 'NEWELL'
+    #polygons['HSR'] = 'PAINTER' # 'PAINTER' or 'NEWELL'
     polygons['HSR'] = 'NEWELL'
     # Hidden to the user for now
     polygons['EXPANSION_TRICK'] = True
@@ -113,12 +113,22 @@ class config:
 
 
 
-# Debug utility function
-print_debug = True
+# Utility functions
+print_debug = False
 def debug(msg):
     if print_debug:
         sys.stderr.write(msg)
 
+EPS = 10e-5
+
+def sign(x):
+    if x < -EPS:
+        return -1
+    elif x > EPS:
+        return 1
+    else:
+        return 0
+
 
 # ---------------------------------------------------------------------
 #
@@ -800,7 +810,8 @@ class SVGVectorWriter(VectorWriter):
             opacity_string = ""
             if color[3] != 255:
                 opacity = float(color[3])/255.0
-                opacity_string = " fill-opacity: %g; stroke-opacity: %g; opacity: 1;" % (opacity, opacity)
+                #opacity_string = " fill-opacity: %g; stroke-opacity: %g; opacity: 1;" % (opacity, opacity)
+                opacity_string = "opacity: %g;" % (opacity)
 
             self.file.write("\tstyle=\"fill:" + str_col + ";")
             self.file.write(opacity_string)
@@ -810,8 +821,9 @@ class SVGVectorWriter(VectorWriter):
             # see http://www.antigrain.com/svg/index.html for more info
             stroke_width = 1.0
 
-            if config.polygons['EXPANSION_TRICK']:
-                str_col = "#000000" # For debug
+            # EXPANSION TRICK is not that useful where there is transparency
+            if config.polygons['EXPANSION_TRICK'] and color[3] == 255:
+                # str_col = "#000000" # For debug
                 self.file.write(" stroke:%s;\n" % str_col)
                 self.file.write(" stroke-width:" + str(stroke_width) + ";\n")
                 self.file.write(" stroke-linecap:round;stroke-linejoin:round")
@@ -1026,8 +1038,18 @@ class Renderer:
 
             self._doBackFaceCulling(mesh)
 
+            # When doing HSR with NEWELL we may want to flip all normals
+            # toward the viewer
+            if config.polygons['HSR'] == "NEWELL":
+                for f in mesh.faces:
+                    f.sel = 1-f.sel
+                mesh.flipNormals()
+                for f in mesh.faces:
+                    f.sel = 1
+
             self._doLighting(mesh)
 
+
             # Do "projection" now so we perform further processing
             # in Normalized View Coordinates
             self._doProjection(mesh, self.proj)
@@ -1144,8 +1166,8 @@ class Renderer:
     def _doConvertGeometricObjsToMesh(self, scene):
         """Convert all "geometric" objects to mesh ones.
         """
-        #geometricObjTypes = ['Mesh', 'Surf', 'Curve', 'Text']
-        geometricObjTypes = ['Mesh', 'Surf', 'Curve']
+        geometricObjTypes = ['Mesh', 'Surf', 'Curve', 'Text']
+        #geometricObjTypes = ['Mesh', 'Surf', 'Curve']
 
         Objects = scene.getChildren()
         objList = [ o for o in Objects if o.getType() in geometricObjTypes ]
@@ -1438,7 +1460,7 @@ class Renderer:
         
 
         by_furthest_z = (lambda f1, f2: progress.update() and
-                cmp(max([v.co[2] for v in f1]), max([v.co[2] for v in f2]))
+                cmp(max([v.co[2] for v in f1]), max([v.co[2] for v in f2])+EPS)
                 )
 
         # FIXME: using NMesh to sort faces. We should avoid that!
@@ -1462,53 +1484,80 @@ class Renderer:
         """Newell's depth sorting.
 
         """
+        global EPS
+
         by_furthest_z = (lambda f1, f2:
-                cmp(max([v.co[2] for v in f1]), max([v.co[2] for v in f2]))
+                cmp(max([v.co[2] for v in f1]), max([v.co[2] for v in f2])+EPS)
                 )
 
+        mesh.quadToTriangle()
 
-        def isOnSegment(v1, v2, p):
+        from split import Distance, isOnSegment
 
-            # when p is at extreme points
-            if p == v1 or p == v2:
-                return False
+        def projectionsOverlap(P, Q):
 
+            for i in range(0, len(P.v)):
 
-            EPS = 10e-7
+                v1 = Vector(P.v[i-1])
+                v1[2] = 0
+                v2 = Vector(P.v[i])
+                v2[2] = 0
 
-            l1 = (v1-p).length
-            l2 = (v2-p).length
-            l = (v1-v2).length
+                EPS = 10e-5
 
-            print "l: ", l, " l1: ", l1, " l2: ", l2, "diff: %.9f" % (l - (l1+l2) )
+                for j in range(0, len(Q.v)):
 
-            if abs(l - (l1+l2)) < EPS:
-                return True
-            else:
-                return False
+                    v3 = Vector(Q.v[j-1])
+                    v3[2] = 0
+                    v4 = Vector(Q.v[j])
+                    v4[2] = 0
+
+                    #print "\n\nTEST if we have coincidence!"
+                    #print v1, v2
+                    #print v3, v4
+                    #print "distances:"
+                    d1 = (v1-v3).length
+                    d2 = (v1-v4).length
+                    d3 = (v2-v3).length
+                    d4 = (v2-v4).length
+                    #print d1, d2, d3, d4
+                    #print "-----------------------\n"
 
+                    if d1 < EPS or d2 < EPS or d3 < EPS or d4 < EPS:
+                        continue
+                    
+                    # TODO: Replace with LineIntersect2D in newer API
+                    ret = LineIntersect(v1, v2, v3, v4)
 
+                    # if line v1-v2 and v3-v4 intersect both return
+                    # values are the same.
+                    if ret and ret[0] == ret[1]  and isOnSegment(v1, v2, ret[0], True) and isOnSegment(v3, v4, ret[1], True):
 
-        def Distance(point, face):
-            """ Calculate the distance between a point and a face.
+                        #l1 = (ret[0] - v1).length
+                        #l2 = (ret[0] - v2).length
 
-            An alternative but more expensive method can be:
+                        #l3 = (ret[1] - v3).length
+                        #l4 = (ret[1] - v4).length
 
-                ip = Intersect(Vector(face[0]), Vector(face[1]), Vector(face[2]),
-                        Vector(face.no), Vector(point), 0)
+                        #print "New DISTACES againt the intersection point:"
+                        #print l1, l2, l3, l4
+                        #print "-----------------------\n"
 
-                d = Vector(ip - point).length
-            """
+                        #if  l1 < EPS or l2 < EPS or l3 < EPS or l4 < EPS:
+                        #    continue
 
-            plNormal = Vector(face.no)
-            plVert0 = Vector(face[0])
+                        debug("Projections OVERLAP!!\n")
+                        debug("line1:"+
+                                " M "+ str(v1[0])+','+str(v1[1]) + ' L ' + str(v2[0])+','+str(v2[1]) + '\n' +
+                                " M "+ str(v3[0])+','+str(v3[1]) + ' L ' + str(v4[0])+','+str(v4[1]) + '\n' +
+                                "\n")
+                        debug("return: "+ str(ret)+"\n")
+                        return True
 
-            #d = abs( (point * plNormal ) - (plVert0 * plNormal) )
-            d = (point * plNormal ) - (plVert0 * plNormal)
-            debug("d: %.10f - sel: %d, %s\n" % (d, face.sel, str(point)) )
+            return False
 
-            return d
 
+        from facesplit import facesplit
 
         # FIXME: using NMesh to sort faces. We should avoid that!
         nmesh = NMesh.GetRaw(mesh.name)
@@ -1527,52 +1576,64 @@ class Renderer:
         facelist = nmesh.faces[:]
         maplist = []
 
-        #EPS = 10e-7
-        EPS = 0
+        EPS = 10e-5
 
         global progress
+
+        # The steps are _at_least_ equal to len(facelist), we do not count the
+        # feces coming out from splitting!!
         progress.setActivity("HSR: Newell", len(facelist))
-        progress.setQuiet(True)
+        #progress.setQuiet(True)
+
+        
+        #split_done = 0
+        #marked_face = 0
 
-        #while len(facelist)-1:
         while len(facelist):
+            debug("\n----------------------\n")
+            debug("len(facelits): %d\n" % len(facelist))
             P = facelist[0]
 
-            pSign = 1
-            if P.sel == 0:
-                pSign = -1
+            pSign = sign(P.normal[2])
+
+            # We can discard faces parallel to the view vector
+            if pSign == 0:
+                facelist.remove(P)
+                continue
+
+            split_done = 0
+            face_marked = 0
 
-            #while False:
             for Q in facelist[1:]:
 
                 debug("P.smooth: " + str(P.smooth) + "\n")
                 debug("Q.smooth: " + str(Q.smooth) + "\n")
                 debug("\n")
 
-                qSign = 1
-                if Q.sel == 0:
-                    qSign = -1
+                qSign = sign(Q.normal[2])
  
                 # We need to test only those Qs whose furthest vertex
                 # is closer to the observer than the closest vertex of P.
 
                 zP = [v.co[2] for v in P.v]
                 zQ = [v.co[2] for v in Q.v]
-                ZOverlap = min(zP) < max(zQ)
+                notZOverlap = min(zP) > max(zQ) + EPS
 
-                if not ZOverlap:
+                if notZOverlap:
                     debug("\nTest 0\n")
                     debug("NOT Z OVERLAP!\n")
-                    if not Q.smooth:
-                        # We can safely print P
+                    if Q.smooth == 0:
+                        # If Q is not marked then we can safely print P
                         break
                     else:
+                        debug("met a marked face\n")
                         continue
                 
                 # Test 1: X extent overlapping
                 xP = [v.co[0] for v in P.v]
                 xQ = [v.co[0] for v in Q.v]
-                notXOverlap = (max(xP) < min(xQ)) or (max(xQ) < min(xP))
+                #notXOverlap = (max(xP) <= min(xQ)) or (max(xQ) <= min(xP))
+                notXOverlap = (min(xQ) >= max(xP)-EPS) or (min(xP) >= max(xQ)-EPS)
 
                 if notXOverlap:
                     debug("\nTest 1\n")
@@ -1582,7 +1643,8 @@ class Renderer:
                 # Test 2: Y extent Overlapping
                 yP = [v.co[1] for v in P.v]
                 yQ = [v.co[1] for v in Q.v]
-                notYOverlap = (max(yP) < min(yQ)) or (max(yQ) < min(yP))
+                #notYOverlap = (max(yP) <= min(yQ)) or (max(yQ) <= min(yP))
+                notYOverlap = (min(yQ) >= max(yP)-EPS) or (min(yP) >= max(yQ)-EPS)
 
                 if notYOverlap:
                     debug("\nTest 2\n")
@@ -1593,9 +1655,8 @@ class Renderer:
                 # Test 3: P vertices are all behind the plane of Q
                 n = 0
                 for Pi in P:
-                    print P.col[0]
                     d = qSign * Distance(Vector(Pi), Q)
-                    if d > EPS:
+                    if d <= EPS:
                         n += 1
                 pVerticesBehindPlaneQ = (n == len(P))
 
@@ -1608,9 +1669,8 @@ class Renderer:
                 # Test 4: Q vertices in front of the plane of P
                 n = 0
                 for Qi in Q:
-                    print Q.col[0]
                     d = pSign * Distance(Vector(Qi), P)
-                    if d <= EPS:
+                    if d >= -EPS:
                         n += 1
                 qVerticesInFrontPlaneP = (n == len(Q))
 
@@ -1621,62 +1681,36 @@ class Renderer:
 
                 # Test 5: Line Intersections... TODO
                 # Check if polygons effectively overlap each other, not only
-                # boundig boxes as dome before.
+                # boundig boxes as done before.
                 # Since we We are working in normalized projection coordinates
                 # we kust check if polygons intersect.
 
-                def projectionsOverlap(P, Q):
-
-                    for i in range(0, len(P.v)):
-
-                        v1 = Vector(P.v[i-1])
-                        v1[2] = 0
-                        v2 = Vector(P.v[i])
-                        v2[2] = 0
-
-                        for j in range(0, len(Q.v)):
-                            v3 = Vector(Q.v[j-1])
-                            v3[2] = 0
-                            v4 = Vector(Q.v[j])
-                            v4[2] = 0
-                            
-                            ret = LineIntersect(v1, v2, v3, v4)
-                            # if line v1-v2 and v3-v4 intersect both return
-                            # values are the same.
-                            if ret and ret[0] == ret[1]  and isOnSegment(v1, v2,
-                                    ret[0]) and isOnSegment(v3, v4, ret[1]):
-                                debug("Projections OVERLAP!!\n")
-                                debug("line1:"+
-                                        " M "+ str(v1[0])+','+str(v1[1]) + ' L ' + str(v2[0])+','+str(v2[1]) + '\n' +
-                                        " M "+ str(v3[0])+','+str(v3[1]) + ' L ' + str(v4[0])+','+str(v4[1]) + '\n' +
-                                        "\n")
-                                debug("return: "+ str(ret)+"\n")
-                                return True
-
-                    return False
-
                 if not projectionsOverlap(P, Q):
                     debug("\nTest 5\n")
                     debug("Projections do not overlap!\n")
                     continue
 
 
-                # We do not know if P obscures Q.
+                # We still do not know if P obscures Q.
+
+                # But if Q is marked we do a split trying to resolve a
+                # difficulty (maybe a visibility cycle).
                 if Q.smooth == 1:
-                    # Split P or Q, TODO
-                    debug("Cycle detected!\n")
+                    # Split P or Q
+                    debug("Possibly a cycle detected!\n")
                     debug("Split here!!\n")
-                    continue
 
+                    facelist = facesplit(P, Q, facelist, nmesh)
+                    split_done = 1
+                    break 
 
                 # The question now is: Does Q obscure P?
 
                 # Test 3bis: Q vertices are all behind the plane of P
                 n = 0
                 for Qi in Q:
-                    print Q.col[0]
                     d = pSign * Distance(Vector(Qi), P)
-                    if d > EPS:
+                    if d <= EPS:
                         n += 1
                 qVerticesBehindPlaneP = (n == len(Q))
 
@@ -1688,9 +1722,8 @@ class Renderer:
                 # Test 4bis: P vertices in front of the plane of Q
                 n = 0
                 for Pi in P:
-                    print P.col[0]
                     d = qSign * Distance(Vector(Pi), Q)
-                    if d <= EPS:
+                    if d >= -EPS:
                         n += 1
                 pVerticesInFrontPlaneQ = (n == len(P))
 
@@ -1698,58 +1731,41 @@ class Renderer:
                     debug("\nTest 4bis\n")
                     debug("P IN FRONT OF Q!\n")
 
-
-                import intersection
-
+                
+                # We don't even know if Q does obscure P, so they should
+                # intersect each other, split one of them in two parts.
                 if not qVerticesBehindPlaneP and not pVerticesInFrontPlaneQ:
                     debug("\nSimple Intersection?\n")
-                    # Split P or Q, TODO
-                    print "Test 3bis or 4bis failed"
-                    print "Split here!!2\n"
-
-                    """newfaces = intersection.splitOn(P, Q, 0)
-                    print newfaces
-                    facelist.remove(Q)
-                    for nf in newfaces:
-                        if nf:
-                            nf.col = Q.col
-                            facelist.append(nf)
-                    """
-
-                    break
-
-                # We do not know
-                if Q.smooth:
-                    # split P or Q
-                    print "Split here!!\n"
-                    """
-                    newfaces = intersection.splitOn(P, Q, 0)
-                    facelist.remove(Q)
-                    for nf in newfaces:
-                        if nf:
-                            nf.col = Q.col
-                            facelist.append(nf)
-
-                    """ 
-                    break
+                    debug("Test 3bis or 4bis failed\n")
+                    debug("Split here!!2\n")
 
-                Q.smooth = 1
+                    facelist = facesplit(P, Q, facelist, nmesh)
+                    split_done = 1
+                    break 
+                    
                 facelist.remove(Q)
                 facelist.insert(0, Q)
+                Q.smooth = 1
+                face_marked = 1
+                debug("Q marked!\n")
+                break
            
             # Write P!                     
-            P = facelist[0]
-            facelist.remove(P)
-            maplist.append(P)
+            if split_done == 0 and face_marked == 0:
+                facelist.remove(P)
+                maplist.append(P)
 
-            progress .update()
+                progress.update()
+
+            # end of while len(facelist)
+         
 
-        
         nmesh.faces = maplist
 
         for f in nmesh.faces:
             f.sel = 1
         nmesh.update()
+        #print nmesh.faces
 
     def _doHiddenSurfaceRemoval(self, mesh):
         """Do HSR for the given mesh.