Changeset 959

Timestamp:

Feb 25, 2005, 2:36:07 PM (20 years ago)

Author:

duncan

Message:

update from vanessa

Location:

inundation/ga/storm_surge/alpha_shape

Files:

• alpha_shape.py (modified) (10 diffs)
• test_alpha_shape.py (modified) (2 diffs)

inundation/ga/storm_surge/alpha_shape/alpha_shape.py

@@ -3 +3 @@
 
 import exceptions
-from Numeric import array, Float, divide_safe, sqrt
+from Numeric import array, Float, divide_safe, sqrt, product
 from load_mesh.loadASCII import load_points_file, export_boundary_file
 import random
@@ -12 +12 @@
 INF = pow(10,9)
 
-def alpha_shape_via_files(point_file, boundary_file, alpha= None):    
+def alpha_shape_via_files(point_file, boundary_file, alpha= None):
     
     point_dict = load_points_file(point_file)
@@ -70 +70 @@
         return self.boundary
 
-    def set_boundary_type(self,flag=0):
+    def set_boundary_type(self,flag=0,small=0.2):
         """
         Use the flag to set constraints on the boundary
@@ -83 +83 @@
         if flag>=1 :
             #remove holes
-            bdry = self._remove_holes(bdry)
+            bdry = self._remove_holes(bdry,small)
         if flag>=2:
             #remove sharp indents
@@ -420 +420 @@
         """
         Given an alpha value,
-        return the edges on the boundary of the alpha-shape
+        return the indices of edges on the boundary of the alpha-shape
         """
         def reg_edge(k):
@@ -437 +437 @@
         return filter(exp_vert, range(len(self.vertexinterval)))        
 
-        
-    def _remove_holes(self,bdry):
+    def _vertices_from_edges(self,elist):
+        """
+        Returns the list of unique vertex labels from edges in elist
+        """
+
+        v1 = [elist[k][0] for k in range(len(elist))]
+        v2 = [elist[k][1] for k in range(len(elist))]
+        v = v1+v2
+        v.sort()
+        vertices = [v[k] for k in range(len(v)) if v[k]!=v[k-1] ]
+        return vertices
+
+        
+    def _remove_holes(self,bdry,small):
         """
         Given the edges in bdry, find the largest exterior components.
@@ -452 +464 @@
 
         # get a list of unique vertex labels:
-        v1 = [bdry[k][0] for k in range(len(bdry))]
-        v2 = [bdry[k][1] for k in range(len(bdry))]
-        v = v1+v2
-        v.sort()
-        verts = [v[k] for k in range(len(v)) if v[k]!=v[k-1] ]
+        verts = self._vertices_from_edges(bdry)
         #print "verts ", verts, "\n"
 
@@ -501 +509 @@
 
         # discard the edges in the little components
-        # (i.e. those components with less than twenty percent of bdry points)
-        cutoff = round(0.2*len(verts))
+        # (i.e. those components with less than 'small' fraction of bdry points)
+        cutoff = round(small*len(verts))
         littleind = [k for k in range(len(vptr)) if (vptr[k]<0 and vptr[k]>-cutoff)] 
         if littleind:
@@ -519 +527 @@
         Given edges in bdry, test for acute-angle indents and remove them.
         """
+        # first find triangles not in alpha-shape
+        def tri_rad_gta(k):
+            return self.triradius[k]>self.alpha
+
+        extrind = filter(tri_rad_gta, range(len(self.triradius)))
+        #print "exterior triangles: ", extrind
+
+        # find triangles that share two edges with bdry
+        # i.e. triangles that have all verts in bdry.
+        verts = self._vertices_from_edges(bdry)
+
+        bdrytri = []
+        for ind in extrind:
+            bvct = 0
+            for j in [0,1,2]:
+                if self.deltri[ind][j] in verts:
+                    bvct +=1
+            if bvct==3:
+                bdrytri.append(ind)
+
+        #print "boundary triangles ", bdrytri
+
+        # test the bdrytri triangles for acute angles 
+        probtri = []
+        for tind in bdrytri:
+            tri = self.deltri[tind]
+            dx = array([self.points[tri[(i+1)%3],0] - self.points[tri[(i+2)%3],0] for i in [0,1,2]])
+            dy = array([self.points[tri[(i+1)%3],1] - self.points[tri[(i+2)%3],1] for i in [0,1,2]])
+            anglesign = array([(-dx[(i+1)%3]*dx[(i+2)%3]-dy[(i+1)%3]*dy[(i+2)%3]) for i in [0,1,2]])
+            if product(anglesign) > 0:
+                probtri.append(tind)
+
+        #print "problem triangles ", probtri 
+
+        # adjust the bdry edges of the probtri triangles
+        for pind in probtri:
+            tri = self.deltri[pind]
+            #print "triangle verts: ", tri
+            for i in [0,1,2]:
+                edga = (tri[(i+1)%3], tri[(i+2)%3])
+                edgb = (tri[(i+2)%3], tri[(i+1)%3])
+                if edga in bdry:
+                    bdry.remove(edga)
+                elif edgb in bdry:
+                    bdry.remove(edgb)
+                else:
+                    bdry.append(edga)       
+        
+        newbdry = bdry
         return newbdry
 
@@ -525 +582 @@
         Given edges in bdry, test for vertices with more than 4 incident edges.
         Expand by adding back in associated triangles... 
-        
-        """
+        """
+
+        #verts = self._vertices_from_edges(bdry)
+
+        v1 = [bdry[k][0] for k in range(len(bdry))]
+        v2 = [bdry[k][1] for k in range(len(bdry))]
+        v = v1+v2
+        v.sort()
+        probv = [v[k] for k in range(len(v)) if (v[k]!=v[k-1] and v.count(k)>2) ]
+
+        # first find triangles not in alpha-shape
+        def tri_rad_gta(k):
+            return self.triradius[k]>self.alpha
+
+        extrind = filter(tri_rad_gta, range(len(self.triradius)))
+        #print "exterior triangles: ", extrind
+
+        # find exterior triangles that have a vertex in probv
+
+        bdrytri = []
+        for ind in extrind:
+            bvct = 0
+            for j in [0,1,2]:
+                if self.deltri[ind][j] in probv:
+                    bvct +=1
+            if bvct>0:
+                bdrytri.append(ind)
+
+        #print "boundary triangles ", bdrytri
+
+
+                
+        newbdry = bdry
         return newbdry
    

inundation/ga/storm_surge/alpha_shape/test_alpha_shape.py

@@ -97 +97 @@
 
 
-    def test_boundary(self):
+    def test_boundary_1(self):
         a = [0.0, 0.0]
         b = [2.0, 0.0]
@@ -111 +111 @@
         #print "result",result
         answer = [(0, 3), (3, 6), (0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 6)]
+        assert allclose(answer, result) 
+        pass
+
+
+    def test_boundary_2(self):
+        a = [0.0, 0.0]
+        b = [2.0, 0.0]
+        c = [4.0, 0.0]
+        cd = [3.0,2.0]
+        d = [4.0, 4.0]
+        e = [2.0, 4.0]
+        f = [0.0, 4.0]
+
+        alpha = Alpha_Shape([a,b,c,cd,d,e,f])
+        alpha.set_boundary_type(2)
+        result = alpha.get_boundary()
+        #print "result",result
+        result [(0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 6), (6, 0)]
         assert allclose(answer, result) 
         pass