pyffi.utils.mathutils

1 """A lightweight library for common vector and matrix operations.""" 2 3 # ***** BEGIN LICENSE BLOCK ***** 4 # 5 # Copyright (c) 2007-2011, Python File Format Interface 6 # All rights reserved. 7 # 8 # Redistribution and use in source and binary forms, with or without 9 # modification, are permitted provided that the following conditions 10 # are met: 11 # 12 # * Redistributions of source code must retain the above copyright 13 # notice, this list of conditions and the following disclaimer. 14 # 15 # * Redistributions in binary form must reproduce the above 16 # copyright notice, this list of conditions and the following 17 # disclaimer in the documentation and/or other materials provided 18 # with the distribution. 19 # 20 # * Neither the name of the Python File Format Interface 21 # project nor the names of its contributors may be used to endorse 22 # or promote products derived from this software without specific 23 # prior written permission. 24 # 25 # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 26 # "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 27 # LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 28 # FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE 29 # COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, 30 # INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 31 # BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 32 # LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 33 # CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 34 # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN 35 # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 36 # POSSIBILITY OF SUCH DAMAGE. 37 # 38 # ***** END LICENSE BLOCK ***** 39 40 from itertools import izip 41 import logging 42 import operator 43

44 -def float_to_int(value):

45 """Convert float to integer, rounding and handling nan and inf 46 gracefully. 47 48 >>> float_to_int(0.4) 49 0 50 >>> float_to_int(-0.4) 51 0 52 >>> float_to_int(0.6) 53 1 54 >>> float_to_int(-0.6) 55 -1 56 >>> float_to_int(float('inf')) 57 pyffi.utils.mathutils:WARNING:float_to_int converted +inf to +2147483648. 58 2147483648 59 >>> float_to_int(-float('inf')) 60 pyffi.utils.mathutils:WARNING:float_to_int converted -inf to -2147483648. 61 -2147483648 62 >>> float_to_int(float('nan')) 63 pyffi.utils.mathutils:WARNING:float_to_int converted nan to 0. 64 0 65 """ 66 try: 67 return int(value + 0.5 if value > 0 else value - 0.5) 68 except ValueError: 69 logging.getLogger("pyffi.utils.mathutils").warn( 70 "float_to_int converted nan to 0.") 71 return 0 72 except OverflowError: 73 if value > 0: 74 logging.getLogger("pyffi.utils.mathutils").warn( 75 "float_to_int converted +inf to +2147483648.") 76 return 2147483648 77 else: 78 logging.getLogger("pyffi.utils.mathutils").warn( 79 "float_to_int converted -inf to -2147483648.") 80 return -2147483648

81

82 -def getBoundingBox(veclist):

83 """Calculate bounding box (pair of vectors with minimum and maximum 84 coordinates). 85 86 >>> getBoundingBox([(0,0,0), (1,1,2), (0.5,0.5,0.5)]) 87 ((0, 0, 0), (1, 1, 2))""" 88 if not veclist: 89 # assume 3 dimensions if veclist is empty 90 return (0,0,0), (0,0,0) 91 92 # find bounding box 93 dim = len(veclist[0]) 94 return ( 95 tuple((min(vec[i] for vec in veclist) for i in xrange(dim))), 96 tuple((max(vec[i] for vec in veclist) for i in xrange(dim))))

97

98 -def getCenterRadius(veclist):

99 """Calculate center and radius of given list of vectors. 100 101 >>> getCenterRadius([(0,0,0), (1,1,2), (0.5,0.5,0.5)]) # doctest: +ELLIPSIS 102 ((0.5, 0.5, 1.0), 1.2247...) 103 """ 104 if not veclist: 105 # assume 3 dimensions if veclist is empty 106 return (0,0,0), 0 107 108 # get bounding box 109 vecmin, vecmax = getBoundingBox(veclist) 110 111 # center is in the center of the bounding box 112 center = tuple((minco + maxco) * 0.5 113 for minco, maxco in izip(vecmin, vecmax)) 114 115 # radius is the largest distance from the center 116 r2 = 0.0 117 for vec in veclist: 118 dist = vecSub(center, vec) 119 r2 = max(r2, vecDotProduct(dist, dist)) 120 radius = r2 ** 0.5 121 122 return center, radius

123

124 -def vecSub(vec1, vec2):

125 """Vector substraction.""" 126 return tuple(x - y for x, y in izip(vec1, vec2))

127

128 -def vecAdd(vec1, vec2):

129 return tuple(x + y for x, y in izip(vec1, vec2))

130

131 -def vecscalarMul(vec, scalar):

132 return tuple(x * scalar for x in vec)

133

134 -def vecDotProduct(vec1, vec2):

135 """The vector dot product (any dimension). 136 137 >>> vecDotProduct((1,2,3),(4,-5,6)) 138 12""" 139 return sum(x1 * x2 for x1, x2 in izip(vec1, vec2))

140

141 -def vecDistance(vec1, vec2):

142 """Return distance between two vectors (any dimension). 143 144 >>> vecDistance((1,2,3),(4,-5,6)) # doctest: +ELLIPSIS 145 8.185... 146 """ 147 return vecNorm(vecSub(vec1, vec2))

148

149 -def vecNormal(vec1, vec2, vec3):

