|
From: Gordon K. <kin...@us...> - 2004-03-25 03:22:34
|
Update of /cvsroot/teem/teem/src/coil In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv12209 Modified Files: coil.h coreCoil.c enumsCoil.c methodsCoil.c realmethods.c scalarCoil.c tensorCoil.c Log Message: some new stuff Index: coil.h =================================================================== RCS file: /cvsroot/teem/teem/src/coil/coil.h,v retrieving revision 1.2 retrieving revision 1.3 diff -C2 -d -r1.2 -r1.3 *** coil.h 24 Mar 2004 06:48:50 -0000 1.2 --- coil.h 25 Mar 2004 03:11:50 -0000 1.3 *************** *** 204,207 **** --- 204,208 ---- /* scalarCoil.c */ TEEM_API const coilKind *coilKindScalar; + TEEM_API const coilKind *coilKindArray[COIL_KIND_TYPE_MAX+1]; /* tensorCoil.c */ Index: tensorCoil.c =================================================================== RCS file: /cvsroot/teem/teem/src/coil/tensorCoil.c,v retrieving revision 1.2 retrieving revision 1.3 diff -C2 -d -r1.2 -r1.3 *** tensorCoil.c 24 Mar 2004 06:48:50 -0000 1.2 --- tensorCoil.c 25 Mar 2004 03:11:50 -0000 1.3 *************** *** 21,24 **** --- 21,45 ---- void + _coilKind7TensorTangents(coil_t traceGrad[6], + coil_t varianceGrad[6], + coil_t skewGrad[6], + coil_t rot0Grad[6], + coil_t rot1Grad[6], + coil_t rot2Grad[6], + coil_t tensor[7]) { + /* + coil_t a, b, c, d, e, f; + + a = tensor[1]; + b = tensor[2]; + c = tensor[3]; + d = tensor[4]; + e = tensor[5]; + f = tensor[6]; + ELL_6V_SET(traceGrad, 1, 0, 0, 1, 0, 1); + */ + } + + void _coilKind7TensorFilterTesting(coil_t *delta, coil_t **iv3, double spacing[3], *************** *** 33,45 **** } void _coilKind7TensorFilterHomogeneous(coil_t *delta, coil_t **iv3, double spacing[3], double parm[COIL_PARMS_NUM]) { ! } void ! _coilKindTensorUpdate(coil_t *val, coil_t *delta) { val[0] += delta[0]; /* WARNING: this could change confidence! */ --- 54,97 ---- } + /* + ** x ----> X + ** \ [0][ 0]-[06] [1][ 0] [2][ 0] + ** | \ [0][ 7]-[13] [1][ 7] [2][ 7] + ** | Y [0][14]-[20] [1][14] [2][14] + ** | + ** | [0][21]-[27] [1][21] [2][21] + ** Z [0][28]-[34] [1][28] [2][28] + ** [0][35]-[41] [1][35] [2][35] + ** + ** [0][42]-[48] [1][42] [2][42] + ** [0][49]-[55] [1][49] [2][49] + ** [0][56]-[62] [1][56] [2][56] + */ + + #define LAPL(iv3, vi, rspX, rspY, rspZ) \ + ( rspX*(iv3[0][vi + 7*4] - 2*iv3[1][vi + 7*4] + iv3[2][vi + 7*4]) \ + + rspY*(iv3[1][vi + 7*3] - 2*iv3[1][vi + 7*4] + iv3[1][vi + 7*5]) \ + + rspZ*(iv3[1][vi + 7*1] - 2*iv3[1][vi + 7*4] + iv3[1][vi + 7*7])) + void _coilKind7TensorFilterHomogeneous(coil_t *delta, coil_t **iv3, double spacing[3], double parm[COIL_PARMS_NUM]) { ! coil_t rspX, rspY, rspZ; ! ! rspX = 1.0/(spacing[0]*spacing[0]); ! rspY = 1.0/(spacing[1]*spacing[1]); ! rspZ = 1.0/(spacing[2]*spacing[2]); ! delta[0] = 0; ! delta[1] = parm[0]*LAPL(iv3, 1, rspX, rspY, rspZ); ! delta[2] = parm[0]*LAPL(iv3, 2, rspX, rspY, rspZ); ! delta[3] = parm[0]*LAPL(iv3, 3, rspX, rspY, rspZ); ! delta[4] = parm[0]*LAPL(iv3, 4, rspX, rspY, rspZ); ! delta[5] = parm[0]*LAPL(iv3, 5, rspX, rspY, rspZ); ! delta[6] = parm[0]*LAPL(iv3, 6, rspX, rspY, rspZ); } void ! _coilKind7TensorUpdate(coil_t *val, coil_t *delta) { val[0] += delta[0]; /* WARNING: this could change confidence! */ *************** *** 53,57 **** const coilKind ! _coilKindTensor = { "tensor", 7, --- 105,109 ---- const coilKind ! _coilKind7Tensor = { "tensor", 7, *************** *** 62,68 **** NULL, NULL}, ! _coilKindTensorUpdate }; const coilKind * ! coilKindTensor = &_coilKindTensor; --- 114,120 ---- NULL, NULL}, ! _coilKind7TensorUpdate }; const coilKind * ! coilKind7Tensor = &_coilKind7Tensor; Index: scalarCoil.c =================================================================== RCS file: /cvsroot/teem/teem/src/coil/scalarCoil.c,v retrieving revision 1.2 retrieving revision 1.3 diff -C2 -d -r1.2 -r1.3 *** scalarCoil.c 24 Mar 2004 06:48:50 -0000 1.2 --- scalarCoil.c 25 Mar 2004 03:11:50 -0000 1.3 *************** *** 20,23 **** --- 20,38 ---- #include "coil.h" + /* + ** x ----> X + ** \ [0][0] [1][0] [2][0] + ** | \ [0][1] [1][1] [2][1] + ** | Y [0][2] [1][2] [2][2] + ** | + ** | [0][3] [1][3] [2][3] + ** Z [0][4] [1][4] [2][4] + ** [0][5] [1][5] [2][5] + ** + ** [0][6] [1][6] [2][6] + ** [0][7] [1][7] [2][7] + ** [0][8] [1][8] [2][8] + */ + coil_t _coilLaplacian3(coil_t **iv3, double spacing[3]) { *************** *** 43,72 **** } - /* - ** x ----> X - ** \ [0][0] [1][0] [2][0] - ** | \ [0][1] [1][1] [2][1] - ** | Y [0][2] [1][2] [2][2] - ** | - ** | [0][3] [1][3] [2][3] - ** Z [0][4] [1][4] [2][4] - ** [0][5] [1][5] [2][5] - ** - ** [0][6] [1][6] [2][6] - ** [0][7] [1][7] [2][7] - ** [0][8] [1][8] [2][8] - */ - void ! _coilKindScalarFilterPeronaMalik(coil_t *delta, coil_t **i, ! double spacing[3], ! double parm[COIL_PARMS_NUM]) { ! coil_t forwX[3], backX[3], forwY[3], backY[3], forwZ[3], backZ[3], ! KK, rspX, rspY, rspZ; ! ! /* reciprocals of spacings in X, Y, and Z */ ! rspX = 1.0/spacing[0]; ! rspY = 1.0/spacing[1]; ! rspZ = 1.0/spacing[2]; /* gradients at forward and backward X */ --- 58,67 ---- } void ! _coilKindScalar3x3x3Gradients(coil_t *forwX, coil_t *backX, ! coil_t *forwY, coil_t *backY, ! coil_t *forwZ, coil_t *backZ, ! coil_t **i, ! coil_t rspX, coil_t rspY, coil_t rspZ) { /* gradients at forward and backward X */ *************** *** 94,113 **** backZ[2] = rspZ*(i[1][4] - i[1][1]); /* compute fluxes */ KK = parm[1]*parm[1]; ! /* ! forwX[0] *= 1.0/(1.0 + ELL_3V_DOT(forwX, forwX)/KK); ! forwY[1] *= 1.0/(1.0 + ELL_3V_DOT(forwY, forwY)/KK); ! forwZ[2] *= 1.0/(1.0 + ELL_3V_DOT(forwZ, forwZ)/KK); ! backX[0] *= 1.0/(1.0 + ELL_3V_DOT(backX, backX)/KK); ! backY[1] *= 1.0/(1.0 + ELL_3V_DOT(backY, backY)/KK); ! backZ[2] *= 1.0/(1.0 + ELL_3V_DOT(backZ, backZ)/KK); ! */ ! forwX[0] *= exp(-0.5*ELL_3V_DOT(forwX, forwX)/KK); ! forwY[1] *= exp(-0.5*ELL_3V_DOT(forwY, forwY)/KK); ! forwZ[2] *= exp(-0.5*ELL_3V_DOT(forwZ, forwZ)/KK); ! backX[0] *= exp(-0.5*ELL_3V_DOT(backX, backX)/KK); ! backY[1] *= exp(-0.5*ELL_3V_DOT(backY, backY)/KK); ! backZ[2] *= exp(-0.5*ELL_3V_DOT(backZ, backZ)/KK); delta[0] = parm[0]*(rspX*(forwX[0] - backX[0]) --- 89,129 ---- backZ[2] = rspZ*(i[1][4] - i[1][1]); + return; + } + + #define _COIL_CONDUCT(LL, KK) \ + (exp(-0.5*(LL)/(KK))) + + /* + #define _COIL_CONDUCT(vec, KK) \ + (1.0/(1.0 + (LL)/(KK))) + */ + + void + _coilKindScalarFilterPeronaMalik(coil_t *delta, coil_t **iv3, + double spacing[3], + double parm[COIL_PARMS_NUM]) { + coil_t forwX[3], backX[3], forwY[3], backY[3], forwZ[3], backZ[3], + KK, rspX, rspY, rspZ; + + /* reciprocals of spacings in X, Y, and Z */ + rspX = 1.0/spacing[0]; + rspY = 1.0/spacing[1]; + rspZ = 1.0/spacing[2]; + + _coilKindScalar3x3x3Gradients(forwX, backX, + forwY, backY, + forwZ, backZ, + iv3, + rspX, rspY, rspZ); + /* compute fluxes */ KK = parm[1]*parm[1]; ! forwX[0] *= _COIL_CONDUCT(ELL_3V_DOT(forwX, forwX), KK); ! forwY[1] *= _COIL_CONDUCT(ELL_3V_DOT(forwY, forwY), KK); ! forwZ[2] *= _COIL_CONDUCT(ELL_3V_DOT(forwZ, forwZ), KK); ! backX[0] *= _COIL_CONDUCT(ELL_3V_DOT(backX, backX), KK); ! backY[1] *= _COIL_CONDUCT(ELL_3V_DOT(backY, backY), KK); ! backZ[2] *= _COIL_CONDUCT(ELL_3V_DOT(backZ, backZ), KK); delta[0] = parm[0]*(rspX*(forwX[0] - backX[0]) *************** *** 117,120 **** --- 133,179 ---- void + _coilKindScalarFilterModifiedCurvature(coil_t *delta, coil_t **iv3, + double spacing[3], + double parm[COIL_PARMS_NUM]) { + coil_t forwX[3], backX[3], forwY[3], backY[3], forwZ[3], backZ[3], + grad[3], gm, eps, KK, LL, rspX, rspY, rspZ; + + /* reciprocals of spacings in X, Y, and Z */ + rspX = 1.0/spacing[0]; + rspY = 1.0/spacing[1]; + rspZ = 1.0/spacing[2]; + + _coilKindScalar3x3x3Gradients(forwX, backX, + forwY, backY, + forwZ, backZ, + iv3, + rspX, rspY, rspZ); + grad[0] = rspX*(iv3[2][4] - iv3[0][4]); + grad[1] = rspY*(iv3[1][5] - iv3[1][3]); + grad[2] = rspZ*(iv3[1][7] - iv3[1][1]); + gm = ELL_3V_LEN(grad); + + /* compute fluxes */ + eps = 0.000001; + KK = parm[1]*parm[1]; + LL = ELL_3V_DOT(forwX, forwX); + forwX[0] *= _COIL_CONDUCT(LL, KK)/(eps + sqrt(LL)); + LL = ELL_3V_DOT(forwY, forwY); + forwY[1] *= _COIL_CONDUCT(LL, KK)/(eps + sqrt(LL)); + LL = ELL_3V_DOT(forwZ, forwZ); + forwZ[2] *= _COIL_CONDUCT(LL, KK)/(eps + sqrt(LL)); + LL = ELL_3V_DOT(backX, backX); + backX[0] *= _COIL_CONDUCT(LL, KK)/(eps + sqrt(LL)); + LL = ELL_3V_DOT(backY, backY); + backY[1] *= _COIL_CONDUCT(LL, KK)/(eps + sqrt(LL)); + LL = ELL_3V_DOT(backZ, backZ); + backZ[2] *= _COIL_CONDUCT(LL, KK)/(eps + sqrt(LL)); + + delta[0] = gm*parm[0]*(rspX*(forwX[0] - backX[0]) + + rspY*(forwY[1] - backY[1]) + + rspZ*(forwZ[2] - backZ[2])); + } + + void _coilKindScalarUpdate(coil_t *val, coil_t *delta) { *************** *** 130,134 **** _coilKindScalarFilterHomogeneous, _coilKindScalarFilterPeronaMalik, ! NULL, NULL}, _coilKindScalarUpdate --- 189,193 ---- _coilKindScalarFilterHomogeneous, _coilKindScalarFilterPeronaMalik, ! _coilKindScalarFilterModifiedCurvature, NULL}, _coilKindScalarUpdate *************** *** 137,138 **** --- 196,209 ---- const coilKind * coilKindScalar = &_coilKindScalar; + + /* ------------------------------------------ */ + + extern const coilKind _coilKind7Tensor; + + const coilKind* + coilKindArray[COIL_KIND_TYPE_MAX+1] = { + NULL, + &_coilKindScalar, + NULL, + &_coilKind7Tensor + }; Index: realmethods.c =================================================================== RCS file: /cvsroot/teem/teem/src/coil/realmethods.c,v retrieving revision 1.3 retrieving revision 1.4 diff -C2 -d -r1.3 -r1.4 *** realmethods.c 24 Mar 2004 06:48:50 -0000 1.3 --- realmethods.c 25 Mar 2004 03:11:50 -0000 1.4 *************** *** 20,23 **** --- 20,24 ---- #include "coil.h" + /* ------------------------------------------ */ const coilMethod _coilMethodTesting = { *************** *** 29,32 **** --- 30,34 ---- coilMethodTesting = &_coilMethodTesting; + /* ------------------------------------------ */ const coilMethod _coilMethodHomogeneous = { *************** *** 38,41 **** --- 40,44 ---- coilMethodHomogeneous = &_coilMethodHomogeneous; + /* ------------------------------------------ */ const coilMethod _coilMethodPeronaMalik = { *************** *** 47,50 **** --- 50,64 ---- coilMethodPeronaMalik = &_coilMethodPeronaMalik; + /* ------------------------------------------ */ + const coilMethod + _coilMethodModifiedCurvature = { + "modified-curvature", + coilMethodTypeModifiedCurvature, + 2 + }; + const coilMethod* + coilMethodModifiedCurvature = &_coilMethodModifiedCurvature; + + /* ------------------------------------------ */ const coilMethod* coilMethodArray[COIL_METHOD_TYPE_MAX+1] = { *************** *** 53,57 **** &_coilMethodHomogeneous, &_coilMethodPeronaMalik, ! NULL, NULL }; --- 67,71 ---- &_coilMethodHomogeneous, &_coilMethodPeronaMalik, ! &_coilMethodModifiedCurvature, NULL }; Index: methodsCoil.c =================================================================== RCS file: /cvsroot/teem/teem/src/coil/methodsCoil.c,v retrieving revision 1.3 retrieving revision 1.4 diff -C2 -d -r1.3 -r1.4 *** methodsCoil.c 24 Mar 2004 06:48:50 -0000 1.3 --- methodsCoil.c 25 Mar 2004 03:11:50 -0000 1.4 *************** *** 191,194 **** --- 191,195 ---- coilOutputGet(Nrrd *nout, coilContext *cctx) { char me[]="coilOutputGet", err[AIR_STRLEN_MED]; + int baseDim; if (!(nout && cctx)) { *************** *** 196,200 **** biffAdd(COIL, err); return 1; } ! if (nrrdSlice(nout, cctx->nvol, 0, 0)) { sprintf(err, "%s: trouble slicing to get output", me); biffMove(COIL, err, NRRD); return 1; --- 197,202 ---- biffAdd(COIL, err); return 1; } ! baseDim = (1 == cctx->kind->valLen ? 0 : 1); ! if (nrrdSlice(nout, cctx->nvol, baseDim, 0)) { sprintf(err, "%s: trouble slicing to get output", me); biffMove(COIL, err, NRRD); return 1; Index: enumsCoil.c =================================================================== RCS file: /cvsroot/teem/teem/src/coil/enumsCoil.c,v retrieving revision 1.2 retrieving revision 1.3 diff -C2 -d -r1.2 -r1.3 *** enumsCoil.c 24 Mar 2004 06:48:50 -0000 1.2 --- enumsCoil.c 25 Mar 2004 03:11:50 -0000 1.3 *************** *** 93,97 **** "scalar", "3color", ! "7tensor", "" }; --- 93,97 ---- "scalar", "3color", ! "7tensor", "tensor", "" }; *************** *** 101,105 **** coilKindTypeScalar, coilKindType3Color, ! coilKindType7Tensor }; --- 101,105 ---- coilKindTypeScalar, coilKindType3Color, ! coilKindType7Tensor, coilKindType7Tensor }; Index: coreCoil.c =================================================================== RCS file: /cvsroot/teem/teem/src/coil/coreCoil.c,v retrieving revision 1.3 retrieving revision 1.4 diff -C2 -d -r1.3 -r1.4 *** coreCoil.c 24 Mar 2004 06:48:50 -0000 1.3 --- coreCoil.c 25 Mar 2004 03:11:50 -0000 1.4 *************** *** 89,92 **** --- 89,104 ---- } + void + _coilIv3Fill_1_7(coil_t **iv3, coil_t *here, int radius, int valLen, + int x0, int y0, int z0, int sizeX, int sizeY, int sizeZ) { + int vi, /* value index */ + xni, yni, zni, /* neighborhood (iv3) indices */ + xvi, yvi, zvi; /* volume indices */ + coil_t *tmp; + + _COIL_IV3_FILL(1, 3, 7); + return; + } + int _coilThisZGet(coilTask *task, int doFilter) { *************** *** 188,191 **** --- 200,205 ---- if (1 == cctx->radius && 1 == cctx->kind->valLen) { task->iv3Fill = _coilIv3Fill_1_1; + } else if (1 == cctx->radius && 7 == cctx->kind->valLen) { + task->iv3Fill = _coilIv3Fill_1_7; } else { task->iv3Fill = _coilIv3Fill_R_L; |