Changeset 1228 for inundation/ga/storm_surge/pmesh/mesh.py

Timestamp:

Apr 15, 2005, 6:38:14 PM (20 years ago)

Author:

prow

Message:

added a new fuction - smooth polygons.
Comments.

File:

• inundation/ga/storm_surge/pmesh/mesh.py (modified) (26 diffs)

inundation/ga/storm_surge/pmesh/mesh.py

@@ -327 +327 @@
     def rotate(self,offset):
         """
+        permute the order of the sides of the triangle
         offset must be 0,1 or 2
         """
@@ -358 +359 @@
         return abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
     
+    def calcP(self):
+        #calculate the perimeter
+        ax = self.vertices[0].x
+        ay = self.vertices[0].y
+        
+        bx = self.vertices[1].x
+        by = self.vertices[1].y
+        
+        cx = self.vertices[2].x
+        cy = self.vertices[2].y
+
+        a =  ((cx-bx)**2+(cy-by)**2)**0.5 
+        b =  ((ax-cx)**2+(ay-cy)**2)**0.5 
+        c =  ((bx-ax)**2+(by-ay)**2)**0.5 
+
+        return a+b+c
+            
     def setNeighbors(self,neighbor1 = None, neighbor2 = None, neighbor3 = None):
         """
@@ -481 +499 @@
         self.sets=[[]]
        
+        self.visualise_graph = True
+
         if userSegments is None:
             self.userSegments=[]
@@ -1441 +1461 @@
         #In [3]: kinds.default_float_kind.M
         #kinds.default_float_kind.MAX         kinds.default_float_kind.MIN
-    #kinds.default_float_kind.MAX_10_EXP  kinds.default_float_kind.MIN_10_EXP
+    #kinds.default_float_kind.MAX_10_EXP  kinds.default_fltesting oat_kind.MIN_10_EXP
         #kinds.default_float_kind.MAX_EXP     kinds.default_float_kind.MIN_EXP
 
@@ -1469 +1489 @@
         Structure: [xmin, ymin, xmax, ymax] 
         """
-        # FIXME dsg!!! large is a hack
+        # FIXME dsg!!! large is a hacktesting 
         #You want the kinds package, part of Numeric:
         #In [2]: import kinds
@@ -2167 +2187 @@
 
     def line2seg(self,line,tag=None):
-        point0 = self.ver2point(line[0])
-        point1 = self.ver2point(line[1])
+        point0 = self.point2ver(line[0])
+        point1 = self.point2ver(line[1])
         return Segment(point0,point1,tag=tag)
 
@@ -2177 +2197 @@
         return Vertex(point[0],point[1])
 
