Changeset 691
 Timestamp:
 Dec 8, 2004, 9:26:14 AM (19 years ago)
 Location:
 inundation/ga/storm_surge/alpha_shape
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

inundation/ga/storm_surge/alpha_shape/alpha_shape.py
r684 r691 5 5 from Numeric import array, Float, divide_safe, sqrt 6 6 from load_mesh.loadASCII import load_xya_file, export_boundary_file 7 import random 7 8 8 9 class PointError(exceptions.Exception): pass … … 44 45 if len (points) <= 2: 45 46 raise PointError, "Too few points to find an alpha shape" 47 if len(points)==3: 48 #check not in a straight line 49 x01 = points[0][0]  points[1][0] 50 y01 = points[0][1]  points[1][1] 51 x12 = points[1][0]  points[2][0] 52 y12 = points[1][1]  points[2][1] 53 crossprod = x01*y12  x12*y01 54 if crossprod==0: 55 raise PointError, "Three points on a straight line" 46 56 47 57 #Convert input to Numeric arrays … … 60 70 """ 61 71 """ 62 #print "getting boundary"63 72 return self.boundary 64 73 … … 67 76 """ 68 77 return self.deltri 78 79 def get_optimum_alpha(self): 80 """ 81 """ 82 return self.optimum_alpha 83 84 def get_alpha(self): 85 """ 86 Return current alpha value 87 """ 88 return self.alpha 89 90 def set_alpha(self,alpha): 91 """ 92 Set alpha and update alphaboundary. 93 """ 94 self.alpha = alpha 95 reg_edge = self.get_regular_edges(alpha) 96 self.boundary = [self.edge[k] for k in reg_edge] 69 97 70 98 … … 76 104 77 105 #print "starting alpha shape algorithm" 106 107 self.alpha = alpha 78 108 79 109 # Build Delaunay triangulation … … 90 120 pointattlist = [ [] for i in range(len(points)) ] 91 121 mode = "Qzcn" 92 print "computing delaunay triangulation ... \n"122 #print "computing delaunay triangulation ... \n" 93 123 tridata = triang.genMesh(points,seglist,holelist,regionlist,pointattlist,segattlist,trilist,mode) 94 124 #print tridata … … 99 129 self.hulledges = tridata['generatedsegmentlist'] 100 130 101 print "Number of delaunay triangles: ", len(self.deltri), "\n"131 #print "Number of delaunay triangles: ", len(self.deltri), "\n" 102 132 #print "deltrinbrs: ", self.deltrinbr, "\n" 103 133 104 134 # Build Alpha table 105 print "Building alpha table ... \n"135 # print "Building alpha table ... \n" 106 136 self._tri_circumradius() 107 print "Largest circumradius ", max(self.triradius)137 # print "Largest circumradius ", max(self.triradius) 108 138 self._edge_intervals() 109 139 self._vertex_intervals() … … 113 143 # Ken Clarkson's hull program uses smallest alpha so that 114 144 # every vertex is nonsingular so... 115 self.optimum_alpha = max([ i nterval[0] for interval in self.vertexinterval])116 print "optimum alpha ", self.optimum_alpha117 145 self.optimum_alpha = max([ iv[0] for iv in self.vertexinterval if iv!=[] ]) 146 # print "optimum alpha ", self.optimum_alpha 147 self.alpha = self.optimum_alpha 118 148 #alpha_tri = self.get_alpha_triangles(self.optimum_alpha) 119 149 #print "alpha shape triangles ", alpha_tri … … 161 191 denom = x21*y31  x31*y21 162 192 #print "denom = ", denom 163 193 delta = 0.00000001 194 random.seed() 195 zeroind = [k for k in range(len(denom)) if denom[k]==0 ] 196 while zeroind!=[]: 197 print "Warning: degenerate triangles found in alpha_shape.py, results may be inaccurate." 198 for d in zeroind: 199 x1[d] = x1[d]+delta*(random.random()0.5) 200 x2[d] = x2[d]+delta*(random.random()0.5) 201 x3[d] = x3[d]+delta*(random.random()0.5) 202 y1[d] = y1[d]+delta*(random.random()0.5) 203 y2[d] = y2[d]+delta*(random.random()0.5) 204 y3[d] = y3[d]+delta*(random.random()0.5) 205 x21 = x2x1 206 x31 = x3x1 207 y21 = y2y1 208 y31 = y3y1 209 dist21 = x21*x21 + y21*y21 210 dist31 = x31*x31 + y31*y31 211 denom = x21*y31  x31*y21 212 zeroind = [k for k in range(len(denom)) if denom[k]==0 ] 213 164 214 # dx/2, dy/2 give circumcenter relative to x1,y1. 165 215 #dx = (y31*dist21  y21*dist31)/denom … … 170 220 dy = divide_safe(x21*dist31  x31*dist21,denom) 171 221 except ZeroDivisionError: 172 print " Found some degenerate triangles."222 print "There are serious problems with the delaunay triangles" 173 223 raise 174 224 # really need to take care of cases when denom = 0. … … 181 231 return 182 232 233 183 234 def _edge_intervals(self): 184 235 # for each edge, find triples … … 254 305 vertexinterval[i][1]=INF 255 306 except ValueError: 256 print "Vertex %i is isolated?"%i257 print "coords: ",self.points[i]258 print "Vertex nbrs: ", vertexnbrs[i]259 print "nbr radii: ",radii260 vertexinterval[i] = [0,INF]307 # print "Vertex %i is isolated?"%i 308 # print "coords: ",self.points[i] 309 # print "Vertex nbrs: ", vertexnbrs[i] 310 # print "nbr radii: ",radii 311 # vertexinterval[i] = [0,INF] 261 312 pass 262 313 … … 278 329 def get_regular_edges(self,alpha): 279 330 """ 280 Given an dalpha value,331 Given an alpha value, 281 332 return the edges on the boundary of the alphashape 282 333 """ … … 313 364 alpha = None 314 365 315 print "about to call alpha shape routine \n"366 #print "about to call alpha shape routine \n" 316 367 alpha_shape_via_files(point_file, boundary_file, alpha) 317 368 
inundation/ga/storm_surge/alpha_shape/test_alpha_shape.py
r674 r691 29 29 result = alpha.get_delaunay() 30 30 answer = [(0, 1, 5), (5, 1, 4), (4, 2, 3), (2, 4, 1)] 31 # print "result",result32 31 assert allclose(answer, result) 32 33 def test_3_points_on_line(self): 34 #print "test_delaunay" 35 a = [0.0, 0.0] 36 b = [1.0, 0.0] 37 c = [2.0, 0.0] 38 39 try: 40 alpha = Alpha_Shape([a,b,c]) 41 except PointError: 42 pass 43 else: 44 self.failUnless(0==1, \ 45 'point list with 2 points did not raise an error!') 46 33 47 34 48
Note: See TracChangeset
for help on using the changeset viewer.