150 """Returns a vector that is orthogonal on C{triangle}.""" 151 return vecCrossProduct(vecSub(vec2, vec1), vecSub(vec3, vec1))

152

153 -def vecDistanceAxis(axis, vec):

154 """Return distance between the axis spanned by axis[0] and axis[1] and the 155 vector v, in 3 dimensions. Raises ZeroDivisionError if the axis points 156 coincide. 157 158 >>> vecDistanceAxis([(0,0,0), (0,0,1)], (0,3.5,0)) 159 3.5 160 >>> vecDistanceAxis([(0,0,0), (1,1,1)], (0,1,0.5)) # doctest: +ELLIPSIS 161 0.70710678... 162 """ 163 return vecNorm(vecNormal(axis[0], axis[1], vec)) / vecDistance(*axis)

164

165 -def vecDistanceTriangle(triangle, vert):

166 """Return (signed) distance between the plane spanned by triangle[0], 167 triangle[1], and triange[2], and the vector v, in 3 dimensions. 168 169 >>> vecDistanceTriangle([(0,0,0),(1,0,0),(0,1,0)], (0,0,1)) 170 1.0 171 >>> vecDistanceTriangle([(0,0,0),(0,1,0),(1,0,0)], (0,0,1)) 172 -1.0 173 """ 174 normal = vecNormal(*triangle) 175 return vecDotProduct(normal, vecSub(vert, triangle[0])) \ 176 / vecNorm(normal)

177

178 -def vecNorm(vec):

179 """Norm of a vector (any dimension). 180 181 >>> vecNorm((2,3,4)) # doctest: +ELLIPSIS 182 5.3851648... 183 """ 184 return vecDotProduct(vec, vec) ** 0.5

185

186 -def vecNormalized(vec):

187 """Normalized version of a vector (any dimension). 188 189 >>> vecNormalized((2,3,4)) # doctest: +ELLIPSIS 190 (0.371..., 0.557..., 0.742...) 191 """ 192 return vecscalarMul(vec, 1.0 / vecNorm(vec))

193

194 -def vecCrossProduct(vec1, vec2):

195 """The vector cross product (in 3d). 196 197 >>> vecCrossProduct((1,0,0),(0,1,0)) 198 (0, 0, 1) 199 >>> vecCrossProduct((1,2,3),(4,5,6)) 200 (-3, 6, -3) 201 """ 202 return (vec1[1] * vec2[2] - vec1[2] * vec2[1], 203 vec1[2] * vec2[0] - vec1[0] * vec2[2], 204 vec1[0] * vec2[1] - vec1[1] * vec2[0])

205

206 -def matTransposed(mat):

207 """Return the transposed of a nxn matrix. 208 209 >>> matTransposed(((1, 2), (3, 4))) 210 ((1, 3), (2, 4))""" 211 dim = len(mat) 212 return tuple( tuple( mat[i][j] 213 for i in xrange(dim) ) 214 for j in xrange(dim) )

215

216 -def matscalarMul(mat, scalar):

217 """Return matrix * scalar.""" 218 dim = len(mat) 219 return tuple( tuple( mat[i][j] * scalar 220 for j in xrange(dim) ) 221 for i in xrange(dim) )

222

223 -def matvecMul(mat, vec):

224 """Return matrix * vector.""" 225 dim = len(mat) 226 return tuple( sum( mat[i][j] * vec[j] for j in xrange(dim) ) 227 for i in xrange(dim) )

228

229 -def matMul(mat1, mat2):

230 """Return matrix * matrix.""" 231 dim = len(mat1) 232 return tuple( tuple( sum( mat1[i][k] * mat2[k][j] 233 for k in xrange(dim) ) 234 for j in xrange(dim) ) 235 for i in xrange(dim) )

236

237 -def matAdd(mat1, mat2):

238 """Return matrix + matrix.""" 239 dim = len(mat1) 240 return tuple( tuple( mat1[i][j] + mat2[i][j] 241 for j in xrange(dim) ) 242 for i in xrange(dim) )

243

244 -def matSub(mat1, mat2):

245 """Return matrix - matrix.""" 246 dim = len(mat1) 247 return tuple( tuple( mat1[i][j] - mat2[i][j] 248 for j in xrange(dim) ) 249 for i in xrange(dim) )

250

251 -def matCofactor(mat, i, j):

252 dim = len(mat) 253 return matDeterminant(tuple( tuple( mat[ii][jj] 254 for jj in xrange(dim) 255 if jj != j ) 256 for ii in xrange(dim) 257 if ii != i ))

258

259 -def matDeterminant(mat):

260 """Calculate determinant. 261 262 >>> matDeterminant( ((1,2,3), (4,5,6), (7,8,9)) ) 263 0 264 >>> matDeterminant( ((1,2,4), (3,0,2), (-3,6,2)) ) 265 36 266 """ 267 dim = len(mat) 268 if dim == 0: return 0 269 elif dim == 1: return mat[0][0] 270 elif dim == 2: return mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1] 271 else: 272 return sum( (-1 if i&1 else 1) * mat[i][0] * matCofactor(mat, i, 0) 273 for i in xrange(dim) )

274 275 if __name__ == "__main__": 276 import doctest 277 doctest.testmod() 278

Source Code for Module pyffi.utils.mathutils