[Teem-commit] teem/src/ten aniso.c,1.32,1.33 chan.c,1.31,1.32 epireg.c,1.28,1.29 glyph.c,1.36,1.37 miscTen.c,1.24,1.25 mod.c,1.10,1.11 ten.h,1.84,1.85 tenGage.c,1.19,1.20 tendEval.c,1.15,1.16 tendEvec.c,1.16,1.17 tendNorm.c,1.4,1.5 tendPoint.c,1.10,1.11 tensor.c,1.30,1.31

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

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

Modified Files:
	aniso.c chan.c epireg.c glyph.c miscTen.c mod.c ten.h 
	tenGage.c tendEval.c tendEvec.c tendNorm.c tendPoint.c 
	tensor.c 
Log Message:
API RENAME: some function renamings, to make way for future move to doubles,
or at least way from float-specificity:
tenEigensolve  --> tenEigensolve_f
tenMakeOne --> tenMakeOne_f
tenSimulateOne --> tenSimulateOne_f
tenAnisoCalc --> tenAnisoCalc_f
tenEvqOne --> tenEvqOne_f
tenEstimateLinearSingle --> tenEstimateLinearSingle_f
API NEW:
tenEigensolve_d
API CHANGE:
Functions xxx which haven't been renamed to xxx_f, and which took
float arguments (not a pointer to float), now take a double instead.
This is not a change that anyone will notice.  Affected functions are:
tenEstimateLinear3D()
tenEstimateLinear4D()
tenExpand()
tenSimulate()
tenAnisoVolume()
tenSizeScale()
tenAnisoScale()
tenEigenvaluePower()
tenEigenvalueClamp()
tenEigenvalueAdd()

API CHANGE:
More of a big deal: this guy took a float array ("weight"), and
  now its a double array:
tenSizeNormalize(Nrrd *nout, Nrrd *nin, double weight[3],
                 double amount, double target);



Index: aniso.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/aniso.c,v
retrieving revision 1.32
retrieving revision 1.33
diff -C2 -d -r1.32 -r1.33
*** aniso.c	7 May 2004 19:09:58 -0000	1.32
--- aniso.c	11 May 2004 02:39:54 -0000	1.33
***************
*** 21,26 ****
  #include "privateTen.h"
  
- float tenAnisoSigma = 0.000001;
- 
  /*
  ** learned: don't take sqrt(FLT_EPSILON) and expect it to still be
--- 21,24 ----
***************
*** 29,33 ****
  
  /*
! ******** tenAnisoCalc
  **
  ** given an array of three SORTED (descending) eigenvalues "e",
--- 27,31 ----
  
  /*
! ******** tenAnisoCalc_f
  **
  ** given an array of three SORTED (descending) eigenvalues "e",
***************
*** 38,42 ****
  */
  void
! tenAnisoCalc(float c[TEN_ANISO_MAX+1], float e[3]) {
    float e0, e1, e2, stdv, mean, sum, cl, cp, ca, ra, fa, vf, denom;
  
--- 36,40 ----
  */
  void
! tenAnisoCalc_f(float c[TEN_ANISO_MAX+1], float e[3]) {
    float e0, e1, e2, stdv, mean, sum, cl, cp, ca, ra, fa, vf, denom;
  
***************
*** 44,53 ****
  
    if (!( e[0] >= e[1] && e[1] >= e[2] )) {
!     fprintf(stderr, "tenAnisoCalc: eigen values not sorted: "
  	    "%g %g %g (%d %d)\n",
  	    e[0], e[1], e[2], e[0] >= e[1], e[1] >= e[2]);
    }
    if (tenVerbose && !( e[0] >= 0 && e[1] >= 0 && e[2] >= 0 )) {
!     fprintf(stderr, "tenAnisoCalc: eigen values not all >= 0: %g %g %g\n",
  	    e[0], e[1], e[2]);
    }
--- 42,51 ----
  
    if (!( e[0] >= e[1] && e[1] >= e[2] )) {
!     fprintf(stderr, "tenAnisoCalc_f: eigen values not sorted: "
  	    "%g %g %g (%d %d)\n",
  	    e[0], e[1], e[2], e[0] >= e[1], e[1] >= e[2]);
    }
    if (tenVerbose && !( e[0] >= 0 && e[1] >= 0 && e[2] >= 0 )) {
!     fprintf(stderr, "tenAnisoCalc_f: eigen values not all >= 0: %g %g %g\n",
  	    e[0], e[1], e[2]);
    }
***************
*** 156,160 ****
        e[2] = c0*m0[2] + c1*m1[2] + c2*m2[2];
        ELL_SORT3(e[0], e[1], e[2], tmp); /* got some warnings w/out this */
!       tenAnisoCalc(c, e);
        out[x + res*y] = c[aniso];
      }
--- 154,158 ----
        e[2] = c0*m0[2] + c1*m1[2] + c2*m2[2];
        ELL_SORT3(e[0], e[1], e[2], tmp); /* got some warnings w/out this */
!       tenAnisoCalc_f(c, e);
        out[x + res*y] = c[aniso];
      }
***************
*** 170,174 ****
  
  int
! tenAnisoVolume(Nrrd *nout, Nrrd *nin, int aniso, float thresh) {
    char me[]="tenAnisoVolume", err[AIR_STRLEN_MED];
    size_t N, I, copyI;
--- 168,172 ----
  
  int
! tenAnisoVolume(Nrrd *nout, Nrrd *nin, int aniso, double thresh) {
    char me[]="tenAnisoVolume", err[AIR_STRLEN_MED];
    size_t N, I, copyI;
***************
*** 203,207 ****
        continue;
      }
!     tenEigensolve(eval, evec, tensor);
      if (!(AIR_EXISTS(eval[0]) && AIR_EXISTS(eval[1]) && AIR_EXISTS(eval[2]))) {
        copyI = I;
--- 201,205 ----
        continue;
      }
!     tenEigensolve_f(eval, evec, tensor);
      if (!(AIR_EXISTS(eval[0]) && AIR_EXISTS(eval[1]) && AIR_EXISTS(eval[2]))) {
        copyI = I;
***************
*** 213,217 ****
        biffAdd(TEN, err); return 1;
      }
!     tenAnisoCalc(c, eval);
      out[I] = c[aniso];
    }
