# source:branches/numpy/anuga/alpha_shape/alpha_shape.py@6304

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

Initial commit of numpy changes. Still a long way to go.

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