Changeset 674

Timestamp:

Dec 6, 2004, 10:01:54 AM (20 years ago)

Author:

duncan

Message:

alpha_shape from vanessa

Location:

inundation/ga/storm_surge/alpha_shape

Files:

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

inundation/ga/storm_surge/alpha_shape/alpha_shape.py

@@ -3 +3 @@
 
 import exceptions
-from Numeric import array, Float
+from Numeric import array, Float, divide_safe, sqrt
 from load_mesh.loadASCII import load_xya_file, export_boundary_file
 
@@ -9 +9 @@
 
 OUTPUT_FILE_TITLE = "# The alpha shape boundary defined by point index pairs of edges"
+INF = pow(10,9)
 
 def alpha_shape_via_files(point_file, boundary_file, alpha= None):
@@ -25 +26 @@
     def __init__(self, points, alpha = None):
 
-        
-        """ Build interpolation matrix mapping from
-        function values at vertices to function values at data points
-
+        """ 
         Inputs:
         
@@ -37 +35 @@
         """
         self._set_points(points)
-        self._alpha_shape_algorithm()
+        self._alpha_shape_algorithm(alpha)
             
 
@@ -43 +41 @@
         """
         """
-        
+        # print "setting points"
         if len (points) <= 2:
             raise PointError, "Too few points to find an alpha shape"
@@ -50 +48 @@
         self.points = array(points).astype(Float)
 
-        if len (points) <= 2:
-            raise PointError, "Too few points to find an alpha shape"
-        
     
     def write_boundary(self,file_name):
@@ -58 +53 @@
         Write the boundary to a file
         """
-        #print " this info will be in the file",boundary
+        #print " this info will be in the file" 
         export_boundary_file(file_name, self.get_boundary(),
                              OUTPUT_FILE_TITLE, delimiter = ',')
@@ -65 +60 @@
         """
         """
+        #print "getting boundary" 
         return self.boundary
-    
-    def _alpha_shape_algorithm(self):
+
+    def get_delaunay(self):
+        """
+        """
+        return self.deltri
+    
+            
+    def _alpha_shape_algorithm(self,alpha):
         
         # A python style guide suggests using _[function name] to
-        # specify interal functions
-        #- this is the first time I've used it though - DSG
-        
+        # specify internal functions
+        # - this is the first time I've used it though - DSG
+
+        #print "starting alpha shape algorithm" 
+
+        # Build Delaunay triangulation
+        import triang
+        points = []
+        seglist = []
+        holelist = []
+        regionlist = []
+        pointattlist = []
+        segattlist = []
+        trilist = []
+
+        points = [(pt[0], pt[1]) for pt in self.points]
+        pointattlist = [ [] for i in range(len(points)) ] 
+        mode = "Qzcn"
+        #print "compute delaunay triangulation" 
+        tridata = triang.genMesh(points,seglist,holelist,regionlist,pointattlist,segattlist,trilist,mode)
+        #print tridata
+        #print "point attribute list: ", tridata['generatedpointattributelist'],"\n"
+        #print "hull segments: ", tridata['generatedsegmentlist'], "\n"
+        self.deltri = tridata['generatedtrianglelist']
+        self.deltrinbr = tridata['generatedtriangleneighborlist']
+        self.hulledges = tridata['generatedsegmentlist']                              
+
+        #print "deltri: ", self.deltri, "\n"
+        #print "deltrinbrs: ", self.deltrinbr, "\n"
+
+        # Build Alpha table
+        self._tri_circumradius()
+        self._edge_intervals()
+        self._vertex_intervals()
+
+        if alpha==None:
+            # Find optimum alpha
+            # Ken Clarkson's hull program uses smallest alpha so that
+            # every vertex is non-singular so... 
+            self.optimum_alpha = max([ interval[0] for interval in self.vertexinterval])
+            #print "optimum alpha ", self.optimum_alpha
+
+            #alpha_tri = self.get_alpha_triangles(self.optimum_alpha)
+            #print "alpha shape triangles ", alpha_tri
+
+            reg_edge = self.get_regular_edges(self.optimum_alpha)
+            #print "alpha boundary edges", reg_edge
+        else:
+            reg_edge = self.get_regular_edges(alpha)
+    
+        self.boundary = [self.edge[k] for k in reg_edge]
+        #print "alpha boundary edges ", self.boundary
+            
         # 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
