|
From: Gordon K. <kin...@us...> - 2004-05-07 19:13:44
|
Update of /cvsroot/teem/teem/src/ten In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv9915 Modified Files: fiber.c tenGage.c tenMacros.h tendHelix.c tensor.c Log Message: changed macros, and started adding shape gradients. Should not break any code Index: fiber.c =================================================================== RCS file: /cvsroot/teem/teem/src/ten/fiber.c,v retrieving revision 1.11 retrieving revision 1.12 diff -C2 -d -r1.11 -r1.12 *** fiber.c 7 Jan 2004 15:34:31 -0000 1.11 --- fiber.c 7 May 2004 19:13:33 -0000 1.12 *************** *** 103,107 **** if (ELL_3V_DOT(tfx->lastDir, tfx->lastDir)) { ELL_3V_COPY(vin, tfx->lastDir); ! TEN_3VLIST_MUL(vout, tfx->dten, tfx->lastDir); ELL_3V_NORM(vout, vout, len); _tenFiberAlign(tfx, vout); /* HEY: is this needed? */ --- 103,107 ---- if (ELL_3V_DOT(tfx->lastDir, tfx->lastDir)) { ELL_3V_COPY(vin, tfx->lastDir); ! TEN_T3V_MUL(vout, tfx->dten, tfx->lastDir); ELL_3V_NORM(vout, vout, len); _tenFiberAlign(tfx, vout); /* HEY: is this needed? */ Index: tensor.c =================================================================== RCS file: /cvsroot/teem/teem/src/ten/tensor.c,v retrieving revision 1.29 retrieving revision 1.30 diff -C2 -d -r1.29 -r1.30 *** tensor.c 13 Mar 2004 20:03:11 -0000 1.29 --- tensor.c 7 May 2004 19:13:34 -0000 1.30 *************** *** 103,107 **** continue; } ! TEN_LIST2MAT(nine, seven); ELL_3M_SCALE(nine, scale, nine); } --- 103,107 ---- continue; } ! TEN_T2M(nine, seven); ELL_3M_SCALE(nine, scale, nine); } *************** *** 163,167 **** N = sx*sy*sz; for (I=0; I<N; I++) { ! TEN_MAT2LIST(seven, nine); seven[0] = conf ? conf[I] : 1.0; seven += 7; --- 163,167 ---- N = sx*sy*sz; for (I=0; I<N; I++) { ! TEN_M2T(seven, nine); seven[0] = conf ? conf[I] : 1.0; seven += 7; *************** *** 198,202 **** int ret; ! TEN_LIST2MAT(m, t); trc = ELL_3M_TRACE(m)/3.0; ELL_3M_IDENTITY_SET(iso); --- 198,202 ---- int ret; ! TEN_T2M(m, t); trc = ELL_3M_TRACE(m)/3.0; ELL_3M_IDENTITY_SET(iso); *************** *** 298,302 **** ELL_3M_MUL(tmpMat2, evecT, tmpMat1); ten[0] = conf; ! TEN_MAT2LIST(ten, tmpMat2); return; } --- 298,302 ---- ELL_3M_MUL(tmpMat2, evecT, tmpMat1); ten[0] = conf; ! TEN_M2T(ten, tmpMat2); return; } *************** *** 473,475 **** --- 473,516 ---- } + void + tenShapeGradients_d(double mean[7], + double var[7], double *varNorm, + double skew[7], double *skewNorm, + double evec[9], double eval[9], int *didEigen, + double ten[7]) { + double meanDot, varDot; + #define Txx (ten[1]) + #define Txy (ten[2]) + #define Txz (ten[3]) + #define Tyy (ten[4]) + #define Tyz (ten[5]) + #define Tzz (ten[6]) + /* 0.57735027 == 1/sqrt(3) */ + TEN_T_SET(mean, ten[0], + 0.57735027, 0, 0, + 0.57735027, 0, + 0.57735027); + TEN_T_SET(var, ten[0], + 2*Txx - Tyy - Tzz, 3*Txy, 3*Txz, + 2*Tyy - Txx - Tzz, 3*Tyz, + 2*Tzz - Txx - Tyy); + *varNorm = TEN_T_NORM(var); + if (*varNorm) { + TEN_T_SCALE(var, 1.0/(*varNorm), var); + } + TEN_T_SET(skew, ten[0], + Tyy*Tzz - Tyz*Tyz, Txz*Tyz - Txy*Tzz, Txy*Tyz - Txz*Tyy, + Txx*Tzz - Txz*Txz, Txy*Txz - Tyz*Txx, + Txx*Tyy - Txy*Txy); + meanDot = TEN_T_DOT(skew, mean); + varDot = TEN_T_DOT(skew, var); + TEN_T_SCALE_INCR2(skew, -meanDot, mean, -varDot, var); + *skewNorm = TEN_T_NORM(skew); + if (*skewNorm) { + TEN_T_SCALE(skew, 1.0/(*skewNorm), skew); + } + + *didEigen = AIR_FALSE; + return; + } Index: tendHelix.c =================================================================== RCS file: /cvsroot/teem/teem/src/ten/tendHelix.c,v retrieving revision 1.5 retrieving revision 1.6 diff -C2 -d -r1.5 -r1.6 *** tendHelix.c 19 Feb 2004 06:45:13 -0000 1.5 --- tendHelix.c 7 May 2004 19:13:34 -0000 1.6 *************** *** 121,125 **** ELL_3M_MUL(mA, H2W, mB); ! TEN_MAT2LIST(out + 7*idx, mA); (out + 7*idx)[0] = 1.0; } --- 121,125 ---- ELL_3M_MUL(mA, H2W, mB); ! TEN_M2T(out + 7*idx, mA); (out + 7*idx)[0] = 1.0; } Index: tenGage.c =================================================================== RCS file: /cvsroot/teem/teem/src/ten/tenGage.c,v retrieving revision 1.17 retrieving revision 1.18 diff -C2 -d -r1.17 -r1.18 *** tenGage.c 22 Apr 2004 07:14:54 -0000 1.17 --- tenGage.c 7 May 2004 19:13:33 -0000 1.18 *************** *** 18,21 **** --- 18,28 ---- */ + extern void + tenShapeGradients_d(double mean[7], + double var[7], double *varNorm, + double skew[7], double *skewNorm, + double evec[9], double eval[9], int *didEigen, + double ten[7]); + #include "ten.h" *************** *** 244,250 **** --- 251,287 ---- /* --- Trace --- */ if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTraceGradVec)) { + /* + double tmpT[7], meanGrad[7], varGrad[7], varNorm, skewGrad[7], skewNorm, + copyT[7]; + int didEigen=0; + + TEN_T_COPY(copyT, tenAns); + tenShapeGradients_d(meanGrad, + varGrad, &varNorm, + skewGrad, &skewNorm, + NULL, NULL, &didEigen, + copyT); + vecTmp = pvl->directAnswer[tenGageTraceGradVec]; + TEN_T_SET(tmpT, tenAns[0], + gradDtA[0], gradDtB[0], gradDtC[0], + gradDtD[0], gradDtE[0], + gradDtF[0]); + vecTmp[0] = TEN_T_DOT(meanGrad, tmpT); + TEN_T_SET(tmpT, tenAns[0], + gradDtA[1], gradDtB[1], gradDtC[1], + gradDtD[1], gradDtE[1], + gradDtF[1]); + vecTmp[1] = TEN_T_DOT(meanGrad, tmpT); + TEN_T_SET(tmpT, tenAns[0], + gradDtA[2], gradDtB[2], gradDtC[2], + gradDtD[2], gradDtE[2], + gradDtF[2]); + vecTmp[2] = TEN_T_DOT(meanGrad, tmpT); + */ + vecTmp = pvl->directAnswer[tenGageTraceGradVec]; ELL_3V_ADD3(vecTmp, gradDtA, gradDtD, gradDtF); ELL_3V_SCALE(gradCbA, -1, vecTmp); + } if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTraceGradMag)) { *************** *** 312,319 **** --- 349,386 ---- /* --- Q --- */ if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageQGradVec)) { + /* + double tmpT[7], meanGrad[7], varGrad[7], varNorm, skewGrad[7], skewNorm, + copyT[7]; + int didEigen=0; + + TEN_T_COPY(copyT, tenAns); + tenShapeGradients_d(meanGrad, + varGrad, &varNorm, + skewGrad, &skewNorm, + NULL, NULL, &didEigen, + copyT); + gradCbQ = vecTmp = pvl->directAnswer[tenGageQGradVec]; + TEN_T_SET(tmpT, tenAns[0], + gradDtA[0], gradDtB[0], gradDtC[0], + gradDtD[0], gradDtE[0], + gradDtF[0]); + vecTmp[0] = TEN_T_DOT(varGrad, tmpT); + TEN_T_SET(tmpT, tenAns[0], + gradDtA[1], gradDtB[1], gradDtC[1], + gradDtD[1], gradDtE[1], + gradDtF[1]); + vecTmp[1] = TEN_T_DOT(varGrad, tmpT); + TEN_T_SET(tmpT, tenAns[0], + gradDtA[2], gradDtB[2], gradDtC[2], + gradDtD[2], gradDtE[2], + gradDtF[2]); + vecTmp[2] = TEN_T_DOT(varGrad, tmpT); + */ + gradCbQ = vecTmp = pvl->directAnswer[tenGageQGradVec]; ELL_3V_SCALE_ADD2(vecTmp, 1.0/9.0, gradCbS, -1.0/9.0, gradCbB); + } if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageQGradMag)) { Index: tenMacros.h =================================================================== RCS file: /cvsroot/teem/teem/src/ten/tenMacros.h,v retrieving revision 1.10 retrieving revision 1.11 diff -C2 -d -r1.10 -r1.11 *** tenMacros.h 7 Jan 2004 15:34:31 -0000 1.10 --- tenMacros.h 7 May 2004 19:13:34 -0000 1.11 *************** *** 26,30 **** /* ! ******** TEN_LIST2MAT, TEN_MAT2LIST ** ** for going between 7-element list and 9-element matrix --- 26,30 ---- /* ! ******** TEN_T2M, TEN_M2T ** ** for going between 7-element list and 9-element matrix *************** *** 40,85 **** ** As in ell, the matrix ordering is given by: ** ! ** 0 3 6 ! ** 1 4 7 ! ** 2 5 8 ** ! ** Note that TEN_MAT2LIST does NOT set the threshold element (index 0), ! ** and that the threshold value plays no role in TEN_LIST2MAT. */ ! #define TEN_LIST2MAT(m, l) ( \ ! (m)[0] = (l)[1], \ ! (m)[1] = (l)[2], \ ! (m)[2] = (l)[3], \ ! (m)[3] = (l)[2], \ ! (m)[4] = (l)[4], \ ! (m)[5] = (l)[5], \ ! (m)[6] = (l)[3], \ ! (m)[7] = (l)[5], \ ! (m)[8] = (l)[6] ) ! #define TEN_MAT2LIST(l, m) ( \ ! (l)[1] = (m)[0], \ ! (l)[2] = (m)[3], \ ! (l)[3] = (m)[6], \ ! (l)[4] = (m)[4], \ ! (l)[5] = (m)[7], \ ! (l)[6] = (m)[8] ) ! #define TEN_LIST_COPY(d, s) ( \ ! (d)[0] = (s)[0], \ ! (d)[1] = (s)[1], \ ! (d)[2] = (s)[2], \ ! (d)[3] = (s)[3], \ ! (d)[4] = (s)[4], \ ! (d)[5] = (s)[5], \ (d)[6] = (s)[6] ) ! #define TEN_LIST_DET(l) ( \ ! (l)[1]*((l)[4]*(l)[6] - (l)[5]*(l)[5]) \ ! + (l)[2]*((l)[5]*(l)[3] - (l)[2]*(l)[6]) \ ! + (l)[3]*((l)[2]*(l)[5] - (l)[3]*(l)[4])) ! #define TEN_LIST_SCALE(a, s, b) ( \ (a)[0] = (b)[0], \ (a)[1] = (s)*(b)[1], \ --- 40,89 ---- ** As in ell, the matrix ordering is given by: ** ! ** 0 1 2 ! ** 3 4 5 ! ** 6 7 8 ** ! ** Note that TEN_M2T does NOT set the threshold element (index 0), ! ** and that the threshold value plays no role in TEN_T2M. */ ! #define TEN_T2M(m, t) ( \ ! (m)[0] = (t)[1], (m)[1] = (t)[2], (m)[2] = (t)[3], \ ! (m)[3] = (t)[2], (m)[4] = (t)[4], (m)[5] = (t)[5], \ ! (m)[6] = (t)[3], (m)[7] = (t)[5], (m)[8] = (t)[6] ) ! #define TEN_M2T(t, m) ( \ ! (t)[1] = (m)[0], (t)[2] = (m)[1], (t)[3] = (m)[2], \ ! (t)[4] = (m)[4], (t)[5] = (m)[5], \ ! (t)[6] = (m)[8] ) ! #define TEN_T_SET(t, conf, a, b, c, d, e, f) ( \ ! (t)[0] = (conf), \ ! (t)[1] = (a), (t)[2] = (b), (t)[3] = (c), \ ! (t)[4] = (d), (t)[5] = (e), \ ! (t)[6] = (f) ) ! ! #define TEN_T_COPY(d, s) ( \ ! (d)[0] = (s)[0], \ ! (d)[1] = (s)[1], \ ! (d)[2] = (s)[2], \ ! (d)[3] = (s)[3], \ ! (d)[4] = (s)[4], \ ! (d)[5] = (s)[5], \ (d)[6] = (s)[6] ) ! #define TEN_T_DET(t) ( \ ! (t)[1]*((t)[4]*(t)[6] - (t)[5]*(t)[5]) \ ! + (t)[2]*((t)[5]*(t)[3] - (t)[2]*(t)[6]) \ ! + (t)[3]*((t)[2]*(t)[5] - (t)[3]*(t)[4])) ! #define TEN_T_DOT(A, B) ( \ ! (A)[1]*(B)[1] + 2*(A)[2]*(B)[2] + 2*(A)[3]*(B)[3] \ ! + (A)[4]*(B)[4] + 2*(A)[5]*(B)[5] \ ! + (A)[6]*(B)[6] ) ! ! #define TEN_T_NORM(A) (sqrt(TEN_T_DOT(A,A))) ! ! #define TEN_T_SCALE(a, s, b) ( \ (a)[0] = (b)[0], \ (a)[1] = (s)*(b)[1], \ *************** *** 90,97 **** (a)[6] = (s)*(b)[6]) ! #define TEN_3VLIST_MUL(b, l, a) ( \ ! (b)[0] = (l)[1]*(a)[0] + (l)[2]*(a)[1] + (l)[3]*(a)[2], \ ! (b)[1] = (l)[2]*(a)[0] + (l)[4]*(a)[1] + (l)[5]*(a)[2], \ ! (b)[2] = (l)[3]*(a)[0] + (l)[5]*(a)[1] + (l)[6]*(a)[2]) #ifdef __cplusplus --- 94,119 ---- (a)[6] = (s)*(b)[6]) ! #define TEN_T_SCALE_INCR(a, s, b) ( \ ! (a)[0] = (b)[0], \ ! (a)[1] += (s)*(b)[1], \ ! (a)[2] += (s)*(b)[2], \ ! (a)[3] += (s)*(b)[3], \ ! (a)[4] += (s)*(b)[4], \ ! (a)[5] += (s)*(b)[5], \ ! (a)[6] += (s)*(b)[6]) ! ! #define TEN_T_SCALE_INCR2(a, s, b, t, c) ( \ ! (a)[0] = AIR_MIN((b)[0], c[0]), \ ! (a)[1] += (s)*(b)[1] + (t)*(c)[1], \ ! (a)[2] += (s)*(b)[2] + (t)*(c)[2], \ ! (a)[3] += (s)*(b)[3] + (t)*(c)[3], \ ! (a)[4] += (s)*(b)[4] + (t)*(c)[4], \ ! (a)[5] += (s)*(b)[5] + (t)*(c)[5], \ ! (a)[6] += (s)*(b)[6] + (t)*(c)[6]) ! ! #define TEN_T3V_MUL(b, t, a) ( \ ! (b)[0] = (t)[1]*(a)[0] + (t)[2]*(a)[1] + (t)[3]*(a)[2], \ ! (b)[1] = (t)[2]*(a)[0] + (t)[4]*(a)[1] + (t)[5]*(a)[2], \ ! (b)[2] = (t)[3]*(a)[0] + (t)[5]*(a)[1] + (t)[6]*(a)[2]) #ifdef __cplusplus |