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_PRIVATE_HAS_BEEN_INCLUDED 25 #define TEN_PRIVATE_HAS_BEEN_INCLUDED 26 27 #ifdef __cplusplus 28 extern "C" { 29 #endif 30 31 #define TEND_CMD(name, info) \ 32 unrrduCmd tend_##name##Cmd = { #name, info, tend_##name##Main } 33 34 /* USAGE, PARSE: both copied verbatim from unrrdu/privateUnrrdu.h, but 35 ** then some hacking was added . . . 36 */ 37 #define USAGE(info) \ 38 if (!argc) { \ 39 hestInfo(stdout, me, (info), hparm); \ 40 hestUsage(stdout, hopt, me, hparm); \ 41 hestGlossary(stdout, hopt, hparm); \ 42 airMopError(mop); \ 43 return 0; \ 44 } 45 46 /* JUSTPARSE is called by the tend functions that do *not* take an 47 ** input 7-component tensor volume 48 */ 49 #define JUSTPARSE() \ 50 if ((pret=hestParse(hopt, argc, argv, &perr, hparm))) { \ 51 if (1 == pret) { \ 52 fprintf(stderr, "%s: %s\n", me, perr); free(perr); \ 53 hestUsage(stderr, hopt, me, hparm); \ 54 airMopError(mop); \ 55 return 2; \ 56 } else { \ 57 exit(1); \ 58 } \ 59 } 60 61 /* 62 ** PARSE is called by tend functions that do take a 7-component tensor 63 ** volume, so that as a hack, we can process 6-component volumes as well, 64 ** by padding on the confidence channel (fixed at 1.0) 65 */ 66 #define PARSE() \ 67 JUSTPARSE(); \ 68 if (4 == nin->dim \ 69 && 6 == nin->axis[0].size \ 70 && nrrdTypeBlock != nin->type) { \ 71 ptrdiff_t padmin[4], padmax[4]; \ 72 Nrrd *npadtmp; \ 73 /* create a confidence channel by padding on 1s */ \ 74 padmin[0] = -1; padmin[1] = padmin[2] = padmin[3] = 0; \ 75 padmax[0] = nin->axis[0].size-1; \ 76 padmax[1] = nin->axis[1].size-1; \ 77 padmax[2] = nin->axis[2].size-1; \ 78 padmax[3] = nin->axis[3].size-1; \ 79 npadtmp = nrrdNew(); \ 80 if (nrrdPad_nva(npadtmp, nin, padmin, padmax, nrrdBoundaryPad, 1.0) \ 81 || nrrdCopy(nin, npadtmp)) { \ 82 char *paderr; \ 83 airMopAdd(mop, paderr = biffGetDone(NRRD), airFree, airMopAlways);\ 84 fprintf(stderr, "%s: can't pad 6-comp tensor:\n%s", me, paderr); \ 85 airMopError(mop); \ 86 nrrdNuke(npadtmp); \ 87 return 2; \ 88 } \ 89 nrrdNuke(npadtmp); \ 90 } 91 92 /* qseg.c: 2-tensor estimation */ 93 extern void _tenQball(const double b, const int gradcount, 94 const double svals[], const double grads[], 95 double qvals[] ); 96 extern void _tenSegsamp2(const int gradcount, const double qvals[], 97 const double grads[], const double qpoints[], 98 unsigned int seg[], double dists[] ); 99 extern void _tenCalcdists(const int centcount, const double centroid[6], 100 const int gradcount, const double qpoints[], 101 double dists[] ); 102 extern void _tenInitcent2(const int gradcount, const double qvals[], 103 const double grads[], double centroids[6] ); 104 extern int _tenCalccent2(const int gradcount, const double qpoints[], 105 const double dists[], double centroid[6], 106 unsigned int seg[] ); 107 extern void _tenSeg2weights(const int gradcount, const int seg[], 108 const int segcount, double weights[] ); 109 extern void _tenQvals2points(const int gradcount, const double qvals[], 110 const double grads[], double qpoints[] ); 111 extern double _tenPldist(const double point[], const double line[] ); 112 113 /* arishFuncs.c: Arish's implementation of Peled's 2-tensor fit */ 114 #define VEC_SIZE 3 115 116 extern void twoTensFunc(double *p, double *x, int m, int n, void *data); 117 extern void formTensor2D(double ten[7], double lam1, double lam3, double phi); 118 119 /* qglox.c: stuff for quaternion geodesic-loxodromes that has dubious 120 utility for the general public */ 121 extern void tenQGLInterpTwoEvalK(double oeval[3], 122 const double evalA[3], 123 const double evalB[3], 124 const double tt); 125 extern void tenQGLInterpTwoEvalR(double oeval[3], 126 const double evalA[3], 127 const double evalB[3], 128 const double tt); 129 extern void tenQGLInterpTwoEvec(double oevec[9], 130 const double evecA[9], const double evecB[9], 131 double tt); 132 extern void tenQGLInterpTwo(double oten[7], 133 const double tenA[7], const double tenB[7], 134 int ptype, double aa, tenInterpParm *tip); 135 extern int _tenQGLInterpNEval(double evalOut[3], 136 const double *evalIn, /* size 3 -by- NN */ 137 const double *wght, /* size NN */ 138 unsigned int NN, 139 int ptype, tenInterpParm *tip); 140 extern int _tenQGLInterpNEvec(double evecOut[9], 141 const double *evecIn, /* size 9 -by- NN */ 142 const double *wght, /* size NN */ 143 unsigned int NN, 144 tenInterpParm *tip); 145 extern int tenQGLInterpN(double tenOut[7], 146 const double *tenIn, 147 const double *wght, 148 unsigned int NN, int ptype, tenInterpParm *tip); 149 150 /* experSpec.c */ 151 TEN_EXPORT double _tenExperSpec_sqe(const double *dwiMeas, 152 const double *dwiSim, 153 const tenExperSpec *espec, 154 int knownB0); 155 TEN_EXPORT double _tenExperSpec_nll(const double *dwiMeas, 156 const double *dwiSim, 157 const tenExperSpec *espec, 158 int rician, double sigma, 159 int knownB0); 160 161 /* tenModel.c */ 162 TEN_EXPORT double _tenModelSqeFitSingle(const tenModel *model, 163 double *testParm, double *grad, 164 double *parm, double *convFrac, 165 unsigned int *itersTaken, 166 const tenExperSpec *espec, 167 double *dwiBuff, 168 const double *dwiMeas, 169 const double *parmInit, 170 int knownB0, 171 unsigned int minIter, 172 unsigned int maxIter, 173 double convEps, 174 int verbose); 175 176 /* model*.c */ 177 /* 178 ** NOTE: these functions rely on per-model information (like PARM_NUM) 179 ** or functions (like "simulate"), which is why these functions don't 180 ** all use some common underlying function: more information would 181 ** have to be passed (annoying for rapid hacking), and the function 182 ** call overhead might be doubled 183 */ 184 185 #define _TEN_PARM_ALLOC \ 186 static double * \ 187 parmAlloc(void) { \ 188 \ 189 return AIR_CAST(double *, calloc(PARM_NUM ? PARM_NUM : 1, sizeof(double))); \ 190 } 191 192 #define _TEN_PARM_RAND \ 193 static void \ 194 parmRand(double *parm, airRandMTState *rng, int knownB0) { \ 195 unsigned int ii; \ 196 \ 197 for (ii=knownB0 ? 1 : 0; ii<PARM_NUM; ii++) { \ 198 if (parmDesc[ii].vec3) { \ 199 /* its unit vector */ \ 200 double xx, yy, zz, theta, rr; \ 201 \ 202 zz = AIR_AFFINE(0.0, airDrandMT_r(rng), 1.0, -1.0, 1.0); \ 203 theta = AIR_AFFINE(0.0, airDrandMT_r(rng), 1.0, 0.0, 2*AIR_PI); \ 204 rr = sqrt(1 - zz*zz); \ 205 xx = rr*cos(theta); \ 206 yy = rr*sin(theta); \ 207 parm[ii + 0] = xx; \ 208 parm[ii + 1] = yy; \ 209 parm[ii + 2] = zz; \ 210 /* bump ii by 2, anticipating end of this for-loop iteration */ \ 211 ii += 2; \ 212 } else { \ 213 parm[ii] = AIR_AFFINE(0.0, airDrandMT_r(rng), 1.0, \ 214 parmDesc[ii].min, parmDesc[ii].max); \ 215 } \ 216 } \ 217 return; \ 218 } 219 220 #define _TEN_PARM_STEP \ 221 static void \ 222 parmStep(double *parm1, const double scl, \ 223 const double *grad, const double *parm0) { \ 224 unsigned int ii; \ 225 \ 226 for (ii=0; ii<PARM_NUM; ii++) { \ 227 parm1[ii] = scl*grad[ii] + parm0[ii]; \ 228 if (parmDesc[ii].cyclic) { \ 229 double delta = parmDesc[ii].max - parmDesc[ii].min; \ 230 while (parm1[ii] > parmDesc[ii].max) { \ 231 parm1[ii] -= delta; \ 232 } \ 233 while (parm1[ii] < parmDesc[ii].min) { \ 234 parm1[ii] += delta; \ 235 } \ 236 } else { \ 237 parm1[ii] = AIR_CLAMP(parmDesc[ii].min, parm1[ii], parmDesc[ii].max); \ 238 } \ 239 if (parmDesc[ii].vec3 && 2 == parmDesc[ii].vecIdx) { \ 240 double len; \ 241 ELL_3V_NORM(parm1 + ii - 2, parm1 + ii - 2, len); \ 242 } \ 243 } \ 244 } 245 246 #define _TEN_PARM_DIST \ 247 static double \ 248 parmDist(const double *parmA, const double *parmB) { \ 249 unsigned int ii; \ 250 double dist; \ 251 \ 252 dist = 0; \ 253 for (ii=0; ii<PARM_NUM; ii++) { \ 254 double dp1, delta; \ 255 delta = parmDesc[ii].max - parmDesc[ii].min; \ 256 dp1 = parmA[ii] - parmB[ii]; \ 257 if (parmDesc[ii].cyclic) { \ 258 double dp2, dp3; \ 259 dp2 = dp1 + delta; \ 260 dp3 = dp1 - delta; \ 261 dp1 *= dp1; \ 262 dp2 *= dp2; \ 263 dp3 *= dp3; \ 264 dp1 = AIR_MIN(dp1, dp2); \ 265 dp1 = AIR_MIN(dp1, dp3); \ 266 dist += dp1/(delta*delta); \ 267 } else { \ 268 dist += dp1*dp1/(delta*delta); \ 269 } \ 270 } \ 271 return sqrt(dist); \ 272 } 273 274 #define _TEN_PARM_COPY \ 275 static void \ 276 parmCopy(double *parmA, const double *parmB) { \ 277 unsigned int ii; \ 278 \ 279 for (ii=0; ii<PARM_NUM; ii++) { \ 280 parmA[ii] = parmB[ii]; \ 281 } \ 282 return; \ 283 } 284 285 #define _TEN_PARM_CONVERT_NOOP \ 286 static int \ 287 parmConvert(double *parmDst, const double *parmSrc, \ 288 const tenModel *modelSrc) { \ 289 unsigned int ii; \ 290 \ 291 AIR_UNUSED(modelSrc); \ 292 parmDst[0] = parmSrc[0]; \ 293 for (ii=1; ii<PARM_NUM; ii++) { \ 294 parmDst[ii] = AIR_NAN; \ 295 } \ 296 return 1; \ 297 } 298 299 #define _TEN_SQE \ 300 static double \ 301 sqe(const double *parm, const tenExperSpec *espec, \ 302 double *dwiBuff, const double *dwiMeas, int knownB0) { \ 303 \ 304 simulate(dwiBuff, parm, espec); \ 305 return _tenExperSpec_sqe(dwiMeas, dwiBuff, espec, knownB0); \ 306 } 307 308 #define _TEN_SQE_GRAD_CENTDIFF \ 309 static void \ 310 sqeGrad(double *grad, const double *parm0, \ 311 const tenExperSpec *espec, \ 312 double *dwiBuff, const double *dwiMeas, \ 313 int knownB0) { \ 314 double parm1[PARM_NUM ? PARM_NUM : 1], sqeForw, sqeBack, dp; \ 315 unsigned int ii, i0; \ 316 \ 317 i0 = knownB0 ? 1 : 0; \ 318 for (ii=0; ii<PARM_NUM; ii++) { \ 319 /* have to copy all parms, even B0 with knownB0, and even if we \ 320 aren't going to diddle these values, because the \ 321 simulation depends on knowing the values */ \ 322 parm1[ii] = parm0[ii]; \ 323 } \ 324 for (ii=i0; ii<PARM_NUM; ii++) { \ 325 dp = (parmDesc[ii].max - parmDesc[ii].min)*TEN_MODEL_PARM_GRAD_EPS; \ 326 parm1[ii] = parm0[ii] + dp; \ 327 sqeForw = sqe(parm1, espec, dwiBuff, dwiMeas, knownB0); \ 328 parm1[ii] = parm0[ii] - dp; \ 329 sqeBack = sqe(parm1, espec, dwiBuff, dwiMeas, knownB0); \ 330 grad[ii] = (sqeForw - sqeBack)/(2*dp); \ 331 parm1[ii] = parm0[ii]; \ 332 if (parmDesc[ii].vec3 && 2 == parmDesc[ii].vecIdx) { \ 333 double dot, *gg; \ 334 const double *vv; \ 335 gg = grad + ii - 2; \ 336 vv = parm0 + ii - 2; \ 337 dot = ELL_3V_DOT(gg, vv); \ 338 ELL_3V_SCALE_INCR(gg, -dot, vv); \ 339 } \ 340 } \ 341 if (knownB0) { \ 342 grad[0] = 0; \ 343 } \ 344 return; \ 345 } 346 347 #define _TEN_SQE_FIT(model) \ 348 static double \ 349 sqeFit(double *parm, double *convFrac, unsigned int *itersTaken, \ 350 const tenExperSpec *espec, \ 351 double *dwiBuff, const double *dwiMeas, \ 352 const double *parmInit, int knownB0, \ 353 unsigned int minIter, unsigned int maxIter, \ 354 double convEps, int verbose) { \ 355 double testparm[PARM_NUM ? PARM_NUM : 1], \ 356 grad[PARM_NUM ? PARM_NUM : 1]; \ 357 \ 358 return _tenModelSqeFitSingle((model), \ 359 testparm, grad, \ 360 parm, convFrac, itersTaken, \ 361 espec, dwiBuff, dwiMeas, \ 362 parmInit, knownB0, \ 363 minIter, maxIter, \ 364 convEps, verbose); \ 365 } 366 367 #define _TEN_SQE_FIT_STUB \ 368 static double \ 369 sqeFit(double *parm, double *convFrac, unsigned int *itersTaken, \ 370 const tenExperSpec *espec, \ 371 double *dwiBuff, const double *dwiMeas, \ 372 const double *parmInit, int knownB0, \ 373 unsigned int minIter, unsigned int maxIter, double convEps) { \ 374 unsigned int pp; \ 375 \ 376 AIR_UNUSED(convFrac); \ 377 AIR_UNUSED(espec); \ 378 AIR_UNUSED(dwiBuff); \ 379 AIR_UNUSED(dwiMeas); \ 380 AIR_UNUSED(knownB0); \ 381 AIR_UNUSED(minIter); \ 382 AIR_UNUSED(maxIter); \ 383 AIR_UNUSED(convEps); \ 384 for (pp=0; pp<PARM_NUM; pp++) { \ 385 parm[pp] = parmInit[pp]; \ 386 } \ 387 return 0; \ 388 } 389 390 #define _TEN_NLL \ 391 static double \ 392 nll(const double *parm, const tenExperSpec *espec, \ 393 double *dwiBuff, const double *dwiMeas, \ 394 int rician, double sigma, int knownB0) { \ 395 \ 396 simulate(dwiBuff, parm, espec); \ 397 return _tenExperSpec_nll(dwiMeas, dwiBuff, espec, \ 398 rician, sigma, knownB0); \ 399 } 400 401 #define _TEN_NLL_GRAD_STUB \ 402 static void \ 403 nllGrad(double *grad, const double *parm, \ 404 const tenExperSpec *espec, \ 405 double *dwiBuff, const double *dwiMeas, \ 406 int rician, double sigma) { \ 407 \ 408 AIR_UNUSED(grad); \ 409 AIR_UNUSED(parm); \ 410 AIR_UNUSED(espec); \ 411 AIR_UNUSED(dwiBuff); \ 412 AIR_UNUSED(dwiMeas); \ 413 AIR_UNUSED(rician); \ 414 AIR_UNUSED(sigma); \ 415 return; \ 416 } 417 418 #define _TEN_NLL_FIT_STUB \ 419 static double \ 420 nllFit(double *parm, const tenExperSpec *espec, \ 421 const double *dwiMeas, const double *parmInit, \ 422 int rician, double sigma, int knownB0) { \ 423 unsigned int pp; \ 424 \ 425 AIR_UNUSED(espec); \ 426 AIR_UNUSED(dwiMeas); \ 427 AIR_UNUSED(rician); \ 428 AIR_UNUSED(sigma); \ 429 AIR_UNUSED(knownB0); \ 430 for (pp=0; pp<PARM_NUM; pp++) { \ 431 parm[pp] = parmInit[pp]; \ 432 } \ 433 return 0; \ 434 } 435 436 #define _TEN_MODEL_FIELDS \ 437 PARM_NUM, parmDesc, \ 438 simulate, \ 439 parmSprint, parmAlloc, parmRand, \ 440 parmStep, parmDist, parmCopy, parmConvert, \ 441 sqe, sqeGrad, sqeFit, \ 442 nll, nllGrad, nllFit 443 444 #ifdef __cplusplus 445 } 446 #endif 447 448 #endif /* TEN_PRIVATE_HAS_BEEN_INCLUDED */ 449