+    def smooth_polySet(self,min_radius=0.05):
+        #for all pairs of connecting segments:
+        #    propose a new segment that replaces the 2
+
+        #    If the difference between the new segment
+        #    and the old lines is small: replace the
+        #    old lines.
+
+        seg2line = self.seg2line
+        ver2point= self.ver2point
+        line2seg = self.line2seg
+        point2ver= self.point2ver
+
+        #create dictionaries of lines -> segments
+        userSegments = self.segs_to_dict(self.userSegments)
+        alphaSegments = self.segs_to_dict(self.alphaUserSegments)
+
+        #lump user and alpha segments
+        for key in alphaSegments.keys():
+            userSegments[key]=alphaSegments[key]
+
+        #point_keys = tuple -> vertex
+        #userVertices = vertex -> [line,line] - lines from that node
+        point_keys = {}
+        userVertices={}
+        for vertex in self.getUserVertices():
+            point = ver2point(vertex)
+            if not point_keys.has_key(point):
+                point_keys[point]=vertex
+                userVertices[vertex]=[]
+        for key in userSegments.keys():
+            line = key
+            point_0 = key[0]
+            point_1 = key[1]
+            userVertices[point_keys[point_0]].append(line)
+            userVertices[point_keys[point_1]].append(line)
+
+        for point in point_keys.keys():
+            try:
+            #removed keys can cause keyerrors
+                vertex = point_keys[point]
+                lines = userVertices[vertex]
+    
+                #if there are 2 lines on the node
+                if len(lines)==2:
+                    line_0 = lines[0]
+                    line_1 = lines[1]
+    
+                    #if the tags are the the same on the 2 lines
+                    if userSegments[line_0].tag == userSegments[line_1].tag:
+                        tag = userSegments[line_0].tag 
+    
+                        #point_a is one of the next nodes, point_b is the other
+                        if point==line_0[0]:
+                            point_a = line_0[1]
+                        if point==line_0[1]:
+                            point_a = line_0[0]
+                        if point==line_1[0]:
+                            point_b = line_1[1]
+                        if point==line_1[1]:
+                            point_b = line_1[0]
+    
+    
+                        #line_2 is proposed
+                        line_2 = (point_a,point_b)
+
+                        #calculate the area of the triangle between
+                        #the two existing segments and the proposed
+                        #new segment
+                        ax = point_a[0]
+                        ay = point_a[1]
+                        bx = point_b[0]
+                        by = point_b[1]
+                        cx = point[0]
+                        cy = point[1]
+                        area=abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
+
+                        #calculate the perimeter
+                        len_a =  ((cx-bx)**2+(cy-by)**2)**0.5 
+                        len_b =  ((ax-cx)**2+(ay-cy)**2)**0.5 
+                        len_c =  ((bx-ax)**2+(by-ay)**2)**0.5 
+                        perimeter = len_a+len_b+len_c
+
+                        #calculate the radius
+                        r = area/(2*perimeter)
+
+                        #if the radius is small: then replace the existing
+                        #segments with the new one
+                        if r < min_radius:
+                            if len_c < min_radius: append = False
+                            else: append = True
+                            #if the new seg is also time, don't add it
+                            if append:
+                                segment = self.line2seg(line_2,tag=tag)
+
+                            list_a=userVertices[point_keys[point_a]]
+                            list_b=userVertices[point_keys[point_b]]
+
+                            if line_0 in list_a:
+                                list_a.remove(line_0)
+                            else:
+                                list_a.remove(line_1)
+
+                            if line_0 in list_b:
+                                list_b.remove(line_0)
+                            else:
+                                list_b.remove(line_1)
+
+                            if append:
+                                list_a.append(line_2)
+                                list_b.append(line_2)
+                            else:
+                                if len(list_a)==0:
+                                    userVertices.pop(point_keys[point_a])
+                                    point_keys.pop(point_a)
+                                if len(list_b)==0:
+                                    userVertices.pop(point_keys[point_b])
+                                    point_keys.pop(point_b)
+
+                            userVertices.pop(point_keys[point])
+                            point_keys.pop(point)
+                            userSegments.pop(line_0)
+                            userSegments.pop(line_1)
+
+                            if append:
+                                userSegments[line_2]=segment
+            except:
+                pass
+
+        #self.userVerticies = userVertices.keys()
+        #self.userSegments = []
+        #for key in userSegments.keys():
+        #    self.userSegments.append(userSegments[key])
+        #self.alphaUserSegments = []
+
+        self.userVerticies = []
+        self.userSegments = []
+        self.alphaUserSegments = []
+
+        return userVertices,userSegments,alphaSegments
+
     def triangles_to_polySet(self,setName):
+        #self.smooth_polySet()
+
         seg2line = self.seg2line
         ver2point= self.ver2point
@@ -2191 +2354 @@
  
 
+        #create a dict of points to vertexes (tuple -> object)
+        #also create a set of vertexes (object -> True)
         point_keys = {}
         userVertices={}
@@ -2199 +2364 @@
                 userVertices[vertex]=True
 
-        #######
-        for vert in userVertices.keys():
-            assert point_keys[(vert.x,vert.y)]==vert
-        assert len(point_keys.keys())==len(userVertices.keys())
-        #######
-
-
+        #create a dict of lines to segments (tuple -> object)
         userSegments = self.segs_to_dict(self.userSegments)
+        #append the userlines in an affine linespace
         affine_lines = Affine_Linespace()
         for line in userSegments.keys():
             affine_lines.append(line)
-
         alphaSegments = self.segs_to_dict(self.alphaUserSegments)
         for line in alphaSegments.keys():
@@ -2228 +2387 @@
                     line = (point_a,point_b)
                     tag = None
+
+
+                    #this bit checks for matching lines
                     possible_lines = affine_lines[line] 
+                    possible_lines = unique(possible_lines)
                     found = 0                            
                     for user_line in possible_lines:
                         if self.point_on_line(midpoint,user_line):
