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

Last change on this file since 5589 was 4955, checked in by steve, 16 years ago

Fixed problems with entries of integer arrays not being integers, but
arrays. So needed to coerce to int in a couple of places.

Also mkstemp was being imported from the wrong module in test_system_tools

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