pyffi.utils.quickhull

1 """A simple implementation of the quick hull algorithm. 2 3 Usually you should only need the L{qhull3d} function, although the module 4 contains some potentially useful helper functions as well. 5 6 Examples 7 ======== 8 9 Tetrahedron 10 ----------- 11 12 >>> import random 13 >>> tetrahedron = [(0,0,0),(1,0,0),(0,1,0),(0,0,1)] 14 >>> for i in range(200): 15 ... alpha = random.random() 16 ... beta = random.random() 17 ... gamma = 1 - alpha - beta 18 ... if gamma >= 0: 19 ... tetrahedron.append((alpha, beta, gamma)) 20 >>> verts, triangles = qhull3d(tetrahedron) 21 >>> (0,0,0) in verts 22 True 23 >>> (1,0,0) in verts 24 True 25 >>> (0,1,0) in verts 26 True 27 >>> (0,0,1) in verts 28 True 29 >>> len(verts) 30 4 31 >>> len(triangles) 32 4 33 34 A double pyramid polyhedron 35 --------------------------- 36 37 >>> poly = [(2,0,0),(0,2,0),(-2,0,0),(0,-2,0),(0,0,2),(0,0,-2)] 38 >>> vertices, triangles = qhull3d(poly) 39 >>> len(vertices) 40 6 41 >>> len(triangles) 42 8 43 >>> for triangle in triangles: # check orientation relative to origin 44 ... verts = [ vertices[i] for i in triangle ] 45 ... assert(vecDotProduct(vecCrossProduct(*verts[:2]), verts[2]) == 8) 46 47 A pyramid 48 --------- 49 50 >>> verts, triangles = qhull3d([(0,0,0),(1,0,0),(0,1,0),(1,1,0),(0.5,0.5,1)]) 51 >>> (0,0,0) in verts 52 True 53 >>> (1,0,0) in verts 54 True 55 >>> (0,1,0) in verts 56 True 57 >>> (1,1,0) in verts 58 True 59 >>> len(verts) 60 5 61 >>> len(triangles) 62 6 63 64 The unit cube 65 ------------- 66 67 >>> import random 68 >>> cube = [(0,0,0),(0,0,1),(0,1,0),(1,0,0),(0,1,1),(1,0,1),(1,1,0),(1,1,1)] 69 >>> for i in range(200): 70 ... cube.append((random.random(), random.random(), random.random())) 71 >>> verts, triangles = qhull3d(cube) 72 >>> len(triangles) # 6 faces, written as 12 triangles 73 12 74 >>> len(verts) 75 8 76 77 A degenerate shape: the unit square 78 ----------------------------------- 79 80 >>> import random 81 >>> plane = [(0,0,0),(1,0,0),(0,1,0),(1,1,0)] 82 >>> for i in range(200): 83 ... plane.append((random.random(), random.random(), 0)) 84 >>> verts, triangles = qhull3d(plane) 85 >>> len(verts) 86 4 87 >>> len(triangles) 88 2 89 90 A random shape 91 -------------- 92 93 >>> import random 94 >>> shape = [] 95 >>> for i in range(2000): 96 ... vert = (random.random(), random.random(), random.random()) 97 ... shape.append(vert) 98 >>> verts, triangles = qhull3d(shape) 99 100 Precision 101 --------- 102 103 >>> plane = [(0,0,0),(1,0,0),(0,1,0),(1,1,0),(1.001, 0.001, 0)] 104 >>> verts, triangles = qhull3d(plane, precision=0.1) 105 >>> len(verts) 106 4 107 >>> len(triangles) 108 2 109 """ 110 111 # ***** BEGIN LICENSE BLOCK ***** 112 # 113 # Copyright (c) 2007-2011, Python File Format Interface 114 # All rights reserved. 115 # 116 # Redistribution and use in source and binary forms, with or without 117 # modification, are permitted provided that the following conditions 118 # are met: 119 # 120 # * Redistributions of source code must retain the above copyright 121 # notice, this list of conditions and the following disclaimer. 122 # 123 # * Redistributions in binary form must reproduce the above 124 # copyright notice, this list of conditions and the following 125 # disclaimer in the documentation and/or other materials provided 126 # with the distribution. 127 # 128 # * Neither the name of the Python File Format Interface 129 # project nor the names of its contributors may be used to endorse 130 # or promote products derived from this software without specific 131 # prior written permission. 132 # 133 # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 134 # "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 135 # LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 136 # FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE 137 # COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, 138 # INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 139 # BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 140 # LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 141 # CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 142 # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN 143 # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 144 # POSSIBILITY OF SUCH DAMAGE. 145 # 146 # ***** END LICENSE BLOCK ***** 147 148 from pyffi.utils.mathutils import * 149 from itertools import izip 150 import operator 151 152 # adapted from 153 # http://en.literateprograms.org/Quickhull_(Python,_arrays)