+                            found+=1
+                            assert found<2
                             if userSegments.has_key(user_line):
                                 parent_segment = userSegments.pop(user_line)
-                            else:
+                            if alphaSegments.has_key(user_line):
                                 parent_segment = alphaSegments.pop(user_line)
                             tag = parent_segment.tag
@@ -2253 +2418 @@
                                 userSegments[newline]=segment
                                 affine_lines.append(newline)
-                            found=1
-                            break
+                            #break
+                    assert found<2
+
+
+
+                    #if no matching lines
                     if not found:
                         line_vertices = []
@@ -2269 +2438 @@
                         affine_lines.append(line)
         
-        #######
-        for vert in userVertices.keys():
-            assert point_keys[(vert.x,vert.y)]==vert
-        assert len(point_keys.keys())==len(userVertices.keys())
-        #######
         self.userVerticies = []
         self.userSegments = []
-        self.alphaSegments = []
-        #for key in userSegments.keys():
-        #    self.userSegments.append(userSegments[key])
-
-        #for key in alphaSegments.keys():
-        #    self.userSegments.append(alphaSegments[key])
-        #FIXME change alpha to user and sync with pmesh
-        #ie self.alphaSegments.append(alphaSegments[key])
-        #print point_keys.keys()
-        #print userVertices.keys()
-        #print userVertices.keys()
+        self.alphaUserSegments = []
+
         return userVertices,userSegments,alphaSegments
 
@@ -2310 +2465 @@
         answer = []
 
