Bug Summary

File:src/ten/fiber.c
Location:line 617, column 5
Description:Value stored to 'fptsIdx' is never read

Annotated Source Code

1/*
2 Teem: Tools to process and visualize scientific data and images .
3 Copyright (C) 2013, 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#include "ten.h"
25#include "privateTen.h"
26
27#define TEN_FIBER_INCR512 512
28
29/*
30** _tenFiberProbe
31**
32** The job here is to probe at (world space) "wPos" and then set:
33** tfx->fiberTen
34** tfx->fiberEval (all 3 evals)
35** tfx->fiberEvec (all 3 eigenvectors)
36** if (tfx->stop & (1 << tenFiberStopAniso): tfx->fiberAnisoStop
37**
38** In the case of non-single-tensor tractography, we do so based on
39** ten2Which (when at the seedpoint) or
40**
41** Note that for performance reasons, a non-zero return value
42** (indicating error) and the associated use of biff, is only possible
43** if seedProbe is non-zero, the reason being that problems can be
44** detected at the seedpoint, and won't arise after the seedpoint.
45**
46** Errors from gage are indicated by *gageRet, which includes leaving
47** the domain of the volume, which is used to terminate fibers.
48**
49** Our use of tfx->seedEvec (shared with _tenFiberAlign), as well as that
50** of tfx->lastDir and tfx->lastDirSet, could stand to have further
51** debugging and documentation ...
52*/
53int
54_tenFiberProbe(tenFiberContext *tfx, int *gageRet,
55 double wPos[3], int seedProbe) {
56 static const char me[]="_tenFiberProbe";
57 double iPos[3];
58 int ret = 0;
59 double tens2[2][7];
60
61 gageShapeWtoI(tfx->gtx->shape, iPos, wPos);
62 *gageRet = gageProbe(tfx->gtx, iPos[0], iPos[1], iPos[2]);
63
64 if (tfx->verbose > 2) {
65 fprintf(stderr__stderrp, "%s(%g,%g,%g, %s): hi ----- %s\n", me,
66 iPos[0], iPos[1], iPos[2], seedProbe ? "***TRUE***" : "false",
67 tfx->useDwi ? "using DWIs" : "");
68 }
69
70 if (!tfx->useDwi) {
71 /* normal single-tensor tracking */
72 TEN_T_COPY(tfx->fiberTen, tfx->gageTen)( (tfx->fiberTen)[0] = (tfx->gageTen)[0], (tfx->fiberTen
)[1] = (tfx->gageTen)[1], (tfx->fiberTen)[2] = (tfx->
gageTen)[2], (tfx->fiberTen)[3] = (tfx->gageTen)[3], (tfx
->fiberTen)[4] = (tfx->gageTen)[4], (tfx->fiberTen)[
5] = (tfx->gageTen)[5], (tfx->fiberTen)[6] = (tfx->gageTen
)[6] );
73 ELL_3V_COPY(tfx->fiberEval, tfx->gageEval)((tfx->fiberEval)[0] = (tfx->gageEval)[0], (tfx->fiberEval
)[1] = (tfx->gageEval)[1], (tfx->fiberEval)[2] = (tfx->
gageEval)[2]);
74 ELL_3M_COPY(tfx->fiberEvec, tfx->gageEvec)((((tfx->fiberEvec)+0)[0] = ((tfx->gageEvec)+0)[0], ((tfx
->fiberEvec)+0)[1] = ((tfx->gageEvec)+0)[1], ((tfx->
fiberEvec)+0)[2] = ((tfx->gageEvec)+0)[2]), (((tfx->fiberEvec
)+3)[0] = ((tfx->gageEvec)+3)[0], ((tfx->fiberEvec)+3)[
1] = ((tfx->gageEvec)+3)[1], ((tfx->fiberEvec)+3)[2] = (
(tfx->gageEvec)+3)[2]), (((tfx->fiberEvec)+6)[0] = ((tfx
->gageEvec)+6)[0], ((tfx->fiberEvec)+6)[1] = ((tfx->
gageEvec)+6)[1], ((tfx->fiberEvec)+6)[2] = ((tfx->gageEvec
)+6)[2]));
75 if (tfx->stop & (1 << tenFiberStopAniso)) {
76 tfx->fiberAnisoStop = tfx->gageAnisoStop[0];
77 }
78 if (seedProbe) {
79 ELL_3V_COPY(tfx->seedEvec, tfx->fiberEvec)((tfx->seedEvec)[0] = (tfx->fiberEvec)[0], (tfx->seedEvec
)[1] = (tfx->fiberEvec)[1], (tfx->seedEvec)[2] = (tfx->
fiberEvec)[2]);
80 }
81 } else { /* tracking in DWIs */
82 if (tfx->verbose > 2 && seedProbe) {
83 fprintf(stderr__stderrp, "%s: fiber type = %s\n", me,
84 airEnumStr(tenDwiFiberType, tfx->fiberType));
85 }
86 switch (tfx->fiberType) {
87 double evec[2][9], eval[2][3];
88 case tenDwiFiberType1Evec0:
89 if (tfx->mframeUse) {
90 double matTmpA[9], matTmpB[9];
91 TEN_T2M(matTmpA, tfx->gageTen)( (matTmpA)[0] = (tfx->gageTen)[1], (matTmpA)[1] = (tfx->
gageTen)[2], (matTmpA)[2] = (tfx->gageTen)[3], (matTmpA)[3
] = (tfx->gageTen)[2], (matTmpA)[4] = (tfx->gageTen)[4]
, (matTmpA)[5] = (tfx->gageTen)[5], (matTmpA)[6] = (tfx->
gageTen)[3], (matTmpA)[7] = (tfx->gageTen)[5], (matTmpA)[8
] = (tfx->gageTen)[6] );
92 ELL_3M_MUL(matTmpB, tfx->mframe, matTmpA)((matTmpB)[0] = (tfx->mframe)[0]*(matTmpA)[0] + (tfx->mframe
)[1]*(matTmpA)[3] + (tfx->mframe)[2]*(matTmpA)[6], (matTmpB
)[1] = (tfx->mframe)[0]*(matTmpA)[1] + (tfx->mframe)[1]
*(matTmpA)[4] + (tfx->mframe)[2]*(matTmpA)[7], (matTmpB)[2
] = (tfx->mframe)[0]*(matTmpA)[2] + (tfx->mframe)[1]*(matTmpA
)[5] + (tfx->mframe)[2]*(matTmpA)[8], (matTmpB)[3] = (tfx->
mframe)[3]*(matTmpA)[0] + (tfx->mframe)[4]*(matTmpA)[3] + (
tfx->mframe)[5]*(matTmpA)[6], (matTmpB)[4] = (tfx->mframe
)[3]*(matTmpA)[1] + (tfx->mframe)[4]*(matTmpA)[4] + (tfx->
mframe)[5]*(matTmpA)[7], (matTmpB)[5] = (tfx->mframe)[3]*(
matTmpA)[2] + (tfx->mframe)[4]*(matTmpA)[5] + (tfx->mframe
)[5]*(matTmpA)[8], (matTmpB)[6] = (tfx->mframe)[6]*(matTmpA
)[0] + (tfx->mframe)[7]*(matTmpA)[3] + (tfx->mframe)[8]
*(matTmpA)[6], (matTmpB)[7] = (tfx->mframe)[6]*(matTmpA)[1
] + (tfx->mframe)[7]*(matTmpA)[4] + (tfx->mframe)[8]*(matTmpA
)[7], (matTmpB)[8] = (tfx->mframe)[6]*(matTmpA)[2] + (tfx->
mframe)[7]*(matTmpA)[5] + (tfx->mframe)[8]*(matTmpA)[8]);
93 ELL_3M_MUL(matTmpA, matTmpB, tfx->mframeT)((matTmpA)[0] = (matTmpB)[0]*(tfx->mframeT)[0] + (matTmpB)
[1]*(tfx->mframeT)[3] + (matTmpB)[2]*(tfx->mframeT)[6],
(matTmpA)[1] = (matTmpB)[0]*(tfx->mframeT)[1] + (matTmpB)
[1]*(tfx->mframeT)[4] + (matTmpB)[2]*(tfx->mframeT)[7],
(matTmpA)[2] = (matTmpB)[0]*(tfx->mframeT)[2] + (matTmpB)
[1]*(tfx->mframeT)[5] + (matTmpB)[2]*(tfx->mframeT)[8],
(matTmpA)[3] = (matTmpB)[3]*(tfx->mframeT)[0] + (matTmpB)
[4]*(tfx->mframeT)[3] + (matTmpB)[5]*(tfx->mframeT)[6],
(matTmpA)[4] = (matTmpB)[3]*(tfx->mframeT)[1] + (matTmpB)
[4]*(tfx->mframeT)[4] + (matTmpB)[5]*(tfx->mframeT)[7],
(matTmpA)[5] = (matTmpB)[3]*(tfx->mframeT)[2] + (matTmpB)
[4]*(tfx->mframeT)[5] + (matTmpB)[5]*(tfx->mframeT)[8],
(matTmpA)[6] = (matTmpB)[6]*(tfx->mframeT)[0] + (matTmpB)
[7]*(tfx->mframeT)[3] + (matTmpB)[8]*(tfx->mframeT)[6],
(matTmpA)[7] = (matTmpB)[6]*(tfx->mframeT)[1] + (matTmpB)
[7]*(tfx->mframeT)[4] + (matTmpB)[8]*(tfx->mframeT)[7],
(matTmpA)[8] = (matTmpB)[6]*(tfx->mframeT)[2] + (matTmpB)
[7]*(tfx->mframeT)[5] + (matTmpB)[8]*(tfx->mframeT)[8]);
94 TEN_M2T(tfx->fiberTen, matTmpA)( (tfx->fiberTen)[1] = (matTmpA)[0], (tfx->fiberTen)[2]
= ((matTmpA)[1]+(matTmpA)[3])/2.0, (tfx->fiberTen)[3] = (
(matTmpA)[2]+(matTmpA)[6])/2.0, (tfx->fiberTen)[4] = (matTmpA
)[4], (tfx->fiberTen)[5] = ((matTmpA)[5]+(matTmpA)[7])/2.0
, (tfx->fiberTen)[6] = (matTmpA)[8] );
95 tfx->fiberTen[0] = tfx->gageTen[0];
96 } else {
97 TEN_T_COPY(tfx->fiberTen, tfx->gageTen)( (tfx->fiberTen)[0] = (tfx->gageTen)[0], (tfx->fiberTen
)[1] = (tfx->gageTen)[1], (tfx->fiberTen)[2] = (tfx->
gageTen)[2], (tfx->fiberTen)[3] = (tfx->gageTen)[3], (tfx
->fiberTen)[4] = (tfx->gageTen)[4], (tfx->fiberTen)[
5] = (tfx->gageTen)[5], (tfx->fiberTen)[6] = (tfx->gageTen
)[6] );
98 }
99 tenEigensolve_d(tfx->fiberEval, tfx->fiberEvec, tfx->fiberTen);
100 if (tfx->stop & (1 << tenFiberStopAniso)) {
101 double tmp;
102 tmp = tenAnisoTen_d(tfx->fiberTen, tfx->anisoStopType);
103 tfx->fiberAnisoStop = AIR_CLAMP(0, tmp, 1)((tmp) < (0) ? (0) : ((tmp) > (1) ? (1) : (tmp)));
104 }
105 if (seedProbe) {
106 ELL_3V_COPY(tfx->seedEvec, tfx->fiberEvec)((tfx->seedEvec)[0] = (tfx->fiberEvec)[0], (tfx->seedEvec
)[1] = (tfx->fiberEvec)[1], (tfx->seedEvec)[2] = (tfx->
fiberEvec)[2]);
107 }
108 break;
109 case tenDwiFiberType2Evec0:
110 /* Estimate principal diffusion direction of each tensor */
111 if (tfx->mframeUse) {
112 /* Transform both the tensors */
113 double matTmpA[9], matTmpB[9];
114
115 TEN_T2M(matTmpA, tfx->gageTen2 + 0)( (matTmpA)[0] = (tfx->gageTen2 + 0)[1], (matTmpA)[1] = (tfx
->gageTen2 + 0)[2], (matTmpA)[2] = (tfx->gageTen2 + 0)[
3], (matTmpA)[3] = (tfx->gageTen2 + 0)[2], (matTmpA)[4] = (
tfx->gageTen2 + 0)[4], (matTmpA)[5] = (tfx->gageTen2 + 0
)[5], (matTmpA)[6] = (tfx->gageTen2 + 0)[3], (matTmpA)[7] =
(tfx->gageTen2 + 0)[5], (matTmpA)[8] = (tfx->gageTen2 +
0)[6] );
116 ELL_3M_MUL(matTmpB, tfx->mframe, matTmpA)((matTmpB)[0] = (tfx->mframe)[0]*(matTmpA)[0] + (tfx->mframe
)[1]*(matTmpA)[3] + (tfx->mframe)[2]*(matTmpA)[6], (matTmpB
)[1] = (tfx->mframe)[0]*(matTmpA)[1] + (tfx->mframe)[1]
*(matTmpA)[4] + (tfx->mframe)[2]*(matTmpA)[7], (matTmpB)[2
] = (tfx->mframe)[0]*(matTmpA)[2] + (tfx->mframe)[1]*(matTmpA
)[5] + (tfx->mframe)[2]*(matTmpA)[8], (matTmpB)[3] = (tfx->
mframe)[3]*(matTmpA)[0] + (tfx->mframe)[4]*(matTmpA)[3] + (
tfx->mframe)[5]*(matTmpA)[6], (matTmpB)[4] = (tfx->mframe
)[3]*(matTmpA)[1] + (tfx->mframe)[4]*(matTmpA)[4] + (tfx->
mframe)[5]*(matTmpA)[7], (matTmpB)[5] = (tfx->mframe)[3]*(
matTmpA)[2] + (tfx->mframe)[4]*(matTmpA)[5] + (tfx->mframe
)[5]*(matTmpA)[8], (matTmpB)[6] = (tfx->mframe)[6]*(matTmpA
)[0] + (tfx->mframe)[7]*(matTmpA)[3] + (tfx->mframe)[8]
*(matTmpA)[6], (matTmpB)[7] = (tfx->mframe)[6]*(matTmpA)[1
] + (tfx->mframe)[7]*(matTmpA)[4] + (tfx->mframe)[8]*(matTmpA
)[7], (matTmpB)[8] = (tfx->mframe)[6]*(matTmpA)[2] + (tfx->
mframe)[7]*(matTmpA)[5] + (tfx->mframe)[8]*(matTmpA)[8]);
117 ELL_3M_MUL(matTmpA, matTmpB, tfx->mframeT)((matTmpA)[0] = (matTmpB)[0]*(tfx->mframeT)[0] + (matTmpB)
[1]*(tfx->mframeT)[3] + (matTmpB)[2]*(tfx->mframeT)[6],
(matTmpA)[1] = (matTmpB)[0]*(tfx->mframeT)[1] + (matTmpB)
[1]*(tfx->mframeT)[4] + (matTmpB)[2]*(tfx->mframeT)[7],
(matTmpA)[2] = (matTmpB)[0]*(tfx->mframeT)[2] + (matTmpB)
[1]*(tfx->mframeT)[5] + (matTmpB)[2]*(tfx->mframeT)[8],
(matTmpA)[3] = (matTmpB)[3]*(tfx->mframeT)[0] + (matTmpB)
[4]*(tfx->mframeT)[3] + (matTmpB)[5]*(tfx->mframeT)[6],
(matTmpA)[4] = (matTmpB)[3]*(tfx->mframeT)[1] + (matTmpB)
[4]*(tfx->mframeT)[4] + (matTmpB)[5]*(tfx->mframeT)[7],
(matTmpA)[5] = (matTmpB)[3]*(tfx->mframeT)[2] + (matTmpB)
[4]*(tfx->mframeT)[5] + (matTmpB)[5]*(tfx->mframeT)[8],
(matTmpA)[6] = (matTmpB)[6]*(tfx->mframeT)[0] + (matTmpB)
[7]*(tfx->mframeT)[3] + (matTmpB)[8]*(tfx->mframeT)[6],
(matTmpA)[7] = (matTmpB)[6]*(tfx->mframeT)[1] + (matTmpB)
[7]*(tfx->mframeT)[4] + (matTmpB)[8]*(tfx->mframeT)[7],
(matTmpA)[8] = (matTmpB)[6]*(tfx->mframeT)[2] + (matTmpB)
[7]*(tfx->mframeT)[5] + (matTmpB)[8]*(tfx->mframeT)[8]);
118 TEN_M2T(tens2[0], matTmpA)( (tens2[0])[1] = (matTmpA)[0], (tens2[0])[2] = ((matTmpA)[1]
+(matTmpA)[3])/2.0, (tens2[0])[3] = ((matTmpA)[2]+(matTmpA)[6
])/2.0, (tens2[0])[4] = (matTmpA)[4], (tens2[0])[5] = ((matTmpA
)[5]+(matTmpA)[7])/2.0, (tens2[0])[6] = (matTmpA)[8] );
119 /* new eigen values and vectors */
120 tenEigensolve_d(eval[0], evec[0], tens2[0]);
121
122 TEN_T2M(matTmpA, tfx->gageTen2 + 7)( (matTmpA)[0] = (tfx->gageTen2 + 7)[1], (matTmpA)[1] = (tfx
->gageTen2 + 7)[2], (matTmpA)[2] = (tfx->gageTen2 + 7)[
3], (matTmpA)[3] = (tfx->gageTen2 + 7)[2], (matTmpA)[4] = (
tfx->gageTen2 + 7)[4], (matTmpA)[5] = (tfx->gageTen2 + 7
)[5], (matTmpA)[6] = (tfx->gageTen2 + 7)[3], (matTmpA)[7] =
(tfx->gageTen2 + 7)[5], (matTmpA)[8] = (tfx->gageTen2 +
7)[6] );
123 ELL_3M_MUL(matTmpB, tfx->mframe, matTmpA)((matTmpB)[0] = (tfx->mframe)[0]*(matTmpA)[0] + (tfx->mframe
)[1]*(matTmpA)[3] + (tfx->mframe)[2]*(matTmpA)[6], (matTmpB
)[1] = (tfx->mframe)[0]*(matTmpA)[1] + (tfx->mframe)[1]
*(matTmpA)[4] + (tfx->mframe)[2]*(matTmpA)[7], (matTmpB)[2
] = (tfx->mframe)[0]*(matTmpA)[2] + (tfx->mframe)[1]*(matTmpA
)[5] + (tfx->mframe)[2]*(matTmpA)[8], (matTmpB)[3] = (tfx->
mframe)[3]*(matTmpA)[0] + (tfx->mframe)[4]*(matTmpA)[3] + (
tfx->mframe)[5]*(matTmpA)[6], (matTmpB)[4] = (tfx->mframe
)[3]*(matTmpA)[1] + (tfx->mframe)[4]*(matTmpA)[4] + (tfx->
mframe)[5]*(matTmpA)[7], (matTmpB)[5] = (tfx->mframe)[3]*(
matTmpA)[2] + (tfx->mframe)[4]*(matTmpA)[5] + (tfx->mframe
)[5]*(matTmpA)[8], (matTmpB)[6] = (tfx->mframe)[6]*(matTmpA
)[0] + (tfx->mframe)[7]*(matTmpA)[3] + (tfx->mframe)[8]
*(matTmpA)[6], (matTmpB)[7] = (tfx->mframe)[6]*(matTmpA)[1
] + (tfx->mframe)[7]*(matTmpA)[4] + (tfx->mframe)[8]*(matTmpA
)[7], (matTmpB)[8] = (tfx->mframe)[6]*(matTmpA)[2] + (tfx->
mframe)[7]*(matTmpA)[5] + (tfx->mframe)[8]*(matTmpA)[8]);
124 ELL_3M_MUL(matTmpA, matTmpB, tfx->mframeT)((matTmpA)[0] = (matTmpB)[0]*(tfx->mframeT)[0] + (matTmpB)
[1]*(tfx->mframeT)[3] + (matTmpB)[2]*(tfx->mframeT)[6],
(matTmpA)[1] = (matTmpB)[0]*(tfx->mframeT)[1] + (matTmpB)
[1]*(tfx->mframeT)[4] + (matTmpB)[2]*(tfx->mframeT)[7],
(matTmpA)[2] = (matTmpB)[0]*(tfx->mframeT)[2] + (matTmpB)
[1]*(tfx->mframeT)[5] + (matTmpB)[2]*(tfx->mframeT)[8],
(matTmpA)[3] = (matTmpB)[3]*(tfx->mframeT)[0] + (matTmpB)
[4]*(tfx->mframeT)[3] + (matTmpB)[5]*(tfx->mframeT)[6],
(matTmpA)[4] = (matTmpB)[3]*(tfx->mframeT)[1] + (matTmpB)
[4]*(tfx->mframeT)[4] + (matTmpB)[5]*(tfx->mframeT)[7],
(matTmpA)[5] = (matTmpB)[3]*(tfx->mframeT)[2] + (matTmpB)
[4]*(tfx->mframeT)[5] + (matTmpB)[5]*(tfx->mframeT)[8],
(matTmpA)[6] = (matTmpB)[6]*(tfx->mframeT)[0] + (matTmpB)
[7]*(tfx->mframeT)[3] + (matTmpB)[8]*(tfx->mframeT)[6],
(matTmpA)[7] = (matTmpB)[6]*(tfx->mframeT)[1] + (matTmpB)
[7]*(tfx->mframeT)[4] + (matTmpB)[8]*(tfx->mframeT)[7],
(matTmpA)[8] = (matTmpB)[6]*(tfx->mframeT)[2] + (matTmpB)
[7]*(tfx->mframeT)[5] + (matTmpB)[8]*(tfx->mframeT)[8]);
125 TEN_M2T(tens2[1], matTmpA)( (tens2[1])[1] = (matTmpA)[0], (tens2[1])[2] = ((matTmpA)[1]
+(matTmpA)[3])/2.0, (tens2[1])[3] = ((matTmpA)[2]+(matTmpA)[6
])/2.0, (tens2[1])[4] = (matTmpA)[4], (tens2[1])[5] = ((matTmpA
)[5]+(matTmpA)[7])/2.0, (tens2[1])[6] = (matTmpA)[8] );
126 tenEigensolve_d(eval[1], evec[1], tens2[1]);
127 } else {
128 tenEigensolve_d(eval[0], evec[0], tfx->gageTen2 + 0);
129 tenEigensolve_d(eval[1], evec[1], tfx->gageTen2 + 7);
130 }
131
132 /* set ten2Use */
133 if (seedProbe) { /* we're on the *very* 1st probe per tract,
134 at the seed pt */
135 ELL_3V_COPY(tfx->seedEvec, evec[tfx->ten2Which])((tfx->seedEvec)[0] = (evec[tfx->ten2Which])[0], (tfx->
seedEvec)[1] = (evec[tfx->ten2Which])[1], (tfx->seedEvec
)[2] = (evec[tfx->ten2Which])[2]);
136 tfx->ten2Use = tfx->ten2Which;
137 if (tfx->verbose > 2) {
138 fprintf(stderr__stderrp, "%s: ** ten2Use == ten2Which == %d\n", me,
139 tfx->ten2Use);
140 }
141 } else {
142 double *lastVec, dot[2];
143
144 if (!tfx->lastDirSet) {
145 /* we're on some probe of the first step */
146 lastVec = tfx->seedEvec;
147 } else {
148 /* we're past the first step */
149 /* Arish says: "Bug len has not been initialized and don't think
150 its needed". The first part is not a problem; "len" is in the
151 *output* argument of ELL_3V_NORM. The second part seems to be
152 true, even though Gordon can't currently see why! */
153 /* ELL_3V_NORM(tfx->lastDir, tfx->lastDir, len); */
154 lastVec = tfx->lastDir;
155 }
156 dot[0] = ELL_3V_DOT(lastVec, evec[0])((lastVec)[0]*(evec[0])[0] + (lastVec)[1]*(evec[0])[1] + (lastVec
)[2]*(evec[0])[2]);
157 dot[1] = ELL_3V_DOT(lastVec, evec[1])((lastVec)[0]*(evec[1])[0] + (lastVec)[1]*(evec[1])[1] + (lastVec
)[2]*(evec[1])[2]);
158 if (dot[0] < 0) {
159 dot[0] *= -1;
160 ELL_3M_SCALE(evec[0], -1, evec[0])((((evec[0])+0)[0] = ((-1))*((evec[0])+0)[0], ((evec[0])+0)[1
] = ((-1))*((evec[0])+0)[1], ((evec[0])+0)[2] = ((-1))*((evec
[0])+0)[2]), (((evec[0])+3)[0] = ((-1))*((evec[0])+3)[0], ((evec
[0])+3)[1] = ((-1))*((evec[0])+3)[1], ((evec[0])+3)[2] = ((-1
))*((evec[0])+3)[2]), (((evec[0])+6)[0] = ((-1))*((evec[0])+6
)[0], ((evec[0])+6)[1] = ((-1))*((evec[0])+6)[1], ((evec[0])+
6)[2] = ((-1))*((evec[0])+6)[2]));
161 }
162 if (dot[1] < 0) {
163 dot[1] *= -1;
164 ELL_3M_SCALE(evec[1], -1, evec[1])((((evec[1])+0)[0] = ((-1))*((evec[1])+0)[0], ((evec[1])+0)[1
] = ((-1))*((evec[1])+0)[1], ((evec[1])+0)[2] = ((-1))*((evec
[1])+0)[2]), (((evec[1])+3)[0] = ((-1))*((evec[1])+3)[0], ((evec
[1])+3)[1] = ((-1))*((evec[1])+3)[1], ((evec[1])+3)[2] = ((-1
))*((evec[1])+3)[2]), (((evec[1])+6)[0] = ((-1))*((evec[1])+6
)[0], ((evec[1])+6)[1] = ((-1))*((evec[1])+6)[1], ((evec[1])+
6)[2] = ((-1))*((evec[1])+6)[2]));
165 }
166 tfx->ten2Use = (dot[0] > dot[1]) ? 0 : 1;
167 if (tfx->verbose > 2) {
168 fprintf(stderr__stderrp, "%s(%g,%g,%g): dot[0,1] = %f, %f -> use %u\n",
169 me, wPos[0], wPos[1], wPos[2], dot[0], dot[1],
170 tfx->ten2Use );
171 }
172 }
173
174 /* based on ten2Use, set the rest of the needed fields */
175 if (tfx->mframeUse) {
176 TEN_T_COPY(tfx->fiberTen, tens2[tfx->ten2Use])( (tfx->fiberTen)[0] = (tens2[tfx->ten2Use])[0], (tfx->
fiberTen)[1] = (tens2[tfx->ten2Use])[1], (tfx->fiberTen
)[2] = (tens2[tfx->ten2Use])[2], (tfx->fiberTen)[3] = (
tens2[tfx->ten2Use])[3], (tfx->fiberTen)[4] = (tens2[tfx
->ten2Use])[4], (tfx->fiberTen)[5] = (tens2[tfx->ten2Use
])[5], (tfx->fiberTen)[6] = (tens2[tfx->ten2Use])[6] );
177 } else {
178 TEN_T_COPY(tfx->fiberTen, tfx->gageTen2 + 7*(tfx->ten2Use))( (tfx->fiberTen)[0] = (tfx->gageTen2 + 7*(tfx->ten2Use
))[0], (tfx->fiberTen)[1] = (tfx->gageTen2 + 7*(tfx->
ten2Use))[1], (tfx->fiberTen)[2] = (tfx->gageTen2 + 7*(
tfx->ten2Use))[2], (tfx->fiberTen)[3] = (tfx->gageTen2
+ 7*(tfx->ten2Use))[3], (tfx->fiberTen)[4] = (tfx->
gageTen2 + 7*(tfx->ten2Use))[4], (tfx->fiberTen)[5] = (
tfx->gageTen2 + 7*(tfx->ten2Use))[5], (tfx->fiberTen
)[6] = (tfx->gageTen2 + 7*(tfx->ten2Use))[6] );
179 }
180 tfx->fiberTen[0] = tfx->gageTen2[0]; /* copy confidence */
181 ELL_3V_COPY(tfx->fiberEval, eval[tfx->ten2Use])((tfx->fiberEval)[0] = (eval[tfx->ten2Use])[0], (tfx->
fiberEval)[1] = (eval[tfx->ten2Use])[1], (tfx->fiberEval
)[2] = (eval[tfx->ten2Use])[2]);
182 ELL_3M_COPY(tfx->fiberEvec, evec[tfx->ten2Use])((((tfx->fiberEvec)+0)[0] = ((evec[tfx->ten2Use])+0)[0]
, ((tfx->fiberEvec)+0)[1] = ((evec[tfx->ten2Use])+0)[1]
, ((tfx->fiberEvec)+0)[2] = ((evec[tfx->ten2Use])+0)[2]
), (((tfx->fiberEvec)+3)[0] = ((evec[tfx->ten2Use])+3)[
0], ((tfx->fiberEvec)+3)[1] = ((evec[tfx->ten2Use])+3)[
1], ((tfx->fiberEvec)+3)[2] = ((evec[tfx->ten2Use])+3)[
2]), (((tfx->fiberEvec)+6)[0] = ((evec[tfx->ten2Use])+6
)[0], ((tfx->fiberEvec)+6)[1] = ((evec[tfx->ten2Use])+6
)[1], ((tfx->fiberEvec)+6)[2] = ((evec[tfx->ten2Use])+6
)[2]));
183 if (tfx->stop & (1 << tenFiberStopAniso)) {
184 double tmp;
185 tmp = tenAnisoEval_d(tfx->fiberEval, tfx->anisoStopType);
186 tfx->fiberAnisoStop = AIR_CLAMP(0, tmp, 1)((tmp) < (0) ? (0) : ((tmp) > (1) ? (1) : (tmp)));
187 /* HEY: what about speed? */
188 } else {
189 tfx->fiberAnisoStop = AIR_NAN(airFloatQNaN.f);
190 }
191 break;
192 default:
193 biffAddf(TENtenBiffKey, "%s: %s %s (%d) unimplemented!", me,
194 tenDwiFiberType->name,
195 airEnumStr(tenDwiFiberType, tfx->fiberType), tfx->fiberType);
196 ret = 1;
197 } /* switch (tfx->fiberType) */
198 }
199 if (tfx->verbose > 2) {
200 fprintf(stderr__stderrp, "%s: fiberEvec = %g %g %g\n", me,
201 tfx->fiberEvec[0], tfx->fiberEvec[1], tfx->fiberEvec[2]);
202 }
203
204 return ret;
205}
206
207int
208_tenFiberStopCheck(tenFiberContext *tfx) {
209 static const char me[]="_tenFiberStopCheck";
210
211 if (tfx->numSteps[tfx->halfIdx] >= TEN_FIBER_NUM_STEPS_MAX10240) {
212 fprintf(stderr__stderrp, "%s: numSteps[%d] exceeded sanity check value of %d!!\n",
213 me, tfx->halfIdx, TEN_FIBER_NUM_STEPS_MAX10240);
214 fprintf(stderr__stderrp, "%s: Check fiber termination conditions, or recompile "
215 "with a larger value for TEN_FIBER_NUM_STEPS_MAX\n", me);
216 return tenFiberStopNumSteps;
217 }
218 if (tfx->stop & (1 << tenFiberStopConfidence)) {
219 if (tfx->fiberTen[0] < tfx->confThresh) {
220 return tenFiberStopConfidence;
221 }
222 }
223 if (tfx->stop & (1 << tenFiberStopRadius)) {
224 if (tfx->radius < tfx->minRadius) {
225 return tenFiberStopRadius;
226 }
227 }
228 if (tfx->stop & (1 << tenFiberStopAniso)) {
229 if (tfx->fiberAnisoStop < tfx->anisoThresh) {
230 return tenFiberStopAniso;
231 }
232 }
233 if (tfx->stop & (1 << tenFiberStopNumSteps)) {
234 if (tfx->numSteps[tfx->halfIdx] > tfx->maxNumSteps) {
235 return tenFiberStopNumSteps;
236 }
237 }
238 if (tfx->stop & (1 << tenFiberStopLength)) {
239 if (tfx->halfLen[tfx->halfIdx] >= tfx->maxHalfLen) {
240 return tenFiberStopLength;
241 }
242 }
243 if (tfx->useDwi
244 && tfx->stop & (1 << tenFiberStopFraction)
245 && tfx->gageTen2) { /* not all DWI fiber types use gageTen2 */
246 double fracUse;
247 fracUse = (0 == tfx->ten2Use
248 ? tfx->gageTen2[7]
249 : 1 - tfx->gageTen2[7]);
250 if (fracUse < tfx->minFraction) {
251 return tenFiberStopFraction;
252 }
253 }
254 return 0;
255}
256
257void
258_tenFiberAlign(tenFiberContext *tfx, double vec[3]) {
259 static const char me[]="_tenFiberAlign";
260 double scale, dot;
261
262 if (tfx->verbose > 2) {
263 fprintf(stderr__stderrp, "%s: hi %s (lds %d):\t%g %g %g\n", me,
264 !tfx->lastDirSet ? "**" : " ",
265 tfx->lastDirSet, vec[0], vec[1], vec[2]);
266 }
267 if (!(tfx->lastDirSet)) {
268 dot = ELL_3V_DOT(tfx->seedEvec, vec)((tfx->seedEvec)[0]*(vec)[0] + (tfx->seedEvec)[1]*(vec)
[1] + (tfx->seedEvec)[2]*(vec)[2]);
269 /* this is the first step (or one of the intermediate steps
270 for RK) in this fiber half; 1st half follows the
271 eigenvector determined at seed point, 2nd goes opposite */
272 if (tfx->verbose > 2) {
273 fprintf(stderr__stderrp, "!%s: dir=%d, dot=%g\n", me, tfx->halfIdx, dot);
274 }
275 if (!tfx->halfIdx) {
276 /* 1st half */
277 scale = dot < 0 ? -1 : 1;
278 } else {
279 /* 2nd half */
280 scale = dot > 0 ? -1 : 1;
281 }
282 } else {
283 dot = ELL_3V_DOT(tfx->lastDir, vec)((tfx->lastDir)[0]*(vec)[0] + (tfx->lastDir)[1]*(vec)[1
] + (tfx->lastDir)[2]*(vec)[2]);
284 /* we have some history in this fiber half */
285 scale = dot < 0 ? -1 : 1;
286 }
287 ELL_3V_SCALE(vec, scale, vec)((vec)[0] = (scale)*(vec)[0], (vec)[1] = (scale)*(vec)[1], (vec
)[2] = (scale)*(vec)[2]);
288 if (tfx->verbose > 2) {
289 fprintf(stderr__stderrp, "!%s: scl = %g -> \t%g %g %g\n",
290 me, scale, vec[0], vec[1], vec[2]);
291 }
292 return;
293}
294
295/*
296** parm[0]: lerp between 1 and the stuff below
297** parm[1]: "t": (parm[1],0) is control point between (0,0) and (1,1)
298** parm[2]: "d": parabolic blend between parm[1]-parm[2] and parm[1]+parm[2]
299*/
300void
301_tenFiberAnisoSpeed(double *step, double xx, double parm[3]) {
302 double aa, dd, tt, yy;
303
304 tt = parm[1];
305 dd = parm[2];
306 aa = 1.0/(DBL_EPSILON2.2204460492503131e-16 + 4*dd*(1.0-tt));
307 yy = xx - tt + dd;
308 xx = (xx < tt - dd
309 ? 0
310 : (xx < tt + dd
311 ? aa*yy*yy
312 : (xx - tt)/(1 - tt)));
313 xx = AIR_LERP(parm[0], 1, xx)((parm[0])*((xx) - (1)) + (1));
314 ELL_3V_SCALE(step, xx, step)((step)[0] = (xx)*(step)[0], (step)[1] = (xx)*(step)[1], (step
)[2] = (xx)*(step)[2]);
315}
316
317/*
318** -------------------------------------------------------------------
319** -------------------------------------------------------------------
320** The _tenFiberStep_* routines are responsible for putting a step into
321** the given step[] vector. Without anisoStepSize, this should be
322** UNIT LENGTH, with anisoStepSize, its scaled by that anisotropy measure
323*/
324void
325_tenFiberStep_Evec(tenFiberContext *tfx, double step[3]) {
326
327 /* fiberEvec points to the correct gage answer based on fiberType */
328 ELL_3V_COPY(step, tfx->fiberEvec + 3*0)((step)[0] = (tfx->fiberEvec + 3*0)[0], (step)[1] = (tfx->
fiberEvec + 3*0)[1], (step)[2] = (tfx->fiberEvec + 3*0)[2]
);
329 _tenFiberAlign(tfx, step);
330 if (tfx->anisoSpeedType) {
331 _tenFiberAnisoSpeed(step, tfx->fiberAnisoSpeed,
332 tfx->anisoSpeedFunc);
333 }
334}
335
336void
337_tenFiberStep_TensorLine(tenFiberContext *tfx, double step[3]) {
338 double cl, evec0[3], vout[3], vin[3], len;
339
340 ELL_3V_COPY(evec0, tfx->fiberEvec + 3*0)((evec0)[0] = (tfx->fiberEvec + 3*0)[0], (evec0)[1] = (tfx
->fiberEvec + 3*0)[1], (evec0)[2] = (tfx->fiberEvec + 3
*0)[2]);
341 _tenFiberAlign(tfx, evec0);
342
343 if (tfx->lastDirSet) {
344 ELL_3V_COPY(vin, tfx->lastDir)((vin)[0] = (tfx->lastDir)[0], (vin)[1] = (tfx->lastDir
)[1], (vin)[2] = (tfx->lastDir)[2]);
345 TEN_T3V_MUL(vout, tfx->fiberTen, tfx->lastDir)( (vout)[0] = (tfx->fiberTen)[1]*(tfx->lastDir)[0] + (tfx
->fiberTen)[2]*(tfx->lastDir)[1] + (tfx->fiberTen)[3
]*(tfx->lastDir)[2], (vout)[1] = (tfx->fiberTen)[2]*(tfx
->lastDir)[0] + (tfx->fiberTen)[4]*(tfx->lastDir)[1]
+ (tfx->fiberTen)[5]*(tfx->lastDir)[2], (vout)[2] = (tfx
->fiberTen)[3]*(tfx->lastDir)[0] + (tfx->fiberTen)[5
]*(tfx->lastDir)[1] + (tfx->fiberTen)[6]*(tfx->lastDir
)[2]);
346 ELL_3V_NORM(vout, vout, len)(len = (sqrt((((vout))[0]*((vout))[0] + ((vout))[1]*((vout))[
1] + ((vout))[2]*((vout))[2]))), ((vout)[0] = (1.0/len)*(vout
)[0], (vout)[1] = (1.0/len)*(vout)[1], (vout)[2] = (1.0/len)*
(vout)[2]));
347 _tenFiberAlign(tfx, vout); /* HEY: is this needed? */
348 } else {
349 ELL_3V_COPY(vin, evec0)((vin)[0] = (evec0)[0], (vin)[1] = (evec0)[1], (vin)[2] = (evec0
)[2]);
350 ELL_3V_COPY(vout, evec0)((vout)[0] = (evec0)[0], (vout)[1] = (evec0)[1], (vout)[2] = (
evec0)[2]);
351 }
352
353 /* HEY: should be using one of the tenAnisoEval[] functions */
354 cl = (tfx->fiberEval[0] - tfx->fiberEval[1])/(tfx->fiberEval[0] + 0.00001);
355
356 ELL_3V_SCALE_ADD3(step,((step)[0] = (cl)*(evec0)[0] + ((1-cl)*(1-tfx->wPunct))*(vin
)[0] + ((1-cl)*tfx->wPunct)*(vout)[0], (step)[1] = (cl)*(evec0
)[1] + ((1-cl)*(1-tfx->wPunct))*(vin)[1] + ((1-cl)*tfx->
wPunct)*(vout)[1], (step)[2] = (cl)*(evec0)[2] + ((1-cl)*(1-tfx
->wPunct))*(vin)[2] + ((1-cl)*tfx->wPunct)*(vout)[2])
357 cl, evec0,((step)[0] = (cl)*(evec0)[0] + ((1-cl)*(1-tfx->wPunct))*(vin
)[0] + ((1-cl)*tfx->wPunct)*(vout)[0], (step)[1] = (cl)*(evec0
)[1] + ((1-cl)*(1-tfx->wPunct))*(vin)[1] + ((1-cl)*tfx->
wPunct)*(vout)[1], (step)[2] = (cl)*(evec0)[2] + ((1-cl)*(1-tfx
->wPunct))*(vin)[2] + ((1-cl)*tfx->wPunct)*(vout)[2])
358 (1-cl)*(1-tfx->wPunct), vin,((step)[0] = (cl)*(evec0)[0] + ((1-cl)*(1-tfx->wPunct))*(vin
)[0] + ((1-cl)*tfx->wPunct)*(vout)[0], (step)[1] = (cl)*(evec0
)[1] + ((1-cl)*(1-tfx->wPunct))*(vin)[1] + ((1-cl)*tfx->
wPunct)*(vout)[1], (step)[2] = (cl)*(evec0)[2] + ((1-cl)*(1-tfx
->wPunct))*(vin)[2] + ((1-cl)*tfx->wPunct)*(vout)[2])
359 (1-cl)*tfx->wPunct, vout)((step)[0] = (cl)*(evec0)[0] + ((1-cl)*(1-tfx->wPunct))*(vin
)[0] + ((1-cl)*tfx->wPunct)*(vout)[0], (step)[1] = (cl)*(evec0
)[1] + ((1-cl)*(1-tfx->wPunct))*(vin)[1] + ((1-cl)*tfx->
wPunct)*(vout)[1], (step)[2] = (cl)*(evec0)[2] + ((1-cl)*(1-tfx
->wPunct))*(vin)[2] + ((1-cl)*tfx->wPunct)*(vout)[2]);
360 /* _tenFiberAlign(tfx, step); */
361 ELL_3V_NORM(step, step, len)(len = (sqrt((((step))[0]*((step))[0] + ((step))[1]*((step))[
1] + ((step))[2]*((step))[2]))), ((step)[0] = (1.0/len)*(step
)[0], (step)[1] = (1.0/len)*(step)[1], (step)[2] = (1.0/len)*
(step)[2]));
362 if (tfx->anisoSpeedType) {
363 _tenFiberAnisoSpeed(step, tfx->fiberAnisoSpeed,
364 tfx->anisoSpeedFunc);
365 }
366}
367
368void
369_tenFiberStep_PureLine(tenFiberContext *tfx, double step[3]) {
370 static const char me[]="_tenFiberStep_PureLine";
371
372 AIR_UNUSED(tfx)(void)(tfx);
373 AIR_UNUSED(step)(void)(step);
374 fprintf(stderr__stderrp, "%s: sorry, unimplemented!\n", me);
375}
376
377void
378_tenFiberStep_Zhukov(tenFiberContext *tfx, double step[3]) {
379 static const char me[]="_tenFiberStep_Zhukov";
380
381 AIR_UNUSED(tfx)(void)(tfx);
382 AIR_UNUSED(step)(void)(step);
383 fprintf(stderr__stderrp, "%s: sorry, unimplemented!\n", me);
384}
385
386void (*
387_tenFiberStep[TEN_FIBER_TYPE_MAX6+1])(tenFiberContext *, double *) = {
388 NULL((void*)0),
389 _tenFiberStep_Evec,
390 _tenFiberStep_Evec,
391 _tenFiberStep_Evec,
392 _tenFiberStep_TensorLine,
393 _tenFiberStep_PureLine,
394 _tenFiberStep_Zhukov
395};
396
397/*
398** -------------------------------------------------------------------
399** -------------------------------------------------------------------
400** The _tenFiberIntegrate_* routines must assume that
401** _tenFiberProbe(tfx, tfx->wPos, AIR_FALSE) has just been called
402*/
403
404int
405_tenFiberIntegrate_Euler(tenFiberContext *tfx, double forwDir[3]) {
406
407 _tenFiberStep[tfx->fiberType](tfx, forwDir);
408 ELL_3V_SCALE(forwDir, tfx->stepSize, forwDir)((forwDir)[0] = (tfx->stepSize)*(forwDir)[0], (forwDir)[1]
= (tfx->stepSize)*(forwDir)[1], (forwDir)[2] = (tfx->stepSize
)*(forwDir)[2]);
409 return 0;
410}
411
412int
413_tenFiberIntegrate_Midpoint(tenFiberContext *tfx, double forwDir[3]) {
414 double loc[3], half[3];
415 int gret;
416
417 _tenFiberStep[tfx->fiberType](tfx, half);
418 ELL_3V_SCALE_ADD2(loc, 1, tfx->wPos, 0.5*tfx->stepSize, half)((loc)[0] = (1)*(tfx->wPos)[0] + (0.5*tfx->stepSize)*(half
)[0], (loc)[1] = (1)*(tfx->wPos)[1] + (0.5*tfx->stepSize
)*(half)[1], (loc)[2] = (1)*(tfx->wPos)[2] + (0.5*tfx->
stepSize)*(half)[2]);
419 _tenFiberProbe(tfx, &gret, loc, AIR_FALSE0); if (gret) return 1;
420 _tenFiberStep[tfx->fiberType](tfx, forwDir);
421 ELL_3V_SCALE(forwDir, tfx->stepSize, forwDir)((forwDir)[0] = (tfx->stepSize)*(forwDir)[0], (forwDir)[1]
= (tfx->stepSize)*(forwDir)[1], (forwDir)[2] = (tfx->stepSize
)*(forwDir)[2]);
422 return 0;
423}
424
425int
426_tenFiberIntegrate_RK4(tenFiberContext *tfx, double forwDir[3]) {
427 double loc[3], k1[3], k2[3], k3[3], k4[3], c1, c2, c3, c4, h;
428 int gret;
429
430 h = tfx->stepSize;
431 c1 = h/6.0; c2 = h/3.0; c3 = h/3.0; c4 = h/6.0;
432
433 _tenFiberStep[tfx->fiberType](tfx, k1);
434 ELL_3V_SCALE_ADD2(loc, 1, tfx->wPos, 0.5*h, k1)((loc)[0] = (1)*(tfx->wPos)[0] + (0.5*h)*(k1)[0], (loc)[1]
= (1)*(tfx->wPos)[1] + (0.5*h)*(k1)[1], (loc)[2] = (1)*(tfx
->wPos)[2] + (0.5*h)*(k1)[2]);
435 _tenFiberProbe(tfx, &gret, loc, AIR_FALSE0); if (gret) return 1;
436 _tenFiberStep[tfx->fiberType](tfx, k2);
437 ELL_3V_SCALE_ADD2(loc, 1, tfx->wPos, 0.5*h, k2)((loc)[0] = (1)*(tfx->wPos)[0] + (0.5*h)*(k2)[0], (loc)[1]
= (1)*(tfx->wPos)[1] + (0.5*h)*(k2)[1], (loc)[2] = (1)*(tfx
->wPos)[2] + (0.5*h)*(k2)[2]);
438 _tenFiberProbe(tfx, &gret, loc, AIR_FALSE0); if (gret) return 1;
439 _tenFiberStep[tfx->fiberType](tfx, k3);
440 ELL_3V_SCALE_ADD2(loc, 1, tfx->wPos, h, k3)((loc)[0] = (1)*(tfx->wPos)[0] + (h)*(k3)[0], (loc)[1] = (
1)*(tfx->wPos)[1] + (h)*(k3)[1], (loc)[2] = (1)*(tfx->wPos
)[2] + (h)*(k3)[2]);
441 _tenFiberProbe(tfx, &gret, loc, AIR_FALSE0); if (gret) return 1;
442 _tenFiberStep[tfx->fiberType](tfx, k4);
443
444 ELL_3V_SET(forwDir,((forwDir)[0] = (c1*k1[0] + c2*k2[0] + c3*k3[0] + c4*k4[0]), (
forwDir)[1] = (c1*k1[1] + c2*k2[1] + c3*k3[1] + c4*k4[1]), (forwDir
)[2] = (c1*k1[2] + c2*k2[2] + c3*k3[2] + c4*k4[2]))
445 c1*k1[0] + c2*k2[0] + c3*k3[0] + c4*k4[0],((forwDir)[0] = (c1*k1[0] + c2*k2[0] + c3*k3[0] + c4*k4[0]), (
forwDir)[1] = (c1*k1[1] + c2*k2[1] + c3*k3[1] + c4*k4[1]), (forwDir
)[2] = (c1*k1[2] + c2*k2[2] + c3*k3[2] + c4*k4[2]))
446 c1*k1[1] + c2*k2[1] + c3*k3[1] + c4*k4[1],((forwDir)[0] = (c1*k1[0] + c2*k2[0] + c3*k3[0] + c4*k4[0]), (
forwDir)[1] = (c1*k1[1] + c2*k2[1] + c3*k3[1] + c4*k4[1]), (forwDir
)[2] = (c1*k1[2] + c2*k2[2] + c3*k3[2] + c4*k4[2]))
447 c1*k1[2] + c2*k2[2] + c3*k3[2] + c4*k4[2])((forwDir)[0] = (c1*k1[0] + c2*k2[0] + c3*k3[0] + c4*k4[0]), (
forwDir)[1] = (c1*k1[1] + c2*k2[1] + c3*k3[1] + c4*k4[1]), (forwDir
)[2] = (c1*k1[2] + c2*k2[2] + c3*k3[2] + c4*k4[2]));
448
449 return 0;
450}
451
452int (*
453_tenFiberIntegrate[TEN_FIBER_INTG_MAX3+1])(tenFiberContext *tfx, double *) = {
454 NULL((void*)0),
455 _tenFiberIntegrate_Euler,
456 _tenFiberIntegrate_Midpoint,
457 _tenFiberIntegrate_RK4
458};
459
460/*
461** modified body of previous tenFiberTraceSet, in order to
462** permit passing the nval for storing desired probed values
463*/
464static int
465_fiberTraceSet(tenFiberContext *tfx, Nrrd *nval, Nrrd *nfiber,
466 double *buff, unsigned int halfBuffLen,
467 unsigned int *startIdxP, unsigned int *endIdxP,
468 double seed[3]) {
469 static const char me[]="_fiberTraceSet";
470 airArray *fptsArr[2], /* airArrays of backward (0) and forward (1)
471 fiber points */
472 *pansArr[2]; /* airArrays of backward (0) and forward (1)
473 probed values */
474 double *fpts[2], /* arrays storing forward and backward
475 fiber points */
476 *pans[2], /* arrays storing forward and backward
477 probed values */
478 tmp[3],
479 iPos[3],
480 currPoint[3],
481 forwDir[3],
482 *fiber, /* array of both forward and backward points,
483 when finished */
484 *valOut; /* same for probed values */
485 const double *pansP; /* pointer to gage's probed values */
486
487 int gret, whyStop, buffIdx, fptsIdx, pansIdx, outIdx, oldStop, keepfiber;
488 unsigned int i, pansLen;
489 airArray *mop;
490 airPtrPtrUnion appu;
491
492 if (!(tfx)) {
493 biffAddf(TENtenBiffKey, "%s: got NULL pointer", me);
494 return 1;
495 }
496 if (nval) {
497 if (!tfx->fiberProbeItem) {
498 biffAddf(TENtenBiffKey, "%s: want to record probed values but no item set", me);
499 return 1;
500 }
501 pansLen = gageAnswerLength(tfx->gtx, tfx->pvl, tfx->fiberProbeItem);
502 pansP = gageAnswerPointer(tfx->gtx, tfx->pvl, tfx->fiberProbeItem);
503 } else {
504 pansLen = 0;
505 pansP = NULL((void*)0);
506 }
507 /*
508 fprintf(stderr, "!%s: =========================== \n", me);
509 fprintf(stderr, "!%s: \n", me);
510 fprintf(stderr, "!%s: item %d -> pansLen = %u\n", me,
511 tfx->fiberProbeItem, pansLen);
512 fprintf(stderr, "!%s: \n", me);
513 fprintf(stderr, "!%s: =========================== \n", me);
514 */
515
516 /* HEY: a hack to preserve the state inside tenFiberContext so that
517 we have fewer side effects (tfx->maxNumSteps may still be set) */
518 oldStop = tfx->stop;
519 if (!nfiber) {
520 if (!( buff && halfBuffLen > 0 && startIdxP && startIdxP )) {
521 biffAddf(TENtenBiffKey, "%s: need either non-NULL nfiber or fpts buffer info", me);
522 return 1;
523 }
524 if (tenFiberStopSet(tfx, tenFiberStopNumSteps, halfBuffLen)) {
525 biffAddf(TENtenBiffKey, "%s: error setting new fiber stop", me);
526 return 1;
527 }
528 }
529
530 /* initialize the quantities which describe the fiber halves */
531 tfx->halfLen[0] = tfx->halfLen[1] = 0.0;
532 tfx->numSteps[0] = tfx->numSteps[1] = 0;
533 tfx->whyStop[0] = tfx->whyStop[1] = tenFiberStopUnknown;
534 /*
535 fprintf(stderr, "!%s: try probing once, at seed %g %g %g\n", me,
536 seed[0], seed[1], seed[2]);
537 */
538 /* try probing once, at seed point */
539 if (tfx->useIndexSpace) {
540 gageShapeItoW(tfx->gtx->shape, tmp, seed);
541 } else {
542 ELL_3V_COPY(tmp, seed)((tmp)[0] = (seed)[0], (tmp)[1] = (seed)[1], (tmp)[2] = (seed
)[2]);
543 }
544 if (_tenFiberProbe(tfx, &gret, tmp, AIR_TRUE1)) {
545 biffAddf(TENtenBiffKey, "%s: first _tenFiberProbe failed", me);
546 return 1;
547 }
548 if (gret) {
549 if (gageErrBoundsSpace != tfx->gtx->errNum) {
550 biffAddf(TENtenBiffKey, "%s: gage problem on first _tenFiberProbe: %s (%d)",
551 me, tfx->gtx->errStr, tfx->gtx->errNum);
552 return 1;
553 } else {
554 /* the problem on the first probe was that it was out of bounds,
555 which is not a catastrophe; its handled the same as below */
556 tfx->whyNowhere = tenFiberStopBounds;
557 if (nval) {
558 nrrdEmpty(nval);
559 }
560 if (nfiber) {
561 nrrdEmpty(nfiber);
562 } else {
563 *startIdxP = *endIdxP = 0;
564 }
565 return 0;
566 }
567 }
568
569 /* see if we're doomed (tract dies before it gets anywhere) */
570 /* have to fake out the possible radius check, since at this point
571 there is no radius of curvature; this will always pass */
572 tfx->radius = DBL_MAX1.7976931348623157e+308;
573 if ((whyStop = _tenFiberStopCheck(tfx))) {
574 /* stopped immediately at seed point, but that's not an error */
575 tfx->whyNowhere = whyStop;
576 if (nval) {
577 nrrdEmpty(nval);
578 }
579 if (nfiber) {
580 nrrdEmpty(nfiber);
581 } else {
582 *startIdxP = *endIdxP = 0;
583 }
584 return 0;
585 } else {
586 /* did not immediately halt */
587 tfx->whyNowhere = tenFiberStopUnknown;
588 }
589
590 /* airMop{Error,Okay}() can safely be called on NULL */
591 mop = (nfiber || nval) ? airMopNew() : NULL((void*)0);
592
593 for (tfx->halfIdx=0; tfx->halfIdx<=1; tfx->halfIdx++) {
594 if (nval) {
595 appu.d = &(pans[tfx->halfIdx]);
596 pansArr[tfx->halfIdx] = airArrayNew(appu.v, NULL((void*)0),
597 pansLen*sizeof(double),
598 TEN_FIBER_INCR512);
599 airMopAdd(mop, pansArr[tfx->halfIdx],
600 (airMopper)airArrayNuke, airMopAlways);
601 } else {
602 pansArr[tfx->halfIdx] = NULL((void*)0);
603 }
604 pansIdx = -1;
605 if (nfiber) {
606 appu.d = &(fpts[tfx->halfIdx]);
607 fptsArr[tfx->halfIdx] = airArrayNew(appu.v, NULL((void*)0),
608 3*sizeof(double), TEN_FIBER_INCR512);
609 airMopAdd(mop, fptsArr[tfx->halfIdx],
610 (airMopper)airArrayNuke, airMopAlways);
611 buffIdx = -1;
612 } else {
613 fptsArr[tfx->halfIdx] = NULL((void*)0);
614 fpts[tfx->halfIdx] = NULL((void*)0);
615 buffIdx = halfBuffLen;
616 }
617 fptsIdx = -1;
Value stored to 'fptsIdx' is never read
618 tfx->halfLen[tfx->halfIdx] = 0;
619 if (tfx->useIndexSpace) {
620 ELL_3V_COPY(iPos, seed)((iPos)[0] = (seed)[0], (iPos)[1] = (seed)[1], (iPos)[2] = (seed
)[2]);
621 gageShapeItoW(tfx->gtx->shape, tfx->wPos, iPos);
622 } else {
623 /*
624 fprintf(stderr, "!%s(A): %p %p %p\n", me,
625 tfx->gtx->shape, iPos, seed);
626 */
627 gageShapeWtoI(tfx->gtx->shape, iPos, seed);
628 ELL_3V_COPY(tfx->wPos, seed)((tfx->wPos)[0] = (seed)[0], (tfx->wPos)[1] = (seed)[1]
, (tfx->wPos)[2] = (seed)[2]);
629 }
630 /* have to initially pass the possible radius check in
631 _tenFiberStopCheck(); this will always pass */
632 tfx->radius = DBL_MAX1.7976931348623157e+308;
633 ELL_3V_SET(tfx->lastDir, 0, 0, 0)((tfx->lastDir)[0] = (0), (tfx->lastDir)[1] = (0), (tfx
->lastDir)[2] = (0));
634 tfx->lastDirSet = AIR_FALSE0;
635 for (tfx->numSteps[tfx->halfIdx] = 0;
636 AIR_TRUE1;
637 tfx->numSteps[tfx->halfIdx]++) {
638 _tenFiberProbe(tfx, &gret, tfx->wPos, AIR_FALSE0);
639 if (gret) {
640 /* even if gageProbe had an error OTHER than going out of bounds,
641 we're not going to report it any differently here, alas */
642 tfx->whyStop[tfx->halfIdx] = tenFiberStopBounds;
643 /*
644 fprintf(stderr, "!%s: A tfx->whyStop[%d] = %s\n", me, tfx->halfIdx,
645 airEnumStr(tenFiberStop, tfx->whyStop[tfx->halfIdx]));
646 */
647 break;
648 }
649 if ((whyStop = _tenFiberStopCheck(tfx))) {
650 if (tenFiberStopNumSteps == whyStop) {
651 /* we stopped along this direction because
652 tfx->numSteps[tfx->halfIdx] exceeded tfx->maxNumSteps.
653 Okay. But tfx->numSteps[tfx->halfIdx] is supposed to be
654 a record of how steps were (successfully) taken. So we
655 need to decrementing before moving on ... */
656 tfx->numSteps[tfx->halfIdx]--;
657 }
658 tfx->whyStop[tfx->halfIdx] = whyStop;
659 /*
660 fprintf(stderr, "!%s: B tfx->whyStop[%d] = %s\n", me, tfx->halfIdx,
661 airEnumStr(tenFiberStop, tfx->whyStop[tfx->halfIdx]));
662 */
663 break;
664 }
665 if (tfx->useIndexSpace) {
666 /*
667 fprintf(stderr, "!%s(B): %p %p %p\n", me,
668 tfx->gtx->shape, iPos, tfx->wPos);
669 */
670 gageShapeWtoI(tfx->gtx->shape, iPos, tfx->wPos);
671 ELL_3V_COPY(currPoint, iPos)((currPoint)[0] = (iPos)[0], (currPoint)[1] = (iPos)[1], (currPoint
)[2] = (iPos)[2]);
672 } else {
673 ELL_3V_COPY(currPoint, tfx->wPos)((currPoint)[0] = (tfx->wPos)[0], (currPoint)[1] = (tfx->
wPos)[1], (currPoint)[2] = (tfx->wPos)[2]);
674 }
675 if (nval) {
676 pansIdx = airArrayLenIncr(pansArr[tfx->halfIdx], 1);
677 /* HEY: speed this up */
678 memcpy(pans[tfx->halfIdx] + pansLen*pansIdx, pansP,__builtin___memcpy_chk (pans[tfx->halfIdx] + pansLen*pansIdx
, pansP, pansLen*sizeof(double), __builtin_object_size (pans[
tfx->halfIdx] + pansLen*pansIdx, 0))
679 pansLen*sizeof(double))__builtin___memcpy_chk (pans[tfx->halfIdx] + pansLen*pansIdx
, pansP, pansLen*sizeof(double), __builtin_object_size (pans[
tfx->halfIdx] + pansLen*pansIdx, 0));
680 /*
681 fprintf(stderr, "!%s: (dir %d) saving to %d: %g @ (%g,%g,%g)\n", me,
682 tfx->halfIdx, pansIdx, pansP[0],
683 currPoint[0], currPoint[1], currPoint[2]);
684 */
685 }
686 if (nfiber) {
687 fptsIdx = airArrayLenIncr(fptsArr[tfx->halfIdx], 1);
688 ELL_3V_COPY(fpts[tfx->halfIdx] + 3*fptsIdx, currPoint)((fpts[tfx->halfIdx] + 3*fptsIdx)[0] = (currPoint)[0], (fpts
[tfx->halfIdx] + 3*fptsIdx)[1] = (currPoint)[1], (fpts[tfx
->halfIdx] + 3*fptsIdx)[2] = (currPoint)[2]);
689 } else {
690 ELL_3V_COPY(buff + 3*buffIdx, currPoint)((buff + 3*buffIdx)[0] = (currPoint)[0], (buff + 3*buffIdx)[1
] = (currPoint)[1], (buff + 3*buffIdx)[2] = (currPoint)[2]);
691 /*
692 fprintf(stderr, "!%s: (dir %d) saving to %d pnt %g %g %g\n", me,
693 tfx->halfIdx, buffIdx,
694 currPoint[0], currPoint[1], currPoint[2]);
695 */
696 buffIdx += !tfx->halfIdx ? -1 : 1;
697 }
698 /* forwDir is set by this to point to the next fiber point */
699 if (_tenFiberIntegrate[tfx->intg](tfx, forwDir)) {
700 tfx->whyStop[tfx->halfIdx] = tenFiberStopBounds;
701 /*
702 fprintf(stderr, "!%s: C tfx->whyStop[%d] = %s\n", me, tfx->halfIdx,
703 airEnumStr(tenFiberStop, tfx->whyStop[tfx->halfIdx]));
704 */
705 break;
706 }
707 /*
708 fprintf(stderr, "!%s: forwDir = %g %g %g\n", me,
709 forwDir[0], forwDir[1], forwDir[2]);
710 */
711 if (tfx->stop & (1 << tenFiberStopRadius)) {
712 /* some more work required to compute radius of curvature */
713 double svec[3], dvec[3], SS, DD, dlen; /* sum,diff length squared */
714 /* tfx->lastDir and forwDir are not normalized to unit-length */
715 if (tfx->lastDirSet) {
716 ELL_3V_ADD2(svec, tfx->lastDir, forwDir)((svec)[0] = (tfx->lastDir)[0] + (forwDir)[0], (svec)[1] =
(tfx->lastDir)[1] + (forwDir)[1], (svec)[2] = (tfx->lastDir
)[2] + (forwDir)[2]);
717 ELL_3V_SUB(dvec, tfx->lastDir, forwDir)((dvec)[0] = (tfx->lastDir)[0] - (forwDir)[0], (dvec)[1] =
(tfx->lastDir)[1] - (forwDir)[1], (dvec)[2] = (tfx->lastDir
)[2] - (forwDir)[2]);
718 SS = ELL_3V_DOT(svec, svec)((svec)[0]*(svec)[0] + (svec)[1]*(svec)[1] + (svec)[2]*(svec)
[2]);
719 DD = ELL_3V_DOT(dvec, dvec)((dvec)[0]*(dvec)[0] + (dvec)[1]*(dvec)[1] + (dvec)[2]*(dvec)
[2]);
720 /* Sun Nov 2 00:04:05 EDT 2008: GLK can't recover how he
721 derived this, and can't see why it would be corrrect,
722 even though it seems to work correctly...
723 tfx->radius = sqrt(SS*(SS+DD)/DD)/4;
724 */
725 dlen = sqrt(DD);
726 tfx->radius = dlen ? (SS + DD)/(4*dlen) : DBL_MAX1.7976931348623157e+308;
727 } else {
728 tfx->radius = DBL_MAX1.7976931348623157e+308;
729 }
730 }
731 /*
732 if (!tfx->lastDirSet) {
733 fprintf(stderr, "!%s: now setting lastDirSet to (%g,%g,%g)\n", me,
734 forwDir[0], forwDir[1], forwDir[2]);
735 }
736 */
737 ELL_3V_COPY(tfx->lastDir, forwDir)((tfx->lastDir)[0] = (forwDir)[0], (tfx->lastDir)[1] = (
forwDir)[1], (tfx->lastDir)[2] = (forwDir)[2]);
738 tfx->lastDirSet = AIR_TRUE1;
739 ELL_3V_ADD2(tfx->wPos, tfx->wPos, forwDir)((tfx->wPos)[0] = (tfx->wPos)[0] + (forwDir)[0], (tfx->
wPos)[1] = (tfx->wPos)[1] + (forwDir)[1], (tfx->wPos)[2
] = (tfx->wPos)[2] + (forwDir)[2]);
740 tfx->halfLen[tfx->halfIdx] += ELL_3V_LEN(forwDir)(sqrt((((forwDir))[0]*((forwDir))[0] + ((forwDir))[1]*((forwDir
))[1] + ((forwDir))[2]*((forwDir))[2])));
741 }
742 }
743
744 keepfiber = AIR_TRUE1;
745 if ((tfx->stop & (1 << tenFiberStopStub))
746 && (2 == fptsArr[0]->len + fptsArr[1]->len)) {
747 /* seed point was actually valid, but neither half got anywhere,
748 and the user has set tenFiberStopStub, so we report this as
749 a non-starter, via tfx->whyNowhere. */
750 tfx->whyNowhere = tenFiberStopStub;
751 keepfiber = AIR_FALSE0;
752 }
753 if ((tfx->stop & (1 << tenFiberStopMinNumSteps))
754 && (fptsArr[0]->len + fptsArr[1]->len < tfx->minNumSteps)) {
755 /* whole fiber didn't have enough steps */
756 tfx->whyNowhere = tenFiberStopMinNumSteps;
757 keepfiber = AIR_FALSE0;
758 }
759 if ((tfx->stop & (1 << tenFiberStopMinLength))
760 && (tfx->halfLen[0] + tfx->halfLen[1] < tfx->minWholeLen)) {
761 /* whole fiber wasn't long enough */
762 tfx->whyNowhere = tenFiberStopMinLength;
763 keepfiber = AIR_FALSE0;
764 }
765 if (!keepfiber) {
766 /* for the curious, tfx->whyStop[0,1], tfx->numSteps[0,1], and
767 tfx->halfLen[1,2] remain set, from above */
768 if (nval) {
769 nrrdEmpty(nval);
770 }
771 if (nfiber) {
772 nrrdEmpty(nfiber);
773 } else {
774 *startIdxP = *endIdxP = 0;
775 }
776 } else {
777 if (nval) {
778 if (nrrdMaybeAlloc_va(nval, nrrdTypeDouble, 2,
779 AIR_CAST(size_t, pansLen)((size_t)(pansLen)),
780 AIR_CAST(size_t, (pansArr[0]->len((size_t)((pansArr[0]->len + pansArr[1]->len - 1)))
781 + pansArr[1]->len - 1))((size_t)((pansArr[0]->len + pansArr[1]->len - 1))))) {
782 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate probed value nrrd", me);
783 airMopError(mop); return 1;
784 }
785 valOut = AIR_CAST(double*, nval->data)((double*)(nval->data));
786 outIdx = 0;
787 /* HEY: speed up memcpy */
788 for (i=pansArr[0]->len-1; i>=1; i--) {
789 memcpy(valOut + pansLen*outIdx, pans[0] + pansLen*i,__builtin___memcpy_chk (valOut + pansLen*outIdx, pans[0] + pansLen
*i, pansLen*sizeof(double), __builtin_object_size (valOut + pansLen
*outIdx, 0))
790 pansLen*sizeof(double))__builtin___memcpy_chk (valOut + pansLen*outIdx, pans[0] + pansLen
*i, pansLen*sizeof(double), __builtin_object_size (valOut + pansLen
*outIdx, 0));
791 outIdx++;
792 }
793 for (i=0; i<=pansArr[1]->len-1; i++) {
794 memcpy(valOut + pansLen*outIdx, pans[1] + pansLen*i,__builtin___memcpy_chk (valOut + pansLen*outIdx, pans[1] + pansLen
*i, pansLen*sizeof(double), __builtin_object_size (valOut + pansLen
*outIdx, 0))
795 pansLen*sizeof(double))__builtin___memcpy_chk (valOut + pansLen*outIdx, pans[1] + pansLen
*i, pansLen*sizeof(double), __builtin_object_size (valOut + pansLen
*outIdx, 0));
796 outIdx++;
797 }
798 }
799 if (nfiber) {
800 if (nrrdMaybeAlloc_va(nfiber, nrrdTypeDouble, 2,
801 AIR_CAST(size_t, 3)((size_t)(3)),
802 AIR_CAST(size_t, (fptsArr[0]->len((size_t)((fptsArr[0]->len + fptsArr[1]->len - 1)))
803 + fptsArr[1]->len - 1))((size_t)((fptsArr[0]->len + fptsArr[1]->len - 1))))) {
804 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate fiber nrrd", me);
805 airMopError(mop); return 1;
806 }
807 fiber = AIR_CAST(double*, nfiber->data)((double*)(nfiber->data));
808 outIdx = 0;
809 for (i=fptsArr[0]->len-1; i>=1; i--) {
810 ELL_3V_COPY(fiber + 3*outIdx, fpts[0] + 3*i)((fiber + 3*outIdx)[0] = (fpts[0] + 3*i)[0], (fiber + 3*outIdx
)[1] = (fpts[0] + 3*i)[1], (fiber + 3*outIdx)[2] = (fpts[0] +
3*i)[2]);
811 outIdx++;
812 }
813 for (i=0; i<=fptsArr[1]->len-1; i++) {
814 ELL_3V_COPY(fiber + 3*outIdx, fpts[1] + 3*i)((fiber + 3*outIdx)[0] = (fpts[1] + 3*i)[0], (fiber + 3*outIdx
)[1] = (fpts[1] + 3*i)[1], (fiber + 3*outIdx)[2] = (fpts[1] +
3*i)[2]);
815 outIdx++;
816 }
817 } else {
818 *startIdxP = halfBuffLen - tfx->numSteps[0];
819 *endIdxP = halfBuffLen + tfx->numSteps[1];
820 }
821 }
822
823 tfx->stop = oldStop;
824 airMopOkay(mop);
825 return 0;
826}
827
828/*
829******** tenFiberTraceSet
830**
831** slightly more flexible API for fiber tracking than tenFiberTrace
832**
833** EITHER: pass a non-NULL nfiber, and NULL, 0, NULL, NULL for
834** the following arguments, and things are the same as with tenFiberTrace:
835** data inside the nfiber is allocated, and the tract vertices are copied
836** into it, having been stored in dynamically allocated airArrays
837**
838** OR: pass a NULL nfiber, and a buff allocated for 3*(2*halfBuffLen + 1)
839** (note the "+ 1" !!!) doubles. The fiber tracking on each half will stop
840** at halfBuffLen points. The given seedpoint will be stored in
841** buff[0,1,2 + 3*halfBuffLen]. The linear (1-D) indices for the end of
842** the first tract half, and the end of the second tract half, will be set in
843** *startIdxP and *endIdxP respectively (this does not include a multiply
844** by 3)
845**
846** it is worth pointing out here that internally, all tractography is done
847** in gage's world space, regardless of tfx->useIndexSpace. The conversion
848** from/to index is space (if tfx->useIndexSpace is non-zero) is only done
849** for seedpoints and when fiber vertices are saved out, respectively.
850**
851** As of Sun Aug 1 20:40:55 CDT 2010 this is just a wrapper around
852** _fiberTraceSet; this will probably change in Teem 2.0
853*/
854int
855tenFiberTraceSet(tenFiberContext *tfx, Nrrd *nfiber,
856 double *buff, unsigned int halfBuffLen,
857 unsigned int *startIdxP, unsigned int *endIdxP,
858 double seed[3]) {
859 static const char me[]="tenFiberTraceSet";
860
861 if (_fiberTraceSet(tfx, NULL((void*)0), nfiber, buff, halfBuffLen,
862 startIdxP, endIdxP, seed)) {
863 biffAddf(TENtenBiffKey, "%s: problem", me);
864 return 1;
865 }
866
867 return 0;
868}
869
870/*
871******** tenFiberTrace
872**
873** takes a starting position in index or world space, depending on the
874** value of tfx->useIndexSpace
875*/
876int
877tenFiberTrace(tenFiberContext *tfx, Nrrd *nfiber, double seed[3]) {
878 static const char me[]="tenFiberTrace";
879
880 if (_fiberTraceSet(tfx, NULL((void*)0), nfiber, NULL((void*)0), 0, NULL((void*)0), NULL((void*)0), seed)) {
881 biffAddf(TENtenBiffKey, "%s: problem computing tract", me);
882 return 1;
883 }
884
885 return 0;
886}
887
888/*
889******** tenFiberDirectionNumber
890**
891** NOTE: for the time being, a return of zero indicates an error, not
892** that we're being clever and detect that the seedpoint is in such
893** isotropy that no directions are possible (though such cleverness
894** will hopefully be implemented soon)
895*/
896unsigned int
897tenFiberDirectionNumber(tenFiberContext *tfx, double seed[3]) {
898 static const char me[]="tenFiberDirectionNumber";
899 unsigned int ret;
900
901 if (!(tfx && seed)) {
902 biffAddf(TENtenBiffKey, "%s: got NULL pointer", me);
903 return 0;
904 }
905
906 /* HEY: eventually this stuff will be specific to the seedpoint ... */
907
908 if (tfx->useDwi) {
909 switch (tfx->fiberType) {
910 case tenDwiFiberType1Evec0:
911 ret = 1;
912 break;
913 case tenDwiFiberType2Evec0:
914 ret = 2;
915 break;
916 case tenDwiFiberType12BlendEvec0:
917 biffAddf(TENtenBiffKey, "%s: sorry, type %s not yet implemented", me,
918 airEnumStr(tenDwiFiberType, tenDwiFiberType12BlendEvec0));
919 ret = 0;
920 break;
921 default:
922 biffAddf(TENtenBiffKey, "%s: type %d unknown!", me, tfx->fiberType);
923 ret = 0;
924 break;
925 }
926 } else {
927 /* not using DWIs */
928 ret = 1;
929 }
930
931 return ret;
932}
933
934/*
935******** tenFiberSingleTrace
936**
937** fiber tracing API that uses new tenFiberSingle, as well as being
938** aware of multi-direction tractography
939**
940** NOTE: this will not try any cleverness in setting "num"
941** according to whether the seedpoint is a non-starter
942*/
943int
944tenFiberSingleTrace(tenFiberContext *tfx, tenFiberSingle *tfbs,
945 double seed[3], unsigned int which) {
946 static const char me[]="tenFiberSingleTrace";
947
948 if (!(tfx && tfbs && seed)) {
949 biffAddf(TENtenBiffKey, "%s: got NULL pointer", me);
950 return 1;
951 }
952
953 /* set input fields in tfbs */
954 ELL_3V_COPY(tfbs->seedPos, seed)((tfbs->seedPos)[0] = (seed)[0], (tfbs->seedPos)[1] = (
seed)[1], (tfbs->seedPos)[2] = (seed)[2]);
955 tfbs->dirIdx = which;
956 /* not our job to set tfbx->dirNum ... */
957
958 /* set tfbs->nvert */
959 /* no harm in setting this even when there are no multiple fibers */
960 tfx->ten2Which = which;
961 if (_fiberTraceSet(tfx, (tfx->fiberProbeItem ? tfbs->nval : NULL((void*)0)),
962 tfbs->nvert, NULL((void*)0), 0, NULL((void*)0), NULL((void*)0), seed)) {
963 biffAddf(TENtenBiffKey, "%s: problem computing tract", me);
964 return 1;
965 }
966
967 /* set other fields based on tfx output */
968 tfbs->halfLen[0] = tfx->halfLen[0];
969 tfbs->halfLen[1] = tfx->halfLen[1];
970 tfbs->seedIdx = tfx->numSteps[0];
971 tfbs->stepNum[0] = tfx->numSteps[0];
972 tfbs->stepNum[1] = tfx->numSteps[1];
973 tfbs->whyStop[0] = tfx->whyStop[0];
974 tfbs->whyStop[1] = tfx->whyStop[1];
975 tfbs->whyNowhere = tfx->whyNowhere;
976
977 return 0;
978}
979
980typedef union {
981 tenFiberSingle **f;
982 void **v;
983} fiberunion;
984
985/* uses biff */
986tenFiberMulti *
987tenFiberMultiNew() {
988 static const char me[]="tenFiberMultiNew";
989 tenFiberMulti *ret;
990 fiberunion tfu;
991
992 ret = AIR_CAST(tenFiberMulti *, calloc(1, sizeof(tenFiberMulti)))((tenFiberMulti *)(calloc(1, sizeof(tenFiberMulti))));
993 if (ret) {
994 ret->fiber = NULL((void*)0);
995 ret->fiberNum = 0;
996 tfu.f = &(ret->fiber);
997 ret->fiberArr = airArrayNew(tfu.v, &(ret->fiberNum),
998 sizeof(tenFiberSingle), 512 /* incr */);
999 if (ret->fiberArr) {
1000 airArrayStructCB(ret->fiberArr,
1001 AIR_CAST(void (*)(void *), tenFiberSingleInit)((void (*)(void *))(tenFiberSingleInit)),
1002 AIR_CAST(void (*)(void *), tenFiberSingleDone)((void (*)(void *))(tenFiberSingleDone)));
1003 } else {
1004 biffAddf(TENtenBiffKey, "%s: couldn't create airArray", me);
1005 return NULL((void*)0);
1006 }
1007 } else {
1008 biffAddf(TENtenBiffKey, "%s: couldn't create tenFiberMulti", me);
1009 return NULL((void*)0);
1010 }
1011 return ret;
1012}
1013
1014int
1015tenFiberMultiCheck(airArray *arr) {
1016 static const char me[]="tenFiberMultiCheck";
1017
1018 if (!arr) {
1019 biffAddf(TENtenBiffKey, "%s: got NULL pointer", me);
1020 return 1;
1021 }
1022 if (sizeof(tenFiberSingle) != arr->unit) {
1023 biffAddf(TENtenBiffKey, "%s: given airArray cannot be for fibers", me);
1024 return 1;
1025 }
1026 if (!(AIR_CAST(void (*)(void *), tenFiberSingleInit)((void (*)(void *))(tenFiberSingleInit)) == arr->initCB
1027 && AIR_CAST(void (*)(void *), tenFiberSingleDone)((void (*)(void *))(tenFiberSingleDone)) == arr->doneCB)) {
1028 biffAddf(TENtenBiffKey, "%s: given airArray not set up with fiber callbacks", me);
1029 return 1;
1030 }
1031 return 0;
1032}
1033
1034tenFiberMulti *
1035tenFiberMultiNix(tenFiberMulti *tfm) {
1036
1037 if (tfm) {
1038 airArrayNuke(tfm->fiberArr);
1039 airFree(tfm);
1040 }
1041 return NULL((void*)0);
1042}
1043
1044/*
1045******** tenFiberMultiTrace
1046**
1047** does tractography for a list of seedpoints
1048**
1049** tfml has been returned from tenFiberMultiNew()
1050*/
1051int
1052tenFiberMultiTrace(tenFiberContext *tfx, tenFiberMulti *tfml,
1053 const Nrrd *_nseed) {
1054 static const char me[]="tenFiberMultiTrace";
1055 airArray *mop;
1056 const double *seedData;
1057 double seed[3];
1058 unsigned int seedNum, seedIdx, fibrNum, dirNum, dirIdx;
1059 Nrrd *nseed;
1060
1061 if (!(tfx && tfml && _nseed)) {
1062 biffAddf(TENtenBiffKey, "%s: got NULL pointer", me);
1063 return 1;
1064 }
1065 if (tenFiberMultiCheck(tfml->fiberArr)) {
1066 biffAddf(TENtenBiffKey, "%s: problem with fiber array", me);
1067 return 1;
1068 }
1069 if (!(2 == _nseed->dim && 3 == _nseed->axis[0].size)) {
1070 biffAddf(TENtenBiffKey, "%s: seed list should be a 2-D (not %u-D) "
1071 "3-by-X (not %u-by-X) array", me, _nseed->dim,
1072 AIR_CAST(unsigned int, _nseed->axis[0].size)((unsigned int)(_nseed->axis[0].size)));
1073 return 1;
1074 }
1075
1076 mop = airMopNew();
1077
1078 seedNum = _nseed->axis[1].size;
1079 if (nrrdTypeDouble == _nseed->type) {
1080 seedData = AIR_CAST(const double *, _nseed->data)((const double *)(_nseed->data));
1081 } else {
1082 nseed = nrrdNew();
1083 airMopAdd(mop, nseed, AIR_CAST(airMopper, nrrdNuke)((airMopper)(nrrdNuke)), airMopAlways);
1084 if (nrrdConvert(nseed, _nseed, nrrdTypeDouble)) {
1085 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: couldn't convert seed list", me);
1086 return 1;
1087 }
1088 seedData = AIR_CAST(const double *, nseed->data)((const double *)(nseed->data));
1089 }
1090
1091 /* HEY: the correctness of the use of the airArray here is quite subtle */
1092 fibrNum = 0;
1093 for (seedIdx=0; seedIdx<seedNum; seedIdx++) {
1094 dirNum = tenFiberDirectionNumber(tfx, seed);
1095 if (!dirNum) {
1096 biffAddf(TENtenBiffKey, "%s: couldn't learn dirNum at seed (%g,%g,%g)", me,
1097 seed[0], seed[1], seed[2]);
1098 return 1;
1099 }
1100 for (dirIdx=0; dirIdx<dirNum; dirIdx++) {
1101 if (tfx->verbose > 1) {
1102 fprintf(stderr__stderrp, "%s: dir %u/%u on seed %u/%u; len %u; # %u\n",
1103 me, dirIdx, dirNum, seedIdx, seedNum,
1104 tfml->fiberArr->len, fibrNum);
1105 }
1106 /* tfml->fiberArr->len can never be < fibrNum */
1107 if (tfml->fiberArr->len == fibrNum) {
1108 airArrayLenIncr(tfml->fiberArr, 1);
1109 }
1110 ELL_3V_COPY(tfml->fiber[fibrNum].seedPos, seedData + 3*seedIdx)((tfml->fiber[fibrNum].seedPos)[0] = (seedData + 3*seedIdx
)[0], (tfml->fiber[fibrNum].seedPos)[1] = (seedData + 3*seedIdx
)[1], (tfml->fiber[fibrNum].seedPos)[2] = (seedData + 3*seedIdx
)[2]);
1111 tfml->fiber[fibrNum].dirIdx = dirIdx;
1112 tfml->fiber[fibrNum].dirNum = dirNum;
1113 ELL_3V_COPY(seed, seedData + 3*seedIdx)((seed)[0] = (seedData + 3*seedIdx)[0], (seed)[1] = (seedData
+ 3*seedIdx)[1], (seed)[2] = (seedData + 3*seedIdx)[2]);
1114 if (tenFiberSingleTrace(tfx, &(tfml->fiber[fibrNum]), seed, dirIdx)) {
1115 biffAddf(TENtenBiffKey, "%s: trouble on seed (%g,%g,%g) %u/%u, dir %u/%u", me,
1116 seed[0], seed[1], seed[2], seedIdx, seedNum, dirIdx, dirNum);
1117 return 1;
1118 }
1119 if (tfx->verbose) {
1120 if (tenFiberStopUnknown == tfml->fiber[fibrNum].whyNowhere) {
1121 fprintf(stderr__stderrp, "%s: (%g,%g,%g) ->\n"
1122 " steps = %u,%u; len = %g,%g; whyStop = %s,%s\n",
1123 me, seed[0], seed[1], seed[2],
1124 tfml->fiber[fibrNum].stepNum[0],
1125 tfml->fiber[fibrNum].stepNum[1],
1126 tfml->fiber[fibrNum].halfLen[0],
1127 tfml->fiber[fibrNum].halfLen[1],
1128 airEnumStr(tenFiberStop, tfml->fiber[fibrNum].whyStop[0]),
1129 airEnumStr(tenFiberStop, tfml->fiber[fibrNum].whyStop[1]));
1130 } else {
1131 fprintf(stderr__stderrp, "%s: (%g,%g,%g) -> whyNowhere: %s\n",
1132 me, seed[0], seed[1], seed[2],
1133 airEnumStr(tenFiberStop, tfml->fiber[fibrNum].whyNowhere));
1134 }
1135 }
1136 fibrNum++;
1137 }
1138 }
1139 /* if the airArray got to be its length only because of the work above,
1140 then the following will be a no-op. Otherwise, via the callbacks,
1141 it will clear out the tenFiberSingle's that we didn't create here */
1142 airArrayLenSet(tfml->fiberArr, fibrNum);
1143
1144 airMopOkay(mop);
1145 return 0;
1146}
1147
1148static int
1149_fiberMultiExtract(tenFiberContext *tfx, Nrrd *nval,
1150 limnPolyData *lpld, tenFiberMulti *tfml) {
1151 static const char me[]="_fiberMultiExtract";
1152 unsigned int seedIdx, vertTotalNum, fiberNum, fiberIdx, vertTotalIdx,
1153 pansLen, pvNum;
1154 double *valOut;
1155
1156 if (!(tfx && (lpld || nval) && tfml)) {
1157 biffAddf(TENtenBiffKey, "%s: got NULL pointer", me);
1158 return 1;
1159 }
1160 if (tenFiberMultiCheck(tfml->fiberArr)) {
1161 biffAddf(TENtenBiffKey, "%s: problem with fiber array", me);
1162 return 1;
1163 }
1164 if (nval) {
1165 if (!tfx->fiberProbeItem) {
1166 biffAddf(TENtenBiffKey, "%s: want probed values but no item set", me);
1167 return 1;
1168 }
1169 pansLen = gageAnswerLength(tfx->gtx, tfx->pvl, tfx->fiberProbeItem);
1170 } else {
1171 pansLen = 0;
1172 }
1173 /*
1174 fprintf(stderr, "!%s: =========================== \n", me);
1175 fprintf(stderr, "!%s: \n", me);
1176 fprintf(stderr, "!%s: item %d -> pansLen = %u\n", me,
1177 tfx->fiberProbeItem, pansLen);
1178 fprintf(stderr, "!%s: \n", me);
1179 fprintf(stderr, "!%s: =========================== \n", me);
1180 */
1181
1182 /* we have to count the real fibers that went somewhere, excluding
1183 fibers that went nowhere (counted in tfml->fiberNum) */
1184 vertTotalNum = 0;
1185 fiberNum = 0;
1186 pvNum = 0;
1187 for (seedIdx=0; seedIdx<tfml->fiberArr->len; seedIdx++) {
1188 tenFiberSingle *tfs;
1189 tfs = tfml->fiber + seedIdx;
1190 if (!(tenFiberStopUnknown == tfs->whyNowhere)) {
1191 continue;
1192 }
1193 if (nval) {
1194 if (tfs->nval) {
1195 if (!(2 == tfs->nval->dim
1196 && pansLen == tfs->nval->axis[0].size
1197 && tfs->nvert->axis[1].size == tfs->nval->axis[1].size)) {
1198 biffAddf(TENtenBiffKey, "%s: fiber[%u]->nval seems wrong", me, seedIdx);
1199 return 1;
1200 }
1201 pvNum++;
1202 }
1203 }
1204 vertTotalNum += tfs->nvert->axis[1].size;
1205 fiberNum++;
1206 }
1207 if (nval && pvNum != fiberNum) {
1208 biffAddf(TENtenBiffKey, "%s: pvNum %u != fiberNum %u", me, pvNum, fiberNum);
1209 return 1;
1210 }
1211
1212 if (nval) {
1213 if (nrrdMaybeAlloc_va(nval, nrrdTypeDouble, 2,
1214 AIR_CAST(size_t, pansLen)((size_t)(pansLen)),
1215 AIR_CAST(size_t, vertTotalNum)((size_t)(vertTotalNum)))) {
1216 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate output", me);
1217 return 1;
1218 }
1219 valOut = AIR_CAST(double *, nval->data)((double *)(nval->data));
1220 } else {
1221 valOut = NULL((void*)0);
1222 }
1223 if (lpld) {
1224 if (limnPolyDataAlloc(lpld, 0, /* no extra per-vertex info */
1225 vertTotalNum, vertTotalNum, fiberNum)) {
1226 biffMovef(TENtenBiffKey, LIMNlimnBiffKey, "%s: couldn't allocate output", me);
1227 return 1;
1228 }
1229 }
1230
1231 fiberIdx = 0;
1232 vertTotalIdx = 0;
1233 for (seedIdx=0; seedIdx<tfml->fiberArr->len; seedIdx++) {
1234 double *vert, *pans;
1235 unsigned int vertIdx, vertNum;
1236 tenFiberSingle *tfs;
1237 tfs = tfml->fiber + seedIdx;
1238 if (!(tenFiberStopUnknown == tfs->whyNowhere)) {
1239 continue;
1240 }
1241 vertNum = tfs->nvert->axis[1].size;
1242 pans = (nval
1243 ? AIR_CAST(double*, tfs->nval->data)((double*)(tfs->nval->data))
1244 : NULL((void*)0));
1245 vert = (lpld
1246 ? AIR_CAST(double*, tfs->nvert->data)((double*)(tfs->nvert->data))
1247 : NULL((void*)0));
1248 for (vertIdx=0; vertIdx<vertNum; vertIdx++) {
1249 if (lpld) {
1250 ELL_3V_COPY_TT(lpld->xyzw + 4*vertTotalIdx, float, vert + 3*vertIdx)((lpld->xyzw + 4*vertTotalIdx)[0] = ((float)((vert + 3*vertIdx
)[0])), (lpld->xyzw + 4*vertTotalIdx)[1] = ((float)((vert +
3*vertIdx)[1])), (lpld->xyzw + 4*vertTotalIdx)[2] = ((float
)((vert + 3*vertIdx)[2])));
1251 (lpld->xyzw + 4*vertTotalIdx)[3] = 1.0;
1252 lpld->indx[vertTotalIdx] = vertTotalIdx;
1253 }
1254 if (nval) {
1255 /* HEY speed up memcpy */
1256 memcpy(valOut + pansLen*vertTotalIdx,__builtin___memcpy_chk (valOut + pansLen*vertTotalIdx, pans +
pansLen*vertIdx, pansLen*sizeof(double), __builtin_object_size
(valOut + pansLen*vertTotalIdx, 0))
1257 pans + pansLen*vertIdx,__builtin___memcpy_chk (valOut + pansLen*vertTotalIdx, pans +
pansLen*vertIdx, pansLen*sizeof(double), __builtin_object_size
(valOut + pansLen*vertTotalIdx, 0))
1258 pansLen*sizeof(double))__builtin___memcpy_chk (valOut + pansLen*vertTotalIdx, pans +
pansLen*vertIdx, pansLen*sizeof(double), __builtin_object_size
(valOut + pansLen*vertTotalIdx, 0));
1259 }
1260 vertTotalIdx++;
1261 }
1262 if (lpld) {
1263 lpld->type[fiberIdx] = limnPrimitiveLineStrip;
1264 lpld->icnt[fiberIdx] = vertNum;
1265 }
1266 fiberIdx++;
1267 }
1268
1269 return 0;
1270}
1271
1272/*
1273******** tenFiberMultiPolyData
1274**
1275** converts tenFiberMulti to polydata.
1276**
1277** currently the tenFiberContext *tfx arg is not used, but it will
1278** probably be needed in the future as the way that parameters to the
1279** polydata creation process are passed.
1280*/
1281int
1282tenFiberMultiPolyData(tenFiberContext *tfx,
1283 limnPolyData *lpld, tenFiberMulti *tfml) {
1284 static const char me[]="tenFiberMultiPolyData";
1285
1286 if (_fiberMultiExtract(tfx, NULL((void*)0), lpld, tfml)) {
1287 biffAddf(TENtenBiffKey, "%s: problem", me);
1288 return 1;
1289 }
1290 return 0;
1291}
1292
1293
1294int
1295tenFiberMultiProbeVals(tenFiberContext *tfx,
1296 Nrrd *nval, tenFiberMulti *tfml) {
1297 static const char me[]="tenFiberMultiProbeVals";
1298
1299 if (_fiberMultiExtract(tfx, nval, NULL((void*)0), tfml)) {
1300 biffAddf(TENtenBiffKey, "%s: problem", me);
1301 return 1;
1302 }
1303 return 0;
1304}