Changeset 1420

Timestamp:

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

Author:

duncan

Message:

dsg added comments and moved some code around

Location:

inundation/ga/storm_surge/alpha_shape

Files:

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

inundation/ga/storm_surge/alpha_shape/alpha_shape.py

@@ -1 @@
-"""Alpha shape   
+"""Alpha shape
+Determine the shape of a set of points.
+
+As mentionned in Edelsbrunner's and Mucke's paper, one can
+intuitively think of an alpha-shape as the following:
+test
+Imagine a huge mass of ice-cream making up the space and containing
+the points S as ``hard'' chocolate pieces. Using one of these
+sphere-formed ice-cream spoons we carve out all parts of the ice-cream
+block we can reach without bumping into chocolate pieces, thereby even
+carving out holes in the inside (eg. parts not reachable by simply
+moving the spoon from the outside). We will eventually end up with a
+(not necessarily convex) object bounded by caps, arcs and points. If
+we now straighten all ``round'' faces to triangles and line segments,
+we have an intuitive description of what is called the alpha-shape.
 """
 
@@ -14 +28 @@
 
 def alpha_shape_via_files(point_file, boundary_file, alpha= None):
-    
+    """
+    Load a point file and return the alpha shape boundary as a boundary file.
+    
+    Inputs:
+    point_file: File location of the input file, points format (.xya or .pts)
+    boundary_file: File location of the generated output file
+    alpha: The alpha value can be optionally specified.  If it is not specified
+    the optimum alpha value will be used.
+"""
     point_dict = import_points_file(point_file)
     points = point_dict['pointlist']
-    #title_string = point_dict['title']
     
     alpha = Alpha_Shape(points)
@@ -25 +46 @@
 
     def __init__(self, points, alpha = None):
-
-        """ 
-        Inputs:
-        
+        """
+        alpha_shape inputs a set of points and can return the alpha
+        shape boundary.
+        
+        Inputs:       
           points: List of coordinate pairs [[x1, y1],[x2, y2]..] 
-
-          alpha: alpha shape parameter
-          
+          alpha: The alpha value can be optionally specified.  If it is
+          not specified the optimum alpha value will be used.
         """
         self._set_points(points)
-        self._alpha_shape_algorithm(alpha)
+        self._alpha_shape_algorithm(alpha=alpha)
 
     def _set_points(self, points):
         """
+        Create self.points array, do Error checking
+        Inputs:       
+          points: List of coordinate pairs [[x1, y1],[x2, y2]..] 
         """
         # print "setting points"
@@ -67 +91 @@
     def get_boundary(self):
         """
+        Return a list of tuples.  Each tuples represents a segment,
+        by giving the index of the segments two end points.
+        The list of tuple represents the alpha shape.
         """
         return self.boundary
@@ -76 +103 @@
                           boundary_points_fraction=0.2):
         """
-        Use the flags to set constraints on the boundary
+        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
@@ -111 +138 @@
     def get_alpha(self):
         """
-        Return current alpha value
+        Return current alpha value.
         """
         return self.alpha
@@ -125 +152 @@
     
             
-    def _alpha_shape_algorithm(self,alpha):
-        
-        # A python style guide suggests using _[function name] to
-        # specify internal functions
-        # - this is the first time I've used it though - DSG
+    def _alpha_shape_algorithm(self, alpha=None):
+        """
+        Given a set of points (self.points) and an optional alpha value
+        determines the alpha shape (stored in self.boundary,
+        accessed by get_boundary).
+        
+        Inputs:
+          alpha: The alpha value can be optionally specified.  If it is
+          not specified the optimum alpha value will be used.
+        """
 
         #print "starting alpha shape algorithm" 
@@ -189 +221 @@
 
     def _tri_circumradius(self):
-
-        # compute circumradii of the delaunay triangles
+        """
+        Compute circumradii of the delaunay triangles
+        """
+        
         x = self.points[:,0]
         y = self.points[:,1]
@@ -244 +278 @@
         except ZeroDivisionError:
             print "There are serious problems with the delaunay triangles"
-            raise
+            raise  #DSG maybe don't catch this error? Is recovery possible?
  
             
         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. 
+        """
+         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 = []
@@ -267 +303 @@
             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]])
+            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]])
+            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"
@@ -291 +330 @@
                     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]]) ])
+                       min(self.triradius[t],self.triradius[trinbr[i]]),\
+                          max(self.triradius[t],self.triradius[trinbr[i]]) ])
                 else:
                     continue
@@ -307 +346 @@
 
     def _vertex_intervals(self):
-        # for each vertex find pairs
-        # (min_adj_triradius, max_adj_triradius)
+        """
+        for each vertex find pairs
+        (min_adj_triradius, max_adj_triradius)
+        """
         nv = len(self.points)
         vertexnbrs = [ [] for i in range(nv)]
@@ -331 +372 @@
                 # print "nbr radii: ",radii
                 # vertexinterval[i] = [0,INF]
+
+                #DSG is it ok to catch this error and continue?
                 pass
 
@@ -354 +397 @@
         """
         def reg_edge(k):
-            return self.edgeinterval[k][1]<=alpha and self.edgeinterval[k][2]>alpha
+            return self.edgeinterval[k][1]<=alpha and \
+                   self.edgeinterval[k][2]>alpha
 
         return filter(reg_edge, range(len(self.edgeinterval)))
@@ -364 +408 @@
         """
         def exp_vert(k):
-            return self.vertexinterval[k][0]<=alpha and self.vertexinterval[k][1]>alpha
+            return self.vertexinterval[k][0]<=alpha and \
+                   self.vertexinterval[k][1]>alpha
 
         return filter(exp_vert, range(len(self.vertexinterval)))        
@@ -562 +607 @@
             # an acute angle spanned by two edges along the boundary..
             # MUST FIX.
+            #DSG - should this be fixed?
             if anglesign[tind[1]] > 0:
                 acutetri.append(tind[0])
@@ -646 +692 @@
     import os, sys
     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:
         print usage

inundation/ga/storm_surge/alpha_shape/test_alpha_shape.py

@@ -346 +346 @@
 (126.437678428,28.0083464873), \
 (133.639824668,50.0149044415), \
-(176.852702105,43.612996673), \alpha = Alpha_Shape(
+(176.852702105,43.612996673), \
 (123.636843779,-24.4072733675), \
 (73.6219393379,35.2104927268), \