Changeset 1435
 Timestamp:
 May 19, 2005, 4:39:30 PM (19 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

inundation/ga/storm_surge/alpha_shape/alpha_shape.py
r1433 r1435 15 15 we now straighten all ``round'' faces to triangles and line segments, 16 16 we have an intuitive description of what is called the alphashape. 17 18 Author: Vanessa Robins, ANU 17 19 """ 18 20 … … 28 30 OUTPUT_FILE_TITLE = "# The alpha shape boundary defined by point index pairs of edges" 29 31 INF = pow(10,9) 32 EPSILON = 1.0e12 30 33 31 34 def alpha_shape_via_files(point_file, boundary_file, alpha= None): … … 50 53 def __init__(self, points, alpha = None): 51 54 """ 52 An Alpha_Shape requires input of a set of points. Other class routines55 An Alpha_Shape requires input of a set of points. Other class routines 53 56 return the alpha shape boundary. 54 57 … … 109 112 """ 110 113 Use the flags to set constraints on the boundary: 111 raw_boundary Return raw boundary i.e. the regular edges of the alpha shape 114 raw_boundary Return raw boundary i.e. the regular edges of the 115 alpha shape. 112 116 remove_holes filter to remove small holes 113 117 (small is defined by boundary_points_fraction ) … … 188 192 mode = "Qzcn" 189 193 #print "computing delaunay triangulation ... \n" 190 tridata = triang.genMesh(points,seglist,holelist,regionlist,pointattlist,segattlist,trilist,mode) 194 tridata = triang.genMesh(points,seglist,holelist,regionlist, 195 pointattlist,segattlist,trilist,mode) 191 196 #print tridata 192 #print "point att ributelist: ", tridata['generatedpointattributelist'],"\n"197 #print "point attlist: ", tridata['generatedpointattributelist'],"\n" 193 198 #print "hull segments: ", tridata['generatedsegmentlist'], "\n" 194 199 self.deltri = tridata['generatedtrianglelist'] 195 200 self.deltrinbr = tridata['generatedtriangleneighborlist'] 196 self.hulledges = tridata['generatedsegmentlist'] 201 self.hulledges = tridata['generatedsegmentlist'] 197 202 198 203 #print "Number of delaunay triangles: ", len(self.deltri), "\n" … … 211 216 # Ken Clarkson's hull program uses smallest alpha so that 212 217 # every vertex is nonsingular so... 213 self.optimum_alpha = max([ iv[0] for iv in self.vertexinterval if iv!=[] ]) 218 self.optimum_alpha = max([iv[0] for iv in self.vertexinterval \ 219 if iv!=[] ]) 214 220 # print "optimum alpha ", self.optimum_alpha 215 self.alpha = self.optimum_alpha 216 #alpha_tri = self.get_alpha_triangles(self.optimum_alpha) 217 #print "alpha shape triangles ", alpha_tri 218 else: 219 self.alpha = alpha 221 alpha = self.optimum_alpha 222 self.alpha = alpha 220 223 reg_edge = self.get_regular_edges(self.alpha) 221 222 224 self.boundary = [self.edge[k] for k in reg_edge] 223 225 #print "alpha boundary edges ", self.boundary 224 self._init_boundary_triangles() 225 226 self._init_boundary_triangles() 226 227 return 227 228 … … 260 261 # first need to check for nearzero values of denom 261 262 delta = 0.00000001 262 # FIXME checking two real numbers 263 zeroind = [k for k in range(len(denom)) if denom[k]==0 ] 264 # if some denom values are zero, we perturb the associated vertices and recalculate 263 zeroind = [k for k in range(len(denom)) if \ 264 (denom[k]< EPSILON and denom[k] > EPSILON)] 265 # if some denom values are close to zero, 266 # we perturb the associated vertices and recalculate 265 267 while zeroind!=[]: 266 268 random.seed() … … 280 282 dist31 = x31*x31 + y31*y31 281 283 denom = x21*y31  x31*y21 282 # FIXME checking two real numbers 283 zeroind = [k for k in range(len(denom)) if denom[k]==0 ] 284 284 zeroind = [k for k in range(len(denom)) if \ 285 (denom[k]< EPSILON and denom[k] > EPSILON)] 285 286 try: 286 287 dx = divide_safe(y31*dist21  y21*dist31,denom) … … 461 462 return i 462 463 k = findroot(vptr[i]) 463 vptr[i] = k # this produces "path compression" in the unionfind tree. 464 vptr[i] = k # this produces "path compression" in the 465 # unionfind tree. 464 466 return k 465 467 … … 511 513 vptr[vr] = rvl 512 514 #print "vptr: ", vptr, "\n" 513 514 515 # end edge loop 515 516 … … 528 529 529 530 # littleind has root indices for small components 530 littleind = [k for k in range(len(vptr)) if (vptr[k]<0 and vptr[k]>cutoff)] 531 littleind = [k for k in range(len(vptr)) if \ 532 (vptr[k]<0 and vptr[k]>cutoff)] 531 533 if littleind: 532 534 # littlecomp has all vptr indices in the small components 533 littlecomp = [k for k in range(len(vptr)) if findroot(k) in littleind] 535 littlecomp = [k for k in range(len(vptr)) if \ 536 findroot(k) in littleind] 534 537 # vdiscard has the vertex indices corresponding to vptr indices 535 538 vdiscard = [verts[k] for k in littlecomp] 536 newbdry = [e for e in bdry if not((e[0] in vdiscard) and (e[1] in vdiscard))] 539 newbdry = [e for e in bdry if \ 540 not((e[0] in vdiscard) and (e[1] in vdiscard))] 537 541 538 542 newverts = self._vertices_from_edges(newbdry) … … 546 550 newbt.append(bt) 547 551 548 self.boundarytriangle = newbt 549 552 self.boundarytriangle = newbt 550 553 else: 551 554 newbdry = bdry … … 556 559 def _smooth_indents(self): 557 560 """ 558 Given edges in bdry, test for acuteangle triangular indents and remove them. 561 Given edges in bdry, test for acuteangle triangular indents 562 and remove them. 559 563 """ 560 564 … … 576 580 for j in [0,1,2]: 577 581 eda = (self.deltri[ind][(j+1)%3], self.deltri[ind][(j+2)%3]) 578 edb = (self.deltri[ind][(j+2)%3], self.deltri[ind][(j+1)%3]) 582 edb = (self.deltri[ind][(j+2)%3], self.deltri[ind][(j+1)%3]) 579 583 if eda in bdry or edb in bdry: 580 584 bect +=1 … … 588 592 tri = self.deltri[tind[0]] 589 593 590 dx = array([self.points[tri[(i+1)%3],0]  self.points[tri[(i+2)%3],0] for i in [0,1,2]]) 591 dy = array([self.points[tri[(i+1)%3],1]  self.points[tri[(i+2)%3],1] for i in [0,1,2]]) 592 anglesign = array([(dx[(i+1)%3]*dx[(i+2)%3]dy[(i+1)%3]*dy[(i+2)%3]) for i in [0,1,2]]) 593 # record any triangle that has an acute angle spanned by two edges along the boundary.. 594 dx = array([self.points[tri[(i+1)%3],0]  \ 595 self.points[tri[(i+2)%3],0] for i in [0,1,2]]) 596 dy = array([self.points[tri[(i+1)%3],1]  \ 597 self.points[tri[(i+2)%3],1] for i in [0,1,2]]) 598 anglesign = array([(dx[(i+1)%3]*dx[(i+2)%3]\ 599 dy[(i+1)%3]*dy[(i+2)%3]) for i in [0,1,2]]) 600 # record any triangle that has an acute angle spanned by 601 #two edges along the boundary.. 594 602 if anglesign[tind[1]] > 0: 595 603 acutetri.append(tind[0]) … … 597 605 #print "acute boundary triangles ", acutetri 598 606 599 # adjust the bdry edges and triangles by adding in the acutetri triangles 607 # adjust the bdry edges and triangles by adding 608 #in the acutetri triangles 600 609 for pind in acutetri: 601 610 bdrytri.remove(pind) … … 627 636 v = v1+v2 628 637 v.sort() 629 probv = [v[k] for k in range(len(v)) if (v[k]!=v[k1] and v.count(v[k])>2) ] 638 probv = [v[k] for k in range(len(v)) \ 639 if (v[k]!=v[k1] and v.count(v[k])>2) ] 630 640 #print "problem vertices: ", probv 631 641
Note: See TracChangeset
for help on using the changeset viewer.