Changeset 1433
 Timestamp:
 May 19, 2005, 4:18:07 PM (20 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

inundation/ga/storm_surge/alpha_shape/alpha_shape.py
r1432 r1433 3 3 4 4 From website by Kaspar Fischer: 5 As mentionned in Edelsbrunner's and Mu cke's paper, one can5 As mentionned in Edelsbrunner's and Muecke's paper, one can 6 6 intuitively think of an alphashape as the following: 7 7 … … 42 42 points = point_dict['pointlist'] 43 43 44 #alpha = Alpha_Shape(points)45 #alpha.write_boundary(boundary_file)46 44 AS = Alpha_Shape(points, alpha) 47 45 AS.write_boundary(boundary_file) … … 61 59 """ 62 60 self._set_points(points) 63 #self._alpha_shape_algorithm(alpha=alpha)64 61 self._alpha_shape_algorithm(alpha) 65 62 … … 75 72 if len(points)==3: 76 73 #check not in a straight line 74 # FIXME check points 1,2,3 if straingt, check if points 2,3,4, ect 77 75 x01 = points[0][0]  points[1][0] 78 76 y01 = points[0][1]  points[1][1] … … 218 216 #alpha_tri = self.get_alpha_triangles(self.optimum_alpha) 219 217 #print "alpha shape triangles ", alpha_tri 220 221 reg_edge = self.get_regular_edges(self.optimum_alpha)222 #print "alpha boundary edges", reg_edge223 218 else: 224 reg_edge = self.get_regular_edges(alpha) 219 self.alpha = alpha 220 reg_edge = self.get_regular_edges(self.alpha) 225 221 226 222 self.boundary = [self.edge[k] for k in reg_edge] … … 264 260 # first need to check for nearzero values of denom 265 261 delta = 0.00000001 266 random.seed()262 # FIXME checking two real numbers 267 263 zeroind = [k for k in range(len(denom)) if denom[k]==0 ] 268 # if some denom values are close tozero, we perturb the associated vertices and recalculate264 # if some denom values are zero, we perturb the associated vertices and recalculate 269 265 while zeroind!=[]: 266 random.seed() 270 267 print "Warning: degenerate triangles found in alpha_shape.py, results may be inaccurate." 271 268 for d in zeroind: … … 283 280 dist31 = x31*x31 + y31*y31 284 281 denom = x21*y31  x31*y21 282 # FIXME checking two real numbers 285 283 zeroind = [k for k in range(len(denom)) if denom[k]==0 ] 286 284 287 288 # from Numeric import divide_safe289 285 try: 290 286 dx = divide_safe(y31*dist21  y21*dist31,denom) 291 287 dy = divide_safe(x21*dist31  x31*dist21,denom) 292 288 except ZeroDivisionError: 293 print "There are serious problems with the delaunay triangles" 294 raise #DSG maybe don't catch this error? Is recovery possible? 295 # VR  we've already tried to deal with these errors above... 296 # but should probably write some better code! 289 raise AlphaError 297 290 298 291 self.triradius = 0.5*sqrt(dx*dx + dy*dy) 299 292 #print "triangle radii", self.triradius 300 301 302 303 293 304 294 def _edge_intervals(self): … … 311 301 312 302 # It should be possible to rewrite this routine in an arrayfriendly 313 # form like _tri_circumradius() if we need to speed things up. 303 # form like _tri_circumradius() if we need to speed things up. 304 # Hard to do though. 314 305 315 306 edges = [] … … 324 315 self.points[tri[(i+2)%3],1] for i in [0,1,2]]) 325 316 elen = sqrt(dx*dx+dy*dy) 326 #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]])327 #angle = arccos(angle)328 317 # really only need sign  not angle value: 329 318 anglesign = array([(dx[(i+1)%3]*dx[(i+2)%3] … … 350 339 else: 351 340 continue 352 #if angle[i]>pi/2:353 341 if anglesign[i] < 0: 354 342 edgeinterval[1][0] = edgeinterval[1][1] … … 383 371 vertexinterval[i][1]=INF 384 372 except ValueError: 385 # print "Vertex %i is isolated?"%i 386 # print "coords: ",self.points[i] 387 # print "Vertex nbrs: ", vertexnbrs[i] 388 # print "nbr radii: ",radii 389 # vertexinterval[i] = [0,INF] 390 391 #DSG is it ok to catch this error and continue? 392 pass 373 raise AlphaError 393 374 394 375 self.vertexnbr = vertexnbrs … … 623 604 bdry.append((tri[(i+1)%3], tri[(i+2)%3])) 624 605 625 self.boundarytriangle = bdrytri626 627 606 newbdry = [] 628 607 for ed in bdry: … … 668 647 for i in [0,1,2]: 669 648 bdry.append((tri[(i+1)%3], tri[(i+2)%3])) 670 671 self.boundarytriangle = bdrytri672 649 673 650 newbdry = []
Note: See TracChangeset
for help on using the changeset viewer.