Changeset 1433

Timestamp:

May 19, 2005, 4:18:07 PM (19 years ago)

Author:

duncan

Message:

changes from vanessa and duncan

File:

• inundation/ga/storm_surge/alpha_shape/alpha_shape.py (modified) (13 diffs)

inundation/ga/storm_surge/alpha_shape/alpha_shape.py

@@ -3 +3 @@
 
 From website by Kaspar Fischer:
-As mentionned in Edelsbrunner's and Mucke's paper, one can
+As mentionned in Edelsbrunner's and Muecke's paper, one can
 intuitively think of an alpha-shape as the following:
 
@@ -42 +42 @@
     points = point_dict['pointlist']
     
-    #alpha = Alpha_Shape(points)
-    #alpha.write_boundary(boundary_file)
     AS = Alpha_Shape(points, alpha)
     AS.write_boundary(boundary_file)
@@ -61 +59 @@
         """
         self._set_points(points)
-        #self._alpha_shape_algorithm(alpha=alpha)
         self._alpha_shape_algorithm(alpha)
 
@@ -75 +72 @@
         if len(points)==3:
             #check not in a straight line
+            # FIXME check points 1,2,3 if straingt, check if points 2,3,4, ect
             x01 = points[0][0] - points[1][0]
             y01 = points[0][1] - points[1][1]
@@ -218 +216 @@
             #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.alpha = alpha
+        reg_edge = self.get_regular_edges(self.alpha)
     
         self.boundary = [self.edge[k] for k in reg_edge]
@@ -264 +260 @@
         # first need to check for near-zero values of denom 
         delta = 0.00000001
-        random.seed()
+        # FIXME checking two real numbers
         zeroind = [k for k in range(len(denom)) if denom[k]==0 ]
-        # if some denom values are close to zero, we perturb the associated vertices and recalculate 
+        # if some denom values are zero, we perturb the associated vertices and recalculate 
         while zeroind!=[]:
+            random.seed()
             print "Warning: degenerate triangles found in alpha_shape.py, results may be inaccurate."
             for d in zeroind:
@@ -283 +280 @@
             dist31 = x31*x31 + y31*y31
             denom = x21*y31 - x31*y21
+            # FIXME checking two real numbers
             zeroind = [k for k in range(len(denom)) if denom[k]==0 ]
 
- 
-        # 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 "There are serious problems with the delaunay triangles"
-            raise  #DSG maybe don't catch this error? Is recovery possible?
-            # VR - we've already tried to deal with these errors above...
-            #      but should probably write some better code!
+            raise  AlphaError
             
         self.triradius = 0.5*sqrt(dx*dx + dy*dy)
         #print "triangle radii", self.triradius
-        
-
-
 
     def _edge_intervals(self):
@@ -311 +301 @@
         
         # It should be possible to rewrite this routine in an array-friendly
-        # form like _tri_circumradius()  if we need to speed things up. 
+        # form like _tri_circumradius()  if we need to speed things up.
+        # Hard to do though.
 
         edges = []
@@ -324 +315 @@
                         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]-
@@ -350 +339 @@
                 else:
                     continue
-                #if angle[i]>pi/2:
                 if anglesign[i] < 0: 
                     edgeinterval[-1][0] = edgeinterval[-1][1]
@@ -383 +371 @@
                     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]
-
-                #DSG is it ok to catch this error and continue?
-                pass
+                raise AlphaError
 
         self.vertexnbr = vertexnbrs
@@ -623 +604 @@
                 bdry.append((tri[(i+1)%3], tri[(i+2)%3]))
 
-        self.boundarytriangle = bdrytri
-
         newbdry = []
         for ed in bdry:
@@ -668 +647 @@
             for i in [0,1,2]:
                 bdry.append((tri[(i+1)%3], tri[(i+2)%3]))
-
-        self.boundarytriangle = bdrytri
 
         newbdry = []