source: branches/source_numpy_conversion/anuga/alpha_shape/alpha_shape.py @ 6982

Last change on this file since 6982 was 5907, checked in by rwilson, 16 years ago

NumPy? conversion.
File size: 25.7 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
22import numpy
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 = numpy.array(points).astype(numpy.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        return
230
231    def _tri_circumradius(self):
232        """
233        Compute circumradii of the delaunay triangles
234        """
235       
236        x = self.points[:,0]
237        y = self.points[:,1]
238        ind1 = [self.deltri[j][0] for j in range(len(self.deltri))]
239        ind2 = [self.deltri[j][1] for j in range(len(self.deltri))]
240        ind3 = [self.deltri[j][2] for j in range(len(self.deltri))]
241
242        x1 = numpy.array([ x[j] for j in ind1 ])
243        y1 = numpy.array([ y[j] for j in ind1 ])
244        x2 = numpy.array([ x[j] for j in ind2 ])
245        y2 = numpy.array([ y[j] for j in ind2 ])
246        x3 = numpy.array([ x[j] for j in ind3 ])
247        y3 = numpy.array([ y[j] for j in ind3 ])
248
249        x21 = x2-x1
250        x31 = x3-x1
251        y21 = y2-y1
252        y31 = y3-y1
253
254        dist21 = x21*x21 + y21*y21
255        dist31 = x31*x31 + y31*y31
256
257        denom = x21*y31 - x31*y21
258        #print "denom = ", denom
259
260        # dx/2, dy/2 give circumcenter relative to x1,y1.
261        # dx = (y31*dist21 - y21*dist31)/denom
262        # dy = (x21*dist31 - x31*dist21)/denom
263        # first need to check for near-zero values of denom
264        delta = 0.00000001
265        zeroind = [k for k in range(len(denom)) if \
266                   (denom[k]< EPSILON and  denom[k] > -EPSILON)]
267        # if some denom values are close to zero,
268        # we perturb the associated vertices and recalculate
269        while zeroind!=[]:
270            random.seed()
271            print "Warning: degenerate triangles found in alpha_shape.py, results may be inaccurate."
272            for d in zeroind:
273                x1[d] = x1[d]+delta*(random.random()-0.5)
274                x2[d] = x2[d]+delta*(random.random()-0.5)
275                x3[d] = x3[d]+delta*(random.random()-0.5)
276                y1[d] = y1[d]+delta*(random.random()-0.5)
277                y2[d] = y2[d]+delta*(random.random()-0.5)
278                y3[d] = y3[d]+delta*(random.random()-0.5)
279            x21 = x2-x1
280            x31 = x3-x1
281            y21 = y2-y1
282            y31 = y3-y1
283            dist21 = x21*x21 + y21*y21
284            dist31 = x31*x31 + y31*y31
285            denom = x21*y31 - x31*y21
286            zeroind = [k for k in range(len(denom)) if \
287                       (denom[k]< EPSILON and  denom[k] > -EPSILON)]
288##NumPy        try:
289##NumPy            dx = divide_safe(y31*dist21 - y21*dist31,denom)
290##NumPy            dy = divide_safe(x21*dist31 - x31*dist21,denom)
291##NumPy        except ZeroDivisionError:
292##NumPy            raise  AlphaError
293# In NumPy there is no divide_safe(), so above code becomes:
294        dx = numpy.divide(y31*dist21 - y21*dist31,denom)
295        dy = numpy.divide(x21*dist31 - x31*dist21,denom)
296        if not numpy.alltrue(numpy.isfinite(dx)) or \
297           not numpy.alltrue(numpy.isfinite(dy)):
298            raise  AlphaError
299           
300        self.triradius = 0.5*numpy.sqrt(dx*dx + dy*dy)
301        #print "triangle radii", self.triradius
302
303    def _edge_intervals(self):
304        """
305         for each edge, find triples
306         (length/2, min_adj_triradius, max_adj_triradius) if unattached
307         (min_adj_triradius, min_adj_triradius, max_adj_triradius) if attached.
308         An edge is attached if it is opposite an obtuse angle
309        """
310       
311        # It should be possible to rewrite this routine in an array-friendly
312        # form like _tri_circumradius()  if we need to speed things up.
313        # Hard to do though.
314
315        edges = []
316        edgenbrs = []
317        edgeinterval = []
318        for t in range(len(self.deltri)):
319            tri = self.deltri[t]
320            trinbr = self.deltrinbr[t]
321            dx = numpy.array([self.points[tri[(i+1)%3],0] -
322                        self.points[tri[(i+2)%3],0] for i in [0,1,2]])
323            dy = numpy.array([self.points[tri[(i+1)%3],1] -
324                        self.points[tri[(i+2)%3],1] for i in [0,1,2]])
325            elen = numpy.sqrt(dx*dx+dy*dy)
326            # really only need sign - not angle value:
327            anglesign = numpy.array([(-dx[(i+1)%3]*dx[(i+2)%3]-
328                                dy[(i+1)%3]*dy[(i+2)%3]) for i in [0,1,2]])
329           
330            #print "dx ",dx,"\n"
331            #print "dy ",dy,"\n"
332            #print "edge lengths of triangle ",t,"\t",elen,"\n"
333            #print "angles ",angle,"\n"
334                         
335            for i in [0,1,2]:
336                j = (i+1)%3
337                k = (i+2)%3
338                if trinbr[i]==-1:
339                    edges.append((tri[j], tri[k]))
340                    edgenbrs.append((t, -1))
341                    edgeinterval.append([0.5*elen[i], self.triradius[t], INF])
342                elif (tri[j]<tri[k]):
343                    edges.append((tri[j], tri[k]))
344                    edgenbrs.append((t, trinbr[i]))
345                    edgeinterval.append([0.5*elen[i],\
346                       min(self.triradius[t],self.triradius[trinbr[i]]),\
347                          max(self.triradius[t],self.triradius[trinbr[i]]) ])
348                else:
349                    continue
350                if anglesign[i] < 0: 
351                    edgeinterval[-1][0] = edgeinterval[-1][1]
352                   
353        self.edge = edges
354        self.edgenbr = edgenbrs
355        self.edgeinterval = edgeinterval
356        #print "edges: ",edges, "\n"
357        #print "edge nbrs:", edgenbrs ,"\n"
358        #print "edge intervals: ",edgeinterval , "\n"
359
360    def _vertex_intervals(self):
361        """
362        for each vertex find pairs
363        (min_adj_triradius, max_adj_triradius)
364        """
365        nv = len(self.points)
366        vertexnbrs = [ [] for i in range(nv)]
367        vertexinterval = [ [] for i in range(nv)] 
368        for t in range(len(self.deltri)):
369            for j in self.deltri[t]:
370                vertexnbrs[int(j)].append(t)
371        for h in range(len(self.hulledges)):
372            for j in self.hulledges[h]:
373                vertexnbrs[int(j)].append(-1)
374       
375        for i in range(nv):
376            radii = [ self.triradius[t] for t in vertexnbrs[i] if t>=0 ]
377            try:
378                vertexinterval[i] = [min(radii), max(radii)]
379                if vertexnbrs[i][-1]==-1:
380                    vertexinterval[i][1]=INF
381            except ValueError:
382                raise AlphaError
383
384        self.vertexnbr = vertexnbrs
385        self.vertexinterval = vertexinterval
386        #print "vertex nbrs ", vertexnbrs, "\n"
387        #print "vertex intervals ",vertexinterval, "\n"
388       
389    def get_alpha_triangles(self,alpha):
390        """
391        Given an alpha value,
392        return indices of triangles in the alpha-shape
393        """
394        def tri_rad_lta(k):
395            return self.triradius[k]<=alpha
396
397        return filter(tri_rad_lta, range(len(self.triradius)))
398
399    def get_regular_edges(self,alpha):
400        """
401        Given an alpha value,
402        return the indices of edges on the boundary of the alpha-shape
403        """
404        def reg_edge(k):
405            return self.edgeinterval[k][1]<=alpha and \
406                   self.edgeinterval[k][2]>alpha
407
408        return filter(reg_edge, range(len(self.edgeinterval)))
409
410    def get_exposed_vertices(self,alpha):
411        """
412        Given an alpha value,
413        return the vertices on the boundary of the alpha-shape
414        """
415        def exp_vert(k):
416            return self.vertexinterval[k][0]<=alpha and \
417                   self.vertexinterval[k][1]>alpha
418
419        return filter(exp_vert, range(len(self.vertexinterval)))       
420
421    def _vertices_from_edges(self,elist):
422        """
423        Returns the list of unique vertex labels from edges in elist
424        """
425
426        v1 = [elist[k][0] for k in range(len(elist))]
427        v2 = [elist[k][1] for k in range(len(elist))]
428        v = v1+v2
429        v.sort()
430        vertices = [v[k] for k in range(len(v)) if v[k]!=v[k-1] ]
431        return vertices
432
433    def _init_boundary_triangles(self):
434        """
435        Creates the initial list of triangle indices
436        exterior to and touching the boundary of the alpha shape
437        """
438        def tri_rad_gta(k):
439            return self.triradius[k]>self.alpha
440
441        extrind = filter(tri_rad_gta, range(len(self.triradius)))
442
443        bv = self._vertices_from_edges(self.boundary)
444       
445        btri = []
446        for et in extrind:
447            v0 = self.deltri[et][0]
448            v1 = self.deltri[et][1]
449            v2 = self.deltri[et][2]
450            if v0 in bv or v1 in bv or v2 in bv:
451                btri.append(et)
452
453        self.boundarytriangle = btri
454 
455        #print "exterior triangles: ", extrind
456
457       
458    def _remove_holes(self,small):
459        """
460        Given the edges in self.boundary, finds the largest components.
461        The routine does this by implementing a union-find algorithm.
462        """
463
464        #print "running _remove_holes \n"
465
466        bdry = self.boundary
467       
468        def findroot(i):
469            if vptr[i] < 0:
470                return i
471            k = findroot(vptr[i])
472            vptr[i] = k    # this produces "path compression" in the
473                           # union-find tree.
474            return k
475
476
477
478        # get a list of unique vertex labels:
479        verts = self._vertices_from_edges(bdry)
480        #print "verts ", verts, "\n"
481
482        # vptr represents the union-find tree.
483        # if vptr[i] = EMPTY < 0, the vertex verts[i] has not been visited yet
484        # if vptr[i] = j > 0, then j verts[j] is the parent of verts[i].
485        # if vptr[i] = n < 0, then verts[i] is a root vertex and
486        #                       represents a connected component of n vertices.
487       
488        #initialise vptr to negative number outside range
489        EMPTY = -max(verts)-len(verts)
490        vptr = [EMPTY for k in range(len(verts))]
491        #print "vptr init: ", vptr, "\n"
492
493        #add edges and maintain union tree
494        for i in range(len(bdry)):
495            #print "edge ",i,"\t",bdry[i]
496            vl = verts.index(bdry[i][0])
497            vr = verts.index(bdry[i][1])
498            rvl = findroot(vl)
499            rvr = findroot(vr)
500            #print "roots:  ",rvl, rvr
501            if not(rvl==rvr):
502                if (vptr[vl]==EMPTY):
503                    if (vptr[vr]==EMPTY):
504                        vptr[vl] = -2
505                        vptr[vr] = vl
506                    else:
507                        vptr[vl] = rvr
508                        vptr[rvr] = vptr[rvr]-1
509                else:
510                    if (vptr[vr]==EMPTY):
511                        vptr[vr] = rvl
512                        vptr[rvl] = vptr[rvl]-1
513                    else:
514                        if vptr[rvl] > vptr[rvr]:
515                            vptr[rvr] = vptr[rvr] + vptr[rvl]
516                            vptr[rvl] = rvr
517                            vptr[vl] = rvr
518                        else:
519                            vptr[rvl] = vptr[rvl] + vptr[rvr]
520                            vptr[rvr] = rvl
521                            vptr[vr] = rvl
522                #print "vptr: ", vptr, "\n"
523        # end edge loop
524
525        if vptr.count(EMPTY):
526            raise FlagError, "We didn't hit all the vertices in the boundary"
527       
528        # print "number of vertices in the connected components: ", [-v for v in vptr if v<0], "\n"
529        # print "largest component has: ", -min(vptr), " points. \n"
530        # discard the edges in the little components
531        # (i.e. those components with less than 'small' fraction of bdry points)
532        cutoff = round(small*len(verts))
533        # print "cutoff component size is ", cutoff, "\n"
534        largest_component = -min(vptr)
535        if cutoff > largest_component:
536            cutoff = round((1-small)*largest_component)
537
538        # littleind has root indices for small components
539        littleind = [k for k in range(len(vptr)) if \
540                     (vptr[k]<0 and vptr[k]>-cutoff)] 
541        if littleind:
542            # littlecomp has all vptr indices in the small components
543            littlecomp = [k for k in range(len(vptr)) if \
544                          findroot(k) in littleind]
545            # vdiscard has the vertex indices corresponding to vptr indices 
546            vdiscard = [verts[k] for k in littlecomp] 
547            newbdry = [e for e in bdry if \
548                       not((e[0] in vdiscard) and (e[1] in vdiscard))]
549
550            newverts = self._vertices_from_edges(newbdry)
551            # recompute the boundary triangles
552            newbt = []
553            for bt in self.boundarytriangle:
554                v0 = self.deltri[bt][0]
555                v1 = self.deltri[bt][1]
556                v2 = self.deltri[bt][2]
557                if (v0 in newverts or v1 in newverts or v2 in newverts):
558                    newbt.append(bt)
559
560            self.boundarytriangle = newbt           
561        else:
562            newbdry = bdry
563   
564        return newbdry
565
566
567    def _smooth_indents(self):
568        """
569        Given edges in bdry, test for acute-angle triangular indents
570        and remove them.
571        """
572       
573        #print "running _smooth_indents \n"
574
575        bdry = self.boundary
576        bdrytri = self.boundarytriangle
577       
578        # find boundary triangles that have two edges in bdry
579        # v2ind has the place index relative to the triangle deltri[ind] 
580        # for the bdry vertex where the two edges meet.
581
582        verts = self._vertices_from_edges(bdry)
583
584        b2etri = []
585        for ind in bdrytri:
586            bect = 0
587            v2ind = [0,1,2]
588            for j in [0,1,2]:
589                eda = (self.deltri[ind][(j+1)%3], self.deltri[ind][(j+2)%3])
590                edb = (self.deltri[ind][(j+2)%3], self.deltri[ind][(j+1)%3])
591                if eda in bdry or edb in bdry:
592                    bect +=1
593                    v2ind.remove(j)
594            if bect==2:
595                b2etri.append((ind,v2ind[0]))
596
597        # test the bdrytri triangles for acute angles
598        acutetri = []
599        for tind in b2etri:
600            tri = self.deltri[tind[0]]
601           
602            dx = numpy.array([self.points[tri[(i+1)%3],0] - \
603                        self.points[tri[(i+2)%3],0] for i in [0,1,2]])
604            dy = numpy.array([self.points[tri[(i+1)%3],1] - \
605                        self.points[tri[(i+2)%3],1] for i in [0,1,2]])
606            anglesign = numpy.array([(-dx[(i+1)%3]*dx[(i+2)%3]-\
607                                dy[(i+1)%3]*dy[(i+2)%3]) for i in [0,1,2]])
608            # record any triangle that has an acute angle spanned by
609            #two edges along the boundary..
610            if anglesign[tind[1]] > 0:
611                acutetri.append(tind[0])
612
613        #print "acute boundary triangles ", acutetri
614
615        # adjust the bdry edges and triangles by adding
616        #in the acutetri triangles
617        for pind in acutetri:
618            bdrytri.remove(pind) 
619            tri = self.deltri[pind]
620            for i in [0,1,2]:
621                bdry.append((tri[(i+1)%3], tri[(i+2)%3]))
622
623        newbdry = []
624        for ed in bdry:
625            numed = bdry.count(ed)+bdry.count((ed[1],ed[0]))
626            if numed%2 == 1:
627                newbdry.append(ed)
628       
629        #print "new boundary ", newbdry
630        return newbdry
631
632    def _expand_pinch(self):
633        """
634        Given edges in bdry, test for vertices with more than 2 incident edges.
635        Expand by adding back in associated triangles.
636        """
637        #print "running _expand_pinch \n"
638
639        bdry = self.boundary
640        bdrytri = self.boundarytriangle
641       
642        v1 = [bdry[k][0] for k in range(len(bdry))]
643        v2 = [bdry[k][1] for k in range(len(bdry))]
644        v = v1+v2
645        v.sort()
646        probv = [v[k] for k in range(len(v)) \
647                 if (v[k]!=v[k-1] and v.count(v[k])>2) ]
648        #print "problem vertices: ", probv 
649
650        # find boundary triangles that have at least one vertex in probv
651        probtri = []
652        for ind in bdrytri:
653            v0 = self.deltri[ind][0]
654            v1 = self.deltri[ind][1]
655            v2 = self.deltri[ind][2]
656            if v0 in probv or v1 in probv or v2 in probv:
657                probtri.append(ind)
658
659        #print "problem boundary triangle indices ", probtri
660
661        # "add in" the problem triangles
662        for pind in probtri:
663            bdrytri.remove(pind) 
664            tri = self.deltri[pind]
665            for i in [0,1,2]:
666                bdry.append((tri[(i+1)%3], tri[(i+2)%3]))
667
668        newbdry = []
669        for ed in bdry:
670            numed = bdry.count(ed)+bdry.count((ed[1],ed[0]))
671            if numed%2 == 1:
672                newbdry.append(ed)
673       
674        #print "new boundary ", newbdry
675        return newbdry
676
677
678#-------------------------------------------------------------
679if __name__ == "__main__":
680    """
681    Load in a data point file.
682    Determine the alpha shape boundary
683    Save the boundary to a file.
684
685    usage: alpha_shape.py point_file.csv boundary_file.bnd [alpha]
686   
687    The alpha value is optional.
688    """
689   
690    import os, sys
691    usage = "usage: %s point_file.csv boundary_file.bnd [alpha]"%os.path.basename(sys.argv[0])
692    if len(sys.argv) < 3:
693        print usage
694    else:
695        point_file = sys.argv[1]
696        boundary_file = sys.argv[2]
697        if len(sys.argv) > 4:
698            alpha = sys.argv[3]
699        else:
700            alpha = None
701
702        #print "about to call alpha shape routine \n"
703        alpha_shape_via_files(point_file, boundary_file, alpha)
704       
Note: See TracBrowser for help on using the repository browser.