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

Last change on this file since 4663 was 4663, checked in by duncan, 17 years ago

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