|
From: <kin...@us...> - 2004-03-07 01:08:19
|
Update of /cvsroot/teem/teem/src/ten In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv18490 Modified Files: chan.c Log Message: tenSimulateOne: - fixed basic bug in handling of off-diagonal entries of B-matrix tenEstimateLinear4D: - fixed basic bug in computation of unknown B0 field - should be able to save out tensor error field with knownB0 either true or false - changed tweaking heuristic in automatic DWI threshold detection Index: chan.c =================================================================== RCS file: /cvsroot/teem/teem/src/ten/chan.c,v retrieving revision 1.30 retrieving revision 1.31 diff -C2 -d -r1.30 -r1.31 *** chan.c 27 Feb 2004 10:23:56 -0000 1.30 --- chan.c 7 Mar 2004 00:52:57 -0000 1.31 *************** *** 141,145 **** ** output: ** ten[0]..ten[6] will be the confidence value followed by the tensor ! ** -------------- !IF knownB0 ------------------------- ** input: ** dwi[0]..dwi[DD-1] are the DD DWI values, --- 141,145 ---- ** output: ** ten[0]..ten[6] will be the confidence value followed by the tensor ! ** -------------- IF !knownB0 ------------------------- ** input: ** dwi[0]..dwi[DD-1] are the DD DWI values, *************** *** 160,163 **** --- 160,167 ---- if (knownB0) { + if (B0P) { + /* we save this as a courtesy */ + *B0P = AIR_MAX(dwi[0], 1); + } logB0 = log(AIR_MAX(dwi[0], 1)); mean = 0; *************** *** 197,202 **** ten[j+1] = tmp; } else { if (B0P) { ! *B0P = tmp; } } --- 201,207 ---- ten[j+1] = tmp; } else { + /* we're on seventh row, for finding B0 */ if (B0P) { ! *B0P = exp(b*tmp); } } *************** *** 266,296 **** float thresh, float soft, float b) { char me[]="tenEstimateLinear4D", err[AIR_STRLEN_MED]; ! Nrrd *nemat, *ncrop, *nhist; airArray *mop; ! int E, DD, d, II, sx, sy, sz, cmin[4], cmax[4], amap[4] = {-1,1,2,3}; float *ten, dwi1[_TEN_MAX_DWI_NUM], *terr=NULL, ! *B0=NULL, (*lup)(const void *, size_t); ! double *emat; NrrdRange *range; - - /* - HEY: uncomment with tenSimulate stuff is working for unknown B0 case - Nrrd *nbmat; - double *bmat; - const char *bk; float dwi2[_TEN_MAX_DWI_NUM], te, d1, d2; - */ if (!(nten && ndwi && _nbmat)) { ! /* nerrP can be NULL */ sprintf(err, "%s: got NULL pointer", me); biffAdd(TEN, err); return 1; } - if (!knownB0 && !nB0P) { - /* For the time being, this needless restriction helps simplify things */ - sprintf(err, "%s: sorry, don't know (and hence will estimate) B=0 image, " - "but didn't get a Nrrd* to put it in", me); - biffAdd(TEN, err); return 1; - } if (!( 4 == ndwi->dim && 7 <= ndwi->axis[0].size )) { sprintf(err, "%s: dwi should be 4-D array with axis 0 size >= 7", me); --- 271,288 ---- float thresh, float soft, float b) { char me[]="tenEstimateLinear4D", err[AIR_STRLEN_MED]; ! Nrrd *nemat, *nbmat, *ncrop, *nhist; airArray *mop; ! int E, DD, d, II, sx, sy, sz, cmin[4], cmax[4], amap[4]; float *ten, dwi1[_TEN_MAX_DWI_NUM], *terr=NULL, ! _B0, *B0, (*lup)(const void *, size_t); ! double *emat, *bmat; NrrdRange *range; float dwi2[_TEN_MAX_DWI_NUM], te, d1, d2; if (!(nten && ndwi && _nbmat)) { ! /* nerrP and _NB0P can be NULL */ sprintf(err, "%s: got NULL pointer", me); biffAdd(TEN, err); return 1; } if (!( 4 == ndwi->dim && 7 <= ndwi->axis[0].size )) { sprintf(err, "%s: dwi should be 4-D array with axis 0 size >= 7", me); *************** *** 318,326 **** mop = airMopNew(); ! /* HEY, right now, nbmat is not used for anything (not tenSimulate stuff) ! because tenSimulate currently can't handle the unknown B0 case */ ! /* airMopAdd(mop, nbmat=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); */ airMopAdd(mop, nemat=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); ! if (tenEMatrixCalc(nemat, _nbmat, knownB0)) { sprintf(err, "%s: trouble computing estimation matrix", me); biffAdd(TEN, err); return 1; --- 310,320 ---- mop = airMopNew(); ! airMopAdd(mop, nbmat=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); ! if (nrrdConvert(nbmat, _nbmat, nrrdTypeDouble)) { ! sprintf(err, "%s: trouble converting to doubles", me); ! biffMove(TEN, err, NRRD); return 1; ! } airMopAdd(mop, nemat=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); ! if (tenEMatrixCalc(nemat, nbmat, knownB0)) { sprintf(err, "%s: trouble computing estimation matrix", me); biffAdd(TEN, err); return 1; *************** *** 347,351 **** biffMove(TEN, err, NRRD); airMopError(mop); return 1; } ! if (_tenFindValley(&thresh, nhist, 0.85)) { sprintf(err, "%s: problem finding DW histogram valley", me); biffAdd(TEN, err); airMopError(mop); return 1; --- 341,345 ---- biffMove(TEN, err, NRRD); airMopError(mop); return 1; } ! if (_tenFindValley(&thresh, nhist, 0.75, AIR_FALSE)) { sprintf(err, "%s: problem finding DW histogram valley", me); biffAdd(TEN, err); airMopError(mop); return 1; *************** *** 358,367 **** } if (nterrP) { - - /* HEY HEY HEY HEY */ - sprintf(err, "%s: sorry, this functionality currently disabled", me); - biffAdd(TEN, err); return 1; - /* HEY HEY HEY HEY */ - if (!(*nterrP)) { *nterrP = nrrdNew(); --- 352,355 ---- *************** *** 386,391 **** airMopAdd(mop, *nB0P, (airMopper)nrrdNuke, airMopOnError); B0 = (float*)((*nB0P)->data); } ! /* bmat = (double*)(nbmat->data); */ emat = (double*)(nemat->data); ten = (float*)(nten->data); --- 374,381 ---- airMopAdd(mop, *nB0P, (airMopper)nrrdNuke, airMopOnError); B0 = (float*)((*nB0P)->data); + } else { + B0 = NULL; } ! bmat = (double*)(nbmat->data); emat = (double*)(nemat->data); ten = (float*)(nten->data); *************** *** 398,421 **** fprintf(stderr, "%s: input dwi1[%d] = %g\n", me, d, dwi1[d]); } ! tenEstimateLinearSingle(ten, B0, dwi1, emat, DD, knownB0, thresh, soft, b); ! if (tenVerbose) fprintf(stderr, "%s: output ten = (%g) %g,%g,%g %g,%g %g\n", me, ten[0], ten[1], ten[2], ten[3], ten[4], ten[5], ten[6]); if (nterrP) { - /* - tenSimulateOne(dwi2, dwi1[0], ten, bmat, DD, b); te = 0; ! for (d=0; d<DD; d++) { ! d1 = AIR_MAX(dwi1[d], 1); ! d2 = AIR_MAX(dwi2[d], 1); ! te += (d1 - d2)*(d1 - d2); ! if (tenVerbose) ! fprintf(stderr, "%s: dwi1,2[%d] = %g,%g --> %g\n", ! me, d, dwi1[d], dwi2[d], te); } *terr = sqrt(te); - */ terr += 1; ! } ten += 7; if (B0) { --- 388,425 ---- fprintf(stderr, "%s: input dwi1[%d] = %g\n", me, d, dwi1[d]); } ! tenEstimateLinearSingle(ten, &_B0, ! dwi1, emat, DD, knownB0, thresh, soft, b); ! if (nB0P) { ! *B0 = _B0; ! } ! if (tenVerbose) { fprintf(stderr, "%s: output ten = (%g) %g,%g,%g %g,%g %g\n", me, ten[0], ten[1], ten[2], ten[3], ten[4], ten[5], ten[6]); + } if (nterrP) { te = 0; ! if (knownB0) { ! tenSimulateOne(dwi2, _B0, ten, bmat, DD, b); ! for (d=1; d<DD; d++) { ! d1 = AIR_MAX(dwi1[d], 1); ! d2 = AIR_MAX(dwi2[d], 1); ! te += (d1 - d2)*(d1 - d2); ! } ! te /= (DD-1); ! } else { ! tenSimulateOne(dwi2, _B0, ten, bmat, DD+1, b); ! for (d=0; d<DD; d++) { ! d1 = AIR_MAX(dwi1[d], 1); ! /* tenSimulateOne always puts the B0 in the beginning of ! the dwi vector, but in this case we didn't have it in ! the input dwi vecctor */ ! d2 = AIR_MAX(dwi2[d+1], 1); ! te += (d1 - d2)*(d1 - d2); ! } ! te /= DD; } *terr = sqrt(te); terr += 1; ! } ten += 7; if (B0) { *************** *** 425,428 **** --- 429,433 ---- /* tenEigenvalueClamp(nten, nten, 0, AIR_NAN); */ + ELL_4V_SET(amap, -1, 1, 2, 3); nrrdAxisInfoCopy(nten, ndwi, amap, NRRD_AXIS_INFO_NONE); *************** *** 434,444 **** ******** tenSimulateOne ** ! ** simulate a diffusion weighted measurement ** */ void ! tenSimulateOne(float *dwi, float B0, float *ten, ! double *bmat, int DD, float b) { ! double v[_TEN_MAX_DWI_NUM]; int i, j; --- 439,456 ---- ******** tenSimulateOne ** ! ** given a tensor, simulate the set of diffusion weighted measurements ! ** represented by the given B matrix ** + ** NOTE: the mindset if this function is very much "knownB0==true": + ** B0 is required as an argument (and its always copied to dwi[0]), + ** and the given bmat is assumed to have DD-1 rows, + ** so dwi[0] through dwi[DD-1] (DD values total) are set in the output. */ void ! tenSimulateOne(float *dwi, ! float B0, float *ten, double *bmat, int DD, float b) { ! double vv; ! /* this is how we multiply the off-diagonal entries by 2 */ ! double matwght[6] = {1, 2, 2, 1, 2, 1}; int i, j; *************** *** 449,459 **** } for (i=0; i<DD-1; i++) { ! v[i] = 0; for (j=0; j<6; j++) { ! v[i] += bmat[j + 6*i]*ten[j+1]; } ! dwi[i+1] = AIR_MAX(B0, 1)*exp(-b*v[i]); if (tenVerbose) { ! fprintf(stderr, "v[%d] = %g --> dwi = %g\n", i, v[i], dwi[i+1]); } } --- 461,471 ---- } for (i=0; i<DD-1; i++) { ! vv = 0; for (j=0; j<6; j++) { ! vv += matwght[j]*bmat[j + 6*i]*ten[j+1]; } ! dwi[i+1] = AIR_MAX(B0, 1)*exp(-b*vv); if (tenVerbose) { ! fprintf(stderr, "v[%d] = %g --> dwi = %g\n", i, vv, dwi[i+1]); } } |