| [674] | 1 | """Alpha shape |
|---|
| 2 | """ |
|---|
| 3 | |
|---|
| 4 | import exceptions |
|---|
| [959] | 5 | from Numeric import array, Float, divide_safe, sqrt, product |
|---|
| [853] | 6 | from load_mesh.loadASCII import load_points_file, export_boundary_file |
|---|
| [691] | 7 | import random |
|---|
| [674] | 8 | |
|---|
| 9 | class PointError(exceptions.Exception): pass |
|---|
| [987] | 10 | class FlagError(exceptions.Exception): pass |
|---|
| [674] | 11 | |
|---|
| 12 | OUTPUT_FILE_TITLE = "# The alpha shape boundary defined by point index pairs of edges" |
|---|
| 13 | INF = pow(10,9) |
|---|
| 14 | |
|---|
| [959] | 15 | def alpha_shape_via_files(point_file, boundary_file, alpha= None): |
|---|
| [674] | 16 | |
|---|
| [853] | 17 | point_dict = load_points_file(point_file) |
|---|
| [674] | 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 | def _set_points(self, points): |
|---|
| 40 | """ |
|---|
| 41 | """ |
|---|
| 42 | # print "setting points" |
|---|
| 43 | if len (points) <= 2: |
|---|
| 44 | raise PointError, "Too few points to find an alpha shape" |
|---|
| [691] | 45 | if len(points)==3: |
|---|
| 46 | #check not in a straight line |
|---|
| 47 | x01 = points[0][0] - points[1][0] |
|---|
| 48 | y01 = points[0][1] - points[1][1] |
|---|
| 49 | x12 = points[1][0] - points[2][0] |
|---|
| 50 | y12 = points[1][1] - points[2][1] |
|---|
| 51 | crossprod = x01*y12 - x12*y01 |
|---|
| 52 | if crossprod==0: |
|---|
| 53 | raise PointError, "Three points on a straight line" |
|---|
| [674] | 54 | |
|---|
| 55 | #Convert input to Numeric arrays |
|---|
| 56 | self.points = array(points).astype(Float) |
|---|
| 57 | |
|---|
| 58 | |
|---|
| 59 | def write_boundary(self,file_name): |
|---|
| 60 | """ |
|---|
| 61 | Write the boundary to a file |
|---|
| 62 | """ |
|---|
| 63 | #print " this info will be in the file" |
|---|
| 64 | export_boundary_file(file_name, self.get_boundary(), |
|---|
| 65 | OUTPUT_FILE_TITLE, delimiter = ',') |
|---|
| 66 | |
|---|
| 67 | def get_boundary(self): |
|---|
| 68 | """ |
|---|
| 69 | """ |
|---|
| 70 | return self.boundary |
|---|
| 71 | |
|---|
| [987] | 72 | def set_boundary_type(self,raw_boundary=True, |
|---|
| 73 | remove_holes=False, |
|---|
| 74 | smooth_indents=False, |
|---|
| 75 | expand_pinch=False, |
|---|
| 76 | boundary_points_fraction=0.2): |
|---|
| [956] | 77 | """ |
|---|
| [987] | 78 | Use the flags to set constraints on the boundary |
|---|
| 79 | raw_boundary Return raw boundary i.e. the regular edges |
|---|
| 80 | remove_holes filter to remove small holes |
|---|
| 81 | smooth_indents remove sharp indents in boundary |
|---|
| 82 | expand_pinch plus test for pinch-off |
|---|
| 83 | i.e. boundary vertex with more than two edges. |
|---|
| [956] | 84 | """ |
|---|
| 85 | |
|---|
| [987] | 86 | if raw_boundary: |
|---|
| 87 | reg_edge = self.get_regular_edges(self.alpha) |
|---|
| 88 | self.boundary = [self.edge[k] for k in reg_edge] |
|---|
| 89 | self._init_boundary_triangles() |
|---|
| 90 | if remove_holes: |
|---|
| [956] | 91 | #remove holes |
|---|
| [987] | 92 | self.boundary = self._remove_holes(boundary_points_fraction) |
|---|
| 93 | if smooth_indents: |
|---|
| [956] | 94 | #remove sharp indents |
|---|
| [987] | 95 | self.boundary = self._smooth_indents() |
|---|
| 96 | if expand_pinch: |
|---|
| [956] | 97 | #deal with pinch-off |
|---|
| [987] | 98 | self.boundary = self._expand_pinch() |
|---|
| [956] | 99 | |
|---|
| 100 | |
|---|
| [674] | 101 | def get_delaunay(self): |
|---|
| 102 | """ |
|---|
| 103 | """ |
|---|
| 104 | return self.deltri |
|---|
| [691] | 105 | |
|---|
| 106 | def get_optimum_alpha(self): |
|---|
| 107 | """ |
|---|
| 108 | """ |
|---|
| 109 | return self.optimum_alpha |
|---|
| 110 | |
|---|
| 111 | def get_alpha(self): |
|---|
| 112 | """ |
|---|
| 113 | Return current alpha value |
|---|
| 114 | """ |
|---|
| 115 | return self.alpha |
|---|
| 116 | |
|---|
| 117 | def set_alpha(self,alpha): |
|---|
| 118 | """ |
|---|
| 119 | Set alpha and update alpha-boundary. |
|---|
| 120 | """ |
|---|
| 121 | self.alpha = alpha |
|---|
| 122 | reg_edge = self.get_regular_edges(alpha) |
|---|
| 123 | self.boundary = [self.edge[k] for k in reg_edge] |
|---|
| [987] | 124 | self._init_boundary_triangles() |
|---|
| [674] | 125 | |
|---|
| 126 | |
|---|
| 127 | def _alpha_shape_algorithm(self,alpha): |
|---|
| 128 | |
|---|
| 129 | # A python style guide suggests using _[function name] to |
|---|
| 130 | # specify internal functions |
|---|
| 131 | # - this is the first time I've used it though - DSG |
|---|
| 132 | |
|---|
| 133 | #print "starting alpha shape algorithm" |
|---|
| 134 | |
|---|
| [691] | 135 | self.alpha = alpha |
|---|
| 136 | |
|---|
| [674] | 137 | # Build Delaunay triangulation |
|---|
| 138 | import triang |
|---|
| 139 | points = [] |
|---|
| 140 | seglist = [] |
|---|
| 141 | holelist = [] |
|---|
| 142 | regionlist = [] |
|---|
| 143 | pointattlist = [] |
|---|
| 144 | segattlist = [] |
|---|
| 145 | trilist = [] |
|---|
| 146 | |
|---|
| 147 | points = [(pt[0], pt[1]) for pt in self.points] |
|---|
| 148 | pointattlist = [ [] for i in range(len(points)) ] |
|---|
| 149 | mode = "Qzcn" |
|---|
| [691] | 150 | #print "computing delaunay triangulation ... \n" |
|---|
| [674] | 151 | tridata = triang.genMesh(points,seglist,holelist,regionlist,pointattlist,segattlist,trilist,mode) |
|---|
| 152 | #print tridata |
|---|
| 153 | #print "point attribute list: ", tridata['generatedpointattributelist'],"\n" |
|---|
| 154 | #print "hull segments: ", tridata['generatedsegmentlist'], "\n" |
|---|
| 155 | self.deltri = tridata['generatedtrianglelist'] |
|---|
| 156 | self.deltrinbr = tridata['generatedtriangleneighborlist'] |
|---|
| 157 | self.hulledges = tridata['generatedsegmentlist'] |
|---|
| 158 | |
|---|
| [691] | 159 | #print "Number of delaunay triangles: ", len(self.deltri), "\n" |
|---|
| [674] | 160 | #print "deltrinbrs: ", self.deltrinbr, "\n" |
|---|
| 161 | |
|---|
| 162 | # Build Alpha table |
|---|
| [691] | 163 | # print "Building alpha table ... \n" |
|---|
| [674] | 164 | self._tri_circumradius() |
|---|
| [691] | 165 | # print "Largest circumradius ", max(self.triradius) |
|---|
| [674] | 166 | self._edge_intervals() |
|---|
| 167 | self._vertex_intervals() |
|---|
| 168 | |
|---|
| 169 | if alpha==None: |
|---|
| 170 | # Find optimum alpha |
|---|
| 171 | # Ken Clarkson's hull program uses smallest alpha so that |
|---|
| 172 | # every vertex is non-singular so... |
|---|
| [691] | 173 | self.optimum_alpha = max([ iv[0] for iv in self.vertexinterval if iv!=[] ]) |
|---|
| 174 | # print "optimum alpha ", self.optimum_alpha |
|---|
| 175 | self.alpha = self.optimum_alpha |
|---|
| [674] | 176 | #alpha_tri = self.get_alpha_triangles(self.optimum_alpha) |
|---|
| 177 | #print "alpha shape triangles ", alpha_tri |
|---|
| 178 | |
|---|
| 179 | reg_edge = self.get_regular_edges(self.optimum_alpha) |
|---|
| 180 | #print "alpha boundary edges", reg_edge |
|---|
| 181 | else: |
|---|
| 182 | reg_edge = self.get_regular_edges(alpha) |
|---|
| 183 | |
|---|
| 184 | self.boundary = [self.edge[k] for k in reg_edge] |
|---|
| 185 | #print "alpha boundary edges ", self.boundary |
|---|
| [987] | 186 | self._init_boundary_triangles() |
|---|
| [674] | 187 | |
|---|
| 188 | return |
|---|
| 189 | |
|---|
| 190 | def _tri_circumradius(self): |
|---|
| 191 | |
|---|
| 192 | # compute circumradii of the delaunay triangles |
|---|
| 193 | x = self.points[:,0] |
|---|
| 194 | y = self.points[:,1] |
|---|
| 195 | ind1 = [self.deltri[j][0] for j in range(len(self.deltri))] |
|---|
| 196 | ind2 = [self.deltri[j][1] for j in range(len(self.deltri))] |
|---|
| 197 | ind3 = [self.deltri[j][2] for j in range(len(self.deltri))] |
|---|
| 198 | |
|---|
| 199 | x1 = array([ x[j] for j in ind1 ]) |
|---|
| 200 | y1 = array([ y[j] for j in ind1 ]) |
|---|
| 201 | x2 = array([ x[j] for j in ind2 ]) |
|---|
| 202 | y2 = array([ y[j] for j in ind2 ]) |
|---|
| 203 | x3 = array([ x[j] for j in ind3 ]) |
|---|
| 204 | y3 = array([ y[j] for j in ind3 ]) |
|---|
| 205 | |
|---|
| 206 | x21 = x2-x1 |
|---|
| 207 | x31 = x3-x1 |
|---|
| 208 | y21 = y2-y1 |
|---|
| 209 | y31 = y3-y1 |
|---|
| 210 | |
|---|
| 211 | dist21 = x21*x21 + y21*y21 |
|---|
| 212 | dist31 = x31*x31 + y31*y31 |
|---|
| 213 | |
|---|
| 214 | denom = x21*y31 - x31*y21 |
|---|
| 215 | #print "denom = ", denom |
|---|
| [691] | 216 | delta = 0.00000001 |
|---|
| 217 | random.seed() |
|---|
| 218 | zeroind = [k for k in range(len(denom)) if denom[k]==0 ] |
|---|
| 219 | while zeroind!=[]: |
|---|
| 220 | print "Warning: degenerate triangles found in alpha_shape.py, results may be inaccurate." |
|---|
| 221 | for d in zeroind: |
|---|
| 222 | x1[d] = x1[d]+delta*(random.random()-0.5) |
|---|
| 223 | x2[d] = x2[d]+delta*(random.random()-0.5) |
|---|
| 224 | x3[d] = x3[d]+delta*(random.random()-0.5) |
|---|
| 225 | y1[d] = y1[d]+delta*(random.random()-0.5) |
|---|
| 226 | y2[d] = y2[d]+delta*(random.random()-0.5) |
|---|
| 227 | y3[d] = y3[d]+delta*(random.random()-0.5) |
|---|
| 228 | x21 = x2-x1 |
|---|
| 229 | x31 = x3-x1 |
|---|
| 230 | y21 = y2-y1 |
|---|
| 231 | y31 = y3-y1 |
|---|
| 232 | dist21 = x21*x21 + y21*y21 |
|---|
| 233 | dist31 = x31*x31 + y31*y31 |
|---|
| 234 | denom = x21*y31 - x31*y21 |
|---|
| 235 | zeroind = [k for k in range(len(denom)) if denom[k]==0 ] |
|---|
| 236 | |
|---|
| [674] | 237 | # dx/2, dy/2 give circumcenter relative to x1,y1. |
|---|
| 238 | #dx = (y31*dist21 - y21*dist31)/denom |
|---|
| 239 | #dy = (x21*dist31 - x31*dist21)/denom |
|---|
| 240 | # from Numeric import divide_safe |
|---|
| 241 | try: |
|---|
| 242 | dx = divide_safe(y31*dist21 - y21*dist31,denom) |
|---|
| 243 | dy = divide_safe(x21*dist31 - x31*dist21,denom) |
|---|
| 244 | except ZeroDivisionError: |
|---|
| [691] | 245 | print "There are serious problems with the delaunay triangles" |
|---|
| [674] | 246 | raise |
|---|
| [956] | 247 | |
|---|
| [674] | 248 | |
|---|
| 249 | self.triradius = 0.5*sqrt(dx*dx + dy*dy) |
|---|
| 250 | #print "triangle radii", self.triradius |
|---|
| 251 | return |
|---|
| 252 | |
|---|
| [710] | 253 | |
|---|
| 254 | |
|---|
| [674] | 255 | def _edge_intervals(self): |
|---|
| 256 | # for each edge, find triples |
|---|
| 257 | # (length/2, min_adj_triradius, max_adj_triradius) if unattached |
|---|
| 258 | # (min_adj_triradius, min_adj_triradius, max_adj_triradius) if attached. |
|---|
| 259 | |
|---|
| 260 | # It should be possible to rewrite this routine in an array-friendly form |
|---|
| 261 | # like _tri_circumradius() if we need to speed things up. |
|---|
| 262 | |
|---|
| 263 | edges = [] |
|---|
| 264 | edgenbrs = [] |
|---|
| 265 | edgeinterval = [] |
|---|
| 266 | for t in range(len(self.deltri)): |
|---|
| 267 | tri = self.deltri[t] |
|---|
| 268 | trinbr = self.deltrinbr[t] |
|---|
| 269 | dx = array([self.points[tri[(i+1)%3],0] - self.points[tri[(i+2)%3],0] for i in [0,1,2]]) |
|---|
| 270 | dy = array([self.points[tri[(i+1)%3],1] - self.points[tri[(i+2)%3],1] for i in [0,1,2]]) |
|---|
| 271 | elen = sqrt(dx*dx+dy*dy) |
|---|
| 272 | #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]]) |
|---|
| 273 | #angle = arccos(angle) |
|---|
| 274 | # really only need sign - not angle value: |
|---|
| 275 | anglesign = array([(-dx[(i+1)%3]*dx[(i+2)%3]-dy[(i+1)%3]*dy[(i+2)%3]) for i in [0,1,2]]) |
|---|
| 276 | |
|---|
| 277 | #print "dx ",dx,"\n" |
|---|
| 278 | #print "dy ",dy,"\n" |
|---|
| 279 | #print "edge lengths of triangle ",t,"\t",elen,"\n" |
|---|
| 280 | #print "angles ",angle,"\n" |
|---|
| 281 | |
|---|
| 282 | for i in [0,1,2]: |
|---|
| 283 | j = (i+1)%3 |
|---|
| 284 | k = (i+2)%3 |
|---|
| 285 | if trinbr[i]==-1: |
|---|
| 286 | edges.append((tri[j], tri[k])) |
|---|
| 287 | edgenbrs.append((t, -1)) |
|---|
| 288 | edgeinterval.append([0.5*elen[i], self.triradius[t], INF]) |
|---|
| 289 | elif (tri[j]<tri[k]): |
|---|
| 290 | edges.append((tri[j], tri[k])) |
|---|
| 291 | edgenbrs.append((t, trinbr[i])) |
|---|
| 292 | edgeinterval.append([0.5*elen[i],\ |
|---|
| 293 | min(self.triradius[t],self.triradius[trinbr[i]]),\ |
|---|
| 294 | max(self.triradius[t],self.triradius[trinbr[i]]) ]) |
|---|
| 295 | else: |
|---|
| 296 | continue |
|---|
| 297 | #if angle[i]>pi/2: |
|---|
| 298 | if anglesign[i] < 0: |
|---|
| 299 | edgeinterval[-1][0] = edgeinterval[-1][1] |
|---|
| 300 | |
|---|
| 301 | self.edge = edges |
|---|
| 302 | self.edgenbr = edgenbrs |
|---|
| 303 | self.edgeinterval = edgeinterval |
|---|
| 304 | #print "edges: ",edges, "\n" |
|---|
| 305 | #print "edge nbrs:", edgenbrs ,"\n" |
|---|
| 306 | #print "edge intervals: ",edgeinterval , "\n" |
|---|
| 307 | |
|---|
| 308 | def _vertex_intervals(self): |
|---|
| 309 | # for each vertex find pairs |
|---|
| 310 | # (min_adj_triradius, max_adj_triradius) |
|---|
| 311 | nv = len(self.points) |
|---|
| 312 | vertexnbrs = [ [] for i in range(nv)] |
|---|
| 313 | vertexinterval = [ [] for i in range(nv)] |
|---|
| 314 | for t in range(len(self.deltri)): |
|---|
| 315 | for j in self.deltri[t]: |
|---|
| 316 | vertexnbrs[j].append(t) |
|---|
| 317 | for h in range(len(self.hulledges)): |
|---|
| 318 | for j in self.hulledges[h]: |
|---|
| 319 | vertexnbrs[j].append(-1) |
|---|
| 320 | |
|---|
| 321 | for i in range(nv): |
|---|
| 322 | radii = [ self.triradius[t] for t in vertexnbrs[i] if t>=0 ] |
|---|
| [684] | 323 | try: |
|---|
| 324 | vertexinterval[i] = [min(radii), max(radii)] |
|---|
| 325 | if vertexnbrs[i][-1]==-1: |
|---|
| 326 | vertexinterval[i][1]=INF |
|---|
| 327 | except ValueError: |
|---|
| [691] | 328 | # print "Vertex %i is isolated?"%i |
|---|
| 329 | # print "coords: ",self.points[i] |
|---|
| 330 | # print "Vertex nbrs: ", vertexnbrs[i] |
|---|
| 331 | # print "nbr radii: ",radii |
|---|
| 332 | # vertexinterval[i] = [0,INF] |
|---|
| [684] | 333 | pass |
|---|
| [674] | 334 | |
|---|
| 335 | self.vertexnbr = vertexnbrs |
|---|
| 336 | self.vertexinterval = vertexinterval |
|---|
| 337 | #print "vertex nbrs ", vertexnbrs, "\n" |
|---|
| 338 | #print "vertex intervals ",vertexinterval, "\n" |
|---|
| 339 | |
|---|
| 340 | def get_alpha_triangles(self,alpha): |
|---|
| 341 | """ |
|---|
| 342 | Given an alpha value, |
|---|
| 343 | return indices of triangles in the alpha-shape |
|---|
| 344 | """ |
|---|
| 345 | def tri_rad_lta(k): |
|---|
| 346 | return self.triradius[k]<=alpha |
|---|
| 347 | |
|---|
| 348 | return filter(tri_rad_lta, range(len(self.triradius))) |
|---|
| 349 | |
|---|
| 350 | def get_regular_edges(self,alpha): |
|---|
| 351 | """ |
|---|
| [691] | 352 | Given an alpha value, |
|---|
| [959] | 353 | return the indices of edges on the boundary of the alpha-shape |
|---|
| [674] | 354 | """ |
|---|
| 355 | def reg_edge(k): |
|---|
| 356 | return self.edgeinterval[k][1]<=alpha and self.edgeinterval[k][2]>alpha |
|---|
| 357 | |
|---|
| 358 | return filter(reg_edge, range(len(self.edgeinterval))) |
|---|
| 359 | |
|---|
| [956] | 360 | def get_exposed_vertices(self,alpha): |
|---|
| 361 | """ |
|---|
| 362 | Given an alpha value, |
|---|
| 363 | return the vertices on the boundary of the alpha-shape |
|---|
| 364 | """ |
|---|
| 365 | def exp_vert(k): |
|---|
| 366 | return self.vertexinterval[k][0]<=alpha and self.vertexinterval[k][1]>alpha |
|---|
| 367 | |
|---|
| 368 | return filter(exp_vert, range(len(self.vertexinterval))) |
|---|
| 369 | |
|---|
| [959] | 370 | def _vertices_from_edges(self,elist): |
|---|
| 371 | """ |
|---|
| 372 | Returns the list of unique vertex labels from edges in elist |
|---|
| 373 | """ |
|---|
| 374 | |
|---|
| 375 | v1 = [elist[k][0] for k in range(len(elist))] |
|---|
| 376 | v2 = [elist[k][1] for k in range(len(elist))] |
|---|
| 377 | v = v1+v2 |
|---|
| 378 | v.sort() |
|---|
| 379 | vertices = [v[k] for k in range(len(v)) if v[k]!=v[k-1] ] |
|---|
| 380 | return vertices |
|---|
| 381 | |
|---|
| [987] | 382 | def _init_boundary_triangles(self): |
|---|
| 383 | """ |
|---|
| 384 | Creates the initial list of triangle indices |
|---|
| 385 | exterior to and touching the boundary of the alpha shape |
|---|
| 386 | """ |
|---|
| 387 | def tri_rad_gta(k): |
|---|
| 388 | return self.triradius[k]>self.alpha |
|---|
| 389 | |
|---|
| 390 | extrind = filter(tri_rad_gta, range(len(self.triradius))) |
|---|
| 391 | |
|---|
| 392 | bv = self._vertices_from_edges(self.boundary) |
|---|
| [674] | 393 | |
|---|
| [987] | 394 | btri = [] |
|---|
| 395 | for et in extrind: |
|---|
| 396 | v0 = self.deltri[et][0] |
|---|
| 397 | v1 = self.deltri[et][1] |
|---|
| 398 | v2 = self.deltri[et][2] |
|---|
| 399 | if v0 in bv or v1 in bv or v2 in bv: |
|---|
| 400 | btri.append(et) |
|---|
| 401 | |
|---|
| 402 | self.boundarytriangle = btri |
|---|
| 403 | |
|---|
| 404 | #print "exterior triangles: ", extrind |
|---|
| 405 | |
|---|
| 406 | |
|---|
| 407 | def _remove_holes(self,small): |
|---|
| [956] | 408 | """ |
|---|
| 409 | Given the edges in bdry, find the largest exterior components. |
|---|
| 410 | """ |
|---|
| [970] | 411 | |
|---|
| [987] | 412 | #print "running _remove_holes \n" |
|---|
| 413 | |
|---|
| 414 | bdry = self.boundary |
|---|
| [970] | 415 | |
|---|
| [956] | 416 | def findroot(i): |
|---|
| 417 | if vptr[i] < 0: |
|---|
| 418 | return i |
|---|
| 419 | k = findroot(vptr[i]) |
|---|
| 420 | vptr[i] = k |
|---|
| 421 | return k |
|---|
| 422 | |
|---|
| 423 | |
|---|
| 424 | |
|---|
| 425 | # get a list of unique vertex labels: |
|---|
| [959] | 426 | verts = self._vertices_from_edges(bdry) |
|---|
| [956] | 427 | #print "verts ", verts, "\n" |
|---|
| 428 | |
|---|
| 429 | #initialise vptr and eptr to negative number outside range |
|---|
| 430 | EMPTY = -max(verts) - len(bdry) |
|---|
| 431 | vptr = [EMPTY for k in range(len(verts))] |
|---|
| 432 | eptr = [EMPTY for k in range(len(bdry))] |
|---|
| 433 | #print "vptr init: ", vptr, "\n" |
|---|
| 434 | |
|---|
| 435 | #add edges and maintain union tree |
|---|
| 436 | for i in range(len(bdry)): |
|---|
| 437 | #print "edge ",i,"\t",bdry[i] |
|---|
| 438 | vl = verts.index(bdry[i][0]) |
|---|
| 439 | vr = verts.index(bdry[i][1]) |
|---|
| 440 | #rvl = vl, rvr = vr |
|---|
| 441 | rvl = findroot(vl) |
|---|
| 442 | rvr = findroot(vr) |
|---|
| 443 | #print "roots: ",rvl, rvr |
|---|
| 444 | if not(rvl==rvr): |
|---|
| 445 | if (vptr[vl]==EMPTY): |
|---|
| 446 | if (vptr[vr]==EMPTY): |
|---|
| 447 | vptr[vl] = -2 |
|---|
| 448 | vptr[vr] = vl |
|---|
| 449 | #eprt[i] = vl |
|---|
| 450 | else: |
|---|
| 451 | vptr[vl] = rvr |
|---|
| 452 | vptr[rvr] -= 1 |
|---|
| 453 | else: |
|---|
| 454 | if (vptr[vr]==EMPTY): |
|---|
| 455 | vptr[vr] = rvl |
|---|
| 456 | vptr[rvl] -= 1 |
|---|
| 457 | else: |
|---|
| 458 | if vptr[rvl] > vptr[rvr]: |
|---|
| [970] | 459 | vptr[rvr] += vptr[rvl] |
|---|
| [956] | 460 | vptr[rvl] = rvr |
|---|
| 461 | else: |
|---|
| [970] | 462 | vptr[rvl] += vptr[rvr] |
|---|
| [956] | 463 | vptr[rvr] = rvl |
|---|
| 464 | #print "vptr: ", vptr, "\n" |
|---|
| 465 | |
|---|
| 466 | # end edge loop |
|---|
| 467 | |
|---|
| 468 | #print "vertex component tree: ", vptr, "\n" |
|---|
| 469 | |
|---|
| 470 | # discard the edges in the little components |
|---|
| [959] | 471 | # (i.e. those components with less than 'small' fraction of bdry points) |
|---|
| 472 | cutoff = round(small*len(verts)) |
|---|
| [956] | 473 | littleind = [k for k in range(len(vptr)) if (vptr[k]<0 and vptr[k]>-cutoff)] |
|---|
| 474 | if littleind: |
|---|
| [970] | 475 | littlecomp = [k for k in range(len(vptr)) if findroot(k) in littleind] |
|---|
| [956] | 476 | vdiscard = [verts[k] for k in littlecomp] |
|---|
| 477 | newbdry = [e for e in bdry if not((e[0] in vdiscard) or (e[1] in vdiscard))] |
|---|
| [987] | 478 | # remove boundary triangles that touch discarded components |
|---|
| 479 | newbt = [] |
|---|
| 480 | for bt in self.boundarytriangle: |
|---|
| 481 | v0 = self.deltri[bt][0] |
|---|
| 482 | v1 = self.deltri[bt][1] |
|---|
| 483 | v2 = self.deltri[bt][2] |
|---|
| 484 | if not (v0 in vdiscard or v1 in vdiscard or v2 in vdiscard): |
|---|
| 485 | newbt.append(bt) |
|---|
| 486 | self.boundarytriangle = newbt |
|---|
| [956] | 487 | else: |
|---|
| 488 | newbdry = bdry |
|---|
| [987] | 489 | |
|---|
| [956] | 490 | return newbdry |
|---|
| 491 | |
|---|
| 492 | |
|---|
| [987] | 493 | def _smooth_indents(self): |
|---|
| [956] | 494 | """ |
|---|
| 495 | Given edges in bdry, test for acute-angle indents and remove them. |
|---|
| 496 | """ |
|---|
| [970] | 497 | |
|---|
| [987] | 498 | #print "running _smooth_indents \n" |
|---|
| [959] | 499 | |
|---|
| [987] | 500 | bdry = self.boundary |
|---|
| 501 | bdrytri = self.boundarytriangle |
|---|
| [959] | 502 | verts = self._vertices_from_edges(bdry) |
|---|
| [987] | 503 | |
|---|
| 504 | # find boundary triangles that have two edges in bdry |
|---|
| 505 | b2etri = [] |
|---|
| 506 | for ind in bdrytri: |
|---|
| 507 | bect = 0 |
|---|
| [959] | 508 | for j in [0,1,2]: |
|---|
| [987] | 509 | eda = (self.deltri[ind][(j+1)%3], self.deltri[ind][(j+2)%3]) |
|---|
| 510 | edb = (self.deltri[ind][(j+2)%3], self.deltri[ind][(j+1)%3]) |
|---|
| 511 | if eda in bdry or edb in bdry: |
|---|
| 512 | bect +=1 |
|---|
| 513 | if bect==2: |
|---|
| 514 | b2etri.append(ind) |
|---|
| [959] | 515 | |
|---|
| 516 | # test the bdrytri triangles for acute angles |
|---|
| [987] | 517 | acutetri = [] |
|---|
| 518 | for tind in b2etri: |
|---|
| [959] | 519 | tri = self.deltri[tind] |
|---|
| 520 | dx = array([self.points[tri[(i+1)%3],0] - self.points[tri[(i+2)%3],0] for i in [0,1,2]]) |
|---|
| 521 | dy = array([self.points[tri[(i+1)%3],1] - self.points[tri[(i+2)%3],1] for i in [0,1,2]]) |
|---|
| 522 | anglesign = array([(-dx[(i+1)%3]*dx[(i+2)%3]-dy[(i+1)%3]*dy[(i+2)%3]) for i in [0,1,2]]) |
|---|
| 523 | if product(anglesign) > 0: |
|---|
| [987] | 524 | acutetri.append(tind) |
|---|
| [959] | 525 | |
|---|
| [987] | 526 | #print "acute boundary triangles ", acutetri |
|---|
| [959] | 527 | |
|---|
| [987] | 528 | # adjust the bdry edges and triangles by adding in the acutetri triangles |
|---|
| 529 | for pind in acutetri: |
|---|
| 530 | bdrytri.remove(pind) |
|---|
| [959] | 531 | tri = self.deltri[pind] |
|---|
| 532 | for i in [0,1,2]: |
|---|
| [987] | 533 | bdry.append((tri[(i+1)%3], tri[(i+2)%3])) |
|---|
| 534 | |
|---|
| 535 | self.boundarytriangle = bdrytri |
|---|
| 536 | |
|---|
| 537 | newbdry = [] |
|---|
| 538 | for ed in bdry: |
|---|
| 539 | numed = bdry.count(ed)+bdry.count((ed[1],ed[0])) |
|---|
| 540 | if numed%2 == 1: |
|---|
| 541 | newbdry.append(ed) |
|---|
| [959] | 542 | |
|---|
| [987] | 543 | #print "new boundary ", newbdry |
|---|
| 544 | return newbdry |
|---|
| [956] | 545 | |
|---|
| [987] | 546 | def _expand_pinch(self): |
|---|
| [956] | 547 | """ |
|---|
| [970] | 548 | Given edges in bdry, test for vertices with more than 2 incident edges. |
|---|
| [956] | 549 | Expand by adding back in associated triangles... |
|---|
| 550 | """ |
|---|
| [987] | 551 | #print "running _expand_pinch \n" |
|---|
| [959] | 552 | |
|---|
| [987] | 553 | bdry = self.boundary |
|---|
| 554 | bdrytri = self.boundarytriangle |
|---|
| 555 | |
|---|
| [959] | 556 | v1 = [bdry[k][0] for k in range(len(bdry))] |
|---|
| 557 | v2 = [bdry[k][1] for k in range(len(bdry))] |
|---|
| 558 | v = v1+v2 |
|---|
| 559 | v.sort() |
|---|
| [987] | 560 | probv = [v[k] for k in range(len(v)) if (v[k]!=v[k-1] and v.count(v[k])>2) ] |
|---|
| 561 | #print "problem vertices: ", probv |
|---|
| [959] | 562 | |
|---|
| [987] | 563 | # find boundary triangles that have at least one vertex in probv |
|---|
| [970] | 564 | probtri = [] |
|---|
| [987] | 565 | for ind in bdrytri: |
|---|
| 566 | v0 = self.deltri[ind][0] |
|---|
| 567 | v1 = self.deltri[ind][1] |
|---|
| 568 | v2 = self.deltri[ind][2] |
|---|
| 569 | if v0 in probv or v1 in probv or v2 in probv: |
|---|
| 570 | probtri.append(ind) |
|---|
| [959] | 571 | |
|---|
| [987] | 572 | #print "problem boundary triangle indices ", probtri |
|---|
| [959] | 573 | |
|---|
| [970] | 574 | # "add in" the problem triangles |
|---|
| [987] | 575 | for pind in probtri: |
|---|
| 576 | bdrytri.remove(pind) |
|---|
| 577 | tri = self.deltri[pind] |
|---|
| 578 | for i in [0,1,2]: |
|---|
| 579 | bdry.append((tri[(i+1)%3], tri[(i+2)%3])) |
|---|
| [959] | 580 | |
|---|
| [987] | 581 | self.boundarytriangle = bdrytri |
|---|
| [970] | 582 | |
|---|
| [987] | 583 | newbdry = [] |
|---|
| 584 | for ed in bdry: |
|---|
| 585 | numed = bdry.count(ed)+bdry.count((ed[1],ed[0])) |
|---|
| 586 | if numed%2 == 1: |
|---|
| 587 | newbdry.append(ed) |
|---|
| 588 | |
|---|
| 589 | #print "new boundary ", newbdry |
|---|
| 590 | return newbdry |
|---|
| 591 | |
|---|
| 592 | |
|---|
| [674] | 593 | #------------------------------------------------------------- |
|---|
| 594 | if __name__ == "__main__": |
|---|
| 595 | """ |
|---|
| 596 | Load in a data point file. |
|---|
| 597 | Determine the alpha shape boundary |
|---|
| 598 | Save the boundary to a file. |
|---|
| 599 | |
|---|
| 600 | usage: alpha_shape.py point_file.xya boundary_file.bnd [alpha] |
|---|
| 601 | |
|---|
| 602 | The alpha value is optional. |
|---|
| 603 | """ |
|---|
| 604 | |
|---|
| 605 | import os, sys |
|---|
| 606 | usage = "usage: %s point_file.xya boundary_file.bnd [alpha]"%os.path.basename(sys.argv[0]) |
|---|
| 607 | # I made up the .bnd affix. Other ideas welcome. -DSG |
|---|
| 608 | if len(sys.argv) < 3: |
|---|
| 609 | print usage |
|---|
| 610 | else: |
|---|
| 611 | point_file = sys.argv[1] |
|---|
| 612 | boundary_file = sys.argv[2] |
|---|
| 613 | if len(sys.argv) > 4: |
|---|
| 614 | alpha = sys.argv[3] |
|---|
| 615 | else: |
|---|
| 616 | alpha = None |
|---|
| 617 | |
|---|
| [691] | 618 | #print "about to call alpha shape routine \n" |
|---|
| [674] | 619 | alpha_shape_via_files(point_file, boundary_file, alpha) |
|---|
| 620 | |
|---|
