[Teem-commit] teem/src/ten chan.c,1.30,1.31

[<kin...@us...>]

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]);
      }
    }