Changeset 691

Timestamp:

Dec 8, 2004, 9:26:14 AM (20 years ago)

Author:

duncan

Message:

updates from vanessa

Location:

inundation/ga/storm_surge/alpha_shape

Files:

• alpha_shape.py (modified) (14 diffs)
• test_alpha_shape.py (modified) (1 diff)

inundation/ga/storm_surge/alpha_shape/alpha_shape.py

@@ -5 +5 @@
 from Numeric import array, Float, divide_safe, sqrt
 from load_mesh.loadASCII import load_xya_file, export_boundary_file
+import random
 
 class PointError(exceptions.Exception): pass
@@ -44 +45 @@
         if len (points) <= 2:
             raise PointError, "Too few points to find an alpha shape"
+        if len(points)==3:
+            #check not in a straight line
+            x01 = points[0][0] - points[1][0]
+            y01 = points[0][1] - points[1][1]
+            x12 = points[1][0] - points[2][0]
+            y12 = points[1][1] - points[2][1]
+            crossprod = x01*y12 - x12*y01
+            if crossprod==0:
+                raise PointError, "Three points on a straight line"
         
         #Convert input to Numeric arrays
@@ -60 +70 @@
         """
         """
-        #print "getting boundary" 
         return self.boundary
 
@@ -67 +76 @@
         """
         return self.deltri
+
+    def get_optimum_alpha(self):
+        """
+        """
+        return self.optimum_alpha
+
+    def get_alpha(self):
+        """
+        Return current alpha value
+        """
+        return self.alpha
+
+    def set_alpha(self,alpha):
+        """
+        Set alpha and update alpha-boundary. 
+        """
+        self.alpha = alpha
+        reg_edge = self.get_regular_edges(alpha)    
+        self.boundary = [self.edge[k] for k in reg_edge]
     
             
@@ -76 +104 @@
 
         #print "starting alpha shape algorithm" 
+
+        self.alpha = alpha
 
         # Build Delaunay triangulation
@@ -90 +120 @@
         pointattlist = [ [] for i in range(len(points)) ] 
         mode = "Qzcn"
-        print "computing delaunay triangulation ... \n" 
+        #print "computing delaunay triangulation ... \n" 
         tridata = triang.genMesh(points,seglist,holelist,regionlist,pointattlist,segattlist,trilist,mode)
         #print tridata
@@ -99 +129 @@
         self.hulledges = tridata['generatedsegmentlist']                              
 
-        print "Number of delaunay triangles: ", len(self.deltri), "\n"
+        #print "Number of delaunay triangles: ", len(self.deltri), "\n"
         #print "deltrinbrs: ", self.deltrinbr, "\n"
 
         # Build Alpha table
-        print "Building alpha table ... \n" 
+        # print "Building alpha table ... \n" 
         self._tri_circumradius()
-        print "Largest circumradius ", max(self.triradius)
+        # print "Largest circumradius ", max(self.triradius)
         self._edge_intervals()
         self._vertex_intervals()
@@ -113 +143 @@
             # 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
-
+            self.optimum_alpha = max([ iv[0] for iv in self.vertexinterval if iv!=[] ])
+            # print "optimum alpha ", self.optimum_alpha
+            self.alpha = self.optimum_alpha
             #alpha_tri = self.get_alpha_triangles(self.optimum_alpha)
             #print "alpha shape triangles ", alpha_tri
@@ -161 +191 @@
         denom = x21*y31 - x31*y21
         #print "denom = ", denom
-                
+        delta = 0.00000001
+        random.seed()
+        zeroind = [k for k in range(len(denom)) if denom[k]==0 ]
+        while zeroind!=[]:
+            print "Warning: degenerate triangles found in alpha_shape.py, results may be inaccurate."
+            for d in zeroind:
+                x1[d] = x1[d]+delta*(random.random()-0.5)
+                x2[d] = x2[d]+delta*(random.random()-0.5)
+                x3[d] = x3[d]+delta*(random.random()-0.5)
+                y1[d] = y1[d]+delta*(random.random()-0.5)
+                y2[d] = y2[d]+delta*(random.random()-0.5)
+                y3[d] = y3[d]+delta*(random.random()-0.5)
+            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
+            zeroind = [k for k in range(len(denom)) if denom[k]==0 ]
+
         # dx/2, dy/2 give circumcenter relative to x1,y1. 
         #dx = (y31*dist21 - y21*dist31)/denom
@@ -170 +220 @@
             dy = divide_safe(x21*dist31 - x31*dist21,denom)
         except ZeroDivisionError:
-            print "Found some degenerate triangles."
+            print "There are serious problems with the delaunay triangles"
             raise
             # really need to take care of cases when denom = 0.
@@ -181 +231 @@
         return
 
+ 
     def _edge_intervals(self):
         # for each edge, find triples
@@ -254 +305 @@
                     vertexinterval[i][1]=INF
             except ValueError:
-                print "Vertex %i is isolated?"%i
-                print "coords: ",self.points[i]
-                print "Vertex nbrs: ", vertexnbrs[i]
-                print "nbr radii: ",radii
-                vertexinterval[i] = [0,INF]
+                # print "Vertex %i is isolated?"%i
+                # print "coords: ",self.points[i]
+                # print "Vertex nbrs: ", vertexnbrs[i]
+                # print "nbr radii: ",radii
+                # vertexinterval[i] = [0,INF]
                 pass
 
@@ -278 +329 @@
     def get_regular_edges(self,alpha):
         """
-        Given and alpha value,
+        Given an alpha value,
         return the edges on the boundary of the alpha-shape
         """
@@ -313 +364 @@
             alpha = None
 
-        print "about to call alpha shape routine \n"
+        #print "about to call alpha shape routine \n"
         alpha_shape_via_files(point_file, boundary_file, alpha)
         

inundation/ga/storm_surge/alpha_shape/test_alpha_shape.py

@@ -29 +29 @@
         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_3_points_on_line(self):
+        #print "test_delaunay"
+        a = [0.0, 0.0]
+        b = [1.0, 0.0]
+        c = [2.0, 0.0]
+
+        try:
+            alpha = Alpha_Shape([a,b,c])
+        except PointError:
+            pass
+        else:
+            self.failUnless(0==1, \
+                        'point list with 2 points did not raise an error!')
+