1 /***************************************************************************** 2 Major portions of this software are copyrighted by the Medical College 3 of Wisconsin, 1994-2000, and are released under the Gnu General Public 4 License, Version 2. See the file README.Copyright for details. 5 ******************************************************************************/ 6 7 #ifndef _MCW_3DVECMAT_ 8 #define _MCW_3DVECMAT_ 9 10 /*------------------------------------------------------------------- 11 06 Feb 2001: modified to remove FLOAT_TYPE, and to create 12 separate float and double data types and macros 13 ---------------------------------------------------------------------*/ 14 15 #include <math.h> 16 17 /*-------------------------------------------------------------------*/ 18 /*----- 3-vector and matrix structures -----*/ 19 /*----- Double precision versions exist at bottom of this file. -----*/ 20 21 /*! 3-vector of integers. */ 22 23 typedef struct { int ijk[3] ; } THD_ivec3 ; 24 25 /*! 3-vector of floats. */ 26 27 typedef struct { float xyz[3] ; } THD_fvec3 ; 28 29 /*! 3x3 matrix of floats. */ 30 31 typedef struct { float mat[3][3] ; } THD_mat33 ; 32 33 /*! 3x3 matrix and a 3-vector (basically an affine transform). */ 34 35 typedef struct { 36 THD_fvec3 vv ; /*!< the vector */ 37 THD_mat33 mm ; /*!< the matrix */ 38 } THD_vecmat ; 39 40 /*-------------------------------------------------------------------*/ 41 /*----- macros that operate on 3 vectors and matrices -----*/ 42 43 /*! Load THD_ivec3 iv with 3 integers. */ 44 45 #define LOAD_IVEC3(iv,i,j,k) ( (iv).ijk[0]=(i), \ 46 (iv).ijk[1]=(j), \ 47 (iv).ijk[2]=(k) ) 48 49 /*! Take 3 integers out of THD_ivec3 iv. */ 50 51 #define UNLOAD_IVEC3(iv,i,j,k) ( (i)=(iv).ijk[0], \ 52 (j)=(iv).ijk[1], \ 53 (k)=(iv).ijk[2] ) 54 55 /*! Load THD_fvec3 fv with 3 floats. */ 56 57 #define LOAD_FVEC3(fv,x,y,z) ( (fv).xyz[0]=(x), \ 58 (fv).xyz[1]=(y), \ 59 (fv).xyz[2]=(z) ) 60 61 /*! Take 3 floats out of THD_fvec3 fv. */ 62 63 #define UNLOAD_FVEC3(fv,x,y,z) ( (x)=(fv).xyz[0], \ 64 (y)=(fv).xyz[1], \ 65 (z)=(fv).xyz[2] ) 66 67 static THD_ivec3 tempA_ivec3 , tempB_ivec3 ; /* temps for macros below */ 68 static THD_fvec3 tempA_fvec3 , tempB_fvec3 ; 69 static THD_mat33 tempA_mat33 , tempB_mat33 ; 70 static float tempRWC ; 71 72 /*! Return a temporary THD_ivec3 from 3 integers. */ 73 74 #define TEMP_IVEC3(i,j,k) ( tempB_ivec3.ijk[0]=(i), \ 75 tempB_ivec3.ijk[1]=(j), \ 76 tempB_ivec3.ijk[2]=(k), tempB_ivec3 ) 77 78 /*! Return a temporary THD_fvec3 from 3 floats. */ 79 80 #define TEMP_FVEC3(x,y,z) ( tempB_fvec3.xyz[0]=(x), \ 81 tempB_fvec3.xyz[1]=(y), \ 82 tempB_fvec3.xyz[2]=(z), tempB_fvec3 ) 83 84 /*! Debug printout of a THD_ivec3. */ 85 86 #define DUMP_IVEC3(str,iv) \ 87 fprintf(stderr,"%s: %d %d %d\n",(str),(iv).ijk[0],(iv).ijk[1],(iv).ijk[2]) 88 89 /*! Debug printout of a THD_fvec3. */ 90 91 #define DUMP_FVEC3(str,fv) \ 92 fprintf(stderr,"%s: %13.6g %13.6g %13.6g\n",(str),(fv).xyz[0],(fv).xyz[1],(fv).xyz[2]) 93 94 /*! Debug printout of a THD_mat33. */ 95 96 #define DUMP_MAT33(str,A) \ 97 fprintf(stderr, \ 98 "%10.10s: [ %13.6g %13.6g %13.6g ]\n" \ 99 " [ %13.6g %13.6g %13.6g ]\n" \ 100 " [ %13.6g %13.6g %13.6g ]\n" , \ 101 str , (A).mat[0][0] , (A).mat[0][1] , (A).mat[0][2] , \ 102 (A).mat[1][0] , (A).mat[1][1] , (A).mat[1][2] , \ 103 (A).mat[2][0] , (A).mat[2][1] , (A).mat[2][2] ) 104 105 /*! Debug printout of a THD_vecmat. */ 106 107 #define DUMP_VECMAT(str,vvmm) ( DUMP_MAT33(str,(vvmm).mm) , DUMP_FVEC3(str,(vvmm).vv) ) 108 109 /*--- macros for operations on floating 3 vectors, 110 with heavy use of the comma operator and structure assignment! ---*/ 111 112 /*! Return negation of THD_fvec3 a. */ 113 114 #define NEGATE_FVEC3(a) ( (a).xyz[0] = -(a).xyz[0] , \ 115 (a).xyz[1] = -(a).xyz[1] , \ 116 (a).xyz[2] = -(a).xyz[2] ) 117 118 /*! Return THD_fvec3 a-b. */ 119 120 #define SUB_FVEC3(a,b) \ 121 ( tempA_fvec3.xyz[0] = (a).xyz[0] - (b).xyz[0] , \ 122 tempA_fvec3.xyz[1] = (a).xyz[1] - (b).xyz[1] , \ 123 tempA_fvec3.xyz[2] = (a).xyz[2] - (b).xyz[2] , tempA_fvec3 ) 124 125 /*! Return THD_fvec3 a+b. */ 126 127 #define ADD_FVEC3(a,b) \ 128 ( tempA_fvec3.xyz[0] = (a).xyz[0] + (b).xyz[0] , \ 129 tempA_fvec3.xyz[1] = (a).xyz[1] + (b).xyz[1] , \ 130 tempA_fvec3.xyz[2] = (a).xyz[2] + (b).xyz[2] , tempA_fvec3 ) 131 132 /*! Return THD_fvec3 elementwise multiplication aj * bj */ 133 134 #define MULT_FVEC3(a,b) \ 135 ( tempA_fvec3.xyz[0] = (a).xyz[0] * (b).xyz[0] , \ 136 tempA_fvec3.xyz[1] = (a).xyz[1] * (b).xyz[1] , \ 137 tempA_fvec3.xyz[2] = (a).xyz[2] * (b).xyz[2] , tempA_fvec3 ) 138 139 /*! Return THD_fvec3 aj* b where b is a scalar */ 140 141 #define SCALE_FVEC3(a,b) \ 142 ( tempA_fvec3.xyz[0] = (a).xyz[0] * (b) , \ 143 tempA_fvec3.xyz[1] = (a).xyz[1] * (b) , \ 144 tempA_fvec3.xyz[2] = (a).xyz[2] * (b) , tempA_fvec3 ) 145 146 147 /*! Return THD_fvec3 a/|a| (unit vector). */ 148 149 #define NORMALIZE_FVEC3(a) \ 150 ( tempRWC = (a).xyz[0] * (a).xyz[0] \ 151 + (a).xyz[1] * (a).xyz[1] \ 152 + (a).xyz[2] * (a).xyz[2] , \ 153 tempRWC = (tempRWC > 0) ? (1.0/sqrt(tempRWC)) : 0 , \ 154 tempA_fvec3.xyz[0] = (a).xyz[0] * tempRWC , \ 155 tempA_fvec3.xyz[1] = (a).xyz[1] * tempRWC , \ 156 tempA_fvec3.xyz[2] = (a).xyz[2] * tempRWC , tempA_fvec3 ) 157 158 /*! Return THD_fvec3 a X b (cross product). */ 159 160 #define CROSS_FVEC3(a,b) \ 161 ( tempA_fvec3.xyz[0] = (a).xyz[1]*(b).xyz[2] - (a).xyz[2]*(b).xyz[1] , \ 162 tempA_fvec3.xyz[1] = (a).xyz[2]*(b).xyz[0] - (a).xyz[0]*(b).xyz[2] , \ 163 tempA_fvec3.xyz[2] = (a).xyz[0]*(b).xyz[1] - (a).xyz[1]*(b).xyz[0] , \ 164 tempA_fvec3 ) 165 166 /*! Return float L2norm(a) from a THD_fvec3. */ 167 168 #define SIZE_FVEC3(a) \ 169 sqrt((a).xyz[0]*(a).xyz[0]+(a).xyz[1]*(a).xyz[1]+(a).xyz[2]*(a).xyz[2]) 170 171 /*! Return float a dot b from THD_fvec3 inputs. */ 172 173 #define DOT_FVEC3(a,b) \ 174 ((a).xyz[0]*(b).xyz[0] + (a).xyz[1]*(b).xyz[1] + (a).xyz[2]*(b).xyz[2]) 175 176 /* scale and add two vectors: fa * a + fb * b */ 177 178 #define SCLADD_FVEC3(fa,a,fb,b) \ 179 ( tempA_fvec3.xyz[0] = (fa)*(a).xyz[0] + (fb)*(b).xyz[0] , \ 180 tempA_fvec3.xyz[1] = (fa)*(a).xyz[1] + (fb)*(b).xyz[1] , \ 181 tempA_fvec3.xyz[2] = (fa)*(a).xyz[2] + (fb)*(b).xyz[2] , tempA_fvec3 ) 182 183 /* round to an integer vector (assume non-negative) */ 184 185 #define INT_FVEC3(a) \ 186 ( tempA_ivec3.ijk[0] = (a).xyz[0] + 0.5 , \ 187 tempA_ivec3.ijk[1] = (a).xyz[1] + 0.5 , \ 188 tempA_ivec3.ijk[1] = (a).xyz[2] + 0.5 , tempA_ivec3 ) ; 189 190 #define SIZE_MAT(A) \ 191 ( fabs((A).mat[0][0]) + fabs((A).mat[0][1]) + fabs((A).mat[0][2]) \ 192 +fabs((A).mat[1][0]) + fabs((A).mat[1][1]) + fabs((A).mat[1][2]) \ 193 +fabs((A).mat[2][0]) + fabs((A).mat[2][1]) + fabs((A).mat[2][2]) ) 194 195 /* matrix-vector multiply */ 196 197 #define MATVEC(A,x) \ 198 ( tempA_fvec3.xyz[0] = (A).mat[0][0] * (x).xyz[0] \ 199 +(A).mat[0][1] * (x).xyz[1] \ 200 +(A).mat[0][2] * (x).xyz[2] ,\ 201 tempA_fvec3.xyz[1] = (A).mat[1][0] * (x).xyz[0] \ 202 +(A).mat[1][1] * (x).xyz[1] \ 203 +(A).mat[1][2] * (x).xyz[2] ,\ 204 tempA_fvec3.xyz[2] = (A).mat[2][0] * (x).xyz[0] \ 205 +(A).mat[2][1] * (x).xyz[1] \ 206 +(A).mat[2][2] * (x).xyz[2] , tempA_fvec3 ) 207 208 /* matrix-vector multiply with subtract: output = A (x-b) */ 209 210 #define VECSUB_MAT(A,x,b) \ 211 ( tempA_fvec3.xyz[0] = (A).mat[0][0] * ((x).xyz[0]-(b).xyz[0]) \ 212 +(A).mat[0][1] * ((x).xyz[1]-(b).xyz[1]) \ 213 +(A).mat[0][2] * ((x).xyz[2]-(b).xyz[2]) ,\ 214 tempA_fvec3.xyz[1] = (A).mat[1][0] * ((x).xyz[0]-(b).xyz[0]) \ 215 +(A).mat[1][1] * ((x).xyz[1]-(b).xyz[1]) \ 216 +(A).mat[1][2] * ((x).xyz[2]-(b).xyz[2]) ,\ 217 tempA_fvec3.xyz[2] = (A).mat[2][0] * ((x).xyz[0]-(b).xyz[0]) \ 218 +(A).mat[2][1] * ((x).xyz[1]-(b).xyz[1]) \ 219 +(A).mat[2][2] * ((x).xyz[2]-(b).xyz[2]) ,\ 220 tempA_fvec3 ) 221 222 /* matrix vector multiply with subtract after: A x - b */ 223 224 #define MATVEC_SUB(A,x,b) \ 225 ( tempA_fvec3.xyz[0] = (A).mat[0][0] * (x).xyz[0] \ 226 +(A).mat[0][1] * (x).xyz[1] \ 227 +(A).mat[0][2] * (x).xyz[2] - (b).xyz[0] ,\ 228 tempA_fvec3.xyz[1] = (A).mat[1][0] * (x).xyz[0] \ 229 +(A).mat[1][1] * (x).xyz[1] \ 230 +(A).mat[1][2] * (x).xyz[2] - (b).xyz[1] ,\ 231 tempA_fvec3.xyz[2] = (A).mat[2][0] * (x).xyz[0] \ 232 +(A).mat[2][1] * (x).xyz[1] \ 233 +(A).mat[2][2] * (x).xyz[2] - (b).xyz[2] ,\ 234 tempA_fvec3 ) 235 236 /* matrix vector multiply with add after: A x + b */ 237 238 #define MATVEC_ADD(A,x,b) \ 239 ( tempA_fvec3.xyz[0] = (A).mat[0][0] * (x).xyz[0] \ 240 +(A).mat[0][1] * (x).xyz[1] \ 241 +(A).mat[0][2] * (x).xyz[2] + (b).xyz[0] ,\ 242 tempA_fvec3.xyz[1] = (A).mat[1][0] * (x).xyz[0] \ 243 +(A).mat[1][1] * (x).xyz[1] \ 244 +(A).mat[1][2] * (x).xyz[2] + (b).xyz[1] ,\ 245 tempA_fvec3.xyz[2] = (A).mat[2][0] * (x).xyz[0] \ 246 +(A).mat[2][1] * (x).xyz[1] \ 247 +(A).mat[2][2] * (x).xyz[2] + (b).xyz[2] ,\ 248 tempA_fvec3 ) 249 250 /* matrix-matrix multiply */ 251 252 #define ROW_DOT_COL(A,B,i,j) ( (A).mat[i][0] * (B).mat[0][j] \ 253 + (A).mat[i][1] * (B).mat[1][j] \ 254 + (A).mat[i][2] * (B).mat[2][j] ) 255 256 #define MAT_MUL(A,B) \ 257 ( tempA_mat33.mat[0][0] = ROW_DOT_COL((A),(B),0,0) , \ 258 tempA_mat33.mat[1][0] = ROW_DOT_COL((A),(B),1,0) , \ 259 tempA_mat33.mat[2][0] = ROW_DOT_COL((A),(B),2,0) , \ 260 tempA_mat33.mat[0][1] = ROW_DOT_COL((A),(B),0,1) , \ 261 tempA_mat33.mat[1][1] = ROW_DOT_COL((A),(B),1,1) , \ 262 tempA_mat33.mat[2][1] = ROW_DOT_COL((A),(B),2,1) , \ 263 tempA_mat33.mat[0][2] = ROW_DOT_COL((A),(B),0,2) , \ 264 tempA_mat33.mat[1][2] = ROW_DOT_COL((A),(B),1,2) , \ 265 tempA_mat33.mat[2][2] = ROW_DOT_COL((A),(B),2,2) , tempA_mat33 ) 266 267 #define MAT_SCALAR(A,c) \ 268 ( tempA_mat33.mat[0][0] = (c)*(A).mat[0][0] , \ 269 tempA_mat33.mat[1][0] = (c)*(A).mat[1][0] , \ 270 tempA_mat33.mat[2][0] = (c)*(A).mat[2][0] , \ 271 tempA_mat33.mat[0][1] = (c)*(A).mat[0][1] , \ 272 tempA_mat33.mat[1][1] = (c)*(A).mat[1][1] , \ 273 tempA_mat33.mat[2][1] = (c)*(A).mat[2][1] , \ 274 tempA_mat33.mat[0][2] = (c)*(A).mat[0][2] , \ 275 tempA_mat33.mat[1][2] = (c)*(A).mat[1][2] , \ 276 tempA_mat33.mat[2][2] = (c)*(A).mat[2][2] , tempA_mat33 ) 277 278 /* matrix determinant */ 279 280 #define MAT_DET(A) \ 281 ( (A).mat[0][0]*(A).mat[1][1]*(A).mat[2][2] \ 282 - (A).mat[0][0]*(A).mat[1][2]*(A).mat[2][1] \ 283 - (A).mat[1][0]*(A).mat[0][1]*(A).mat[2][2] \ 284 + (A).mat[1][0]*(A).mat[0][2]*(A).mat[2][1] \ 285 + (A).mat[2][0]*(A).mat[0][1]*(A).mat[1][2] \ 286 - (A).mat[2][0]*(A).mat[0][2]*(A).mat[1][1] ) 287 288 /* matrix norm [28 Aug 2002] */ 289 290 #define MAT_FNORM(A) \ 291 sqrt( (A).mat[0][0]*(A).mat[0][0] + \ 292 (A).mat[0][1]*(A).mat[0][1] + \ 293 (A).mat[0][2]*(A).mat[0][2] + \ 294 (A).mat[1][0]*(A).mat[1][0] + \ 295 (A).mat[1][1]*(A).mat[1][1] + \ 296 (A).mat[1][2]*(A).mat[1][2] + \ 297 (A).mat[2][0]*(A).mat[2][0] + \ 298 (A).mat[2][1]*(A).mat[2][1] + \ 299 (A).mat[2][2]*(A).mat[2][2] ) 300 301 /* matrix trace [5 Oct 1998] */ 302 303 #define MAT_TRACE(A) ( (A).mat[0][0] + (A).mat[1][1] + (A).mat[2][2] ) 304 305 /* matrix inverse */ 306 307 #define MAT_INV(A) \ 308 ( tempRWC = 1.0 / MAT_DET(A) , \ 309 tempA_mat33.mat[1][1] = \ 310 ( (A).mat[0][0]*(A).mat[2][2] - (A).mat[0][2]*(A).mat[2][0]) * tempRWC,\ 311 tempA_mat33.mat[2][2] = \ 312 ( (A).mat[0][0]*(A).mat[1][1] - (A).mat[0][1]*(A).mat[1][0]) * tempRWC,\ 313 tempA_mat33.mat[2][0] = \ 314 ( (A).mat[1][0]*(A).mat[2][1] - (A).mat[1][1]*(A).mat[2][0]) * tempRWC,\ 315 tempA_mat33.mat[1][2] = \ 316 (-(A).mat[0][0]*(A).mat[1][2] + (A).mat[0][2]*(A).mat[1][0]) * tempRWC,\ 317 tempA_mat33.mat[0][1] = \ 318 (-(A).mat[0][1]*(A).mat[2][2] + (A).mat[0][2]*(A).mat[2][1]) * tempRWC,\ 319 tempA_mat33.mat[0][0] = \ 320 ( (A).mat[1][1]*(A).mat[2][2] - (A).mat[1][2]*(A).mat[2][1]) * tempRWC,\ 321 tempA_mat33.mat[2][1] = \ 322 (-(A).mat[0][0]*(A).mat[2][1] + (A).mat[0][1]*(A).mat[2][0]) * tempRWC,\ 323 tempA_mat33.mat[1][0] = \ 324 (-(A).mat[1][0]*(A).mat[2][2] + (A).mat[1][2]*(A).mat[2][0]) * tempRWC,\ 325 tempA_mat33.mat[0][2] = \ 326 ( (A).mat[0][1]*(A).mat[1][2] - (A).mat[0][2]*(A).mat[1][1]) * tempRWC,\ 327 tempA_mat33 ) 328 329 /* load a matrix from scalars [3 Oct 1998] */ 330 331 #define LOAD_MAT(A,a11,a12,a13,a21,a22,a23,a31,a32,a33) \ 332 ( (A).mat[0][0] = (a11) , (A).mat[0][1] = (a12) , \ 333 (A).mat[0][2] = (a13) , (A).mat[1][0] = (a21) , \ 334 (A).mat[1][1] = (a22) , (A).mat[1][2] = (a23) , \ 335 (A).mat[2][0] = (a31) , (A).mat[2][1] = (a32) , (A).mat[2][2] = (a33) ) 336 337 /* unload a matrix into scalars [3 Oct 1998] */ 338 339 #define UNLOAD_MAT(A,a11,a12,a13,a21,a22,a23,a31,a32,a33) \ 340 ( (a11) = (A).mat[0][0] , (a12) = (A).mat[0][1] , \ 341 (a13) = (A).mat[0][2] , (a21) = (A).mat[1][0] , \ 342 (a22) = (A).mat[1][1] , (a23) = (A).mat[1][2] , \ 343 (a31) = (A).mat[2][0] , (a32) = (A).mat[2][1] , (a33) = (A).mat[2][2] ) 344 345 /* diagonal matrix */ 346 347 #define LOAD_DIAG_MAT(A,x,y,z) \ 348 ( (A).mat[0][0] = (x) , \ 349 (A).mat[1][1] = (y) , \ 350 (A).mat[2][2] = (z) , \ 351 (A).mat[0][1] = (A).mat[0][2] = (A).mat[1][0] = \ 352 (A).mat[1][2] = (A).mat[2][0] = (A).mat[2][1] = 0.0 ) 353 354 /* zero matrix */ 355 356 #define LOAD_ZERO_MAT(A) LOAD_DIAG_MAT((A),0,0,0) 357 358 /* elementary rotation matrices: 359 rotate about axis #ff, from axis #aa toward #bb, 360 where ff, aa, and bb are a permutation of {0,1,2} */ 361 362 #define LOAD_ROTGEN_MAT(A,th,ff,aa,bb) \ 363 ( (A).mat[aa][aa] = (A).mat[bb][bb] = cos((th)) , \ 364 (A).mat[aa][bb] = sin((th)) , \ 365 (A).mat[bb][aa] = -(A).mat[aa][bb] , \ 366 (A).mat[ff][ff] = 1.0 , \ 367 (A).mat[aa][ff] = (A).mat[bb][ff] = (A).mat[ff][aa] = (A).mat[ff][bb] = 0.0 ) 368 369 #define LOAD_ROTX_MAT(A,th) LOAD_ROTGEN_MAT(A,th,0,1,2) 370 #define LOAD_ROTY_MAT(A,th) LOAD_ROTGEN_MAT(A,th,1,2,0) 371 #define LOAD_ROTZ_MAT(A,th) LOAD_ROTGEN_MAT(A,th,2,0,1) 372 373 #define LOAD_ROT_MAT(A,th,i) \ 374 do{ switch( (i) ){ \ 375 case 0: LOAD_ROTX_MAT(A,th) ; break ; \ 376 case 1: LOAD_ROTY_MAT(A,th) ; break ; \ 377 case 2: LOAD_ROTZ_MAT(A,th) ; break ; \ 378 default: LOAD_DIAG_MAT(A,1,1,1); break ; \ 379 } } while(0) 380 381 /* shear matrices [3 Oct 1998] */ 382 383 #define LOAD_SHEARX_MAT(A,f,b,c) ( LOAD_DIAG_MAT(A,(f),1,1) , \ 384 (A).mat[0][1] = (b) , (A).mat[0][2] = (c) ) 385 386 #define LOAD_SHEARY_MAT(A,f,a,c) ( LOAD_DIAG_MAT(A,1,(f),1) , \ 387 (A).mat[1][0] = (a) , (A).mat[1][2] = (c) ) 388 389 #define LOAD_SHEARZ_MAT(A,f,a,b) ( LOAD_DIAG_MAT(A,1,1,(f)) , \ 390 (A).mat[2][0] = (a) , (A).mat[2][1] = (b) ) 391 392 /* matrix transpose */ 393 394 #define TRANSPOSE_MAT(A) \ 395 ( tempA_mat33.mat[0][0] = (A).mat[0][0] , \ 396 tempA_mat33.mat[1][0] = (A).mat[0][1] , \ 397 tempA_mat33.mat[2][0] = (A).mat[0][2] , \ 398 tempA_mat33.mat[0][1] = (A).mat[1][0] , \ 399 tempA_mat33.mat[1][1] = (A).mat[1][1] , \ 400 tempA_mat33.mat[2][1] = (A).mat[1][2] , \ 401 tempA_mat33.mat[0][2] = (A).mat[2][0] , \ 402 tempA_mat33.mat[1][2] = (A).mat[2][1] , \ 403 tempA_mat33.mat[2][2] = (A).mat[2][2] , tempA_mat33 ) 404 405 /* matrix add & subtract */ 406 407 #define ADD_MAT(A,B) \ 408 ( tempA_mat33.mat[0][0] = (A).mat[0][0] + (B).mat[0][0] , \ 409 tempA_mat33.mat[1][0] = (A).mat[1][0] + (B).mat[1][0] , \ 410 tempA_mat33.mat[2][0] = (A).mat[2][0] + (B).mat[2][0] , \ 411 tempA_mat33.mat[0][1] = (A).mat[0][1] + (B).mat[0][1] , \ 412 tempA_mat33.mat[1][1] = (A).mat[1][1] + (B).mat[1][1] , \ 413 tempA_mat33.mat[2][1] = (A).mat[2][1] + (B).mat[2][1] , \ 414 tempA_mat33.mat[0][2] = (A).mat[0][2] + (B).mat[0][2] , \ 415 tempA_mat33.mat[1][2] = (A).mat[1][2] + (B).mat[1][2] , \ 416 tempA_mat33.mat[2][2] = (A).mat[2][2] + (B).mat[2][2] , tempA_mat33 ) 417 418 #define SUB_MAT(A,B) \ 419 ( tempA_mat33.mat[0][0] = (A).mat[0][0] - (B).mat[0][0] , \ 420 tempA_mat33.mat[1][0] = (A).mat[1][0] - (B).mat[1][0] , \ 421 tempA_mat33.mat[2][0] = (A).mat[2][0] - (B).mat[2][0] , \ 422 tempA_mat33.mat[0][1] = (A).mat[0][1] - (B).mat[0][1] , \ 423 tempA_mat33.mat[1][1] = (A).mat[1][1] - (B).mat[1][1] , \ 424 tempA_mat33.mat[2][1] = (A).mat[2][1] - (B).mat[2][1] , \ 425 tempA_mat33.mat[0][2] = (A).mat[0][2] - (B).mat[0][2] , \ 426 tempA_mat33.mat[1][2] = (A).mat[1][2] - (B).mat[1][2] , \ 427 tempA_mat33.mat[2][2] = (A).mat[2][2] - (B).mat[2][2] , tempA_mat33 ) 428 429 /* component-wise min and max of 3-vectors */ 430 431 #define MAX_FVEC3(a,b) \ 432 ( tempA_fvec3.xyz[0] = (((a).xyz[0] > (b).xyz[0]) ? (a).xyz[0] : (b).xyz[0]) ,\ 433 tempA_fvec3.xyz[1] = (((a).xyz[1] > (b).xyz[1]) ? (a).xyz[1] : (b).xyz[1]) ,\ 434 tempA_fvec3.xyz[2] = (((a).xyz[2] > (b).xyz[2]) ? (a).xyz[2] : (b).xyz[2]) ,\ 435 tempA_fvec3 ) 436 437 #define MIN_FVEC3(a,b) \ 438 ( tempA_fvec3.xyz[0] = (((a).xyz[0] < (b).xyz[0]) ? (a).xyz[0] : (b).xyz[0]) ,\ 439 tempA_fvec3.xyz[1] = (((a).xyz[1] < (b).xyz[1]) ? (a).xyz[1] : (b).xyz[1]) ,\ 440 tempA_fvec3.xyz[2] = (((a).xyz[2] < (b).xyz[2]) ? (a).xyz[2] : (b).xyz[2]) ,\ 441 tempA_fvec3 ) 442 443 /*-------------------------------------------------------------------*/ 444 /*----- double precision versions of the above -----*/ 445 446 typedef struct { double xyz[3] ; } THD_dfvec3 ; 447 typedef struct { double mat[3][3] ; } THD_dmat33 ; 448 449 typedef struct { /* 3x3 matrix + 3-vector [16 Jul 2000] */ 450 THD_dfvec3 vv ; 451 THD_dmat33 mm ; 452 } THD_dvecmat ; 453 454 #define LOAD_DFVEC3(fv,x,y,z) ( (fv).xyz[0]=(x), \ 455 (fv).xyz[1]=(y), \ 456 (fv).xyz[2]=(z) ) 457 458 #define UNLOAD_DFVEC3(fv,x,y,z) ( (x)=(fv).xyz[0], \ 459 (y)=(fv).xyz[1], \ 460 (z)=(fv).xyz[2] ) 461 462 #define DFVEC3_TO_FVEC3(dv,fv) LOAD_FVEC3(fv,(dv).xyz[0],(dv).xyz[1],(dv).xyz[2]) 463 #define FVEC3_TO_DFVEC3(fv,dv) LOAD_DFVEC3(dv,(fv).xyz[0],(fv).xyz[1],(fv).xyz[2]) 464 465 static THD_dfvec3 tempA_dfvec3 , tempB_dfvec3 ; 466 static THD_dmat33 tempA_dmat33 , tempB_dmat33 ; 467 static double dtempRWC ; 468 469 #define TEMP_DFVEC3(x,y,z) ( tempB_dfvec3.xyz[0]=(x), \ 470 tempB_dfvec3.xyz[1]=(y), \ 471 tempB_dfvec3.xyz[2]=(z), tempB_dfvec3 ) 472 473 #define DUMP_DFVEC3(str,fv) \ 474 fprintf(stderr,"%s: %13.6g %13.6g %13.6g\n",(str),(fv).xyz[0],(fv).xyz[1],(fv).xyz[2]) 475 476 #define DUMP_DMAT33(str,A) \ 477 fprintf(stderr, \ 478 "%10.10s: [ %13.6g %13.6g %13.6g ]\n" \ 479 " [ %13.6g %13.6g %13.6g ]\n" \ 480 " [ %13.6g %13.6g %13.6g ]\n" , \ 481 str , (A).mat[0][0] , (A).mat[0][1] , (A).mat[0][2] , \ 482 (A).mat[1][0] , (A).mat[1][1] , (A).mat[1][2] , \ 483 (A).mat[2][0] , (A).mat[2][1] , (A).mat[2][2] ) 484 485 #define DUMP_DVECMAT(str,vvmm) ( DUMP_DMAT33(str,(vvmm).mm) , DUMP_DFVEC3(str,(vvmm).vv) ) 486 487 /*! Convert from dmat33 to mat33 (double to single precision) */ 488 489 #define DMAT_TO_MAT(D,M) \ 490 LOAD_MAT(M,(D).mat[0][0] , (D).mat[0][1] , (D).mat[0][2] , \ 491 (D).mat[1][0] , (D).mat[1][1] , (D).mat[1][2] , \ 492 (D).mat[2][0] , (D).mat[2][1] , (D).mat[2][2] ) 493 494 /*--- macros for operations on floating 3 vectors, 495 with heavy use of the comma operator and structure assignment! ---*/ 496 497 /* negation */ 498 499 #define NEGATE_DFVEC3(a) ( (a).xyz[0] = -(a).xyz[0] , \ 500 (a).xyz[1] = -(a).xyz[1] , \ 501 (a).xyz[2] = -(a).xyz[2] ) 502 503 /* subtraction */ 504 505 #define SUB_DFVEC3(a,b) \ 506 ( tempA_dfvec3.xyz[0] = (a).xyz[0] - (b).xyz[0] , \ 507 tempA_dfvec3.xyz[1] = (a).xyz[1] - (b).xyz[1] , \ 508 tempA_dfvec3.xyz[2] = (a).xyz[2] - (b).xyz[2] , tempA_dfvec3 ) 509 510 /* addition */ 511 512 #define ADD_DFVEC3(a,b) \ 513 ( tempA_dfvec3.xyz[0] = (a).xyz[0] + (b).xyz[0] , \ 514 tempA_dfvec3.xyz[1] = (a).xyz[1] + (b).xyz[1] , \ 515 tempA_dfvec3.xyz[2] = (a).xyz[2] + (b).xyz[2] , tempA_dfvec3 ) 516 517 /* make into a unit vector */ 518 519 #define NORMALIZE_DFVEC3(a) \ 520 ( dtempRWC = (a).xyz[0] * (a).xyz[0] \ 521 + (a).xyz[1] * (a).xyz[1] \ 522 + (a).xyz[2] * (a).xyz[2] , \ 523 dtempRWC = (dtempRWC > 0) ? (1.0/sqrt(dtempRWC)) : 0 , \ 524 tempA_dfvec3.xyz[0] = (a).xyz[0] * dtempRWC , \ 525 tempA_dfvec3.xyz[1] = (a).xyz[1] * dtempRWC , \ 526 tempA_dfvec3.xyz[2] = (a).xyz[2] * dtempRWC , tempA_dfvec3 ) 527 528 /* cross product */ 529 530 #define CROSS_DFVEC3(a,b) \ 531 ( tempA_dfvec3.xyz[0] = (a).xyz[1]*(b).xyz[2] - (a).xyz[2]*(b).xyz[1] , \ 532 tempA_dfvec3.xyz[1] = (a).xyz[2]*(b).xyz[0] - (a).xyz[0]*(b).xyz[2] , \ 533 tempA_dfvec3.xyz[2] = (a).xyz[0]*(b).xyz[1] - (a).xyz[1]*(b).xyz[0] , \ 534 tempA_dfvec3 ) 535 536 /* L2 norm */ 537 538 #define SIZE_DFVEC3(a) \ 539 sqrt((a).xyz[0]*(a).xyz[0]+(a).xyz[1]*(a).xyz[1]+(a).xyz[2]*(a).xyz[2]) 540 541 /* dot product */ 542 543 #define DOT_DFVEC3(a,b) \ 544 ((a).xyz[0]*(b).xyz[0] + (a).xyz[1]*(b).xyz[1] + (a).xyz[2]*(b).xyz[2]) 545 546 /* scale and add two vectors: fa * a + fb * b */ 547 548 #define SCLADD_DFVEC3(fa,a,fb,b) \ 549 ( tempA_dfvec3.xyz[0] = (fa)*(a).xyz[0] + (fb)*(b).xyz[0] , \ 550 tempA_dfvec3.xyz[1] = (fa)*(a).xyz[1] + (fb)*(b).xyz[1] , \ 551 tempA_dfvec3.xyz[2] = (fa)*(a).xyz[2] + (fb)*(b).xyz[2] , tempA_dfvec3 ) 552 553 /* round to an integer vector (assume non-negative) */ 554 555 #define INT_DFVEC3(a) \ 556 ( tempA_ivec3.ijk[0] = (a).xyz[0] + 0.5 , \ 557 tempA_ivec3.ijk[1] = (a).xyz[1] + 0.5 , \ 558 tempA_ivec3.ijk[1] = (a).xyz[2] + 0.5 , tempA_ivec3 ) ; 559 560 #define SIZE_DMAT(A) \ 561 ( fabs((A).mat[0][0]) + fabs((A).mat[0][1]) + fabs((A).mat[0][2]) \ 562 +fabs((A).mat[1][0]) + fabs((A).mat[1][1]) + fabs((A).mat[1][2]) \ 563 +fabs((A).mat[2][0]) + fabs((A).mat[2][1]) + fabs((A).mat[2][2]) ) 564 565 /* matrix-vector multiply */ 566 567 #define DMATVEC(A,x) \ 568 ( tempA_dfvec3.xyz[0] = (A).mat[0][0] * (x).xyz[0] \ 569 +(A).mat[0][1] * (x).xyz[1] \ 570 +(A).mat[0][2] * (x).xyz[2] ,\ 571 tempA_dfvec3.xyz[1] = (A).mat[1][0] * (x).xyz[0] \ 572 +(A).mat[1][1] * (x).xyz[1] \ 573 +(A).mat[1][2] * (x).xyz[2] ,\ 574 tempA_dfvec3.xyz[2] = (A).mat[2][0] * (x).xyz[0] \ 575 +(A).mat[2][1] * (x).xyz[1] \ 576 +(A).mat[2][2] * (x).xyz[2] , tempA_dfvec3 ) 577 578 /* matrix-vector multiply with subtract: output = A (x-b) */ 579 580 #define VECSUB_DMAT(A,x,b) \ 581 ( tempA_dfvec3.xyz[0] = (A).mat[0][0] * ((x).xyz[0]-(b).xyz[0]) \ 582 +(A).mat[0][1] * ((x).xyz[1]-(b).xyz[1]) \ 583 +(A).mat[0][2] * ((x).xyz[2]-(b).xyz[2]) ,\ 584 tempA_dfvec3.xyz[1] = (A).mat[1][0] * ((x).xyz[0]-(b).xyz[0]) \ 585 +(A).mat[1][1] * ((x).xyz[1]-(b).xyz[1]) \ 586 +(A).mat[1][2] * ((x).xyz[2]-(b).xyz[2]) ,\ 587 tempA_dfvec3.xyz[2] = (A).mat[2][0] * ((x).xyz[0]-(b).xyz[0]) \ 588 +(A).mat[2][1] * ((x).xyz[1]-(b).xyz[1]) \ 589 +(A).mat[2][2] * ((x).xyz[2]-(b).xyz[2]) ,\ 590 tempA_dfvec3 ) 591 592 /* matrix vector multiply with subtract after: A x - b */ 593 594 #define DMATVEC_SUB(A,x,b) \ 595 ( tempA_dfvec3.xyz[0] = (A).mat[0][0] * (x).xyz[0] \ 596 +(A).mat[0][1] * (x).xyz[1] \ 597 +(A).mat[0][2] * (x).xyz[2] - (b).xyz[0] ,\ 598 tempA_dfvec3.xyz[1] = (A).mat[1][0] * (x).xyz[0] \ 599 +(A).mat[1][1] * (x).xyz[1] \ 600 +(A).mat[1][2] * (x).xyz[2] - (b).xyz[1] ,\ 601 tempA_dfvec3.xyz[2] = (A).mat[2][0] * (x).xyz[0] \ 602 +(A).mat[2][1] * (x).xyz[1] \ 603 +(A).mat[2][2] * (x).xyz[2] - (b).xyz[2] ,\ 604 tempA_dfvec3 ) 605 606 /* matrix vector multiply with add after: A x + b */ 607 608 #define DMATVEC_ADD(A,x,b) \ 609 ( tempA_dfvec3.xyz[0] = (A).mat[0][0] * (x).xyz[0] \ 610 +(A).mat[0][1] * (x).xyz[1] \ 611 +(A).mat[0][2] * (x).xyz[2] + (b).xyz[0] ,\ 612 tempA_dfvec3.xyz[1] = (A).mat[1][0] * (x).xyz[0] \ 613 +(A).mat[1][1] * (x).xyz[1] \ 614 +(A).mat[1][2] * (x).xyz[2] + (b).xyz[1] ,\ 615 tempA_dfvec3.xyz[2] = (A).mat[2][0] * (x).xyz[0] \ 616 +(A).mat[2][1] * (x).xyz[1] \ 617 +(A).mat[2][2] * (x).xyz[2] + (b).xyz[2] ,\ 618 tempA_dfvec3 ) 619 620 /* matrix-matrix multiply */ 621 622 #define ROW_DOT_COL(A,B,i,j) ( (A).mat[i][0] * (B).mat[0][j] \ 623 + (A).mat[i][1] * (B).mat[1][j] \ 624 + (A).mat[i][2] * (B).mat[2][j] ) 625 626 #define DMAT_MUL(A,B) \ 627 ( tempA_dmat33.mat[0][0] = ROW_DOT_COL((A),(B),0,0) , \ 628 tempA_dmat33.mat[1][0] = ROW_DOT_COL((A),(B),1,0) , \ 629 tempA_dmat33.mat[2][0] = ROW_DOT_COL((A),(B),2,0) , \ 630 tempA_dmat33.mat[0][1] = ROW_DOT_COL((A),(B),0,1) , \ 631 tempA_dmat33.mat[1][1] = ROW_DOT_COL((A),(B),1,1) , \ 632 tempA_dmat33.mat[2][1] = ROW_DOT_COL((A),(B),2,1) , \ 633 tempA_dmat33.mat[0][2] = ROW_DOT_COL((A),(B),0,2) , \ 634 tempA_dmat33.mat[1][2] = ROW_DOT_COL((A),(B),1,2) , \ 635 tempA_dmat33.mat[2][2] = ROW_DOT_COL((A),(B),2,2) , tempA_dmat33 ) 636 637 #define DMAT_SCALAR(A,c) \ 638 ( tempA_dmat33.mat[0][0] = (c)*(A).mat[0][0] , \ 639 tempA_dmat33.mat[1][0] = (c)*(A).mat[1][0] , \ 640 tempA_dmat33.mat[2][0] = (c)*(A).mat[2][0] , \ 641 tempA_dmat33.mat[0][1] = (c)*(A).mat[0][1] , \ 642 tempA_dmat33.mat[1][1] = (c)*(A).mat[1][1] , \ 643 tempA_dmat33.mat[2][1] = (c)*(A).mat[2][1] , \ 644 tempA_dmat33.mat[0][2] = (c)*(A).mat[0][2] , \ 645 tempA_dmat33.mat[1][2] = (c)*(A).mat[1][2] , \ 646 tempA_dmat33.mat[2][2] = (c)*(A).mat[2][2] , tempA_dmat33 ) 647 648 /* matrix determinant */ 649 650 #define DMAT_DET(A) \ 651 ( (A).mat[0][0]*(A).mat[1][1]*(A).mat[2][2] \ 652 - (A).mat[0][0]*(A).mat[1][2]*(A).mat[2][1] \ 653 - (A).mat[1][0]*(A).mat[0][1]*(A).mat[2][2] \ 654 + (A).mat[1][0]*(A).mat[0][2]*(A).mat[2][1] \ 655 + (A).mat[2][0]*(A).mat[0][1]*(A).mat[1][2] \ 656 - (A).mat[2][0]*(A).mat[0][2]*(A).mat[1][1] ) 657 658 /* matrix trace [5 Oct 1998] */ 659 660 #define DMAT_TRACE(A) ( (A).mat[0][0] + (A).mat[1][1] + (A).mat[2][2] ) 661 662 /* matrix inverse */ 663 664 #define DMAT_INV(A) \ 665 ( dtempRWC = 1.0 / DMAT_DET(A) , \ 666 tempA_dmat33.mat[1][1] = \ 667 ( (A).mat[0][0]*(A).mat[2][2] - (A).mat[0][2]*(A).mat[2][0]) * dtempRWC,\ 668 tempA_dmat33.mat[2][2] = \ 669 ( (A).mat[0][0]*(A).mat[1][1] - (A).mat[0][1]*(A).mat[1][0]) * dtempRWC,\ 670 tempA_dmat33.mat[2][0] = \ 671 ( (A).mat[1][0]*(A).mat[2][1] - (A).mat[1][1]*(A).mat[2][0]) * dtempRWC,\ 672 tempA_dmat33.mat[1][2] = \ 673 (-(A).mat[0][0]*(A).mat[1][2] + (A).mat[0][2]*(A).mat[1][0]) * dtempRWC,\ 674 tempA_dmat33.mat[0][1] = \ 675 (-(A).mat[0][1]*(A).mat[2][2] + (A).mat[0][2]*(A).mat[2][1]) * dtempRWC,\ 676 tempA_dmat33.mat[0][0] = \ 677 ( (A).mat[1][1]*(A).mat[2][2] - (A).mat[1][2]*(A).mat[2][1]) * dtempRWC,\ 678 tempA_dmat33.mat[2][1] = \ 679 (-(A).mat[0][0]*(A).mat[2][1] + (A).mat[0][1]*(A).mat[2][0]) * dtempRWC,\ 680 tempA_dmat33.mat[1][0] = \ 681 (-(A).mat[1][0]*(A).mat[2][2] + (A).mat[1][2]*(A).mat[2][0]) * dtempRWC,\ 682 tempA_dmat33.mat[0][2] = \ 683 ( (A).mat[0][1]*(A).mat[1][2] - (A).mat[0][2]*(A).mat[1][1]) * dtempRWC,\ 684 tempA_dmat33 ) 685 686 /* load a matrix from scalars [3 Oct 1998] */ 687 688 #define LOAD_DMAT(A,a11,a12,a13,a21,a22,a23,a31,a32,a33) \ 689 ( (A).mat[0][0] = (a11) , (A).mat[0][1] = (a12) , \ 690 (A).mat[0][2] = (a13) , (A).mat[1][0] = (a21) , \ 691 (A).mat[1][1] = (a22) , (A).mat[1][2] = (a23) , \ 692 (A).mat[2][0] = (a31) , (A).mat[2][1] = (a32) , (A).mat[2][2] = (a33) ) 693 694 /* unload a matrix into scalars [3 Oct 1998] */ 695 696 #define UNLOAD_DMAT(A,a11,a12,a13,a21,a22,a23,a31,a32,a33) \ 697 ( (a11) = (A).mat[0][0] , (a12) = (A).mat[0][1] , \ 698 (a13) = (A).mat[0][2] , (a21) = (A).mat[1][0] , \ 699 (a22) = (A).mat[1][1] , (a23) = (A).mat[1][2] , \ 700 (a31) = (A).mat[2][0] , (a32) = (A).mat[2][1] , (a33) = (A).mat[2][2] ) 701 702 /* diagonal matrix */ 703 704 #define LOAD_DIAG_DMAT(A,x,y,z) \ 705 ( (A).mat[0][0] = (x) , \ 706 (A).mat[1][1] = (y) , \ 707 (A).mat[2][2] = (z) , \ 708 (A).mat[0][1] = (A).mat[0][2] = (A).mat[1][0] = \ 709 (A).mat[1][2] = (A).mat[2][0] = (A).mat[2][1] = 0.0 ) 710 711 /* zero matrix */ 712 713 #define LOAD_ZERO_DMAT(A) LOAD_DIAG_DMAT((A),0,0,0) 714 715 /* elementary rotation matrices: 716 rotate about axis #ff, from axis #aa toward #bb, 717 where ff, aa, and bb are a permutation of {0,1,2} */ 718 719 #define LOAD_ROTGEN_DMAT(A,th,ff,aa,bb) \ 720 ( (A).mat[aa][aa] = (A).mat[bb][bb] = cos((th)) , \ 721 (A).mat[aa][bb] = sin((th)) , \ 722 (A).mat[bb][aa] = -(A).mat[aa][bb] , \ 723 (A).mat[ff][ff] = 1.0 , \ 724 (A).mat[aa][ff] = (A).mat[bb][ff] = (A).mat[ff][aa] = (A).mat[ff][bb] = 0.0 ) 725 726 #define LOAD_ROTX_DMAT(A,th) LOAD_ROTGEN_DMAT(A,th,0,1,2) 727 #define LOAD_ROTY_DMAT(A,th) LOAD_ROTGEN_DMAT(A,th,1,2,0) 728 #define LOAD_ROTZ_DMAT(A,th) LOAD_ROTGEN_DMAT(A,th,2,0,1) 729 730 #define LOAD_ROT_DMAT(A,th,i) \ 731 do{ switch( (i) ){ \ 732 case 0: LOAD_ROTX_DMAT(A,th) ; break ; \ 733 case 1: LOAD_ROTY_DMAT(A,th) ; break ; \ 734 case 2: LOAD_ROTZ_DMAT(A,th) ; break ; \ 735 default: LOAD_ZERO_DMAT(A) ; break ; \ 736 } } while(0) 737 738 /* shear matrices [3 Oct 1998] */ 739 740 #define LOAD_SHEARX_DMAT(A,f,b,c) ( LOAD_DIAG_DMAT(A,(f),1,1) , \ 741 (A).mat[0][1] = (b) , (A).mat[0][2] = (c) ) 742 743 #define LOAD_SHEARY_DMAT(A,f,a,c) ( LOAD_DIAG_DMAT(A,1,(f),1) , \ 744 (A).mat[1][0] = (a) , (A).mat[1][2] = (c) ) 745 746 #define LOAD_SHEARZ_DMAT(A,f,a,b) ( LOAD_DIAG_DMAT(A,1,1,(f)) , \ 747 (A).mat[2][0] = (a) , (A).mat[2][1] = (b) ) 748 749 /* matrix transpose */ 750 751 #define TRANSPOSE_DMAT(A) \ 752 ( tempA_dmat33.mat[0][0] = (A).mat[0][0] , \ 753 tempA_dmat33.mat[1][0] = (A).mat[0][1] , \ 754 tempA_dmat33.mat[2][0] = (A).mat[0][2] , \ 755 tempA_dmat33.mat[0][1] = (A).mat[1][0] , \ 756 tempA_dmat33.mat[1][1] = (A).mat[1][1] , \ 757 tempA_dmat33.mat[2][1] = (A).mat[1][2] , \ 758 tempA_dmat33.mat[0][2] = (A).mat[2][0] , \ 759 tempA_dmat33.mat[1][2] = (A).mat[2][1] , \ 760 tempA_dmat33.mat[2][2] = (A).mat[2][2] , tempA_dmat33 ) 761 762 /* matrix add & subtract */ 763 764 #define ADD_DMAT(A,B) \ 765 ( tempA_dmat33.mat[0][0] = (A).mat[0][0] + (B).mat[0][0] , \ 766 tempA_dmat33.mat[1][0] = (A).mat[1][0] + (B).mat[1][0] , \ 767 tempA_dmat33.mat[2][0] = (A).mat[2][0] + (B).mat[2][0] , \ 768 tempA_dmat33.mat[0][1] = (A).mat[0][1] + (B).mat[0][1] , \ 769 tempA_dmat33.mat[1][1] = (A).mat[1][1] + (B).mat[1][1] , \ 770 tempA_dmat33.mat[2][1] = (A).mat[2][1] + (B).mat[2][1] , \ 771 tempA_dmat33.mat[0][2] = (A).mat[0][2] + (B).mat[0][2] , \ 772 tempA_dmat33.mat[1][2] = (A).mat[1][2] + (B).mat[1][2] , \ 773 tempA_dmat33.mat[2][2] = (A).mat[2][2] + (B).mat[2][2] , tempA_dmat33 ) 774 775 #define SUB_DMAT(A,B) \ 776 ( tempA_dmat33.mat[0][0] = (A).mat[0][0] - (B).mat[0][0] , \ 777 tempA_dmat33.mat[1][0] = (A).mat[1][0] - (B).mat[1][0] , \ 778 tempA_dmat33.mat[2][0] = (A).mat[2][0] - (B).mat[2][0] , \ 779 tempA_dmat33.mat[0][1] = (A).mat[0][1] - (B).mat[0][1] , \ 780 tempA_dmat33.mat[1][1] = (A).mat[1][1] - (B).mat[1][1] , \ 781 tempA_dmat33.mat[2][1] = (A).mat[2][1] - (B).mat[2][1] , \ 782 tempA_dmat33.mat[0][2] = (A).mat[0][2] - (B).mat[0][2] , \ 783 tempA_dmat33.mat[1][2] = (A).mat[1][2] - (B).mat[1][2] , \ 784 tempA_dmat33.mat[2][2] = (A).mat[2][2] - (B).mat[2][2] , tempA_dmat33 ) 785 786 /* component-wise min and max of 3-vectors */ 787 788 #define MAX_DFVEC3(a,b) \ 789 ( tempA_dfvec3.xyz[0] = (((a).xyz[0] > (b).xyz[0]) ? (a).xyz[0] : (b).xyz[0]) ,\ 790 tempA_dfvec3.xyz[1] = (((a).xyz[1] > (b).xyz[1]) ? (a).xyz[1] : (b).xyz[1]) ,\ 791 tempA_dfvec3.xyz[2] = (((a).xyz[2] > (b).xyz[2]) ? (a).xyz[2] : (b).xyz[2]) ,\ 792 tempA_dfvec3 ) 793 794 #define MIN_DFVEC3(a,b) \ 795 ( tempA_dfvec3.xyz[0] = (((a).xyz[0] < (b).xyz[0]) ? (a).xyz[0] : (b).xyz[0]) ,\ 796 tempA_dfvec3.xyz[1] = (((a).xyz[1] < (b).xyz[1]) ? (a).xyz[1] : (b).xyz[1]) ,\ 797 tempA_dfvec3.xyz[2] = (((a).xyz[2] < (b).xyz[2]) ? (a).xyz[2] : (b).xyz[2]) ,\ 798 tempA_dfvec3 ) 799 800 /* multiply two vecmat structures */ 801 802 static THD_dvecmat tempA_dvm33 ; 803 804 #define MUL_DVECMAT(A,B) \ 805 ( tempA_dvm33.mm = DMAT_MUL((A).mm,(B).mm) , \ 806 tempB_dfvec3 = DMATVEC((A).mm,(B).vv) , \ 807 tempA_dvm33.vv = ADD_DFVEC3(tempB_dfvec3,(A).vv) , tempA_dvm33 ) 808 809 /* invert a vecmat structure: */ 810 /* [y] = [R][x] + [v] is transformation, so */ 811 /* [x] = inv[R] - inv[R][v] is inverse transformation */ 812 813 #define INV_DVECMAT(A) ( tempA_dvm33.mm = DMAT_INV((A).mm) , \ 814 tempA_dvm33.vv = DMATVEC(tempA_dvm33.mm,(A).vv) , \ 815 NEGATE_DFVEC3(tempA_dvm33.vv) , tempA_dvm33 ) 816 817 /* same for single precision stuctures */ 818 819 static THD_vecmat tempA_vm33 ; 820 821 #define MUL_VECMAT(A,B) \ 822 ( tempA_vm33.mm = MAT_MUL((A).mm,(B).mm) , \ 823 tempB_fvec3 = MATVEC((A).mm,(B).vv) , \ 824 tempA_vm33.vv = ADD_FVEC3(tempB_fvec3,(A).vv) , tempA_vm33 ) 825 826 #define INV_VECMAT(A) ( tempA_vm33.mm = MAT_INV((A).mm) , \ 827 tempA_vm33.vv = MATVEC(tempA_vm33.mm,(A).vv) , \ 828 NEGATE_FVEC3(tempA_vm33.vv) , tempA_vm33 ) 829 830 /* apply a vecmat to a vector */ 831 832 #define VECMAT_VEC(A,x) MATVEC_ADD( (A).mm , (x) , (A).vv ) 833 #define DVECMAT_VEC(A,x) DMATVEC_ADD( (A).mm , (x) , (A).vv ) 834 835 #endif /* _MCW_3DVECMAT_ */ 836