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

Last change on this file since 6174 was 6174, checked in by rwilson, 15 years ago

Changed .array(A) to .array(A, num.Int) where appropriate. Helps when going to numpy.

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