1 #ifndef SUMA_MACROSm_INCLUDED 2 #define SUMA_MACROSm_INCLUDED 3 4 /* A lazy macro to free a DList */ 5 #define SUMA_FREE_DLIST(striplist) { \ 6 if (striplist) { dlist_destroy(striplist); SUMA_free(striplist); } \ 7 striplist = NULL; \ 8 } 9 10 11 #ifdef isfinite 12 # define SUMA_IS_GOOD_FLOAT(x) isfinite(x) /* stolen from Bob's thd_floatscan.c */ 13 #else 14 # define SUMA_IS_GOOD_FLOAT(x) finite(x) 15 #endif 16 17 /*! sqrt(2 pi) and 1/sqrt(2 pi) */ 18 #define SQ_2PI 2.50662827463100024161 19 #define iSQ_2PI 0.39894228040143270286 20 21 22 /*! Return a random number from a gaussian distribution of mean m and std s 23 Based on parser.c and parser_int.c functions unif_ and gran_1*/ 24 #define SUMA_GRAN(m,s) \ 25 ( m + \ 26 (drand48() + drand48() + drand48() + drand48() + \ 27 drand48() + drand48() + drand48() + drand48() + \ 28 drand48() + drand48() + drand48() + drand48() - 6.0) * s) 29 30 31 /* return integer between 0 and t-1 */ 32 #define SUMA_IRAN(t) (lrand48() % (t)) 33 34 /*! get a random color */ 35 #define SUMA_RAND_COL(i,r,g,b,a) {\ 36 if ((i)>=0) srand((unsigned)(i)); \ 37 (r) = (float) (double)rand()/(double)RAND_MAX; \ 38 (g) = (float) (double)rand()/(double)RAND_MAX; \ 39 (b) = (float) (double)rand()/(double)RAND_MAX; \ 40 (a) = (float) (double)rand()/(double)RAND_MAX; \ 41 } 42 43 /*! 44 Macro to create a new ID code at random (when strn = NULL) or a hash of a string 45 written mostly to allow for the allocation of newcode with SUMA's functions 46 */ 47 #define SUMA_NEW_ID(newcode, istrn) { \ 48 char *strn=(char*)istrn; /* a trick to quiet compiler warnings like \ 49 warning: the address of 'stmp' will always evaluate as 'true'\ 50 when 'stmp' is a string of the type char [] */\ 51 if ((newcode)) { SUMA_S_Err("newcode pointer must be null"); } \ 52 else if (!(strn)) { (newcode) = (char*)SUMA_calloc(SUMA_IDCODE_LENGTH, sizeof(char)); UNIQ_idcode_fill((newcode)); } \ 53 else { char *m_tmp; m_tmp = UNIQ_hashcode((strn)); (newcode) = SUMA_copy_string(m_tmp); free(m_tmp); m_tmp = NULL; } \ 54 } 55 56 #define SUMA_PUT_ID_ATTR(nel,name,strn) { \ 57 if (strn) { \ 58 char *m_tmp; m_tmp = UNIQ_hashcode((strn)); \ 59 NI_set_attribute(nel, name,m_tmp); free(m_tmp); \ 60 } else { \ 61 char m_sbuf[SUMA_IDCODE_LENGTH]; \ 62 UNIQ_idcode_fill(m_sbuf); \ 63 NI_set_attribute(nel, name, m_sbuf); \ 64 } \ 65 } 66 67 68 #define SUMA_WHAT_ENDIAN(End){ \ 69 int m_one = 1; \ 70 /* From RickR's Imon */ \ 71 End = (*(char *)&m_one == 1) ? LSB_FIRST : MSB_FIRST; \ 72 } 73 74 #define SUMA_OTHER_ENDIAN(End){ \ 75 End = (End == LSB_FIRST) ? MSB_FIRST : LSB_FIRST; \ 76 } 77 78 /* nonzero if not power of 2 (from GLUT-3.7's advanced97/volume.c)*/ 79 #define SUMA_NOTPOW2(num) ((num) & (num - 1)) 80 81 #define SUMA_SWAP_THIS(nip,chnk){ \ 82 if (chnk == 4) SUMA_swap_4( nip ) ; \ 83 else if (chnk == 8) SUMA_swap_8( nip ) ; \ 84 else if (chnk == 2) SUMA_swap_2( nip ) ; \ 85 else { SUMA_SL_Err ("No swapping performed.") } \ 86 } 87 /*! 88 \brief a macro for reading one number at a time from a file 89 \param nip (void *) where values go 90 \param fp (FILE *) 91 \param ex (int) value returned by fread 92 \param chnk (int) sizeof(TYPE_YOU_READ) 93 94 \sa SUMA_READ_NUM_BS 95 */ 96 #define SUMA_READ_NUM(nip, fp, ex, chnk) \ 97 { \ 98 ex = fread (nip, chnk, 1, fp); \ 99 } 100 /*! 101 SUMA_READ_NUM with swapping 102 */ 103 #define SUMA_READ_NUM_BS(nip, fp, ex, chnk) \ 104 { \ 105 SUMA_READ_NUM(nip, fp, ex, chnk); \ 106 if (chnk == 4) SUMA_swap_4( nip ) ; \ 107 else if (chnk == 8) SUMA_swap_8( nip ) ; \ 108 else if (chnk == 2) SUMA_swap_2( nip ) ; \ 109 else { SUMA_SL_Err ("No swapping performed.") } \ 110 } 111 /*! 112 \brief a macro for reading one integer from a file pointer. 113 114 \param nip (int *) 115 \param bs (int) 0: no swap, 1 swap 116 \param fp (FILE *) 117 \param ex (int) value returned by fread 118 */ 119 #define SUMA_READ_INT(nip, bs, fp, ex) \ 120 { static int m_chnk = sizeof(int);\ 121 ex = fread (nip, m_chnk, 1, fp); \ 122 if (bs) { \ 123 if (m_chnk == 4) SUMA_swap_4( nip ) ; \ 124 else if (m_chnk == 8) SUMA_swap_8( nip ) ; \ 125 else if (m_chnk == 2) SUMA_swap_2( nip ) ; /* keep compiler quiet about using swap_2 */ \ 126 else { SUMA_SL_Err ("No swapping performed.") } \ 127 } \ 128 } 129 #define SUMA_READ_FLOAT(nip, bs, fp, ex) \ 130 { static int m_chnk = sizeof(float);\ 131 ex = fread (nip, m_chnk, 1, fp); \ 132 if (bs) { \ 133 if (m_chnk == 4) SUMA_swap_4( nip ) ; \ 134 else if (m_chnk == 8) SUMA_swap_8( nip ) ; \ 135 else { SUMA_SL_Err ("No swapping performed.") } \ 136 } \ 137 } 138 139 #define SUMA_SWAP_VEC(vec,N_alloc,chnk) {\ 140 int m_i; \ 141 if (chnk == 4) { for (m_i=0; m_i<N_alloc; ++m_i) SUMA_swap_4(&(vec[m_i])); } \ 142 else if (chnk == 2) { for (m_i=0; m_i<N_alloc; ++m_i) SUMA_swap_2(&(vec[m_i])); } \ 143 else if (chnk == 8) { for (m_i=0; m_i<N_alloc; ++m_i) SUMA_swap_8(&(vec[m_i])); } \ 144 else { SUMA_SL_Err("Bad chnk"); } \ 145 } 146 147 #define SUMA_NODELIST_XY_NEGATE(v, n) { \ 148 int m_n3 = 3*n, m_n=0; \ 149 while (m_n < m_n3) {\ 150 v[m_n] = -v[m_n]; m_n++;\ 151 v[m_n] = -v[m_n]; m_n += 2;\ 152 }\ 153 } 154 155 #define SUMA_CHECK_NULL_STR(str) ((str) ? (str) : "(NULL)") 156 #define CNS SUMA_CHECK_NULL_STR 157 158 #define SUMA_IS_STRICT_POS(a) ( ((a) > 0) ? 1 : 0 ) 159 160 #define SUMA_IS_POS(a) ( ((a) >= 0) ? 1 : 0 ) 161 162 #define SUMA_IS_STRICT_NEG(a) ( ((a) < 0) ? 1 : 0 ) 163 164 #define SUMA_IS_NEG(a) ( ((a) <= 0) ? 1 : 0 ) 165 166 #define SUMA_SIGN(a) ( ((a) < 0) ? -1 : 1 ) 167 168 #define SUMA_SIGN_CHAR(p) ( ( (p) < 0 ) ? '-' : '+' ) 169 170 #define SUMA_MIN_PAIR(a,b) ( ((a) <= (b)) ? (a) : (b) ) 171 172 #define SUMA_MAX_PAIR(a,b) ( ((a) <= (b)) ? (b) : (a) ) 173 174 #define SUMA_ABS(a) ( ((a) < 0 ) ? (-(a)) : (a) ) 175 176 #define SUMA_COMPLEX_ADD(a,b,m) { \ 177 (m).r = (a).r+(b).r; \ 178 (m).i = (a).i+(b).i; } 179 #define SUMA_COMPLEX_REAL_ADD(a,b,m) { \ 180 (m).r = (a).r+(b); \ 181 (m).i = (a).i; } 182 #define SUMA_COMPLEX_MULT(a,b,m) { \ 183 (m).r = (a).r*(b).r - (a).i*(b).i; \ 184 (m).i = (a).i*(b).r + (a).r*(b).i; } 185 #define SUMA_COMPLEX_SCALE(a,f,m) {\ 186 (m).r = (a).r*(f); \ 187 (m).i = (a).i*(f); } 188 #define SUMA_COMPLEX_S2(a) ((a).r*(a).r+(a).i*(a).i) 189 #define SUMA_COMPLEX_ABS(a) sqrt(SUMA_COMPLEX_S2(a)) 190 191 #define SUMA_COMPLEX_INV(a,ai) { \ 192 double m_a=SUMA_COMPLEX_S2(a); \ 193 ai.r = a.r/m_a; ai.i = -a.i/m_a; \ 194 } 195 196 /*! Transpose a 4x4 matrix */ 197 #define SUMA_TRANSP_4MATRIX(mm){\ 198 double tt;\ 199 tt=mm[1][0]; mm[1][0] = mm[0][1]; mm[0][1]=tt; \ 200 tt=mm[2][0]; mm[2][0] = mm[0][2]; mm[0][2]=tt; \ 201 tt=mm[3][0]; mm[3][0] = mm[0][3]; mm[0][3]=tt; \ 202 tt=mm[1][2]; mm[1][2] = mm[2][1]; mm[2][1]=tt; \ 203 tt=mm[1][3]; mm[1][3] = mm[3][1]; mm[3][1]=tt; \ 204 tt=mm[3][2]; mm[3][2] = mm[2][3]; mm[2][3]=tt; \ 205 } 206 207 /*! invert a 4x4 affine xform matrix */ 208 #define SUMA_INV_44ROTATIONMATRIX(xform) \ 209 { \ 210 static mat44 A, A0; \ 211 LOAD_MAT44( A0, \ 212 xform[0][0], xform[0][1], xform[0][2], xform[0][3], \ 213 xform[1][0], xform[1][1], xform[1][2], xform[1][3], \ 214 xform[2][0], xform[2][1], xform[2][2], xform[2][3] ); \ 215 A = nifti_mat44_inverse(A0); \ 216 UNLOAD_MAT44(A, \ 217 xform[0][0], xform[0][1], xform[0][2], xform[0][3], \ 218 xform[1][0], xform[1][1], xform[1][2], xform[1][3], \ 219 xform[2][0], xform[2][1], xform[2][2], xform[2][3] ); \ 220 } 221 222 #define SUMA_INV_V12MATRIX(vxform) \ 223 { \ 224 static mat44 A, A0; \ 225 LOAD_MAT44( A0, \ 226 V[0], V[1], V[2], V[3], \ 227 V[4], V[5], V[6], V[7], \ 228 V[8], V[9], V[10], V[11] ); \ 229 A = nifti_mat44_inverse(A0); \ 230 UNLOAD_MAT44(A, \ 231 V[0], V[1], V[2], V[3], \ 232 V[4], V[5], V[6], V[7], \ 233 V[8], V[9], V[10], V[11] ); \ 234 } 235 236 #define SUMA_POW2(a) ((a)*(a)) 237 238 #define SUMA_POW3(a) ((a)*(a)*(a)) 239 240 #define SUMA_R2D(a) ( (a)*180.0/SUMA_PI ) 241 242 #define SUMA_ROUND(a) ( ((a) < 0 ) ? (int)((a)-0.5) : (int)((a)+0.5) ) 243 244 #define SUMA_CEIL_POS(a) ( ( ((a) - (int)(a)) == 0.0 ) ? (int)(a) : ((int)(a)+1) ) 245 #define SUMA_CEIL(a) ( ((a) < 0 ) ? (int)(a) : SUMA_CEIL_POS(a) ) 246 247 #define SUMA_3D_2_1D_index AFNI_3D_to_1D_index 248 249 #define SUMA_1D_2_3D_index AFNI_1D_to_3D_index 250 251 /*! 252 \brief Returns the two points that are at a distance d from P1 along the direction of U 253 SUMA_POINT_AT_DISTANCE(U, P1, d, P2) 254 Input paramters : 255 \param U (float *) 3x1 vector specifying directions along x, y, z axis 256 U does not have to be normalized 257 \param P1 (float *) 3x1 vector containing the XYZ of P1 258 \param d (float) distance from P1 259 \param P2 (float 2x3) 2D array to hold the two points equidistant from P1 260 along U (first row) and -U (second row). 261 I recommend you initialize P2 to all zeros. 262 263 \sa SUMA_POINT_AT_DISTANCE_NORM 264 265 { 266 float P1[3] = { 54.255009, 85.570129, -4.534704 }; 267 float Un[3] = { -0.525731, -0.850651, 0.000000 }; 268 float U[3] = { -2.6287 , -4.2533 , 0 }; 269 float d = 25, P2[2][3]; 270 271 SUMA_POINT_AT_DISTANCE(U, P1, d, P2); 272 fprintf(SUMA_STDERR,"P2 [%f %f %f]\n [%f %f %f]\n", P2[0][0], P2[0][1], P2[0][2], P2[1][0], P2[1][1], P2[1][2]); 273 SUMA_POINT_AT_DISTANCE_NORM(Un, P1, d, P2); use this macro if you have a normalized direction vector ... 274 fprintf(SUMA_STDERR,"P2 [%f %f %f]\n [%f %f %f]\n", P2[0][0], P2[0][1], P2[0][2], P2[1][0], P2[1][1], P2[1][2]); 275 } 276 277 */ 278 #define SUMA_POINT_AT_DISTANCE(U, P1, d, P2){ \ 279 float m_n, m_Un[3]; \ 280 SUMA_NORM_VEC(U, 3, m_n); \ 281 if (m_n) { \ 282 if (m_n != 1.0) { \ 283 m_Un[0] = (U[0]) / (m_n); m_Un[1] = (U[1]) / (m_n); m_Un[2] = (U[2]) / (m_n); \ 284 SUMA_POINT_AT_DISTANCE_NORM(m_Un, P1, d, P2); \ 285 } else { SUMA_POINT_AT_DISTANCE_NORM(U, P1, d, P2); } \ 286 } \ 287 } 288 /*! 289 \brief Returns the two points that are at a distance d from P1 along the direction of U 290 SUMA_POINT_AT_DISTANCE_NORM(U, P1, d, P2) 291 Input paramters : 292 \param U (float *) 3x1 vector specifying directions along x, y, z axis 293 U MUST BE A UNIT VECTOR 294 \param P1 (float *) 3x1 vector containing the XYZ of P1 295 \param d (float) distance from P1 296 \param P2 (float 2x3) 2D array to hold the two points equidistant from P1 297 along U (first row) and -U (second row). 298 I recommend you initialize P2 to all zeros. 299 300 \sa SUMA_POINT_AT_DISTANCE 301 302 { 303 float P1[3] = { 54.255009, 85.570129, -4.534704 }; 304 float Un[3] = { -0.525731, -0.850651, 0.000000 }; 305 float U[3] = { -2.6287 , -4.2533 , 0 }; 306 float d = 25, P2[2][3]; 307 308 SUMA_POINT_AT_DISTANCE(U, P1, d, P2); 309 fprintf(SUMA_STDERR,"P2 [%f %f %f]\n [%f %f %f]\n", 310 P2[0][0], P2[0][1], P2[0][2], P2[1][0], P2[1][1], P2[1][2]); 311 SUMA_POINT_AT_DISTANCE_NORM(Un, P1, d, P2); use this macro if you have a normalized direction vector ... 312 fprintf(SUMA_STDERR,"P2 [%f %f %f]\n [%f %f %f]\n", 313 P2[0][0], P2[0][1], P2[0][2], P2[1][0], P2[1][1], P2[1][2]); 314 } 315 316 */ 317 #define SUMA_POINT_AT_DISTANCE_NORM(U, P1, d, P2){ \ 318 P2[0][0] = (d) * U[0]; P2[1][0] = -(d) * U[0]; P2[0][0] += P1[0]; P2[1][0] += P1[0]; \ 319 P2[0][1] = (d) * U[1]; P2[1][1] = -(d) * U[1]; P2[0][1] += P1[1]; P2[1][1] += P1[1]; \ 320 P2[0][2] = (d) * U[2]; P2[1][2] = -(d) * U[2]; P2[0][2] += P1[2]; P2[1][2] += P1[2]; \ 321 } 322 323 /*! 324 \brief SUMA_FROM_BARYCENTRIC(u, v, p1, p2, p3, p) 325 change from barycentric coordinates. 326 */ 327 #define SUMA_FROM_BARYCENTRIC(u, v, p1, p2, p3, p){ \ 328 (p)[0] = (p1)[0] + (u) * ((p2)[0] - (p1)[0] ) + (v) * ((p3)[0] - (p1)[0]); \ 329 (p)[1] = (p1)[1] + (u) * ((p2)[1] - (p1)[1] ) + (v) * ((p3)[1] - (p1)[1]); \ 330 (p)[2] = (p1)[2] + (u) * ((p2)[2] - (p1)[2] ) + (v) * ((p3)[2] - (p1)[2]); \ 331 } 332 333 /*! 334 \brief calculate spherical coordinates from cartesian. 335 Assuming coords are centered on 0 0 0. 336 c XYZ coordinates in cartesian 337 s rtp Rho, theta (azimuth), phi (elevation) in spherical 338 \sa SUMA_Cart2Sph 339 \sa SUMA_SPH_2_CART 340 */ 341 #define SUMA_CART_2_SPH(c,s){\ 342 SUMA_NORM_VEC(c, 3, s[0]); \ 343 s[1] = atan2(c[1], c[0]); /* theta (azimuth) */ \ 344 s[2] = atan2(c[2],sqrt(c[0]*c[0]+c[1]*c[1])); /* phi (elevation)*/ \ 345 } 346 347 /*! 348 \brief calculate cartesian coordinates from spherical. 349 Assuming coords are centered on 0 0 0. 350 s rtp Rho, theta (azimuth), phi (elevation) in spherical 351 c XYZ coordinates in cartesian 352 \sa SUMA_CART_2_SPH 353 */ 354 #define SUMA_SPH_2_CART(s,c){\ 355 c[0] = s[0] * cos(s[2]) * cos(s[1]) ; \ 356 c[1] = s[0] * cos(s[2]) * sin(s[1]) ; \ 357 c[2] = s[0] * sin(s[2]); \ 358 } 359 /* largest absolute value */ 360 #define SUMA_LARG_ABS(a, b) ( ( fabs((double)(a)) > fabs((double)(b)) ) ? fabs((double)(a)) : fabs((double)(b)) ) 361 362 /*! the 1D index of element 363 [r][c] in a row major matrix 364 of nc columns */ 365 #define SRM(r,c,nc) ((nc)*(r)+(c)) 366 /*! the 1D index of element 367 [r][c] in a column major matrix 368 of nr rows */ 369 #define SCM(r,c,nr) ((r)+(nr)*(c)) 370 371 372 #define SUMA_WRAP_VALUE(v, min, max) \ 373 { \ 374 if (v > max) v = min; \ 375 else if (v < min) v = max; \ 376 } 377 378 #define SUMA_CLIP_VALUE(v, min, max) \ 379 { \ 380 if (v > max) v = max; \ 381 else if (v < min) v = min; \ 382 } 383 384 #define SUMA_CLIP_UB(v, max) \ 385 { \ 386 if (v > max) v = max; \ 387 } 388 389 #define SUMA_CLIP_LB(v, min) \ 390 { \ 391 if (v < min) v = min; \ 392 } 393 394 /*! determines the distance and side of the plane that a point is located 395 see also SUMA_Surf_Plane_Intersect 396 if Dist = 0, point on plane, if Dist > 0 point above plane (along normal), if Dist < 0 point is below plane 397 */ 398 #define SUMA_DIST_FROM_PLANE(P1, P2, P3, P, Dist){ \ 399 static float m_Eq[4]; \ 400 SUMA_Plane_Equation ( P1, P2, P3, m_Eq); \ 401 Dist = m_Eq[0] * P[0] + m_Eq[1] * P[1] + m_Eq[2] * P[2] + m_Eq[3] ; \ 402 } 403 404 #define SUMA_DIST_FROM_PLANE_EQ(Eq, P) ( Eq[0] * P[0] + Eq[1] * P[1] + \ 405 Eq[2] * P[2] + Eq[3] ) 406 407 408 #define SUMA_DIST_FROM_PLANE2(P1, P2, P3, P, Dist, InAirSpace, dot){ \ 409 static float m_Eq[4], m_P2[3], m_U[3], m_Un; \ 410 SUMA_Plane_Equation ( P1, P2, P3, m_Eq); \ 411 Dist = m_Eq[0] * P[0] + m_Eq[1] * P[1] + m_Eq[2] * P[2] + m_Eq[3] ; \ 412 m_P2[0] = m_Eq[0]+P[0]; m_P2[1] = m_Eq[1]+P[1]; m_P2[2] = m_Eq[2]+P[2]; \ 413 InAirSpace = SUMA_MT_isIntersect_Triangle (\ 414 P, m_P2, P1, P2, P3, NULL, NULL, NULL);\ 415 m_P2[0]=(P1[0]+P2[0]+P3[0])/3.0; \ 416 m_P2[1]=(P1[1]+P2[1]+P3[1])/3.0; \ 417 m_P2[2]=(P1[2]+P2[2]+P3[2])/3.0; \ 418 SUMA_UNIT_VEC(m_P2, P, m_U, m_Un); \ 419 SUMA_NORM(m_Un, m_Eq); \ 420 dot = m_Eq[0]*m_U[0]/m_Un + m_Eq[1]*m_U[1]/m_Un + m_Eq[2]*m_U[2]/m_Un;\ 421 } 422 423 /*! 424 Equation of a plane given its normal and a point 425 See also SUMA_Plane_Equation 426 */ 427 #define SUMA_PLANE_NORMAL_POINT(N, P, Eq) { \ 428 Eq[0] = N[0]; Eq[1] = N[1]; Eq[2] = N[2]; \ 429 Eq[3] = -(Eq[0]*P[0] + Eq[1]*P[1] + Eq[2]*P[2]); \ 430 } 431 432 #define SUMA_MARK_PLANE_NOT_SET(P) { (P)[0] = 0.0; (P)[1] = 0.0; \ 433 (P)[2] = 0.0; (P)[2] = 0.0; } 434 /* Is equation of plane marked as not set? */ 435 #define SUMA_PLANE_NOT_SET(P) (!(P) || \ 436 ((P)[0] == 0.0 && (P)[1] == 0.0 && \ 437 (P)[2] == 0.0 && (P)[2] == 0.0)) 438 #define SUMA_PLANE_IS_SET(P) (!(SUMA_PLANE_NOT_SET((P)))) 439 440 /*! 441 Shift a plane to pass through P 442 See also SUMA_Plane_Equation 443 */ 444 #define SUMA_SHIFT_PLANE_TO_P(Eq, P) { \ 445 Eq[3] = -(Eq[0]*P[0] + Eq[1]*P[1] + Eq[2]*P[2]); \ 446 } 447 448 /*! 449 determines the bounding box for a triangle 450 */ 451 #define SUMA_TRIANGLE_BOUNDING_BOX(n1, n2, n3, min_v, max_v){ \ 452 min_v[0] = SUMA_MIN_PAIR( (n1)[0], (n2)[0]); min_v[0] = SUMA_MIN_PAIR( (n3)[0], min_v[0]); \ 453 min_v[1] = SUMA_MIN_PAIR( (n1)[1], (n2)[1]); min_v[1] = SUMA_MIN_PAIR( (n3)[1], min_v[1]); \ 454 min_v[2] = SUMA_MIN_PAIR( (n1)[2], (n2)[2]); min_v[2] = SUMA_MIN_PAIR( (n3)[2], min_v[2]); \ 455 max_v[0] = SUMA_MAX_PAIR( (n1)[0], (n2)[0]); max_v[0] = SUMA_MAX_PAIR( (n3)[0], max_v[0]); \ 456 max_v[1] = SUMA_MAX_PAIR( (n1)[1], (n2)[1]); max_v[1] = SUMA_MAX_PAIR( (n3)[1], max_v[1]); \ 457 max_v[2] = SUMA_MAX_PAIR( (n1)[2], (n2)[2]); max_v[2] = SUMA_MAX_PAIR( (n3)[2], max_v[2]); \ 458 } 459 460 /*! 461 \brief find out if a point p on the line formed by points p0 and p1 is between p0 and p1 462 ans = SUMA_IS_POINT_IN_SEGMENT(p, p0, p1); 463 if ans then p is between p0 and p1 464 p, p0, p1 (float *) xyz of points 465 466 - NOTE: macro does not check that three points are colinear (as they should be)! 467 */ 468 #define SUMA_IS_POINT_IN_SEGMENT(p, p0, p1) ( ( \ 469 ( \ 470 (p[0] > p0[0] && p[0] < p1[0]) || \ 471 (p[0] < p0[0] && p[0] > p1[0]) || \ 472 (SUMA_ABS(p[0] - p0[0]) < 0.00000001 || \ 473 SUMA_ABS(p[0] - p1[0]) < 0.00000001 ) \ 474 )\ 475 && \ 476 ( \ 477 (p[1] > p0[1] && p[1] < p1[1]) || \ 478 (p[1] < p0[1] && p[1] > p1[1]) || \ 479 (SUMA_ABS(p[1] - p0[1]) < 0.00000001 || \ 480 SUMA_ABS(p[1] - p1[1]) < 0.00000001 ) \ 481 )\ 482 && \ 483 ( \ 484 (p[2] > p0[2] && p[2] < p1[2]) || \ 485 (p[2] < p0[2] && p[2] > p1[2]) || \ 486 (SUMA_ABS(p[2] - p0[2]) < 0.00000001 || \ 487 SUMA_ABS(p[2] - p1[2]) < 0.00000001 ) \ 488 )\ 489 ) ? 1 : 0 ) 490 491 /*! 492 \brief Project point N onto plane Eq 493 */ 494 #define SUMA_PROJECT_ONTO_PLANE(Eq,N,P) {\ 495 double cn = (Eq)[0]*(Eq)[0]+(Eq)[1]*(Eq)[1]+(Eq)[2]*(Eq)[2]; \ 496 if (cn != 0.0f) { \ 497 cn = ((Eq)[0]*(N)[0]+(Eq)[1]*(N)[1]+(Eq)[2]*(N)[2]+(Eq)[3])/cn; \ 498 (P)[0] = (N)[0] - (Eq)[0] * cn; \ 499 (P)[1] = (N)[1] - (Eq)[1] * cn; \ 500 (P)[2] = (N)[2] - (Eq)[2] * cn; \ 501 } else { \ 502 (P)[0] = (P)[1] = (P)[2] = 0.0; \ 503 } \ 504 } 505 506 /*! 507 \brief Project point C onto line with direction U (unit vector)and passing 508 through the origin, projection is in P. It is like using 509 SUMA_PROJECT_C_ONTO_AB with A being the origin, and B being U 510 */ 511 #define SUMA_PROJECT_ONTO_DIR(U,C,P) {\ 512 double f; \ 513 f = (U[0]*C[0]+U[1]*C[1]+U[2]*C[2]);\ 514 P[0] = f*U[0]; P[1]=f*U[1]; P[2]=f*U[2];\ 515 } 516 517 /*! 518 \brief Project point C onto vector AB, projection is in P, fractional distance 519 of projection from A is in f 520 */ 521 #define SUMA_PROJECT_C_ONTO_AB(C, A, B, P, f) {\ 522 static double AB[3], AC[3], AP[3], ABdAC, ABdAB; \ 523 AB[0]=B[0]-A[0]; AB[1]=B[1]-A[1]; AB[2]=B[2]-A[2]; \ 524 AC[0]=C[0]-A[0]; AC[1]=C[1]-A[1]; AC[2]=C[2]-A[2]; \ 525 ABdAC = (AB[0]*AC[0]+AB[1]*AC[1]+AB[2]*AC[2]);\ 526 ABdAB = (AB[0]*AB[0]+AB[1]*AB[1]+AB[2]*AB[2]);\ 527 if (ABdAB) {\ 528 f = ABdAC/ABdAB; \ 529 AP[0] = f*AB[0]; AP[1]=f*AB[1]; AP[2]=f*AB[2];\ 530 P[0] = AP[0]+A[0]; P[1] = AP[1]+A[1]; P[2] = AP[2]+A[2]; \ 531 } else {\ 532 f=0; P[0] = P[1] = P[2] = 0.0; \ 533 }\ 534 } 535 /* Same as SUMA_PROJECT_C_ONTO_AB but external variable to contain vector AB.*/ 536 #define SUMA_PROJECT_C_ONTO_ABv(C, A, B, P, f, AB) {\ 537 static double AC[3], AP[3], ABdAC, ABdAB; \ 538 AB[0]=B[0]-A[0]; AB[1]=B[1]-A[1]; AB[2]=B[2]-A[2]; \ 539 AC[0]=C[0]-A[0]; AC[1]=C[1]-A[1]; AC[2]=C[2]-A[2]; \ 540 ABdAC = (AB[0]*AC[0]+AB[1]*AC[1]+AB[2]*AC[2]);\ 541 ABdAB = (AB[0]*AB[0]+AB[1]*AB[1]+AB[2]*AB[2]);\ 542 if (ABdAB) {\ 543 f = ABdAC/ABdAB; \ 544 AP[0] = f*AB[0]; AP[1]=f*AB[1]; AP[2]=f*AB[2];\ 545 P[0] = AP[0]+A[0]; P[1] = AP[1]+A[1]; P[2] = AP[2]+A[2]; \ 546 } else {\ 547 f=0; P[0] = P[1] = P[2] = 0.0; \ 548 }\ 549 } 550 551 552 /*! 553 \brief intersection of line defined by point N and direction U 554 and plane Eq 555 P is the resultant intersection point and ispar = 1 if 556 line is parallel to the plane. 557 558 MACRO IS UNUSED, NEEDS TESTING 559 */ 560 #define SUMA_LINE_PLANE_INTERSECT(U,N,Eq,P,ispar) {\ 561 double m_dot, m_fr; \ 562 SUMA_S_Warn("UNTESTED-2"); \ 563 dot = SUMA_MT_DOT(U,Eq); \ 564 ispar = 0; \ 565 if (dot != 0.0f) m_fr = -(Eq[0]*N[0]+Eq[1]*N[1]+Eq[2]*N[2]+Eq[3])/dot; \ 566 else ispar = 1; \ 567 P[0] = N[0]+m_fr*U[0]; P[1] = N[1]+m_fr*U[1]; P[2] = N[2]+m_fr*U[2]; \ 568 } 569 570 /*! 571 \brief Reset quaternion values for default 3D view 572 */ 573 #define SUMA_HOME_QUAT(vw,qq) { \ 574 float m_a[3]; \ 575 switch (vw) {\ 576 case SUMA_2D_Z0: \ 577 m_a[0] = 1.0; m_a[1] = 0.0; m_a[2] = 0.0; \ 578 axis_to_quat(m_a, 0, qq); \ 579 break; \ 580 case SUMA_2D_Z0L: \ 581 m_a[0] = 1.0; m_a[1] = 0.0; m_a[2] = 0.0; \ 582 axis_to_quat(m_a, SUMA_PI, qq); \ 583 break; \ 584 case SUMA_3D_Z0: \ 585 m_a[0] = 1.0; m_a[1] = 0.0; m_a[2] = 0.0; \ 586 axis_to_quat(m_a, 0, qq); \ 587 break; \ 588 default: \ 589 case SUMA_3D: \ 590 *((qq) ) = 0.252199; *((qq)+1) = -0.129341; \ 591 *((qq)+2) = -0.016295; *((qq)+3) = 0.958854; \ 592 break; \ 593 } \ 594 } 595 596 /*! 597 \brief Intersection of a segment with a plane 598 \param p1 x y z of point 1 599 \param p2 x y z of point 2 600 \param Eq equation of plane 601 \param Hit: On exit: 1 segment interects plane 602 0 segment does not intersect plane 603 \param pinter: On exit: x y z of intersection point (0,0,0) if no intersection 604 605 Note: Macro is not efficient to use when intersection if between entire surface 606 and plane. For that use SUMA_Surf_Plane_Intersect 607 */ 608 #define SUMA_SEGMENT_PLANE_INTERSECT(p1, p2, Eq, Hit, pinter) { \ 609 double m_p1, m_p2, m_u; \ 610 m_p1 = Eq[0] * p1[0] + Eq[1] * p1[1] + Eq[2] * p1[2] + Eq[3]; \ 611 m_p2 = Eq[0] * p2[0] + Eq[1] * p2[1] + Eq[2] * p2[2] + Eq[3]; \ 612 /* fprintf(SUMA_STDERR,"m_p1=%f, m_p2 = %f\n", m_p1, m_p2); */\ 613 if ((SUMA_SIGN(m_p1)) != (SUMA_SIGN(m_p2)) ) { \ 614 Hit = 1; \ 615 m_u = -m_p1 / (m_p2 - m_p1); \ 616 pinter[0] = p1[0] + m_u * ( p2[0] - p1[0] ); \ 617 pinter[1] = p1[1] + m_u * ( p2[1] - p1[1] ); \ 618 pinter[2] = p1[2] + m_u * ( p2[2] - p1[2] ); \ 619 } else { \ 620 Hit = 0; \ 621 pinter[0] = pinter[1] = pinter[2] = 0.0; \ 622 } \ 623 } 624 625 #define SUMA_SET_GL_RENDER_MODE(m_PolyMode) \ 626 { \ 627 switch (m_PolyMode) { \ 628 case SRM_Fill: \ 629 glEnable (GL_POLYGON_OFFSET_FILL); \ 630 glPolygonOffset(1.0, 1.0); \ 631 /* Polygon offset is needed to deal with rendering 632 lines that are coplanar with filled polygons. 633 Without polygon offset, lines can get stripy under 634 certain angles, quite ugly. The effect is known as 635 stitching, bleeding, or Z fighting. 636 http://www.opengl.org/resources/faq/technical/polygonoffset.htm */ \ 637 glPolygonMode(GL_FRONT_AND_BACK, GL_FILL); \ 638 break; \ 639 case SRM_Line: \ 640 glDisable (GL_POLYGON_OFFSET_FILL); \ 641 glPolygonMode(GL_FRONT_AND_BACK, GL_LINE); \ 642 break; \ 643 case SRM_Points: \ 644 glDisable (GL_POLYGON_OFFSET_FILL); \ 645 glPolygonMode(GL_FRONT_AND_BACK, GL_POINT); \ 646 break; \ 647 case SRM_ViewerDefault: \ 648 break; \ 649 case SRM_Hide: \ 650 break; \ 651 default: \ 652 fprintf (SUMA_STDERR, "Error %s: Wrong Rendering Mode.\n", FuncName); \ 653 break; \ 654 } \ 655 } 656 657 #define SUMA_SET_GL_RENDER_MODE_TRACK(m_PolyMode, m_ust) \ 658 { float m_fv[4]={0.0, 0.0, 0.0, 0.0};\ 659 DList *m_st=m_ust; \ 660 switch (m_PolyMode) { \ 661 case SRM_Fill: \ 662 SUMA_GLStateTrack("set", &m_st, FuncName, \ 663 "GL_POLYGON_OFFSET_FILL", (void *)1); \ 664 /* glEnable (GL_POLYGON_OFFSET_FILL); */\ 665 m_fv[0] = 1.0; m_fv[1] = 1.0; \ 666 SUMA_GLStateTrack("set", &m_st, FuncName, \ 667 "glPolygonOffset", (void *)m_fv); \ 668 /* glPolygonOffset(1.0, 1.0); */\ 669 /* Polygon offset is needed to deal with rendering 670 lines that are coplanar with filled polygons. 671 Without polygon offset, lines can get stripy under 672 certain angles, quite ugly. The effect is known as 673 stitching, bleeding, or Z fighting. 674 http://www.opengl.org/resources/faq/technical/polygonoffset.htm */ \ 675 m_fv[0] = GL_FRONT_AND_BACK; m_fv[1] = GL_FILL; \ 676 SUMA_GLStateTrack("set", &m_st, FuncName, \ 677 "glPolygonMode", (void *)m_fv); \ 678 /* glPolygonMode(GL_FRONT_AND_BACK, GL_FILL); */ \ 679 break; \ 680 case SRM_Line: \ 681 SUMA_GLStateTrack("set", &m_st, FuncName, \ 682 "GL_POLYGON_OFFSET_FILL", (void *)0); \ 683 /* glDisable (GL_POLYGON_OFFSET_FILL); */ \ 684 m_fv[0] = GL_FRONT_AND_BACK; m_fv[1] = GL_LINE; \ 685 SUMA_GLStateTrack("set", &m_st, FuncName, \ 686 "glPolygonMode", (void *)m_fv); \ 687 /* glPolygonMode(GL_FRONT_AND_BACK, GL_LINE); */ \ 688 break; \ 689 case SRM_Points: \ 690 SUMA_GLStateTrack("set", &m_st, FuncName, \ 691 "GL_POLYGON_OFFSET_FILL", (void *)0); \ 692 /* glDisable (GL_POLYGON_OFFSET_FILL); */ \ 693 m_fv[0] = GL_FRONT_AND_BACK; m_fv[1] = GL_POINT; \ 694 SUMA_GLStateTrack("set", &m_st, FuncName, \ 695 "glPolygonMode", (void *)m_fv); \ 696 /* glPolygonMode(GL_FRONT_AND_BACK, GL_POINT); */ \ 697 break; \ 698 case SRM_ViewerDefault: \ 699 break; \ 700 case SRM_Hide: \ 701 break; \ 702 default: \ 703 fprintf (SUMA_STDERR, "Error %s: Wrong Rendering Mode.\n", FuncName); \ 704 break; \ 705 } \ 706 } 707 708 #define SUMA_SET_GL_TRANS_MODE(m_TransMode, m_ust, btp) \ 709 { DList *m_st=m_ust; \ 710 if (m_TransMode == STM_0) { \ 711 /* classic rendering, no polygon stippling*/ \ 712 SUMA_LH("Setting GL_POLYGON_STIPPLE off");\ 713 SUMA_GLStateTrack("set", &m_st, FuncName, \ 714 "GL_POLYGON_STIPPLE", (void *)0); \ 715 /* glDisable(GL_POLYGON_STIPPLE); */ \ 716 if (0) {/* ZSS March 11 2014 */ \ 717 SUMA_LH("Setting GL_DEPTH_TEST on"); \ 718 SUMA_GLStateTrack("set", &m_st, FuncName, \ 719 "GL_DEPTH_TEST", (void *)1); \ 720 /* glEnable(GL_DEPTH_TEST); */ } \ 721 } else { \ 722 SUMA_LH("Setting GL_POLYGON_STIPPLE on");\ 723 SUMA_GLStateTrack("set", &m_st, FuncName, \ 724 "GL_POLYGON_STIPPLE", (void *)1); \ 725 glPolygonStipple(SUMA_StippleMask_shift(m_TransMode, btp)); \ 726 if (0) { SUMA_LH("Setting GL_DEPTH_TEST off"); \ 727 SUMA_GLStateTrack("set", &m_st, \ 728 FuncName, "GL_DEPTH_TEST", (void *)0); }\ 729 } \ 730 } 731 732 #define SUMA_GLX_BUF_SWAP(sv) {\ 733 if ((sv)->X->DOUBLEBUFFER) { \ 734 if ((sv)->DO_PickMode) { \ 735 glFlush(); \ 736 } else { \ 737 glXSwapBuffers((sv)->X->DPY, XtWindow((sv)->X->GLXAREA));\ 738 } \ 739 } else { \ 740 glFlush();\ 741 } \ 742 } 743 744 #define SUMA_NEW_MASKSTATE() (SUMAg_CF->X ? (++SUMAg_CF->X->MaskStateID):0) 745 746 /*! 747 \brief Get the pointer for the last visited viewer 748 */ 749 #define SUMA_LAST_VIEWER (SUMAg_CF->PointerLastInViewer >=0 && SUMAg_N_SVv >1 ? \ 750 &(SUMAg_SVv[SUMAg_CF->PointerLastInViewer]) : \ 751 &(SUMAg_SVv[0])) 752 753 /*! 754 \brief SO->Show is not quite not the end of the story 755 */ 756 #define SO_SHOWING(SO,sv) ( SO->Show && SO->PolyMode != SRM_Hide && (SO->PolyMode != SRM_ViewerDefault || sv->PolyMode != SRM_Hide) && (SO->TransMode != STM_ViewerDefault || sv->TransMode != STM_16)) 757 758 #define MDO_SHOWING(MDO,sv) ( (MDO) && (MDO)->SO && (sv) && \ 759 (SO_SHOWING((MDO)->SO,sv)) ) 760 /*! 761 \brief set polymode 762 */ 763 #define SUMA_SET_SO_POLYMODE(SO,i){ \ 764 if (i < 0 || i >= SRM_N_RenderModes) { SO->PolyMode = SRM_ViewerDefault; }\ 765 else { SO->PolyMode = i; } \ 766 if (SO->PolyMode == SRM_Hide) { SO->Show = NOPE; } \ 767 else { SO->Show = YUP; }\ 768 } 769 770 /*! 771 \brief calculates the average 'radius' of a surface. 772 avg(dist(node_i,center)); 773 */ 774 #define SUMA_SO_RADIUS(SO, r){ \ 775 int m_i, m_i3; \ 776 float m_dx, m_dy, m_dz; \ 777 r = 0.0; \ 778 for (m_i=0; m_i<SO->N_Node; ++m_i) { \ 779 m_i3 = 3 * m_i; \ 780 m_dx = SO->NodeList[m_i3 ] - SO->Center[0]; \ 781 m_dy = SO->NodeList[m_i3+1] - SO->Center[1]; \ 782 m_dz = SO->NodeList[m_i3+2] - SO->Center[2]; \ 783 r += sqrt( (m_dx * m_dx) + (m_dy * m_dy) + (m_dz * m_dz)); \ 784 } \ 785 r /= (float)SO->N_Node; \ 786 } 787 788 /*! 789 \brief calculates the farthest point on a surface from its center. 790 max(dist(node_i,center)); 791 */ 792 #define SUMA_SO_MAX_MIN_DIST(SO, DD, NN, dd, nn){ \ 793 int m_i, m_i3=0; \ 794 float m_dx, m_dy, m_dz; double m_Dm, m_d, m_dm; \ 795 dd = 0.0; DD = 0.0; m_dm=11111111110.0; m_Dm=0.0; nn = -1; NN = -1;\ 796 for (m_i=0; m_i<SO->N_Node; ++m_i) { \ 797 m_dx = SO->NodeList[m_i3++] - SO->Center[0]; \ 798 m_dy = SO->NodeList[m_i3++] - SO->Center[1]; \ 799 m_dz = SO->NodeList[m_i3++] - SO->Center[2]; \ 800 m_d = (m_dx * m_dx) + (m_dy * m_dy) + (m_dz * m_dz); \ 801 if (m_d > m_Dm) { m_Dm = m_d; NN = m_i; } \ 802 else if (m_d < m_dm) { m_dm = m_d; nn = m_i; } \ 803 } \ 804 if (NN > -1) DD = sqrt(m_Dm); \ 805 if (nn > -1) dd = sqrt(m_dm); \ 806 } 807 808 /*! 809 A macro to recalculate a surface object's normals 810 */ 811 #define SUMA_RECOMPUTE_NORMALS(SO){ \ 812 SUMA_SURF_NORM m_SN; \ 813 if (SO->NodeNormList) SUMA_free(SO->NodeNormList); \ 814 SO->NodeNormList = NULL; \ 815 if (SO->FaceNormList) SUMA_free(SO->FaceNormList); \ 816 SO->FaceNormList = NULL; \ 817 set_surf_norm_quiet(1); \ 818 m_SN = SUMA_SurfNorm(SO->NodeList, SO->N_Node, \ 819 SO->FaceSetList, SO->N_FaceSet ); \ 820 SO->NodeNormList = m_SN.NodeNormList; \ 821 SO->FaceNormList = m_SN.FaceNormList; \ 822 SO->glar_NodeNormList = (GLfloat *) SO->NodeNormList; \ 823 /* just copy the pointer, not the data */\ 824 SO->glar_FaceNormList = (GLfloat *) SO->FaceNormList; \ 825 } 826 827 /*! 828 A macro to recalculate a surface object's normals 829 */ 830 #define SUMA_RECOMPUTE_POLYGON_AREAS(SO){ \ 831 if (SO->PolyArea) SUMA_free(SO->PolyArea); SO->PolyArea=NULL; \ 832 if (!SUMA_SurfaceMetrics_eng(SO, "PolyArea", NULL, 0, \ 833 SUMAg_CF->DsetList)) { \ 834 SUMA_S_Err("Failed to recompute poly areas"); \ 835 } \ 836 }\ 837 838 #define SUMA_RECOMPUTE_NORMALS_and_AREAS(SO) {\ 839 SUMA_RECOMPUTE_NORMALS(SO); \ 840 SUMA_RECOMPUTE_POLYGON_AREAS(SO); \ 841 } 842 843 /*! 844 Pause prompt, stdin 845 */ 846 #define SUMA_PAUSE_PROMPT_STDIN(s) { int m_jnk; fprintf(SUMA_STDOUT,"Pausing: %s ...", s); fflush(SUMA_STDOUT); m_jnk = getchar(); fprintf(SUMA_STDOUT,"\n"); fflush(SUMA_STDOUT);} 847 848 #define SUMA_PAUSE_PROMPT(s) { \ 849 XtAppContext m_app; \ 850 Widget m_w; \ 851 XEvent m_ev; \ 852 XtInputMask m_pp; \ 853 int m_ii=0, m_wkill=1; \ 854 if (SUMAg_SVv && SUMAg_SVv[0].X->TOPLEVEL){\ 855 m_w = SUMAg_SVv[0].X->TOPLEVEL; \ 856 m_wkill = 0; \ 857 } else { \ 858 m_w = XtOpenApplication(&m_app, FuncName, \ 859 NULL, 0, &m_ii, NULL, \ 860 SUMA_get_fallbackResources(), \ 861 topLevelShellWidgetClass, NULL, 0); \ 862 } \ 863 if (m_w) { \ 864 SUMA_PauseForUser(m_w, s, \ 865 SWP_POINTER_OFF, &m_app, 0, -1.0); \ 866 if (m_wkill) { \ 867 /* This is set when the top level widget is created */ \ 868 /* and there is no XtAppMainLoop to process events and */\ 869 /* ensure that the destroyed widget gets off of the screen */\ 870 /* so you'll need to do the event processing yourself */\ 871 /* or you'll end up with ghostly windows if events come in */\ 872 /* too rapidly */ \ 873 /* A simpler call to XtAppNextEvent(app, &ev); was enough */ \ 874 /* to cause the widget to go away, but not when events are */ \ 875 /* piling up. Not too surprising since events are not getting */ \ 876 /* processed. The while statement seems to do the trick */ \ 877 while ((m_pp = XtAppPending(m_app))) { \ 878 XtAppProcessEvent(m_app, m_pp); \ 879 } \ 880 XtDestroyWidget(m_w); \ 881 }; \ 882 } else {\ 883 SUMA_PAUSE_PROMPT_STDIN(s); \ 884 } \ 885 } 886 887 /*! 888 A macro to recalculate a surface's center and its bounding box 889 */ 890 #define SUMA_DIM_CENTER(SO){ \ 891 SUMA_MIN_MAX_SUM_VECMAT_COL (SO->NodeList, SO->N_Node, SO->NodeDim, SO->MinDims, SO->MaxDims, SO->Center); \ 892 SO->Center[0] /= SO->N_Node; \ 893 SO->Center[1] /= SO->N_Node; \ 894 SO->Center[2] /= SO->N_Node; \ 895 SUMA_MIN_VEC (SO->MinDims, 3, SO->aMinDims ); \ 896 SUMA_MAX_VEC (SO->MaxDims, 3, SO->aMaxDims); \ 897 } 898 899 /*! 900 A macro to scale the surface so that the largest bounding box size is SZ 901 */ 902 #define SUMA_LARGEST_SIZE_SCALE(SO, SZ){ \ 903 float m_lg , m_lg1 , m_lg2 ; \ 904 double scl = 1.0; \ 905 int m_i_mx, m_i, m_itop; \ 906 m_i_mx = 0; \ 907 m_lg = SO->MaxDims[0] - SO->Center[0]; \ 908 m_lg1 = SO->MaxDims[1] - SO->Center[1], \ 909 m_lg2 = SO->MaxDims[2] - SO->Center[2]; \ 910 m_itop = SO->NodeDim * SO->N_Node; \ 911 if (m_lg < m_lg1) { m_lg = m_lg1; m_i_mx = 1; } \ 912 if (m_lg < m_lg2) { m_lg = m_lg2; m_i_mx = 2; } \ 913 if (m_lg > 0.0) { \ 914 scl = (double)SZ / 2.0 / (double)m_lg; \ 915 for (m_i=0; m_i < m_itop; ++m_i) { SO->NodeList[m_i] = (float)(scl*(double)SO->NodeList[m_i]); } \ 916 for (m_i=0; m_i < 3; ++m_i) { \ 917 SO->MinDims[m_i] = (float)(scl*(double)SO->MinDims[m_i]); \ 918 SO->MaxDims[m_i] = (float)(scl*(double)SO->MaxDims[m_i]); \ 919 SO->Center[m_i] = (float)(scl*(double)SO->Center[m_i]); \ 920 } \ 921 SO->aMinDims = (float)(scl*(double)SO->aMinDims); \ 922 SO->aMaxDims = (float)(scl*(double)SO->aMaxDims); \ 923 }\ 924 } 925 926 /*! calculate the centroid of a triangular faceset */ 927 #define SUMA_FACE_CENTROID(SO, ifc, c){ \ 928 static int m_n1, m_n2, m_n3; \ 929 m_n1 = SO->FaceSetList[3*ifc]; m_n2 = SO->FaceSetList[3*ifc+1]; m_n3 = SO->FaceSetList[3*ifc+2]; \ 930 c[0] = (SO->NodeList[3*m_n1] + SO->NodeList[3*m_n2] + SO->NodeList[3*m_n3] )/3; \ 931 c[1] = (SO->NodeList[3*m_n1+1] + SO->NodeList[3*m_n2+1] + SO->NodeList[3*m_n3+1])/3; \ 932 c[2] = (SO->NodeList[3*m_n1+2] + SO->NodeList[3*m_n2+2] + SO->NodeList[3*m_n3+2])/3; \ 933 } 934 935 /*! 936 A macro to find the third node forming a triangle 937 */ 938 #define SUMA_THIRD_TRIANGLE_NODE(n1,n2,t,facelist,n3){ \ 939 static int m_t3; \ 940 m_t3 = 3 * t; \ 941 n3 = -1; \ 942 do { \ 943 if (facelist[m_t3] != n1 && facelist[m_t3] != n2) n3 = facelist[m_t3]; \ 944 else ++m_t3; \ 945 } while (n3 < 0); \ 946 } 947 948 /*! 949 A macro to find the two other nodes forming a triangle 950 */ 951 #define SUMA_TWO_OTHER_TRIANGLE_NODES(n1,t,facelist,n2,n3){ \ 952 static int m_t3; \ 953 m_t3 = 3 * t; \ 954 n2 = n3 = -1; \ 955 if (facelist[m_t3 ] == n1) { \ 956 n2 = facelist[m_t3+1]; \ 957 n3 = facelist[m_t3+2]; \ 958 } else if (facelist[m_t3+1] == n1) { \ 959 n2 = facelist[m_t3 ]; \ 960 n3 = facelist[m_t3+2]; \ 961 } else if (facelist[m_t3+2] == n1) { \ 962 n2 = facelist[m_t3 ]; \ 963 n3 = facelist[m_t3+1]; \ 964 } \ 965 } 966 967 968 969 /*! 970 A macro version of SUMA_FindEdge 971 Use function for robust error checking. 972 Here you are on your own. 973 974 you should make sure n1 < n2 975 you should initialize iseg to -1 before calling the macro 976 977 */ 978 /* NEW VERSION: Bug found in previous one (_OLD), 979 it is possible that m_eloc is -1 at initialization: m_eloc = m_EL->ELloc[m_n1]; 980 see below for fix 981 */ 982 #define SUMA_FIND_EDGE(m_EL, m_n1, m_n2, m_iseg) \ 983 { int m_eloc = m_EL->ELloc[m_n1]; \ 984 if (m_eloc >= 0) { /* It is possible that m_n1 has m_EL->ELloc[m_n1] = -1, this happens when m_n1 is surrounded by nodes of lower indices */\ 985 for (m_eloc=m_EL->ELloc[m_n1]; m_eloc<m_EL->N_EL; ++m_eloc) { \ 986 /* fprintf(stderr, "%d/%d\n", m_eloc, m_EL->N_EL); */\ 987 if (m_EL->EL[m_eloc][0] != m_n1) break; \ 988 if (m_EL->EL[m_eloc][1] == m_n2) m_iseg = m_eloc; \ 989 } \ 990 } \ 991 } 992 #define SUMA_FIND_EDGE_OLD(m_EL, m_n1, m_n2, m_iseg) \ 993 { int m_eloc ; \ 994 m_eloc = m_EL->ELloc[m_n1]; \ 995 do { \ 996 if (m_EL->EL[m_eloc][1] == m_n2) m_iseg = m_eloc; \ 997 ++m_eloc; \ 998 } while (m_EL->EL[m_eloc][0] == m_n1 && m_eloc < m_EL->N_EL); \ 999 } 1000 1001 /*! 1002 Macro that sets Ntail to be the index of the last node of the last ROIdatum in a 1003 DrawnROI. Ntail is -1 if the macro fails 1004 */ 1005 #define SUMA_DRAWN_ROI_TAIL_NODE(D_ROI, Ntail) \ 1006 {\ 1007 DListElmt *m_Tail=NULL; \ 1008 SUMA_ROI_DATUM *m_ROId=NULL; \ 1009 Ntail = -1; \ 1010 m_Tail = dlist_tail(D_ROI->ROIstrokelist); \ 1011 if (m_Tail) { \ 1012 m_ROId = (SUMA_ROI_DATUM *)m_Tail->data; \ 1013 if (m_ROId->N_n) Ntail = m_ROId->nPath[m_ROId->N_n-1]; \ 1014 } \ 1015 } 1016 1017 /*! 1018 Macro that sets Nhead to be the index of the 1st node of the 1st ROIdatum in a 1019 DrawnROI. Nhead is -1 if the macro fails 1020 */ 1021 #define SUMA_DRAWN_ROI_HEAD_NODE(D_ROI, Nhead) \ 1022 {\ 1023 DListElmt *m_Head=NULL; \ 1024 SUMA_ROI_DATUM *m_ROId=NULL; \ 1025 Nhead = -1; \ 1026 m_Head = dlist_head(D_ROI->ROIstrokelist); \ 1027 if (m_Head) { \ 1028 m_ROId = (SUMA_ROI_DATUM *)m_Head->data; \ 1029 if (m_ROId->N_n) Nhead = m_ROId->nPath[0]; \ 1030 } \ 1031 } 1032 1033 /*! is 3x4 xform an indentity matrix? */ 1034 #define SUMA_IS_XFORM_IDENTITY(Xform) \ 1035 ( Xform[0][0] == 1.0 && Xform[1][1] == 1.0 && Xform[2][2] == 1.0 && \ 1036 Xform[0][3] == 0.0 && Xform[1][3] == 0.0 && Xform[2][3] == 0.0 && \ 1037 Xform[0][1] == 0.0 && Xform[0][2] == 0.0 && \ 1038 Xform[1][0] == 0.0 && Xform[1][2] == 0.0 && \ 1039 Xform[2][0] == 0.0 && Xform[2][1] == 0.0 ) ? 1:0 1040 1041 /* definitions for SUMA_MT_intersect */ 1042 #define SUMA_MT_CROSS(m_MTCR_dest,m_MTCR_v1,m_MTCR_v2) \ 1043 m_MTCR_dest[0]=m_MTCR_v1[1]*m_MTCR_v2[2]-m_MTCR_v1[2]*m_MTCR_v2[1]; \ 1044 m_MTCR_dest[1]=m_MTCR_v1[2]*m_MTCR_v2[0]-m_MTCR_v1[0]*m_MTCR_v2[2]; \ 1045 m_MTCR_dest[2]=m_MTCR_v1[0]*m_MTCR_v2[1]-m_MTCR_v1[1]*m_MTCR_v2[0]; 1046 #define SUMA_MT_DOT(m_MTDOT_v1,m_MTDOT_v2) (m_MTDOT_v1[0]*m_MTDOT_v2[0]+m_MTDOT_v1[1]*m_MTDOT_v2[1]+m_MTDOT_v1[2]*m_MTDOT_v2[2]) 1047 #define SUMA_MT_SUB(m_MTSUB_dest,m_MTSUB_v1,m_MTSUB_v2) \ 1048 m_MTSUB_dest[0]=m_MTSUB_v1[0]-m_MTSUB_v2[0]; \ 1049 m_MTSUB_dest[1]=m_MTSUB_v1[1]-m_MTSUB_v2[1]; \ 1050 m_MTSUB_dest[2]=m_MTSUB_v1[2]-m_MTSUB_v2[2]; 1051 #define SUMA_NORM(m_NORM_dest, m_NORM_v1) \ 1052 m_NORM_dest= sqrt(m_NORM_v1[0]*m_NORM_v1[0]+\ 1053 m_NORM_v1[1]*m_NORM_v1[1]+\ 1054 m_NORM_v1[2]*m_NORM_v1[2]); 1055 #define SUMA_NORM3(m_NORM_v1) ( sqrt(*((m_NORM_v1) ) * *((m_NORM_v1) )+\ 1056 *((m_NORM_v1)+1) * *((m_NORM_v1)+1)+\ 1057 *((m_NORM_v1)+2) * *((m_NORM_v1)+2) ) ) 1058 #define SUMA_TRI_AREA(m_TRIAREA_n0,m_TRIAREA_n1,m_TRIAREA_n2, m_TRIAREA_A) {\ 1059 float m_TRIAREA_dv[3], m_TRIAREA_dw[3], m_TRIAREA_cross[3]; \ 1060 SUMA_MT_SUB (m_TRIAREA_dv, m_TRIAREA_n1, m_TRIAREA_n0); \ 1061 SUMA_MT_SUB (m_TRIAREA_dw, m_TRIAREA_n2, m_TRIAREA_n0); \ 1062 SUMA_MT_CROSS(m_TRIAREA_cross,m_TRIAREA_dv,m_TRIAREA_dw); \ 1063 SUMA_NORM(m_TRIAREA_A, m_TRIAREA_cross); \ 1064 m_TRIAREA_A *= 0.5; \ 1065 } 1066 1067 /*! 1068 XYZ: vector of coordinates 1069 nrm: unit normal vector of axis of rotaion 1070 phi: angle of rotation, counterclockwise, in radians 1071 XYZr: vector to contain rotated coordinates 1072 Equation from: http://mathworld.wolfram.com/RotationFormula.html 1073 Note that this equation subtracts the last term, the cross product. 1074 The website equation adds the last term. 1075 This change permits counterclockwise angle of rotation and the proper axis of rotation. 1076 See also RotateAboutAxis.m 1077 */ 1078 1079 #define SUMA_ROTATE_ABOUT_AXIS(xyz, nrm, phi, xyzr) { \ 1080 double m_cop, m_sip, m_dop, m_cro[3], m_1cop; \ 1081 m_cop = cos(phi); m_1cop = 1 - m_cop; m_sip = sin(phi); \ 1082 SUMA_MT_CROSS(m_cro, xyz, nrm); \ 1083 m_dop = SUMA_MT_DOT(nrm,xyz); \ 1084 xyzr[0] = xyz[0]*m_cop + nrm[0] * m_dop * m_1cop - m_cro[0] * m_sip; \ 1085 xyzr[1] = xyz[1]*m_cop + nrm[1] * m_dop * m_1cop - m_cro[1] * m_sip; \ 1086 xyzr[2] = xyz[2]*m_cop + nrm[2] * m_dop * m_1cop - m_cro[2] * m_sip; \ 1087 } 1088 1089 /*! 1090 SUMA_3D_Rotation_Matrix(P1, P2, M, alpha, nrm) 1091 Generate the rotation matrix for moving from P1 to P2. 1092 For rotation about a specified origin and direction. 1093 Rotation through angle alpha about an axis m_u. Alpha is always positive and counterclockwise. 1094 P1, P2: intial and final vector location 1095 m_alpha: angle between P1 and P2 1096 m_u: axis of rotation and cross product of P1, P2 1097 mn_u: magnitude of axis of rotation needed to normalize the axis of rotation 1098 m_M: 3x3 rotation matrix 1099 Equations based on matlab script AxisRotate3D.m 1100 */ 1101 1102 #define SUMA_3D_Rotation_Matrix(P1, P2, M, m_alpha, m_u) { \ 1103 static double m_x, m_y, m_z, m_mag_u; \ 1104 static double m_xy_vera, m_xz_vera, m_yz_vera; \ 1105 static double m_cosa, m_sina, m_vera; \ 1106 SUMA_ANGLE_DIST_NC(P2 , P1, m_alpha, m_u); \ 1107 m_mag_u = sqrt( m_u[0]*m_u[0] + m_u[1]*m_u[1] + m_u[2]*m_u[2] ); \ 1108 if( m_mag_u > 0.000000001 ) { m_u[0] = m_u[0]/m_mag_u ; m_u[1] = m_u[1]/m_mag_u ; m_u[2] = m_u[2]/m_mag_u ;}; \ 1109 m_x = m_u[0]; m_y = m_u[1]; m_z = m_u[2]; \ 1110 m_cosa = cos(m_alpha); m_sina = sin(m_alpha); m_vera = 1 - m_cosa; \ 1111 m_xy_vera = m_x * m_y * m_vera; m_xz_vera = m_x * m_z * m_vera; m_yz_vera = m_y * m_z * m_vera;\ 1112 \ 1113 M[0][0] = m_cosa + m_x * m_x * m_vera; \ 1114 M[0][1] = m_xy_vera - m_z * m_sina; \ 1115 M[0][2] = m_xz_vera + m_y * m_sina; \ 1116 M[1][0] = m_xy_vera + m_z * m_sina; \ 1117 M[1][1] = m_cosa + m_y * m_y * m_vera; \ 1118 M[1][2] = m_yz_vera - m_x * m_sina; \ 1119 M[2][0] = m_xz_vera - m_y*m_sina; \ 1120 M[2][1] = m_yz_vera + m_x*m_sina; \ 1121 M[2][2] = m_cosa + m_z * m_z * m_vera; \ 1122 } 1123 1124 /*! 1125 SUMA_3Dax_Rotation_Matrix(Ax, Alpha, M) 1126 Generate the rotation matrix for rotating an angle Alpha 1127 around unit length axis Ax. 1128 M: 3x3 rotation matrix 1129 */ 1130 1131 #define SUMA_3Dax_Rotation_Matrix(AX, m_alpha, M) { \ 1132 static double m_x, m_y, m_z, m_mag_u; \ 1133 static double m_xy_vera, m_xz_vera, m_yz_vera; \ 1134 static double m_cosa, m_sina, m_vera; \ 1135 m_x = AX[0]; m_y = AX[1]; m_z = AX[2]; \ 1136 m_cosa = cos(m_alpha); m_sina = sin(m_alpha); m_vera = 1 - m_cosa; \ 1137 m_xy_vera = m_x * m_y * m_vera; \ 1138 m_xz_vera = m_x * m_z * m_vera; m_yz_vera = m_y * m_z * m_vera;\ 1139 \ 1140 M[0][0] = m_cosa + m_x * m_x * m_vera; \ 1141 M[0][1] = m_xy_vera - m_z * m_sina; \ 1142 M[0][2] = m_xz_vera + m_y * m_sina; \ 1143 M[1][0] = m_xy_vera + m_z * m_sina; \ 1144 M[1][1] = m_cosa + m_y * m_y * m_vera; \ 1145 M[1][2] = m_yz_vera - m_x * m_sina; \ 1146 M[2][0] = m_xz_vera - m_y*m_sina; \ 1147 M[2][1] = m_yz_vera + m_x*m_sina; \ 1148 M[2][2] = m_cosa + m_z * m_z * m_vera; \ 1149 } 1150 1151 1152 /*! 1153 SUMA_ANGLE_DIST(p2,p1,cent,a,nrm) 1154 SUMA_ANGLE_DIST_NC(p2,p1,a, nrm) 1155 calculate angular distance between two points 1156 on a sphere. 1157 For the _NC version, the sphere is centered on 1158 0 0 0 and has a radius of 1 1159 p1 and p2 are the XYZ of the two points 1160 cent is the center of the sphere 1161 and 'a' is the angle in radians between them. 1162 m_cr contains the cross product p1 cross p2 1163 Tx to tip from JHU's Applied Physics Laboratory web page 1164 */ 1165 #define SUMA_ANGLE_DIST_NC(m_p2,m_p1,a, m_cr) \ 1166 {\ 1167 double m_p2r[3]; \ 1168 static double m_mag; \ 1169 SUMA_MT_CROSS(m_cr, m_p1, m_p2); \ 1170 a = atan2(sqrt(m_cr[0]*m_cr[0]+m_cr[1]*m_cr[1]+m_cr[2]*m_cr[2]),SUMA_MT_DOT(m_p2, m_p1)); \ 1171 } 1172 1173 #define SUMA_ANGLE_DIST(p2,p1,cent,a,m_cr) \ 1174 {\ 1175 double m_p1[3],m_p2[3],m_np1, m_np2; \ 1176 SUMA_MT_SUB(m_p1, p1, cent); SUMA_NORM(m_np1, m_p1); if (m_np1) { m_p1[0] /= m_np1; m_p1[1] /= m_np1; m_p1[2] /= m_np1; }\ 1177 SUMA_MT_SUB(m_p2, p2, cent); SUMA_NORM(m_np2, m_p2); if (m_np2) { m_p2[0] /= m_np2; m_p2[1] /= m_np2; m_p2[2] /= m_np2; }\ 1178 SUMA_MT_CROSS(m_cr, m_p2, m_p1); \ 1179 a = atan2(sqrt(m_cr[0]*m_cr[0]+m_cr[1]*m_cr[1]+m_cr[2]*m_cr[2]),SUMA_MT_DOT(m_p2, m_p1)); \ 1180 } 1181 1182 /*! 1183 SUMA_SINC(alpha) 1184 calculate the sinc function. 1185 alpha is an angle in radians 1186 sinc is the value of the function at alpha 1187 */ 1188 #define SUMA_SINC(alpha, sinc) \ 1189 {\ 1190 if( alpha < 0.0000001 ) { sinc = 1.0; } \ 1191 else { sinc = sin(alpha) / alpha ; } \ 1192 } 1193 1194 /*! 1195 \brief SUMA_EULER_SO (SO, eu) 1196 computes the euler number = N - E + F 1197 eu = SO->N_Node - SO->EL->N_Distinct_Edges + SO->N_FaceSet 1198 eu = 2 for closed surfaces 1199 -1000 --> NULL SO 1200 -1001 --> NULL SO->EL 1201 */ 1202 1203 #define SUMA_EULER_SO(SO, eu) { \ 1204 if (!SO) { eu = -1000; } \ 1205 else if (!SO->EL) { eu = -1001; } \ 1206 else eu = SO->N_Node - SO->EL->N_Distinct_Edges + SO->N_FaceSet; \ 1207 } 1208 1209 /*! 1210 1211 \brief SUMA_IS_IN_VEC(vec, nel, val, loc); 1212 \param vec pointer to vector of at least nel elements 1213 \param nel (int) number of elements to check in vec 1214 \param val value to look for in vec 1215 \param loc (int) index in vec where val was found. -1 if nothing was found 1216 */ 1217 #define SUMA_IS_IN_VEC(m_vec, m_nel, m_val, m_loc) { \ 1218 int m_i=0;\ 1219 m_loc = -1;\ 1220 while (m_i < m_nel) { \ 1221 if (m_vec[m_i] == m_val) { \ 1222 m_loc = m_i; \ 1223 m_i = m_nel; \ 1224 } else { \ 1225 ++ m_i; \ 1226 } \ 1227 } \ 1228 } 1229 1230 #define SUMA_DBG_IN_NOTIFY(m_fname) { \ 1231 int m_i;\ 1232 ++SUMAg_CF->InOut_Level; \ 1233 for (m_i=0; m_i < SUMAg_CF->InOut_Level; ++m_i) fprintf (SUMA_STDERR," ");\ 1234 fprintf (SUMA_STDERR,"--dbg: Entered %s (lvl %d).\n", m_fname, SUMAg_CF->InOut_Level); \ 1235 } 1236 1237 #define SUMA_DBG_OUT_NOTIFY(m_fname) { \ 1238 int m_i;\ 1239 for (m_i=0; m_i < SUMAg_CF->InOut_Level; ++m_i) fprintf (SUMA_STDERR," ");\ 1240 fprintf (SUMA_STDERR,"--dbg: Left %s (lvl %d).\n", m_fname, SUMAg_CF->InOut_Level); \ 1241 --SUMAg_CF->InOut_Level; \ 1242 } 1243 1244 1245 /*! \def SUMA_REPORT_WICH_WIDGET_SV(m_w) 1246 \brief SUMA_REPORT_WICH_WIDGET_SV macro for determining what type of widget m_w is and which sv it belongs to 1247 */ 1248 #define SUMA_REPORT_WICH_WIDGET_SV(m_w) {\ 1249 int m_i = 0; \ 1250 SUMA_SurfaceViewer *m_sv = NULL; \ 1251 int m_svi = -1; \ 1252 while (m_i < SUMA_MAX_SURF_VIEWERS) { \ 1253 /* the series of if statements should be else-ifs but I left them like that for debugging purposes */ \ 1254 if (SUMAg_SVv[m_i].X->GLXAREA == m_w) { \ 1255 fprintf (SUMA_STDERR,"SUMA_REPORT_WICH_WIDGET_SV: %p is GLXAREA widget in Surface Viewer %d.\n", m_w, m_i); \ 1256 m_svi = m_i; \ 1257 } \ 1258 if (SUMAg_SVv[m_i].X->TOPLEVEL == m_w) { \ 1259 fprintf (SUMA_STDERR,"SUMA_REPORT_WICH_WIDGET_SV: %p is TOPLEVEL widget in Surface Viewer %d.\n", m_w, m_i); \ 1260 m_svi = m_i; \ 1261 } \ 1262 if (SUMAg_SVv[m_i].X->FORM == m_w) { \ 1263 fprintf (SUMA_STDERR,"SUMA_REPORT_WICH_WIDGET_SV: %p is FORM widget in Surface Viewer %d.\n", m_w, m_i); \ 1264 m_svi = m_i; \ 1265 } \ 1266 if (SUMAg_SVv[m_i].X->FRAME == m_w) { \ 1267 fprintf (SUMA_STDERR,"SUMA_REPORT_WICH_WIDGET_SV: %p is FRAME widget in Surface Viewer %d.\n", m_w, m_i); \ 1268 m_svi = m_i; \ 1269 } \ 1270 ++m_i; \ 1271 } \ 1272 } 1273 1274 1275 /*! \def SUMA_ANY_WIDGET2SV(m_w, m_sv, m_svi) 1276 \brief SUMA_ANY_WIDGET2SV macro for determining the SurfaceViewer structure containing any of the following widgets: GLXAREA, TOPLEVEL, FORM, FRAME. The macro searches all the SurfaceViewer structures in SUMAg_SVv. 1277 1278 m_w the widget in question 1279 m_sv a pointer to SUMA_SurfaceViewer structure. This pointer is NULL if no matching SurfaceViewer structure is found in SUMAg_SVv 1280 m_svi the index of m_sv in SUMAg_SVv vector of Surface Viewer structures. m_sv = &(SUMAg_SVv[m_svi]). -1 if no match was found 1281 */ 1282 1283 #define SUMA_ANY_WIDGET2SV(m_w, m_sv, m_svi) {\ 1284 int m_i = 0; \ 1285 m_sv = NULL; \ 1286 m_svi = -1; \ 1287 while (m_i < SUMA_MAX_SURF_VIEWERS) { \ 1288 if (SUMAg_SVv[m_i].X->GLXAREA == m_w || SUMAg_SVv[m_i].X->TOPLEVEL == m_w || SUMAg_SVv[m_i].X->FORM == m_w || SUMAg_SVv[m_i].X->FRAME == m_w) { \ 1289 m_svi = m_i; \ 1290 m_sv = &(SUMAg_SVv[m_i]); \ 1291 m_i = SUMA_MAX_SURF_VIEWERS; \ 1292 } else { \ 1293 ++m_i; \ 1294 } \ 1295 } \ 1296 } 1297 1298 /*! \def SUMA_GLXAREA_WIDGET2SV(m_w, m_sv, m_svi) 1299 \brief SUMA_GLXAREA_WIDGET2SV macro for determining the SurfaceViewer structure containing the widget: GLXAREA. The macro searches all the SurfaceViewer structures in SUMAg_SVv. 1300 1301 m_w the widget in question 1302 m_sv a pointer to SUMA_SurfaceViewer structure. This pointer is NULL if no matching SurfaceViewer structure is found in SUMAg_SVv 1303 m_svi the index of m_sv in SUMAg_SVv vector of Surface Viewer structures. m_sv = &(SUMAg_SVv[m_svi]). -1 if no match was found 1304 */ 1305 1306 #define SUMA_GLXAREA_WIDGET2SV(m_w, m_sv, m_svi) {\ 1307 int m_i = 0; \ 1308 m_sv = NULL; \ 1309 m_svi = -1; \ 1310 while (m_i < SUMA_MAX_SURF_VIEWERS) { \ 1311 if (SUMAg_SVv[m_i].X->GLXAREA == m_w) { \ 1312 m_svi = m_i; \ 1313 m_sv = &(SUMAg_SVv[m_i]); \ 1314 m_i = SUMA_MAX_SURF_VIEWERS; \ 1315 } else { \ 1316 ++m_i; \ 1317 } \ 1318 } \ 1319 } 1320 1321 /*! 1322 Find the closest node in SO to coordinate xyz 1323 distance^2 is stored in d and nc is the index of the 1324 the closest node 1325 */ 1326 #define SUMA_CLOSEST_NODE(SO, xyz, nc, d) { \ 1327 double m_dxyz = 1023734552736672366372.0; \ 1328 float *m_p; \ 1329 int m_n=0; \ 1330 nc = -1; \ 1331 d = m_dxyz; \ 1332 for (m_n=0; m_n<SO->N_Node; ++m_n) { \ 1333 m_p = &(SO->NodeList[SO->NodeDim*m_n]); \ 1334 SUMA_SEG_LENGTH_SQ(m_p, xyz, m_dxyz); \ 1335 if (d < m_dxyz) { \ 1336 d = m_dxyz; nc = m_n; \ 1337 } \ 1338 } \ 1339 } 1340 #define SUMA_CLOSEST_NODE_VEC(SO, xyz, nc, d, dv) { \ 1341 double m_dxyz = 1023734552736672366372.0; \ 1342 double m_dv[3];\ 1343 float *m_p; \ 1344 int m_n=0; \ 1345 nc = -1; \ 1346 d = m_dxyz; \ 1347 for (m_n=0; m_n<SO->N_Node; ++m_n) { \ 1348 m_p = &(SO->NodeList[SO->NodeDim*m_n]); \ 1349 m_dv[0] = m_p[0] - xyz[0]; \ 1350 m_dv[1] = m_p[1] - xyz[1]; \ 1351 m_dv[2] = m_p[2] - xyz[2]; \ 1352 m_dxyz = m_dv[0]*m_dv[0]+\ 1353 m_dv[1]*m_dv[1]+\ 1354 m_dv[2]*m_dv[2]; \ 1355 if (d > m_dxyz) { \ 1356 d=m_dxyz; nc = m_n; \ 1357 dv[0]=m_dv[0]; dv[1]=m_dv[1]; dv[2]=m_dv[2];\ 1358 } \ 1359 } \ 1360 } 1361 1362 /* Some macros are taken from DSP_in_C examples in 1363 C Language Algorithms for Digital Signal Processing 1364 by 1365 Bruce Kimball, Paul Embree and Bruce Kimble 1366 1991, Prentice Hall 1367 */ 1368 1369 /*! \def SUMA_SEG_LENGTH(a,b, dist) 1370 \brief SUMA_SEG_LENGTH macro for a segment's length 1371 a pointer to xyz coordinates 1372 b pointer to xyz coordinates 1373 dist (float) sqrt( (m_b[0] - m_a[0]) * (m_b[0] - m_a[0]) 1374 +(m_b[1] - m_a[1]) * (m_b[1] - m_a[1]) 1375 +(m_b[2] - m_a[2]) * (m_b[2] - m_a[2]) ); 1376 */ 1377 #define SUMA_SEG_LENGTH(m_a, m_b, m_dist) { \ 1378 m_dist = sqrt( (m_b[0] - m_a[0]) * (m_b[0] - m_a[0]) \ 1379 +(m_b[1] - m_a[1]) * (m_b[1] - m_a[1]) \ 1380 +(m_b[2] - m_a[2]) * (m_b[2] - m_a[2]) ); \ 1381 } 1382 1383 /*! \def SUMA_SEG_LENGTH_SQ(a,b, dist) 1384 \brief SUMA_SEG_LENGTH_SQ macro for a segment's squared length 1385 a pointer to xyz coordinates 1386 b pointer to xyz coordinates 1387 dist (float) ( (m_b[0] - m_a[0]) * (m_b[0] - m_a[0]) 1388 +(m_b[1] - m_a[1]) * (m_b[1] - m_a[1]) 1389 +(m_b[2] - m_a[2]) * (m_b[2] - m_a[2]) ); 1390 */ 1391 #define SUMA_SEG_LENGTH_SQ(m_a, m_b, m_dist) { \ 1392 m_dist = (m_b[0] - m_a[0]) * (m_b[0] - m_a[0]) \ 1393 +(m_b[1] - m_a[1]) * (m_b[1] - m_a[1]) \ 1394 +(m_b[2] - m_a[2]) * (m_b[2] - m_a[2]) ; \ 1395 } 1396 1397 /*! \def SUMA_NORM_VEC(a,nel, norm) 1398 \brief SUMA_NORM_VEC macro for vectors's norm (sqrt of sum of squares) 1399 a pointer to vector 1400 nel number of elements in vector 1401 norm (float) norm of a 1402 */ 1403 #define SUMA_NORM_VEC(a,nel,norm) { \ 1404 int m_I; \ 1405 norm = 0.0; \ 1406 for (m_I = 0; m_I < nel; m_I++) { \ 1407 norm += a[m_I]*a[m_I]; \ 1408 } \ 1409 norm = sqrt(norm); \ 1410 } 1411 1412 /*! 1413 Make vector have unit norm 1414 */ 1415 #define SUMA_UNITIZE_VEC(a,nel) {\ 1416 double m_norm=0.0; \ 1417 int m_I=0; \ 1418 for (m_I = 0; m_I < nel; m_I++) { \ 1419 m_norm += a[m_I]*a[m_I]; \ 1420 } \ 1421 m_norm = sqrt(m_norm); \ 1422 if (m_norm != 0.0f) { \ 1423 for (m_I = 0; m_I < nel; m_I++) {\ 1424 a[m_I] /= m_norm; \ 1425 } \ 1426 } \ 1427 } 1428 1429 /*! \def SUMA_MIN_VEC(a,nel, amin) 1430 \brief SUMA_MIN_VEC macro for minimum 1431 a pointer to vector 1432 nel number of elements in vector 1433 amin minimum of a (make sure types of a and amin match) 1434 */ 1435 #define SUMA_MIN_VEC(a,nel,amin) { \ 1436 int m_I; \ 1437 amin = a[0]; \ 1438 for (m_I = 1; m_I < nel; m_I++) { \ 1439 if (a[m_I] < amin) amin = a[m_I]; \ 1440 } \ 1441 } 1442 1443 /*! \def SUMA_MIN_LOC_VEC(a,nel, amin, minloc) 1444 \brief SUMA_MIN_VEC macro for minimum identification and location 1445 a pointer to vector 1446 nel (int) number of elements in vector 1447 amin minimum of a (make sure types of a and amin match) 1448 minloc (int) index into a where the minimum was found 1449 */ 1450 #define SUMA_MIN_LOC_VEC(a,nel,amin, minloc) { \ 1451 int m_I; \ 1452 amin = a[0]; \ 1453 minloc = 0; \ 1454 for (m_I = 1; m_I < nel; m_I++) { \ 1455 if (a[m_I] < amin) { \ 1456 amin = a[m_I]; \ 1457 minloc = m_I; \ 1458 } \ 1459 } \ 1460 } 1461 1462 /*! \def SUMA_MAX_VEC(a,nel,amax) 1463 \brief SUMA_MAX_VEC macro for minimum 1464 a pointer to vector 1465 nel number of elements in vector 1466 amax maximum of a (make sure types of a and amax match) 1467 */ 1468 #define SUMA_MAX_VEC(a,nel,amax) { \ 1469 int m_I; \ 1470 amax = a[0]; \ 1471 for (m_I = 1; m_I < nel; m_I++) { \ 1472 if (a[m_I] > amax) amax = a[m_I]; \ 1473 } \ 1474 } 1475 1476 /* 1477 Replace value in a with forward sum of w values 1478 if (avg) then divide resultant vector by w 1479 Only nel - w +1 values are set in the output 1480 a[k] = S(a[k]+...a[k+w-1]) 1481 */ 1482 #define SUMA_MOVING_SUM(a,nel,w,avg, seal) { \ 1483 int m_I=w; double mbuf0, mbuf1; \ 1484 if (w < nel) { \ 1485 mbuf0 = a[0]; for (m_I=1; m_I<w; ++m_I) a[0]+=a[m_I]; \ 1486 m_I = 1; \ 1487 while (m_I+w <= nel) { \ 1488 mbuf1 = a[m_I]; \ 1489 a[m_I] = a[m_I-1] - mbuf0 +a[m_I+w-1]; mbuf0 = mbuf1; \ 1490 ++m_I; \ 1491 } \ 1492 if (avg) { for (m_I=0; m_I<nel-w+1; ++m_I) a[m_I] /= (double)w; } \ 1493 if (seal) for (m_I=nel-w+1; m_I<nel; ++m_I) a[m_I] = 0; \ 1494 }\ 1495 } 1496 1497 #define SUMA_MEAN_VEC(a,nel,amean,nozero) { \ 1498 int m_I, m_c=0; double m_mean=0.0; \ 1499 amean = 0; \ 1500 if (nozero) {\ 1501 for (m_I = 0; m_I < nel; m_I++) { \ 1502 if (a[m_I]!=0.0f) { ++m_c; m_mean += a[m_I]; }\ 1503 }\ 1504 if (m_c) amean = m_mean/m_c; \ 1505 } else { \ 1506 for (m_I = 0; m_I < nel; m_I++) { \ 1507 m_mean += a[m_I]; \ 1508 }\ 1509 if (nel) amean = m_mean/nel; \ 1510 }\ 1511 } 1512 1513 #define SUMA_MEAN_STD_VEC(a,nel,amean, astd, nozero) { \ 1514 int m_I, m_c=0; double m_std=0.0, m_mmmm=0.0, dd=0.0;\ 1515 SUMA_MEAN_VEC(a,nel,m_mmmm,nozero); \ 1516 amean=m_mmmm; astd=0; \ 1517 if (nozero) {\ 1518 for (m_I = 0; m_I < nel; m_I++) { \ 1519 if (a[m_I]!=0.0f) { ++m_c; dd=(a[m_I]-m_mmmm); m_std += dd*dd; }\ 1520 }\ 1521 if (m_c > 1) astd = sqrt(m_std/(m_c-1)); \ 1522 } else { \ 1523 for (m_I = 0; m_I < nel; m_I++) { \ 1524 dd=(a[m_I]-m_mmmm); m_std += dd*dd; \ 1525 }\ 1526 if (nel > 1) astd = sqrt(m_std/(nel-1)); \ 1527 }\ 1528 } 1529 1530 1531 1532 /*! \def SUMA_MIN_MAX_VEC(a,nel,amin, amax, aminloc, amaxloc) 1533 \brief SUMA_MIN_MAX_VEC macro for minimum and maximum 1534 a pointer to vector 1535 nel number of elements in vector 1536 amin minimum of a (make sure types of a and amin match) 1537 amax maximum of a (make sure types of a and amax match) 1538 aminloc index where minimum is found 1539 amaxloc index where maximum is found 1540 */ 1541 #define SUMA_MIN_MAX_VEC(a,nel, amin, amax, aminloc, amaxloc) { \ 1542 int m_I; \ 1543 amaxloc = 0; \ 1544 amax = a[0]; \ 1545 aminloc = 0; \ 1546 amin = a[0];\ 1547 for (m_I = 1; m_I < nel; m_I++) { \ 1548 if (a[m_I] > amax) { amax = a[m_I]; amaxloc = m_I; } \ 1549 else { if (a[m_I] < amin) { amin = a[m_I]; aminloc = m_I; } } \ 1550 } \ 1551 } 1552 #define SUMA_MIN_MAX_CVEC(a,nel, amin, amax, aminloc, amaxloc, phase) { \ 1553 int m_I; \ 1554 double m_abs; \ 1555 amaxloc = 0; \ 1556 aminloc = 0; \ 1557 if (phase) { \ 1558 amax = CARG(a[0]); \ 1559 amin = CARG(a[0]); \ 1560 } else { \ 1561 amax = CABS(a[0]); \ 1562 amin = CABS(a[0]); \ 1563 } \ 1564 for (m_I = 1; m_I < nel; m_I++) { \ 1565 if (phase) { \ 1566 m_abs = CARG((a[m_I])); \ 1567 } else { \ 1568 m_abs = CABS((a[m_I])); \ 1569 }\ 1570 if (m_abs > amax) { amax = m_abs; amaxloc = m_I; } \ 1571 else { if (m_abs < amin) { amin = m_abs; aminloc = m_I; } } \ 1572 } \ 1573 } 1574 1575 #define SUMA_MIN_MAX_VEC_STRIDE(a,nel, amin, amax, aminloc, amaxloc, stride) { \ 1576 int m_I; \ 1577 amaxloc = 0; \ 1578 amax = a[0]; \ 1579 aminloc = 0; \ 1580 amin = a[0];\ 1581 for (m_I = stride; m_I < nel; m_I = m_I+stride) { \ 1582 if (a[m_I] > amax) { amax = a[m_I]; amaxloc = m_I; } \ 1583 else { if (a[m_I] < amin) { amin = a[m_I]; aminloc = m_I; } } \ 1584 } \ 1585 } 1586 1587 #define SUMA_MIN_MAX_CVEC_STRIDE(a,nel, amin, amax, aminloc, amaxloc, stride, phase) { \ 1588 int m_I; \ 1589 double m_abs; \ 1590 complex m_c; \ 1591 amaxloc = 0; \ 1592 aminloc = 0; \ 1593 if (phase) { \ 1594 amax = CARG(a[0]); \ 1595 amin = CARG(a[0]); \ 1596 } else { \ 1597 amax = CABS(a[0]); \ 1598 amin = CABS(a[0]); \ 1599 } \ 1600 for (m_I = stride; m_I < nel; m_I = m_I+stride) { \ 1601 if (phase) { \ 1602 m_abs = CARG((a[m_I])); \ 1603 } else { \ 1604 m_abs = CABS((a[m_I])); \ 1605 }\ 1606 if ( m_abs > amax) { amax = m_abs; amaxloc = m_I; } \ 1607 else { if (m_abs < amin) { amin = m_abs; aminloc = m_I; } } \ 1608 } \ 1609 } 1610 1611 /* 1612 SUMA_ADD_VEC macro: 1613 1614 ADDS TWO VECTORS (a,b) POINT BY POINT (PROMOTING AS REQUIRED) AND 1615 PUTS THE RESULT IN THE c VECTOR (DEMOTING IF REQUIRED). 1616 1617 SUMA_ADD_VEC(a,b,c,len,typea,typeb,typec) 1618 1619 a pointer to first vector. 1620 b pointer to second vector. 1621 c pointer to result vector. 1622 len length of vectors (integer). 1623 typea legal C type describing the type of a data. 1624 typeb legal C type describing the type of b data. 1625 typec legal C type describing the type of c data. 1626 */ 1627 1628 #ifndef SUMA_ADD_VEC 1629 #define SUMA_ADD_VEC(a,b,c,len,typea,typeb,typec) { \ 1630 typea *_PTA = a; \ 1631 typeb *_PTB = b; \ 1632 typec *_PTC = c; \ 1633 int m_IX; \ 1634 for(m_IX = 0 ; m_IX < (len) ; m_IX++) \ 1635 *_PTC++ = (typec)((*_PTA++) + (*_PTB++)); \ 1636 } 1637 1638 #endif 1639 1640 /* 1641 SUMA_SUB_VEC macro: 1642 1643 SUBTRACTS TWO VECTORS (a,b) POINT BY POINT (PROMOTING AS REQUIRED) AND 1644 PUTS THE RESULT IN THE c VECTOR (DEMOTING IF REQUIRED). 1645 1646 SUMA_SUB_VEC(a,b,c,len,typea,typeb,typec) 1647 1648 a pointer to first vector. 1649 b pointer to second vector. 1650 c pointer to result vector. 1651 len length of vectors (integer). 1652 typea legal C type describing the type of a data. 1653 typeb legal C type describing the type of b data. 1654 typec legal C type describing the type of c data. 1655 1656 */ 1657 1658 #ifndef SUMA_SUB_VEC 1659 #define SUMA_SUB_VEC(a,b,c,len,typea,typeb,typec) { \ 1660 typea *_PTA = a; \ 1661 typeb *_PTB = b; \ 1662 typec *_PTC = c; \ 1663 int m_IX; \ 1664 for(m_IX = 0 ; m_IX < (len) ; m_IX++) \ 1665 *_PTC++ = (typec)((*_PTA++) - (*_PTB++)); \ 1666 } 1667 #endif 1668 /* 1669 SUMA_MULT_VEC macro: 1670 1671 MULTIPLIES TWO VECTORS (a,b) POINT BY POINT (PROMOTING AS REQUIRED) AND 1672 PUTS THE RESULT IN THE c VECTOR (DEMOTING IF REQUIRED). 1673 1674 SUMA_MULT_VEC(a,b,c,len,typea,typeb,typec) 1675 1676 a pointer to first vector. 1677 b pointer to second vector. 1678 c pointer to result vector. 1679 len length of vectors (integer). 1680 typea legal C type describing the type of a data. 1681 typeb legal C type describing the type of b data. 1682 typec legal C type describing the type of c data. 1683 1684 WARNING: The input data vectors are not cast to the type of c. 1685 This means that at least one of the input types must 1686 be able to represent the individual products without 1687 overflow. 1688 1689 */ 1690 1691 #ifndef SUMA_MULT_VEC 1692 #define SUMA_MULT_VEC(a,b,c,len,typea,typeb,typec) { \ 1693 typea *_PTA = a; \ 1694 typeb *_PTB = b; \ 1695 typec *_PTC = c; \ 1696 int m_IX; \ 1697 for(m_IX = 0 ; m_IX < (len) ; m_IX++) \ 1698 *_PTC++ = (typec)((*_PTA++) * (*_PTB++)); \ 1699 } 1700 #endif 1701 1702 /* 1703 SUMA_SUM_VEC macro: 1704 1705 FORMS THE SUM THE VECTOR a AND PUT THE RESULT IN THE 1706 PREVIOUSLY DEFINED VARIABLE s. 1707 1708 SUMA_SUM_VEC(a,s,len,typea) 1709 1710 a pointer to first vector. 1711 s variable used to store result (not a pointer). 1712 len length of vector (integer). 1713 typea legal C type describing the type of a data. 1714 1715 */ 1716 1717 #ifndef SUMA_SUM_VEC 1718 #define SUMA_SUM_VEC(a,s,len,typea) { \ 1719 typea *_PTA = a; \ 1720 int m_IX; \ 1721 s = (*_PTA++); \ 1722 for(m_IX = 1 ; m_IX < (len) ; m_IX++) \ 1723 s += (*_PTA++); \ 1724 } 1725 #endif 1726 1727 /* 1728 SUMA_SCALE_VEC macro: 1729 1730 SCALES AND/OR CONVERTS (PROMOTES OR DEMOTES) THE INPUT VECTOR a 1731 (of typea) AND COPIES THE SCALED VECTOR INTO ANOTHER VECTOR b 1732 (of typeb). 1733 1734 SUMA_SCALE_VEC(a,b,s,len,typea,typeb) 1735 1736 a pointer to input vector. 1737 b pointer to output vector. 1738 s variable used to scale output vector (not a pointer). 1739 len length of vectors (integer). 1740 typea legal C type describing the type of a data. 1741 typeb legal C type describing the type of b data. 1742 1743 */ 1744 1745 #ifndef SUMA_SCALE_VEC 1746 #define SUMA_SCALE_VEC(a,b,s,len,typea,typeb) { \ 1747 typea *_PTA = (typea *)a; \ 1748 typeb *_PTB = (typeb *)b; \ 1749 int m_IX; \ 1750 for(m_IX = 0 ; m_IX < (len) ; m_IX++) \ 1751 *(_PTB)++ = (typeb)(s * (*(_PTA)++)); \ 1752 } 1753 #endif 1754 1755 /* SUMA_EXTRACT_VEC macro: 1756 1757 copies a[ind] into b where ind is a vector of indices 1758 1759 SUMA_EXTRACT_VEC (a,b,ind,len) 1760 1761 a pointer to input vector 1762 b pointer to output vector 1763 ind vector containing the indices of the values in a to copy to b 1764 len the length of ind (which is that of b too) 1765 1766 DO NOT SEND TWO POINTERS THAT POINT TO THE SAME LOCATION ! 1767 */ 1768 1769 1770 #define SUMA_EXTRACT_VEC(a,b,ind,len,typea,typeb) { \ 1771 typea *_PTA = (typea *)a; \ 1772 typeb *_PTB = (typeb *)b; \ 1773 int m_IX; \ 1774 for(m_IX = 0 ; m_IX < (len) ; m_IX++) \ 1775 _PTB[m_IX] = _PTA[ind[m_IX]]; \ 1776 } 1777 1778 1779 /* SUMA_CAT_VEC macro : 1780 concatenates two vectors and b together such that a = [a b]; 1781 1782 SUMA_CAT_VEC(a,b, catata, lenb,typea,typeb) 1783 a pointer to first vector 1784 b pointer to second vector 1785 catata index indicating where b is to be concatenated 1786 to a at. if catata = 6, then a[6] = b[0] 1787 and a[7] = b[1] etc ... 1788 make sure that you allocate enough space for a 1789 lenb number of values in b 1790 1791 typea legal C type describing the type of a data. 1792 typeb legal C type describing the type of b data. 1793 1794 */ 1795 1796 #define SUMA_CAT_VEC(a,b, catata, lenb,typea,typeb) { \ 1797 typea *_PTA = (typea *)a; \ 1798 typeb *_PTB = (typeb *)b; \ 1799 int m_IX; \ 1800 _PTA = _PTA + catata; \ 1801 for(m_IX = 0 ; m_IX < (lenb) ; m_IX++) \ 1802 *(_PTA)++ = (typea)(*(_PTB)++); \ 1803 } 1804 1805 1806 /*! 1807 SUMA_GET_MAT_ROW MACRO FOR GETTING A ROW FROM A MATRIX: 1808 1809 SUMA_GET_MAT_ROW(a,b,row,cols,typea,typeb) 1810 1811 a pointer to input 2D matrix. 1812 b pointer to resultant 1D vector. 1813 row index of row to extract from a 1814 cols number of columns in matrix a 1815 typea legal C type describing the type of a 1816 typeb legal C type describing the type of b 1817 1818 */ 1819 #define SUMA_GET_MAT_ROW(a,b,row,cols,typea,typeb) { \ 1820 typea **_AMX = (typea **)a; \ 1821 typeb *_PTB = (typeb *)b; \ 1822 typea *_PTA; \ 1823 int _JX; \ 1824 _PTA = _AMX[row]; \ 1825 for(_JX = 0 ; _JX < cols ; _JX++) \ 1826 *_PTB++ = (typeb) (*_PTA++); \ 1827 } 1828 1829 /*! 1830 SUMA_GET_MAT_COL MACRO FOR GETTING A COLUMN FROM A MATRIX: 1831 1832 SUMA_GET_MAT_COL(a,b, col, rows,typea,typeb) 1833 1834 a pointer to input 2D matrix. 1835 b pointer to resultant 1D vector. 1836 col index of column in matrix a to extract to b 1837 rows number of rows to in a 1838 typea legal C type describing the type of a 1839 typeb legal C type describing the type of b 1840 1841 */ 1842 1843 #define SUMA_GET_MAT_COL(a,b,col , rows, typea,typeb) { \ 1844 typea **_AMX = (typea **)a; \ 1845 typeb *_PTB = (typeb *)b; \ 1846 typea *_PTA; \ 1847 int m_IX,_JX; \ 1848 for(m_IX = 0 ; m_IX < rows ; m_IX++) { \ 1849 _PTA = _AMX[m_IX ] ; \ 1850 for(_JX = 0 ; _JX < col ; _JX++) \ 1851 _PTA++; \ 1852 *_PTB++ = (typeb)(*_PTA++); \ 1853 } \ 1854 } 1855 1856 1857 /*! \def SUMA_MIN_MAT_COL(a, rows, cols, amin) 1858 \brief SUMA_MIN_MAT_COL macro for minimum of each column in a matrix 1859 a pointer to matrix (**) 1860 rows number of rows 1861 cols number of cols 1862 amin minimum of each column in a (make sure types of a and amin match) 1863 */ 1864 #define SUMA_MIN_MAT_COL(a, rows, cols, amin) { \ 1865 int m_IX, _JX; \ 1866 for (m_IX = 0; m_IX < cols ; m_IX++) { \ 1867 amin[m_IX]=a[0][m_IX]; \ 1868 for (_JX = 1 ; _JX < rows ; _JX++) \ 1869 if (a[_JX][m_IX] < amin[m_IX]) amin[m_IX] = a[_JX][m_IX];\ 1870 }\ 1871 } 1872 1873 1874 /*! \def SUMA_MAX_MAT_COL(a, rows, cols, amax) 1875 \brief SUMA_MAX_MAT_COL macro for maximum of each column in a matrix 1876 a pointer to matrix (**) 1877 rows number of rows 1878 cols number of cols 1879 amax maximum of each column in a (make sure types of a and amin match) 1880 */ 1881 #define SUMA_MAX_MAT_COL(a, rows, cols, amax) { \ 1882 int m_IX, _JX; \ 1883 for (m_IX = 0; m_IX < cols ; m_IX++) { \ 1884 amax[m_IX]=a[0][m_IX]; \ 1885 for (_JX = 1 ; _JX < rows ; _JX++) \ 1886 if (a[_JX][m_IX] > amax[m_IX]) amax[m_IX] = a[_JX][m_IX];\ 1887 }\ 1888 } 1889 1890 /*! \def SUMA_MIN_MAX_SUM_MAT_COL(a, rows, cols, amin, amax, asum) 1891 \brief SUMA_MIN_MAX_SUM_MAT_COL macro for minimum, maximum and sum of each column in a matrix 1892 a pointer to matrix (**) 1893 rows number of rows 1894 cols number of cols 1895 amin minimum of each column in a (make sure types of a and amin match) 1896 amax maximum of each column in a (make sure types of a and amin match) 1897 asum sum of each column in a (the mean is not computed because the / operation would then depend on the type of a) 1898 1899 */ 1900 #define SUMA_MIN_MAX_SUM_MAT_COL(a, rows, cols, amin, amax, asum) { \ 1901 int m_IX, _JX; \ 1902 for (m_IX = 0; m_IX < cols ; m_IX++) { \ 1903 amax[m_IX]=a[0][m_IX]; \ 1904 amin[m_IX]=a[0][m_IX]; \ 1905 asum[m_IX]=a[0][m_IX]; \ 1906 for (_JX = 1 ; _JX < rows ; _JX++) { \ 1907 if (a[_JX][m_IX] > amax[m_IX]) amax[m_IX] = a[_JX][m_IX];\ 1908 if (a[_JX][m_IX] < amin[m_IX]) amin[m_IX] = a[_JX][m_IX];\ 1909 asum[m_IX] += a[_JX][m_IX]; \ 1910 } \ 1911 }\ 1912 } 1913 1914 /*! \def SUMA_MIN_MAX_SUM_VECMAT_COL(a, rows, cols, amin, amax, asum) 1915 \brief SUMA_MIN_MAX_SUM_VECMAT_COL macro for minimum, maximum and sum of each column in a matrix stored in vector format 1916 matrix 1 2 3 1917 4 5 6 1918 is stored as 1 2 3 4 5 6 ... 1919 1920 a pointer to vector containing rwos x cols elements 1921 rows number of rows 1922 cols number of cols 1923 amin minimum of each column in a (make sure types of a and amin match) 1924 amax maximum of each column in a (make sure types of a and amin match) 1925 asum sum of each column in a (the mean is not computed because the / operation would then depend on the type of a) 1926 1927 \sa SUMA_MIN_MAX_SUM_VECMAT_MASK_COL 1928 */ 1929 #define SUMA_MIN_MAX_SUM_VECMAT_COL(a, rows, cols, amin, amax, asum) { \ 1930 int m_IX, m_JX, m_id; \ 1931 for (m_IX = 0; m_IX < (cols) ; m_IX++) { \ 1932 (amax)[m_IX]=(a)[m_IX]; \ 1933 (amin)[m_IX]=(a)[m_IX]; \ 1934 (asum)[m_IX]=(a)[m_IX]; \ 1935 for (m_JX = 1 ; m_JX < (rows) ; m_JX++) { \ 1936 m_id = (cols) * m_JX + m_IX; \ 1937 if ((a)[m_id] > (amax)[m_IX]) (amax)[m_IX] = (a)[m_id];\ 1938 if ((a)[m_id] < (amin)[m_IX]) (amin)[m_IX] = (a)[m_id];\ 1939 (asum)[m_IX] += (a)[m_id]; \ 1940 } \ 1941 } \ 1942 } 1943 /*! \def SUMA_MIN_MAX_SUM_VECMAT_MASK_COL(a, rows, cols, rowmask, amin, amax, asum) 1944 \brief SUMA_MIN_MAX_SUM_VECMAT_MASK_COL macro for minimum, maximum and sum of each column in a matrix stored in vector format 1945 ONLY rows n where rowmask[n] is not 0 are used 1946 matrix 1 2 3 1947 4 5 6 1948 7 8 9 1949 is stored as 1 2 3 4 5 6 7 8 9... 1950 1951 a pointer to vector containing rwos x cols elements 1952 rows number of rows 1953 cols number of cols 1954 rowmask pointer to vector containing rows x 1 mask values 1955 amin minimum of each column in a (make sure types of a and amin match) 1956 amax maximum of each column in a (make sure types of a and amin match) 1957 asum sum of each column in a (the mean is not computed because the / operation would then depend on the type of a) 1958 1959 */ 1960 #define SUMA_MIN_MAX_SUM_VECMAT_MASK_COL(a, rows, cols, rowmask, amin, amax, asum) { \ 1961 int m_IX, m_JX, m_id, m_start_row=0; \ 1962 m_JX = 0; \ 1963 while (!m_start_row && m_JX < rows) { \ 1964 if (rowmask[m_JX]) m_start_row = m_JX+1; \ 1965 ++m_JX; \ 1966 } \ 1967 --m_start_row; \ 1968 for (m_IX = 0; m_IX < cols ; m_IX++) { \ 1969 amax[m_IX]=a[cols * m_start_row + m_IX]; \ 1970 amin[m_IX]=a[cols * m_start_row + m_IX]; \ 1971 asum[m_IX]=a[cols * m_start_row + m_IX]; \ 1972 ++m_start_row; \ 1973 for (m_JX = m_start_row ; m_JX < rows ; m_JX++) { \ 1974 if (rowmask[m_JX]) { \ 1975 m_id = cols * m_JX + m_IX; \ 1976 if (a[m_id] > amax[m_IX]) amax[m_IX] = a[m_id];\ 1977 if (a[m_id] < amin[m_IX]) amin[m_IX] = a[m_id];\ 1978 asum[m_IX] += a[m_id]; \ 1979 } \ 1980 } \ 1981 } \ 1982 } 1983 1984 /*! 1985 SUMA_MAT_TO_VEC MACRO rearranges a matrix into a vector 1986 one row after the other: 1987 1988 SUMA_MAT_TO_VEC(a,b,rows,cols,typea,typeb) 1989 1990 a pointer to input 2D matrix. 1991 b pointer to resultant D matrix. 1992 rows number of rows in matrix a 1993 cols number of columns in matrix a 1994 typea legal C type describing the type of a 1995 typeb legal C type describing the type of b 1996 1997 */ 1998 #define SUMA_MAT_TO_VEC(a,b,rows,cols,typea,typeb) { \ 1999 typea **_AMX = (typea **)a; \ 2000 typeb *_PTB = (typeb *)b; \ 2001 typea *_PTA; \ 2002 int m_IX,_JX; \ 2003 for(m_IX = 0 ; m_IX < rows ; m_IX++) { \ 2004 _PTA = _AMX[m_IX ]; \ 2005 for(_JX = 0 ; _JX < cols ; _JX++) \ 2006 *_PTB++ = (typeb) (*_PTA++); \ 2007 } \ 2008 } 2009 2010 2011 /* 2012 SUMA_COPY_VEC macro: copies the contents of vector a into vector b 2013 2014 SUMA_COPY_VEC(a,b,len,typea,typeb) 2015 2016 a pointer to input vector. 2017 b pointer to output vector. 2018 len length of vectors (integer). 2019 typea legal C type describing the type of a data. 2020 typeb legal C type describing the type of b data. 2021 2022 */ 2023 2024 #define SUMA_COPY_VEC(a,b,len,typea,typeb) { \ 2025 typea *_PTA = (typea *)a; \ 2026 typeb *_PTB = (typeb *)b; \ 2027 int m_IX; \ 2028 for(m_IX = 0 ; m_IX < (len) ; m_IX++) \ 2029 *(_PTB)++ = (typeb)(*(_PTA)++); \ 2030 } 2031 2032 #define SUMA_VEC_DIFF3(a, b) (( ((a)[0] != (b)[0]) || ((a)[1] != (b)[1]) \ 2033 || ((a)[2] != (b)[2]) ) ? 1:0 ) 2034 #define SUMA_VEC_DIFF4(a, b) (( ((a)[0] != (b)[0]) || ((a)[1] != (b)[1]) \ 2035 || ((a)[2] != (b)[2]) || ((a)[3] != (b)[3]))\ 2036 ? 1:0 ) 2037 2038 /*! 2039 SUMA_INIT_VEC macro: initializes values in a vector 2040 SUMA_INIT_VEC(a,len,val,typea) 2041 */ 2042 #define SUMA_INIT_VEC(a,len,val,typea) { \ 2043 int m_i; \ 2044 for (m_i = 0; m_i < (len) ; m_i ++) \ 2045 a[m_i] = (typea)val; \ 2046 } 2047 2048 #define SUMA_COPY_VALUE_IN_VEC(a, b, ia, ib, typea, typeb) { \ 2049 typea *_PTA = (typea *)a; \ 2050 typeb *_PTB = (typeb *)b; \ 2051 _PTB[ib] = (typeb)_PTA[ia]; \ 2052 } 2053 2054 #define SUMA_ASSIGN_VALUE_IN_VEC(a, ia, typea, val) { \ 2055 typea *_PTA = (typea *)a; \ 2056 _PTA[ia] = (typea)val; \ 2057 } 2058 2059 2060 /*! 2061 SUMA_DOTP_VEC macro: 2062 2063 FORMS THE SUM OF PRODUCTS OF TWO VECTORS (a,b) AND 2064 PUTS THE RESULT IN THE PREVIOUSLY DEFINED VARIABLE s. 2065 2066 SUMA_DOTP_VEC(a,b,s,len,typea,typeb) 2067 2068 a pointer to first vector. 2069 b pointer to second vector. 2070 s variable used to store result (not a pointer). 2071 len length of vectors (integer). 2072 typea legal C type describing the type of a data. 2073 typeb legal C type describing the type of b data. 2074 2075 WARNING: The input data vectors are not cast to the type of s. 2076 This means that at least one of the input types must 2077 be able to represent the individual products without 2078 overflow. 2079 2080 */ 2081 2082 #define SUMA_DOTP_VEC(a,b,s,len,typea,typeb) { \ 2083 typea *_PTA = a; \ 2084 typeb *_PTB = b; \ 2085 int m_IX; \ 2086 s = (*_PTA++) * (*_PTB++); \ 2087 for(m_IX = 1 ; m_IX < (len) ; m_IX++) \ 2088 s += (*_PTA++) * (*_PTB++); \ 2089 } 2090 2091 /*! 2092 \brief Macro to calculate the unit direction vector U from P1-->P2 and 2093 the distance Un between P1 and P2. 2094 If Un is 0, U is all zeros 2095 \param P1/P2 (float *) 3-elements arrays containing XYZ of P1 and P2 2096 \param U (float *) 3-elements array to contain unit direction vector 2097 \param Un (float) the norm of |P1--P2| 2098 */ 2099 #define SUMA_UNIT_VEC(P1, P2, U, Un){ \ 2100 /* Calculate normalized unit vector of line formed by P1, P2 */ \ 2101 U[0] = P2[0] - P1[0]; \ 2102 U[1] = P2[1] - P1[1]; \ 2103 U[2] = P2[2] - P1[2]; \ 2104 Un = sqrt(U[0]*U[0] + U[1]*U[1] + U[2]*U[2]); \ 2105 if (Un) { \ 2106 U[0] /= Un; U[1] /= Un; U[2] /= Un; \ 2107 }else { \ 2108 U[0] = U[1] = U[2] = 0; \ 2109 } \ 2110 } \ 2111 2112 /*! 2113 Pick a color based on the direction of the unit vector 2114 between pa and pb. Direction sign ignored 2115 */ 2116 #define SUMA_SEG_DELTA_COL(pa,pb,U, Un) {\ 2117 SUMA_UNIT_VEC(pa, pb, U, Un); \ 2118 U[0] = SUMA_ABS(U[0]); \ 2119 U[1] = SUMA_ABS(U[1]); \ 2120 U[2] = SUMA_ABS(U[2]); U[3]=1.0; \ 2121 } 2122 2123 /*! 2124 \brief Macro to calculate normal of a triangle 2125 \params P1, P2, P3: XYZ coordinates forming triangle 2126 \param norm: Contains UN-NORMALIZED normal upon return 2127 use SUMA_TRI_NORM_NORM for normalized normals 2128 */ 2129 #define SUMA_TRI_NORM(P1, P2, P3, norm){ \ 2130 double m_d1[3], m_d2[3]; \ 2131 m_d1[0] = P1[0] - P2[0]; \ 2132 m_d2[0] = P2[0] - P3[0]; \ 2133 m_d1[1] = P1[1] - P2[1]; \ 2134 m_d2[1] = P2[1] - P3[1]; \ 2135 m_d1[2] = P1[2] - P2[2]; \ 2136 m_d2[2] = P2[2] - P3[2]; \ 2137 norm[0] = m_d1[1]*m_d2[2] - m_d1[2]*m_d2[1]; \ 2138 norm[1] = m_d1[2]*m_d2[0] - m_d1[0]*m_d2[2]; \ 2139 norm[2] = m_d1[0]*m_d2[1] - m_d1[1]*m_d2[0]; \ 2140 } 2141 #define SUMA_TRI_NORM_NORM(P1, P2, P3, norm){ \ 2142 double m_d; \ 2143 SUMA_TRI_NORM(P1, P2, P3, norm); \ 2144 m_d = sqrt(norm[0]*norm[0]+norm[1]*norm[1]+norm[2]*norm[2]); \ 2145 if (m_d == 0.0f) { norm[0] = norm[1] = norm[2] = 0.0; } \ 2146 else { norm[0] /= m_d; norm[1] /= m_d; norm[2] /= m_d;} \ 2147 } 2148 2149 2150 /*! 2151 \brief Macro to calculate the distance Un from P1-->P2 a 2152 \param P1/P2 (float *) 3-elements arrays containing XYZ of P1 and P2 2153 \param Un (float) the norm of |P1--P2| 2154 */ 2155 #define SUMA_SEG_NORM(P1, P2, Un){ \ 2156 static double m_dx, m_dy,m_dz; \ 2157 /* Calculate normalized unit vector of line formed by P1, P2 */ \ 2158 m_dx = P2[0] - P1[0]; \ 2159 m_dy = P2[1] - P1[1]; \ 2160 m_dz = P2[2] - P1[2]; \ 2161 Un = sqrt(m_dx*m_dx + m_dy*m_dy + m_dz*m_dz); \ 2162 } \ 2163 2164 /*! 2165 SUMA_MULT_MAT MACRO FOR MATRIX MULTIPLICATION: 2166 2167 SUMA_MULT_MAT(a,b,c,rowsa,colsa,colsb,typea,typeb,typec) 2168 2169 a pointer to first matirx. 2170 b pointer to second matrix. 2171 c pointer to result matrix. 2172 rowsa number of rows in matrix a 2173 colsa number of columns in matrix a 2174 colsb number of columns in matrix b 2175 typea legal C type describing the type of a 2176 typeb legal C type describing the type of b 2177 typec legal C type describing the type of c 2178 2179 */ 2180 2181 #define SUMA_MULT_MAT(a,b,c,rowsa,colsa,colsb,typea,typeb,typec) { \ 2182 typea **_AMX = (typea **)a; \ 2183 typeb **_BMX = (typeb **)b; \ 2184 typec **_CMX = (typec **)c; \ 2185 typea *_PTA; \ 2186 typeb *_PTB; \ 2187 typec *_PTC; \ 2188 int m_IX,_JX,_KX; \ 2189 for(m_IX = 0 ; m_IX < rowsa ; m_IX++) { \ 2190 _PTC = _CMX[m_IX]; \ 2191 _PTB = _BMX[0]; \ 2192 for(_JX = 0 ; _JX < colsb ; _JX++) { \ 2193 _PTA = _AMX[m_IX]; \ 2194 *_PTC = (*_PTA++) * (*_PTB++); \ 2195 for(_KX = 1 ; _KX < colsa ; _KX++) \ 2196 *_PTC += (*_PTA++)* _BMX[_KX][_JX]; \ 2197 _PTC++; \ 2198 } \ 2199 } \ 2200 } 2201 2202 2203 /*! 2204 SUMA_ADD_MAT MACRO FOR MATRIX ADDITION: 2205 2206 SUMA_ADD_MAT(a,b,c,rowsa,colsa,typea,typeb,typec) 2207 2208 a pointer to first matirx. 2209 b pointer to second matrix. 2210 c pointer to result matrix. 2211 rowsa number of rows in matrix a 2212 colsa number of columns in matrix a 2213 typea legal C type describing the type of a 2214 typeb legal C type describing the type of b 2215 typec legal C type describing the type of c 2216 2217 */ 2218 2219 #define SUMA_ADD_MAT(a,b,c,rowsa,colsa,typea,typeb,typec) { \ 2220 typea **_AMX = (typea **)a; \ 2221 typeb **_BMX = (typeb **)b; \ 2222 typec **_CMX = (typec **)c; \ 2223 int m_IX,_JX; \ 2224 for(m_IX = 0 ; m_IX < rowsa ; m_IX++) { \ 2225 for(_JX = 0 ; _JX < colsa ; _JX++) { \ 2226 _CMX[m_IX][_JX] = _AMX[m_IX][_JX] + _BMX[m_IX][_JX]; \ 2227 } \ 2228 } \ 2229 } 2230 2231 /*! 2232 SUMA_SUB_MAT MACRO FOR MATRIX SUBTRACTION: 2233 2234 SUMA_SUB_MAT(a,b,c,rowsa,colsa,typea,typeb,typec) 2235 2236 a pointer to first matirx. 2237 b pointer to second matrix. 2238 c pointer to result matrix. 2239 rowsa number of rows in matrix a 2240 colsa number of columns in matrix a 2241 typea legal C type describing the type of a 2242 typeb legal C type describing the type of b 2243 typec legal C type describing the type of c 2244 2245 */ 2246 2247 #define SUMA_SUB_MAT(a,b,c,rowsa,colsa,typea,typeb,typec) { \ 2248 typea **_AMX = (typea **)a; \ 2249 typeb **_BMX = (typeb **)b; \ 2250 typec **_CMX = (typec **)c; \ 2251 int m_IX,_JX; \ 2252 for(m_IX = 0 ; m_IX < rowsa ; m_IX++) { \ 2253 for(_JX = 0 ; _JX < colsa ; _JX++) { \ 2254 _CMX[m_IX][_JX] = _AMX[m_IX][_JX] - _BMX[m_IX][_JX]; \ 2255 } \ 2256 } \ 2257 } 2258 2259 2260 /*! 2261 SUMA_TRANSP_MAT MACRO FOR MATRIX TRANSPOSE: 2262 2263 SUMA_TRANSP_MAT(a,b,rowsa,colsa,typea,typeb) 2264 2265 a pointer to first matirx. 2266 b pointer to result matrix. 2267 rowsa number of rows in matrix a 2268 colsa number of columns in matrix a 2269 typea legal C type describing the type of a 2270 typeb legal C type describing the type of b 2271 2272 */ 2273 2274 #define SUMA_TRANSP_MAT(a,b,rowsa,colsa,typea,typeb) { \ 2275 typea **_AMX = (typea **)a; \ 2276 typeb **_BMX = (typeb **)b; \ 2277 int m_IX,_JX; \ 2278 for(m_IX = 0 ; m_IX < rowsa ; m_IX++) { \ 2279 for(_JX = 0 ; _JX < colsa ; _JX++) { \ 2280 _BMX[_JX][m_IX] = _AMX[m_IX][_JX]; \ 2281 } \ 2282 } \ 2283 } 2284 2285 /*! 2286 SUMA_RGBmat_2_GLCOLAR4 copies an N x 3 RGB matrix into a 4N x 1 GL color array format 2287 SUMA_RGBmat_2_GLCOLAR4(RGBmat, glcolar, N) 2288 RGBmat (float **) N x 3 matrix of RGB values 2289 glcolar (GLfloat *) (4 N) x 1 vector 2290 *** Mar 17 04: 2291 Added the SUMA_RGBvec versions using RGBvec instead of RGBmat. 2292 */ 2293 2294 #define SUMA_RGBmat_2_GLCOLAR4(RGBmat, glcolar, nrgb) {\ 2295 int m_I, m_I4 = 0; \ 2296 for (m_I=0; m_I < nrgb; ++m_I) {\ 2297 glcolar[m_I4] = RGBmat[m_I][0]; ++m_I4;\ 2298 glcolar[m_I4] = RGBmat[m_I][1]; ++m_I4;\ 2299 glcolar[m_I4] = RGBmat[m_I][2]; ++m_I4;\ 2300 ++m_I4;\ 2301 }\ 2302 } 2303 #define SUMA_RGBvec_2_GLCOLAR4(RGBvec, glcolar, nrgb) {\ 2304 int m_I, m_I4 = 0, m_I3=0; \ 2305 for (m_I=0; m_I < nrgb; ++m_I) {\ 2306 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2307 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2308 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2309 ++m_I4;\ 2310 }\ 2311 } 2312 2313 /*! 2314 SUMA_GLCOLAR4_2_RGBmat copies 4N x 1 GL color array into an N x 3 RGB matrix format 2315 2316 SUMA_GLCOLAR4_2_RGBmat (glcolar, RGBmat, N) 2317 glcolar (GLfloat *) (4 N) x 1 vector 2318 RGBmat (float **) N x 3 matrix of RGB values 2319 *** Mar 17 04: 2320 Added the SUMA_RGBvec versions using RGBvec instead of RGBmat. 2321 */ 2322 2323 #define SUMA_GLCOLAR4_2_RGBmat(glcolar, RGBmat, nrgb) {\ 2324 int m_I, m_I4 = 0; \ 2325 for (m_I=0; m_I < nrgb; ++m_I) {\ 2326 RGBmat[m_I][0] = glcolar[m_I4]; ++m_I4;\ 2327 RGBmat[m_I][1] = glcolar[m_I4]; ++m_I4;\ 2328 RGBmat[m_I][2] = glcolar[m_I4]; ++m_I4;\ 2329 ++m_I4;\ 2330 }\ 2331 } 2332 #define SUMA_GLCOLAR4_2_RGBvec(glcolar, RGBvec, nrgb) {\ 2333 int m_I, m_I4 = 0, m_I3=0; \ 2334 for (m_I=0; m_I < nrgb; ++m_I) {\ 2335 RGBvec[m_I3] = glcolar[m_I4]; ++m_I4; ++m_I3;\ 2336 RGBvec[m_I3] = glcolar[m_I4]; ++m_I4; ++m_I3;\ 2337 RGBvec[m_I3] = glcolar[m_I4]; ++m_I4; ++m_I3;\ 2338 ++m_I4;\ 2339 }\ 2340 } 2341 2342 2343 /*! 2344 SUMA_FillBlanks_GLCOLAR4 2345 fills nodes that received no color with the Nocolor Color 2346 SUMA_FillBlanks_GLCOLAR4(isColored, N_Nodes, R, G, B, glcolar) 2347 isColored (SUMA_Boolean) N_Nodes x 1 vector indicating that a node was colored. 2348 N_Nodes (int) total number of nodes 2349 R, G, B (float 0..1) RGB values of NoColor color 2350 glcolar (GLfloat *) (4 N_Nodes) x 1 vector 2351 */ 2352 2353 #define SUMA_FillBlanks_GLCOLAR4(isColored, N_Nodes, R, G, B, glcolar) {\ 2354 int m_I, m_I4; \ 2355 for (m_I=0; m_I < N_Nodes; ++m_I) {\ 2356 if (!isColored[m_I]) {\ 2357 m_I4 = 4*m_I; \ 2358 glcolar[m_I4] = R; ++m_I4;\ 2359 glcolar[m_I4] = G; ++m_I4;\ 2360 glcolar[m_I4] = B; ++m_I4;\ 2361 }\ 2362 }\ 2363 } 2364 2365 /*! 2366 This macro used to be called: SUMA_RGBmat_FullNoGlobNoLoc2_GLCOLAR4_opacity 2367 SUMA_RGB_FnGnL_AR4op_opacity 2368 copies an N x 3 RGB matrix into a 4N x 1 GL color array format 2369 F (Full) means that N is equal to all the nodes in the surface 2370 nG (NoGlob) means no opacity is applied to the color values (fully opaque) 2371 nL (NoLocal) means no local gain (per node) is applied to the color values 2372 2373 SUMA_RGB_FnGnL_AR4op(RGBmat, glcolar, N, add) 2374 RGBmat (float **) N x 3 matrix of RGB values 2375 glcolar (GLfloat *) (4 N) x 1 vector 2376 isColored (SUMA_Boolean) N_Nodes x 1 vector indicating that a node was colored. ONLY YUP/1 are placed 2377 *** Mar 17 04: 2378 Added the SUMA_RGBv versions using RGBvec instead of RGBmat. 2379 2380 */ 2381 2382 #define SUMA_RGB_FnGnL_AR4op(RGBmat, glcolar, nrgb, isColored) {\ 2383 int m_I, m_I4 = 0; \ 2384 for (m_I=0; m_I < nrgb; ++m_I) {\ 2385 isColored[m_I] = YUP;\ 2386 glcolar[m_I4] = RGBmat[m_I][0]; ++m_I4;\ 2387 glcolar[m_I4] = RGBmat[m_I][1]; ++m_I4;\ 2388 glcolar[m_I4] = RGBmat[m_I][2]; ++m_I4;\ 2389 ++m_I4;\ 2390 }\ 2391 } 2392 2393 #define SUMA_RGBv_FnGnL_AR4op(RGBvec, glcolar, nrgb, isColored) {\ 2394 int m_I, m_I4 = 0, m_I3=0; \ 2395 for (m_I=0; m_I < nrgb; ++m_I) {\ 2396 isColored[m_I] = YUP;\ 2397 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2398 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2399 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2400 ++m_I4;\ 2401 }\ 2402 } 2403 /*! 2404 This macro used to be called: SUMA_RGBmat_FullGlobNoLoc2_GLCOLAR4_opacity 2405 SUMA_RGB_FGnL_AR4op 2406 copies an N x 3 RGB matrix into a 4N x 1 GL color array format 2407 F (Full) means that N is equal to all the nodes in the surface 2408 G (Glob) means an opacity is applied to the color values 2409 nL (NoLocal) means no local gain (per node) is applied to the color values 2410 2411 SUMA_RGB_FGnL_AR4op(RGBmat, glcolar, N, opacity, add) 2412 RGBmat (float **) N x 3 matrix of RGB values 2413 glcolar (GLfloat *) (4 N) x 1 vector 2414 opacity (float) an opacity factor applied to each color R, G, B values in the entire list before adding it 2415 to a pre-exising color opacity is not applied to the color of nodes that had not been colored thus far. 2416 isColored (SUMA_Boolean) N_Nodes x 1 vector indicating that a node was colored. ONLY YUP/1 are placed 2417 2418 *** Dec 04 03: 2419 Variant SUMA_RGB_FGnL_AR4op2 2420 removed opacity scaling of RGBmat[m_I] and added clipping to 1.0 2421 *** Mar 17 04: 2422 Added the SUMA_RGBv versions using RGBvec instead of RGBmat. 2423 */ 2424 2425 #define SUMA_RGB_FGnL_AR4op(RGBmat, glcolar, nrgb, opacity, isColored) {\ 2426 int m_I, m_I4 = 0; \ 2427 float m_of;\ 2428 m_of = 1-opacity;\ 2429 for (m_I=0; m_I < nrgb; ++m_I) {\ 2430 if (isColored[m_I]) { \ 2431 glcolar[m_I4] = m_of * glcolar[m_I4] + opacity * RGBmat[m_I][0]; ++m_I4;\ 2432 glcolar[m_I4] = m_of * glcolar[m_I4] + opacity * RGBmat[m_I][1]; ++m_I4;\ 2433 glcolar[m_I4] = m_of * glcolar[m_I4] + opacity * RGBmat[m_I][2]; ++m_I4;\ 2434 } else { /* not yet colored */\ 2435 glcolar[m_I4] = RGBmat[m_I][0]; ++m_I4;\ 2436 glcolar[m_I4] = RGBmat[m_I][1]; ++m_I4;\ 2437 glcolar[m_I4] = RGBmat[m_I][2]; ++m_I4;\ 2438 isColored[m_I] = YUP;\ 2439 } \ 2440 ++m_I4;\ 2441 }\ 2442 } 2443 #define SUMA_RGBv_FGnL_AR4op(RGBvec, glcolar, nrgb, opacity, isColored) {\ 2444 int m_I, m_I4 = 0, m_I3 = 0; \ 2445 float m_of;\ 2446 m_of = 1-opacity;\ 2447 for (m_I=0; m_I < nrgb; ++m_I) {\ 2448 if (isColored[m_I]) { \ 2449 glcolar[m_I4] = m_of * glcolar[m_I4] + opacity * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2450 glcolar[m_I4] = m_of * glcolar[m_I4] + opacity * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2451 glcolar[m_I4] = m_of * glcolar[m_I4] + opacity * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2452 } else { /* not yet colored */\ 2453 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2454 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2455 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2456 isColored[m_I] = YUP;\ 2457 } \ 2458 ++m_I4;\ 2459 }\ 2460 } 2461 #define SUMA_RGB_FGnL_AR4op2(RGBmat, glcolar, nrgb, opacity, isColored) {\ 2462 int m_I, m_I4 = 0; \ 2463 float m_of;\ 2464 m_of = 1-opacity;\ 2465 for (m_I=0; m_I < nrgb; ++m_I) {\ 2466 if (isColored[m_I]) { \ 2467 glcolar[m_I4] = m_of * glcolar[m_I4] + RGBmat[m_I][0]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4;\ 2468 glcolar[m_I4] = m_of * glcolar[m_I4] + RGBmat[m_I][1]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4;\ 2469 glcolar[m_I4] = m_of * glcolar[m_I4] + RGBmat[m_I][2]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4;\ 2470 } else { /* not yet colored */\ 2471 glcolar[m_I4] = RGBmat[m_I][0]; ++m_I4;\ 2472 glcolar[m_I4] = RGBmat[m_I][1]; ++m_I4;\ 2473 glcolar[m_I4] = RGBmat[m_I][2]; ++m_I4;\ 2474 isColored[m_I] = YUP;\ 2475 } \ 2476 ++m_I4;\ 2477 }\ 2478 } 2479 #define SUMA_RGBv_FGnL_AR4op2(RGBvec, glcolar, nrgb, opacity, isColored) {\ 2480 int m_I, m_I4 = 0, m_I3 = 0; \ 2481 float m_of;\ 2482 m_of = 1-opacity;\ 2483 for (m_I=0; m_I < nrgb; ++m_I) {\ 2484 if (isColored[m_I]) { \ 2485 glcolar[m_I4] = m_of * glcolar[m_I4] + RGBvec[m_I3]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4; ++m_I3;\ 2486 glcolar[m_I4] = m_of * glcolar[m_I4] + RGBvec[m_I3]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4; ++m_I3;\ 2487 glcolar[m_I4] = m_of * glcolar[m_I4] + RGBvec[m_I3]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4; ++m_I3;\ 2488 } else { /* not yet colored */\ 2489 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2490 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2491 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2492 isColored[m_I] = YUP;\ 2493 } \ 2494 ++m_I4;\ 2495 }\ 2496 } 2497 /*! 2498 This macro used to be called: SUMA_RGBmat_FullGlobLoc2_GLCOLAR4_opacity 2499 but name was too long for some compilers 2500 2501 SUMA_RGB_FGL_AR4op 2502 copies an N x 3 RGB matrix into a 4N x 1 GL color array format 2503 F (Full) means that N is equal to all the nodes in the surface 2504 G (Glob) means an opacity is applied to the color values 2505 L (Local) means a local gain (per node) is applied to the color values 2506 2507 SUMA_RGB_FGL_AR4op(RGBmat, glcolar, N, opacity, locgain, add) 2508 RGBmat (float **) N x 3 matrix of RGB values 2509 glcolar (GLfloat *) (4 N) x 1 vector 2510 opacity (float) an opacity factor applied to each color R, G, B values in the entire list before adding it 2511 to a pre-exising color opacity is not applied to the color of nodes that had not been colored thus far. 2512 locgain (float *) a N x 1 vector of gains applied to their respective nodes 2513 isColored (SUMA_Boolean) N_Nodes x 1 vector indicating that a node was colored. ONLY YUP/1 are placed 2514 2515 *** Dec 04 03: 2516 Variant SUMA_RGB_FGL_AR4op2 2517 removed opacity from m_of and added clipping to 1.0 2518 *** Mar 17 04: 2519 Added the SUMA_RGBv versions using RGBvec instead of RGBmat. 2520 */ 2521 2522 #define SUMA_RGB_FGL_AR4op(RGBmat, glcolar, nrgb, opacity, locgain, isColored) {\ 2523 int m_I, m_I4 = 0; \ 2524 float m_of, m_of2;\ 2525 for (m_I=0; m_I < nrgb; ++m_I) {\ 2526 if (isColored[m_I]) { \ 2527 m_of = locgain[m_I] * opacity; /* Dec 04 03 changed from locgain[m_I] * opacity; */\ 2528 m_of2 = (1-opacity);\ 2529 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBmat[m_I][0]; ++m_I4;\ 2530 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBmat[m_I][1]; ++m_I4;\ 2531 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBmat[m_I][2]; ++m_I4;\ 2532 } else { /* not yet colored */\ 2533 glcolar[m_I4] = locgain[m_I] * RGBmat[m_I][0]; ++m_I4;\ 2534 glcolar[m_I4] = locgain[m_I] * RGBmat[m_I][1]; ++m_I4;\ 2535 glcolar[m_I4] = locgain[m_I] * RGBmat[m_I][2]; ++m_I4;\ 2536 isColored[m_I] = YUP;\ 2537 } \ 2538 ++m_I4;\ 2539 }\ 2540 } 2541 #define SUMA_RGBv_FGL_AR4op(RGBvec, glcolar, nrgb, opacity, locgain, isColored) {\ 2542 int m_I, m_I4 = 0, m_I3 = 0; \ 2543 float m_of, m_of2;\ 2544 for (m_I=0; m_I < nrgb; ++m_I) {\ 2545 if (isColored[m_I]) { \ 2546 m_of = locgain[m_I] * opacity; /* Dec 04 03 changed from locgain[m_I] * opacity; */\ 2547 m_of2 = (1-opacity);\ 2548 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2549 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2550 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2551 } else { /* not yet colored */\ 2552 glcolar[m_I4] = locgain[m_I] * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2553 glcolar[m_I4] = locgain[m_I] * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2554 glcolar[m_I4] = locgain[m_I] * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2555 isColored[m_I] = YUP;\ 2556 } \ 2557 ++m_I4;\ 2558 }\ 2559 } 2560 #define SUMA_RGB_FGL_AR4op2(RGBmat, glcolar, nrgb, opacity, locgain, isColored) {\ 2561 int m_I, m_I4 = 0; \ 2562 float m_of, m_of2;\ 2563 for (m_I=0; m_I < nrgb; ++m_I) {\ 2564 if (isColored[m_I]) { \ 2565 m_of = locgain[m_I]; /* Dec 04 03 changed from locgain[m_I] * opacity; */\ 2566 m_of2 = (1-opacity);\ 2567 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBmat[m_I][0]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4;\ 2568 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBmat[m_I][1]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4;\ 2569 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBmat[m_I][2]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4;\ 2570 } else { /* not yet colored */\ 2571 glcolar[m_I4] = locgain[m_I] * RGBmat[m_I][0]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4;\ 2572 glcolar[m_I4] = locgain[m_I] * RGBmat[m_I][1]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4;\ 2573 glcolar[m_I4] = locgain[m_I] * RGBmat[m_I][2]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4;\ 2574 isColored[m_I] = YUP;\ 2575 } \ 2576 ++m_I4;\ 2577 }\ 2578 } 2579 #define SUMA_RGBv_FGL_AR4op2(RGBvec, glcolar, nrgb, opacity, locgain, isColored) {\ 2580 int m_I, m_I4 = 0, m_I3=0; \ 2581 float m_of, m_of2;\ 2582 for (m_I=0; m_I < nrgb; ++m_I) {\ 2583 if (isColored[m_I]) { \ 2584 m_of = locgain[m_I]; /* Dec 04 03 changed from locgain[m_I] * opacity; */\ 2585 m_of2 = (1-opacity);\ 2586 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBvec[m_I3]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4; ++m_I3;\ 2587 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBvec[m_I3]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4; ++m_I3;\ 2588 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBvec[m_I3]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4; ++m_I3;\ 2589 } else { /* not yet colored */\ 2590 glcolar[m_I4] = locgain[m_I] * RGBvec[m_I3]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4; ++m_I3;\ 2591 glcolar[m_I4] = locgain[m_I] * RGBvec[m_I3]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4; ++m_I3;\ 2592 glcolar[m_I4] = locgain[m_I] * RGBvec[m_I3]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4; ++m_I3;\ 2593 isColored[m_I] = YUP;\ 2594 } \ 2595 ++m_I4;\ 2596 }\ 2597 } 2598 /*! 2599 This macro used to be called: SUMA_RGBmat_FullNoGlobLoc2_GLCOLAR4_opacity 2600 2601 SUMA_RGB_FnGL_AR4op 2602 copies an N x 3 RGB matrix into a 4N x 1 GL color array format 2603 F (Full) means that N is equal to all the nodes in the surface 2604 nG (NoGlob) means no opacity is applied to the color values (fully opaque) 2605 L (Local) means a local gain (per node) is applied to the color values 2606 2607 SUMA_RGB_FnGL_AR4op(RGBmat, glcolar, N, locgain, add) 2608 RGBmat (float **) N x 3 matrix of RGB values 2609 glcolar (GLfloat *) (4 N) x 1 vector 2610 locgain (float *) a N x 1 vector of gains applied to their respective nodes 2611 isColored (SUMA_Boolean) N_Nodes x 1 vector indicating that a node was colored. ONLY YUP/1 are placed 2612 *** Mar 17 04: 2613 Added the SUMA_RGBv versions using RGBvec instead of RGBmat. 2614 */ 2615 2616 #define SUMA_RGB_FnGL_AR4op(RGBmat, glcolar, nrgb, locgain, isColored) {\ 2617 int m_I, m_I4 = 0; \ 2618 float m_of;\ 2619 for (m_I=0; m_I < nrgb; ++m_I) {\ 2620 glcolar[m_I4] = locgain[m_I] * RGBmat[m_I][0]; ++m_I4;\ 2621 glcolar[m_I4] = locgain[m_I] * RGBmat[m_I][1]; ++m_I4;\ 2622 glcolar[m_I4] = locgain[m_I] * RGBmat[m_I][2]; ++m_I4;\ 2623 isColored[m_I] = YUP;\ 2624 ++m_I4;\ 2625 }\ 2626 } 2627 #define SUMA_RGBv_FnGL_AR4op(RGBvec, glcolar, nrgb, locgain, isColored) {\ 2628 int m_I, m_I4 = 0, m_I3=0; \ 2629 float m_of;\ 2630 for (m_I=0; m_I < nrgb; ++m_I) {\ 2631 glcolar[m_I4] = locgain[m_I] * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2632 glcolar[m_I4] = locgain[m_I] * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2633 glcolar[m_I4] = locgain[m_I] * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2634 isColored[m_I] = YUP;\ 2635 ++m_I4;\ 2636 }\ 2637 } 2638 /*! 2639 This macro used to be called: SUMA_RGBmat_PartNoGlobNoLoc2_GLCOLAR4_opacity 2640 SUMA_RGB_PnGnL_AR4op 2641 copies an N x 3 RGB matrix into a 4N x 1 GL color array format 2642 P (Part) means that colors are specified for some of the nodes only N < N_Nodes 2643 nG (NoGlob) means no opacity is applied to the color values (fully opaque) 2644 nL (NoLocal) means no local gain (per node) is applied to the color values 2645 2646 SUMA_RGB_PnGnL_AR4op(RGBmat, NodeId, glcolar, N, isColored, N_Nodes) 2647 RGBmat (float **) N x 3 matrix of RGB values 2648 NodeId (int *) N x 1 vector containing indices of nodes for wich color is specified in RGBmat 2649 glcolar (GLfloat *) (4 N_Nodes) x 1 vector 2650 isColored (SUMA_Boolean) N_Nodes x 1 vector indicating that a node was colored. ONLY YUP/1 are placed 2651 when a node is assigned a color. Values of isColored for nodes that have not been visited remain unchanged 2652 *** Mar 17 04: 2653 Added the SUMA_RGBv versions using RGBvec instead of RGBmat. 2654 */ 2655 2656 #define SUMA_RGB_PnGnL_AR4op(RGBmat, NodeId, glcolar, nrgb, isColored, N_Node) {\ 2657 int m_I, m_I4 = 0; \ 2658 for (m_I=0; m_I < nrgb; ++m_I) {\ 2659 if (NodeId[m_I] < N_Node) {\ 2660 m_I4 = 4*NodeId[m_I]; \ 2661 glcolar[m_I4] = RGBmat[m_I][0]; ++m_I4;\ 2662 glcolar[m_I4] = RGBmat[m_I][1]; ++m_I4;\ 2663 glcolar[m_I4] = RGBmat[m_I][2]; ++m_I4;\ 2664 isColored[NodeId[m_I]] = YUP;\ 2665 }\ 2666 }\ 2667 } 2668 2669 #define SUMA_RGBv_PnGnL_AR4op(RGBvec, NodeId, glcolar, nrgb, isColored, N_Node) {\ 2670 int m_I, m_I4 = 0, m_I3=0; \ 2671 for (m_I=0; m_I < nrgb; ++m_I) {\ 2672 if (NodeId[m_I] < N_Node) {\ 2673 m_I4 = 4*NodeId[m_I]; \ 2674 m_I3 = 3*m_I; \ 2675 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2676 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2677 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2678 isColored[NodeId[m_I]] = YUP;\ 2679 }\ 2680 }\ 2681 } 2682 /*! 2683 This macro used to be called: SUMA_RGBmat_PartGlobNoLoc2_GLCOLAR4_opacity 2684 SUMA_RGB_PGnL_AR4op 2685 copies an N x 3 RGB matrix into a 4N x 1 GL color array format 2686 P (Part) means that colors are specified for some of the nodes only N < N_Nodes 2687 G (Glob) means an opacity is applied to the color values 2688 nL (NoLocal) means no local gain (per node) is applied to the color values 2689 2690 SUMA_RGB_PGnL_AR4op(RGBmat, NodeId, glcolar, N, isColored, opacity, N_Nodes) 2691 RGBmat (float **) N x 3 matrix of RGB values 2692 NodeId (int *) N x 1 vector containing indices of nodes for wich color is specified in RGBmat 2693 glcolar (GLfloat *) (4 N_Nodes) x 1 vector 2694 isColored (SUMA_Boolean) N_Nodes x 1 vector indicating that a node was colored. ONLY YUP/1 are placed 2695 when a node is assigned a color. Values of isColored for nodes that have not been visited remain unchanged 2696 opacity (float) an opacity factor applied to each color R, G, B values in the entire list before adding it 2697 to a pre-exising color opacity is not applied to the color of nodes that had not been colored thus far. 2698 2699 *** Dec. 04 03: 2700 Variant SUMA_RGB_PGnL_AR4op2 2701 Instead of mixing like this: Col_new = (1 - opacity)*Col_1 + opacity *Col_2; 2702 I am using Col_new = (1 - opacity)*Col_1 + Col_2; if (Col_new > 1) Col_new = 1 2703 *** Mar 17 04: 2704 Added the SUMA_RGBv versions using RGBvec instead of RGBmat. 2705 2706 */ 2707 2708 #define SUMA_RGB_PGnL_AR4op(RGBmat, NodeId, glcolar, nrgb, isColored, opacity, N_Node) {\ 2709 int m_I, m_I4 = 0, m_II; \ 2710 float m_of;\ 2711 m_of = (1 - opacity); \ 2712 for (m_I=0; m_I < nrgb; ++m_I) {\ 2713 if (NodeId[m_I] < N_Node) {\ 2714 if (!isColored[NodeId[m_I]]) { /* a new color, put it down as it is, Opacity gain should not be applied Wed Apr 2 17:31:33 EST 2003*/\ 2715 m_I4 = 4*NodeId[m_I]; \ 2716 glcolar[m_I4] = RGBmat[m_I][0]; ++m_I4;\ 2717 glcolar[m_I4] = RGBmat[m_I][1]; ++m_I4;\ 2718 glcolar[m_I4] = RGBmat[m_I][2]; ++m_I4;\ 2719 isColored[NodeId[m_I]] = YUP;\ 2720 }else {/* mixing to be done */\ 2721 m_I4 = 4*NodeId[m_I]; \ 2722 glcolar[m_I4] = m_of * glcolar[m_I4] + opacity * RGBmat[m_I][0]; ++m_I4;\ 2723 glcolar[m_I4] = m_of * glcolar[m_I4] + opacity * RGBmat[m_I][1]; ++m_I4;\ 2724 glcolar[m_I4] = m_of * glcolar[m_I4] + opacity * RGBmat[m_I][2]; ++m_I4;\ 2725 }\ 2726 }\ 2727 }\ 2728 } 2729 #define SUMA_RGBv_PGnL_AR4op(RGBvec, NodeId, glcolar, nrgb, isColored, opacity, N_Node) {\ 2730 int m_I, m_I4 = 0, m_II, m_I3=0; \ 2731 float m_of;\ 2732 m_of = (1 - opacity); \ 2733 for (m_I=0; m_I < nrgb; ++m_I) {\ 2734 if (NodeId[m_I] < N_Node) {\ 2735 if (!isColored[NodeId[m_I]]) { /* a new color, put it down as it is, Opacity gain should not be applied Wed Apr 2 17:31:33 EST 2003*/\ 2736 m_I4 = 4*NodeId[m_I]; \ 2737 m_I3 = 3*m_I; \ 2738 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2739 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2740 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2741 isColored[NodeId[m_I]] = YUP;\ 2742 }else {/* mixing to be done */\ 2743 m_I4 = 4*NodeId[m_I]; \ 2744 m_I3 = 3*m_I; \ 2745 glcolar[m_I4] = m_of * glcolar[m_I4] + opacity * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2746 glcolar[m_I4] = m_of * glcolar[m_I4] + opacity * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2747 glcolar[m_I4] = m_of * glcolar[m_I4] + opacity * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2748 }\ 2749 }\ 2750 }\ 2751 } 2752 2753 #define SUMA_RGB_PGnL_AR4op2(RGBmat, NodeId, glcolar, nrgb, isColored, opacity, N_Node) {\ 2754 int m_I, m_I4 = 0, m_II; \ 2755 float m_of;\ 2756 m_of = (1 - opacity); \ 2757 for (m_I=0; m_I < nrgb; ++m_I) {\ 2758 if (NodeId[m_I] < N_Node) {\ 2759 if (!isColored[NodeId[m_I]]) { /* a new color, put it down as it is, Opacity gain should not be applied Wed Apr 2 17:31:33 EST 2003*/\ 2760 m_I4 = 4*NodeId[m_I]; \ 2761 glcolar[m_I4] = RGBmat[m_I][0]; ++m_I4;\ 2762 glcolar[m_I4] = RGBmat[m_I][1]; ++m_I4;\ 2763 glcolar[m_I4] = RGBmat[m_I][2]; ++m_I4;\ 2764 isColored[NodeId[m_I]] = YUP;\ 2765 }else {/* mixing to be done */\ 2766 m_I4 = 4*NodeId[m_I]; \ 2767 glcolar[m_I4] = m_of * glcolar[m_I4] + RGBmat[m_I][0]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4;\ 2768 glcolar[m_I4] = m_of * glcolar[m_I4] + RGBmat[m_I][1]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4;\ 2769 glcolar[m_I4] = m_of * glcolar[m_I4] + RGBmat[m_I][2]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4;\ 2770 }\ 2771 }\ 2772 }\ 2773 } 2774 #define SUMA_RGBv_PGnL_AR4op2(RGBvec, NodeId, glcolar, nrgb, isColored, opacity, N_Node) {\ 2775 int m_I, m_I4 = 0, m_II, m_I3=0; \ 2776 float m_of;\ 2777 m_of = (1 - opacity); \ 2778 for (m_I=0; m_I < nrgb; ++m_I) {\ 2779 if (NodeId[m_I] < N_Node) {\ 2780 if (!isColored[NodeId[m_I]]) { /* a new color, put it down as it is, Opacity gain should not be applied Wed Apr 2 17:31:33 EST 2003*/\ 2781 m_I4 = 4*NodeId[m_I]; \ 2782 m_I3 = 3*m_I; \ 2783 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2784 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2785 glcolar[m_I4] = RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2786 isColored[NodeId[m_I]] = YUP;\ 2787 }else {/* mixing to be done */\ 2788 m_I4 = 4*NodeId[m_I]; \ 2789 m_I3 = 3*m_I; \ 2790 glcolar[m_I4] = m_of * glcolar[m_I4] + RGBvec[m_I3]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4; ++m_I3;\ 2791 glcolar[m_I4] = m_of * glcolar[m_I4] + RGBvec[m_I3]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4; ++m_I3;\ 2792 glcolar[m_I4] = m_of * glcolar[m_I4] + RGBvec[m_I3]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4; ++m_I3;\ 2793 }\ 2794 }\ 2795 }\ 2796 } 2797 /*! 2798 This macro used to be called: SUMA_RGBmat_PartGlobLoc2_GLCOLAR4_opacity 2799 2800 SUMA_RGB_PGL_AR4op 2801 copies an N x 3 RGB matrix into a 4N x 1 GL color array format 2802 P (Part) means that colors are specified for some of the nodes only N < N_Nodes 2803 G (Glob) means an opacity is applied to the color values 2804 L (Local) means a local gain (per node) is applied to the color values 2805 2806 SUMA_RGB_PGL_AR4op(RGBmat, NodeId, glcolar, N, isColored, opacity, locgain, add, N_Nodes) 2807 RGBmat (float **) N x 3 matrix of RGB values 2808 NodeId (int *) N x 1 vector containing indices of nodes for wich color is specified in RGBmat 2809 glcolar (GLfloat *) (4 N_Nodes) x 1 vector 2810 isColored (SUMA_Boolean) N_Nodes x 1 vector indicating that a node was colored. ONLY YUP/1 are placed 2811 when a node is assigned a color. Values of isColored for nodes that have not been visited remain unchanged 2812 opacity (float) an opacity factor applied to each color R, G, B values in the entire list before adding it 2813 to a pre-exising color opacity is not applied to the color of nodes that had not been colored thus far. 2814 locgain (float *) N x 1 vector of gains applied to their respective nodes 2815 2816 *** Dec. 04 03: 2817 Variant SUMA_RGB_PGL_AR4op2 2818 Instead of mixing like this: Col_new = (1 - opacity) *Col_1 + opacity * locgain *Col_2; 2819 I am using Col_new = (1 - opacity)*Col_1 + locgain * Col_2; if (Col_new > 1) Col_new = 1 2820 *** Mar 17 04: 2821 Added the SUMA_RGBv versions using RGBvec instead of RGBmat. 2822 2823 */ 2824 2825 #define SUMA_RGB_PGL_AR4op(RGBmat, NodeId, glcolar, nrgb, isColored, opacity, locgain, N_Node) {\ 2826 int m_I, m_I4 = 0; \ 2827 float m_of, m_of2;\ 2828 for (m_I=0; m_I < nrgb; ++m_I) {\ 2829 if (NodeId[m_I] < N_Node) {\ 2830 if (!isColored[NodeId[m_I]]) { /* a new color, put it down as it is */\ 2831 m_I4 = 4*NodeId[m_I]; \ 2832 glcolar[m_I4] = locgain[m_I] * RGBmat[m_I][0]; ++m_I4;\ 2833 glcolar[m_I4] = locgain[m_I] * RGBmat[m_I][1]; ++m_I4;\ 2834 glcolar[m_I4] = locgain[m_I] * RGBmat[m_I][2]; ++m_I4;\ 2835 isColored[NodeId[m_I]] = YUP;\ 2836 }else { /* mixing to be done */\ 2837 m_I4 = 4*NodeId[m_I]; \ 2838 m_of = (locgain[m_I] * opacity); /* Dec. 04 03, changed from locgain[m_I] * opacity */\ 2839 m_of2 = (1 - opacity);\ 2840 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBmat[m_I][0]; ++m_I4;\ 2841 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBmat[m_I][1]; ++m_I4;\ 2842 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBmat[m_I][2]; ++m_I4;\ 2843 }\ 2844 }\ 2845 }\ 2846 } 2847 #define SUMA_RGBv_PGL_AR4op(RGBvec, NodeId, glcolar, nrgb, isColored, opacity, locgain, N_Node) {\ 2848 int m_I, m_I4 = 0, m_I3=0; \ 2849 float m_of, m_of2;\ 2850 for (m_I=0; m_I < nrgb; ++m_I) {\ 2851 if (NodeId[m_I] < N_Node) {\ 2852 if (!isColored[NodeId[m_I]]) { /* a new color, put it down as it is */\ 2853 m_I4 = 4*NodeId[m_I]; \ 2854 m_I3 = 3*m_I; \ 2855 glcolar[m_I4] = locgain[m_I] * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2856 glcolar[m_I4] = locgain[m_I] * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2857 glcolar[m_I4] = locgain[m_I] * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2858 isColored[NodeId[m_I]] = YUP;\ 2859 }else { /* mixing to be done */\ 2860 m_I4 = 4*NodeId[m_I]; \ 2861 m_I3 = 3*m_I; \ 2862 m_of = (locgain[m_I] * opacity); /* Dec. 04 03, changed from locgain[m_I] * opacity */\ 2863 m_of2 = (1 - opacity);\ 2864 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2865 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2866 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2867 }\ 2868 }\ 2869 }\ 2870 } 2871 2872 #define SUMA_RGB_PGL_AR4op2(RGBmat, NodeId, glcolar, nrgb, isColored, opacity, locgain, N_Node) {\ 2873 int m_I, m_I4 = 0; \ 2874 float m_of, m_of2;\ 2875 for (m_I=0; m_I < nrgb; ++m_I) {\ 2876 if (NodeId[m_I] < N_Node) {\ 2877 if (!isColored[NodeId[m_I]]) { /* a new color, put it down as it is */\ 2878 m_I4 = 4*NodeId[m_I]; \ 2879 glcolar[m_I4] = locgain[m_I] * RGBmat[m_I][0]; ++m_I4;\ 2880 glcolar[m_I4] = locgain[m_I] * RGBmat[m_I][1]; ++m_I4;\ 2881 glcolar[m_I4] = locgain[m_I] * RGBmat[m_I][2]; ++m_I4;\ 2882 isColored[NodeId[m_I]] = YUP;\ 2883 }else { /* mixing to be done */\ 2884 m_I4 = 4*NodeId[m_I]; \ 2885 m_of = (locgain[m_I]); /* Dec. 04 03, changed from locgain[m_I] * opacity */\ 2886 m_of2 = (1 - opacity);\ 2887 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBmat[m_I][0]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4;\ 2888 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBmat[m_I][1]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4;\ 2889 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBmat[m_I][2]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4;\ 2890 }\ 2891 }\ 2892 }\ 2893 } 2894 #define SUMA_RGBv_PGL_AR4op2(RGBvec, NodeId, glcolar, nrgb, isColored, opacity, locgain, N_Node) {\ 2895 int m_I, m_I4 = 0, m_I3=0; \ 2896 float m_of, m_of2;\ 2897 for (m_I=0; m_I < nrgb; ++m_I) {\ 2898 if (NodeId[m_I] < N_Node) {\ 2899 if (!isColored[NodeId[m_I]]) { /* a new color, put it down as it is */\ 2900 m_I4 = 4*NodeId[m_I]; \ 2901 m_I3 = 3*m_I; \ 2902 glcolar[m_I4] = locgain[m_I] * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2903 glcolar[m_I4] = locgain[m_I] * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2904 glcolar[m_I4] = locgain[m_I] * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2905 isColored[NodeId[m_I]] = YUP;\ 2906 }else { /* mixing to be done */\ 2907 m_I4 = 4*NodeId[m_I]; \ 2908 m_I3 = 3*m_I; \ 2909 m_of = (locgain[m_I]); /* Dec. 04 03, changed from locgain[m_I] * opacity */\ 2910 m_of2 = (1 - opacity);\ 2911 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBvec[m_I3]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4; ++m_I3;\ 2912 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBvec[m_I3]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4; ++m_I3;\ 2913 glcolar[m_I4] = m_of2 * glcolar[m_I4] + m_of * RGBvec[m_I3]; SUMA_CLIP_UB(glcolar[m_I4], 1.0); ++m_I4; ++m_I3;\ 2914 }\ 2915 }\ 2916 }\ 2917 } 2918 2919 /*! 2920 This macro used to be called: SUMA_RGBmat_PartNoGlobLoc2_GLCOLAR4_opacity 2921 SUMA_RGB_PnGL_AR4op 2922 copies an N x 3 RGB matrix into a 4N x 1 GL color array format 2923 P (Part) means that colors are specified for some of the nodes only N < N_Nodes 2924 nG (NoGlob) means no opacity is applied to the color values (fully opaque) 2925 L (Local) means a local gain (per node) is applied to the color values 2926 2927 SUMA_RGB_PnGL_AR4op(RGBmat, NodeId, glcolar, N, isColored, locgain, add, N_Nodes) 2928 RGBmat (float **) N x 3 matrix of RGB values 2929 NodeId (int *) N x 1 vector containing indices of nodes for wich color is specified in RGBmat 2930 glcolar (GLfloat *) (4 N_Nodes) x 1 vector 2931 isColored (SUMA_Boolean) N_Nodes x 1 vector indicating that a node was colored. ONLY YUP/1 are placed 2932 when a node is assigned a color. Values of isColored for nodes that have not been visited remain unchanged 2933 locgain (float *) N x 1 vector of gains applied to their respective nodes 2934 2935 *** Mar 17 04: 2936 Added the SUMA_RGBv versions using RGBvec instead of RGBmat. 2937 */ 2938 #define SUMA_RGB_PnGL_AR4op(RGBmat, NodeId, glcolar, nrgb, isColored, locgain, N_Node) {\ 2939 int m_I, m_I4 = 0; \ 2940 float m_of;\ 2941 for (m_I=0; m_I < nrgb; ++m_I) {\ 2942 if (NodeId[m_I] < N_Node) {\ 2943 m_I4 = 4*NodeId[m_I]; \ 2944 glcolar[m_I4] = locgain[m_I] * RGBmat[m_I][0]; ++m_I4;\ 2945 glcolar[m_I4] = locgain[m_I] * RGBmat[m_I][1]; ++m_I4;\ 2946 glcolar[m_I4] = locgain[m_I] * RGBmat[m_I][2]; ++m_I4;\ 2947 isColored[NodeId[m_I]] = YUP;\ 2948 }\ 2949 }\ 2950 } 2951 #define SUMA_RGBv_PnGL_AR4op(RGBvec, NodeId, glcolar, nrgb, isColored, locgain, N_Node) {\ 2952 int m_I, m_I4 = 0, m_I3=0; \ 2953 float m_of;\ 2954 for (m_I=0; m_I < nrgb; ++m_I) {\ 2955 if (NodeId[m_I] < N_Node) {\ 2956 m_I4 = 4*NodeId[m_I]; \ 2957 m_I3 = 3*m_I; \ 2958 glcolar[m_I4] = locgain[m_I] * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2959 glcolar[m_I4] = locgain[m_I] * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2960 glcolar[m_I4] = locgain[m_I] * RGBvec[m_I3]; ++m_I4; ++m_I3;\ 2961 isColored[NodeId[m_I]] = YUP;\ 2962 }\ 2963 }\ 2964 } 2965 2966 /*! 2967 Change a 12 element 1D affine transform vector to a 3x4 double array format 2968 */ 2969 #define SUMA_Xform1Dto2D(v,m){\ 2970 m[0][0] = v[0]; m[0][1] = v[1]; m[0][2] = v[2]; m[0][3] = v[3]; \ 2971 m[1][0] = v[4]; m[1][1] = v[5]; m[1][2] = v[6]; m[1][3] = v[7]; \ 2972 m[2][0] = v[8]; m[2][1] = v[9]; m[2][2] = v[10]; m[2][3] = v[11]; \ 2973 } 2974 2975 2976 2977 #endif 2978 2979