1 /*
2   Teem: Tools to process and visualize scientific data and images             .
3   Copyright (C) 2012, 2011, 2010, 2009  University of Chicago
4   Copyright (C) 2008, 2007, 2006, 2005  Gordon Kindlmann
5   Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998  University of Utah
6 
7   This library is free software; you can redistribute it and/or
8   modify it under the terms of the GNU Lesser General Public License
9   (LGPL) as published by the Free Software Foundation; either
10   version 2.1 of the License, or (at your option) any later version.
11   The terms of redistributing and/or modifying this software also
12   include exceptions to the LGPL that facilitate static linking.
13 
14   This library is distributed in the hope that it will be useful,
15   but WITHOUT ANY WARRANTY; without even the implied warranty of
16   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17   Lesser General Public License for more details.
18 
19   You should have received a copy of the GNU Lesser General Public License
20   along with this library; if not, write to Free Software Foundation, Inc.,
21   51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
22 */
23 
24 #ifndef TEN_HAS_BEEN_INCLUDED
25 #define TEN_HAS_BEEN_INCLUDED
26 
27 #include <math.h>
28 
29 #include <teem/air.h>
30 #include <teem/biff.h>
31 #include <teem/ell.h>
32 #include <teem/nrrd.h>
33 #include <teem/unrrdu.h>
34 #include <teem/dye.h>
35 #include <teem/gage.h>
36 #include <teem/limn.h>
37 #include <teem/echo.h>
38 
39 #include "tenMacros.h"
40 
41 #if defined(_WIN32) && !defined(__CYGWIN__) && !defined(TEEM_STATIC)
42 #  if defined(TEEM_BUILD) || defined(ten_EXPORTS) || defined(teem_EXPORTS)
43 #    define TEN_EXPORT extern __declspec(dllexport)
44 #  else
45 #    define TEN_EXPORT extern __declspec(dllimport)
46 #  endif
47 #else /* TEEM_STATIC || UNIX */
48 #  define TEN_EXPORT extern
49 #endif
50 
51 #ifdef __cplusplus
52 extern "C" {
53 #endif
54 
55 #define TEN tenBiffKey
56 
57 /*
58 ****** tenAniso* enum
59 **
60 ** the different scalar values that can describe a tensor.  Nearly
61 ** all of these are anisotropy metrics, but with time this has become
62 ** a convenient way to present any scalar invariant (such as trace),
63 ** and even the individual eigenvalues
64 **
65 ** keep in sync:
66 ** aniso.c: _tenAnisoEval_X_f(), _tenAnisoEval_f[]
67 ** aniso.c: _tenAnisoTen_X_d(), _tenAnisoTen_d[]
68 ** aniso.c: tenAnisoCalc_f()
69 ** enumsTen.c: tenAniso
70 */
71 enum {
72   tenAnisoUnknown,  /*  0: nobody knows */
73   tenAniso_Conf,    /*  1: not an anisotropy, but enables some useful hacks */
74   tenAniso_Cl1,     /*  2: Westin's linear (first version) */
75   tenAniso_Cp1,     /*  3: Westin's planar (first version) */
76   tenAniso_Ca1,     /*  4: Westin's linear + planar (first version) */
77   tenAniso_Clpmin1, /*  5: minimum of Cl and Cp (first version) */
78   tenAniso_Cs1,     /*  6: Westin's spherical (first version) */
79   tenAniso_Ct1,     /*  7: gk's anisotropy type (first version) */
80   tenAniso_Cl2,     /*  8: Westin's linear (second version) */
81   tenAniso_Cp2,     /*  9: Westin's planar (second version) */
82   tenAniso_Ca2,     /* 10: Westin's linear + planar (second version) */
83   tenAniso_Clpmin2, /* 11: minimum of Cl and Cp (second version) */
84   tenAniso_Cs2,     /* 12: Westin's spherical (second version) */
85   tenAniso_Ct2,     /* 13: gk's anisotropy type (second version) */
86   tenAniso_RA,      /* 14: Bass+Pier's relative anisotropy */
87   tenAniso_FA,      /* 15: (Bass+Pier's fractional anisotropy)/sqrt(2) */
88   tenAniso_VF,      /* 16: volume fraction = 1-(Bass+Pier's volume ratio) */
89   tenAniso_B,       /* 17: linear term in cubic characteristic polynomial */
90   tenAniso_Q,       /* 18: radius of root circle is 2*sqrt(Q) */
91   tenAniso_R,       /* 19: half of third moment of eigenvalues */
92   tenAniso_S,       /* 20: frobenius norm, squared */
93   tenAniso_Skew,    /* 21: R/sqrt(2*Q^3) */
94   tenAniso_Mode,    /* 22: 3*sqrt(6)*det(dev)/norm(dev) = sqrt(2)*skew */
95   tenAniso_Th,      /* 23: acos(sqrt(2)*skew)/3 */
96   tenAniso_Omega,   /* 24: FA*(1+mode)/2 */
97   tenAniso_Det,     /* 25: plain old determinant */
98   tenAniso_Tr,      /* 26: plain old trace */
99   tenAniso_eval0,   /* 27: largest eigenvalue */
100   tenAniso_eval1,   /* 28: middle eigenvalue */
101   tenAniso_eval2,   /* 29: smallest eigenvalue */
102   tenAnisoLast
103 };
104 #define TEN_ANISO_MAX  29
105 
106 /*
107 ******** tenInterpType* enum
108 **
109 ** different kinds of interpolations paths between tensors
110 */
111 enum {
112   tenInterpTypeUnknown,         /*  0: nobody knows */
113   tenInterpTypeLinear,          /*  1: simple per-coefficient linear */
114   tenInterpTypeLogLinear,       /*  2: linear on logs (Log-Euclidean) */
115   tenInterpTypeAffineInvariant, /*  3: Riemannian approach of many authors */
116   tenInterpTypeWang,            /*  4: affine-invariant of Z Wang & B Vemuri */
117   tenInterpTypeGeoLoxK,         /*  5: geodesic-loxodrome on K_i invariants */
118   tenInterpTypeGeoLoxR,         /*  6: geodesic-loxodrome on R_i invariants */
119   tenInterpTypeLoxK,            /*  7: total loxodrome on K_i invariants */
120   tenInterpTypeLoxR,            /*  8: total loxodrome on R_i invariants */
121   tenInterpTypeQuatGeoLoxK,     /*  9: geodesic-loxodrome on K_i invariants */
122   tenInterpTypeQuatGeoLoxR,     /* 10: geodesic-loxodrome on R_i invariants */
123   tenInterpTypeRThetaPhiLinear, /* 11: linear interpolation of (R,theta,phi),
124                                    oriented by eigensystem of linear avg */
125   tenInterpTypeLast
126 };
127 #define TEN_INTERP_TYPE_MAX        11
128 
129 /*
130 ******** tenGlyphType* enum
131 **
132 ** the different types of glyphs that may be used for tensor viz
133 */
134 enum {
135   tenGlyphTypeUnknown,    /* 0: nobody knows */
136   tenGlyphTypeBox,        /* 1 */
137   tenGlyphTypeSphere,     /* 2 */
138   tenGlyphTypeCylinder,   /* 3 */
139   tenGlyphTypeSuperquad,  /* 4 */
140   tenGlyphTypeBetterquad, /* 5: for T Schultz, GL Kindlmann.
141                              Superquadric Glyphs for Symmetric
142                              Second-Order Tensors. IEEE TVCG
143                              Nov/Dec 2010, 16(6):1595-1604 */
144   tenGlyphTypePolarPlot,  /* 6 */
145   tenGlyphTypeLast
146 };
147 #define TEN_GLYPH_TYPE_MAX   6
148 
149 /*
150 ******** tenGlyphParm struct
151 **
152 ** all input parameters to tenGlyphGen
153 */
154 typedef struct {
155   int verbose;
156 
157   /* glyphs will be shown at samples that have confidence >= confThresh,
158      and anisotropy anisoType >= anisoThresh, and if nmask is non-NULL,
159      then the corresponding mask value must be >= maskThresh.  If
160      onlyPositive, then samples with a non-positive eigenvalue will
161      be skipped, regardless of their purported anisotropy */
162   Nrrd *nmask;
163   int anisoType, onlyPositive;
164   float confThresh, anisoThresh, maskThresh;
165 
166   /* glyphs have shape glyphType and size glyphScale. Superquadrics
167      are tuned by sqdSharp, and things that must polygonalize do so
168      according to facetRes.  Postscript rendering of glyph edges is
169      governed by edgeWidth[] */
170   int glyphType, facetRes;
171   float glyphScale, sqdSharp;
172   float edgeWidth[5];  /* 0: contour, 1: front crease, 2: front non-crease */
173 
174   /* glyphs are colored by eigenvector colEvec with the standard XYZ-RGB
175      colormapping, with maximal saturation colMaxSat (use 0.0 to turn off
176      coloring).  Saturation is modulated by anisotropy colAnisoType, to a
177      degree controlled by colAnisoModulate (if 0, saturation is not at all
178      modulated by anistropy).  Post-saturation, there is a per-channel
179      gamma of colGamma. */
180   int colEvec, colAnisoType;
181   float colMaxSat, colIsoGray, colGamma, colAnisoModulate, ADSP[4];
182 
183   /* if doSlice, a slice of anisotropy sliceAnisoType will be depicted
184      in grayscale as a sheet of grayscale squares, one per sample. As
185      with glyphs, these are thresholded by confThresh and maskThresh
186      (but not anisoThresh).  Things can be lightened up with a
187      sliceGamma > 1, after the slice's gray values have been mapped from
188      [0,1] to [sliceBias,1]. The squares will be at their corresponding
189      sample locations, but offset by sliceOffset */
190   unsigned int sliceAxis;
191   size_t slicePos;
192   int doSlice, sliceAnisoType;
193   float sliceOffset, sliceBias, sliceGamma;
194 } tenGlyphParm;
195 
196 #define TEN_ANISO_DESC \
197   "All the Westin metrics come in two versions.  Currently supported:\n " \
198   "\b\bo \"cl1\", \"cl2\": Westin's linear\n " \
199   "\b\bo \"cp1\", \"cp2\": Westin's planar\n " \
200   "\b\bo \"ca1\", \"ca2\": Westin's linear + planar\n " \
201   "\b\bo \"cs1\", \"cs2\": Westin's spherical (1-ca)\n " \
202   "\b\bo \"ct1\", \"ct2\": GK's anisotropy type (cp/ca)\n " \
203   "\b\bo \"ra\": Basser/Pierpaoli relative anisotropy/sqrt(2)\n " \
204   "\b\bo \"fa\": Basser/Pierpaoli fractional anisotropy\n " \
205   "\b\bo \"vf\": volume fraction = 1-(Basser/Pierpaoli volume ratio)\n " \
206   "\b\bo \"tr\": trace"
207 
208 /*
209 ******** tenGage* enum
210 **
211 ** all the possible queries supported in the tenGage gage kind
212 ** various properties of the quantities below (eigenvalues = v1, v2, v3):
213 ** eigenvalue cubic equation: v^3 + A*v^2 + B*v + C = 0
214 ** Trace = v1 + v2 + v3 = -A
215 ** B = v1*v2 + v1*v3 + v2*v3
216 ** Det = v1*v2*v3 = -C
217 ** S = v1*v1 + v2*v2 + v3*v3
218 ** Q = (S-B)/9 = variance({v1,v2,v3})/2 = (RootRadius/2)^2
219 ** FA = 3*sqrt(Q/S)
220 ** R = (9*A*B - 2*A^3 - 27*C)/54
221      = (5*A*B - 2*A*S - 27*C)/54 = thirdmoment({v1,v2,v3})/2
222 ** P = arccos(R/sqrt(Q)^3)/3 = phase angle of cubic solution
223 **
224 ** NOTE: currently tenGage knows *nothing* about nrrd->measurementFrame!
225 ** You probably want to call tenMeasurementFrameReduce() first.
226 **
227 ** Hey- the problem with adding the RGB eigenvector coloring to the
228 ** capability of tenGage is that because this is visualization, you
229 ** can't easily control whether the measurement frame is applied, if
230 ** known- in that sense the RGB info is uniquely different from the
231 ** other vector and tensor items that can be queried . . . so after a
232 ** brief appearance here the RGB evec coloring was removed.  The
233 ** gagePerVolume->data field that it motivated has rightly remained.
234 **
235 ** !!! Changes to this list need to be propogated to:
236 ** !!! tenGage.c: _tenGageTable[], _tenGageAnswer(),
237 ** !!! enumsTen.c: tenGage airEnum.
238 **
239 */
240 enum {
241   tenGageUnknown,          /*   0: nobody knows */
242 
243   tenGageTensor,           /*   1: "t", the reconstructed tensor: [7] */
244   tenGageConfidence,       /*   2: "c", first of seven tensor values: [1] */
245   tenGageTrace,            /*   3: "tr", trace of tensor: [1] */
246   tenGageNorm,             /*   4: "n", frobenius norm  of tensor: [1] */
247   tenGageB,                /*   5: "b": [1] */
248   tenGageDet,              /*   6: "det", determinant of tensor: [1] */
249   tenGageS,                /*   7: "s", square of frobenius norm: [1] */
250   tenGageQ,                /*   8: "q", (S - B)/9: [1] */
251   tenGageFA,               /*   9: "fa", fractional anisotropy: [1] */
252   tenGageR,                /*  10: "r", 9*A*B - 2*A^3 - 27*C: [1] */
253   tenGageMode,             /*  11: "mode", sqrt(2)*R/sqrt(Q^3): [1] */
254   tenGageTheta,            /*  12: "th", arccos(-mode)/AIR_PI: [1] */
255   tenGageModeWarp,         /*  13: "modew", mode warped for better contrast:
256                                    cos((1-m)*pi/2): [1] */
257   tenGageOmega,            /*  14: "om", fa*(mode+1)/2: [1] */
258 
259   tenGageEval,             /*  15: "eval", all eigenvals of tensor : [3] */
260   tenGageEval0,            /*  16: "eval0", major eigenval of tensor : [1] */
261   tenGageEval1,            /*  17: "eval1", medium eigenval of tensor : [1] */
262   tenGageEval2,            /*  18: "eval2", minor eigenval of tensor : [1] */
263   tenGageEvec,             /*  19: "evec", major eigenvects of tensor: [9] */
264   tenGageEvec0,            /*  20: "evec0", major eigenvect of tensor: [3] */
265   tenGageEvec1,            /*  21: "evec1", medium eigenvect of tensor: [3] */
266   tenGageEvec2,            /*  22: "evec2", minor eigenvect of tensor: [3] */
267 
268   tenGageDelNormK2,        /*  23: "delnk2": normalized gradient tensor of
269                                    K2 = eval variance: [7] */
270   tenGageDelNormK3,        /*  24: "delnk3": normal gradient tensor of
271                                    K3 = R3 = eval skewness: [7] */
272   tenGageDelNormR1,        /*  25: "delnr1": normalized gradient tensor of
273                                    R1 = tensor norm: [7] */
274   tenGageDelNormR2,        /*  26: "delnr2": normalized gradient tensor of
275                                    R2 = FA: [7] */
276   tenGageDelNormPhi1,      /*  27: "delnphi1": normalized rotation tangent
277                                     around principal evector: [7] */
278   tenGageDelNormPhi2,      /*  28: "delnphi2": normalized rotation tangent
279                                     rotation around medium evector: [7] */
280   tenGageDelNormPhi3,      /*  29: "delnphi3": normalized rotation tangent
281                                     rotation around minor evector: [7] */
282 
283   tenGageTensorGrad,       /*  30: "tg", all tensor component gradients,
284                                    starting with confidence gradient: [21] */
285   tenGageTensorGradMag,    /*  31: "tgm", actually a 3-vector of tensor
286                                    gradient norms, one for each axis: [3] */
287   tenGageTensorGradMagMag, /*  32: "tgmm", single scalar magnitude: [1] */
288 
289   tenGageTraceGradVec,     /*  33: "trgv": gradient (vector) of trace: [3] */
290   tenGageTraceGradMag,     /*  34: "trgm": gradient magnitude of trace: [1] */
291   tenGageTraceNormal,      /*  35: "trn": normal of trace: [3] */
292 
293   tenGageNormGradVec,      /*  36: "ngv", gradient (vector) of norm: [3] */
294   tenGageNormGradMag,      /*  37: "ngm", gradient magnitude of norm: [1] */
295   tenGageNormNormal,       /*  38: "nn", normal of norm: [3] */
296 
297   tenGageBGradVec,         /*  39: "bgv", gradient (vector) of B: [3] */
298   tenGageBGradMag,         /*  40: "bgm", gradient magnitude of B: [1] */
299   tenGageBNormal,          /*  41: "bn", normal of B: [3] */
300 
301   tenGageDetGradVec,       /*  42: "detgv", gradient (vector) of Det: [3] */
302   tenGageDetGradMag,       /*  43: "detgm", gradient magnitude of Det: [1] */
303   tenGageDetNormal,        /*  44: "detn", normal of Det: [3] */
304 
305   tenGageSGradVec,         /*  45: "sgv", gradient (vector) of S: [3] */
306   tenGageSGradMag,         /*  46: "sgm", gradient magnitude of S: [1] */
307   tenGageSNormal,          /*  47: "sn", normal of S: [3] */
308 
309   tenGageQGradVec,         /*  48: "qgv", gradient vector of Q: [3] */
310   tenGageQGradMag,         /*  49: "qgm", gradient magnitude of Q: [1] */
311   tenGageQNormal,          /*  50: "qn", normalized gradient of Q: [3] */
312 
313   tenGageFAGradVec,        /*  51: "fagv", gradient vector of FA: [3] */
314   tenGageFAGradMag,        /*  52: "fagm", gradient magnitude of FA: [1] */
315   tenGageFANormal,         /*  53: "fan", normalized gradient of FA: [3] */
316 
317   tenGageRGradVec,         /*  54: "rgv", gradient vector of Q: [3] */
318   tenGageRGradMag,         /*  55: "rgm", gradient magnitude of Q: [1] */
319   tenGageRNormal,          /*  56: "rn", normalized gradient of Q: [3] */
320 
321   tenGageModeGradVec,      /*  57: "mgv", gradient vector of Mode: [3] */
322   tenGageModeGradMag,      /*  58: "mgm", gradient magnitude of Mode: [1] */
323   tenGageModeNormal,       /*  59: "mn", normalized gradient of Mode: [3] */
324 
325   tenGageThetaGradVec,     /*  60: "thgv", gradient vector of Th: [3] */
326   tenGageThetaGradMag,     /*  61: "thgm", gradient magnitude of Th: [1] */
327   tenGageThetaNormal,      /*  62: "thn", normalized gradient of Th: [3] */
328 
329   tenGageOmegaGradVec,     /*  63: "omgv", gradient vector of Omega: [3] */
330   tenGageOmegaGradMag,     /*  64: "omgm", gradient magnitude of Omega: [1] */
331   tenGageOmegaNormal,      /*  65: "omn", normalized gradient of Omega: [3] */
332 
333   tenGageInvarKGrads,      /*  66: "ikgs", projections of tensor gradient onto
334                                   the normalized shape gradients: eval mean,
335                                   variance, skew, in that order: [9]  */
336   tenGageInvarKGradMags,   /*  67: "ikgms", vector magnitude of the spatial
337                                   invariant gradients (above): [3] */
338   tenGageInvarRGrads,      /*  68: "irgs", projections of tensor gradient onto
339                                   the normalized shape gradients: eval mean,
340                                   variance, skew, in that order: [9]  */
341   tenGageInvarRGradMags,   /*  69: "irgms", vector magnitude of the spatial
342                                   invariant gradients (above): [3] */
343   tenGageRotTans,          /*  70: "rts", projections of the tensor grad onto
344                                    *normalized* rotation tangents: [9] */
345   tenGageRotTanMags,       /*  71: "rtms", mags of vectors above: [3] */
346   tenGageEvalGrads,        /*  72: "evgs", projections of tensor gradient onto
347                                    gradients of eigenvalues: [9] */
348   tenGageCl1,              /*  73: same as tenAniso_Cl1, but faster */
349   tenGageCp1,              /*  74: same as tenAniso_Cp1, but faster */
350   tenGageCa1,              /*  75: same as tenAniso_Ca1, but faster */
351   tenGageClpmin1,          /*  76: min(cl1,cp1) */
352   tenGageCl2,              /*  77: same as tenAniso_Cl2, but faster */
353   tenGageCp2,              /*  78: same as tenAniso_Cp2, but faster */
354   tenGageCa2,              /*  79: same as tenAniso_Ca2, but faster */
355   tenGageClpmin2,          /*  80: min(cl2,cp2) */
356 
357   tenGageHessian,          /*  81: "hess", all hessians of tensor
358                                    components: [63] */
359   tenGageTraceHessian,     /*  82: "trhess", hessian(trace): [9] */
360   tenGageTraceHessianEval,  /*  83 */
361   tenGageTraceHessianEval0, /*  84 */
362   tenGageTraceHessianEval1, /*  85 */
363   tenGageTraceHessianEval2, /*  86 */
364   tenGageTraceHessianEvec,  /*  87 */
365   tenGageTraceHessianEvec0, /*  88 */
366   tenGageTraceHessianEvec1, /*  89 */
367   tenGageTraceHessianEvec2, /*  90 */
368   tenGageTraceHessianFrob,  /*  91 */
369   tenGageBHessian,         /*  92: "bhess": [9] */
370   tenGageDetHessian,       /*  93: "dethess": [9] */
371   tenGageSHessian,         /*  94: "shess": [9] */
372   tenGageQHessian,         /*  95: "qhess": [9] */
373 
374   tenGageFAHessian,        /*  96: "fahess": [9] */
375   tenGageFAHessianEval,    /*  97: "fahesseval": [3] */
376   tenGageFAHessianEval0,   /*  98: "fahesseval0": [1] */
377   tenGageFAHessianEval1,   /*  99: "fahesseval1": [1] */
378   tenGageFAHessianEval2,   /* 100: "fahesseval2": [1] */
379   tenGageFAHessianEvec,    /* 101: "fahessevec": [9] */
380   tenGageFAHessianEvec0,   /* 102: "fahessevec0": [3] */
381   tenGageFAHessianEvec1,   /* 103: "fahessevec1": [3] */
382   tenGageFAHessianEvec2,   /* 104: "fahessevec2": [3] */
383   tenGageFAHessianFrob,    /* 105 */
384   tenGageFARidgeSurfaceStrength,  /* 106: "farsurfstrn": [1] */
385   tenGageFAValleySurfaceStrength, /* 107: "favsurfstrn": [1] */
386   tenGageFALaplacian,      /* 108: "falapl": [1] */
387   tenGageFAHessianEvalMode,/* 109: "fahessevalmode": [1] */
388   tenGageFARidgeLineAlignment,    /* 110: "farlinealn": [1] */
389   tenGageFARidgeSurfaceAlignment, /* 111: "farsurfaln": [1] */
390   tenGageFA2ndDD,          /* 112: "fa2d": [1] */
391 
392   tenGageFAGeomTens,       /* 113: "fagten", sym. matx w/ evals {0, K1, K2}
393                                    and evecs {grad, cdir0, cdir1}: [9] */
394   tenGageFAKappa1,         /* 114: "fak1", 1st princ curv: [1] */
395   tenGageFAKappa2,         /* 115: "fak2", 2nd princ curv (k2 <= k1): [1] */
396   tenGageFATotalCurv,      /* 116: "fatc", L2 norm(K1,K2): [1] */
397   tenGageFAShapeIndex,     /* 117: "fasi", Koen.'s shape index, ("S"): [1] */
398   tenGageFAMeanCurv,       /* 118: "famc", mean curvature (K1 + K2)/2: [1] */
399   tenGageFAGaussCurv,      /* 119: "fagc", gaussian curvature K1*K2: [1] */
400   tenGageFACurvDir1,       /* 120: "facdir1", 1st princ curv direction: [3] */
401   tenGageFACurvDir2,       /* 121: "facdir2", 2nd princ curv direction: [3] */
402   tenGageFAFlowlineCurv,   /* 122: "fafc", curv of normal streamline: [1] */
403 
404   tenGageRHessian,         /* 123: "rhess": [9] */
405 
406   tenGageModeHessian,      /* 124: "mhess": [9] */
407   tenGageModeHessianEval,  /* 125: "mhesseval": [3] */
408   tenGageModeHessianEval0, /* 126: "mhesseval0": [1] */
409   tenGageModeHessianEval1, /* 127: "mhesseval1": [1] */
410   tenGageModeHessianEval2, /* 128: "mhesseval2": [1] */
411   tenGageModeHessianEvec,  /* 129: "mhessevec": [9] */
412   tenGageModeHessianEvec0, /* 130: "mhessevec0": [3] */
413   tenGageModeHessianEvec1, /* 131: "mhessevec1": [3] */
414   tenGageModeHessianEvec2, /* 132: "mhessevec2": [3] */
415   tenGageModeHessianFrob,  /* 133 */
416 
417   tenGageOmegaHessian,     /* 134: "omhess": [9] */
418   tenGageOmegaHessianEval, /* 135: "omhesseval": [3] */
419   tenGageOmegaHessianEval0,/* 136: "omhesseval0": [1] */
420   tenGageOmegaHessianEval1,/* 137: "omhesseval1": [1] */
421   tenGageOmegaHessianEval2,/* 138: "omhesseval2": [1] */
422   tenGageOmegaHessianEvec, /* 139: "omhessevec": [9] */
423   tenGageOmegaHessianEvec0,/* 140: "omhessevec0": [3] */
424   tenGageOmegaHessianEvec1,/* 141: "omhessevec1": [3] */
425   tenGageOmegaHessianEvec2,/* 142: "omhessevec2": [3] */
426   tenGageOmegaLaplacian,   /* 143: "omlapl": [1] */
427   tenGageOmega2ndDD,       /* 144: "om2d": [1] */
428   tenGageOmegaHessianContrTenEvec0, /* 145: "omhesscte0": [1] */
429   tenGageOmegaHessianContrTenEvec1, /* 146: "omhesscte1": [1] */
430   tenGageOmegaHessianContrTenEvec2, /* 147: "omhesscte2": [1] */
431 
432   tenGageTraceGradVecDotEvec0,   /* 148: "trgvdotevec0": [1] */
433   tenGageTraceDiffusionAlign,    /* 149: "datr": [1] */
434   tenGageTraceDiffusionFraction, /* 150: "dftr": [1] */
435   tenGageFAGradVecDotEvec0,      /* 151: "fagvdotevec0": [1] */
436   tenGageFADiffusionAlign,       /* 152: "dafa": [1] */
437   tenGageFADiffusionFraction,    /* 153: "dffa": [1] */
438   tenGageOmegaGradVecDotEvec0,   /* 154: "omgvdotevec0": [1] */
439   tenGageOmegaDiffusionAlign,    /* 155: "daom": [1] */
440   tenGageOmegaDiffusionFraction, /* 156: "dfom": [1] */
441   tenGageConfGradVecDotEvec0,    /* 157: "confgvdotevec0": [1] */
442   tenGageConfDiffusionAlign,     /* 158: "daconf": [1] */
443   tenGageConfDiffusionFraction,  /* 159: "dfconf": [1] */
444 
445   tenGageCovariance, /* 160: "cov" 4rth order covariance tensor: [21]
446                         in order of appearance:
447                         0:xxxx  1:xxxy  2:xxxz  3:xxyy  4:xxyz  5:xxzz
448                                 6:xyxy  7:xyxz  8:xyyy  9:xyyz 10:xyzz
449                                        11:xzxz 12:xzyy 13:xzyz 14:xzzz
450                                                15:yyyy 16:yyyz 17:yyzz
451                                                        18:yzyz 19:yzzz
452                                                                20:zzzz */
453   tenGageCovarianceRGRT, /* 161: "covr" covariance tensor expressed in frame
454                             of R invariant gradients and rotation tangents */
455   tenGageCovarianceKGRT, /* 162: "covk" covariance tensor expressed in frame
456                             of K invariant gradients and rotation tangents */
457   tenGageTensorLogEuclidean,     /* 163: "logeuc" log-euclidean interp */
458   tenGageTensorQuatGeoLoxK,      /* 164: "qglk" QGL-K interpolation */
459   tenGageTensorQuatGeoLoxR,      /* 165: "qglr" QGL-R interpolation */
460   tenGageTensorRThetaPhiLinear,  /* 166: "rtpl" RThetaPhiLinear interp */
461 
462   tenGageCl1GradVec,       /* 167: "cl1gv" gradient vector of cl1: [3] */
463   tenGageCl1GradMag,       /* 168: "cl1gm" gradient magnitude of cl1: [1] */
464   tenGageCl1Normal,        /* 169: "cl1gn" normal of cl1: [3] */
465   tenGageCp1GradVec,       /* 170: "cp1gv" gradient vector of cp1: [3] */
466   tenGageCp1GradMag,       /* 171: "cp1gm" gradient magnitude of cp1: [1] */
467   tenGageCp1Normal,        /* 172: "cl1gn" normal of cp1: [3] */
468   tenGageCa1GradVec,       /* 173: "ca1gv" gradient vector of ca1: [3] */
469   tenGageCa1GradMag,       /* 174: "ca1gm" gradient magnitude of ca1: [1] */
470   tenGageCa1Normal,        /* 175: "cl1gn" normal of ca1: [3] */
471   tenGageTensorGradRotE,   /* 176: "tgrote" all tensor component gradients,
472                               starting with confidence gradient.
473                               Rotated such that eigenvalue
474                               derivatives are on the diagonal: [21] */
475   tenGageEvalHessian,    /* 177: "evalhess" Hessian of the eigenvalues: [27] */
476   tenGageCl1Hessian,     /* 178: "cl1hess" Hessian of cl1: [9] */
477   tenGageCl1HessianEval, /* 179: "cl1hesseval" Hessian evals of cl1: [3] */
478   tenGageCl1HessianEval0,/* 180: "cl1hesseval0" First Hess eval of cl1: [1] */
479   tenGageCl1HessianEval1,/* 181: "cl1hesseval1" Second Hess eval of cl1: [1] */
480   tenGageCl1HessianEval2,/* 182: "cl1hesseval2" Third Hess eval of cl1: [1] */
481   tenGageCl1HessianEvec, /* 183: "cl1hessevec" Hessian evecs of cl1: [9] */
482   tenGageCl1HessianEvec0,/* 184: "cl1hessevec0" First Hess evec of cl1: [3] */
483   tenGageCl1HessianEvec1,/* 185: "cl1hessevec1" Second Hess evec of cl1: [3] */
484   tenGageCl1HessianEvec2,/* 186: "cl1hessevec2" Third Hess evec of cl1: [3] */
485   tenGageCp1Hessian,     /* 187: "cp1hess" Hessian of cp1: [9] */
486   tenGageCp1HessianEval, /* 188: "cp1hesseval" Hessian evals of cp1: [3] */
487   tenGageCp1HessianEval0,/* 189: "cp1hesseval0" First Hess eval of cp1: [1] */
488   tenGageCp1HessianEval1,/* 190: "cp1hesseval1" Second Hess eval of cp1: [1] */
489   tenGageCp1HessianEval2,/* 191: "cp1hesseval2" Third Hess eval of cp1: [1] */
490   tenGageCp1HessianEvec, /* 192: "cp1hessevec" Hessian evecs of cp1: [9] */
491   tenGageCp1HessianEvec0,/* 193: "cp1hessevec0" First Hess evec of cp1: [3] */
492   tenGageCp1HessianEvec1,/* 194: "cp1hessevec1" Second Hess evec of cp1: [3] */
493   tenGageCp1HessianEvec2,/* 195: "cp1hessevec2" Third Hess evec of cp1: [3] */
494   tenGageCa1Hessian,     /* 196: "ca1hess" Hessian of ca1: [9] */
495   tenGageCa1HessianEval, /* 197: "ca1hesseval" Hessian evals of ca1: [3] */
496   tenGageCa1HessianEval0,/* 198: "ca1hesseval0" First Hess eval of ca1: [1] */
497   tenGageCa1HessianEval1,/* 199: "ca1hesseval1" Second Hess eval of ca1: [1] */
498   tenGageCa1HessianEval2,/* 200: "ca1hesseval2" Third Hess eval of ca1: [1] */
499   tenGageCa1HessianEvec, /* 201: "ca1hessevec" Hessian evecs of ca1: [9] */
500   tenGageCa1HessianEvec0,/* 202: "ca1hessevec0" First Hess evec of ca1: [3] */
501   tenGageCa1HessianEvec1,/* 203: "ca1hessevec1" Second Hess evec of ca1: [3] */
502   tenGageCa1HessianEvec2,/* 204: "ca1hessevec2" Third Hess evec of ca1: [3] */
503   tenGageFiberCurving,   /* 205: "fibcurv" Savadjiev et al. fiber curving */
504   tenGageFiberDispersion,/* 206: "fibdisp" Savadjiev et al. fiber dispersion */
505   tenGageAniso,          /* 207: "an", all anisos: [TEN_ANISO_MAX+1] */
506   tenGageLast
507 };
508 #define TEN_GAGE_ITEM_MAX     207
509 
510 /*
511 ******** tenDwiGage* enum
512 **
513 ** all things that can be measured in the diffusion weighted images that
514 ** underly diffusion tensor imaging
515 */
516 enum {
517   tenDwiGageUnknown,        /* 0: nobody knows */
518 
519   /*  1: "all", all the measured values, both baseline and diffusion
520       weighted: [N], where N is the number of DWIs */
521   tenDwiGageAll,
522 
523   /*  2: "b0", the non-Dwi image value, either by direct measurement
524       or by estimation: [1]
525       HEY: currently a hack, because it assumes a single known B0 */
526   tenDwiGageB0,
527 
528   /*  3: "jdwi", just the DWIs, no B0: [N-1] (HEY same hack) */
529   tenDwiGageJustDWI,
530 
531   /*  4: "adc", ADCs from the DWIs: [N-1] (HEY same hack) */
532   tenDwiGageADC,
533 
534   /*  5: "mdwi", the average Dwi image value, which is thresholded to
535       create the confidence mask: [1] */
536   tenDwiGageMeanDWIValue,
537 
538   /*  6: "tlls": [7],
539       7: "tllserr": [1],
540       8: "tllserrlog": [1],
541       9: "tllslike": [1],
542      linear least squares fit of tensor value to log(Dwi)s */
543   tenDwiGageTensorLLS,
544   tenDwiGageTensorLLSError,      /* RMS error w/ Dwis */
545   tenDwiGageTensorLLSErrorLog,   /* RMS error w/ log(Dwi)s */
546   tenDwiGageTensorLLSLikelihood,
547 
548   /* 10: "twls": [7],
549      11: "twlserr": [1],
550      12: "twlserrlog": [1],
551      13: "twlslike": [1],
552      weighted least squares fit of tensor value to log(Dwi)s */
553   tenDwiGageTensorWLS,
554   tenDwiGageTensorWLSError,
555   tenDwiGageTensorWLSErrorLog,
556   tenDwiGageTensorWLSLikelihood,
557 
558   /* 14: "tnls": [7],
559      15: "tnlserr": [1],
560      16: "tnlserrlog": [1],
561      17: "tnlslike": [1],
562      non-linear least squares fit of tensor value to Dwis (not log) */
563   tenDwiGageTensorNLS,
564   tenDwiGageTensorNLSError,
565   tenDwiGageTensorNLSErrorLog,
566   tenDwiGageTensorNLSLikelihood,
567 
568   /* 18: "tmle": [7],
569      19: "tmleerr": [1],
570      20: "tmleerrlog": [1],
571      21: "tmlelike": [1],
572      maximum-likelihood fit of tensor value to Dwis */
573   tenDwiGageTensorMLE,
574   tenDwiGageTensorMLEError,
575   tenDwiGageTensorMLEErrorLog,
576   tenDwiGageTensorMLELikelihood,
577 
578   /* 22: "t": [7],
579      23: "terr": [1],
580      24: "terrlog": [1],
581      25: "tlike": [1],
582      one of the above tensors and its errors, depending on settings */
583   tenDwiGageTensor,
584   tenDwiGageTensorError,
585   tenDwiGageTensorErrorLog,
586   tenDwiGageTensorLikelihood,
587 
588   /* 26: "c", first of seven tensor values: [1] */
589   tenDwiGageConfidence,
590 
591   /* 27: "fa", FA computed from the single tensor: [1] */
592   tenDwiGageFA,
593 
594   /* 28: "adwie", all errors between measured and predicted DWIs
595      [N-1] (HEY same hack) */
596   tenDwiGageTensorAllDWIError,
597 
598   /* 29: "2qserr": [1]
599      30: "2qs", two tensor fitting by q-ball segmentation: [14]
600      31: "2qsnerr": [15] */
601   tenDwiGage2TensorQSeg,
602   tenDwiGage2TensorQSegError,
603   tenDwiGage2TensorQSegAndError,
604 
605   /* 32: "2pelederr": [1]
606      33: "2peled", two tensor fitting by q-ball segmentation: [14]
607      34: "2pelednerr": [15] */
608   tenDwiGage2TensorPeled,
609   tenDwiGage2TensorPeledError,
610   tenDwiGage2TensorPeledAndError,
611 
612   /* 35: "2peledlminfo", levmar output info vector: [9]
613      note: length 9 being correct is checked in _tenDwiGagePvlDataNew() */
614   tenDwiGage2TensorPeledLevmarInfo,
615 
616   tenDwiGageLast
617 };
618 #define TEN_DWI_GAGE_ITEM_MAX 35
619 
620 /*
621 ******** tenEstimate1Method* enum
622 **
623 ** the different ways of doing single tensor estimation
624 */
625 enum {
626   tenEstimate1MethodUnknown,   /* 0 */
627   tenEstimate1MethodLLS,       /* 1 */
628   tenEstimate1MethodWLS,       /* 2 */
629   tenEstimate1MethodNLS,       /* 3 */
630   tenEstimate1MethodMLE,       /* 4 */
631   tenEstimate1MethodLast
632 };
633 #define TEN_ESTIMATE_1_METHOD_MAX 4
634 
635 /*
636 ******** tenEstimate2Method* enum
637 **
638 ** the different ways of doing two-tensor estimation
639 */
640 enum {
641   tenEstimate2MethodUnknown,   /* 0 */
642   tenEstimate2MethodQSegLLS,   /* 1 */
643   tenEstimate2MethodPeled,     /* 2 */
644   tenEstimate2MethodLast
645 };
646 #define TEN_ESTIMATE_2_METHOD_MAX 2
647 
648 /*
649 ******** tenEvecRGBParm struct
650 **
651 ** dumb little bag for the parameters relating to how to do the
652 ** eigenvector -> RGB mapping, since its needed by various things in
653 ** various contexts.  Note that you may need two of these, one for
654 ** doing rgb(evec0) (the linear part) and one for doing rgb(evec2)
655 ** (the planar part).  This used to have "aniso0" and "aniso2", but
656 ** the associated methods were clumsy and redundant.
657 */
658 typedef struct {
659   unsigned int which; /* when the eigenvector hasn't already been computed,
660                          which eigenvector to map:
661                          0 for linear, 2 or planar, 1 for orthotropic */
662   int aniso;          /* which anisotropy metric modulates saturation */
663   double confThresh,  /* confidence threshold */
664     anisoGamma,       /* gamma on aniso, pre-mapping */
665     gamma,            /* per RGB component gamma, post-mapping */
666     bgGray,           /* gray-value for low confidence samples */
667     isoGray,          /* gray-value for isotropic samples */
668     maxSat;           /* maximum saturation */
669   int typeOut,        /* when output type is flexible, and if this is
670                          nrrdTypeUChar or nrrdTypeUShort, then output will
671                          be quantized to those types (range [0,255] and
672                          [0,65535] respectively); otherwise values are
673                          copied directly to output */
674     genAlpha;         /* when output value set is flexible, create RGBA
675                          values instead of just RGB */
676 } tenEvecRGBParm;
677 
678 /*
679 ******** tenFiberType* enum
680 **
681 ** the different kinds of fiber tractography that we do
682 */
683 enum {
684   tenFiberTypeUnknown,    /* 0: nobody knows */
685   tenFiberTypeEvec0,      /* 1: standard following of principal eigenvector */
686   tenFiberTypeEvec1,      /* 2: following medium eigenvector */
687   tenFiberTypeEvec2,      /* 3: following minor eigenvector */
688   tenFiberTypeTensorLine, /* 4: Weinstein-Kindlmann tensorlines */
689   tenFiberTypePureLine,   /* 5: "pure" tensorlines- multiplication only */
690   tenFiberTypeZhukov,     /* 6: Zhukov's oriented tensor reconstruction */
691   tenFiberTypeLast
692 };
693 #define TEN_FIBER_TYPE_MAX   6
694 
695 /*
696 ******** tenDwiFiberType* enum
697 **
698 ** how tractography is done in DWI volumes.  This is orthogonal to
699 ** how single- or two-tensor estimation is done; it describes what we
700 ** do with the model(s) once estimated
701 */
702 enum {
703   tenDwiFiberTypeUnknown,      /* 0: nobody knows */
704   tenDwiFiberType1Evec0,       /* 1: like old-fashioned tractography */
705   tenDwiFiberType2Evec0,       /* 2: only using 2-tensor fits */
706   tenDwiFiberType12BlendEvec0, /* 3: blend between 1- and 2-ten evec0 methods,
707                                   based on something else */
708   tenDwiFiberTypeLast
709 };
710 #define TEN_DWI_FIBER_TYPE_MAX    3
711 
712 /*
713 ******** tenFiberIntg* enum
714 **
715 ** the different integration styles supported.  Obviously, this is more
716 ** general purpose than fiber tracking, so this will be moved (elsewhere
717 ** in Teem) as needed
718 */
719 enum {
720   tenFiberIntgUnknown,   /* 0: nobody knows */
721   tenFiberIntgEuler,     /* 1: dumb but fast */
722   tenFiberIntgMidpoint,  /* 2: 2nd order Runge-Kutta */
723   tenFiberIntgRK4,       /* 3: 4rth order Runge-Kutta */
724   tenFiberIntgLast
725 };
726 #define TEN_FIBER_INTG_MAX  3
727 
728 /*
729 ******** tenFiberStop* enum
730 **
731 ** the different reasons why fibers stop going (as stored in
732 ** tenFiberSingle->whyStop[]), or never got started
733 ** (tenFiberSingle->whyNowhere), or never went far enough (also
734 ** tenFiberSingle->whyNowhere).
735 **
736 ** The addition of tenFiberStopMinLength and tenFiberStopMinNumSteps
737 ** really stretch the meaningfulness of "tenFiberStop", but its the
738 ** only logical place for such constraints to go.
739 ** NOTE: tenFiberStopMinLength and tenFiberStopMinNumSteps only make
740 ** sense as a value for whyNowhere, not whyStop, despite the name.
741 */
742 enum {
743   tenFiberStopUnknown,     /*  0: nobody knows,
744                                   or, for tfx->whyNowhere: no, actually,
745                                   we *did* get somewhere with this fiber */
746   tenFiberStopAniso,       /*  1: specified aniso got below specified level */
747   tenFiberStopLength,      /*  2: fiber length in world space got too long */
748   tenFiberStopNumSteps,    /*  3: took too many steps along fiber */
749   tenFiberStopConfidence,  /*  4: tensor "confidence" value went too low */
750   tenFiberStopRadius,      /*  5: radius of curvature got too small */
751   tenFiberStopBounds,      /*  6: fiber position stepped outside volume */
752   tenFiberStopFraction,    /*  7: during multi-tensor tracking, fractional
753                                   constituency of the tracked tensor got
754                                   too small */
755   tenFiberStopStub,        /*  8: treat single vertex fibers as non-starters */
756   tenFiberStopMinLength,   /*  9: fibers with total (both halves) small length
757                                   are discarded */
758   tenFiberStopMinNumSteps, /* 10: fibers with total (both halves) small # steps
759                                   are discarded (more general-purpose than
760                                   tenFiberStopStub) */
761   tenFiberStopLast
762 };
763 #define TEN_FIBER_STOP_MAX    10
764 
765 /*
766 ******** #define TEN_FIBER_NUM_STEPS_MAX
767 **
768 ** whatever the stop criteria are for fiber tracing, no fiber half can
769 ** have more points than this- a useful sanity check against fibers
770 ** done amok.
771 */
772 #define TEN_FIBER_NUM_STEPS_MAX 10240
773 
774 enum {
775   tenFiberParmUnknown,         /* 0: nobody knows */
776   tenFiberParmStepSize,        /* 1: base step size */
777   tenFiberParmUseIndexSpace,   /* 2: non-zero iff output of fiber should be
778                                   seeded in and output in index space,
779                                   instead of default world */
780   tenFiberParmWPunct,          /* 3: tensor-line parameter */
781   tenFiberParmVerbose,         /* 4: verbosity */
782   tenFiberParmLast
783 };
784 #define TEN_FIBER_PARM_MAX        4
785 
786 enum {
787   tenTripleTypeUnknown,    /* 0: nobody knows */
788   tenTripleTypeEigenvalue, /* 1: eigenvalues sorted in descending order */
789   tenTripleTypeMoment,     /* 2: (mu1,mu2,mu3) */
790   tenTripleTypeXYZ,        /* 3: eval rotation, after Bahn'99 JMR:141(68-77) */
791   tenTripleTypeRThetaZ,    /* 4: cylindrical coords of rotated evals */
792   tenTripleTypeRThetaPhi,  /* 5: spherical coords of rotated evals */
793   tenTripleTypeJ,          /* 6: (J1,J2,J3) principal invariants */
794   tenTripleTypeK,          /* 7: (K1,K2,K3) cylindrical invariants */
795   tenTripleTypeR,          /* 8: (R1,R2,R3) spherical invariants */
796   tenTripleTypeWheelParm,  /* 9: eigenvalue wheel (center,radius,angle) */
797   tenTripleTypeLast
798 };
799 #define TEN_TRIPLE_TYPE_MAX   9
800 
801 /*
802 ******** tenFiberContext
803 **
804 ** catch-all for input, state, and output of fiber tracing.  Luckily, like
805 ** in a gageContext, NOTHING in here is set directly by the user; everything
806 ** should be through the tenFiber* calls
807 */
808 typedef struct {
809   /* ---- input -------- */
810   const Nrrd *nin;      /* the tensor OR DWI volume being analyzed */
811   NrrdKernelSpec *ksp;  /* reconstruction kernel for tensors or DWIs */
812   int useDwi,           /* we're working in a DWI, not a tensor, volume */
813     fiberType,          /* from tenFiberType* OR tenDwiFiberType* enum */
814     fiberProbeItem,     /* item to probe along fiber and possibly save in
815                            tenFiberSingle->nval */
816     intg,               /* from tenFiberIntg* enum */
817     anisoStopType,      /* which aniso we do a threshold on */
818     anisoSpeedType,     /* base step size is function of this anisotropy */
819     stop,               /* BITFLAG for different reasons to stop a fiber */
820     useIndexSpace,      /* output in index space, not world space */
821     verbose;            /* blah blah blah */
822   double anisoThresh,   /* anisotropy threshold */
823     anisoSpeedFunc[3];  /* parameters of mapping aniso to speed */
824   unsigned int maxNumSteps, /* max # steps allowed on one fiber *half* */
825     minNumSteps;        /* min signficiant # steps on *whole* fiber */
826   double stepSize,      /* step size in world space */
827     maxHalfLen,         /* longest propagation (forward or backward) allowed
828                            from midpoint */
829     minWholeLen,        /* minimum significant length of whole fiber */
830     confThresh,         /* confidence threshold */
831     minRadius,          /* minimum radius of curvature of path */
832     minFraction;        /* minimum fractional constituency in multi-tensor */
833   double wPunct;        /* knob for tensor lines */
834   unsigned int ten2Which;  /* which path to follow in 2-tensor tracking */
835   /* ---- internal ----- */
836   gageQuery query;      /* query we'll send to gageQuerySet */
837   int halfIdx,          /* current fiber half being computed (0 or 1) */
838     mframeUse;          /* need to use mframe[] and mframeT[] */
839   double mframe[9],     /* measurement frame in normal matrix lay-out */
840     mframeT[9],         /* transpose of mframe[] */
841     wPos[3],            /* current world space location */
842     wDir[3],            /* difference between this and last world space pos */
843     lastDir[3],         /* previous value of wDir */
844     seedEvec[3];        /* principal eigenvector first found at seed point */
845   int lastDirSet,       /* lastDir[] is usefully set */
846     lastTenSet;         /* lastTen[] is usefully set */
847   unsigned int ten2Use; /* which of the 2-tensors was last used */
848   gageContext *gtx;     /* wrapped around pvl */
849   gagePerVolume *pvl;   /* wrapped around dtvol */
850 
851   const double *gageTen,    /* gageAnswerPointer(pvl, tenGageTensor) */
852     *gageEval,              /* gageAnswerPointer(pvl, tenGageEval) */
853     *gageEvec,              /* gageAnswerPointer(pvl, tenGageEvec) */
854     *gageAnisoStop,         /* gageAnswerPointer(pvl, tenGage<anisoStop>) */
855     *gageAnisoSpeed,        /* gageAnswerPointer(pvl, tenGage<anisoSpeed>) */
856     *gageTen2;              /* gageAnswerPointer(pvl, tenDwiGage..2Tensor..) */
857   double ten2AnisoStop;
858   double fiberTen[7], fiberEval[3], fiberEvec[9],
859     fiberAnisoStop, fiberAnisoSpeed;
860   double radius;        /* current radius of curvature */
861   /* ---- output ------- */
862   double halfLen[2];    /* length of each fiber half in world space */
863   unsigned int numSteps[2]; /* how many samples are used for each fiber half */
864   int whyStop[2],       /* why backward/forward (0/1) tracing stopped
865                            (from tenFiberStop* enum) */
866     whyNowhere;         /* why fiber never got started (from tenFiberStop*) */
867 
868 } tenFiberContext;
869 
870 /*
871 ******** tenFiberSingle
872 **
873 ** experimental struct for holding results from a single tracing
874 */
875 typedef struct {
876   /* ------- available for recording for reference, not used by ten */
877   double seedPos[3];    /* where was the seed point */
878   unsigned int dirIdx;  /* which direction at seedpoint to follow */
879   unsigned int dirNum;  /* how many directions at seedpnt could be followed */
880   /* ------- output ------- */
881   Nrrd *nvert;          /* locations of tract vertices */
882   double halfLen[2];    /* (same as in tenFiberContext) */
883   unsigned int seedIdx, /* which index in nvert is for seedpoint */
884     stepNum[2];         /* (same as in tenFiberContext) */
885   int whyStop[2],       /* (same as in tenFiberContext) */
886     whyNowhere;         /* (same as in tenFiberContext) */
887   Nrrd *nval;           /* results of probing at vertices */
888   double measr[NRRD_MEASURE_MAX+1];  /* a controlled mess */
889 } tenFiberSingle;
890 
891 /*
892 ******** tenFiberMulti
893 **
894 ** container for multiple fibers
895 */
896 typedef struct {
897   tenFiberSingle *fiber;
898   unsigned int fiberNum;
899   airArray *fiberArr;
900 } tenFiberMulti;
901 
902 /*
903 ******** struct tenEmBimodalParm
904 **
905 ** input and output parameters for tenEMBimodal (for fitting two
906 ** gaussians to a histogram).  Currently, all fields are directly
907 ** set/read; no API help here.
908 **
909 ** "fraction" means prior probability
910 **
911 ** In the output, material #1 is the one with the lower mean
912 */
913 typedef struct {
914   /* ----- input -------- */
915   double minProb,        /* threshold for negligible posterior prob. values */
916     minProb2,            /* minProb for 2nd stage fitting */
917     minDelta,            /* convergence test for maximization */
918     minFraction,         /* smallest fraction (in 0.0 to 1.0) that material
919                             1 or 2 can legitimately have */
920     minConfidence,       /* smallest confidence value that the model fitting
921                             is allowed to have */
922     twoStage,            /* wacky two-stage fitting */
923     verbose;             /* output messages and/or progress images */
924   unsigned int maxIteration; /* cap on # of non-convergent iters allowed */
925   /* ----- internal ----- */
926   double *histo,         /* double version of histogram */
927     *pp1, *pp2,          /* pre-computed posterior probabilities for the
928                             current iteration */
929     vmin, vmax,          /* value range represented by histogram. This is
930                             saved from given histogram, and used to inform
931                             final output values, but it is not used for
932                             any intermediate histogram calculations, all of
933                             which are done entirely in index space */
934     delta;               /* some measure of model change between iters */
935   int N,                 /* number of bins in histogram */
936     stage;               /* current stage (1 or 2) */
937   unsigned int iteration;  /* current iteration */
938   /* ----- output ------- */
939   double mean1, stdv1,   /* material 1 mean and  standard dev */
940     mean2, stdv2,        /* same for material 2 */
941     fraction1,           /* fraction of material 1 (== 1 - fraction2) */
942     confidence,          /* (mean2 - mean1)/(stdv1 + stdv2) */
943     threshold;           /* minimum-error threshold */
944 } tenEMBimodalParm;
945 
946 /*
947 ******** struct tenGradientParm
948 **
949 ** all parameters for repulsion-based generation of gradient directions
950 **
951 ** the old physics-based point-repulsion code (RK2 integration of
952 ** equations of motion, with drag force) is gone; this is only the
953 ** fast gradient descent with some strategies for adaptive step size
954 */
955 typedef struct {
956   /* ----------------------- INPUT */
957   double initStep,        /* initial step size for gradient descent */
958     jitter,               /* amount by which distribution is jittered
959                              when starting with a given input set, as
960                              a fraction of the ideal edge length */
961     minVelocity,          /* minimum mean gradient velocity that signifies
962                              end of first distribution phase */
963     minPotentialChange,   /* minimum change in potential that signifies
964                              end of first distribution phase */
965     minMean,              /* mean gradient length that signifies end of
966                              secondary balancing phase */
967     minMeanImprovement;   /* magnitude of improvement (reduction) of mean
968                              gradient length that signifies end of
969                              secondary balancing phase */
970   int single,             /* distribute single points, instead of
971                              anti-podal pairs of points */
972     insertZeroVec,        /* when computing output in
973                              tenGradientDistribute (and only there, though
974                              this is called by tenGradientGenerate), insert
975                              the zero vector at the beginning of the list,
976                              corresponding to a non-DWI B0 image */
977     verbose;              /* verbosity level; 0 turns off everything */
978   unsigned int snap,      /* interval of interations at which to save
979                              snapshats of distribution */
980     report,               /* interval of interations at which to report
981                              on progress */
982     expo,                 /* the exponent N that defines the potential
983                              energy profile 1/r^N (coulomb: N=1) */
984     seed,                 /* seed value for random number generator */
985     maxEdgeShrink,        /* max number of times we try to compute
986                              an update with smaller edge normalization */
987     minIteration,         /* run for at least this many iterations,
988                              which can be useful for high exponents,
989                              for which potential measurements can
990                              easily go to infinity */
991     maxIteration;         /* bail if we haven't converged by this number
992                              of iterations */
993   double expo_d;          /* floating point exponent.  If expo is zero,
994                              this is the value that matters */
995   /* ----------------------- INTERNAL */
996   double step,            /* actual current step size (adjusted during
997                              the algorithm depending on progress) */
998     nudge;                /* how to increase realDT with each iteration */
999   /* ----------------------- OUTPUT */
1000   unsigned int itersUsed; /* total number of iterations */
1001   double potential,       /* potential, without edge normalization */
1002     potentialNorm,        /* potential, with edge normalization */
1003     angle,                /* minimum angle */
1004     edge;                 /* minimum edge length */
1005 } tenGradientParm;
1006 
1007 /*
1008 ******** struct tenEstimateContext
1009 **
1010 ** for handling estimation of diffusion models
1011 */
1012 typedef struct {
1013   /* input ----------- */
1014   double bValue,           /* scalar b value */
1015     valueMin,              /* smallest sensible for input Dwi value,
1016                               must be > 0.0 (for taking log) */
1017     sigma,                 /* noise parameter */
1018     dwiConfThresh,         /* mean Dwi threshold for confidence mask */
1019     dwiConfSoft;           /* softness in confidence mask */
1020                            /* NOTE: for both _ngrad and _nbmat:
1021                               1) only one can be non-NULL, and axis[1].size
1022                               is the total # values, both Dwi and non-Dwi
1023                               2) NO additional re-normalization is done on
1024                               the grads/bmats, UNLIKE the normalization
1025                               performed by tenDWMRIKeyValueParse(). */
1026   const Nrrd *_ngrad,      /* caller's 3-by-allNum gradient list */
1027     *_nbmat;               /* caller's 6-by-allNum B-matrix list,
1028                               off-diagonals are *NOT* pre-multiplied by 2 */
1029   unsigned int *skipList;  /* list of value indices that we are to skip */
1030   airArray *skipListArr;   /* airArray around skipList */
1031   const float *all_f;      /* caller's list of all values (length allNum) */
1032   const double *all_d;     /* caller's list of all values (length allNum) */
1033   int simulate,            /* if non-zero, we're being used for simulation,
1034                               not estimation: be tolerant of unset parms */
1035     estimate1Method,       /* what kind of single-tensor estimation to do */
1036     estimateB0,            /* if non-zero, B0 should be estimated along with
1037                               rest of model. Otherwise, B0 is found by simply
1038                               taking average of non-Dwi images */
1039     recordTime,            /* if non-zero, record estimation time */
1040     recordErrorDwi,
1041     recordErrorLogDwi,
1042     recordLikelihoodDwi,
1043     verbose,               /* blah blah blah */
1044     negEvalShift,          /* if non-zero, shift eigenvalues upwards so that
1045                               smallest one is non-negative */
1046     progress;              /* progress indication for volume processing */
1047   unsigned int WLSIterNum; /* number of iterations for WLS */
1048   /* internal -------- */
1049   /* a "dwi" in here is basically any value (diffusion-weighted or not)
1050      that varies as a function of the model parameters being estimated */
1051   int flag[128];           /* flags for state management */
1052   unsigned int allNum,     /* total number of images (Dwi and non-Dwi) */
1053     dwiNum;                /* number of Dwis */
1054   Nrrd *nbmat,             /* B-matrices (dwiNum of them) for the Dwis, with
1055                               off-diagonals (*YES*) pre-multiplied by 2,
1056                               and with a 7th column of -1.0 if estimateB0 */
1057     *nwght,                /* dwiNum x dwiNum matrix of weights */
1058     *nemat;                /* pseudo-inverse of nbmat */
1059   double knownB0,          /* B0 known from DWI, only if !estimateB0 */
1060     *all,                  /* (copy of) vector of input values,
1061                               allocated for allNum */
1062     *bnorm,                /* frob norm of B-matrix, allocated for allNum */
1063     *allTmp, *dwiTmp,      /* for storing intermediate values,
1064                               allocated for allNum and dwiNum respectively */
1065     *dwi;                  /* the Dwi values, allocated for dwiNum */
1066   unsigned char *skipLut;  /* skipLut[i] non-zero if we should ignore all[i],
1067                               allocated for allNum */
1068   /* output ---------- */
1069   double estimatedB0,      /* estimated non-Dwi value, only if estimateB0 */
1070     ten[7],                /* the estimated single tensor */
1071     conf,                  /* the "confidence" mask value (i.e. ten[0]) */
1072     mdwi,                  /* mean Dwi value (used for conf mask calc) */
1073     time,                  /* time required for estimation */
1074     errorDwi,              /* error in Dwi of estimate */
1075     errorLogDwi,           /* error in log(Dwi) of estimate */
1076     likelihoodDwi;         /* the maximized likelihood */
1077 } tenEstimateContext;
1078 
1079 /*
1080 ******** struct tenDwiGageKindData
1081 **
1082 ** the kind data has static info that is set by the user
1083 ** and then not altered-- or else it wouldn't really be per-kind,
1084 ** and it wouldn't be thread-safe
1085 */
1086 typedef struct {
1087   Nrrd *ngrad, *nbmat;       /* owned by us */
1088   double thresh, soft, bval, valueMin;
1089   int est1Method, est2Method;
1090   unsigned int randSeed;
1091 } tenDwiGageKindData;
1092 
1093 /*
1094 ******** struct tenDwiGagePvlData
1095 **
1096 ** the pvl data is the dynamic stuff (buffers, mostly)
1097 ** that is potentially modified per query.
1098 ** being part of the pvl, such modifications are thread-safe
1099 **
1100 ** the reason for having two tenEstimateContexts is that within
1101 ** one volume you will sometimes want to measure both one and
1102 ** two tensors, and it would be crazy to incur the overhead of
1103 ** switching between the two *per-query*.
1104 */
1105 typedef struct {
1106   tenEstimateContext *tec1, /* for estimating single tensors */
1107     *tec2;                  /* for estimating two-tensors */
1108   double *vbuf;
1109   unsigned int *wght;
1110   double *qvals;
1111   double *qpoints;
1112   double *dists;
1113   double *weights;
1114   Nrrd *nten1EigenGrads;
1115   airRandMTState *randState;
1116   unsigned int randSeed;
1117   double ten1[7], ten1Evec[9], ten1Eval[3];
1118   int levmarUseFastExp;
1119   unsigned int levmarMaxIter;
1120   double levmarTau, levmarEps1, levmarEps2, levmarEps3,
1121     levmarDelta, levmarMinCp;
1122   double levmarInfo[9]; /* output */
1123 } tenDwiGagePvlData;
1124 
1125 typedef struct {
1126   /* input ------------- */
1127   int verbose;
1128   double convStep, minNorm, convEps, wghtSumEps;
1129   int enableRecurse;
1130   unsigned int maxIter, numSteps;
1131   int lengthFancy;
1132   /* internal ------------ */
1133   unsigned int allocLen;
1134   double *eval, *evec, *rtIn, *rtLog, *qIn, *qBuff, *qInter;
1135   /* output ------------ */
1136   unsigned int numIter;
1137   double convFinal;
1138   double lengthShape, lengthOrient;
1139 } tenInterpParm;
1140 
1141 /* the idea is that there should be a uniform way of describing
1142    DWI experiments, and as a result of calling one of the Set()
1143    functions below, the information is set in tenExperSpec so that
1144    combined with a diffusion model spec, DWIs can be simulated */
1145 typedef struct {
1146   int set;               /* has been set */
1147   unsigned imgNum;       /* total number of images, dwi or not */
1148   double *bval,          /* all b values: imgNum doubles.
1149                             NOTE: if dwi[i] is to be considered a "b0" value
1150                             (for the various functions that have a knownB0 or
1151                             similarly named flag), bval[i] must be 0.0;
1152                             regardless of the gradient length */
1153   /* *wght,                 all weights, but not yet used */
1154     *grad;               /* all gradients: 3 x imgNum doubles */
1155 } tenExperSpec;
1156 
1157 /* description of *one* parameter in the parameter vector that
1158    defines one instance of a given model */
1159 typedef struct {
1160   char name[AIR_STRLEN_SMALL]; /* name */
1161   double min, max;             /* bounds */
1162   int cyclic;                  /* non-zero if effectively min == max */
1163   int vec3;                    /* non-zero if this is a coefficient
1164                                   of a UNIT-LENGTH 3-vector */
1165   unsigned int vecIdx;         /* *if* this param is part of vector,
1166                                   the index into it, i.e. 0, 1, or 2 */
1167 } tenModelParmDesc;
1168 
1169 #define TEN_MODEL_B0_MAX 65500    /* HEY: fairly arbitrary, but is set to be
1170                                      a little below the max storable value
1171                                      of the unsigned short in which DWI
1172                                      values might be saved */
1173 #define TEN_MODEL_DIFF_MAX 0.006  /* in units of mm^2/sec; diffusivity of
1174                                      water is about 0.003 mm^2/sec */
1175 #define TEN_MODEL_PARM_GRAD_EPS 0.000005 /* for gradient calculations */
1176 
1177 /*
1178 ******** struct tenModel
1179 **
1180 ** Container for *information* about how DWI values arise (i.e. model
1181 ** parameters), with functions that help in fitting and evaluating
1182 ** models.  A tenModel does not contain any vector of model parameter
1183 ** values, it is only a description of those values (especially its
1184 ** tenModelParmDesc vector), and utilities for computing with those
1185 ** parameter values.
1186 **
1187 ** NOTE: the current convention is to *always* have the non-DW T2
1188 ** image value ("B0") be the first value in the parameter vector that
1189 ** a tenModel is used to do describe.  The B0 value will be found either
1190 ** by trivial means (i.e copied from the image data), or by a different
1191 ** method than the rest of the model parameters, but it is stored along with
1192 ** the parameters because
1193 ** (1) its very handy to have in one place all the information that, when
1194 ** combined with the tenExperSpec, gives you a DWI value, and
1195 ** (2) sometimes B0 does have to be estimated at the same time as the rest
1196 ** of the model, so we thereby avoid doubling the number of models
1197 ** to be able to support this.
1198 ** On the other hand, the B0 value may not be stored along with the rest of
1199 ** the parm vec in the case of saving out whole nrrd of parm vecs.
1200 ** The "parmNum" field below therefore always includes the B0 value
1201 */
1202 typedef struct tenModel_t {
1203   char name[AIR_STRLEN_SMALL];
1204   unsigned int parmNum;
1205   const tenModelParmDesc *parmDesc;
1206   /* noise free simulation */
1207   void (*simulate)(double *dwiSim, const double *parm,
1208                    const tenExperSpec *espec);
1209   /* parameter vector operations */
1210   char *(*sprint)(char str[AIR_STRLEN_MED], const double *parm);
1211   double *(*alloc)(void);
1212   void (*rand)(double *parm, airRandMTState *rng, int knownB0);
1213   void (*step)(double *parm1, const double scl,
1214                const double *grad, const double *parm0);
1215   double (*dist)(const double *parmA, const double *parmB);
1216   void (*copy)(double *parmA, const double *parmB);
1217   int (*convert)(double *parmDst, const double *parmSrc,
1218                  const struct tenModel_t *modelSrc);
1219 
1220   /* "sqe" == squared error in DWI values */
1221   double (*sqe)(const double *parm, const tenExperSpec *espec,
1222                 double *dwiBuff, const double *dwiMeas, int knownB0);
1223   void (*sqeGrad)(double *grad, const double *parm,
1224                   const tenExperSpec *espec,
1225                   double *dwiBuff, const double *dwiMeas,
1226                   int knownB0);
1227   double (*sqeFit)(double *parm, double *convFrac, unsigned int *itersTaken,
1228                    const tenExperSpec *espec,
1229                    double *dwiBuff, const double *dwiMeas,
1230                    const double *parmInit, int knownB0,
1231                    unsigned int minIter, unsigned int maxIter,
1232                    double convEps, int verbose);
1233   /* "nll" == negative log likelihood */
1234   double (*nll)(const double *parm, const tenExperSpec *espec,
1235                 double *dwiBuff, const double *dwiMeas,
1236                 int rician, double sigma, int knownB0);
1237   void (*nllGrad)(double *grad, const double *parm,
1238                   const tenExperSpec *espec,
1239                   double *dwiBuff, const double *dwiMeas,
1240                   int rician, double sigma);
1241   double (*nllFit)(double *parm, const tenExperSpec *espec,
1242                    const double *dwiMeas, const double *parmInit,
1243                    int rician, double sigma, int knownB0);
1244 } tenModel;
1245 
1246 /* defaultsTen.c */
1247 TEN_EXPORT const int tenPresent;
1248 TEN_EXPORT const char *tenBiffKey;
1249 TEN_EXPORT const char tenDefFiberKernel[];
1250 TEN_EXPORT double tenDefFiberStepSize;
1251 TEN_EXPORT int tenDefFiberUseIndexSpace;
1252 TEN_EXPORT int tenDefFiberMaxNumSteps;
1253 TEN_EXPORT double tenDefFiberMaxHalfLen;
1254 TEN_EXPORT int tenDefFiberAnisoStopType;
1255 TEN_EXPORT double tenDefFiberAnisoThresh;
1256 TEN_EXPORT int tenDefFiberIntg;
1257 TEN_EXPORT double tenDefFiberWPunct;
1258 
1259 /* triple.c */
1260 TEN_EXPORT void tenTripleConvertSingle_d(double dst[3],
1261                                          int dstType,
1262                                          const double src[3],
1263                                          const int srcType);
1264 TEN_EXPORT void tenTripleConvertSingle_f(float dst[3],
1265                                          int dstType,
1266                                          const float src[3],
1267                                          const int srcType);
1268 TEN_EXPORT void tenTripleCalcSingle_d(double dst[3],
1269                                       int ttype, double ten[7]);
1270 TEN_EXPORT void tenTripleCalcSingle_f(float dst[3],
1271                                       int ttype, float ten[7]);
1272 TEN_EXPORT int tenTripleCalc(Nrrd *nout, int ttype, const Nrrd *nten);
1273 TEN_EXPORT int tenTripleConvert(Nrrd *nout, int dstType,
1274                                 const Nrrd *nin, int srcType);
1275 
1276 /* grads.c */
1277 TEN_EXPORT tenGradientParm *tenGradientParmNew(void);
1278 TEN_EXPORT tenGradientParm *tenGradientParmNix(tenGradientParm *tgparm);
1279 TEN_EXPORT int tenGradientCheck(const Nrrd *ngrad, int type,
1280                                 unsigned int minnum);
1281 TEN_EXPORT int tenGradientRandom(Nrrd *ngrad, unsigned int num,
1282                                  unsigned int seed);
1283 TEN_EXPORT double tenGradientIdealEdge(unsigned int N, int single);
1284 TEN_EXPORT int tenGradientJitter(Nrrd *nout, const Nrrd *nin, double dist);
1285 TEN_EXPORT int tenGradientBalance(Nrrd *nout, const Nrrd *nin,
1286                                   tenGradientParm *tgparm);
1287 TEN_EXPORT void tenGradientMeasure(double *pot, double *minAngle,
1288                                    double *minEdge,
1289                                    const Nrrd *npos, tenGradientParm *tgparm,
1290                                    int edgeNormalize);
1291 TEN_EXPORT int tenGradientDistribute(Nrrd *nout, const Nrrd *nin,
1292                                      tenGradientParm *tgparm);
1293 TEN_EXPORT int tenGradientGenerate(Nrrd *nout, unsigned int num,
1294                                    tenGradientParm *tgparm);
1295 
1296 /* enumsTen.c */
1297 TEN_EXPORT const airEnum *const tenAniso;
1298 TEN_EXPORT const airEnum *const tenInterpType;
1299 TEN_EXPORT const airEnum _tenGage;
1300 TEN_EXPORT const airEnum *const tenGage;
1301 TEN_EXPORT const airEnum *const tenFiberType;
1302 TEN_EXPORT const airEnum *const tenDwiFiberType;
1303 TEN_EXPORT const airEnum *const tenFiberStop;
1304 TEN_EXPORT const airEnum *const tenFiberIntg;
1305 TEN_EXPORT const airEnum *const tenGlyphType;
1306 TEN_EXPORT const airEnum *const tenEstimate1Method;
1307 TEN_EXPORT const airEnum *const tenEstimate2Method;
1308 TEN_EXPORT const airEnum *const tenTripleType;
1309 
1310 /* path.c */
1311 TEN_EXPORT tenInterpParm *tenInterpParmNew(void);
1312 TEN_EXPORT tenInterpParm *tenInterpParmCopy(tenInterpParm *tip);
1313 TEN_EXPORT int tenInterpParmBufferAlloc(tenInterpParm *tip,
1314                                         unsigned int num);
1315 TEN_EXPORT tenInterpParm *tenInterpParmNix(tenInterpParm *tip);
1316 TEN_EXPORT void tenInterpTwo_d(double oten[7],
1317                                const double tenA[7],
1318                                const double tenB[7],
1319                                int ptype, double aa,
1320                                tenInterpParm *tip);
1321 TEN_EXPORT int tenInterpN_d(double tenOut[7],
1322                             const double *tenIn,
1323                             const double *wght,
1324                             unsigned int num, int ptype, tenInterpParm *tip);
1325 TEN_EXPORT double tenInterpPathLength(Nrrd *npath, int doubleVerts,
1326                                       int fancy, int shape);
1327 TEN_EXPORT int tenInterpTwoDiscrete_d(Nrrd *nout,
1328                                       const double tenA[7],
1329                                       const double tenB[7],
1330                                       int ptype, unsigned int num,
1331                                       tenInterpParm *tip);
1332 TEN_EXPORT double tenInterpDistanceTwo_d(const double tenA[7],
1333                                          const double tenB[7],
1334                                          int ptype, tenInterpParm *tip);
1335 TEN_EXPORT int tenInterpMulti3D(Nrrd *nout, const Nrrd *const *nin,
1336                                 const double *wght,
1337                                 unsigned int ninNum,
1338                                 int ptype, tenInterpParm *tip);
1339 
1340 /* glyph.c */
1341 TEN_EXPORT tenGlyphParm *tenGlyphParmNew(void);
1342 TEN_EXPORT tenGlyphParm *tenGlyphParmNix(tenGlyphParm *parm);
1343 TEN_EXPORT int tenGlyphParmCheck(tenGlyphParm *parm,
1344                                  const Nrrd *nten, const Nrrd *npos,
1345                                  const Nrrd *nslc);
1346 TEN_EXPORT int tenGlyphGen(limnObject *glyphs, echoScene *scene,
1347                            tenGlyphParm *parm,
1348                            const Nrrd *nten, const Nrrd *npos,
1349                            const Nrrd *nslc);
1350 TEN_EXPORT unsigned int tenGlyphBqdZoneEval(const double eval[3]);
1351 TEN_EXPORT void tenGlyphBqdUvEval(double uv[2], const double eval[3]);
1352 TEN_EXPORT void tenGlyphBqdEvalUv(double eval[3], const double uv[2]);
1353 TEN_EXPORT unsigned int tenGlyphBqdZoneUv(const double uv[2]);
1354 TEN_EXPORT void tenGlyphBqdAbcUv(double abc[3], const double uv[2],
1355                                  double betaMax);
1356 
1357 /* tensor.c */
1358 TEN_EXPORT int tenVerbose;
1359 TEN_EXPORT void tenRotateSingle_f(float tenOut[7],
1360                                   const float rot[9], const float tenIn[7]);
1361 TEN_EXPORT int tenTensorCheck(const Nrrd *nin,
1362                               int wantType, int want4D, int useBiff);
1363 TEN_EXPORT int tenMeasurementFrameReduce(Nrrd *nout, const Nrrd *nin);
1364 TEN_EXPORT int tenExpand2D(Nrrd *nout, const Nrrd *nin,
1365                            double scale, double thresh);
1366 TEN_EXPORT int tenExpand(Nrrd *tnine, const Nrrd *tseven,
1367                          double scale, double thresh);
1368 TEN_EXPORT int tenShrink(Nrrd *tseven, const Nrrd *nconf, const Nrrd *tnine);
1369 TEN_EXPORT int tenEigensolve_f(float eval[3], float evec[9],
1370                                const float ten[7]);
1371 TEN_EXPORT int tenEigensolve_d(double eval[3], double evec[9],
1372                                const double ten[7]);
1373 TEN_EXPORT void tenMakeSingle_f(float ten[7],
1374                                 float conf, const float eval[3],
1375                                 const float evec[9]);
1376 TEN_EXPORT void tenMakeSingle_d(double ten[7],
1377                                 double conf, const double eval[3],
1378                                 const double evec[9]);
1379 TEN_EXPORT int tenMake(Nrrd *nout, const Nrrd *nconf,
1380                        const Nrrd *neval, const Nrrd *nevec);
1381 TEN_EXPORT int tenSlice(Nrrd *nout, const Nrrd *nten,
1382                         unsigned int axis, size_t pos, unsigned int dim);
1383 TEN_EXPORT void tenInvariantGradientsK_d(double K1[7],
1384                                          double K2[7],
1385                                          double K3[7],
1386                                          const double ten[7],
1387                                          const double minnorm);
1388 TEN_EXPORT void tenInvariantGradientsR_d(double R1[7],
1389                                          double R2[7],
1390                                          double R3[7],
1391                                          const double ten[7],
1392                                          const double minnorm);
1393 TEN_EXPORT void tenRotationTangents_d(double phi1[7],
1394                                       double phi2[7],
1395                                       double phi3[7],
1396                                       const double evec[9]);
1397 TEN_EXPORT void tenLogSingle_d(double logten[7], const double ten[7]);
1398 TEN_EXPORT void tenLogSingle_f(float logten[7], const float ten[7]);
1399 TEN_EXPORT void tenExpSingle_d(double expten[7], const double ten[7]);
1400 TEN_EXPORT void tenExpSingle_f(float expten[7], const float ten[7]);
1401 TEN_EXPORT void tenSqrtSingle_d(double sqrtten[7], const double ten[7]);
1402 TEN_EXPORT void tenSqrtSingle_f(float sqrtten[7], const float ten[7]);
1403 TEN_EXPORT void tenPowSingle_d(double powten[7], const double ten[7],
1404                                double power);
1405 TEN_EXPORT void tenPowSingle_f(float powten[7], const float ten[7],
1406                                float power);
1407 /* should rename to tenInvSingle_x */
1408 TEN_EXPORT void tenInv_f(float inv[7], const float ten[7]);
1409 TEN_EXPORT void tenInv_d(double inv[7], const double ten[7]);
1410 TEN_EXPORT double tenDoubleContract_d(double a[7], double T[21], double b[7]);
1411 
1412 /* chan.c */
1413 /* old tenCalc* functions superceded/deprecated by new tenEstimate* code */
1414 TEN_EXPORT const char *tenDWMRIModalityKey;
1415 TEN_EXPORT const char *tenDWMRIModalityVal;
1416 TEN_EXPORT const char *tenDWMRINAVal;
1417 TEN_EXPORT const char *tenDWMRIBValueKey;
1418 TEN_EXPORT const char *tenDWMRIGradKeyFmt;
1419 TEN_EXPORT const char *tenDWMRIBmatKeyFmt;
1420 TEN_EXPORT const char *tenDWMRINexKeyFmt;
1421 TEN_EXPORT const char *tenDWMRISkipKeyFmt;
1422 TEN_EXPORT int tenDWMRIKeyValueParse(Nrrd **ngradP, Nrrd **nbmatP, double *bP,
1423                                      unsigned int **skip,
1424                                      unsigned int *skipNum,
1425                                      const Nrrd *ndwi);
1426 TEN_EXPORT int tenBMatrixCalc(Nrrd *nbmat, const Nrrd *ngrad);
1427 TEN_EXPORT int tenEMatrixCalc(Nrrd *nemat, const Nrrd *nbmat, int knownB0);
1428 TEN_EXPORT void tenEstimateLinearSingle_f(float *ten, float *B0P,
1429                                           const float *dwi, const double *emat,
1430                                           double *vbuf, unsigned int DD,
1431                                           int knownB0, float thresh,
1432                                           float soft, float b);
1433 TEN_EXPORT void tenEstimateLinearSingle_d(double *ten, double *B0P,
1434                                           const double *dwi, const double*emat,
1435                                           double *vbuf, unsigned int DD,
1436                                           int knownB0, double thresh,
1437                                           double soft, double b);
1438 TEN_EXPORT int tenEstimateLinear3D(Nrrd *nten, Nrrd **nterrP, Nrrd **nB0P,
1439                                    const Nrrd *const *ndwi,
1440                                    unsigned int dwiLen,
1441                                    const Nrrd *nbmat, int knownB0,
1442                                    double thresh, double soft, double b);
1443 TEN_EXPORT int tenEstimateLinear4D(Nrrd *nten, Nrrd **nterrP, Nrrd **nB0P,
1444                                    const Nrrd *ndwi, const Nrrd *_nbmat,
1445                                    int knownB0,
1446                                    double thresh, double soft, double b);
1447 TEN_EXPORT void tenSimulateSingle_f(float *dwi, float B0, const float *ten,
1448                                     const double *bmat, unsigned int DD,
1449                                     float b);
1450 TEN_EXPORT int tenSimulate(Nrrd *ndwi, const Nrrd *nT2, const Nrrd *nten,
1451                            const Nrrd *nbmat, double b);
1452 
1453 /* estimate.c */
1454 TEN_EXPORT tenEstimateContext *tenEstimateContextNew(void);
1455 TEN_EXPORT void tenEstimateVerboseSet(tenEstimateContext *tec,
1456                                       int verbose);
1457 TEN_EXPORT void tenEstimateNegEvalShiftSet(tenEstimateContext *tec,
1458                                            int doit);
1459 TEN_EXPORT int tenEstimateMethodSet(tenEstimateContext *tec,
1460                                     int estMethod);
1461 TEN_EXPORT int tenEstimateSigmaSet(tenEstimateContext *tec,
1462                                    double sigma);
1463 TEN_EXPORT int tenEstimateValueMinSet(tenEstimateContext *tec,
1464                                       double valueMin);
1465 TEN_EXPORT int tenEstimateGradientsSet(tenEstimateContext *tec,
1466                                        const Nrrd *ngrad,
1467                                        double bValue, int estimateB0);
1468 TEN_EXPORT int tenEstimateBMatricesSet(tenEstimateContext *tec,
1469                                        const Nrrd *nbmat,
1470                                        double bValue, int estimateB0);
1471 TEN_EXPORT int tenEstimateSkipSet(tenEstimateContext *tec,
1472                                   unsigned int valIdx,
1473                                   int doSkip);
1474 TEN_EXPORT int tenEstimateSkipReset(tenEstimateContext *tec);
1475 TEN_EXPORT int tenEstimateThresholdSet(tenEstimateContext *tec,
1476                                        double thresh, double soft);
1477 TEN_EXPORT int tenEstimateUpdate(tenEstimateContext *tec);
1478 TEN_EXPORT int tenEstimate1TensorSimulateSingle_f(tenEstimateContext *tec,
1479                                                   float *simval,
1480                                                   float sigma,
1481                                                   float bValue, float B0,
1482                                                   const float _ten[7]);
1483 TEN_EXPORT int tenEstimate1TensorSimulateSingle_d(tenEstimateContext *tec,
1484                                                   double *simval,
1485                                                   double sigma,
1486                                                   double bValue, double B0,
1487                                                   const double ten[7]);
1488 TEN_EXPORT int tenEstimate1TensorSimulateVolume(tenEstimateContext *tec,
1489                                                 Nrrd *ndwi,
1490                                                 double sigma, double bValue,
1491                                                 const Nrrd *nB0,
1492                                                 const Nrrd *nten,
1493                                                 int outType,
1494                                                 int keyValueSet);
1495 TEN_EXPORT int tenEstimate1TensorSingle_f(tenEstimateContext *tec,
1496                                           float ten[7], const float *all);
1497 TEN_EXPORT int tenEstimate1TensorSingle_d(tenEstimateContext *tec,
1498                                           double ten[7], const double *all);
1499 TEN_EXPORT int tenEstimate1TensorVolume4D(tenEstimateContext *tec,
1500                                           Nrrd *nten,
1501                                           Nrrd **nB0P, Nrrd **nterrP,
1502                                           const Nrrd *ndwi, int outType);
1503 TEN_EXPORT tenEstimateContext *tenEstimateContextNix(tenEstimateContext *tec);
1504 
1505 /* aniso.c */
1506 TEN_EXPORT float (*_tenAnisoEval_f[TEN_ANISO_MAX+1])(const float eval[3]);
1507 TEN_EXPORT float tenAnisoEval_f(const float eval[3], int aniso);
1508 TEN_EXPORT double (*_tenAnisoEval_d[TEN_ANISO_MAX+1])(const double eval[3]);
1509 TEN_EXPORT double tenAnisoEval_d(const double eval[3], int aniso);
1510 
1511 TEN_EXPORT float (*_tenAnisoTen_f[TEN_ANISO_MAX+1])(const float ten[7]);
1512 TEN_EXPORT float tenAnisoTen_f(const float ten[7], int aniso);
1513 TEN_EXPORT double (*_tenAnisoTen_d[TEN_ANISO_MAX+1])(const double ten[7]);
1514 TEN_EXPORT double tenAnisoTen_d(const double ten[7], int aniso);
1515 
1516 TEN_EXPORT int tenAnisoPlot(Nrrd *nout, int aniso, unsigned int res,
1517                             int hflip, int whole, int nanout);
1518 TEN_EXPORT int tenAnisoVolume(Nrrd *nout, const Nrrd *nin,
1519                               int aniso, double confThresh);
1520 TEN_EXPORT int tenAnisoHistogram(Nrrd *nout, const Nrrd *nin,
1521                                  const Nrrd *nwght, int right,
1522                                  int version, unsigned int resolution);
1523 TEN_EXPORT tenEvecRGBParm *tenEvecRGBParmNew(void);
1524 TEN_EXPORT tenEvecRGBParm *tenEvecRGBParmNix(tenEvecRGBParm *rgbp);
1525 TEN_EXPORT int tenEvecRGBParmCheck(const tenEvecRGBParm *rgbp);
1526 TEN_EXPORT void tenEvecRGBSingle_f(float RGB[3], float conf,
1527                                    const float eval[3], const float evec[3],
1528                                    const tenEvecRGBParm *rgbp);
1529 TEN_EXPORT void tenEvecRGBSingle_d(double RGB[3], double conf,
1530                                    const double eval[3], const double evec[3],
1531                                    const tenEvecRGBParm *rgbp);
1532 
1533 /* miscTen.c */
1534 TEN_EXPORT int tenEvecRGB(Nrrd *nout, const Nrrd *nin,
1535                           const tenEvecRGBParm *rgbp);
1536 TEN_EXPORT short tenEvqSingle_f(float vec[3], float scl);
1537 TEN_EXPORT int tenEvqVolume(Nrrd *nout, const Nrrd *nin, int which,
1538                             int aniso, int scaleByAniso);
1539 TEN_EXPORT int tenBMatrixCheck(const Nrrd *nbmat,
1540                                int type, unsigned int minnum);
1541 TEN_EXPORT int _tenFindValley(double *valP, const Nrrd *nhist,
1542                               double tweak, int save);
1543 
1544 /* fiberMethods.c */
1545 TEN_EXPORT void tenFiberSingleInit(tenFiberSingle *tfbs);
1546 TEN_EXPORT void tenFiberSingleDone(tenFiberSingle *tfbs);
1547 TEN_EXPORT tenFiberSingle *tenFiberSingleNew(void);
1548 TEN_EXPORT tenFiberSingle *tenFiberSingleNix(tenFiberSingle *tfbs);
1549 TEN_EXPORT tenFiberContext *tenFiberContextNew(const Nrrd *dtvol);
1550 TEN_EXPORT tenFiberContext *tenFiberContextDwiNew(const Nrrd *dwivol,
1551                                                   double thresh,
1552                                                   double soft,
1553                                                   double valueMin,
1554                                                   int ten1method,
1555                                                   int ten2method);
1556 TEN_EXPORT void tenFiberVerboseSet(tenFiberContext *tfx, int verbose);
1557 TEN_EXPORT int tenFiberTypeSet(tenFiberContext *tfx, int type);
1558 TEN_EXPORT int tenFiberKernelSet(tenFiberContext *tfx,
1559                                  const NrrdKernel *kern,
1560                                  const double parm[NRRD_KERNEL_PARMS_NUM]);
1561 TEN_EXPORT int tenFiberProbeItemSet(tenFiberContext *tfx, int item);
1562 TEN_EXPORT int tenFiberIntgSet(tenFiberContext *tfx, int intg);
1563 TEN_EXPORT int tenFiberStopSet(tenFiberContext *tfx, int stop, ...);
1564 TEN_EXPORT int tenFiberStopAnisoSet(tenFiberContext *tfx,
1565                                     int anisoType, double anisoThresh);
1566 TEN_EXPORT int tenFiberStopDoubleSet(tenFiberContext *tfx,
1567                                      int stop, double val);
1568 TEN_EXPORT int tenFiberStopUIntSet(tenFiberContext *tfx,
1569                                    int stop, unsigned int val);
1570 TEN_EXPORT void tenFiberStopOn(tenFiberContext *tfx, int stop);
1571 TEN_EXPORT void tenFiberStopOff(tenFiberContext *tfx, int stop);
1572 TEN_EXPORT void tenFiberStopReset(tenFiberContext *tfx);
1573 TEN_EXPORT int tenFiberAnisoSpeedSet(tenFiberContext *tfx, int aniso,
1574                                      double lerp, double thresh, double soft);
1575 TEN_EXPORT int tenFiberAnisoSpeedReset(tenFiberContext *tfx);
1576 TEN_EXPORT int tenFiberParmSet(tenFiberContext *tfx, int parm, double val);
1577 TEN_EXPORT int tenFiberUpdate(tenFiberContext *tfx);
1578 TEN_EXPORT tenFiberContext *tenFiberContextCopy(tenFiberContext *tfx);
1579 TEN_EXPORT tenFiberContext *tenFiberContextNix(tenFiberContext *tfx);
1580 
1581 /* fiber.c */
1582 TEN_EXPORT int tenFiberTraceSet(tenFiberContext *tfx, Nrrd *nfiber,
1583                                 double *buff, unsigned int halfBuffLen,
1584                                 unsigned int *startIdxP, unsigned int *endIdxP,
1585                                 double seed[3]);
1586 TEN_EXPORT int tenFiberTrace(tenFiberContext *tfx,
1587                              Nrrd *nfiber, double seed[3]);
1588 TEN_EXPORT unsigned int tenFiberDirectionNumber(tenFiberContext *tfx,
1589                                                 double seed[3]);
1590 TEN_EXPORT int tenFiberSingleTrace(tenFiberContext *tfx, tenFiberSingle *tfbs,
1591                                    double seed[3], unsigned int which);
1592 TEN_EXPORT tenFiberMulti *tenFiberMultiNew(void);
1593 TEN_EXPORT tenFiberMulti *tenFiberMultiNix(tenFiberMulti *tfm);
1594 TEN_EXPORT int tenFiberMultiTrace(tenFiberContext *tfx, tenFiberMulti *tfml,
1595                                   const Nrrd *nseed);
1596 TEN_EXPORT int tenFiberMultiPolyData(tenFiberContext *tfx,
1597                                      limnPolyData *lpld, tenFiberMulti *tfml);
1598 TEN_EXPORT int tenFiberMultiProbeVals(tenFiberContext *tfx,
1599                                       Nrrd *nval, tenFiberMulti *tfml);
1600 
1601 /* epireg.c */
1602 TEN_EXPORT int _tenEpiRegThresholdFind(double *DWthrP, Nrrd **nin,
1603                                        int ninLen, int save, double expo);
1604 TEN_EXPORT int tenEpiRegister3D(Nrrd **nout, Nrrd **ndwi,
1605                                 unsigned int dwiLen, Nrrd *ngrad,
1606                                 int reference,
1607                                 double bwX, double bwY,
1608                                 double fitFrac, double DWthr,
1609                                 int doCC,
1610                                 const NrrdKernel *kern, double *kparm,
1611                                 int progress, int verbose);
1612 TEN_EXPORT int tenEpiRegister4D(Nrrd *nout, Nrrd *nin, Nrrd *ngrad,
1613                                 int reference,
1614                                 double bwX, double bwY,
1615                                 double fitFrac, double DWthr,
1616                                 int doCC,
1617                                 const NrrdKernel *kern, double *kparm,
1618                                 int progress, int verbose);
1619 
1620 /* experSpec.c */
1621 TEN_EXPORT tenExperSpec* tenExperSpecNew(void);
1622 TEN_EXPORT int tenExperSpecGradSingleBValSet(tenExperSpec *espec,
1623                                              int insertB0,
1624                                              double bval,
1625                                              const double *grad,
1626                                              unsigned int gradNum);
1627 TEN_EXPORT int tenExperSpecGradBValSet(tenExperSpec *espec,
1628                                        int insertB0,
1629                                        const double *bval,
1630                                        const double *grad,
1631                                        unsigned int bgNum);
1632 /*
1633 TEN_EXPORT int tenExperSpecGradBValWghtSet(tenExperSpec *espec,
1634                                            int insertB0,
1635                                            const double *bval,
1636                                            const double *grad,
1637                                            const double *wght,
1638                                            unsigned int bgwNum);
1639 */
1640 TEN_EXPORT int tenExperSpecFromKeyValueSet(tenExperSpec *espec,
1641                                            const Nrrd *ndwi);
1642 TEN_EXPORT tenExperSpec* tenExperSpecNix(tenExperSpec *espec);
1643 TEN_EXPORT double tenExperSpecKnownB0Get(const tenExperSpec *espec,
1644                                          const double *dwi);
1645 TEN_EXPORT double tenExperSpecMaxBGet(const tenExperSpec *espec);
1646 TEN_EXPORT int tenDWMRIKeyValueFromExperSpecSet(Nrrd *ndwi,
1647                                                 const tenExperSpec *espec);
1648 
1649 /* tenModel.c */
1650 TEN_EXPORT const char *tenModelPrefixStr;
1651 TEN_EXPORT int tenModelParse(const tenModel **model, int *plusB0,
1652                              int requirePrefix, const char *str);
1653 TEN_EXPORT int tenModelFromAxisLearnPossible(const NrrdAxisInfo *axinfo);
1654 TEN_EXPORT int tenModelFromAxisLearn(const tenModel **model, int *plusB0,
1655                                      const NrrdAxisInfo *axinfo);
1656 TEN_EXPORT int tenModelSimulate(Nrrd *ndwi, int typeOut,
1657                                 tenExperSpec *espec,
1658                                 const tenModel *model,
1659                                 const Nrrd *nB0, /* maybe NULL */
1660                                 const Nrrd *nparm,
1661                                 int keyValueSet);
1662 TEN_EXPORT int tenModelSqeFit(Nrrd *nparm,
1663                               Nrrd **nsqeP, Nrrd **nconvP, Nrrd **niterP,
1664                               const tenModel *model,
1665                               const tenExperSpec *espec, const Nrrd *ndwi,
1666                               int knownB0, int saveB0, int typeOut,
1667                               unsigned int minIter, unsigned int maxIter,
1668                               unsigned int starts, double convEps,
1669                               airRandMTState *rng,
1670                               int verbose);
1671 TEN_EXPORT int tenModelNllFit(Nrrd *nparm, Nrrd **nnllP,
1672                               const tenModel *model,
1673                               const tenExperSpec *espec, const Nrrd *ndwi,
1674                               int rician, double sigma, int knownB0);
1675 TEN_EXPORT int tenModelConvert(Nrrd *nparmDst, int *convRet,
1676                                const tenModel *modelDst,
1677                                const Nrrd *nparmSrc,
1678                                const tenModel *modelSrc);
1679 
1680 /*
1681 ** Have to keep this list of model declarations in sync with:
1682 ** * tenModel.c/str2model()
1683 ** * all model*.c/parmConvert()  (fallen behind here)
1684 ** * sources.cmake, GNUmakefile
1685 **
1686 ** Note that the TEN_MODEL_STR_* strings should be all lowercase,
1687 ** as long as we want the model parsing to be case insensitive, and
1688 ** as long as the logic in tenModel.c/str2model() is so simplistic
1689 */
1690 /* modelZero.c */
1691 #define TEN_MODEL_STR_ZERO "zero"
1692 TEN_EXPORT const tenModel *const tenModelZero;
1693 /* modelB0.c */
1694 #define TEN_MODEL_STR_B0 "b0"
1695 TEN_EXPORT const tenModel *const tenModelB0;
1696 /* modelBall.c */
1697 #define TEN_MODEL_STR_BALL "ball"
1698 TEN_EXPORT const tenModel *const tenModelBall;
1699 /* model1Vector2D.c */
1700 #define TEN_MODEL_STR_1VECTOR2D "1vector2d"
1701 TEN_EXPORT const tenModel *const tenModel1Vector2D;
1702 /* model1Unit2D.c */
1703 #define TEN_MODEL_STR_1UNIT2D "1unit2d"
1704 TEN_EXPORT const tenModel *const tenModel1Unit2D;
1705 /* model2Unit2D.c */
1706 #define TEN_MODEL_STR_2UNIT2D "2unit2d"
1707 TEN_EXPORT const tenModel *const tenModel2Unit2D;
1708 /* model1Stick.c */
1709 #define TEN_MODEL_STR_1STICK "1stick"
1710 TEN_EXPORT const tenModel *const tenModel1Stick;
1711 /* modelBall1StickEMD.c */
1712 #define TEN_MODEL_STR_BALL1STICKEMD "ball1stickemd"
1713 TEN_EXPORT const tenModel *const tenModelBall1StickEMD;
1714 /* modelBall1Stick.c */
1715 #define TEN_MODEL_STR_BALL1STICK "ball1stick"
1716 TEN_EXPORT const tenModel *const tenModelBall1Stick;
1717 /* modelBall1Cylinder.c */
1718 #define TEN_MODEL_STR_BALL1CYLINDER "ball1cylinder"
1719 TEN_EXPORT const tenModel *const tenModelBall1Cylinder;
1720 /* model1Cylinder.c */
1721 #define TEN_MODEL_STR_1CYLINDER "1cylinder"
1722 TEN_EXPORT const tenModel *const tenModel1Cylinder;
1723 /* model1Tensor2.c: 2nd-order tensor (one of them), not two-tensor */
1724 #define TEN_MODEL_STR_1TENSOR2 "1tensor2"
1725 TEN_EXPORT const tenModel *const tenModel1Tensor2;
1726 
1727 /* mod.c */
1728 TEN_EXPORT int tenSizeNormalize(Nrrd *nout, const Nrrd *nin, double weight[3],
1729                                 double amount, double target);
1730 TEN_EXPORT int tenSizeScale(Nrrd *nout, const Nrrd *nin, double amount);
1731 TEN_EXPORT int tenAnisoScale(Nrrd *nout, const Nrrd *nin, double scale,
1732                              int fixDet, int makePositive);
1733 TEN_EXPORT int tenEigenvaluePower(Nrrd *nout, const Nrrd *nin, double expo);
1734 TEN_EXPORT int tenEigenvalueClamp(Nrrd *nout, const Nrrd *nin,
1735                                   double min, double max);
1736 TEN_EXPORT int tenEigenvalueAdd(Nrrd *nout, const Nrrd *nin, double val);
1737 TEN_EXPORT int tenEigenvalueMultiply(Nrrd *nout, const Nrrd *nin, double val);
1738 TEN_EXPORT int tenLog(Nrrd *nout, const Nrrd *nin);
1739 TEN_EXPORT int tenExp(Nrrd *nout, const Nrrd *nin);
1740 
1741 /* bvec.c */
1742 TEN_EXPORT int tenBVecNonLinearFit(Nrrd *nout, const Nrrd *nin,
1743                                    double *bb, double *ww,
1744                                    int iterMax, double eps);
1745 
1746 /* tenGage.c */
1747 TEN_EXPORT gageKind *tenGageKind;
1748 
1749 /* tenDwiGage.c */
1750 /* we can't declare or define a tenDwiGageKind->name (analogous to
1751    tenGageKind->name or gageSclKind->name) because the DWI kind is
1752    dynamically allocated, but at least we can declare a cannonical
1753    name (HEY: ugly) */
1754 #define TEN_DWI_GAGE_KIND_NAME "dwi"
1755 TEN_EXPORT const airEnum _tenDwiGage;
1756 TEN_EXPORT const airEnum *const tenDwiGage;
1757 TEN_EXPORT gageKind *tenDwiGageKindNew(void);
1758 TEN_EXPORT gageKind *tenDwiGageKindNix(gageKind *dwiKind);
1759 /* warning: this function will likely change its arguments in the future */
1760 TEN_EXPORT int tenDwiGageKindSet(gageKind *dwiKind,
1761                                  double thresh, double soft,
1762                                  double bval, double valueMin,
1763                                  const Nrrd *ngrad,
1764                                  const Nrrd *nbmat,
1765                                  int emethod1, int emethod2,
1766                                  unsigned int randSeed);
1767 TEN_EXPORT int tenDwiGageKindCheck(const gageKind *kind);
1768 
1769 /* bimod.c */
1770 TEN_EXPORT tenEMBimodalParm* tenEMBimodalParmNew(void);
1771 TEN_EXPORT tenEMBimodalParm* tenEMBimodalParmNix(tenEMBimodalParm *biparm);
1772 TEN_EXPORT int tenEMBimodal(tenEMBimodalParm *biparm, const Nrrd *nhisto);
1773 
1774 /* tend{Flotsam,Anplot,Anvol,Evec,Eval,. . .}.c */
1775 #define TEND_DECLARE(C) TEN_EXPORT unrrduCmd tend_##C##Cmd;
1776 #define TEND_LIST(C) &tend_##C##Cmd,
1777 /* removed from below (superseded by estim): F(calc) \ */
1778 #define TEND_MAP(F) \
1779 F(about) \
1780 F(grads) \
1781 F(epireg) \
1782 F(bmat) \
1783 F(estim) \
1784 F(sim) \
1785 F(mfit) \
1786 F(mconv) \
1787 F(msim) \
1788 F(make) \
1789 F(avg) \
1790 F(helix) \
1791 F(sten) \
1792 F(glyph) \
1793 F(ellipse) \
1794 F(anplot) \
1795 F(anvol) \
1796 F(anscale) \
1797 F(anhist) \
1798 F(triple) \
1799 F(tconv) \
1800 F(point) \
1801 F(slice) \
1802 F(norm) \
1803 F(fiber) \
1804 F(eval) \
1805 F(evalpow) \
1806 F(evalclamp) \
1807 F(evaladd) \
1808 F(evalmult) \
1809 F(log) \
1810 F(exp) \
1811 F(evec) \
1812 F(evecrgb) \
1813 F(evq) \
1814 F(unmf) \
1815 F(expand) \
1816 F(shrink) \
1817 F(bfit) \
1818 F(satin)
1819 TEND_MAP(TEND_DECLARE)
1820 TEN_EXPORT unrrduCmd *tendCmdList[];
1821 TEN_EXPORT hestCB *tendFiberStopCB;
1822 TEN_EXPORT const char *tendTitle;
1823 
1824 #ifdef __cplusplus
1825 }
1826 #endif
1827 
1828 #endif /* TEN_HAS_BEEN_INCLUDED */
1829