Changeset 956

Timestamp:

Feb 25, 2005, 11:59:06 AM (20 years ago)

Author:

duncan

Message:

update from vanessa

Location:

inundation/ga/storm_surge/alpha_shape

Files:

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

inundation/ga/storm_surge/alpha_shape/alpha_shape.py

@@ -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)
@@ -69 +69 @@
         """
         return self.boundary
+
+    def set_boundary_type(self,flag=0):
+        """
+        Use the flag to set constraints on the boundary
+        0   Return raw boundary i.e. the regular edges
+        1   filter to remove small holes
+        2   (includes 1) plus remove sharp indents in boundary
+        3   (includes 2) plus test for pinch-off i.e. boundary vertex with more than two edges.
+        """
+        reg_edge = self.get_regular_edges(self.alpha)
+        bdry = [self.edge[k] for k in reg_edge]
+
+        if flag>=1 :
+            #remove holes
+            bdry = self._remove_holes(bdry)
+        if flag>=2:
+            #remove sharp indents
+            bdry = self._smooth_indents(bdry)
+        if flag>=3:
+            #deal with pinch-off
+            bdry = self._expand_pinch(bdry)
+        
+        self.boundary = bdry
+        
 
     def get_delaunay(self):
@@ -220 +244 @@
             print "There are serious problems with the delaunay triangles"
             raise
-            # really need to take care of cases when denom = 0.
-            # eg: something like:
-            # print "Handling degenerate triangles."
-            # ... add some random displacments to trouble points? 
+ 
             
         self.triradius = 0.5*sqrt(dx*dx + dy*dy)
@@ -406 +427 @@
         return filter(reg_edge, range(len(self.edgeinterval)))
 
-        
+    def get_exposed_vertices(self,alpha):
+        """
+        Given an alpha value,
+        return the vertices on the boundary of the alpha-shape
+        """
+        def exp_vert(k):
+            return self.vertexinterval[k][0]<=alpha and self.vertexinterval[k][1]>alpha
+
+        return filter(exp_vert, range(len(self.vertexinterval)))        
+
+        
+    def _remove_holes(self,bdry):
+        """
+        Given the edges in bdry, find the largest exterior components.
+        """
+        def findroot(i):
+            if vptr[i] < 0:
+                return i
+            k = findroot(vptr[i])
+            vptr[i] = k
+            return k
+
+
+
+        # 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] ]
+        #print "verts ", verts, "\n"
+
+        #initialise vptr and eptr to negative number outside range 
+        EMPTY = -max(verts) - len(bdry)
+        vptr = [EMPTY for k in range(len(verts))]
+        eptr = [EMPTY for k in range(len(bdry))]
+        #print "vptr init: ", vptr, "\n"
+
+        #add edges and maintain union tree
+        for i in range(len(bdry)):
+            #print "edge ",i,"\t",bdry[i] 
+            vl = verts.index(bdry[i][0])
+            vr = verts.index(bdry[i][1])
+            #rvl = vl, rvr = vr
+            rvl = findroot(vl)
+            rvr = findroot(vr)
+            #print "roots:  ",rvl, rvr
+            if not(rvl==rvr):
+                if (vptr[vl]==EMPTY):
+                    if (vptr[vr]==EMPTY):
+                        vptr[vl] = -2
+                        vptr[vr] = vl
+                        #eprt[i] = vl
+                    else:
+                        vptr[vl] = rvr
+                        vptr[rvr] -= 1
+                else:
+                    if (vptr[vr]==EMPTY):
+                        vptr[vr] = rvl
+                        vptr[rvl] -= 1
+                    else:
+                        if vptr[rvl] > vptr[rvr]:
+                            vptr[rvl] = rvr
+                            vptr[rvr] += vptr[rvl]
+                        else:
+                            vptr[rvr] = rvl
+                            vptr[rvl] += vptr[rvr]
+                #print "vptr: ", vptr, "\n"
+
+        # end edge loop
+
+        #print "vertex component tree: ", vptr, "\n"
+
+        # 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))
+        littleind = [k for k in range(len(vptr)) if (vptr[k]<0 and vptr[k]>-cutoff)] 
+        if littleind:
+            littlecomp = [k for k in range(len(vptr)) if findroot(k) in littlecomp]
+            vdiscard = [verts[k] for k in littlecomp] 
+            newbdry = [e for e in bdry if not((e[0] in vdiscard) or (e[1] in vdiscard))]
+        else:
+            newbdry = bdry
+
+        return newbdry
+
+
+
+    def _smooth_indents(self,bdry):
+        """
+        Given edges in bdry, test for acute-angle indents and remove them.
+        """
+        return newbdry
+
+    def _expand_pinch(self,bdry):
+        """
+        Given edges in bdry, test for vertices with more than 4 incident edges.
+        Expand by adding back in associated triangles... 
+        
+        """
+        return newbdry
    
 #-------------------------------------------------------------

inundation/ga/storm_surge/alpha_shape/test_alpha_shape.py

@@ -97 +97 @@
 
 
-    
-    def test_alpha_optional_alpha(self):
-        #print "test_alpha" 
+    def test_boundary(self):
         a = [0.0, 0.0]
-        b = [1.0, 0.0]
-        c = [2.0, 0.0]
-        d = [2.0, 2.0]
-        e = [1.0, 2.0]
-        f = [0.0, 2.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]
 
-        shape = Alpha_Shape([a,b,c,d,e,f])
-        result = shape.get_boundary()
+        alpha = Alpha_Shape([a,b,c,cd,d,e,f])
+        alpha.set_boundary_type(1)
+        result = alpha.get_boundary()
         #print "result",result
-        answer = [(5, 0), (0, 1), (4, 5), (2, 3), (3, 4), (1, 2)]
-        assert allclose(answer, result)
+        answer = [(0, 3), (3, 6), (0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 6)]
+        assert allclose(answer, result) 
+        pass
+
         
-        optimum_alpha = shape.get_optimum_alpha()
-        alpha = shape.get_alpha()
-        
-        self.failUnless( optimum_alpha == alpha,
-                        'alpha is wrong!')
-
-        #New instance, setting optional alpha
-        shape = Alpha_Shape([a,b,c,d,e,f], alpha = optimum_alpha)
-        result = shape.get_boundary()
-        #print "result",result
-        answer = [(5, 0), (0, 1), (4, 5), (2, 3), (3, 4), (1, 2)]
-        assert allclose(answer, result)
-
-        shape.set_alpha(optimum_alpha)
-        result = shape.get_boundary()
-        #print "result",result
-        answer = [(5, 0), (0, 1), (4, 5), (2, 3), (3, 4), (1, 2)]
-        assert allclose(answer, result)
-        
-        
-            
     def test_not_enough_points(self):
         #print "test_not_enought_points"
@@ -156 +137 @@
         #(h,fileName) = tempfile.mkstemp(".xya") 
         file = open(fileName,"w")
-        file.write("title row \n\
+        file.write("\n\
 0.0, 0.0\n\
 1.0, 0.0\n\