[Teem-commit] teem/src/ten epireg.c,1.26,1.27 tendEpireg.c,1.14,1.15

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

Update of /cvsroot/teem/teem/src/ten
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv29863

Modified Files:
	epireg.c tendEpireg.c 
Log Message:
NO API CHANGES: attempts at refining registration method, and other inconsequential changes

Index: epireg.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/epireg.c,v
retrieving revision 1.26
retrieving revision 1.27
diff -C2 -d -r1.26 -r1.27
*** epireg.c	7 Mar 2004 00:49:06 -0000	1.26
--- epireg.c	8 Mar 2004 14:41:24 -0000	1.27
***************
*** 333,336 ****
--- 333,337 ----
         and merge up (to bright) all small dark pieces, where
         (currently) small is big/2 */
+     big = -1;
      if (nrrdCCFind(ncc, &nval, nthr[ni], nrrdTypeDefault, conny)
  	|| nrrdCCSize(nsize, ncc)
***************
*** 338,345 ****
  	|| nrrdCCMerge(ncc, ncc, nval, -1, big-1, 0, conny)
  	|| nrrdCCRevalue(nthr[ni], ncc, nval)) {
!       sprintf(err, "%s: trouble with 3-D processing nthr[%d]", me, ni);
!       biffMove(TEN, err, NRRD); return 1;
      }
      for (z=0; z<sz; z++) {
        if ( nrrdSlice(nslc, nthr[ni], 2, z)
  	   || nrrdCCFind(ncc, &nval, nslc, nrrdTypeDefault, conny)
--- 339,353 ----
  	|| nrrdCCMerge(ncc, ncc, nval, -1, big-1, 0, conny)
  	|| nrrdCCRevalue(nthr[ni], ncc, nval)) {
!       if (big) {
! 	sprintf(err, "%s: trouble with 3-D processing nthr[%d]", me, ni);
! 	biffMove(TEN, err, NRRD); return 1;
!       } else {
! 	sprintf(err, "%s: got size 0 for biggest bright CC of nthr[%d]",
! 		me, ni);
! 	biffAdd(TEN, err); return 1;
!       }
      }
      for (z=0; z<sz; z++) {
+       big = -1;
        if ( nrrdSlice(nslc, nthr[ni], 2, z)
  	   || nrrdCCFind(ncc, &nval, nslc, nrrdTypeDefault, conny)
***************
*** 349,354 ****
  	   || nrrdCCRevalue(nslc, ncc, nval)
  	   || nrrdSplice(nthr[ni], nthr[ni], nslc, 2, z) ) {
! 	sprintf(err, "%s: trouble processing slice %d of nthr[%d]", me, z, ni);
! 	biffMove(TEN, err, NRRD); return 1;
        }
      }
--- 357,368 ----
  	   || nrrdCCRevalue(nslc, ncc, nval)
  	   || nrrdSplice(nthr[ni], nthr[ni], nslc, 2, z) ) {
! 	if (big) {
! 	  sprintf(err, "%s: trouble processing slice %d of nthr[%d]",
! 		  me, z, ni);
! 	  biffMove(TEN, err, NRRD); return 1;
! 	} else {
! 	  /* biggest bright CC on this slice had size 0
! 	     <==> there was no bright CC on this slice, move on */
! 	}
        }
      }
***************
*** 414,422 ****
  	}
        }
-       if (!N) {
- 	sprintf(err, "%s: saw no non-zero pixels in nthresh[%d]; "
- 		"DWI threshold too high?", me, ni);
- 	biffAdd(TEN, err); return 1;
-       }
        if (N == sx*sy) {
  	sprintf(err, "%s: saw only non-zero pixels in nthresh[%d]; "
--- 428,431 ----
***************
*** 424,452 ****
  	biffAdd(TEN, err); return 1;
        }
!       mx /= N;
!       my /= N;
!       cx = sx/2.0;
!       cy = sy/2.0;
!       /* ------ find M02, M11, M20 */
!       M02 = M11 = M20 = 0.0;
!       for (yi=0; yi<sy; yi++) {
! 	for (xi=0; xi<sx; xi++) {
! 	  val = thr[xi + sx*yi];
! 	  x = xi - cx;
! 	  y = yi - cy;
! 	  M02 += y*y*val;
! 	  M11 += x*y*val;
! 	  M20 += x*x*val;
  	}
        }
-       M02 /= N;
-       M11 /= N;
-       M20 /= N;
-       /* ------ set output */
-       mom[MEAN_X] = mx;
-       mom[MEAN_Y] = my;
-       mom[M_02] = M02;
-       mom[M_11] = M11;
-       mom[M_20] = M20;
        thr += sx*sy;
        mom += 5;
--- 433,471 ----
  	biffAdd(TEN, err); return 1;
        }
