From: <kin...@us...> - 2008-12-09 11:02:03
|
Revision: 4035 http://teem.svn.sourceforge.net/teem/?rev=4035&view=rev Author: kindlmann Date: 2008-12-09 11:02:01 +0000 (Tue, 09 Dec 2008) Log Message: ----------- Integrating last fixes and changes from Arish Qazi based on his two-tensor tractography work: * If the single tensor is too planar (based on pvlData->levmarMinCp), then the two-tensor fit will not be computed * removed "ELL_3V_NORM(tfx->lastDir, tfx->lastDir, len);" from _tenFiberProbe in the case of two-tensor work; the reason why this works is not yet clear.... * the measurement frame is applied to the whole tensor, instead of just its eigenvectors Modified Paths: -------------- teem/trunk/src/ten/fiber.c teem/trunk/src/ten/tenDwiGage.c Modified: teem/trunk/src/ten/fiber.c =================================================================== --- teem/trunk/src/ten/fiber.c 2008-12-08 17:38:33 UTC (rev 4034) +++ teem/trunk/src/ten/fiber.c 2008-12-09 11:02:01 UTC (rev 4035) @@ -55,6 +55,7 @@ char me[]="_tenFiberProbe", err[BIFF_STRLEN]; double iPos[3]; int ret = 0; + double tens2[2][7]; gageShapeWtoI(tfx->gtx->shape, iPos, wPos); *gageRet = gageProbe(tfx->gtx, iPos[0], iPos[1], iPos[2]); @@ -82,7 +83,7 @@ airEnumStr(tenDwiFiberType, tfx->fiberType)); } switch (tfx->fiberType) { - double evec[2][9], eval[2][3], len; + double evec[2][9], eval[2][3]; case tenDwiFiberType1Evec0: if (tfx->mframeUse) { double matTmpA[9], matTmpB[9]; @@ -106,18 +107,27 @@ break; case tenDwiFiberType2Evec0: /* Estimate principal diffusion direction of each tensor */ - tenEigensolve_d(eval[0], evec[0], tfx->gageTen2 + 0); - tenEigensolve_d(eval[1], evec[1], tfx->gageTen2 + 7); + if (tfx->mframeUse) { + /* Transform both the tensors */ + double matTmpA[9], matTmpB[9]; - if (tfx->mframeUse) { - double vectmp[3]; - /* NOTE: only fixing first (of three) eigenvector, for both tensors */ - ELL_3MV_MUL(vectmp, tfx->mframe, evec[0]); - ELL_3V_COPY(evec[0], vectmp); - ELL_3MV_MUL(vectmp, tfx->mframe, evec[1]); - ELL_3V_COPY(evec[1], vectmp); + TEN_T2M(matTmpA, tfx->gageTen2 + 0); + ELL_3M_MUL(matTmpB, tfx->mframe, matTmpA); + ELL_3M_MUL(matTmpA, matTmpB, tfx->mframeT); + TEN_M2T(tens2[0], matTmpA); + /* new eigen values and vectors */ + tenEigensolve_d(eval[0], evec[0], tens2[0]); + + TEN_T2M(matTmpA, tfx->gageTen2 + 7); + ELL_3M_MUL(matTmpB, tfx->mframe, matTmpA); + ELL_3M_MUL(matTmpA, matTmpB, tfx->mframeT); + TEN_M2T(tens2[1], matTmpA); + tenEigensolve_d(eval[1], evec[1], tens2[1]); + } else { + tenEigensolve_d(eval[0], evec[0], tfx->gageTen2 + 0); + tenEigensolve_d(eval[1], evec[1], tfx->gageTen2 + 7); } - + /* set ten2Use */ if (seedProbe) { /* we're on the *very* 1st probe per tract, at the seed pt */ @@ -135,7 +145,11 @@ lastVec = tfx->seedEvec; } else { /* we're past the first step */ - ELL_3V_NORM(tfx->lastDir, tfx->lastDir, len); + /* Arish says: "Bug len has not been initialized and don't think + its needed". The first part is not a problem; "len" is in the + *output* argument of ELL_3V_NORM. The second part seems to be + true, even though Gordon can't currently see why! */ + /* ELL_3V_NORM(tfx->lastDir, tfx->lastDir, len); */ lastVec = tfx->lastDir; } dot[0] = ELL_3V_DOT(lastVec, evec[0]); @@ -158,11 +172,7 @@ /* based on ten2Use, set the rest of the needed fields */ if (tfx->mframeUse) { - double matTmpA[9], matTmpB[9]; - TEN_T2M(matTmpA, tfx->gageTen2 + 7*(tfx->ten2Use)); - ELL_3M_MUL(matTmpB, tfx->mframe, matTmpA); - ELL_3M_MUL(matTmpA, matTmpB, tfx->mframeT); - TEN_M2T(tfx->fiberTen, matTmpA); + TEN_T_COPY(tfx->fiberTen, tens2[tfx->ten2Use]); } else { TEN_T_COPY(tfx->fiberTen, tfx->gageTen2 + 7*(tfx->ten2Use)); } Modified: teem/trunk/src/ten/tenDwiGage.c =================================================================== --- teem/trunk/src/ten/tenDwiGage.c 2008-12-08 17:38:33 UTC (rev 4034) +++ teem/trunk/src/ten/tenDwiGage.c 2008-12-09 11:02:01 UTC (rev 4035) @@ -306,8 +306,8 @@ egrad += 3; for (ii=0; ii<nn; ii++) { double argA, argB, sigA, sigB; - argA = -pvlData->tec2->bValue*TEN_T3V_CONTR(tenA, egrad + 3*ii); - argB = -pvlData->tec2->bValue*TEN_T3V_CONTR(tenB, egrad + 3*ii); + argA = -pvlData->tec1->bValue*TEN_T3V_CONTR(tenA, egrad + 3*ii); + argB = -pvlData->tec1->bValue*TEN_T3V_CONTR(tenB, egrad + 3*ii); if (pvlData->levmarUseFastExp) { sigA = airFastExp(argA); sigB = airFastExp(argB); @@ -315,7 +315,7 @@ sigA = exp(argA); sigB = exp(argB); } - xx[ii] = pvlData->tec2->knownB0*(pp[1]*sigA + (1-pp[1])*sigB); + xx[ii] = pvlData->tec1->knownB0*(pp[1]*sigA + (1-pp[1])*sigB); } return; } @@ -580,106 +580,114 @@ /* Pointer to the location where the two tensor will be written */ twoTen = pvl->directAnswer[tenDwiGage2TensorPeled]; /* Estimate the DWI error, error is given as standard deviation */ - pvlData->tec2->recordErrorDwi = AIR_FALSE; + pvlData->tec1->recordErrorDwi = AIR_FALSE; /* Estimate the single tensor */ - tenEstimate1TensorSingle_d(pvlData->tec2, pvlData->ten1, dwiAll); + tenEstimate1TensorSingle_d(pvlData->tec1, pvlData->ten1, dwiAll); /* Get the eigenValues and eigen vectors for this tensor */ tenEigensolve_d(pvlData->ten1Eval, pvlData->ten1Evec, pvlData->ten1); /* Get westins Cp */ Cp = tenAnisoEval_d(pvlData->ten1Eval, tenAniso_Cp1); - /* Calculate the residual, need the variance to sqr it */ - /* residual = pvlData->tec2->errorDwi*pvlData->tec2->errorDwi; */ - /* Calculate the AIC for single tensor fit */ - /* AICSingFit = _tenComputeAIC(residual, pvlData->tec2->dwiNum, 6); */ + /* Only do two-tensor fitting if CP is greater or equal to than a + user-defined threshold */ + if (Cp >= pvlData->levmarMinCp) { + /* Calculate the residual, need the variance to sqr it */ + /* residual = pvlData->tec1->errorDwi*pvlData->tec1->errorDwi; */ + /* Calculate the AIC for single tensor fit */ + /* AICSingFit = _tenComputeAIC(residual, pvlData->tec1->dwiNum, 6); */ - /* the CP-based test is gone; caller's responsibility */ + /* the CP-based test is gone; caller's responsibility */ - /* rotate DW gradients by inverse of eigenvector column matrix - and place into pvlData->nten1EigenGrads (which has been - allocated by _tenDwiGagePvlDataNew()) */ - grad = AIR_CAST(double *, kindData->ngrad->data); - egrad = AIR_CAST(double *, pvlData->nten1EigenGrads->data); - for (gi=0; gi<kindData->ngrad->axis[1].size; gi++) { - /* yes, this is also transforming some zero-length (B0) gradients; - that's harmless */ - ELL_3MV_MUL(egrad, pvlData->ten1Evec, grad); - grad += 3; - egrad += 3; - } + /* rotate DW gradients by inverse of eigenvector column matrix + and place into pvlData->nten1EigenGrads (which has been + allocated by _tenDwiGagePvlDataNew()) */ + grad = AIR_CAST(double *, kindData->ngrad->data); + egrad = AIR_CAST(double *, pvlData->nten1EigenGrads->data); + for (gi=0; gi<kindData->ngrad->axis[1].size; gi++) { + /* yes, this is also transforming some zero-length (B0) gradients; + that's harmless */ + ELL_3MV_MUL(egrad, pvlData->ten1Evec, grad); + grad += 3; + egrad += 3; + } - /* Lower and upper bounds for the NLLS routine */ - loBnd[0] = 0.0; - loBnd[1] = 0.0; - loBnd[2] = -AIR_PI/2; - loBnd[3] = -AIR_PI/2; - upBnd[0] = pvlData->ten1Eval[0]*5; - upBnd[1] = 1.0; - upBnd[2] = AIR_PI/2; - upBnd[3] = AIR_PI/2; - /* Starting point for the NLLS */ - guess[0] = pvlData->ten1Eval[0]; - guess[1] = 0.5; + /* Lower and upper bounds for the NLLS routine */ + loBnd[0] = 0.0; + loBnd[1] = 0.0; + loBnd[2] = -AIR_PI/2; + loBnd[3] = -AIR_PI/2; + upBnd[0] = pvlData->ten1Eval[0]*5; + upBnd[1] = 1.0; + upBnd[2] = AIR_PI/2; + upBnd[3] = AIR_PI/2; + /* Starting point for the NLLS */ + guess[0] = pvlData->ten1Eval[0]; + guess[1] = 0.5; - guess[2] = AIR_PI/4; - guess[3] = -AIR_PI/4; - /* - guess[2] = AIR_AFFINE(0, airDrandMT_r(pvlData->randState), 1, - AIR_PI/6, AIR_PI/3); - guess[3] = AIR_AFFINE(0, airDrandMT_r(pvlData->randState), 1, - -AIR_PI/6, -AIR_PI/3); - */ - /* Fill in the constraints for the LM optimization, - the threshold of error difference */ - opts[0] = pvlData->levmarTau; - opts[1] = pvlData->levmarEps1; - opts[2] = pvlData->levmarEps2; - opts[3] = pvlData->levmarEps3; - /* Very imp to set this opt, note that only forward - differences are used to approx Jacobian */ - opts[4] = pvlData->levmarDelta; - - /* run NLLS, results are stored back into guess[] */ - pvlData->levmarUseFastExp = AIR_FALSE; - lmret = dlevmar_bc_dif(_tenLevmarPeledCB, guess, pvlData->tec2->dwi, - PARAMS, pvlData->tec2->dwiNum, loBnd, upBnd, - pvlData->levmarMaxIter, opts, - pvlData->levmarInfo, - NULL, NULL, pvlData); - if (-1 == lmret) { - ctx->errNum = 1; - sprintf(ctx->errStr, "%s: dlevmar_bc_dif() failed!", me); - } else { - /* Get the AIC for the two tensor fit, use the levmarinfo - to get the residual */ + guess[2] = AIR_PI/4; + guess[3] = -AIR_PI/4; /* - residual = pvlData->levmarInfo[1]/pvlData->tec2->dwiNum; - AICTwoFit = _tenComputeAIC(residual, pvlData->tec2->dwiNum, 12); + guess[2] = AIR_AFFINE(0, airDrandMT_r(pvlData->randState), 1, + AIR_PI/6, AIR_PI/3); + guess[3] = AIR_AFFINE(0, airDrandMT_r(pvlData->randState), 1, + -AIR_PI/6, -AIR_PI/3); */ - /* Form the tensors using the estimated pp, returned in guess */ - _tenPeledRotate2D(tenA, guess[0], pvlData->ten1Eval[2], guess[2]); - _tenPeledRotate2D(tenB, guess[0], pvlData->ten1Eval[2], guess[3]); - TEN_T2M(matA, tenA); - TEN_T2M(matB, tenB); + /* Fill in the constraints for the LM optimization, + the threshold of error difference */ + opts[0] = pvlData->levmarTau; + opts[1] = pvlData->levmarEps1; + opts[2] = pvlData->levmarEps2; + opts[3] = pvlData->levmarEps3; + /* Very imp to set this opt, note that only forward + differences are used to approx Jacobian */ + opts[4] = pvlData->levmarDelta; + + /* run NLLS, results are stored back into guess[] */ + pvlData->levmarUseFastExp = AIR_FALSE; + lmret = dlevmar_bc_dif(_tenLevmarPeledCB, guess, pvlData->tec1->dwi, + PARAMS, pvlData->tec1->dwiNum, loBnd, upBnd, + pvlData->levmarMaxIter, opts, + pvlData->levmarInfo, + NULL, NULL, pvlData); + if (-1 == lmret) { + ctx->errNum = 1; + sprintf(ctx->errStr, "%s: dlevmar_bc_dif() failed!", me); + } else { + /* Get the AIC for the two tensor fit, use the levmarinfo + to get the residual */ + /* + residual = pvlData->levmarInfo[1]/pvlData->tec1->dwiNum; + AICTwoFit = _tenComputeAIC(residual, pvlData->tec1->dwiNum, 12); + */ + /* Form the tensors using the estimated pp, returned in guess */ + _tenPeledRotate2D(tenA, guess[0], pvlData->ten1Eval[2], guess[2]); + _tenPeledRotate2D(tenB, guess[0], pvlData->ten1Eval[2], guess[3]); + TEN_T2M(matA, tenA); + TEN_T2M(matB, tenB); - ELL_3M_TRANSPOSE(rott, pvlData->ten1Evec); - ELL_3M_MUL(matTmp, matA, pvlData->ten1Evec); - ELL_3M_MUL(matA, rott, matTmp); - ELL_3M_MUL(matTmp, matB, pvlData->ten1Evec); - ELL_3M_MUL(matB, rott, matTmp); - - /* Copy two two tensors */ - /* guess[1] is population fraction of first tensor */ - if (guess[1] > 0.5) { - twoTen[7] = guess[1]; - TEN_M2T(twoTen + 0, matA); - TEN_M2T(twoTen + 7, matB); - } else { - twoTen[7] = 1 - guess[1]; - TEN_M2T(twoTen + 0, matB); - TEN_M2T(twoTen + 7, matA); + ELL_3M_TRANSPOSE(rott, pvlData->ten1Evec); + ELL_3M_MUL(matTmp, matA, pvlData->ten1Evec); + ELL_3M_MUL(matA, rott, matTmp); + ELL_3M_MUL(matTmp, matB, pvlData->ten1Evec); + ELL_3M_MUL(matB, rott, matTmp); + + /* Copy two two tensors */ + /* guess[1] is population fraction of first tensor */ + if (guess[1] > 0.5) { + twoTen[7] = guess[1]; + TEN_M2T(twoTen + 0, matA); + TEN_M2T(twoTen + 7, matB); + } else { + twoTen[7] = 1 - guess[1]; + TEN_M2T(twoTen + 0, matB); + TEN_M2T(twoTen + 7, matA); + } + twoTen[0] = 1; } - twoTen[0] = 1; + } else { + /* its too planar- just do single tensor fit */ + TEN_T_COPY(twoTen, pvlData->ten1); + TEN_T_SET(twoTen + 7, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0); } #undef PARAMS #else @@ -700,7 +708,7 @@ if (info[1] > 0) { /* Returning the standard deviation */ pvl->directAnswer[tenDwiGage2TensorPeledError][0] = - sqrt(info[1]/pvlData->tec2->dwiNum); + sqrt(info[1]/pvlData->tec1->dwiNum); } } if (GAGE_QUERY_ITEM_TEST(pvl->query, tenDwiGage2TensorPeledAndError)) { @@ -833,7 +841,7 @@ pvlData->levmarEps2 = 1E-8; pvlData->levmarEps3 = 1E-8; pvlData->levmarDelta = 1E-8; - pvlData->levmarMinCp = 0.08; + pvlData->levmarMinCp = 0.1; /* pvlData->levmarInfo[] is output; not initialized */ This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |