source: anuga_core/source/anuga/alpha_shape/alpha_shape.py @ 4668

Last change on this file since 4668 was 4668, checked in by duncan, 16 years ago

ticket#189 - hacky fix to unproven memory leak
File size: 25.8 KB
Line 
1"""Alpha shape
2Determine the shape of a set of points.
3
4From website by Kaspar Fischer:
5As mentionned in Edelsbrunner's and Muecke's paper, one can
6intuitively think of an alpha-shape as the following:
7
8Imagine a huge mass of ice-cream making up the space and containing
9the points S as ``hard'' chocolate pieces. Using one of these
10sphere-formed ice-cream spoons we carve out all parts of the ice-cream
11block we can reach without bumping into chocolate pieces, even
12carving out holes in the inside (eg. parts not reachable by simply
13moving the spoon from the outside). We will eventually end up with a
14(not necessarily convex) object bounded by caps, arcs and points. If
15we now straighten all ``round'' faces to triangles and line segments,
16we have an intuitive description of what is called the alpha-shape.
17
18Author: Vanessa Robins, ANU
19"""
20
21import exceptions
22from Numeric import array, Float, divide_safe, sqrt, product
23import random
24
25from load_mesh.loadASCII import export_boundary_file
26from anuga.geospatial_data.geospatial_data import Geospatial_data
27
28class AlphaError(exceptions.Exception):pass
29class PointError(AlphaError): pass
30class FlagError(AlphaError): pass
31
32OUTPUT_FILE_TITLE = "# The alpha shape boundary defined by point index pairs of edges"
33INF = pow(10,9)
34EPSILON = 1.0e-12
35
36def alpha_shape_via_files(point_file, boundary_file, alpha= None):
37    """
38    Load a point file and return the alpha shape boundary as a boundary file.
39   
40    Inputs:
41    point_file: File location of the input file, points format (.csv or .pts)
42    boundary_file: File location of the generated output file
43    alpha: The alpha value can be optionally specified.  If it is not specified
44    the optimum alpha value will be used.
45    """
46    geospatial = Geospatial_data(point_file)
47    points = geospatial.get_data_points(absolute=False)
48   
49    AS = Alpha_Shape(points, alpha)
50    AS.write_boundary(boundary_file)
51   
52
53class Alpha_Shape:
54
55    def __init__(self, points, alpha = None):
56        """
57        An Alpha_Shape requires input of a set of points. Other class routines
58        return the alpha shape boundary.
59       
60        Inputs:       
61          points: List of coordinate pairs [[x1, y1],[x2, y2]..]
62          alpha: The alpha value can be optionally specified.  If it is
63          not specified the optimum alpha value will be used.
64        """
65        self._set_points(points)
66        self._alpha_shape_algorithm(alpha)
67
68    def _set_points(self, points):
69        """
70        Create self.points array, do Error checking
71        Inputs:       
72          points: List of coordinate pairs [[x1, y1],[x2, y2]..]
73        """
74        # print "setting points"
75        if len (points) <= 2:
76            raise PointError, "Too few points to find an alpha shape"
77        if len(points)==3:
78            #check not in a straight line
79            # FIXME check points 1,2,3 if straingt, check if points 2,3,4, ect
80            x01 = points[0][0] - points[1][0]
81            y01 = points[0][1] - points[1][1]
82            x12 = points[1][0] - points[2][0]
83            y12 = points[1][1] - points[2][1]
84            crossprod = x01*y12 - x12*y01
85            if crossprod==0:
86                raise PointError, "Three points on a straight line"
87       
88        #Convert input to Numeric arrays
89        self.points = array(points).astype(Float)
90
91   
92    def write_boundary(self,file_name):
93        """
94        Write the boundary to a file
95        """
96        #print " this info will be in the file"
97        export_boundary_file(file_name, self.get_boundary(),
98                             OUTPUT_FILE_TITLE, delimiter = ',')
99   
100    def get_boundary(self):
101        """
102        Return a list of tuples.
103        Each tuple represents a segment in the boundary
104        by the index of its two end points.
105        The list of tuples represents the alpha shape boundary.
106        """
107        return self.boundary
108
109    def set_boundary_type(self,raw_boundary=True,
110                          remove_holes=False,
111                          smooth_indents=False,
112                          expand_pinch=False,
113                          boundary_points_fraction=0.2):
114        """
115        Use the flags to set constraints on the boundary:
116        raw_boundary    Return raw boundary i.e. the regular edges of the
117                        alpha shape.
118        remove_holes    filter to remove small holes
119                        (small is defined by  boundary_points_fraction )
120        smooth_indents  remove sharp triangular indents in boundary
121        expand_pinch    test for pinch-off and correct
122                           i.e. a boundary vertex with more than two edges.
123        """
124
125        if raw_boundary:
126            # reset alpha shape boundary
127            reg_edge = self.get_regular_edges(self.alpha)
128            self.boundary = [self.edge[k] for k in reg_edge]
129            self._init_boundary_triangles()
130        if remove_holes:
131            #remove small holes
132            self.boundary = self._remove_holes(boundary_points_fraction)
133        if smooth_indents:
134            #remove sharp triangular indents
135            self.boundary = self._smooth_indents()
136        if expand_pinch:
137            #deal with pinch-off
138            self.boundary = self._expand_pinch()
139       
140
141    def get_delaunay(self):
142        """
143        """
144        return self.deltri
145
146    def get_optimum_alpha(self):
147        """
148        """
149        return self.optimum_alpha
150
151    def get_alpha(self):
152        """
153        Return current alpha value.
154        """
155        return self.alpha
156
157    def set_alpha(self,alpha):
158        """
159        Set alpha and update alpha-boundary.
160        """
161        self.alpha = alpha
162        reg_edge = self.get_regular_edges(alpha)   
163        self.boundary = [self.edge[k] for k in reg_edge]
164        self._init_boundary_triangles()
165   
166           
167    def _alpha_shape_algorithm(self, alpha=None):
168        """
169        Given a set of points (self.points) and an optional alpha value
170        determines the alpha shape boundary (stored in self.boundary,
171        accessed by get_boundary).
172       
173        Inputs:
174          alpha: The alpha value can be optionally specified.  If it is
175          not specified the optimum alpha value will be used.
176        """
177
178        #print "starting alpha shape algorithm"
179
180        self.alpha = alpha
181
182        ## Build Delaunay triangulation
183        from anuga.mesh_engine.mesh_engine import generate_mesh
184        points = []
185        seglist = []
186        holelist = []
187        regionlist = []
188        pointattlist = []
189        segattlist = []
190
191        points = [(pt[0], pt[1]) for pt in self.points]
192        pointattlist = [ [] for i in range(len(points)) ] 
193        mode = "Qzcn"
194        #print "computing delaunay triangulation ... \n"
195        tridata = generate_mesh(points,seglist,holelist,regionlist,
196                                 pointattlist,segattlist,mode)
197        #print tridata
198        #print "point attlist: ", tridata['generatedpointattributelist'],"\n"
199        #print "hull segments: ", tridata['generatedsegmentlist'], "\n"
200        self.deltri = tridata['generatedtrianglelist']
201        self.deltrinbr = tridata['generatedtriangleneighborlist']
202        self.hulledges = tridata['generatedsegmentlist'] 
203
204        #print "Number of delaunay triangles: ", len(self.deltri), "\n"
205        #print "deltrinbrs: ", self.deltrinbr, "\n"
206
207        ## Build Alpha table
208        ## the following routines determine alpha thresholds for the
209        ## triangles, edges, and vertices of the delaunay triangulation
210        self._tri_circumradius()
211        # print "Largest circumradius ", max(self.triradius)
212        self._edge_intervals()
213        self._vertex_intervals()
214
215        if alpha==None:
216            # Find optimum alpha
217            # Ken Clarkson's hull program uses smallest alpha so that
218            # every vertex is non-singular so...
219            self.optimum_alpha = max([iv[0] for iv in self.vertexinterval \
220                                      if iv!=[] ])
221            # print "optimum alpha ", self.optimum_alpha
222            alpha = self.optimum_alpha
223        self.alpha = alpha
224        reg_edge = self.get_regular_edges(self.alpha)
225        self.boundary = [self.edge[k] for k in reg_edge]
226        #print "alpha boundary edges ", self.boundary
227        self._init_boundary_triangles()
228       
229        del tridata['generatedtriangleneighborlist']
230        del tridata['generatedpointmarkerlist']
231        del tridata['generatedpointlist']
232        #del tridata['generatedpointattributetitlelist']
233        del tridata['generatedsegmentlist']
234        del tridata['generatedsegmentmarkerlist']
235        del tridata['generatedtrianglelist']
236        del tridata['generatedtriangleattributelist']
237        del tridata['generatedpointattributelist']
238        return
239
240    def _tri_circumradius(self):
241        """
242        Compute circumradii of the delaunay triangles
243        """
244       
245        x = self.points[:,0]
246        y = self.points[:,1]
247        ind1 = [self.deltri[j][0] for j in range(len(self.deltri))]
248        ind2 = [self.deltri[j][1] for j in range(len(self.deltri))]
249        ind3 = [self.deltri[j][2] for j in range(len(self.deltri))]
250
251        x1 = array([ x[j] for j in ind1 ])
252        y1 = array([ y[j] for j in ind1 ])
253        x2 = array([ x[j] for j in ind2 ])
254        y2 = array([ y[j] for j in ind2 ])
255        x3 = array([ x[j] for j in ind3 ])
256        y3 = array([ y[j] for j in ind3 ])
257
258        x21 = x2-x1
259        x31 = x3-x1
260        y21 = y2-y1
261        y31 = y3-y1
262
263        dist21 = x21*x21 + y21*y21
264        dist31 = x31*x31 + y31*y31
265
266        denom = x21*y31 - x31*y21
267        #print "denom = ", denom
268
269        # dx/2, dy/2 give circumcenter relative to x1,y1.
270        # dx = (y31*dist21 - y21*dist31)/denom
271        # dy = (x21*dist31 - x31*dist21)/denom
272        # first need to check for near-zero values of denom
273        delta = 0.00000001
274        zeroind = [k for k in range(len(denom)) if \
275                   (denom[k]< EPSILON and  denom[k] > -EPSILON)]
276        # if some denom values are close to zero,
277        # we perturb the associated vertices and recalculate
278        while zeroind!=[]:
279            random.seed()
280            print "Warning: degenerate triangles found in alpha_shape.py, results may be inaccurate."
281            for d in zeroind:
282                x1[d] = x1[d]+delta*(random.random()-0.5)
283                x2[d] = x2[d]+delta*(random.random()-0.5)
284                x3[d] = x3[d]+delta*(random.random()-0.5)
285                y1[d] = y1[d]+delta*(random.random()-0.5)
286                y2[d] = y2[d]+delta*(random.random()-0.5)
287                y3[d] = y3[d]+delta*(random.random()-0.5)
288            x21 = x2-x1
289            x31 = x3-x1
290            y21 = y2-y1
291            y31 = y3-y1
292            dist21 = x21*x21 + y21*y21
293            dist31 = x31*x31 + y31*y31
294            denom = x21*y31 - x31*y21
295            zeroind = [k for k in range(len(denom)) if \
296                       (denom[k]< EPSILON and  denom[k] > -EPSILON)]
297        try:
298            dx = divide_safe(y31*dist21 - y21*dist31,denom)
299            dy = divide_safe(x21*dist31 - x31*dist21,denom)
300        except ZeroDivisionError:
301            raise  AlphaError
302           
303        self.triradius = 0.5*sqrt(dx*dx + dy*dy)
304        #print "triangle radii", self.triradius
305
306    def _edge_intervals(self):
307        """
308         for each edge, find triples
309         (length/2, min_adj_triradius, max_adj_triradius) if unattached
310         (min_adj_triradius, min_adj_triradius, max_adj_triradius) if attached.
311         An edge is attached if it is opposite an obtuse angle
312        """
313       
314        # It should be possible to rewrite this routine in an array-friendly
315        # form like _tri_circumradius()  if we need to speed things up.
316        # Hard to do though.
317
318        edges = []
319        edgenbrs = []
320        edgeinterval = []
321        for t in range(len(self.deltri)):
322            tri = self.deltri[t]
323            trinbr = self.deltrinbr[t]
324            dx = array([self.points[tri[(i+1)%3],0] -
325                        self.points[tri[(i+2)%3],0] for i in [0,1,2]])
326            dy = array([self.points[tri[(i+1)%3],1] -
327                        self.points[tri[(i+2)%3],1] for i in [0,1,2]])
328            elen = sqrt(dx*dx+dy*dy)
329            # really only need sign - not angle value:
330            anglesign = array([(-dx[(i+1)%3]*dx[(i+2)%3]-
331                                dy[(i+1)%3]*dy[(i+2)%3]) for i in [0,1,2]])
332           
333            #print "dx ",dx,"\n"
334            #print "dy ",dy,"\n"
335            #print "edge lengths of triangle ",t,"\t",elen,"\n"
336            #print "angles ",angle,"\n"
337                         
338            for i in [0,1,2]:
339                j = (i+1)%3
340                k = (i+2)%3
341                if trinbr[i]==-1:
342                    edges.append((tri[j], tri[k]))
343                    edgenbrs.append((t, -1))
344                    edgeinterval.append([0.5*elen[i], self.triradius[t], INF])
345                elif (tri[j]<tri[k]):
346                    edges.append((tri[j], tri[k]))
347                    edgenbrs.append((t, trinbr[i]))
348                    edgeinterval.append([0.5*elen[i],\
349                       min(self.triradius[t],self.triradius[trinbr[i]]),\
350                          max(self.triradius[t],self.triradius[trinbr[i]]) ])
351                else:
352                    continue
353                if anglesign[i] < 0: 
354                    edgeinterval[-1][0] = edgeinterval[-1][1]
355                   
356        self.edge = edges
357        self.edgenbr = edgenbrs
358        self.edgeinterval = edgeinterval
359        #print "edges: ",edges, "\n"
360        #print "edge nbrs:", edgenbrs ,"\n"
361        #print "edge intervals: ",edgeinterval , "\n"
362
363    def _vertex_intervals(self):
364        """
365        for each vertex find pairs
366        (min_adj_triradius, max_adj_triradius)
367        """
368        nv = len(self.points)
369        vertexnbrs = [ [] for i in range(nv)]
370        vertexinterval = [ [] for i in range(nv)] 
371        for t in range(len(self.deltri)):
372            for j in self.deltri[t]:
373                vertexnbrs[j].append(t)
374        for h in range(len(self.hulledges)):
375            for j in self.hulledges[h]:
376                vertexnbrs[j].append(-1)
377       
378        for i in range(nv):
379            radii = [ self.triradius[t] for t in vertexnbrs[i] if t>=0 ]
380            try:
381                vertexinterval[i] = [min(radii), max(radii)]
382                if vertexnbrs[i][-1]==-1:
383                    vertexinterval[i][1]=INF
384            except ValueError:
385                raise AlphaError
386
387        self.vertexnbr = vertexnbrs
388        self.vertexinterval = vertexinterval
389        #print "vertex nbrs ", vertexnbrs, "\n"
390        #print "vertex intervals ",vertexinterval, "\n"
391       
392    def get_alpha_triangles(self,alpha):
393        """
394        Given an alpha value,
395        return indices of triangles in the alpha-shape
396        """
397        def tri_rad_lta(k):
398            return self.triradius[k]<=alpha
399
400        return filter(tri_rad_lta, range(len(self.triradius)))
401
402    def get_regular_edges(self,alpha):
403        """
404        Given an alpha value,
405        return the indices of edges on the boundary of the alpha-shape
406        """
407        def reg_edge(k):
408            return self.edgeinterval[k][1]<=alpha and \
409                   self.edgeinterval[k][2]>alpha
410
411        return filter(reg_edge, range(len(self.edgeinterval)))
412
413    def get_exposed_vertices(self,alpha):
414        """
415        Given an alpha value,
416        return the vertices on the boundary of the alpha-shape
417        """
418        def exp_vert(k):
419            return self.vertexinterval[k][0]<=alpha and \
420                   self.vertexinterval[k][1]>alpha
421
422        return filter(exp_vert, range(len(self.vertexinterval)))       
423
424    def _vertices_from_edges(self,elist):
425        """
426        Returns the list of unique vertex labels from edges in elist
427        """
428
429        v1 = [elist[k][0] for k in range(len(elist))]
430        v2 = [elist[k][1] for k in range(len(elist))]
431        v = v1+v2
432        v.sort()
433        vertices = [v[k] for k in range(len(v)) if v[k]!=v[k-1] ]
434        return vertices
435
436    def _init_boundary_triangles(self):
437        """
438        Creates the initial list of triangle indices
439        exterior to and touching the boundary of the alpha shape
440        """
441        def tri_rad_gta(k):
442            return self.triradius[k]>self.alpha
443
444        extrind = filter(tri_rad_gta, range(len(self.triradius)))
445
446        bv = self._vertices_from_edges(self.boundary)
447       
448        btri = []
449        for et in extrind:
450            v0 = self.deltri[et][0]
451            v1 = self.deltri[et][1]
452            v2 = self.deltri[et][2]
453            if v0 in bv or v1 in bv or v2 in bv:
454                btri.append(et)
455
456        self.boundarytriangle = btri
457 
458        #print "exterior triangles: ", extrind
459
460       
461    def _remove_holes(self,small):
462        """
463        Given the edges in self.boundary, finds the largest components.
464        The routine does this by implementing a union-find algorithm.
465        """
466
467        #print "running _remove_holes \n"
468
469        bdry = self.boundary
470       
471        def findroot(i):
472            if vptr[i] < 0:
473                return i
474            k = findroot(vptr[i])
475            vptr[i] = k    # this produces "path compression" in the
476                           # union-find tree.
477            return k
478
479
480
481        # get a list of unique vertex labels:
482        verts = self._vertices_from_edges(bdry)
483        #print "verts ", verts, "\n"
484
485        # vptr represents the union-find tree.
486        # if vptr[i] = EMPTY < 0, the vertex verts[i] has not been visited yet
487        # if vptr[i] = j > 0, then j verts[j] is the parent of verts[i].
488        # if vptr[i] = n < 0, then verts[i] is a root vertex and
489        #                       represents a connected component of n vertices.
490       
491        #initialise vptr to negative number outside range
492        EMPTY = -max(verts)-len(verts)
493        vptr = [EMPTY for k in range(len(verts))]
494        #print "vptr init: ", vptr, "\n"
495
496        #add edges and maintain union tree
497        for i in range(len(bdry)):
498            #print "edge ",i,"\t",bdry[i]
499            vl = verts.index(bdry[i][0])
500            vr = verts.index(bdry[i][1])
501            rvl = findroot(vl)
502            rvr = findroot(vr)
503            #print "roots:  ",rvl, rvr
504            if not(rvl==rvr):
505                if (vptr[vl]==EMPTY):
506                    if (vptr[vr]==EMPTY):
507                        vptr[vl] = -2
508                        vptr[vr] = vl
509                    else:
510                        vptr[vl] = rvr
511                        vptr[rvr] = vptr[rvr]-1
512                else:
513                    if (vptr[vr]==EMPTY):
514                        vptr[vr] = rvl
515                        vptr[rvl] = vptr[rvl]-1
516                    else:
517                        if vptr[rvl] > vptr[rvr]:
518                            vptr[rvr] = vptr[rvr] + vptr[rvl]
519                            vptr[rvl] = rvr
520                            vptr[vl] = rvr
521                        else:
522                            vptr[rvl] = vptr[rvl] + vptr[rvr]
523                            vptr[rvr] = rvl
524                            vptr[vr] = rvl
525                #print "vptr: ", vptr, "\n"
526        # end edge loop
527
528        if vptr.count(EMPTY):
529            raise FlagError, "We didn't hit all the vertices in the boundary"
530       
531        # print "number of vertices in the connected components: ", [-v for v in vptr if v<0], "\n"
532        # print "largest component has: ", -min(vptr), " points. \n"
533        # discard the edges in the little components
534        # (i.e. those components with less than 'small' fraction of bdry points)
535        cutoff = round(small*len(verts))
536        # print "cutoff component size is ", cutoff, "\n"
537        largest_component = -min(vptr)
538        if cutoff > largest_component:
539            cutoff = round((1-small)*largest_component)
540
541        # littleind has root indices for small components
542        littleind = [k for k in range(len(vptr)) if \
543                     (vptr[k]<0 and vptr[k]>-cutoff)] 
544        if littleind:
545            # littlecomp has all vptr indices in the small components
546            littlecomp = [k for k in range(len(vptr)) if \
547                          findroot(k) in littleind]
548            # vdiscard has the vertex indices corresponding to vptr indices 
549            vdiscard = [verts[k] for k in littlecomp] 
550            newbdry = [e for e in bdry if \
551                       not((e[0] in vdiscard) and (e[1] in vdiscard))]
552
553            newverts = self._vertices_from_edges(newbdry)
554            # recompute the boundary triangles
555            newbt = []
556            for bt in self.boundarytriangle:
557                v0 = self.deltri[bt][0]
558                v1 = self.deltri[bt][1]
559                v2 = self.deltri[bt][2]
560                if (v0 in newverts or v1 in newverts or v2 in newverts):
561                    newbt.append(bt)
562
563            self.boundarytriangle = newbt           
564        else:
565            newbdry = bdry
566   
567        return newbdry
568
569
570    def _smooth_indents(self):
571        """
572        Given edges in bdry, test for acute-angle triangular indents
573        and remove them.
574        """
575       
576        #print "running _smooth_indents \n"
577
578        bdry = self.boundary
579        bdrytri = self.boundarytriangle
580       
581        # find boundary triangles that have two edges in bdry
582        # v2ind has the place index relative to the triangle deltri[ind] 
583        # for the bdry vertex where the two edges meet.
584
585        verts = self._vertices_from_edges(bdry)
586
587        b2etri = []
588        for ind in bdrytri:
589            bect = 0
590            v2ind = [0,1,2]
591            for j in [0,1,2]:
592                eda = (self.deltri[ind][(j+1)%3], self.deltri[ind][(j+2)%3])
593                edb = (self.deltri[ind][(j+2)%3], self.deltri[ind][(j+1)%3])
594                if eda in bdry or edb in bdry:
595                    bect +=1
596                    v2ind.remove(j)
597            if bect==2:
598                b2etri.append((ind,v2ind[0]))
599
600        # test the bdrytri triangles for acute angles
601        acutetri = []
602        for tind in b2etri:
603            tri = self.deltri[tind[0]]
604           
605            dx = array([self.points[tri[(i+1)%3],0] - \
606                        self.points[tri[(i+2)%3],0] for i in [0,1,2]])
607            dy = array([self.points[tri[(i+1)%3],1] - \
608                        self.points[tri[(i+2)%3],1] for i in [0,1,2]])
609            anglesign = array([(-dx[(i+1)%3]*dx[(i+2)%3]-\
610                                dy[(i+1)%3]*dy[(i+2)%3]) for i in [0,1,2]])
611            # record any triangle that has an acute angle spanned by
612            #two edges along the boundary..
613            if anglesign[tind[1]] > 0:
614                acutetri.append(tind[0])
615
616        #print "acute boundary triangles ", acutetri
617
618        # adjust the bdry edges and triangles by adding
619        #in the acutetri triangles
620        for pind in acutetri:
621            bdrytri.remove(pind) 
622            tri = self.deltri[pind]
623            for i in [0,1,2]:
624                bdry.append((tri[(i+1)%3], tri[(i+2)%3]))
625
626        newbdry = []
627        for ed in bdry:
628            numed = bdry.count(ed)+bdry.count((ed[1],ed[0]))
629            if numed%2 == 1:
630                newbdry.append(ed)
631       
632        #print "new boundary ", newbdry
633        return newbdry
634
635    def _expand_pinch(self):
636        """
637        Given edges in bdry, test for vertices with more than 2 incident edges.
638        Expand by adding back in associated triangles.
639        """
640        #print "running _expand_pinch \n"
641
642        bdry = self.boundary
643        bdrytri = self.boundarytriangle
644       
645        v1 = [bdry[k][0] for k in range(len(bdry))]
646        v2 = [bdry[k][1] for k in range(len(bdry))]
647        v = v1+v2
648        v.sort()
649        probv = [v[k] for k in range(len(v)) \
650                 if (v[k]!=v[k-1] and v.count(v[k])>2) ]
651        #print "problem vertices: ", probv 
652
653        # find boundary triangles that have at least one vertex in probv
654        probtri = []
655        for ind in bdrytri:
656            v0 = self.deltri[ind][0]
657            v1 = self.deltri[ind][1]
658            v2 = self.deltri[ind][2]
659            if v0 in probv or v1 in probv or v2 in probv:
660                probtri.append(ind)
661
662        #print "problem boundary triangle indices ", probtri
663
664        # "add in" the problem triangles
665        for pind in probtri:
666            bdrytri.remove(pind) 
667            tri = self.deltri[pind]
668            for i in [0,1,2]:
669                bdry.append((tri[(i+1)%3], tri[(i+2)%3]))
670
671        newbdry = []
672        for ed in bdry:
673            numed = bdry.count(ed)+bdry.count((ed[1],ed[0]))
674            if numed%2 == 1:
675                newbdry.append(ed)
676       
677        #print "new boundary ", newbdry
678        return newbdry
679
680
681#-------------------------------------------------------------
682if __name__ == "__main__":
683    """
684    Load in a data point file.
685    Determine the alpha shape boundary
686    Save the boundary to a file.
687
688    usage: alpha_shape.py point_file.csv boundary_file.bnd [alpha]
689   
690    The alpha value is optional.
691    """
692   
693    import os, sys
694    usage = "usage: %s point_file.csv boundary_file.bnd [alpha]"%os.path.basename(sys.argv[0])
695    if len(sys.argv) < 3:
696        print usage
697    else:
698        point_file = sys.argv[1]
699        boundary_file = sys.argv[2]
700        if len(sys.argv) > 4:
701            alpha = sys.argv[3]
702        else:
703            alpha = None
704
705        #print "about to call alpha shape routine \n"
706        alpha_shape_via_files(point_file, boundary_file, alpha)
707       
Note: See TracBrowser for help on using the repository browser.