--- 211,215 ----
        biffAdd(TEN, err); return 1;
      }
!     tenAnisoCalc_f(c, eval);
      out[I] = c[aniso];
    }
***************
*** 262,267 ****
    N = nrrdElementNumber(nin)/7;
    for (I=0; I<=N-1; I++) {
!     tenEigensolve(eval, evec, tdata);
!     tenAnisoCalc(c, eval);
      cl = c[clIdx];
      cp = c[cpIdx];
--- 260,265 ----
    N = nrrdElementNumber(nin)/7;
    for (I=0; I<=N-1; I++) {
!     tenEigensolve_f(eval, evec, tdata);
!     tenAnisoCalc_f(c, eval);
      cl = c[clIdx];
      cp = c[cpIdx];

Index: glyph.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/glyph.c,v
retrieving revision 1.36
retrieving revision 1.37
diff -C2 -d -r1.36 -r1.37
*** glyph.c	1 Mar 2004 21:59:03 -0000	1.36
--- glyph.c	11 May 2004 02:39:54 -0000	1.37
***************
*** 264,270 ****
        }
      }
!     tenEigensolve(eval, evec, tdata);
      ELL_3M_TRANSPOSE(rotEvec, evec);
!     tenAnisoCalc(aniso, eval);
      if (parm->doSlice
  	&& pI[parm->sliceAxis] == parm->slicePos) {
--- 264,270 ----
        }
      }
!     tenEigensolve_f(eval, evec, tdata);
      ELL_3M_TRANSPOSE(rotEvec, evec);
