Changeset 1364
- Timestamp:
- May 11, 2005, 11:49:06 AM (20 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
r1227 r1364 79 79 raw_boundary Return raw boundary i.e. the regular edges 80 80 remove_holes filter to remove small holes 81 smooth_indents remove sharp indents in boundary81 smooth_indents remove sharp triangular indents in boundary 82 82 expand_pinch plus test for pinch-off 83 83 i.e. boundary vertex with more than two edges. … … 407 407 def _remove_holes(self,small): 408 408 """ 409 Given the edges in bdry, find the largest exterior components. 409 Given the edges in self.boundary, finds the largest components. 410 The routine does this by implementing a union-find algorithm. 410 411 """ 411 412 … … 418 419 return i 419 420 k = findroot(vptr[i]) 420 vptr[i] = k 421 vptr[i] = k # this produces "path compression" in the union-find tree. 421 422 return k 422 423 … … 427 428 #print "verts ", verts, "\n" 428 429 429 #initialise vptr and eptr to negative number outside range 430 # vptr represents the union-find tree. 431 # if vptr[i] = EMPTY < 0, the vertex verts[i] has not been visited yet 432 # if vptr[i] = j > 0, then j verts[j] is the parent of verts[i]. 433 # if vptr[i] = n < 0, then verts[i] is a root vertex and 434 # represents a connected component of n vertices. 435 436 #initialise vptr to negative number outside range 430 437 EMPTY = -max(verts)-len(verts) 431 438 vptr = [EMPTY for k in range(len(verts))] … … 437 444 vl = verts.index(bdry[i][0]) 438 445 vr = verts.index(bdry[i][1]) 439 #rvl = vl, rvr = vr440 446 rvl = findroot(vl) 441 447 rvr = findroot(vr) … … 466 472 # end edge loop 467 473 468 # print "vertex component tree: ", [-v for v in vptr if v<0], "\n" 474 if vptr.count(EMPTY): 475 raise FlagError, "We didn't hit all the vertices in the boundary" 476 477 # print "number of vertices in the connected components: ", [-v for v in vptr if v<0], "\n" 469 478 # print "largest component has: ", -min(vptr), " points. \n" 470 479 # discard the edges in the little components … … 475 484 if cutoff > largest_component: 476 485 cutoff = round((1-small)*largest_component) 477 486 487 # littleind has root indices for small components 478 488 littleind = [k for k in range(len(vptr)) if (vptr[k]<0 and vptr[k]>-cutoff)] 479 489 if littleind: 490 # littlecomp has all vptr indices in the small components 480 491 littlecomp = [k for k in range(len(vptr)) if findroot(k) in littleind] 492 # vdiscard has the vertex indices corresponding to vptr indices 481 493 vdiscard = [verts[k] for k in littlecomp] 482 494 newbdry = [e for e in bdry if not((e[0] in vdiscard) and (e[1] in vdiscard))] 483 # remove boundary triangles that touch discarded components 495 496 newverts = self._vertices_from_edges(newbdry) 497 # recompute the boundary triangles 484 498 newbt = [] 485 499 for bt in self.boundarytriangle: … … 487 501 v1 = self.deltri[bt][1] 488 502 v2 = self.deltri[bt][2] 489 if not (v0 in vdiscard or v1 in vdiscard or v2 in vdiscard):503 if (v0 in newverts or v1 in newverts or v2 in newverts): 490 504 newbt.append(bt) 505 491 506 self.boundarytriangle = newbt 507 508 ## # remove boundary triangles that touch discarded components 509 ## newbt = [] 510 ## for bt in self.boundarytriangle: 511 ## v0 = self.deltri[bt][0] 512 ## v1 = self.deltri[bt][1] 513 ## v2 = self.deltri[bt][2] 514 ## if not (v0 in vdiscard or v1 in vdiscard or v2 in vdiscard): 515 ## newbt.append(bt) 516 ## self.boundarytriangle = newbt 492 517 else: 493 518 newbdry = bdry … … 505 530 bdry = self.boundary 506 531 bdrytri = self.boundarytriangle 532 533 # find boundary triangles that have two edges in bdry 534 # v2ind has the place index relative to the triangle deltri[ind] 535 # for the bdry vertex where the two edges meet. 536 507 537 verts = self._vertices_from_edges(bdry) 508 509 # find boundary triangles that have two edges in bdry 538 510 539 b2etri = [] 511 540 for ind in bdrytri: 512 541 bect = 0 542 v2ind = [0,1,2] 513 543 for j in [0,1,2]: 514 544 eda = (self.deltri[ind][(j+1)%3], self.deltri[ind][(j+2)%3]) … … 516 546 if eda in bdry or edb in bdry: 517 547 bect +=1 548 v2ind.remove(j) 518 549 if bect==2: 519 b2etri.append( ind)550 b2etri.append((ind,v2ind[0])) 520 551 521 552 # test the bdrytri triangles for acute angles 522 553 acutetri = [] 523 554 for tind in b2etri: 524 tri = self.deltri[tind] 555 tri = self.deltri[tind[0]] 556 525 557 dx = array([self.points[tri[(i+1)%3],0] - self.points[tri[(i+2)%3],0] for i in [0,1,2]]) 526 558 dy = array([self.points[tri[(i+1)%3],1] - self.points[tri[(i+2)%3],1] for i in [0,1,2]]) 527 559 anglesign = array([(-dx[(i+1)%3]*dx[(i+2)%3]-dy[(i+1)%3]*dy[(i+2)%3]) for i in [0,1,2]]) 528 if product(anglesign) > 0: 529 acutetri.append(tind) 560 # if product(anglesign) > 0: tests if all angles are acute. 561 # but we really want to remove any triangle that has 562 # an acute angle spanned by two edges along the boundary.. 563 # MUST FIX. 564 if anglesign[tind[1]] > 0: 565 acutetri.append(tind[0]) 530 566 531 567 #print "acute boundary triangles ", acutetri -
inundation/ga/storm_surge/alpha_shape/test_alpha_shape.py
r1219 r1364 219 219 assert allclose(answer, result) 220 220 221 def ztest_sharp_indents(self):222 a = [ 1.0, 1.0]223 b = [ 1.0, 5.0]224 c = [4.0, 3.0]225 d = [ 8.0, 5.0]226 e = [ 8.0, 1.0]221 def test_sharp_indents(self): 222 a = [3.0, 1.0] 223 b = [5.0, 3.0] 224 c = [4.0, 4.0] 225 d = [3.0, 5.0] 226 e = [1.0, 3.0] 227 227 228 228 alpha = Alpha_Shape([a,b,c,d,e]) 229 229 alpha.set_boundary_type(raw_boundary=True, 230 230 remove_holes=False, 231 smooth_indents=False, 231 smooth_indents=True, 232 expand_pinch=False) 233 result = alpha.get_boundary() 234 #print "result",result 235 answer = [(3, 4), (2, 3), (0, 1), (1, 2), (4, 0)] 236 assert allclose(answer, result) 237 238 239 def test_small_islands(self): 240 """ 241 I couldn't find a small data set that could test this feature... 242 """ 243 alpha = Alpha_Shape([(191.0,92.0), \ 244 (166.0,79.0), \ 245 (142.0,63.0), \ 246 (124.0,44.0), \ 247 (134.0,19.0), \ 248 (154.0,-9.0), \ 249 (182.0,-21.0), \ 250 (200.0,-3.0), \ 251 (227.0,13.0), \ 252 (237.0,59.0), \ 253 (208.0,86.0), \ 254 (243.0,34.0), \ 255 (209.0,35.0), \ 256 (171.0,15.0), \ 257 (159.0,41.0), \ 258 (189.0,56.0), \ 259 (235.0,-6.0), \ 260 (221.0,-17.0), \ 261 (204.0,-30.0), \ 262 (184.0,-41.0), \ 263 (160.0,-37.0), \ 264 (144.0,-24.0), \ 265 (128.0,-8.0), \ 266 (113.0,9.0), \ 267 (102.0,29.0), \ 268 (95.0,46.0), \ 269 (112.0,61.0), \ 270 (127.0,72.0), \ 271 (145.0,85.0), \ 272 (164.0,100.0), \ 273 (187.0,112.0), \ 274 (210.0,107.0), \ 275 (227.0,94.0), \ 276 (244.0,74.0), \ 277 (261.0,49.0), \ 278 (266.0,20.0), \ 279 (256.0,1.0), \ 280 (244.0,-28.0), \ 281 (228.0,-40.0), \ 282 (208.0,-50.0), \ 283 (192.0,-56.0), \ 284 (169.0,-64.0), \ 285 (155.0,-58.0), \ 286 (141.0,-46.0), \ 287 (129.0,-35.0), \ 288 (113.0,-17.0), \ 289 (103.0,-3.0), \ 290 (90.0,9.0), \ 291 (78.0,26.0), \ 292 (70.0,46.0), \ 293 (86.0,62.0), \ 294 (101.0,79.0), \ 295 (118.0,89.0), \ 296 (135.0,102.0), \ 297 (152.0,115.0), \ 298 (177.0,125.0), \ 299 (192.0,130.0), \ 300 (209.0,125.0), \ 301 (227.0,116.0), \ 302 (234.0,104.0), \ 303 (249.0,92.0), \ 304 (258.0,74.0), \ 305 (264.0,60.0), \ 306 (280.0,42.0), \ 307 (279.0,26.0), \ 308 (274.0,3.0), \ 309 (270.0,-19.0), \ 310 (259.0,-34.0), \ 311 (244.0,-45.0), \ 312 (229.0,-60.0), \ 313 (186.0,4.0), \ 314 (203.0,20.0), \ 315 (215.0,75.0), \ 316 (228.0,49.0), \ 317 (151.0,21.0), \ 318 (163.0,60.0), \ 319 (178.0,70.0), \ 320 (196.0,45.0), \ 321 (261.0,37.0), \ 322 (113.0,36.0), \ 323 (241.0,12.0), \ 324 (196.0,75.0), \ 325 (210.0,58.0), \ 326 (163.248648097,111.633266713), \ 327 (150.844951796,94.8282588211), \ 328 (130.438870783,84.0250394618), \ 329 (112.033385949,72.8217008669), \ 330 (98.8294511765,61.2182430364), \ 331 (184.855086816,-50.4150236771), \ 332 (171.251032808,-52.0155006192), \ 333 (156.446621093,-49.2146659705), \ 334 (148.044117147,-42.8127582019), \ 335 (264.878933922,-2.80083464873), \ 336 (260.077503096,-12.4036963015), \ 337 (255.27607227,-28.0083464873), \ 338 (163.248648097,-22.0065579543), \ 339 (128.03815537,6.00178853298), \ 340 (221.265937249,0.0), \ 341 (202.060213944,-13.6040540081), \ 342 (166.449601981,-2.80083464873), \ 343 (149.244474854,6.80202700405), \ 344 (182.454371403,29.2087041938), \ 345 (147.243878676,30.4090619004), \ 346 (126.437678428,28.0083464873), \ 347 (133.639824668,50.0149044415), \ 348 (176.852702105,43.612996673), \alpha = Alpha_Shape( 349 (123.636843779,-24.4072733675), \ 350 (73.6219393379,35.2104927268), \ 351 (212.463314068,-39.2116850822), \ 352 (254.875953034,66.4197930983), \ 353 (231.669037373,71.6213431603), \ 354 (85.6255164039,23.6070348964), \ 355 (98.4293319409,16.4048886568), \ 356 (223.266533427,-52.0155006192), \ 357 (208.062002477,-57.2170506811), \ 358 (243.672614439,86.425754875), \ 359 (214.06379101,113.633862891), \ 360 (205.261167828,115.234339833), \ 361 (194.057829233,120.836009131), \ 362 (85.2253971684,50.8151429126), \ 363 (84.0250394618,34.0101350202), \ 364 (194.457948469,-47.6141890283), \ 365 (232.86939508,-20.4060810121), \ 366 (252.875356856,23.2069156609), \ 367 (232.86939508,27.2081080162), \ 368 (214.463910245,13.6040540081), \ 369 (182.054252167,18.4054848345), \ 370 (163.248648097,29.6088234294), \ 371 (218.865221836,42.8127582019), \ 372 (176.45258287,97.2289742343), \ 373 (154.04590568,70.8211046892), \ 374 (200.059617766,110.833028242)]) 375 alpha.set_boundary_type(raw_boundary=False , 376 remove_holes=True, 377 smooth_indents=True, 232 378 expand_pinch=True) 233 379 result = alpha.get_boundary() 234 380 #print "result",result 235 answer = [(1, 0), (4, 3), (0, 4), (3, 1)] 236 assert allclose(answer, result) 381 answer = [(45, 106), (44, 106), (46, 47), (44, 43), (48, 111), \ 382 (47, 111), (43, 42), (42, 41), (41, 88), (40, 88), \ 383 (107, 48), (49, 107), (49, 50), (50, 51), (52, 53), \ 384 (51, 52), (53, 54), (55, 83), (54, 83), (40, 114), \ 385 (67, 68), (68, 69), (67, 66), (65, 66), (64, 65), \ 386 (69, 114), (57, 56), (56, 55), (58, 57), (64, 63), \ 387 (61, 62), (62, 63), (59, 60), (58, 59), (60, 61), \ 388 (46, 45)] 389 assert allclose(answer, result) 237 390 #------------------------------------------------------------- 238 391 if __name__ == "__main__":
Note: See TracChangeset
for help on using the changeset viewer.