| 1 | """Alpha shape |
|---|
| 2 | """ |
|---|
| 3 | |
|---|
| 4 | import exceptions |
|---|
| 5 | from Numeric import array, Float, divide_safe, sqrt |
|---|
| 6 | from load_mesh.loadASCII import load_xya_file, export_boundary_file |
|---|
| 7 | |
|---|
| 8 | class PointError(exceptions.Exception): pass |
|---|
| 9 | |
|---|
| 10 | OUTPUT_FILE_TITLE = "# The alpha shape boundary defined by point index pairs of edges" |
|---|
| 11 | INF = pow(10,9) |
|---|
| 12 | |
|---|
| 13 | def alpha_shape_via_files(point_file, boundary_file, alpha= None): |
|---|
| 14 | |
|---|
| 15 | from load_mesh.loadASCII import load_xya_file |
|---|
| 16 | |
|---|
| 17 | point_dict = load_xya_file(point_file) |
|---|
| 18 | points = point_dict['pointlist'] |
|---|
| 19 | #title_string = point_dict['title'] |
|---|
| 20 | |
|---|
| 21 | alpha = Alpha_Shape(points) |
|---|
| 22 | alpha.write_boundary(boundary_file) |
|---|
| 23 | |
|---|
| 24 | class Alpha_Shape: |
|---|
| 25 | |
|---|
| 26 | def __init__(self, points, alpha = None): |
|---|
| 27 | |
|---|
| 28 | """ |
|---|
| 29 | Inputs: |
|---|
| 30 | |
|---|
| 31 | points: List of coordinate pairs [[x1, y1],[x2, y2]..] |
|---|
| 32 | |
|---|
| 33 | alpha: alpha shape parameter |
|---|
| 34 | |
|---|
| 35 | """ |
|---|
| 36 | self._set_points(points) |
|---|
| 37 | self._alpha_shape_algorithm(alpha) |
|---|
| 38 | |
|---|
| 39 | |
|---|
| 40 | def _set_points(self, points): |
|---|
| 41 | """ |
|---|
| 42 | """ |
|---|
| 43 | # print "setting points" |
|---|
| 44 | if len (points) <= 2: |
|---|
| 45 | raise PointError, "Too few points to find an alpha shape" |
|---|
| 46 | |
|---|
| 47 | #Convert input to Numeric arrays |
|---|
| 48 | self.points = array(points).astype(Float) |
|---|
| 49 | |
|---|
| 50 | |
|---|
| 51 | def write_boundary(self,file_name): |
|---|
| 52 | """ |
|---|
| 53 | Write the boundary to a file |
|---|
| 54 | """ |
|---|
| 55 | #print " this info will be in the file" |
|---|
| 56 | export_boundary_file(file_name, self.get_boundary(), |
|---|
| 57 | OUTPUT_FILE_TITLE, delimiter = ',') |
|---|
| 58 | |
|---|
| 59 | def get_boundary(self): |
|---|
| 60 | """ |
|---|
| 61 | """ |
|---|
| 62 | #print "getting boundary" |
|---|
| 63 | return self.boundary |
|---|
| 64 | |
|---|
| 65 | def get_delaunay(self): |
|---|
| 66 | """ |
|---|
| 67 | """ |
|---|
| 68 | return self.deltri |
|---|
| 69 | |
|---|
| 70 | |
|---|
| 71 | def _alpha_shape_algorithm(self,alpha): |
|---|
| 72 | |
|---|
| 73 | # A python style guide suggests using _[function name] to |
|---|
| 74 | # specify internal functions |
|---|
| 75 | # - this is the first time I've used it though - DSG |
|---|
| 76 | |
|---|
| 77 | #print "starting alpha shape algorithm" |
|---|
| 78 | |
|---|
| 79 | # Build Delaunay triangulation |
|---|
| 80 | import triang |
|---|
| 81 | points = [] |
|---|
| 82 | seglist = [] |
|---|
| 83 | holelist = [] |
|---|
| 84 | regionlist = [] |
|---|
| 85 | pointattlist = [] |
|---|
| 86 | segattlist = [] |
|---|
| 87 | trilist = [] |
|---|
| 88 | |
|---|
| 89 | points = [(pt[0], pt[1]) for pt in self.points] |
|---|
| 90 | pointattlist = [ [] for i in range(len(points)) ] |
|---|
| 91 | mode = "Qzcn" |
|---|
| 92 | #print "compute delaunay triangulation" |
|---|
| 93 | tridata = triang.genMesh(points,seglist,holelist,regionlist,pointattlist,segattlist,trilist,mode) |
|---|
| 94 | #print tridata |
|---|
| 95 | #print "point attribute list: ", tridata['generatedpointattributelist'],"\n" |
|---|
| 96 | #print "hull segments: ", tridata['generatedsegmentlist'], "\n" |
|---|
| 97 | self.deltri = tridata['generatedtrianglelist'] |
|---|
| 98 | self.deltrinbr = tridata['generatedtriangleneighborlist'] |
|---|
| 99 | self.hulledges = tridata['generatedsegmentlist'] |
|---|
| 100 | |
|---|
| 101 | #print "deltri: ", self.deltri, "\n" |
|---|
| 102 | #print "deltrinbrs: ", self.deltrinbr, "\n" |
|---|
| 103 | |
|---|
| 104 | # Build Alpha table |
|---|
| 105 | self._tri_circumradius() |
|---|
| 106 | self._edge_intervals() |
|---|
| 107 | self._vertex_intervals() |
|---|
| 108 | |
|---|
| 109 | if alpha==None: |
|---|
| 110 | # Find optimum alpha |
|---|
| 111 | # Ken Clarkson's hull program uses smallest alpha so that |
|---|
| 112 | # every vertex is non-singular so... |
|---|
| 113 | self.optimum_alpha = max([ interval[0] for interval in self.vertexinterval]) |
|---|
| 114 | #print "optimum alpha ", self.optimum_alpha |
|---|
| 115 | |
|---|
| 116 | #alpha_tri = self.get_alpha_triangles(self.optimum_alpha) |
|---|
| 117 | #print "alpha shape triangles ", alpha_tri |
|---|
| 118 | |
|---|
| 119 | reg_edge = self.get_regular_edges(self.optimum_alpha) |
|---|
| 120 | #print "alpha boundary edges", reg_edge |
|---|
| 121 | else: |
|---|
| 122 | reg_edge = self.get_regular_edges(alpha) |
|---|
| 123 | |
|---|
| 124 | self.boundary = [self.edge[k] for k in reg_edge] |
|---|
| 125 | #print "alpha boundary edges ", self.boundary |
|---|
| 126 | |
|---|
| 127 | # this produces a baaad boundary |
|---|
| 128 | #boundary = [] |
|---|
| 129 | #for point_index in range(len(self.points)-1): |
|---|
| 130 | # boundary.append([point_index,point_index +1]) |
|---|
| 131 | #boundary.append([len(self.points)-1,0]) |
|---|
| 132 | #self.boundary = boundary |
|---|
| 133 | return |
|---|
| 134 | |
|---|
| 135 | def _tri_circumradius(self): |
|---|
| 136 | |
|---|
| 137 | # compute circumradii of the delaunay triangles |
|---|
| 138 | x = self.points[:,0] |
|---|
| 139 | y = self.points[:,1] |
|---|
| 140 | ind1 = [self.deltri[j][0] for j in range(len(self.deltri))] |
|---|
| 141 | ind2 = [self.deltri[j][1] for j in range(len(self.deltri))] |
|---|
| 142 | ind3 = [self.deltri[j][2] for j in range(len(self.deltri))] |
|---|
| 143 | |
|---|
| 144 | x1 = array([ x[j] for j in ind1 ]) |
|---|
| 145 | y1 = array([ y[j] for j in ind1 ]) |
|---|
| 146 | x2 = array([ x[j] for j in ind2 ]) |
|---|
| 147 | y2 = array([ y[j] for j in ind2 ]) |
|---|
| 148 | x3 = array([ x[j] for j in ind3 ]) |
|---|
| 149 | y3 = array([ y[j] for j in ind3 ]) |
|---|
| 150 | |
|---|
| 151 | x21 = x2-x1 |
|---|
| 152 | x31 = x3-x1 |
|---|
| 153 | y21 = y2-y1 |
|---|
| 154 | y31 = y3-y1 |
|---|
| 155 | |
|---|
| 156 | dist21 = x21*x21 + y21*y21 |
|---|
| 157 | dist31 = x31*x31 + y31*y31 |
|---|
| 158 | |
|---|
| 159 | denom = x21*y31 - x31*y21 |
|---|
| 160 | #print "denom = ", denom |
|---|
| 161 | |
|---|
| 162 | # dx/2, dy/2 give circumcenter relative to x1,y1. |
|---|
| 163 | #dx = (y31*dist21 - y21*dist31)/denom |
|---|
| 164 | #dy = (x21*dist31 - x31*dist21)/denom |
|---|
| 165 | # from Numeric import divide_safe |
|---|
| 166 | try: |
|---|
| 167 | dx = divide_safe(y31*dist21 - y21*dist31,denom) |
|---|
| 168 | dy = divide_safe(x21*dist31 - x31*dist21,denom) |
|---|
| 169 | except ZeroDivisionError: |
|---|
| 170 | print "Found some degenerate triangles." |
|---|
| 171 | raise |
|---|
| 172 | # really need to take care of cases when denom = 0. |
|---|
| 173 | # eg: something like: |
|---|
| 174 | # print "Handling degenerate triangles." |
|---|
| 175 | # ... add some random displacments to trouble points? |
|---|
| 176 | |
|---|
| 177 | self.triradius = 0.5*sqrt(dx*dx + dy*dy) |
|---|
| 178 | #print "triangle radii", self.triradius |
|---|
| 179 | return |
|---|
| 180 | |
|---|
| 181 | def _edge_intervals(self): |
|---|
| 182 | # for each edge, find triples |
|---|
| 183 | # (length/2, min_adj_triradius, max_adj_triradius) if unattached |
|---|
| 184 | # (min_adj_triradius, min_adj_triradius, max_adj_triradius) if attached. |
|---|
| 185 | |
|---|
| 186 | # It should be possible to rewrite this routine in an array-friendly form |
|---|
| 187 | # like _tri_circumradius() if we need to speed things up. |
|---|
| 188 | |
|---|
| 189 | edges = [] |
|---|
| 190 | edgenbrs = [] |
|---|
| 191 | edgeinterval = [] |
|---|
| 192 | for t in range(len(self.deltri)): |
|---|
| 193 | tri = self.deltri[t] |
|---|
| 194 | trinbr = self.deltrinbr[t] |
|---|
| 195 | dx = array([self.points[tri[(i+1)%3],0] - self.points[tri[(i+2)%3],0] for i in [0,1,2]]) |
|---|
| 196 | dy = array([self.points[tri[(i+1)%3],1] - self.points[tri[(i+2)%3],1] for i in [0,1,2]]) |
|---|
| 197 | elen = sqrt(dx*dx+dy*dy) |
|---|
| 198 | #angle = array([(-dx[(i+1)%3]*dx[(i+2)%3]-dy[(i+1)%3]*dy[(i+2)%3])/(elen[(i+1)%3]*elen[(i+2)%3]) for i in [0,1,2]]) |
|---|
| 199 | #angle = arccos(angle) |
|---|
| 200 | # really only need sign - not angle value: |
|---|
| 201 | anglesign = array([(-dx[(i+1)%3]*dx[(i+2)%3]-dy[(i+1)%3]*dy[(i+2)%3]) for i in [0,1,2]]) |
|---|
| 202 | |
|---|
| 203 | #print "dx ",dx,"\n" |
|---|
| 204 | #print "dy ",dy,"\n" |
|---|
| 205 | #print "edge lengths of triangle ",t,"\t",elen,"\n" |
|---|
| 206 | #print "angles ",angle,"\n" |
|---|
| 207 | |
|---|
| 208 | for i in [0,1,2]: |
|---|
| 209 | j = (i+1)%3 |
|---|
| 210 | k = (i+2)%3 |
|---|
| 211 | if trinbr[i]==-1: |
|---|
| 212 | edges.append((tri[j], tri[k])) |
|---|
| 213 | edgenbrs.append((t, -1)) |
|---|
| 214 | edgeinterval.append([0.5*elen[i], self.triradius[t], INF]) |
|---|
| 215 | elif (tri[j]<tri[k]): |
|---|
| 216 | edges.append((tri[j], tri[k])) |
|---|
| 217 | edgenbrs.append((t, trinbr[i])) |
|---|
| 218 | edgeinterval.append([0.5*elen[i],\ |
|---|
| 219 | min(self.triradius[t],self.triradius[trinbr[i]]),\ |
|---|
| 220 | max(self.triradius[t],self.triradius[trinbr[i]]) ]) |
|---|
| 221 | else: |
|---|
| 222 | continue |
|---|
| 223 | #if angle[i]>pi/2: |
|---|
| 224 | if anglesign[i] < 0: |
|---|
| 225 | edgeinterval[-1][0] = edgeinterval[-1][1] |
|---|
| 226 | |
|---|
| 227 | self.edge = edges |
|---|
| 228 | self.edgenbr = edgenbrs |
|---|
| 229 | self.edgeinterval = edgeinterval |
|---|
| 230 | #print "edges: ",edges, "\n" |
|---|
| 231 | #print "edge nbrs:", edgenbrs ,"\n" |
|---|
| 232 | #print "edge intervals: ",edgeinterval , "\n" |
|---|
| 233 | |
|---|
| 234 | def _vertex_intervals(self): |
|---|
| 235 | # for each vertex find pairs |
|---|
| 236 | # (min_adj_triradius, max_adj_triradius) |
|---|
| 237 | nv = len(self.points) |
|---|
| 238 | vertexnbrs = [ [] for i in range(nv)] |
|---|
| 239 | vertexinterval = [ [] for i in range(nv)] |
|---|
| 240 | for t in range(len(self.deltri)): |
|---|
| 241 | for j in self.deltri[t]: |
|---|
| 242 | vertexnbrs[j].append(t) |
|---|
| 243 | for h in range(len(self.hulledges)): |
|---|
| 244 | for j in self.hulledges[h]: |
|---|
| 245 | vertexnbrs[j].append(-1) |
|---|
| 246 | |
|---|
| 247 | for i in range(nv): |
|---|
| 248 | radii = [ self.triradius[t] for t in vertexnbrs[i] if t>=0 ] |
|---|
| 249 | vertexinterval[i] = [min(radii), max(radii)] |
|---|
| 250 | if vertexnbrs[i][-1]==-1: |
|---|
| 251 | vertexinterval[i][1]=INF |
|---|
| 252 | |
|---|
| 253 | self.vertexnbr = vertexnbrs |
|---|
| 254 | self.vertexinterval = vertexinterval |
|---|
| 255 | #print "vertex nbrs ", vertexnbrs, "\n" |
|---|
| 256 | #print "vertex intervals ",vertexinterval, "\n" |
|---|
| 257 | |
|---|
| 258 | def get_alpha_triangles(self,alpha): |
|---|
| 259 | """ |
|---|
| 260 | Given an alpha value, |
|---|
| 261 | return indices of triangles in the alpha-shape |
|---|
| 262 | """ |
|---|
| 263 | def tri_rad_lta(k): |
|---|
| 264 | return self.triradius[k]<=alpha |
|---|
| 265 | |
|---|
| 266 | return filter(tri_rad_lta, range(len(self.triradius))) |
|---|
| 267 | |
|---|
| 268 | def get_regular_edges(self,alpha): |
|---|
| 269 | """ |
|---|
| 270 | Given and alpha value, |
|---|
| 271 | return the edges on the boundary of the alpha-shape |
|---|
| 272 | """ |
|---|
| 273 | def reg_edge(k): |
|---|
| 274 | return self.edgeinterval[k][1]<=alpha and self.edgeinterval[k][2]>alpha |
|---|
| 275 | |
|---|
| 276 | return filter(reg_edge, range(len(self.edgeinterval))) |
|---|
| 277 | |
|---|
| 278 | |
|---|
| 279 | |
|---|
| 280 | #------------------------------------------------------------- |
|---|
| 281 | if __name__ == "__main__": |
|---|
| 282 | """ |
|---|
| 283 | Load in a data point file. |
|---|
| 284 | Determine the alpha shape boundary |
|---|
| 285 | Save the boundary to a file. |
|---|
| 286 | |
|---|
| 287 | usage: alpha_shape.py point_file.xya boundary_file.bnd [alpha] |
|---|
| 288 | |
|---|
| 289 | The alpha value is optional. |
|---|
| 290 | """ |
|---|
| 291 | |
|---|
| 292 | import os, sys |
|---|
| 293 | usage = "usage: %s point_file.xya boundary_file.bnd [alpha]"%os.path.basename(sys.argv[0]) |
|---|
| 294 | # I made up the .bnd affix. Other ideas welcome. -DSG |
|---|
| 295 | if len(sys.argv) < 3: |
|---|
| 296 | print usage |
|---|
| 297 | else: |
|---|
| 298 | point_file = sys.argv[1] |
|---|
| 299 | boundary_file = sys.argv[2] |
|---|
| 300 | if len(sys.argv) > 4: |
|---|
| 301 | alpha = sys.argv[3] |
|---|
| 302 | else: |
|---|
| 303 | alpha = None |
|---|
| 304 | |
|---|
| 305 | alpha_shape_via_files(point_file, boundary_file, alpha) |
|---|
| 306 | |
|---|