!       if (N) {
! 	/* there were non-zero pixels */
! 	mx /= N;
! 	my /= N;
! 	cx = sx/2.0;
! 	cy = sy/2.0;
! 	/* ------ find M02, M11, M20 */
! 	M02 = M11 = M20 = 0.0;
! 	for (yi=0; yi<sy; yi++) {
! 	  for (xi=0; xi<sx; xi++) {
! 	    val = thr[xi + sx*yi];
! 	    x = xi - cx;
! 	    y = yi - cy;
! 	    M02 += y*y*val;
! 	    M11 += x*y*val;
! 	    M20 += x*x*val;
! 	  }
  	}
+ 	M02 /= N;
+ 	M11 /= N;
+ 	M20 /= N;
+ 	/* ------ set output */
+ 	mom[MEAN_X] = mx;
+ 	mom[MEAN_Y] = my;
+ 	mom[M_02] = M02;
+ 	mom[M_11] = M11;
+ 	mom[M_20] = M20;
+       } else {
+ 	/* there were no non-zero pixels */
+ 	mom[MEAN_X] = 0;
+ 	mom[MEAN_Y] = 0;
+ 	mom[M_02] = 0;
+ 	mom[M_11] = 0;
+ 	mom[M_20] = 0;
        }
        thr += sx*sy;
        mom += 5;
