Changeset 1364


Ignore:
Timestamp:
May 11, 2005, 11:49:06 AM (20 years ago)
Author:
duncan
Message:

updated tests

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  
    7979        raw_boundary   Return raw boundary i.e. the regular edges
    8080        remove_holes   filter to remove small holes
    81         smooth_indents   remove sharp indents in boundary
     81        smooth_indents   remove sharp triangular indents in boundary
    8282        expand_pinch   plus test for pinch-off
    8383                       i.e. boundary vertex with more than two edges.
     
    407407    def _remove_holes(self,small):
    408408        """
    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.
    410411        """
    411412
     
    418419                return i
    419420            k = findroot(vptr[i])
    420             vptr[i] = k
     421            vptr[i] = k    # this produces "path compression" in the union-find tree.
    421422            return k
    422423
     
    427428        #print "verts ", verts, "\n"
    428429
    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
    430437        EMPTY = -max(verts)-len(verts)
    431438        vptr = [EMPTY for k in range(len(verts))]
     
    437444            vl = verts.index(bdry[i][0])
    438445            vr = verts.index(bdry[i][1])
    439             #rvl = vl, rvr = vr
    440446            rvl = findroot(vl)
    441447            rvr = findroot(vr)
     
    466472        # end edge loop
    467473
    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"
    469478        # print "largest component has: ", -min(vptr), " points. \n"
    470479        # discard the edges in the little components
     
    475484        if cutoff > largest_component:
    476485            cutoff = round((1-small)*largest_component)
    477        
     486
     487        # littleind has root indices for small components
    478488        littleind = [k for k in range(len(vptr)) if (vptr[k]<0 and vptr[k]>-cutoff)]
    479489        if littleind:
     490            # littlecomp has all vptr indices in the small components
    480491            littlecomp = [k for k in range(len(vptr)) if findroot(k) in littleind]
     492            # vdiscard has the vertex indices corresponding to vptr indices 
    481493            vdiscard = [verts[k] for k in littlecomp]
    482494            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
    484498            newbt = []
    485499            for bt in self.boundarytriangle:
     
    487501                v1 = self.deltri[bt][1]
    488502                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):
    490504                    newbt.append(bt)
     505
    491506            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
    492517        else:
    493518            newbdry = bdry
     
    505530        bdry = self.boundary
    506531        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
    507537        verts = self._vertices_from_edges(bdry)
    508        
    509         # find boundary triangles that have two edges in bdry
     538
    510539        b2etri = []
    511540        for ind in bdrytri:
    512541            bect = 0
     542            v2ind = [0,1,2]
    513543            for j in [0,1,2]:
    514544                eda = (self.deltri[ind][(j+1)%3], self.deltri[ind][(j+2)%3])
     
    516546                if eda in bdry or edb in bdry:
    517547                    bect +=1
     548                    v2ind.remove(j)
    518549            if bect==2:
    519                 b2etri.append(ind)
     550                b2etri.append((ind,v2ind[0]))
    520551
    521552        # test the bdrytri triangles for acute angles
    522553        acutetri = []
    523554        for tind in b2etri:
    524             tri = self.deltri[tind]
     555            tri = self.deltri[tind[0]]
     556           
    525557            dx = array([self.points[tri[(i+1)%3],0] - self.points[tri[(i+2)%3],0] for i in [0,1,2]])
    526558            dy = array([self.points[tri[(i+1)%3],1] - self.points[tri[(i+2)%3],1] for i in [0,1,2]])
    527559            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])
    530566
    531567        #print "acute boundary triangles ", acutetri
  • inundation/ga/storm_surge/alpha_shape/test_alpha_shape.py

    r1219 r1364  
    219219        assert allclose(answer, result)
    220220 
    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]
    227227
    228228        alpha = Alpha_Shape([a,b,c,d,e])
    229229        alpha.set_boundary_type(raw_boundary=True,
    230230                          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,
    232378                          expand_pinch=True)
    233379        result = alpha.get_boundary()
    234380        #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)   
    237390#-------------------------------------------------------------
    238391if __name__ == "__main__":
Note: See TracChangeset for help on using the changeset viewer.