+        #if the new line does not share a 
+        #vertex with the old one
         if not (allclose(A_array,a_array)\
              or allclose(B_array,b_array)\
@@ -2344 +2501 @@
             return [sibling]
 
-
     def point_on_line(self,point,line):
+        #returns true within a tolerance of 3 degrees
         x=point[0]
         y=point[1]
@@ -2354 +2511 @@
         from Numeric import array, dot, allclose
         from math import sqrt
+        tol = 3. #DEGREES
+        tol = tol*3.1415/180
 
         a = array([x - x0, y - y0]) 
-        a_normal = array([a[1], -a[0]])
-
-        b = array([x1 - x0, y1 - y0])
-    
-        if allclose(dot(a_normal, b),0):
+        a_normal = array([a[1], -a[0]])      
+        len_a_normal = sqrt(sum(a_normal**2)) 
+
+        b = array([x1 - x0, y1 - y0])          
+        len_b = sqrt(sum(b**2)) 
+    
+        if abs(dot(a_normal, b)/(len_b*len_a_normal))< tol:
             #Point is somewhere on the infinite extension of the line
 
-            len_a = sqrt(sum(a**2))            
-            len_b = sqrt(sum(b**2))                             
+            len_a = sqrt(sum(a**2))                                     
             if dot(a, b) >= 0 and len_a <= len_b:
                return True
@@ -2370 +2530 @@
                return False
         else:
-          return False  
+          return False
 
     def line_length(self,line):
@@ -2379 +2539 @@
         return ((x1-x0)**2-(y1-y0)**2)**0.5     
 
-    def triangles_to_polySet2(self,setName):
-        from Numeric import array,allclose
-
-        #turn the triangles into a set
-        Triangles = self.sets[self.setID[setName]]
-        Triangles_dict = {}
-        for triangle in Triangles:
-            Triangles_dict[triangle]=None
-
-        userVertices = {}
-        userSegments = {}
-        point_keys = {}
-        affine_lines = {}
-        for vertex in self.getUserVertices():
-            point = (vertex.x,vertex.y)
-            point_keys[point]=vertex
-
-        line_keys = {}
-        for segment in self.getUserSegments():
-        #inlined would looks very ugly
-            vertex1 = segment.vertices[0]
-            vertex2 = segment.vertices[1]
-            point1 = (vertex1.x,vertex1.y)
-            point2 = (vertex2.x,vertex2.y)
-            #segment.vertices[0]=point_keys[point1]
-            #segment.vertices[1]=point_keys[point2]
-            #vertex1 = segment.vertices[0]
-            #vertex2 = segment.vertices[1]
-            #point1 = (vertex1.x,vertex1.y)
-            #point2 = (vertex2.x,vertex2.y)
-            line1 = (point1,point2)
-            line2 = (point2,point1)
-            AB = affine_line(point1,point2)
-            flat_AB_up = ceilAB[0]+AB[1]+AB[2]
-            if not (line_keys.has_key(line1) \
-                 or line_keys.has_key(line2)):
-                 line_keys[line1]=segment
-
-        for triangle in Triangles:
-            for i in (0,1,2):
-                #for every triangles neighbour:
-
-                if not Triangles_dict.has_key(triangle.neighbors[i]):
-                #if the neighbour is not in the set:
-                    a = triangle.vertices[i-1]
-                    b = triangle.vertices[i-2]
-                    if not point_keys.has_key((a.x,a.y)):
-                        #if point a does not already exist
-                        #then add it to the points.
-                        userVertices[a]=True
-                        point_keys[(a.x,a.y)]=a
-                    else: 
-                        a=point_keys[(a.x,a.y)]
-                        userVertices[a]=point_keys[(a.x,a.y)]
-                    assert userVertices.has_key(a)
-
-                    if not point_keys.has_key((b.x,b.y)):
-                        #if point b does not already exist
-                        #then add it to the points.
-                        userVertices[b]=True
-                        point_keys[(b.x,b.y)]=b
-                    else: 
-                        b=point_keys[(b.x,b.y)]
-                        userVertices[b]=point_keys[(b.x,b.y)]
-                    assert userVertices.has_key(b)
-
-                    if not (line_keys.has_key(((a.x,a.y),(b.x,b.y)))\
-                         or line_keys.has_key(((b.x,b.y),(a.x,a.y)))):
-                        #if the segment does not already exist then
-                        #add it to the segments
-                        assert ((a.x,a.y)!=(b.x,b.y))
-                        assert a!=b
-                        AB = affine_line(a,b)
-                        flat_AB = AB[0]+AB[1]+AB[2]
-                        if affine_lines.has_key(flat_AB):
-                            userSegments[Segment(a,b)]=None
-                        line_keys[((a.x,a.y),(b.x,b.y))]=None
-
-        userVertices,userSegments = self.weed(userVertices.keys(),userSegments.keys())
-        self.userVertices.extend(userVertices)
-        self.userSegments.extend(userSegments)
-        self.userVertices,self.userSegments = \
-           self.weed(self.userVertices,self.userSegments)
-
     def threshold(self,setName,min=None,max=None,attribute_name = 'elevation'):
         """
@@ -2472 +2548 @@
         if attribute_name in self.attributeTitles:
             i = self.attributeTitles.index(attribute_name)
+        else: i = -1#no attribute
         if not max == None:
             for t in triangles:
@@ -2482 +2559 @@
         self.sets[self.setID[setName]] = A
 
-    def general_threshold(self,setName,min=None,max=None,attribute_name = 'elevation',function = None):
-        """
-        threshold using  d
-        """
+    def general_threshold(self,setName,min=None,max=None\
+              ,attribute_name = 'elevation',function=None):
+        """
+        Thresholds the triangles 
+        """
+        from visual.graph import arange,ghistogram,color as colour
         triangles = self.sets[self.setID[setName]]
         A = []
+        data=[]
+        #data is for the graph
 
         if attribute_name in self.attributeTitles:
             i = self.attributeTitles.index(attribute_name)
+        else: i = -1
         if not max == None:
             for t in triangles:
-                if (min<function(t,i)<max):
+                value=function(t,i)
+                if (min<value<max):
                     A.append(t)
+                data.append(value)
         else:
             for t in triangles:
-                if (min<self.function(t,i)):
+                value=function(t,i)
+                if (min<value):
                     A.append(t)
+                data.append(value)
         self.sets[self.setID[setName]] = A
 
+        if self.visualise_graph:
+            if len(data)>0:
+                max=data[0]
+                min=data[0]
+                for value in data:
+                    if value > max:
+                        max = value
+                    if value < min:
+                        min = value
+
+                inc = (max-min)/100
+
+                histogram = ghistogram(bins=arange(min,max,inc),\
+                             color = colour.red)
+                histogram.plot(data=data)
+        
     def av_att(self,triangle,i):
-    #evaluates the average attribute of the vertices of a triangle.
-        V = triangle.getVertices()
-        a0 = (V[0].attributes[i])
-        a1 = (V[1].attributes[i])
-        a2 = (V[2].attributes[i])
-        return (a0+a1+a2)/3
+        if i==-1: return 1
+        else:
+            #evaluates the average attribute of the vertices of a triangle.
+            V = triangle.getVertices()
+            a0 = (V[0].attributes[i])
+            a1 = (V[1].attributes[i])
+            a2 = (V[2].attributes[i])
+            return (a0+a1+a2)/3
 
     def Courant_ratio(self,triangle,index):
         """
-        Not the true Courant ratio, just elevation on area
+        Uses the courant threshold
         """
         e = self.av_att(triangle,index)
         A = triangle.calcArea()
-        return e/A
+        P = triangle.calcP()
+        r = A/(2*P)
+        e = max(0.1,abs(e))
+        return r/e**0.5
 
     def Gradient(self,triangle,index):
@@ -3187 +3294 @@
         from math import frexp
         self.frexp = frexp
+        roundings = [self.round_up_rel,\
+        self.round_down_rel,self.round_flat_rel,\
+        self.round_down_abs,self.round_up_abs,\
+        self.round_flat_abs]#
+
+        self.roundings = roundings
 
     def __repr__(self):
@@ -3202 +3315 @@
         rounded_values=list(key)
 
-        roundings = [self.round_up_rel,\
-        self.round_down_rel,self.round_flat_rel,\
-        self.round_down_abs,self.round_up_abs,\
-        self.round_flat_abs]
-
-
-        roundings = \
-        [self.round_up_rel, self.round_down_rel, \
-         self.round_up_rel2,self.round_down_rel2,\
-         self.round_up_abs, self.round_down_abs]
+        roundings = list(self.roundings)
 
         #initialise rounded_values
@@ -3371 +3475 @@
         return self.__items__.keys()
 
-    
-
 
 class Mapped_Discretised_Tuple_Set(Discretised_Tuple_Set):
@@ -3412 +3514 @@
     Nearly.
     """
-    def __init__(self,p_rel = 3,t_rel=0.1):
+    def __init__(self,p_rel=4,t_rel=0.2):
         Mapped_Discretised_Tuple_Set.__init__\
             (self,self.affine_line,\
             p_rel=p_rel,t_rel=t_rel)
+
+        roundings = \
+        [self.round_down_rel,self.round_up_rel,self.round_flat_rel]
+        self.roundings = roundings
+        #roundings = \
+        #[self.round_down_abs,self.round_up_abs,self.round_flat_abs]
+        #self.roundings = roundings
 
     def affine_line(self,line):
@@ -3445 +3554 @@
             alpha = (-dif_y)/dif_x
             #a = alpha * b
-            b = 1.
+            b = -1.
             c = (x1*alpha+x2*alpha+y1+y2)/2.
             a = alpha*b
@@ -3454 +3563 @@
             c = (x1+x2+y1*beta+y2*beta)/2.
             b = beta*a
-
-        mag = (a**2+b**2+c**2)**(0.5)
-
-        if a == 0:
-            sign_a = 1.
-        else:
-            sign_a = a/((a**2)**0.5)
+        mag = abs(a)+abs(b)
         #This does not change the mathematical
-        #properties, but it makes comparism
-        if b == 0:
-            sign_b = 1.
-        else:
-            sign_b = b/((b**2)**0.5)
-        if c == 0:
-            sign_c = 1.
-        else:
-            sign_c = c/((c**2)**0.5)
-        a = a/mag*sign_a
-        b = b/mag*sign_b
-        c = c/mag*sign_c
+        #properties, but it makes comparism possible.
+
+        #note that the gradient is b/a, or (a/b)**-1.
+        #so 
+
+        #if a == 0:
+        #    sign_a = 1.
+        #else:
+        #    sign_a = a/((a**2)**0.5)
+        #if b == 0:
+        #    sign_b = 1.
+        #else:
+        #    sign_b = b/((b**2)**0.5)
+        #if c == 0:
+        #    sign_c = 1.
+        #else:
+        #    sign_c = c/((c**2)**0.5)
+        #a = a/mag*sign_a
+        #b = b/mag*sign_b
+        #c = c/mag*sign_c
+        a = a/mag
+        b = b/mag
+        c = c/mag
         return a,b,c
 
 
-    
+
 if __name__ == "__main__":
     #from mesh import *
@@ -3484 +3598 @@
     dict = importUngenerateFile("ungen_test.txt")
     m.addVertsSegs(dict)
-    print m
+    print m3