!     tenAnisoCalc_f(aniso, eval);
      if (parm->doSlice
  	&& pI[parm->sliceAxis] == parm->slicePos) {

Index: tendEvec.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/tendEvec.c,v
retrieving revision 1.16
retrieving revision 1.17
diff -C2 -d -r1.16 -r1.17
*** tendEvec.c	13 Mar 2004 20:03:11 -0000	1.16
--- tendEvec.c	11 May 2004 02:39:54 -0000	1.17
***************
*** 88,92 ****
    if (1 == compLen) {
      for (I=0; I<N; I++) {
!       tenEigensolve(eval, evec, tdata);
        scl = tdata[0] >= thresh;
        ELL_3V_SCALE(edata, scl, evec+3*comp[0]);
--- 88,92 ----
    if (1 == compLen) {
      for (I=0; I<N; I++) {
!       tenEigensolve_f(eval, evec, tdata);
        scl = tdata[0] >= thresh;
        ELL_3V_SCALE(edata, scl, evec+3*comp[0]);
***************
*** 96,100 ****
    } else {
      for (I=0; I<N; I++) {
!       tenEigensolve(eval, evec, tdata);
        scl = tdata[0] >= thresh;
        for (cc=0; cc<compLen; cc++) {
--- 96,100 ----
    } else {
      for (I=0; I<N; I++) {
!       tenEigensolve_f(eval, evec, tdata);
        scl = tdata[0] >= thresh;
        for (cc=0; cc<compLen; cc++) {

Index: mod.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/mod.c,v
retrieving revision 1.10
retrieving revision 1.11
diff -C2 -d -r1.10 -r1.11
*** mod.c	7 Jan 2004 15:34:31 -0000	1.10
--- mod.c	11 May 2004 02:39:54 -0000	1.11
***************
*** 23,28 ****
  
  int
! tenSizeNormalize(Nrrd *nout, Nrrd *nin, float _weight[3],
! 		 float amount, float target) {
    char me[]="tenSizeNormalize", err[AIR_STRLEN_MED];
    float *tin, *tout, eval[3], evec[9], size, weight[3];
--- 23,28 ----
  
  int
! tenSizeNormalize(Nrrd *nout, Nrrd *nin, double _weight[3],
! 		 double amount, double target) {
    char me[]="tenSizeNormalize", err[AIR_STRLEN_MED];
    float *tin, *tout, eval[3], evec[9], size, weight[3];
***************
*** 59,63 ****
    N = nrrdElementNumber(nin)/7;
    for (I=0; I<=N-1; I++) {
!     tenEigensolve(eval, evec, tin);
      size = (weight[0]*AIR_ABS(eval[0])
  	    + weight[1]*AIR_ABS(eval[1])
--- 59,63 ----
    N = nrrdElementNumber(nin)/7;
    for (I=0; I<=N-1; I++) {
!     tenEigensolve_f(eval, evec, tin);
      size = (weight[0]*AIR_ABS(eval[0])
  	    + weight[1]*AIR_ABS(eval[1])
***************
*** 73,77 ****
      fprintf(stderr, "%g %g %g\n", eval[0], eval[1], eval[2]);
      */
!     tenMakeOne(tout, tin[0], eval, evec);
      tin += 7;
      tout += 7;
--- 73,77 ----
      fprintf(stderr, "%g %g %g\n", eval[0], eval[1], eval[2]);
      */
!     tenMakeOne_f(tout, tin[0], eval, evec);
      tin += 7;
      tout += 7;
***************
*** 82,86 ****
  
  int
! tenSizeScale(Nrrd *nout, Nrrd *nin, float amount) {
    char me[]="tenSizeScale", err[AIR_STRLEN_MED];
    size_t I, N;
--- 82,86 ----
  
  int
! tenSizeScale(Nrrd *nout, Nrrd *nin, double amount) {
    char me[]="tenSizeScale", err[AIR_STRLEN_MED];
    size_t I, N;
***************
*** 125,129 ****
  */
  int
! tenAnisoScale(Nrrd *nout, Nrrd *nin, float scale,
  	      int fixDet, int makePositive) {
    char me[]="tenAnisoScale", err[AIR_STRLEN_MED];
--- 125,129 ----
  */
  int
! tenAnisoScale(Nrrd *nout, Nrrd *nin, double scale,
  	      int fixDet, int makePositive) {
    char me[]="tenAnisoScale", err[AIR_STRLEN_MED];
***************
*** 149,153 ****
    N = nrrdElementNumber(nin)/7;
    for (I=0; I<N; I++) {
!     tenEigensolve(eval, evec, tin);
      if (fixDet) {
        eval[0] = AIR_MAX(eval[0], 0.00001);
--- 149,153 ----
    N = nrrdElementNumber(nin)/7;
    for (I=0; I<N; I++) {
!     tenEigensolve_f(eval, evec, tin);
      if (fixDet) {
        eval[0] = AIR_MAX(eval[0], 0.00001);
***************
*** 172,176 ****
        eval[2] = AIR_MAX(eval[2], 0.0);
      }
!     tenMakeOne(tout, tin[0], eval, evec);
      tin += 7;
      tout += 7;
--- 172,176 ----
        eval[2] = AIR_MAX(eval[2], 0.0);
      }
!     tenMakeOne_f(tout, tin[0], eval, evec);
      tin += 7;
      tout += 7;
***************
*** 185,189 ****
  */
  int
! tenEigenvalueClamp(Nrrd *nout, Nrrd *nin, float min, float max) {
    char me[]="tenEigenvalueClamp", err[AIR_STRLEN_MED];
    size_t I, N;
--- 185,189 ----
  */
  int
! tenEigenvalueClamp(Nrrd *nout, Nrrd *nin, double min, double max) {
    char me[]="tenEigenvalueClamp", err[AIR_STRLEN_MED];
    size_t I, N;
***************
*** 208,212 ****
    N = nrrdElementNumber(nin)/7;
    for (I=0; I<N; I++) {
!     tenEigensolve(eval, evec, tin);
      if (AIR_EXISTS(min)) {
        eval[0] = AIR_MAX(eval[0], min);
--- 208,212 ----
    N = nrrdElementNumber(nin)/7;
    for (I=0; I<N; I++) {
!     tenEigensolve_f(eval, evec, tin);
      if (AIR_EXISTS(min)) {
        eval[0] = AIR_MAX(eval[0], min);
***************
*** 219,223 ****
        eval[2] = AIR_MIN(eval[2], max);
      }
!     tenMakeOne(tout, tin[0], eval, evec);
      tin += 7;
      tout += 7;
--- 219,223 ----
        eval[2] = AIR_MIN(eval[2], max);
      }
!     tenMakeOne_f(tout, tin[0], eval, evec);
      tin += 7;
      tout += 7;
***************
*** 232,236 ****
  */
  int
! tenEigenvaluePower(Nrrd *nout, Nrrd *nin, float expo) {
    char me[]="tenEigenvaluePower", err[AIR_STRLEN_MED];
    size_t I, N;
--- 232,236 ----
  */
  int
! tenEigenvaluePower(Nrrd *nout, Nrrd *nin, double expo) {
    char me[]="tenEigenvaluePower", err[AIR_STRLEN_MED];
    size_t I, N;
***************
*** 255,263 ****
    N = nrrdElementNumber(nin)/7;
    for (I=0; I<N; I++) {
!     tenEigensolve(eval, evec, tin);
      eval[0] = pow(eval[0], expo);
      eval[1] = pow(eval[1], expo);
      eval[2] = pow(eval[2], expo);
!     tenMakeOne(tout, tin[0], eval, evec);
      tin += 7;
      tout += 7;
--- 255,263 ----
    N = nrrdElementNumber(nin)/7;
    for (I=0; I<N; I++) {
!     tenEigensolve_f(eval, evec, tin);
      eval[0] = pow(eval[0], expo);
      eval[1] = pow(eval[1], expo);
      eval[2] = pow(eval[2], expo);
!     tenMakeOne_f(tout, tin[0], eval, evec);
      tin += 7;
      tout += 7;
***************
*** 272,276 ****
  */
  int
! tenEigenvalueAdd(Nrrd *nout, Nrrd *nin, float val) {
    char me[]="tenEigenvalueAdd", err[AIR_STRLEN_MED];
    size_t I, N;
--- 272,276 ----
  */
  int
! tenEigenvalueAdd(Nrrd *nout, Nrrd *nin, double val) {
    char me[]="tenEigenvalueAdd", err[AIR_STRLEN_MED];
    size_t I, N;
***************
*** 295,303 ****
    N = nrrdElementNumber(nin)/7;
    for (I=0; I<N; I++) {
!     tenEigensolve(eval, evec, tin);
      eval[0] += val;
      eval[1] += val;
      eval[2] += val;
!     tenMakeOne(tout, tin[0], eval, evec);
      tin += 7;
      tout += 7;
--- 295,303 ----
    N = nrrdElementNumber(nin)/7;
    for (I=0; I<N; I++) {
!     tenEigensolve_f(eval, evec, tin);
      eval[0] += val;
      eval[1] += val;
      eval[2] += val;
!     tenMakeOne_f(tout, tin[0], eval, evec);
      tin += 7;
      tout += 7;

Index: chan.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/chan.c,v
retrieving revision 1.31
retrieving revision 1.32
diff -C2 -d -r1.31 -r1.32
*** chan.c	7 Mar 2004 00:52:57 -0000	1.31
--- chan.c	11 May 2004 02:39:54 -0000	1.32
***************
*** 129,133 ****
  
  /*
! ******** tenEstimateSingle
  **
  ** estimate one single tensor
--- 129,133 ----
  
  /*
! ******** tenEstimateSingle_f
  **
  ** estimate one single tensor
***************
*** 153,159 ****
  */
  void
! tenEstimateLinearSingle(float *ten, float *B0P, float *dwi, double *emat,
! 			int DD, int knownB0,
! 			float thresh, float soft, float b) {
    double logB0, v[_TEN_MAX_DWI_NUM], tmp, mean;
    int i, j;
--- 153,159 ----
  */
  void
! tenEstimateLinearSingle_f(float *ten, float *B0P, float *dwi, double *emat,
! 			  int DD, int knownB0,
! 			  float thresh, float soft, float b) {
    double logB0, v[_TEN_MAX_DWI_NUM], tmp, mean;
    int i, j;
***************
*** 223,227 ****
  		    Nrrd **_ndwi, int dwiLen, 
  		    Nrrd *_nbmat, int knownB0, 
! 		    float thresh, float soft, float b) {
    char me[]="tenEstimateLinear3D", err[AIR_STRLEN_MED];
    Nrrd *ndwi;
--- 223,227 ----
  		    Nrrd **_ndwi, int dwiLen, 
  		    Nrrd *_nbmat, int knownB0, 
! 		    double thresh, double soft, double b) {
    char me[]="tenEstimateLinear3D", err[AIR_STRLEN_MED];
    Nrrd *ndwi;
***************
*** 269,273 ****
  tenEstimateLinear4D(Nrrd *nten, Nrrd **nterrP, Nrrd **nB0P,
  		    Nrrd *ndwi, Nrrd *_nbmat, int knownB0,
! 		    float thresh, float soft, float b) {
    char me[]="tenEstimateLinear4D", err[AIR_STRLEN_MED];
    Nrrd *nemat, *nbmat, *ncrop, *nhist;
--- 269,273 ----
  tenEstimateLinear4D(Nrrd *nten, Nrrd **nterrP, Nrrd **nB0P,
  		    Nrrd *ndwi, Nrrd *_nbmat, int knownB0,
! 		    double thresh, double soft, double b) {
    char me[]="tenEstimateLinear4D", err[AIR_STRLEN_MED];
    Nrrd *nemat, *nbmat, *ncrop, *nhist;
***************
*** 388,393 ****
  	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;
--- 388,393 ----
  	fprintf(stderr, "%s: input dwi1[%d] = %g\n", me, d, dwi1[d]);
      }
!     tenEstimateLinearSingle_f(ten, &_B0,
! 			      dwi1, emat, DD, knownB0, thresh, soft, b);
      if (nB0P) {
        *B0 = _B0;
***************
*** 400,404 ****
        te = 0;
        if (knownB0) {
! 	tenSimulateOne(dwi2, _B0, ten, bmat, DD, b);
  	for (d=1; d<DD; d++) {
  	  d1 = AIR_MAX(dwi1[d], 1);
--- 400,404 ----
        te = 0;
        if (knownB0) {
! 	tenSimulateOne_f(dwi2, _B0, ten, bmat, DD, b);
  	for (d=1; d<DD; d++) {
  	  d1 = AIR_MAX(dwi1[d], 1);
***************
*** 408,415 ****
  	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 */
--- 408,415 ----
  	te /= (DD-1);
        } else {
! 	tenSimulateOne_f(dwi2, _B0, ten, bmat, DD+1, b);
  	for (d=0; d<DD; d++) {
  	  d1 = AIR_MAX(dwi1[d], 1);
! 	  /* tenSimulateOne_f 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 */
***************
*** 437,441 ****
  
  /*
! ******** tenSimulateOne
  **
  ** given a tensor, simulate the set of diffusion weighted measurements
--- 437,441 ----
  
  /*
! ******** tenSimulateOne_f
  **
  ** given a tensor, simulate the set of diffusion weighted measurements
***************
*** 448,453 ****
  */
  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 */
--- 448,453 ----
  */
  void
! tenSimulateOne_f(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 */
***************
*** 475,479 ****
  
  int
! tenSimulate(Nrrd *ndwi, Nrrd *nT2, Nrrd *nten, Nrrd *_nbmat, float b) {
    char me[]="tenSimulate", err[AIR_STRLEN_MED];
    size_t II;
--- 475,479 ----
  
  int
! tenSimulate(Nrrd *ndwi, Nrrd *nT2, Nrrd *nten, Nrrd *_nbmat, double b) {
    char me[]="tenSimulate", err[AIR_STRLEN_MED];
    size_t II;
***************
*** 520,524 ****
    for (II=0; II<sx*sy*sz; II++) {
      /* tenVerbose = (II == 42 + 190*(96 + 196*0)); */
!     tenSimulateOne(dwi, lup(nT2->data, II), ten, bmat, DD, b);
      dwi += DD;
      ten += 7;
--- 520,524 ----
    for (II=0; II<sx*sy*sz; II++) {
      /* tenVerbose = (II == 42 + 190*(96 + 196*0)); */
!     tenSimulateOne_f(dwi, lup(nT2->data, II), ten, bmat, DD, b);
      dwi += DD;
      ten += 7;

Index: miscTen.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/miscTen.c,v
retrieving revision 1.24
retrieving revision 1.25
diff -C2 -d -r1.24 -r1.25
*** miscTen.c	30 Mar 2004 19:25:14 -0000	1.24
--- miscTen.c	11 May 2004 02:39:54 -0000	1.25
***************
*** 24,28 ****
  int
  tenEvecRGB(Nrrd *nout, Nrrd *nin, int which, int aniso,
! 	   float gamma, float bgGray, float isoGray) {
    char me[]="tenEvecRGB", err[AIR_STRLEN_MED];
    int size[NRRD_DIM_MAX];
--- 24,28 ----
  int
  tenEvecRGB(Nrrd *nout, Nrrd *nin, int which, int aniso,
! 	   double gamma, double bgGray, double isoGray) {
    char me[]="tenEvecRGB", err[AIR_STRLEN_MED];
    int size[NRRD_DIM_MAX];
***************
*** 58,63 ****
    for (II=0; II<NN; II++) {
      /* tenVerbose = (II == (50 + 64*(32 + 64*0))); */
!     tenEigensolve(eval, evec, tdata);
!     tenAnisoCalc(an, eval);
      R = AIR_ABS(evec[0 + 3*which]);
      G = AIR_ABS(evec[1 + 3*which]);
--- 58,63 ----
    for (II=0; II<NN; II++) {
      /* tenVerbose = (II == (50 + 64*(32 + 64*0))); */
!     tenEigensolve_f(eval, evec, tdata);
!     tenAnisoCalc_f(an, eval);
      R = AIR_ABS(evec[0 + 3*which]);
      G = AIR_ABS(evec[1 + 3*which]);
***************
*** 163,169 ****
    qdata = nout->data;
    for (I=0; I<N; I++) {
!     tenEigensolve(eval, evec, tdata);
      if (scaleByAniso) {
!       tenAnisoCalc(c, eval);
        an = c[aniso];
      } else {
--- 163,169 ----
    qdata = nout->data;
    for (I=0; I<N; I++) {
!     tenEigensolve_f(eval, evec, tdata);
      if (scaleByAniso) {
!       tenAnisoCalc_f(c, eval);
        an = c[aniso];
      } else {
***************
*** 215,219 ****
  */
  int
! _tenFindValley(float *valP, Nrrd *nhist, float tweak, int save) {
    char me[]="_tenFindValley", err[AIR_STRLEN_MED];
    double gparm[NRRD_KERNEL_PARMS_NUM], dparm[NRRD_KERNEL_PARMS_NUM];
--- 215,219 ----
  */
  int
! _tenFindValley(double *valP, Nrrd *nhist, double tweak, int save) {
    char me[]="_tenFindValley", err[AIR_STRLEN_MED];
    double gparm[NRRD_KERNEL_PARMS_NUM], dparm[NRRD_KERNEL_PARMS_NUM];

Index: tenGage.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/tenGage.c,v
retrieving revision 1.19
retrieving revision 1.20
diff -C2 -d -r1.19 -r1.20
*** tenGage.c	9 May 2004 22:55:16 -0000	1.19
--- tenGage.c	11 May 2004 02:39:54 -0000	1.20
***************
*** 163,167 ****
  #if !GAGE_TYPE_FLOAT
    int ci;
!   float tenAnsF[7], evalAnsF[3], evecAnsF[9], aniso[TEN_ANISO_MAX+1];
  #endif
  
--- 163,167 ----
  #if !GAGE_TYPE_FLOAT
    int ci;
!   float evalAnsF[3], aniso[TEN_ANISO_MAX+1];
  #endif
  
***************
*** 221,239 ****
         we always find the eigenvalues, whether or not they were asked for */
  #if GAGE_TYPE_FLOAT
!     tenEigensolve(evalAns, evecAns, tenAns);
  #else
!     TEN_T_COPY(tenAnsF, tenAns);
!     tenEigensolve(evalAnsF, evecAnsF, tenAnsF);
!     ELL_3V_COPY(evalAns, evalAnsF);
!     ELL_3M_COPY(evecAns, evecAnsF);
  #endif
    } else if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageEval)) {
      /* else eigenvectors are NOT needed, but eigenvalues ARE needed */
  #if GAGE_TYPE_FLOAT
!     tenEigensolve(evalAns, NULL, tenAns);
  #else
!     TEN_T_COPY(tenAnsF, tenAns);
!     tenEigensolve(evalAnsF, NULL, tenAnsF);
!     ELL_3V_COPY(evalAns, evalAnsF);
  #endif
    }
--- 221,234 ----
         we always find the eigenvalues, whether or not they were asked for */
  #if GAGE_TYPE_FLOAT
!     tenEigensolve_f(evalAns, evecAns, tenAns);
  #else
!     tenEigensolve_d(evalAns, evecAns, tenAns);
  #endif
    } else if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageEval)) {
      /* else eigenvectors are NOT needed, but eigenvalues ARE needed */
  #if GAGE_TYPE_FLOAT
!     tenEigensolve_f(evalAns, NULL, tenAns);
  #else
!     tenEigensolve_d(evalAns, NULL, tenAns);
  #endif
    }
***************
*** 445,451 ****
    if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageAniso)) {
  #if GAGE_TYPE_FLOAT
!     tenAnisoCalc(pvl->directAnswer[tenGageAniso], evalAns);
  #else
!     tenAnisoCalc(aniso, evalAnsF);
      for (ci=0; ci<=TEN_ANISO_MAX; ci++) {
        pvl->directAnswer[tenGageAniso][ci] = aniso[ci];
--- 440,446 ----
    if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageAniso)) {
  #if GAGE_TYPE_FLOAT
!     tenAnisoCalc_f(pvl->directAnswer[tenGageAniso], evalAns);
  #else
!     tenAnisoCalc_f(aniso, evalAnsF);
      for (ci=0; ci<=TEN_ANISO_MAX; ci++) {
        pvl->directAnswer[tenGageAniso][ci] = aniso[ci];

Index: ten.h
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/ten.h,v
retrieving revision 1.84
retrieving revision 1.85
diff -C2 -d -r1.84 -r1.85
*** ten.h	22 Apr 2004 07:14:54 -0000	1.84
--- ten.h	11 May 2004 02:39:54 -0000	1.85
***************
*** 435,443 ****
  TEEM_API int tenVerbose;
  TEEM_API int tenTensorCheck(Nrrd *nin, int wantType, int want4D, int useBiff);
! TEEM_API int tenExpand(Nrrd *tnine, Nrrd *tseven, float scale, float thresh);
  TEEM_API int tenShrink(Nrrd *tseven, Nrrd *nconf, Nrrd *tnine);
! TEEM_API int tenEigensolve(float eval[3], float evec[9], float ten[7]);
! TEEM_API void tenMakeOne(float ten[7],
! 			 float conf, float eval[3], float evec[9]);
  TEEM_API int tenMake(Nrrd *nout, Nrrd *nconf, Nrrd *neval, Nrrd *nevec);
  TEEM_API int tenSlice(Nrrd *nout, Nrrd *nten, int axis, int pos, int dim);
--- 435,444 ----
  TEEM_API int tenVerbose;
  TEEM_API int tenTensorCheck(Nrrd *nin, int wantType, int want4D, int useBiff);
! TEEM_API int tenExpand(Nrrd *tnine, Nrrd *tseven, double scale, double thresh);
  TEEM_API int tenShrink(Nrrd *tseven, Nrrd *nconf, Nrrd *tnine);
! TEEM_API int tenEigensolve_f(float eval[3], float evec[9], float ten[7]);
! TEEM_API int tenEigensolve_d(double eval[3], double evec[9], double ten[7]);
! TEEM_API void tenMakeOne_f(float ten[7],
! 			   float conf, float eval[3], float evec[9]);
  TEEM_API int tenMake(Nrrd *nout, Nrrd *nconf, Nrrd *neval, Nrrd *nevec);
  TEEM_API int tenSlice(Nrrd *nout, Nrrd *nten, int axis, int pos, int dim);
***************
*** 447,471 ****
  TEEM_API int tenBMatrixCalc(Nrrd *nbmat, Nrrd *ngrad);
  TEEM_API int tenEMatrixCalc(Nrrd *nemat, Nrrd *nbmat, int knownB0);
! TEEM_API void tenEstimateLinearSingle(float *ten, float *B0P, float *dwi,
! 				      double *emat, int DD, int knownB0,
! 				      float thresh, float soft, float b);
  TEEM_API int tenEstimateLinear3D(Nrrd *nten, Nrrd **nterrP, Nrrd **nB0P,
  				 Nrrd **ndwi, int dwiLen, 
  				 Nrrd *nbmat, int knownB0, 
! 				 float thresh, float soft, float b);
  TEEM_API int tenEstimateLinear4D(Nrrd *nten, Nrrd **nterrP, Nrrd **nB0P,
  				 Nrrd *ndwi, Nrrd *_nbmat, int knownB0,
! 				 float thresh, float soft, float b);
! TEEM_API void tenSimulateOne(float *dwi, float B0, float *ten,
! 			     double *bmat, int DD, float b);
  TEEM_API int tenSimulate(Nrrd *ndwi, Nrrd *nT2, Nrrd *nten,
! 			 Nrrd *nbmat, float b);
  
  /* aniso.c */
! TEEM_API float tenAnisoSigma;  /* added to denominator in Westin anisos */
! TEEM_API void tenAnisoCalc(float c[TEN_ANISO_MAX+1], float eval[3]);
  TEEM_API int tenAnisoPlot(Nrrd *nout, int aniso, int res,
  			  int whole, int nanout);
! TEEM_API int tenAnisoVolume(Nrrd *nout, Nrrd *nin, int aniso, float thresh);
  TEEM_API int tenAnisoHistogram(Nrrd *nout, Nrrd *nin,
  			       int version, int resolution);
--- 448,471 ----
  TEEM_API int tenBMatrixCalc(Nrrd *nbmat, Nrrd *ngrad);
  TEEM_API int tenEMatrixCalc(Nrrd *nemat, Nrrd *nbmat, int knownB0);
! TEEM_API void tenEstimateLinearSingle_f(float *ten, float *B0P, float *dwi,
! 					double *emat, int DD, int knownB0,
! 					float thresh, float soft, float b);
  TEEM_API int tenEstimateLinear3D(Nrrd *nten, Nrrd **nterrP, Nrrd **nB0P,
  				 Nrrd **ndwi, int dwiLen, 
  				 Nrrd *nbmat, int knownB0, 
! 				 double thresh, double soft, double b);
  TEEM_API int tenEstimateLinear4D(Nrrd *nten, Nrrd **nterrP, Nrrd **nB0P,
  				 Nrrd *ndwi, Nrrd *_nbmat, int knownB0,
! 				 double thresh, double soft, double b);
! TEEM_API void tenSimulateOne_f(float *dwi, float B0, float *ten,
! 			       double *bmat, int DD, float b);
  TEEM_API int tenSimulate(Nrrd *ndwi, Nrrd *nT2, Nrrd *nten,
! 			 Nrrd *nbmat, double b);
  
  /* aniso.c */
! TEEM_API void tenAnisoCalc_f(float c[TEN_ANISO_MAX+1], float eval[3]);
  TEEM_API int tenAnisoPlot(Nrrd *nout, int aniso, int res,
  			  int whole, int nanout);
! TEEM_API int tenAnisoVolume(Nrrd *nout, Nrrd *nin, int aniso, double thresh);
  TEEM_API int tenAnisoHistogram(Nrrd *nout, Nrrd *nin,
  			       int version, int resolution);
***************
*** 473,483 ****
  /* miscTen.c */
  TEEM_API int tenEvecRGB(Nrrd *nout, Nrrd *nin, int which, int aniso,
! 			float gamma, float bgGray, float isoGray);
! TEEM_API short tenEvqOne(float vec[3], float scl);
  TEEM_API int tenEvqVolume(Nrrd *nout, Nrrd *nin, int which,
  			  int aniso, int scaleByAniso);
  TEEM_API int tenBMatrixCheck(Nrrd *nbmat);
! TEEM_API int _tenFindValley(float *valP, Nrrd *nhist,
! 			    float tweak, int save);
  
  /* fiberMethods.c */
--- 473,483 ----
  /* miscTen.c */
  TEEM_API int tenEvecRGB(Nrrd *nout, Nrrd *nin, int which, int aniso,
! 			double gamma, double bgGray, double isoGray);
! TEEM_API short tenEvqOne_f(float vec[3], float scl);
  TEEM_API int tenEvqVolume(Nrrd *nout, Nrrd *nin, int which,
  			  int aniso, int scaleByAniso);
  TEEM_API int tenBMatrixCheck(Nrrd *nbmat);
! TEEM_API int _tenFindValley(double *valP, Nrrd *nhist,
! 			    double tweak, int save);
  
  /* fiberMethods.c */
***************
*** 502,506 ****
  			      int dwiLen, Nrrd *ngrad,
  			      int reference,
! 			      float bwX, float bwY, float fitFrac, float DWthr,
  			      int doCC,
  			      const NrrdKernel *kern, double *kparm,
--- 502,507 ----
  			      int dwiLen, Nrrd *ngrad,
  			      int reference,
! 			      double bwX, double bwY,
! 			      double fitFrac, double DWthr,
  			      int doCC,
  			      const NrrdKernel *kern, double *kparm,
***************
*** 508,525 ****
  TEEM_API int tenEpiRegister4D(Nrrd *nout, Nrrd *nin, Nrrd *ngrad,
  			      int reference,
! 			      float bwX, float bwY, float fitFrac,
! 			      float DWthr, int doCC, 
  			      const NrrdKernel *kern, double *kparm,
  			      int progress, int verbose);
  
  /* mod.c */
! TEEM_API int tenSizeNormalize(Nrrd *nout, Nrrd *nin, float weight[3],
! 			      float amount, float target);
! TEEM_API int tenSizeScale(Nrrd *nout, Nrrd *nin, float amount);
! TEEM_API int tenAnisoScale(Nrrd *nout, Nrrd *nin, float amount,
  			   int fixDet, int makePositive);
! TEEM_API int tenEigenvaluePower(Nrrd *nout, Nrrd *nin, float expo);
! TEEM_API int tenEigenvalueClamp(Nrrd *nout, Nrrd *nin, float min, float max);
! TEEM_API int tenEigenvalueAdd(Nrrd *nout, Nrrd *nin, float val);
  
  /* tenGage.c */
--- 509,527 ----
  TEEM_API int tenEpiRegister4D(Nrrd *nout, Nrrd *nin, Nrrd *ngrad,
  			      int reference,
! 			      double bwX, double bwY,
! 			      double fitFrac, double DWthr,
! 			      int doCC, 
  			      const NrrdKernel *kern, double *kparm,
  			      int progress, int verbose);
  
  /* mod.c */
! TEEM_API int tenSizeNormalize(Nrrd *nout, Nrrd *nin, double weight[3],
! 			      double amount, double target);
! TEEM_API int tenSizeScale(Nrrd *nout, Nrrd *nin, double amount);
! TEEM_API int tenAnisoScale(Nrrd *nout, Nrrd *nin, double scale,
  			   int fixDet, int makePositive);
! TEEM_API int tenEigenvaluePower(Nrrd *nout, Nrrd *nin, double expo);
! TEEM_API int tenEigenvalueClamp(Nrrd *nout, Nrrd *nin, double min, double max);
! TEEM_API int tenEigenvalueAdd(Nrrd *nout, Nrrd *nin, double val);
  
  /* tenGage.c */

Index: tendPoint.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/tendPoint.c,v
retrieving revision 1.10
retrieving revision 1.11
diff -C2 -d -r1.10 -r1.11
*** tendPoint.c	22 Apr 2004 07:46:52 -0000	1.10
--- tendPoint.c	11 May 2004 02:39:54 -0000	1.11
***************
*** 77,81 ****
    fprintf(stderr, "% 15.7f % 15.7f % 15.7f\n", tdata[2], tdata[4], tdata[5]);
    fprintf(stderr, "% 15.7f % 15.7f % 15.7f\n", tdata[3], tdata[5], tdata[6]);
!   tenEigensolve(eval, evec, tdata);
    fprintf(stderr, "eigensystem = (<eigenvalue> : <eigenvector>):\n");
    fprintf(stderr, "% 15.7f : % 15.7f % 15.7f % 15.7f\n",
--- 77,81 ----
    fprintf(stderr, "% 15.7f % 15.7f % 15.7f\n", tdata[2], tdata[4], tdata[5]);
    fprintf(stderr, "% 15.7f % 15.7f % 15.7f\n", tdata[3], tdata[5], tdata[6]);
!   tenEigensolve_f(eval, evec, tdata);
    fprintf(stderr, "eigensystem = (<eigenvalue> : <eigenvector>):\n");
    fprintf(stderr, "% 15.7f : % 15.7f % 15.7f % 15.7f\n",
***************
*** 95,99 ****
    fprintf(stderr, "% 15.7f % 15.7f % 15.7f\n",
  	  mat[6], mat[7], mat[8]);
!   tenAnisoCalc(c, eval);
    fprintf(stderr, "anisotropies = \n");
    for (i=1; i<=TEN_ANISO_MAX; i++) {
--- 95,99 ----
    fprintf(stderr, "% 15.7f % 15.7f % 15.7f\n",
  	  mat[6], mat[7], mat[8]);
!   tenAnisoCalc_f(c, eval);
    fprintf(stderr, "anisotropies = \n");
    for (i=1; i<=TEN_ANISO_MAX; i++) {

Index: tendNorm.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/tendNorm.c,v
retrieving revision 1.4
retrieving revision 1.5
diff -C2 -d -r1.4 -r1.5
*** tendNorm.c	7 Jan 2004 15:34:31 -0000	1.4
--- tendNorm.c	11 May 2004 02:39:54 -0000	1.5
***************
*** 37,43 ****
    Nrrd *nin, *nout;
    char *outS;
!   float weight[3], amount, target;
  
!   hestOptAdd(&hopt, "w", "w0 w1 w2", airTypeFloat, 3, 3, weight, NULL,
  	     "relative weights to put on major, medium, and minor "
  	     "eigenvalue when performing normalization (internally "
--- 37,44 ----
    Nrrd *nin, *nout;
    char *outS;
!   float amount, target;
!   double weight[3];
  
!   hestOptAdd(&hopt, "w", "w0 w1 w2", airTypeDouble, 3, 3, weight, NULL,
  	     "relative weights to put on major, medium, and minor "
  	     "eigenvalue when performing normalization (internally "

Index: epireg.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/epireg.c,v
retrieving revision 1.28
retrieving revision 1.29
diff -C2 -d -r1.28 -r1.29
*** epireg.c	13 Mar 2004 20:03:11 -0000	1.28
--- epireg.c	11 May 2004 02:39:54 -0000	1.29
***************
*** 191,195 ****
  
  int
! _tenEpiRegFindThresh(float *DWthrP, Nrrd **nin, int ninLen, int save) {
    char me[]="_tenEpiRegFindThresh", err[AIR_STRLEN_MED];
    Nrrd *nhist, *ntmp;
--- 191,195 ----
  
  int
! _tenEpiRegFindThresh(double *DWthrP, Nrrd **nin, int ninLen, int save) {
    char me[]="_tenEpiRegFindThresh", err[AIR_STRLEN_MED];
    Nrrd *nhist, *ntmp;
***************
*** 247,251 ****
  int
  _tenEpiRegThreshold(Nrrd **nthresh, Nrrd **nblur, int ninLen,
! 		    float DWthr, int verb, int progress) {
    char me[]="_tenEpiRegThreshold", err[AIR_STRLEN_MED];
    airArray *mop;
--- 247,251 ----
  int
  _tenEpiRegThreshold(Nrrd **nthresh, Nrrd **nblur, int ninLen,
! 		    double DWthr, int verb, int progress) {
    char me[]="_tenEpiRegThreshold", err[AIR_STRLEN_MED];
    airArray *mop;
***************
*** 999,1004 ****
  tenEpiRegister3D(Nrrd **nout, Nrrd **nin, int ninLen, Nrrd *_ngrad,
  		 int reference,
! 		 float bwX, float bwY, float fitFrac,
! 		 float DWthr, int doCC, 
  		 const NrrdKernel *kern, double *kparm,
  		 int progress, int verbose) {
--- 999,1004 ----
  tenEpiRegister3D(Nrrd **nout, Nrrd **nin, int ninLen, Nrrd *_ngrad,
  		 int reference,
! 		 double bwX, double bwY, double fitFrac,
! 		 double DWthr, int doCC, 
  		 const NrrdKernel *kern, double *kparm,
  		 int progress, int verbose) {
***************
*** 1135,1140 ****
  tenEpiRegister4D(Nrrd *_nout, Nrrd *_nin, Nrrd *ngrad,
  		 int reference,
! 		 float bwX, float bwY, float fitFrac,
! 		 float DWthr, int doCC, 
  		 const NrrdKernel *kern, double *kparm,
  		 int progress, int verbose) {
--- 1135,1140 ----
  tenEpiRegister4D(Nrrd *_nout, Nrrd *_nin, Nrrd *ngrad,
  		 int reference,
! 		 double bwX, double bwY, double fitFrac,
! 		 double DWthr, int doCC, 
  		 const NrrdKernel *kern, double *kparm,
  		 int progress, int verbose) {

Index: tendEval.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/tendEval.c,v
retrieving revision 1.15
retrieving revision 1.16
diff -C2 -d -r1.15 -r1.16
*** tendEval.c	13 Mar 2004 20:03:11 -0000	1.15
--- tendEval.c	11 May 2004 02:39:54 -0000	1.16
***************
*** 93,97 ****
      ELL_3V_SET(map, 1, 2, 3);
      for (I=0; I<N; I++) {
!       tenEigensolve(eval, evec, tdata);
        edata[I] = (tdata[0] >= thresh)*eval[comp[0]];
        tdata += 7;
--- 93,97 ----
      ELL_3V_SET(map, 1, 2, 3);
      for (I=0; I<N; I++) {
!       tenEigensolve_f(eval, evec, tdata);
        edata[I] = (tdata[0] >= thresh)*eval[comp[0]];
        tdata += 7;
***************
*** 100,104 ****
      ELL_4V_SET(map, 0, 1, 2, 3);
      for (I=0; I<N; I++) {
!       tenEigensolve(eval, evec, tdata);
        for (cc=0; cc<compLen; cc++)
  	edata[cc] = (tdata[0] >= thresh)*eval[comp[cc]];
--- 100,104 ----
      ELL_4V_SET(map, 0, 1, 2, 3);
      for (I=0; I<N; I++) {
!       tenEigensolve_f(eval, evec, tdata);
        for (cc=0; cc<compLen; cc++)
  	edata[cc] = (tdata[0] >= thresh)*eval[comp[cc]];

Index: tensor.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/tensor.c,v
retrieving revision 1.30
retrieving revision 1.31
diff -C2 -d -r1.30 -r1.31
*** tensor.c	7 May 2004 19:13:34 -0000	1.30
--- tensor.c	11 May 2004 02:39:54 -0000	1.31
***************
*** 73,77 ****
  
  int
! tenExpand(Nrrd *nout, Nrrd *nin, float scale, float thresh) {
    char me[]="tenExpand", err[AIR_STRLEN_MED];
    size_t N, I;
--- 73,77 ----
  
  int
! tenExpand(Nrrd *nout, Nrrd *nin, double scale, double thresh) {
    char me[]="tenExpand", err[AIR_STRLEN_MED];
    size_t N, I;
***************
*** 179,183 ****
  
  /*
! ******** tenEigensolve
  **
  ** uses ell_3m_eigensolve_d to get the eigensystem of a single tensor
--- 179,183 ----
  
  /*
! ******** tenEigensolve_f
  **
  ** uses ell_3m_eigensolve_d to get the eigensystem of a single tensor
***************
*** 194,198 ****
  */
  int
! tenEigensolve(float _eval[3], float _evec[9], float t[7]) {
    double m[9], eval[3], evec[9], trc, iso[9];
    int ret;
--- 194,198 ----
  */
  int
! tenEigensolve_f(float _eval[3], float _evec[9], float t[7]) {
    double m[9], eval[3], evec[9], trc, iso[9];
    int ret;
***************
*** 201,206 ****
    trc = ELL_3M_TRACE(m)/3.0;
    ELL_3M_IDENTITY_SET(iso);
!   ELL_3M_SCALE(iso, trc, iso);
!   ELL_3M_SUB(m, m, iso);
    if (_evec) {
      ret = ell_3m_eigensolve_d(eval, evec, m, AIR_TRUE);
--- 201,206 ----
    trc = ELL_3M_TRACE(m)/3.0;
    ELL_3M_IDENTITY_SET(iso);
!   ELL_3M_SCALE_SET(iso, -trc, -trc, -trc);
!   ELL_3M_ADD2(m, m, iso);
    if (_evec) {
      ret = ell_3m_eigensolve_d(eval, evec, m, AIR_TRUE);
***************
*** 249,252 ****
--- 249,283 ----
  }
  
+ int
+ tenEigensolve_d(double _eval[3], double evec[9], double t[7]) {
+   double m[9], eval[3], trc, iso[9];
+   int ret;
+   
+   TEN_T2M(m, t);
+   trc = ELL_3M_TRACE(m)/3.0;
+   ELL_3M_SCALE_SET(iso, -trc, -trc, -trc);
+   ELL_3M_ADD2(m, m, iso);
+   if (evec) {
+     ret = ell_3m_eigensolve_d(eval, evec, m, AIR_TRUE);
+     ELL_3V_SET(_eval, eval[0] + trc, eval[1] + trc, eval[2] + trc);
+     if (ell_cubic_root_single_double == ret) {
+       /* this was added to fix a stupid problem with very nearly
+ 	 isotropic glyphs, used for demonstration figures */
+       if (eval[0] == eval[1]) {
+ 	ELL_3V_CROSS(evec+6, evec+0, evec+3);
+       } else {
+ 	ELL_3V_CROSS(evec+0, evec+3, evec+6);
+       }
+     }
+   } else {
+     /* caller only wants eigenvalues */
+     ret = ell_3m_eigenvalues_d(eval, m, AIR_TRUE);
+     ELL_3V_SET(_eval, eval[0] + trc, eval[1] + trc, eval[2] + trc);
+   }    
+   return ret;
+ }
+ 
+ 
+ 
  /*  lop A
      fprintf(stderr, "###################################  I = %d\n", (int)I);
***************
*** 289,293 ****
  
  void
! tenMakeOne(float ten[7], float conf, float eval[3], float evec[9]) {
    double tmpMat1[9], tmpMat2[9], diag[9], evecT[9];
  
--- 320,324 ----
  
  void
! tenMakeOne_f(float ten[7], float conf, float eval[3], float evec[9]) {
    double tmpMat1[9], tmpMat2[9], diag[9], evecT[9];
  
***************
*** 378,382 ****
    out = (float *)nout->data;
    for (I=0; I<N; I++) {
!     tenMakeOne(out, conf[I], eval, evec);
      /* lop A */
      out += 7;
--- 409,413 ----
    out = (float *)nout->data;
    for (I=0; I<N; I++) {
!     tenMakeOne_f(out, conf[I], eval, evec);
      /* lop A */
      out += 7;