154 -def qdome2d(vertices, base, normal, precision = 0.0001):

155 """Build a convex dome from C{vertices} on top of the two C{base} vertices, 156 in the plane with normal C{normal}. This is a helper function for 157 L{qhull2d}, and should usually not be called directly. 158 159 :param vertices: The vertices to construct the dome from. 160 :param base: Two vertices that serve as a base for the dome. 161 :param normal: Orientation of the projection plane used for calculating 162 distances. 163 :param precision: Distance used to decide whether points lie outside of 164 the hull or not. 165 :return: A list of vertices that make up a fan of the dome.""" 166 167 vert0, vert1 = base 168 outer = [ (dist, vert) 169 for dist, vert 170 in izip( ( vecDotProduct(vecCrossProduct(normal, 171 vecSub(vert1, vert0)), 172 vecSub(vert, vert0)) 173 for vert in vertices ), 174 vertices ) 175 if dist > precision ] 176 177 if outer: 178 pivot = max(outer)[1] 179 outer_verts = map(operator.itemgetter(1), outer) 180 return qdome2d(outer_verts, [vert0, pivot], normal, precision) \ 181 + qdome2d(outer_verts, [pivot, vert1], normal, precision)[1:] 182 else: 183 return base

184

185 -def qhull2d(vertices, normal, precision = 0.0001):

186 """Simple implementation of the 2d quickhull algorithm in 3 dimensions for 187 vertices viewed from the direction of C{normal}. 188 Returns a fan of vertices that make up the surface. Called by 189 L{qhull3d} to convexify coplanar vertices. 190 191 >>> import random 192 >>> import math 193 >>> plane = [(0,0,0),(1,0,0),(0,1,0),(1,1,0)] 194 >>> for i in range(200): 195 ... plane.append((random.random(), random.random(), 0)) 196 >>> verts = qhull2d(plane, (0,0,1)) 197 >>> len(verts) 198 4 199 >>> disc = [] 200 >>> for i in range(50): 201 ... theta = (2 * math.pi * i) / 50 202 ... disc.append((0, math.sin(theta), math.cos(theta))) 203 >>> verts = qhull2d(disc, (1,0,0)) 204 >>> len(verts) 205 50 206 >>> for i in range(400): 207 ... disc.append((0, 1.4 * random.random() - 0.7, 1.4 * random.random() - 0.7)) 208 >>> verts = qhull2d(disc, (1,0,0)) 209 >>> len(verts) 210 50 211 >>> dist = 2 * math.pi / 50 212 >>> for i in range(len(verts) - 1): 213 ... assert(abs(vecDistance(verts[i], verts[i+1]) - dist) < 0.001) 214 215 :param vertices: The vertices to construct the hull from. 216 :param normal: Orientation of the projection plane used for calculating 217 distances. 218 :param precision: Distance used to decide whether points lie outside of 219 the hull or not. 220 :return: A list of vertices that make up a fan of extreme points. 221 """ 222 base = basesimplex3d(vertices, precision) 223 if len(base) >= 2: 224 vert0, vert1 = base[:2] 225 return qdome2d(vertices, [vert0, vert1], normal, precision) \ 226 + qdome2d(vertices, [vert1, vert0], normal, precision)[1:-1] 227 else: 228 return base