***************
*** 509,526 ****
  _tenEpiRegEstimHST(Nrrd *nhst, Nrrd *npxfr, int ninLen, Nrrd *ngrad) {
    char me[]="_tenEpiRegEstimHST", err[AIR_STRLEN_MED];
!   double *hst, *grad, *mat1, *vec, *ans, *pxfr, *gA, *gB;
    int z, sz, A, B, npairs, ri;
!   Nrrd *nmat1, *nvec, *ninv, *nans;
    airArray *mop;
  
    mop = airMopNew();
!   airMopAdd(mop, nmat1=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
!   airMopAdd(mop, ninv=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
    airMopAdd(mop, nvec=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
    airMopAdd(mop, nans=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
!   npairs = ninLen*(ninLen-1);
!   sz = npxfr->axis[1].size;
!   if (nrrdMaybeAlloc(nhst, nrrdTypeDouble, 2, 9, sz)
!       || nrrdMaybeAlloc(nmat1, nrrdTypeDouble, 2, 3, npairs)
        || nrrdMaybeAlloc(nvec, nrrdTypeDouble, 2, 1, npairs)) {
      sprintf(err, "%s: couldn't allocate HST nrrd", me);
--- 528,560 ----
  _tenEpiRegEstimHST(Nrrd *nhst, Nrrd *npxfr, int ninLen, Nrrd *ngrad) {
    char me[]="_tenEpiRegEstimHST", err[AIR_STRLEN_MED];
!   double *hst, *grad, *mat, *vec, *ans, *pxfr, *gA, *gB;
    int z, sz, A, B, npairs, ri;
!   Nrrd **nmat, *nvec, **ninv, *nans;
    airArray *mop;
+   int order;
+ 
+   order = 1;
+ 
+   sz = npxfr->axis[1].size;
+   npairs = ninLen*(ninLen-1);
  
    mop = airMopNew();
!   nmat = (Nrrd**)calloc(sz, sizeof(Nrrd*));
!   ninv = (Nrrd**)calloc(sz, sizeof(Nrrd*));
!   airMopAdd(mop, nmat, airFree, airMopAlways);
!   airMopAdd(mop, ninv, airFree, airMopAlways);
!   for (z=0; z<sz; z++) {
!     airMopAdd(mop, nmat[z]=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
!     if (nrrdMaybeAlloc(nmat[z], nrrdTypeDouble, 2,
! 		       (1 == order ? 3 : 9), npairs)) {
!       sprintf(err, "%s: couldn't allocate fitting matrices", me);
!       biffMove(TEN, err, NRRD); airMopError(mop); return 1;
!     }
!     airMopAdd(mop, ninv[z]=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
!   }
    airMopAdd(mop, nvec=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
    airMopAdd(mop, nans=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
!   if (nrrdMaybeAlloc(nhst, nrrdTypeDouble, 2,
! 		     (1 == order ? 9 : 27), sz)
        || nrrdMaybeAlloc(nvec, nrrdTypeDouble, 2, 1, npairs)) {
      sprintf(err, "%s: couldn't allocate HST nrrd", me);
***************
*** 528,539 ****
    }
    nrrdAxisInfoSet(nhst, nrrdAxisInfoLabel,
! 		  "Hx,Hy,Hz,Sx,Sy,Sz,Tx,Ty,Tz", "z");
!   grad = (double *)(ngrad->data);
!   mat1 = (double *)(nmat1->data);
!   vec = (double *)(nvec->data);
  
!   /* ------ find Sx, Sy, Sz per slice */
    for (z=0; z<sz; z++) {
!     hst = (double *)(nhst->data) + 0 + 9*z;
      ri = 0;
      for (A=0; A<ninLen; A++) {
--- 562,574 ----
    }
    nrrdAxisInfoSet(nhst, nrrdAxisInfoLabel,
! 		  (1 == order 
! 		   ? "Hx,Hy,Hz,Sx,Sy,Sz,Tx,Ty,Tz"
! 		   : "HST parms"), "z");
  
!   /* ------ per-slice: compute model fitting matrix and its pseudo-inverse */
!   grad = (double *)(ngrad->data);
    for (z=0; z<sz; z++) {
!     hst = (double *)(nhst->data) + (1 == order ? 9 : 27)*z;
!     mat = (double *)(nmat[z]->data);
      ri = 0;
      for (A=0; A<ninLen; A++) {
***************
*** 543,565 ****
  	gA = grad + 0 + 3*A;
  	gB = grad + 0 + 3*B;
! 	ELL_3V_SET(mat1 + 3*ri,
! 		   pxfr[SCALE]*gA[0] - gB[0],
! 		   pxfr[SCALE]*gA[1] - gB[1],
! 		   pxfr[SCALE]*gA[2] - gB[2]);
! 	vec[ri] = 1 - pxfr[SCALE];
  	ri += 1;
        }
      }
!     ell_Nm_pseudo_inv(ninv, nmat1);
!     ell_Nm_mul(nans, ninv, nvec);
!     ans = (double *)(nans->data);
!     hst[3] = ans[0];
!     hst[4] = ans[1];
!     hst[5] = ans[2];
    }
  
    /* ------ find Hx, Hy, Hz per slice */
    for (z=0; z<sz; z++) {
!     hst = (double *)(nhst->data) + 0 + 9*z;
      ri = 0;
      for (A=0; A<ninLen; A++) {
--- 578,636 ----
  	gA = grad + 0 + 3*A;
  	gB = grad + 0 + 3*B;
! 	if (1 == order) {
! 	  ELL_3V_SET(mat + 3*ri,
! 		     gB[0] - pxfr[SCALE]*gA[0],
! 		     gB[1] - pxfr[SCALE]*gA[1],
! 		     gB[2] - pxfr[SCALE]*gA[2]);
! 	} else {
! 	  ELL_9V_SET(mat + 9*ri,
! 		     gB[0] - pxfr[SCALE]*gA[0],
! 		     gB[1] - pxfr[SCALE]*gA[1],
! 		     gB[2] - pxfr[SCALE]*gA[2],
! 		     gB[0] - pxfr[SCALE]*gA[0],
! 		     gB[1] - pxfr[SCALE]*gA[1],
! 		     gB[2] - pxfr[SCALE]*gA[2],
! 		     gB[0] - pxfr[SCALE]*gA[0],
! 		     gB[1] - pxfr[SCALE]*gA[1],
! 		     gB[2] - pxfr[SCALE]*gA[2]);
! 		     /*
! 		     gB[0]*gB[0] - pxfr[SCALE]*gA[0]*gA[0],
! 		     gB[1]*gB[1] - pxfr[SCALE]*gA[1]*gA[1],
! 		     gB[2]*gB[2] - pxfr[SCALE]*gA[2]*gA[2],
! 		     gB[0]*gB[1] - pxfr[SCALE]*gA[0]*gA[1],
! 		     gB[0]*gB[2] - pxfr[SCALE]*gA[0]*gA[2],
! 		     gB[1]*gB[2] - pxfr[SCALE]*gA[1]*gA[2]);
! 		     */
! 	}
  	ri += 1;
        }
      }
!     if (nrrdHasNonExist(nmat[z])) {
!       /* as happens if there were zero slices in the segmentation output */
!       if (1 == order) {
! 	ELL_3V_SET(hst + 0*3, 0, 0, 0);
! 	ELL_3V_SET(hst + 1*3, 0, 0, 0);
! 	ELL_3V_SET(hst + 2*3, 0, 0, 0);
!       } else {
! 	ELL_9V_SET(hst + 0*9, 0, 0, 0, 0, 0, 0, 0, 0, 0);
! 	ELL_9V_SET(hst + 1*9, 0, 0, 0, 0, 0, 0, 0, 0, 0);
! 	ELL_9V_SET(hst + 2*9, 0, 0, 0, 0, 0, 0, 0, 0, 0);
!       }
!     } else {
!       if (ell_Nm_pseudo_inv(ninv[z], nmat[z])) {
! 	sprintf(err, "%s: trouble estimating model (slice %d)", me, z);
! 	biffMove(TEN, err, ELL); airMopError(mop); return 1;
!       }
!     }
    }
  
    /* ------ find Hx, Hy, Hz per slice */
+   vec = (double *)(nvec->data);
    for (z=0; z<sz; z++) {
!     if (nrrdHasNonExist(nmat[z])) {
!       /* we've already zero-ed out this row in the HST nrrd */
!       continue;
!     }
!     hst = (double *)(nhst->data) + (1 == order ? 9 : 27)*z;
      ri = 0;
      for (A=0; A<ninLen; A++) {
***************
*** 567,591 ****
  	if (A == B) continue;
  	pxfr = (double *)(npxfr->data) + 0 + 5*(z + sz*(A + ninLen*B));
- 	gA = grad + 0 + 3*A;
- 	gB = grad + 0 + 3*B;
- 	ELL_3V_SET(mat1 + 3*ri,
- 		   gB[0] - pxfr[SCALE]*gA[0],
- 		   gB[1] - pxfr[SCALE]*gA[1],
- 		   gB[2] - pxfr[SCALE]*gA[2]);
  	vec[ri] = pxfr[SHEAR];
  	ri += 1;
        }
      }
!     ell_Nm_pseudo_inv(ninv, nmat1);
!     ell_Nm_mul(nans, ninv, nvec);
      ans = (double *)(nans->data);
!     hst[0] = ans[0];
!     hst[1] = ans[1];
!     hst[2] = ans[2];
    }
  
    /* ------ find Tx, Ty, Tz per slice */
    for (z=0; z<sz; z++) {
!     hst = (double *)(nhst->data) + 0 + 9*z;
      ri = 0;
      for (A=0; A<ninLen; A++) {
--- 638,694 ----
  	if (A == B) continue;
  	pxfr = (double *)(npxfr->data) + 0 + 5*(z + sz*(A + ninLen*B));
  	vec[ri] = pxfr[SHEAR];
  	ri += 1;
        }
      }
!     if (ell_Nm_mul(nans, ninv[z], nvec)) {
!       sprintf(err, "%s: trouble estimating model (slice %d): Hx, Hy, Hz",
! 	      me, z);
!       biffMove(TEN, err, ELL); airMopError(mop); return 1;
!     }
      ans = (double *)(nans->data);
!     if (1 == order) {
!       ELL_3V_COPY(hst + 0*3, ans);
!     } else {
!       ELL_9V_COPY(hst + 0*9, ans);
!     }
!   }
! 
!   /* ------ find Sx, Sy, Sz per slice */
!   for (z=0; z<sz; z++) {
!     if (nrrdHasNonExist(nmat[z])) {
!       /* we've already zero-ed out this row in the HST nrrd */
!       continue;
!     }
!     hst = (double *)(nhst->data) + (1 == order ? 9 : 27)*z;
!     ri = 0;
!     for (A=0; A<ninLen; A++) {
!       for (B=0; B<ninLen; B++) {
! 	if (A == B) continue;
! 	pxfr = (double *)(npxfr->data) + 0 + 5*(z + sz*(A + ninLen*B));
! 	vec[ri] = pxfr[SCALE] - 1;
! 	ri += 1;
!       }
!     }
!     if (ell_Nm_mul(nans, ninv[z], nvec)) {
!       sprintf(err, "%s: trouble estimating model (slice %d): Sx, Sy, Sz",
! 	      me, z);
!       biffMove(TEN, err, ELL); airMopError(mop); return 1;
!     }
!     ans = (double *)(nans->data);
!     if (1 == order) {
!       ELL_3V_COPY(hst + 1*3, ans);
!     } else {
!       ELL_9V_COPY(hst + 1*9, ans);
!     }
    }
  
    /* ------ find Tx, Ty, Tz per slice */
    for (z=0; z<sz; z++) {
!     if (nrrdHasNonExist(nmat[z])) {
!       /* we've already zero-ed out this row in the HST nrrd */
!       continue;
!     }
!     hst = (double *)(nhst->data) + (1 == order ? 9 : 27)*z;
      ri = 0;
      for (A=0; A<ninLen; A++) {
***************
*** 593,612 ****
  	if (A == B) continue;
  	pxfr = (double *)(npxfr->data) + 0 + 5*(z + sz*(A + ninLen*B));
- 	gA = grad + 0 + 3*A;
- 	gB = grad + 0 + 3*B;
- 	ELL_3V_SET(mat1 + 3*ri,
- 		   gB[0] - pxfr[SCALE]*gA[0],
- 		   gB[1] - pxfr[SCALE]*gA[1],
- 		   gB[2] - pxfr[SCALE]*gA[2]);
  	vec[ri] = pxfr[TRAN];
  	ri += 1;
        }
      }
!     ell_Nm_pseudo_inv(ninv, nmat1);
!     ell_Nm_mul(nans, ninv, nvec);
      ans = (double *)(nans->data);
!     hst[6] = ans[0];
!     hst[7] = ans[1];
!     hst[8] = ans[2];
    }
  
--- 696,714 ----
  	if (A == B) continue;
  	pxfr = (double *)(npxfr->data) + 0 + 5*(z + sz*(A + ninLen*B));
  	vec[ri] = pxfr[TRAN];
  	ri += 1;
        }
      }
!     if (ell_Nm_mul(nans, ninv[z], nvec)) {
!       sprintf(err, "%s: trouble estimating model (slice %d): Tx, Ty, Tz",
! 	      me, z);
!       biffMove(TEN, err, ELL); airMopError(mop); return 1;
!     }
      ans = (double *)(nans->data);
!     if (1 == order) {
!       ELL_3V_COPY(hst + 2*3, ans);
!     } else {
!       ELL_9V_COPY(hst + 2*9, ans);
!     }
    }
  
***************
*** 670,674 ****
    /* initial ordering: messiness, slice index */
    for (zi=0; zi<sz; zi++) {
!     two[0 + 2*zi] = mess[zi];
      two[1 + 2*zi] = zi;
    }
--- 772,778 ----
    /* initial ordering: messiness, slice index */
    for (zi=0; zi<sz; zi++) {
!     two[0 + 2*zi] = (AIR_EXISTS(mess[zi])
! 		     ? mess[zi]
! 		     : 666);  /* don't use empty slices */
      two[1 + 2*zi] = zi;
    }
***************
*** 693,697 ****
    }
  
!   /* perform fitting for each column in hst */
    hst = (double*)(nhst->data);
    sh = nhst->axis[0].size;
--- 797,802 ----
    }
  
!   /* perform fitting for each column in hst (regardless of 
!      whether we're using a 1st or 2nd order model */
    hst = (double*)(nhst->data);
    sh = nhst->axis[0].size;
***************
*** 724,730 ****
  		 int reference, int ni, int zi,
  		 Nrrd *npxfr, Nrrd *nhst, Nrrd *ngrad) {
!   double *xfr, *hst, *grad;
    int sz, ninLen;
  
    /* these could also have been passed to us, but we can also discover them */
    sz = npxfr->axis[1].size;
--- 829,839 ----
  		 int reference, int ni, int zi,
  		 Nrrd *npxfr, Nrrd *nhst, Nrrd *ngrad) {
!   double *xfr, *hst, *_grad, grad[9]; /* big enough for 2nd order */
    int sz, ninLen;
  
+   int order;
+   
+   order = 1;
+ 
    /* these could also have been passed to us, but we can also discover them */
    sz = npxfr->axis[1].size;
***************
*** 734,742 ****
      /* we use the estimated H,S,T vectors to determine distortion
         as a function of gradient direction, and then invert this */
!     hst = (double*)(nhst->data) + 0 + 9*zi;
!     grad = (double*)(ngrad->data) + 0 + 3*ni;
!     *hhP = ELL_3V_DOT(grad, hst + 0*3);
!     *ssP = 1 + ELL_3V_DOT(grad, hst + 1*3);
!     *ttP = ELL_3V_DOT(grad, hst + 2*3);
    } else {
      /* we register against a specific DWI */
--- 843,865 ----
      /* we use the estimated H,S,T vectors to determine distortion
         as a function of gradient direction, and then invert this */
!     _grad = (double*)(ngrad->data) + 0 + 3*ni;
!     if (1 == order) {
!       hst = (double*)(nhst->data) + 0 + 9*zi;
!       *hhP = ELL_3V_DOT(_grad, hst + 0*3);
!       *ssP = 1 + ELL_3V_DOT(_grad, hst + 1*3);
!       *ttP = ELL_3V_DOT(_grad, hst + 2*3);
!     } else {
!       hst = (double*)(nhst->data) + 0 + 27*zi;
!       ELL_9V_SET(grad, _grad[0], _grad[1], _grad[2],
! 		 _grad[0], _grad[1], _grad[2],
! 		 _grad[0], _grad[1], _grad[2]);
! 		 /*
! 		 _grad[0]*_grad[0], _grad[1]*_grad[1], _grad[2]*_grad[2],
! 		 _grad[0]*_grad[1], _grad[0]*_grad[2], _grad[1]*_grad[2]);
! 		 */
!       *hhP = ELL_9V_DOT(grad, hst + 0*9);
!       *ssP = 1 + ELL_9V_DOT(grad, hst + 1*9);
!       *ttP = ELL_9V_DOT(grad, hst + 2*9);
!     }
    } else {
      /* we register against a specific DWI */

Index: tendEpireg.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/tendEpireg.c,v
retrieving revision 1.14
retrieving revision 1.15
diff -C2 -d -r1.14 -r1.15
*** tendEpireg.c	7 Jan 2004 15:34:31 -0000	1.14
--- tendEpireg.c	8 Mar 2004 14:41:24 -0000	1.15
***************
*** 140,144 ****
    if (rret) {
      airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways);
!     fprintf(stderr, "%s: trouble doing epireg \"%s\":\n%s\n", me, buff, err);
      airMopError(mop); return 1;
    }
--- 140,144 ----
    if (rret) {
      airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways);
!     fprintf(stderr, "%s: trouble doing epireg:\n%s\n", me, err);
      airMopError(mop); return 1;
    }