-        
+        #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
+
+    def _tri_circumradius(self):
+
+        # compute circumradii of the delaunay triangles
+        x = self.points[:,0]
+        y = self.points[:,1]
+        ind1 = [self.deltri[j][0] for j in range(len(self.deltri))]
+        ind2 = [self.deltri[j][1] for j in range(len(self.deltri))]
+        ind3 = [self.deltri[j][2] for j in range(len(self.deltri))]
+
+        x1 = array([ x[j] for j in ind1 ])
+        y1 = array([ y[j] for j in ind1 ])
+        x2 = array([ x[j] for j in ind2 ])
+        y2 = array([ y[j] for j in ind2 ])
+        x3 = array([ x[j] for j in ind3 ])
+        y3 = array([ y[j] for j in ind3 ])
+
+        x21 = x2-x1
+        x31 = x3-x1
+        y21 = y2-y1
+        y31 = y3-y1
+
+        dist21 = x21*x21 + y21*y21
+        dist31 = x31*x31 + y31*y31
+
+        denom = x21*y31 - x31*y21
+        #print "denom = ", denom
+                
+        # dx/2, dy/2 give circumcenter relative to x1,y1. 
+        #dx = (y31*dist21 - y21*dist31)/denom
+        #dy = (x21*dist31 - x31*dist21)/denom
+        # from Numeric import divide_safe
+        try:
+            dx = divide_safe(y31*dist21 - y21*dist31,denom)
+            dy = divide_safe(x21*dist31 - x31*dist21,denom)
+        except ZeroDivisionError:
+            print "Found some degenerate 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)
+        #print "triangle radii", self.triradius
+        return
+
+    def _edge_intervals(self):
+        # 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.
+
+        # It should be possible to rewrite this routine in an array-friendly form
+        # like _tri_circumradius()  if we need to speed things up. 
+
+        edges = []
+        edgenbrs = []
+        edgeinterval = []
+        for t in range(len(self.deltri)):
+            tri = self.deltri[t]
+            trinbr = self.deltrinbr[t]
+            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]])
+            elen = sqrt(dx*dx+dy*dy)
+            #angle = array([(-dx[(i+1)%3]*dx[(i+2)%3]-dy[(i+1)%3]*dy[(i+2)%3])/(elen[(i+1)%3]*elen[(i+2)%3]) for i in [0,1,2]])
+            #angle = arccos(angle)
+            # really only need sign - not angle value:
+            anglesign = array([(-dx[(i+1)%3]*dx[(i+2)%3]-dy[(i+1)%3]*dy[(i+2)%3]) for i in [0,1,2]])
+            
+            #print "dx ",dx,"\n"
+            #print "dy ",dy,"\n"
+            #print "edge lengths of triangle ",t,"\t",elen,"\n"
+            #print "angles ",angle,"\n"
+                          
+            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" 
+
+    def _vertex_intervals(self):
+        # for each vertex find pairs
+        # (min_adj_triradius, max_adj_triradius)
+        nv = len(self.points)
+        vertexnbrs = [ [] for i in range(nv)]
+        vertexinterval = [ [] for i in range(nv)] 
+        for t in range(len(self.deltri)):
+            for j in self.deltri[t]:
+                vertexnbrs[j].append(t)
+        for h in range(len(self.hulledges)):
+            for j in self.hulledges[h]:
+                vertexnbrs[j].append(-1)
+        
+        for i in range(nv):
+            radii = [ self.triradius[t] for t in vertexnbrs[i] if t>=0 ]
+            vertexinterval[i] = [min(radii), max(radii)]
+            if vertexnbrs[i][-1]==-1:
+                vertexinterval[i][1]=INF
+
+        self.vertexnbr = vertexnbrs
+        self.vertexinterval = vertexinterval
+        #print "vertex nbrs ", vertexnbrs, "\n"
+        #print "vertex intervals ",vertexinterval, "\n"
+        
+    def get_alpha_triangles(self,alpha):
+        """
+        Given an alpha value,
+        return indices of triangles in the alpha-shape
+        """
+        def tri_rad_lta(k):
+            return self.triradius[k]<=alpha
+
+        return filter(tri_rad_lta, range(len(self.triradius)))
+
+    def get_regular_edges(self,alpha):
+        """
+        Given and alpha value,
+        return the edges on the boundary of the alpha-shape
+        """
+        def reg_edge(k):
+            return self.edgeinterval[k][1]<=alpha and self.edgeinterval[k][2]>alpha
+
+        return filter(reg_edge, range(len(self.edgeinterval)))
+
+        
+   
 #-------------------------------------------------------------
 if __name__ == "__main__":
