Changeset 987

Timestamp:

Mar 2, 2005, 1:22:13 PM (20 years ago)

Author:

duncan

Message:

vanessa's changes, my changes

Location:

inundation/ga/storm_surge/alpha_shape

Files:

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

inundation/ga/storm_surge/alpha_shape/alpha_shape.py

@@ -8 +8 @@
 
 class PointError(exceptions.Exception): pass
+class FlagError(exceptions.Exception): pass
 
 OUTPUT_FILE_TITLE = "# The alpha shape boundary defined by point index pairs of edges"
@@ -35 +36 @@
         self._set_points(points)
         self._alpha_shape_algorithm(alpha)
-            
 
     def _set_points(self, points):
@@ -70 +70 @@
         return self.boundary
 
-    def set_boundary_type(self,flag=0,small=0.2):
-        """
-        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 :
+    def set_boundary_type(self,raw_boundary=True,
+                          remove_holes=False,
+                          smooth_indents=False,
+                          expand_pinch=False,
+                          boundary_points_fraction=0.2):
+        """
+        Use the flags to set constraints on the boundary
+        raw_boundary   Return raw boundary i.e. the regular edges
+        remove_holes   filter to remove small holes
+        smooth_indents   remove sharp indents in boundary
+        expand_pinch   plus test for pinch-off
+                       i.e. boundary vertex with more than two edges.
+        """
+
+        if raw_boundary:
+            reg_edge = self.get_regular_edges(self.alpha)
+            self.boundary = [self.edge[k] for k in reg_edge]
+            self._init_boundary_triangles()
+        if remove_holes:
             #remove holes
-            bdry = self._remove_holes(bdry,small)
-        if flag>=2:
+            self.boundary = self._remove_holes(boundary_points_fraction)
+        if smooth_indents:
             #remove sharp indents
-            bdry = self._smooth_indents(bdry)
-        if flag>=3:
+            self.boundary = self._smooth_indents()
+        if expand_pinch:
             #deal with pinch-off
-            bdry = self._expand_pinch(bdry)
-        
-        self.boundary = bdry
+            self.boundary = self._expand_pinch()
         
 
@@ -117 +122 @@
         reg_edge = self.get_regular_edges(alpha)    
         self.boundary = [self.edge[k] for k in reg_edge]
+        self._init_boundary_triangles()
     
             
@@ -178 +184 @@
         self.boundary = [self.edge[k] for k in reg_edge]
         #print "alpha boundary edges ", self.boundary
+        self._init_boundary_triangles()
             
-        # this produces a baaad boundary
-        #boundary = []
-        #for point_index in range(len(self.points)-1):
-        #    boundary.append([point_index,point_index +1])
-        #boundary.append([len(self.points)-1,0])
-        #self.boundary = boundary
         return
 
@@ -250 +251 @@
         return
 
-    def _edge_intervals_array_based(self):
-        # THIS DOESN'T WORK YET 
-
-        # for each edge, find triples
-        # (length/2, min_adj_triradius, max_adj_triradius) if unattached
-        # (min_adj_triradius, min_adj_triradius, max_adj_triradius) if attached.
-
-        edges = []
-        edgenbrs = []
-        edgeinterval = []
-
-        x = self.points[:,0]
-        y = self.points[:,1]
-        ind0 = [self.deltri[j][0] for j in range(len(self.deltri))]
-        ind1 = [self.deltri[j][1] for j in range(len(self.deltri))]
-        ind2 = [self.deltri[j][2] for j in range(len(self.deltri))]
-
-        x0 = array([ x[j] for j in ind1 ])
-        y0 = array([ y[j] for j in ind1 ])
-        x1 = array([ x[j] for j in ind2 ])
-        y1 = array([ y[j] for j in ind2 ])
-        x2 = array([ x[j] for j in ind3 ])
-        y2 = array([ y[j] for j in ind3 ])
-
-        dx0 = x1-x2
-        dx1 = x2-x0
-        dx2 = x0-x1
-        dy0 = y1-y2
-        dy1 = y2-y0
-        dy2 = y0-y1
-
-        elen0 = sqrt(dx0*dx0 + dy0*dy0)
-        elen1 = sqrt(dx1*dx1 + dy1*dy1)
-        elen2 = sqrt(dx2*dx2 + dy2*dy2)
-
-        anglesgn0 = -dx1*dx2-dy1*dy2
-        anglesgn1 = -dx2*dx0-dy2*dy0
-        anglesgn2 = -dx0*dx1-dy0*dy1
-
-        # not sure how to convert the rest of the algorithm.... 
-        
-        for t in range(len(self.deltri)):
-            tri = self.deltri[t]
-            trinbr = self.deltrinbr[t]
-                          
-            for i in [0,1,2]:
-                j = (i+1)%3
-                k = (i+2)%3
-                if trinbr[i]==-1:
-                    edges.append((tri[j], tri[k]))
-                    edgenbrs.append((t, -1))
-                    edgeinterval.append([0.5*elen[i], self.triradius[t], INF])
-                elif (tri[j]<tri[k]):
-                    edges.append((tri[j], tri[k]))
-                    edgenbrs.append((t, trinbr[i]))
-                    edgeinterval.append([0.5*elen[i],\
-                                         min(self.triradius[t],self.triradius[trinbr[i]]),\
-                                             max(self.triradius[t],self.triradius[trinbr[i]]) ])
-                else:
-                    continue
-                #if angle[i]>pi/2:
-                if anglesign[i] < 0: 
-                    edgeinterval[-1][0] = edgeinterval[-1][1]
-                    
-        self.edge = edges
-        self.edgenbr = edgenbrs
-        self.edgeinterval = edgeinterval
-        #print "edges: ",edges, "\n"
-        #print "edge nbrs:", edgenbrs ,"\n"
-        #print "edge intervals: ",edgeinterval , "\n" 
 
 
@@ -449 +380 @@
         return vertices
 
-        
-    def _remove_holes(self,bdry,small):
+    def _init_boundary_triangles(self):
+        """
+        Creates the initial list of triangle indices
+        exterior to and touching the boundary of the alpha shape
+        """
+        def tri_rad_gta(k):
+            return self.triradius[k]>self.alpha
+
+        extrind = filter(tri_rad_gta, range(len(self.triradius)))
+
+        bv = self._vertices_from_edges(self.boundary)
+        
+        btri = []
+        for et in extrind:
+            v0 = self.deltri[et][0]
+            v1 = self.deltri[et][1]
+            v2 = self.deltri[et][2]
+            if v0 in bv or v1 in bv or v2 in bv:
+                btri.append(et)
+
+        self.boundarytriangle = btri
+  
+        #print "exterior triangles: ", extrind
+
+       
+    def _remove_holes(self,small):
         """
         Given the edges in bdry, find the largest exterior components.
         """
 
-        print "running _remove_holes \n"
+        #print "running _remove_holes \n"
+
+        bdry = self.boundary
         
         def findroot(i):
@@ -519 +476 @@
             vdiscard = [verts[k] for k in littlecomp] 
             newbdry = [e for e in bdry if not((e[0] in vdiscard) or (e[1] in vdiscard))]
+            # remove boundary triangles that touch discarded components
+            newbt = []
+            for bt in self.boundarytriangle:
+                v0 = self.deltri[bt][0]
+                v1 = self.deltri[bt][1]
+                v2 = self.deltri[bt][2]
+                if not (v0 in vdiscard or v1 in vdiscard or v2 in vdiscard):
+                    newbt.append(bt)
+            self.boundarytriangle = newbt
         else:
             newbdry = bdry
-
+    
         return newbdry
 
 
-
-    def _smooth_indents(self,bdry):
+    def _smooth_indents(self):
         """
         Given edges in bdry, test for acute-angle indents and remove them.
         """
-
-        print "running _smooth_indents \n"
-        
-        # 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.
+        
+        #print "running _smooth_indents \n"
+
+        bdry = self.boundary
+        bdrytri = self.boundarytriangle
         verts = self._vertices_from_edges(bdry)
-
-        bdrytri = []
-        for ind in extrind:
-            bvct = 0
+        
+        # find boundary triangles that have two edges in bdry
+        b2etri = []
+        for ind in bdrytri:
+            bect = 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
+                eda = (self.deltri[ind][(j+1)%3], self.deltri[ind][(j+2)%3])
+                edb = (self.deltri[ind][(j+2)%3], self.deltri[ind][(j+1)%3])               
+                if eda in bdry or edb in bdry:
+                    bect +=1
+            if bect==2:
+                b2etri.append(ind)
 
         # test the bdrytri triangles for acute angles 
-        probtri = []
-        for tind in bdrytri:
+        acutetri = []
+        for tind in b2etri:
             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]])
@@ -563 +522 @@
             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:
+                acutetri.append(tind)
+
+        #print "acute boundary triangles ", acutetri
+
+        # adjust the bdry edges and triangles by adding in the acutetri triangles
+        for pind in acutetri:
+            bdrytri.remove(pind) 
             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)       
-        
-        return bdry
-
-    def _expand_pinch(self,bdry):
+                bdry.append((tri[(i+1)%3], tri[(i+2)%3]))
+
+        self.boundarytriangle = bdrytri
+
+        newbdry = []
+        for ed in bdry:
+            numed = bdry.count(ed)+bdry.count((ed[1],ed[0]))
+            if numed%2 == 1:
+                newbdry.append(ed)
+        
+        #print "new boundary ", newbdry
+        return newbdry
+
+    def _expand_pinch(self):
         """
         Given edges in bdry, test for vertices with more than 2 incident edges.
         Expand by adding back in associated triangles... 
         """
-        print "running _expand_pinch \n"
-
-        #verts = self._vertices_from_edges(bdry)
-
+        #print "running _expand_pinch \n"
+
+        bdry = self.boundary
+        bdrytri = self.boundarytriangle
+        
         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
+        probv = [v[k] for k in range(len(v)) if (v[k]!=v[k-1] and v.count(v[k])>2) ]
+        #print "problem vertices: ", probv  
+
+        # find boundary triangles that have at least one vertex in probv
         probtri = []
-        for ind in extrind:
-            for j in [0,1,2]:
-                if (self.deltri[ind][j] in probv):
-                    tr = (ind,j)
-                    if probtri.count(tr)==0:
-                        probtri.append(tr)
-
-        # print "problem boundary triangle indices ", probtri
+        for ind in bdrytri:
+            v0 = self.deltri[ind][0]
+            v1 = self.deltri[ind][1]
+            v2 = self.deltri[ind][2]
+            if v0 in probv or v1 in probv or v2 in probv:
+                probtri.append(ind)
+
+        #print "problem boundary triangle indices ", probtri
 
         # "add in" the problem triangles
-        for t in probtri:
-            tri = self.deltri[t[0]] 
-            opeda = ( tri[(t[1]+1)%3], tri[(t[1]+2)%3] )
-            opedb = ( tri[(t[1]+2)%3], tri[(t[1]+1)%3] )
-            # oped is the edge opposite the problem vertex
-            # if oped is in bdry  remove it and the other triangle edges
-            # otherwise add oped to the bdry and remove the other triangle edges
-
-            if opeda in bdry:
-                bdry.remove(opeda)
-            elif opedb in bdry:
-                bdry.remove(opedb)
-            else:
-                bdry.append(opeda)
-
-            if (tri[t[1]], tri[(t[1]+1)%3]) in bdry:
-                bdry.remove((tri[t[1]], tri[(t[1]+1)%3]))
-            if (tri[(t[1]+1)%3], tri[t[1]]) in bdry:
-                bdry.remove((tri[(t[1]+1)%3], tri[t[1]]))
-            if (tri[t[1]], tri[(t[1]+2)%3]) in bdry:
-                bdry.remove((tri[t[1]], tri[(t[1]+2)%3]))
-            if (tri[(t[1]+2)%3], tri[t[1]]) in bdry:
-                bdry.remove((tri[(t[1]+2)%3], tri[t[1]]))
-              
-        return bdry
-   
+        for pind in probtri:
+            bdrytri.remove(pind) 
+            tri = self.deltri[pind]
+            for i in [0,1,2]:
+                bdry.append((tri[(i+1)%3], tri[(i+2)%3]))
+
+        self.boundarytriangle = bdrytri
+
+        newbdry = []
+        for ed in bdry:
+            numed = bdry.count(ed)+bdry.count((ed[1],ed[0]))
+            if numed%2 == 1:
+                newbdry.append(ed)
+        
+        #print "new boundary ", newbdry
+        return newbdry
+
+
 #-------------------------------------------------------------
 if __name__ == "__main__":

inundation/ga/storm_surge/alpha_shape/test_alpha_shape.py

@@ -107 +107 @@
 
         alpha = Alpha_Shape([a,b,c,cd,d,e,f])
-        alpha.set_boundary_type(1)
+        alpha.set_boundary_type(0,1,0,0)
         result = alpha.get_boundary()
         #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
+   
 
 
@@ -125 +125 @@
 
         alpha = Alpha_Shape([a,b,c,cd,d,e,f])
-        alpha.set_boundary_type(2)
+        alpha.set_boundary_type(0,0,1,0)
         result = alpha.get_boundary()
         #print "result",result
         answer = [(0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 6), (6, 0)]
         assert allclose(answer, result) 
-        pass
+
+
+    def test_boundary_3(self):
+        a = [452.0, -300.0]
+        b = [637.0, -322.0]
+        c = [844.0, -474.0]
+        d = [900.0, -585.0]
+        e = [726.0, -492.0]
+        f = [481.0, -354.0]
+        g = [744.0, -388.0]
+
+        alpha = Alpha_Shape([a,b,c,d,e,f,g])
+        alpha.set_boundary_type(0,0,0,1)
+        result = alpha.get_boundary()
+        answer = [(1, 0), (0, 5), (3, 2), (4, 3), (2, 6), (6, 1), (5, 4)]
+        assert allclose(answer, result)
+