# Changeset 1420

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

dsg added comments and moved some code around

Location:
inundation/ga/storm_surge/alpha_shape
Files:
2 edited

Unmodified
Added
Removed
• ## inundation/ga/storm_surge/alpha_shape/alpha_shape.py

 r1381 """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. """ 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) 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" 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 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 def get_alpha(self): """ Return current alpha value Return current alpha value. """ return self.alpha 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" def _tri_circumradius(self): # compute circumradii of the delaunay triangles """ Compute circumradii of the delaunay triangles """ x = self.points[:,0] y = self.points[:,1] 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 = [] 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" 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 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)] # print "nbr radii: ",radii # vertexinterval[i] = [0,INF] #DSG is it ok to catch this error and continue? pass """ 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))) """ 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))) # 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]) 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

 r1364 (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), \
Note: See TracChangeset for help on using the changeset viewer.