229

230 -def basesimplex3d(vertices, precision = 0.0001):

231 """Find four extreme points, to be used as a starting base for the 232 quick hull algorithm L{qhull3d}. 233 234 The algorithm tries to find four points that are 235 as far apart as possible, because that speeds up the quick hull 236 algorithm. The vertices are ordered so their signed volume is positive. 237 238 If the volume zero up to C{precision} then only three vertices are 239 returned. If the vertices are colinear up to C{precision} then only two 240 vertices are returned. Finally, if the vertices are equal up to C{precision} 241 then just one vertex is returned. 242 243 >>> import random 244 >>> cube = [(0,0,0),(0,0,1),(0,1,0),(1,0,0),(0,1,1),(1,0,1),(1,1,0),(1,1,1)] 245 >>> for i in range(200): 246 ... cube.append((random.random(), random.random(), random.random())) 247 >>> base = basesimplex3d(cube) 248 >>> len(base) 249 4 250 >>> (0,0,0) in base 251 True 252 >>> (1,1,1) in base 253 True 254 255 :param vertices: The vertices to construct extreme points from. 256 :param precision: Distance used to decide whether points coincide, 257 are colinear, or coplanar. 258 :return: A list of one, two, three, or four vertices, depending on the 259 the configuration of the vertices. 260 """ 261 # sort axes by their extent in vertices 262 extents = sorted(range(3), 263 key=lambda i: 264 max(vert[i] for vert in vertices) 265 - min(vert[i] for vert in vertices)) 266 # extents[0] has the index with largest extent etc. 267 # so let us minimize and maximize vertices with key 268 # (vert[extents[0]], vert[extents[1]], vert[extents[2]]) 269 # which we can write as operator.itemgetter(*extents)(vert) 270 vert0 = min(vertices, key=operator.itemgetter(*extents)) 271 vert1 = max(vertices, key=operator.itemgetter(*extents)) 272 # check if all vertices coincide 273 if vecDistance(vert0, vert1) < precision: 274 return [ vert0 ] 275 # as a third extreme point select that one which maximizes the distance 276 # from the vert0 - vert1 axis 277 vert2 = max(vertices, 278 key=lambda vert: vecDistanceAxis((vert0, vert1), vert)) 279 #check if all vertices are colinear 280 if vecDistanceAxis((vert0, vert1), vert2) < precision: 281 return [ vert0, vert1 ] 282 # as a fourth extreme point select one which maximizes the distance from 283 # the v0, v1, v2 triangle 284 vert3 = max(vertices, 285 key=lambda vert: abs(vecDistanceTriangle((vert0, vert1, vert2), 286 vert))) 287 # ensure positive orientation and check if all vertices are coplanar 288 orientation = vecDistanceTriangle((vert0, vert1, vert2), vert3) 289 if orientation > precision: 290 return [ vert0, vert1, vert2, vert3 ] 291 elif orientation < -precision: 292 return [ vert1, vert0, vert2, vert3 ] 293 else: 294 # coplanar 295 return [ vert0, vert1, vert2 ]

296

297 -def qhull3d(vertices, precision = 0.0001, verbose = False):