@@ -88 +286 @@
 
     usage: alpha_shape.py point_file.xya boundary_file.bnd [alpha]
-
+    
     The alpha value is optional.
     """
     
     import os, sys
-    usage = "usage: %s point_file.xya boundary_file.bnd [alpha]"\
-            %os.path.basename(sys.argv[0])
+    usage = "usage: %s point_file.xya boundary_file.bnd [alpha]"%os.path.basename(sys.argv[0])
     # I made up the .bnd affix. Other ideas welcome. -DSG 
     if len(sys.argv) < 3:
@@ -105 +302 @@
         else:
             alpha = None
+
         alpha_shape_via_files(point_file, boundary_file, alpha)
         

inundation/ga/storm_surge/alpha_shape/test_alpha_shape.py

@@ -17 +17 @@
         pass
 
+    def test_delaunay(self):
+        #print "test_delaunay"
+        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]
 
-    def test_alpha(self):
+        alpha = Alpha_Shape([a,b,c,d,e,f])
+        result = alpha.get_delaunay()
+        answer = [(0, 1, 5), (5, 1, 4), (4, 2, 3), (2, 4, 1)]
+        # print "result",result
+        assert allclose(answer, result) 
+
+
+    def test_alpha_1(self):
+        #print "test_alpha" 
         a = [0.0, 0.0]
         b = [1.0, 0.0]
@@ -28 +44 @@
         alpha = Alpha_Shape([a,b,c,d,e,f])
         result = alpha.get_boundary()
-        answer = [[0, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 0]]
         #print "result",result
+        answer = [(5, 0), (0, 1), (4, 5), (2, 3), (3, 4), (1, 2)]
         assert allclose(answer, result) 
+
+
+    def test_alpha_2(self):
+        #print "test_alpha" 
+        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])
+        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) 
+
+
+    def test_alpha_3(self):
+        #print "test_alpha" 
+        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]
+        alf = 1.5
+
+        alpha = Alpha_Shape([a,b,c,d,e,f],alf)
+        result = alpha.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"
         a = [0.0, 0.0]
         b = [1.0, 0.0]
@@ -47 +99 @@
 
     def test_alpha_stand_alone(self):
+        #print "test_alpha_stand_alone"
         import os
         import tempfile
-       
+
         fileName = tempfile.mktemp(".xya")
+        #(h,fileName) = tempfile.mkstemp(".xya") 
         file = open(fileName,"w")
         file.write("title row \n\
@@ -62 +116 @@
 
         output_file_name = tempfile.mktemp(".bnd")
-        command = 'python .\\alpha_shape.py ' + fileName + ' ' + output_file_name
+        #(ho, output_file_name) = tempfile.mkstemp(".bnd")
+        command = 'python alpha_shape.py ' + fileName + ' ' + output_file_name
         #print "command", command
         os.system(command)
@@ -71 +126 @@
         file.close()
         os.remove(output_file_name)
-        self.failUnless(lFile[1] == "0,1" and
-                        lFile[2] == "1,2" and 
-                        lFile[3] == "2,3" and 
-                        lFile[4] == "3,4" and 
-                        lFile[5] == "4,5" and 
-                        lFile[6] == "5,0" and 
+        self.failUnless(lFile[1] == "5,0" and
+                        lFile[2] == "0,1" and 
+                        lFile[3] == "4,5" and 
+                        lFile[4] == "2,3" and 
+                        lFile[5] == "3,4" and 
+                        lFile[6] == "1,2" and 
                         lFile[7] == ""
                         ,
@@ -84 +139 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
+    #print "starting tests \n"
     suite = unittest.makeSuite(TestCase,'test')
     runner = unittest.TextTestRunner(verbosity=1)
     runner.run(suite)
+    #print "finished tests \n"
+