Changeset 1435 for inundation/ga/storm_surge/alpha_shape/alpha_shape.py

Timestamp:

May 19, 2005, 4:39:30 PM (19 years ago)

Author:

duncan

Message:

code format changes

File:

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

inundation/ga/storm_surge/alpha_shape/alpha_shape.py

@@ -15 +15 @@
 we now straighten all ``round'' faces to triangles and line segments,
 we have an intuitive description of what is called the alpha-shape.
+
+Author: Vanessa Robins, ANU
 """
 
@@ -28 +30 @@
 OUTPUT_FILE_TITLE = "# The alpha shape boundary defined by point index pairs of edges"
 INF = pow(10,9)
+EPSILON = 1.0e-12
 
 def alpha_shape_via_files(point_file, boundary_file, alpha= None):
@@ -50 +53 @@
     def __init__(self, points, alpha = None):
         """
-        An Alpha_Shape requires input of a set of points.  Other class routines 
+        An Alpha_Shape requires input of a set of points. Other class routines
         return the alpha shape boundary.
         
@@ -109 +112 @@
         """
         Use the flags to set constraints on the boundary:
-        raw_boundary    Return raw boundary i.e. the regular edges of the  alpha shape 
+        raw_boundary    Return raw boundary i.e. the regular edges of the
+                        alpha shape. 
         remove_holes    filter to remove small holes
                         (small is defined by  boundary_points_fraction ) 
@@ -188 +192 @@
         mode = "Qzcn"
         #print "computing delaunay triangulation ... \n" 
-        tridata = triang.genMesh(points,seglist,holelist,regionlist,pointattlist,segattlist,trilist,mode)
+        tridata = triang.genMesh(points,seglist,holelist,regionlist,
+                                 pointattlist,segattlist,trilist,mode)
         #print tridata
-        #print "point attribute list: ", tridata['generatedpointattributelist'],"\n"
+        #print "point attlist: ", tridata['generatedpointattributelist'],"\n"
         #print "hull segments: ", tridata['generatedsegmentlist'], "\n"
         self.deltri = tridata['generatedtrianglelist']
         self.deltrinbr = tridata['generatedtriangleneighborlist']
-        self.hulledges = tridata['generatedsegmentlist']                              
+        self.hulledges = tridata['generatedsegmentlist'] 
 
         #print "Number of delaunay triangles: ", len(self.deltri), "\n"
@@ -211 +216 @@
             # Ken Clarkson's hull program uses smallest alpha so that
             # every vertex is non-singular so... 
-            self.optimum_alpha = max([ iv[0] for iv in self.vertexinterval if iv!=[] ])
+            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
-        else:
-            self.alpha = alpha
+            alpha = self.optimum_alpha
+        self.alpha = alpha
         reg_edge = self.get_regular_edges(self.alpha)
-    
         self.boundary = [self.edge[k] for k in reg_edge]
         #print "alpha boundary edges ", self.boundary
-        self._init_boundary_triangles()
-            
+        self._init_boundary_triangles()            
         return
 
@@ -260 +261 @@
         # first need to check for near-zero values of denom 
         delta = 0.00000001
-        # FIXME checking two real numbers
-        zeroind = [k for k in range(len(denom)) if denom[k]==0 ]
-        # if some denom values are zero, we perturb the associated vertices and recalculate 
+        zeroind = [k for k in range(len(denom)) if \
+                   (denom[k]< EPSILON and  denom[k] > -EPSILON)]
+        # if some denom values are close to zero,
+        # we perturb the associated vertices and recalculate 
         while zeroind!=[]:
             random.seed()
@@ -280 +282 @@
             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 ]
-
+            zeroind = [k for k in range(len(denom)) if \
+                       (denom[k]< EPSILON and  denom[k] > -EPSILON)]
         try:
             dx = divide_safe(y31*dist21 - y21*dist31,denom)
@@ -461 +462 @@
                 return i
             k = findroot(vptr[i])
-            vptr[i] = k    # this produces "path compression" in the union-find tree. 
+            vptr[i] = k    # this produces "path compression" in the
+                           # union-find tree. 
             return k
 
@@ -511 +513 @@
                             vptr[vr] = rvl 
                 #print "vptr: ", vptr, "\n"
-
         # end edge loop
 
@@ -528 +529 @@
 
         # littleind has root indices for small components 
-        littleind = [k for k in range(len(vptr)) if (vptr[k]<0 and vptr[k]>-cutoff)] 
+        littleind = [k for k in range(len(vptr)) if \
+                     (vptr[k]<0 and vptr[k]>-cutoff)] 
         if littleind:
             # littlecomp has all vptr indices in the small components
-            littlecomp = [k for k in range(len(vptr)) if findroot(k) in littleind]
+            littlecomp = [k for k in range(len(vptr)) if \
+                          findroot(k) in littleind]
             # vdiscard has the vertex indices corresponding to vptr indices  
             vdiscard = [verts[k] for k in littlecomp] 
-            newbdry = [e for e in bdry if not((e[0] in vdiscard) and (e[1] in vdiscard))]
+            newbdry = [e for e in bdry if \
+                       not((e[0] in vdiscard) and (e[1] in vdiscard))]
 
             newverts = self._vertices_from_edges(newbdry)
@@ -546 +550 @@
                     newbt.append(bt)
 
-            self.boundarytriangle = newbt
-            
+            self.boundarytriangle = newbt           
         else:
             newbdry = bdry
@@ -556 +559 @@
     def _smooth_indents(self):
         """
-        Given edges in bdry, test for acute-angle triangular indents and remove them.
+        Given edges in bdry, test for acute-angle triangular indents
+        and remove them.
         """
         
@@ -576 +580 @@
             for j in [0,1,2]:
                 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])               
+                edb = (self.deltri[ind][(j+2)%3], self.deltri[ind][(j+1)%3])
                 if eda in bdry or edb in bdry:
                     bect +=1
@@ -588 +592 @@
             tri = self.deltri[tind[0]]
             
-            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]])
-            anglesign = array([(-dx[(i+1)%3]*dx[(i+2)%3]-dy[(i+1)%3]*dy[(i+2)%3]) for i in [0,1,2]])
-            # record any triangle that has an acute angle spanned by two edges along the boundary..
+            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]])
+            anglesign = array([(-dx[(i+1)%3]*dx[(i+2)%3]-\
+                                dy[(i+1)%3]*dy[(i+2)%3]) for i in [0,1,2]])
+            # record any triangle that has an acute angle spanned by
+            #two edges along the boundary..
             if anglesign[tind[1]] > 0:
                 acutetri.append(tind[0])
@@ -597 +605 @@
         #print "acute boundary triangles ", acutetri
 
-        # adjust the bdry edges and triangles by adding in the acutetri triangles
+        # adjust the bdry edges and triangles by adding
+        #in the acutetri triangles
         for pind in acutetri:
             bdrytri.remove(pind) 
@@ -627 +636 @@
         v = v1+v2
         v.sort()
-        probv = [v[k] for k in range(len(v)) if (v[k]!=v[k-1] and v.count(v[k])>2) ]
+        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