298 """Return the triangles making up the convex hull of C{vertices}. 299 Considers distances less than C{precision} to be zero (useful to simplify 300 the hull of a complex mesh, at the expense of exactness of the hull). 301 302 :param vertices: The vertices to find the hull of. 303 :param precision: Distance used to decide whether points lie outside of 304 the hull or not. Larger numbers mean fewer triangles, but some vertices 305 may then end up outside the hull, at a distance of no more than 306 C{precision}. 307 :param verbose: Print information about what the algorithm is doing. Only 308 useful for debugging. 309 :return: A list cointaining the extreme points of C{vertices}, and 310 a list of triangle indices containing the triangles that connect 311 all extreme points. 312 """ 313 # find a simplex to start from 314 hull_vertices = basesimplex3d(vertices, precision) 315 316 # handle degenerate cases 317 if len(hull_vertices) == 3: 318 # coplanar 319 hull_vertices = qhull2d(vertices, vecNormal(*hull_vertices), precision) 320 return hull_vertices, [ (0, i+1, i+2) 321 for i in xrange(len(hull_vertices) - 2) ] 322 elif len(hull_vertices) <= 2: 323 # colinear or singular 324 # no triangles for these cases 325 return hull_vertices, [] 326 327 # construct list of triangles of this simplex 328 hull_triangles = set([ operator.itemgetter(i,j,k)(hull_vertices) 329 for i, j, k in ((1,0,2), (0,1,3), (0,3,2), (3,1,2)) ]) 330 331 if verbose: 332 print("starting set", hull_vertices) 333 334 # construct list of outer vertices for each triangle 335 outer_vertices = {} 336 for triangle in hull_triangles: 337 outer = \ 338 [ (dist, vert) 339 for dist, vert 340 in izip( ( vecDistanceTriangle(triangle, vert) 341 for vert in vertices ), 342 vertices ) 343 if dist > precision ] 344 if outer: 345 outer_vertices[triangle] = outer 346 347 # as long as there are triangles with outer vertices 348 while outer_vertices: 349 # grab a triangle and its outer vertices 350 tmp_iter = outer_vertices.iteritems() 351 triangle, outer = tmp_iter.next() # tmp_iter trick to make 2to3 work 352 # calculate pivot point 353 pivot = max(outer)[1] 354 if verbose: 355 print("pivot", pivot) 356 # add it to the list of extreme vertices 357 hull_vertices.append(pivot) 358 # and update the list of triangles: 359 # 1. calculate visibility of triangles to pivot point 360 visibility = [ vecDistanceTriangle(othertriangle, pivot) > precision 361 for othertriangle in outer_vertices.iterkeys() ] 362 # 2. get list of visible triangles 363 visible_triangles = [ othertriangle 364 for othertriangle, visible 365 in izip(outer_vertices.iterkeys(), visibility) 366 if visible ] 367 # 3. find all edges of visible triangles 368 visible_edges = [] 369 for visible_triangle in visible_triangles: 370 visible_edges += [operator.itemgetter(i,j)(visible_triangle) 371 for i, j in ((0,1),(1,2),(2,0))] 372 if verbose: 373 print("visible edges", visible_edges) 374 # 4. construct horizon: edges that are not shared with another triangle 375 horizon_edges = [ edge for edge in visible_edges 376 if not tuple(reversed(edge)) in visible_edges ] 377 # 5. remove visible triangles from list 378 # this puts a hole inside the triangle list 379 visible_outer = set() 380 for outer_verts in outer_vertices.itervalues(): 381 visible_outer |= set(map(operator.itemgetter(1), outer_verts)) 382 for triangle in visible_triangles: 383 if verbose: 384 print("removing", triangle) 385 hull_triangles.remove(triangle) 386 del outer_vertices[triangle] 387 # 6. close triangle list by adding cone from horizon to pivot 388 # also update the outer triangle list as we go 389 for edge in horizon_edges: 390 newtriangle = edge + ( pivot, ) 391 newouter = \ 392 [ (dist, vert) 393 for dist, vert in izip( ( vecDistanceTriangle(newtriangle, 394 vert) 395 for vert in visible_outer ), 396 visible_outer ) 397 if dist > precision ] 398 hull_triangles.add(newtriangle) 399 if newouter: 400 outer_vertices[newtriangle] = newouter 401 if verbose: 402 print("adding", newtriangle, newouter) 403 404 # no triangle has outer vertices anymore 405 # so the convex hull is complete! 406 # remap the triangles to indices that point into hull_vertices 407 return hull_vertices, [ tuple(hull_vertices.index(vert) 408 for vert in triangle) 409 for triangle in hull_triangles ]

410 411 if __name__ == "__main__": 412 import doctest 413 doctest.testmod() 414

Source Code for Module pyffi.utils.quickhull