Bug Summary

File:src/pull/trace.c
Location:line 1171, column 7
Description:Potential leak of memory pointed to by 'trc'

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
25#include "pull.h"
26#include "privatePull.h"
27
28pullTrace *
29pullTraceNew(void) {
30 pullTrace *ret;
31
32 ret = AIR_CALLOC(1, pullTrace)(pullTrace*)(calloc((1), sizeof(pullTrace)));
11
←
Within the expansion of the macro 'AIR_CALLOC':
→
a
Memory is allocated
33 if (ret) {
12
←
Assuming 'ret' is non-null
→
13
←
Taking true branch
→
34 ret->seedPos[0] = ret->seedPos[1] = AIR_NAN(airFloatQNaN.f);
35 ret->seedPos[2] = ret->seedPos[3] = AIR_NAN(airFloatQNaN.f);
36 ret->nvert = nrrdNew();
37 ret->nstrn = nrrdNew();
38 ret->nstab = nrrdNew();
39 ret->norin = nrrdNew();
40 ret->seedIdx = 0;
41 ret->whyStop[0] = ret->whyStop[1] = pullTraceStopUnknown;
42 ret->whyStopCFail[0] = ret->whyStopCFail[1] = pullConstraintFailUnknown;
43 ret->whyNowhere = pullTraceStopUnknown;
44 ret->whyNowhereCFail = pullConstraintFailUnknown;
45 }
46 return ret;
47}
48
49pullTrace *
50pullTraceNix(pullTrace *pts) {
51
52 if (pts) {
53 nrrdNuke(pts->nvert);
54 nrrdNuke(pts->nstrn);
55 nrrdNuke(pts->nstab);
56 nrrdNuke(pts->norin);
57 free(pts);
58 }
59 return NULL((void*)0);
60}
61
62void
63pullTraceStability(double *spcStab,
64 double *oriStab,
65 const double pos0[4],
66 const double pos1[4],
67 const double ori0[3],
68 const double ori1[3],
69 double sigma0,
70 const pullContext *pctx) {
71 double sc, stb, dx, ds, diff[4];
72
73 ELL_4V_SUB(diff, pos1, pos0)((diff)[0] = (pos1)[0] - (pos0)[0], (diff)[1] = (pos1)[1] - (
pos0)[1], (diff)[2] = (pos1)[2] - (pos0)[2], (diff)[3] = (pos1
)[3] - (pos0)[3]);
74 dx = ELL_3V_LEN(diff)(sqrt((((diff))[0]*((diff))[0] + ((diff))[1]*((diff))[1] + ((
diff))[2]*((diff))[2])))/(pctx->voxelSizeSpace);
75 sc = (pos0[3] + pos1[3])/2;
76 if (pctx->flag.scaleIsTau) {
77 sc = gageSigOfTau(sc)airSigmaOfTau(sc);
78 }
79 sc += sigma0;
80 dx /= sc;
81 ds = diff[3];
82 stb = atan2(ds, dx)/(AIR_PI3.14159265358979323846/2); /* [-1,1]: 0 means least stable */
83 stb = AIR_ABS(stb)((stb) > 0.0f ? (stb) : -(stb)); /* [0,1]: 0 means least stable */
84 *spcStab = AIR_CLAMP(0, stb, 1)((stb) < (0) ? (0) : ((stb) > (1) ? (1) : (stb)));
85 /*
86 static double maxvelo = 0;
87 if (1 && vv > maxvelo) {
88 char me[]="pullTraceStability";
89 printf("\n%s: [%g,%g,%g,%g] - [%g,%g,%g,%g] = [%g,%g,%g,%g]\n", me,
90 pos0[0], pos0[1], pos0[2], pos0[3],
91 pos1[0], pos1[1], pos1[2], pos1[3],
92 diff[0], diff[1], diff[2], diff[3]);
93 printf("%s: dx = %g -> %g; ds = %g\n", me, dx0, dx, ds);
94 printf("%s: vv = atan2(%g,%g)/pi = %g -> %g > %g\n", me, ds, dx, vv0, vv, maxvelo);
95 maxvelo = vv;
96 }
97 */
98 if (ori0 && ori1) {
99 /* dori = delta in orientation */
100 double dori = ell_3v_angle_d(ori0, ori1);
101 /* dori in [0,pi]; 0 and pi mean no change */
102 if (dori > AIR_PI3.14159265358979323846/2) {
103 dori -= AIR_PI3.14159265358979323846;
104 }
105 /* dori in [-pi/2,pi/2]; 0 means no change; +-pi/2 means most change */
106 dori = AIR_ABS(dori)((dori) > 0.0f ? (dori) : -(dori));
107 /* dori in [0,pi/2]; 0 means no change; pi/2 means most change */
108 *oriStab = atan2(ds, dori)/(AIR_PI3.14159265358979323846/2);
109 /* *oriStab in [0,1]; 0 means stable, 1 means not stable */
110 } else {
111 *oriStab = 0;
112 }
113 return;
114}
115
116int
117_pullConstrTanSlap(pullContext *pctx, pullPoint *point,
118 double tlen,
119 /* input+output */ double toff[3],
120 /* output */ int *constrFailP) {
121 static const char me[]="_pullConstrTanSlap";
122 double pos0[4], tt;
123
124 if (!(pctx && point && toff && constrFailP)) {
125 biffAddf(PULLpullBiffKey, "%s: got NULL pointer", me);
126 return 1;
127 }
128 if (pctx->flag.zeroZ) {
129 toff[2] = 0;
130 }
131 ELL_3V_NORM(toff, toff, tt)(tt = (sqrt((((toff))[0]*((toff))[0] + ((toff))[1]*((toff))[1
] + ((toff))[2]*((toff))[2]))), ((toff)[0] = (1.0/tt)*(toff)[
0], (toff)[1] = (1.0/tt)*(toff)[1], (toff)[2] = (1.0/tt)*(toff
)[2]));
132 if (!(tlen > 0 && tt > 0)) {
133 biffAddf(PULLpullBiffKey, "%s: tlen %g or |toff| %g not positive", me, tlen, tt);
134 return 1;
135 }
136 ELL_4V_COPY(pos0, point->pos)((pos0)[0] = (point->pos)[0], (pos0)[1] = (point->pos)[
1], (pos0)[2] = (point->pos)[2], (pos0)[3] = (point->pos
)[3]); /* pos0[3] should stay put; will test at end */
137
138 ELL_3V_SCALE_ADD2(point->pos, 1, pos0, tlen, toff)((point->pos)[0] = (1)*(pos0)[0] + (tlen)*(toff)[0], (point
->pos)[1] = (1)*(pos0)[1] + (tlen)*(toff)[1], (point->pos
)[2] = (1)*(pos0)[2] + (tlen)*(toff)[2]);
139 if (_pullConstraintSatisfy(pctx->task[0], point, 0 /* travmax */, constrFailP)) {
140 biffAddf(PULLpullBiffKey, "%s: trouble", me);
141 ELL_4V_COPY(point->pos, pos0)((point->pos)[0] = (pos0)[0], (point->pos)[1] = (pos0)[
1], (point->pos)[2] = (pos0)[2], (point->pos)[3] = (pos0
)[3]); return 1;
142 }
143 /* save offset to new position */
144 ELL_3V_SUB(toff, point->pos, pos0)((toff)[0] = (point->pos)[0] - (pos0)[0], (toff)[1] = (point
->pos)[1] - (pos0)[1], (toff)[2] = (point->pos)[2] - (pos0
)[2]);
145
146 if (pos0[3] != point->pos[3]) {
147 biffAddf(PULLpullBiffKey, "%s: point->pos[3] %g was changed (from %g)",
148 me, point->pos[3], pos0[3]);
149 ELL_4V_COPY(point->pos, pos0)((point->pos)[0] = (pos0)[0], (point->pos)[1] = (pos0)[
1], (point->pos)[2] = (pos0)[2], (point->pos)[3] = (pos0
)[3]); return 1;
150 }
151 ELL_4V_COPY(point->pos, pos0)((point->pos)[0] = (pos0)[0], (point->pos)[1] = (pos0)[
1], (point->pos)[2] = (pos0)[2], (point->pos)[3] = (pos0
)[3]);
152 return 0;
153}
154
155int
156_pullConstrOrientFind(pullContext *pctx, pullPoint *point,
157 int normalfind, /* find normal two 2D surface,
158 else find tangent to 1D curve */
159 double tlen,
160 const double tdir0[3], /* if non-NULL, try using this direction */
161 /* output */
162 double dir[3],
163 int *constrFailP) {
164 static const char me[]="_pullConstrOrientFind";
165 double tt;
166
167#define SLAP(LEN, DIR)if (_pullConstrTanSlap(pctx, point, (LEN), (DIR), constrFailP
)) { biffAddf(pullBiffKey, "%s: looking for tangent, starting with (%g,%g,%g)"
, me, (DIR)[0], (DIR)[1], (DIR)[2]); return 1; } if (*constrFailP
) { return 0; } \
168 /* fprintf(stderr, "!%s: SLAP %g %g %g -->", me, (DIR)[0], (DIR)[1], (DIR)[2]); */ \
169 if (_pullConstrTanSlap(pctx, point, (LEN), (DIR), constrFailP)) { \
170 biffAddf(PULLpullBiffKey, "%s: looking for tangent, starting with (%g,%g,%g)", \
171 me, (DIR)[0], (DIR)[1], (DIR)[2]); \
172 return 1; \
173 } \
174 if (*constrFailP) { \
175 /* unsuccessful in finding tangent, but not a biff error */ \
176 return 0; \
177 } \
178 /* fprintf(stderr, " %g %g %g\n", (DIR)[0], (DIR)[1], (DIR)[2]); */
179
180 if (!(pctx && point && constrFailP)) {
181 biffAddf(PULLpullBiffKey, "%s: got NULL pointer", me);
182 return 1;
183 }
184 *constrFailP = pullConstraintFailUnknown;
185 if (normalfind) {
186 double tan0[3], tan1[3];
187 /* looking for a surface normal */
188 if (tdir0) { /* have something to start with */
189 ell_3v_perp_d(tan0, tdir0);
190 ELL_3V_CROSS(tan1, tan0, tdir0)((tan1)[0] = (tan0)[1]*(tdir0)[2] - (tan0)[2]*(tdir0)[1], (tan1
)[1] = (tan0)[2]*(tdir0)[0] - (tan0)[0]*(tdir0)[2], (tan1)[2]
= (tan0)[0]*(tdir0)[1] - (tan0)[1]*(tdir0)[0]);
191 SLAP(tlen, tan1)if (_pullConstrTanSlap(pctx, point, (tlen), (tan1), constrFailP
)) { biffAddf(pullBiffKey, "%s: looking for tangent, starting with (%g,%g,%g)"
, me, (tan1)[0], (tan1)[1], (tan1)[2]); return 1; } if (*constrFailP
) { return 0; };
192 SLAP(tlen, tan0)if (_pullConstrTanSlap(pctx, point, (tlen), (tan0), constrFailP
)) { biffAddf(pullBiffKey, "%s: looking for tangent, starting with (%g,%g,%g)"
, me, (tan0)[0], (tan0)[1], (tan0)[2]); return 1; } if (*constrFailP
) { return 0; };
193 ELL_3V_CROSS(dir, tan1, tan0)((dir)[0] = (tan1)[1]*(tan0)[2] - (tan1)[2]*(tan0)[1], (dir)[
1] = (tan1)[2]*(tan0)[0] - (tan1)[0]*(tan0)[2], (dir)[2] = (tan1
)[0]*(tan0)[1] - (tan1)[1]*(tan0)[0]);
194 } else { /* have to start from scratch */
195 double tns[9] = {1,0,0, 0,1,0, 0,0,1};
196 SLAP(tlen, tns+0)if (_pullConstrTanSlap(pctx, point, (tlen), (tns+0), constrFailP
)) { biffAddf(pullBiffKey, "%s: looking for tangent, starting with (%g,%g,%g)"
, me, (tns+0)[0], (tns+0)[1], (tns+0)[2]); return 1; } if (*constrFailP
) { return 0; };
197 SLAP(tlen, tns+3)if (_pullConstrTanSlap(pctx, point, (tlen), (tns+3), constrFailP
)) { biffAddf(pullBiffKey, "%s: looking for tangent, starting with (%g,%g,%g)"
, me, (tns+3)[0], (tns+3)[1], (tns+3)[2]); return 1; } if (*constrFailP
) { return 0; };
198 SLAP(tlen, tns+6)if (_pullConstrTanSlap(pctx, point, (tlen), (tns+6), constrFailP
)) { biffAddf(pullBiffKey, "%s: looking for tangent, starting with (%g,%g,%g)"
, me, (tns+6)[0], (tns+6)[1], (tns+6)[2]); return 1; } if (*constrFailP
) { return 0; };
199 ell_3m_1d_nullspace_d(dir, tns);
200 }
201 } else {
202 /* looking for a curve tangent */
203 if (tdir0) { /* have something to start with */
204 ELL_3V_COPY(dir, tdir0)((dir)[0] = (tdir0)[0], (dir)[1] = (tdir0)[1], (dir)[2] = (tdir0
)[2]);
205 SLAP(tlen, dir)if (_pullConstrTanSlap(pctx, point, (tlen), (dir), constrFailP
)) { biffAddf(pullBiffKey, "%s: looking for tangent, starting with (%g,%g,%g)"
, me, (dir)[0], (dir)[1], (dir)[2]); return 1; } if (*constrFailP
) { return 0; };
206 /* SLAP(tlen, dir); (didn't have much effect, in one test) */
207 } else { /* have to start from scratch */
208 double tX[3] = {1,0,0}, tY[3] = {0,1,0}, tZ[3] = {0,0,1};
209 SLAP(tlen, tX)if (_pullConstrTanSlap(pctx, point, (tlen), (tX), constrFailP
)) { biffAddf(pullBiffKey, "%s: looking for tangent, starting with (%g,%g,%g)"
, me, (tX)[0], (tX)[1], (tX)[2]); return 1; } if (*constrFailP
) { return 0; };
210 SLAP(tlen, tY)if (_pullConstrTanSlap(pctx, point, (tlen), (tY), constrFailP
)) { biffAddf(pullBiffKey, "%s: looking for tangent, starting with (%g,%g,%g)"
, me, (tY)[0], (tY)[1], (tY)[2]); return 1; } if (*constrFailP
) { return 0; };
211 if (ELL_3V_DOT(tX, tY)((tX)[0]*(tY)[0] + (tX)[1]*(tY)[1] + (tX)[2]*(tY)[2]) < 0) { ELL_3V_SCALE(tY, -1, tY)((tY)[0] = (-1)*(tY)[0], (tY)[1] = (-1)*(tY)[1], (tY)[2] = (-
1)*(tY)[2]); }
212 if (pctx->flag.zeroZ) {
213 ELL_3V_SET(tZ, 0, 0, 0)((tZ)[0] = (0), (tZ)[1] = (0), (tZ)[2] = (0));
214 } else {
215 SLAP(tlen, tZ)if (_pullConstrTanSlap(pctx, point, (tlen), (tZ), constrFailP
)) { biffAddf(pullBiffKey, "%s: looking for tangent, starting with (%g,%g,%g)"
, me, (tZ)[0], (tZ)[1], (tZ)[2]); return 1; } if (*constrFailP
) { return 0; };
216 if (ELL_3V_DOT(tX, tZ)((tX)[0]*(tZ)[0] + (tX)[1]*(tZ)[1] + (tX)[2]*(tZ)[2]) < 0) { ELL_3V_SCALE(tZ, -1, tZ)((tZ)[0] = (-1)*(tZ)[0], (tZ)[1] = (-1)*(tZ)[1], (tZ)[2] = (-
1)*(tZ)[2]); }
217 if (ELL_3V_DOT(tY, tZ)((tY)[0]*(tZ)[0] + (tY)[1]*(tZ)[1] + (tY)[2]*(tZ)[2]) < 0) { ELL_3V_SCALE(tY, -1, tZ)((tY)[0] = (-1)*(tZ)[0], (tY)[1] = (-1)*(tZ)[1], (tY)[2] = (-
1)*(tZ)[2]); }
218 }
219 ELL_3V_ADD3(dir, tX, tY, tZ)((dir)[0] = (tX)[0] + (tY)[0] + (tZ)[0], (dir)[1] = (tX)[1] +
(tY)[1] + (tZ)[1], (dir)[2] = (tX)[2] + (tY)[2] + (tZ)[2]);
220 }
221 }
222 ELL_3V_NORM(dir, dir, tt)(tt = (sqrt((((dir))[0]*((dir))[0] + ((dir))[1]*((dir))[1] + (
(dir))[2]*((dir))[2]))), ((dir)[0] = (1.0/tt)*(dir)[0], (dir)
[1] = (1.0/tt)*(dir)[1], (dir)[2] = (1.0/tt)*(dir)[2]));
223 if (!(tt > 0)) {
224 biffAddf(PULLpullBiffKey, "%s: computed direction is zero (%g)?", me, tt);
225 return 1;
226 }
227 return 0;
228}
229
230/*
231******** pullTraceSet
232**
233** computes a single trace, according to the given parameters,
234** and store it in the pullTrace
235int recordStrength: should strength be recorded along trace
236double scaleDelta: discrete step along scale
237double halfScaleWin: how far, along scale, trace should go in each direction
238unsigned int arrIncr: increment for storing position (maybe strength)
239const double _seedPos[4]: starting position
240*/
241int
242pullTraceSet(pullContext *pctx, pullTrace *pts,
243 int recordStrength,
244 double scaleDelta, double halfScaleWin,
245 double orientTestLen,
246 unsigned int arrIncr,
247 const double _seedPos[4]) {
248 static const char me[]="pullTraceSet";
249 pullPoint *point;
250 airArray *mop, *trceArr[2], *hstrnArr[2], *horinArr[2];
251 double *trce[2], ssrange[2], *vert, *hstrn[2], *horin[2], *strn, *orin, *stab,
252 seedPos[4], polen, porin[3];
253 int constrFail;
254 unsigned int dirIdx, lentmp, tidx, oidx, vertNum;
255
256 if (!( pctx && pts && _seedPos )) {
257 biffAddf(PULLpullBiffKey, "%s: got NULL pointer", me);
258 return 1;
259 }
260 if (!( AIR_EXISTS(scaleDelta)(((int)(!((scaleDelta) - (scaleDelta))))) && scaleDelta > 0.0 )) {
261 biffAddf(PULLpullBiffKey, "%s: need existing scaleDelta > 0 (not %g)",
262 me, scaleDelta);
263 return 1;
264 }
265 if (!( halfScaleWin > 0 )) {
266 biffAddf(PULLpullBiffKey, "%s: need halfScaleWin > 0", me);
267 return 1;
268 }
269 if (!(pctx->constraint)) {
270 biffAddf(PULLpullBiffKey, "%s: given context doesn't have constraint set", me);
271 return 1;
272 }
273 pts->fdim = _pullConstraintDim(pctx);
274 if (0 > pts->fdim) {
275 biffAddf(PULLpullBiffKey, "%s: couldn't learn dimension of feature", me);
276 return 1;
277 }
278 if (pts->fdim == 2 && pctx->flag.zeroZ) {
279 biffAddf(PULLpullBiffKey, "%s: can't have feature dim 2 with zeroZ", me);
280 return 1;
281 }
282 if (!( AIR_EXISTS(orientTestLen)(((int)(!((orientTestLen) - (orientTestLen))))) && orientTestLen >= 0 )) {
283 biffAddf(PULLpullBiffKey, "%s: need non-negative orientTestLen (not %g)\n",
284 me, orientTestLen);
285 return 1;
286 }
287 if (pts->fdim && !orientTestLen) {
288 /* not really an error */
289 /*
290 fprintf(stderr, "\n\n%s: WARNING: have a %d-D feature, but not "
291 "measuring its orientation\n\n\n", me, pts->fdim);
292 */
293 }
294 /* The test for "should I measure orientation"
295 is "if (pts->fdim && orientTestLen)" */
296 if (recordStrength && !pctx->ispec[pullInfoStrength]) {
297 biffAddf(PULLpullBiffKey, "%s: want to record strength but %s not set in context",
298 me, airEnumStr(pullInfo, pullInfoStrength));
299 return 1;
300 }
301 if (pullConstraintScaleRange(pctx, ssrange)) {
302 biffAddf(PULLpullBiffKey, "%s: trouble getting scale range", me);
303 return 1;
304 }
305
306 /* re-initialize termination descriptions (in case of trace re-use) */
307 pts->whyStop[0] = pts->whyStop[1] = pullTraceStopUnknown;
308 pts->whyStopCFail[0] = pts->whyStopCFail[1] = pullConstraintFailUnknown;
309 pts->whyNowhere = pullTraceStopUnknown;
310 pts->whyNowhereCFail = pullConstraintFailUnknown;
311
312 /* enforce zeroZ */
313 ELL_4V_COPY(seedPos, _seedPos)((seedPos)[0] = (_seedPos)[0], (seedPos)[1] = (_seedPos)[1], (
seedPos)[2] = (_seedPos)[2], (seedPos)[3] = (_seedPos)[3]);
314 if (pctx->flag.zeroZ) {
315 seedPos[2] = 0.0;
316 }
317
318 /* save seedPos in any case */
319 ELL_4V_COPY(pts->seedPos, seedPos)((pts->seedPos)[0] = (seedPos)[0], (pts->seedPos)[1] = (
seedPos)[1], (pts->seedPos)[2] = (seedPos)[2], (pts->seedPos
)[3] = (seedPos)[3]);
320
321 mop = airMopNew();
322 point = pullPointNew(pctx); /* we'll want to decrement idtagNext later */
323 airMopAdd(mop, point, (airMopper)pullPointNix, airMopAlways);
324
325 ELL_4V_COPY(point->pos, seedPos)((point->pos)[0] = (seedPos)[0], (point->pos)[1] = (seedPos
)[1], (point->pos)[2] = (seedPos)[2], (point->pos)[3] =
(seedPos)[3]);
326 /*
327 if (pctx->verbose) {
328 fprintf(stderr, "%s: trying at seed=(%g,%g,%g,%g)\n", me,
329 seedPos[0], seedPos[1], seedPos[2], seedPos[3]);
330 }
331 */
332 /* The termination of the trace due to low stablity is handled here
333 (or could be), not by _pullConstraintSatisfy, so we set 0 for the
334 travelMax arg of _pullConstraintSatisfy (no travel limit) */
335 if (_pullConstraintSatisfy(pctx->task[0], point,
336 0 /* travmax */, &constrFail)) {
337 biffAddf(PULLpullBiffKey, "%s: constraint sat on seed point", me);
338 airMopError(mop);
339 return 1;
340 }
341 /*
342 if (pctx->verbose) {
343 fprintf(stderr, "%s: seed=(%g,%g,%g,%g) -> %s (%g,%g,%g,%g)\n", me,
344 seedPos[0], seedPos[1], seedPos[2], seedPos[3],
345 constrFail ? "!NO!" : "(yes)",
346 point->pos[0] - seedPos[0], point->pos[1] - seedPos[1],
347 point->pos[2] - seedPos[2], point->pos[3] - seedPos[3]);
348 }
349 */
350 if (constrFail) {
351 pts->whyNowhere = pullTraceStopConstrFail;
352 airMopOkay(mop);
353 /* pctx->idtagNext -= 1; / * HACK * / */
354 return 0;
355 }
356 if (point->status & PULL_STATUS_EDGE_BIT(1<< 4)) {
357 pts->whyNowhere = pullTraceStopVolumeEdge;
358 airMopOkay(mop);
359 /* pctx->idtagNext -= 1; / * HACK * / */
360 return 0;
361 }
362 if (pctx->flag.zeroZ && point->pos[2] != 0) {
363 biffAddf(PULLpullBiffKey, "%s: zeroZ violated (a)", me);
364 airMopError(mop);
365 return 1;
366 }
367
368 /* else constraint sat worked at seed point; we have work to do */
369 if (pts->fdim && orientTestLen) {
370 /* learn orientation at seed point */
371 polen = (orientTestLen
372 *pctx->voxelSizeSpace
373 /* if used, the effect of this this last (unprincipled) factor
374 is gradually increase the test distance with scale
375 *(1 + gageTauOfSig(_pullSigma(pctx, point->pos))) */ );
376 double pos0[4], dp[4];
377 int cf;
378 ELL_4V_COPY(pos0, point->pos)((pos0)[0] = (point->pos)[0], (pos0)[1] = (point->pos)[
1], (pos0)[2] = (point->pos)[2], (pos0)[3] = (point->pos
)[3]);
379
380 if (_pullConstrOrientFind(pctx, point, pts->fdim == 2,
381 polen, NULL((void*)0) /* no initial guess */, porin, &cf)) {
382 biffAddf(PULLpullBiffKey, "%s: trying to find orientation at seed", me);
383 airMopError(mop);
384 return 1;
385 }
386 ELL_4V_SUB(dp, pos0, point->pos)((dp)[0] = (pos0)[0] - (point->pos)[0], (dp)[1] = (pos0)[1
] - (point->pos)[1], (dp)[2] = (pos0)[2] - (point->pos)
[2], (dp)[3] = (pos0)[3] - (point->pos)[3]);
387 /*
388 fprintf(stderr, "!%s: cf = %d (%s)\n", me, cf, airEnumStr(pullConstraintFail, cf));
389 fprintf(stderr, "!%s: (fdim=%u) pos=[%g,%g,%g,%g] polen=%g porin=[%g,%g,%g] |%g|\n",
390 me, pts->fdim,
391 point->pos[0], point->pos[1], point->pos[2], point->pos[3],
392 polen, porin[0], porin[1], porin[2], ELL_4V_LEN(dp));
393 */
394 if (cf) {
395 pts->whyNowhere = pullTraceStopConstrFail;
396 pts->whyNowhereCFail = cf;
397 airMopOkay(mop);
398 /* pctx->idtagNext -= 1; / * HACK * / */
399 return 0;
400 }
401 } else {
402 /* either feature is 0D points, or don't care about orientation */
403 polen = AIR_NAN(airFloatQNaN.f);
404 ELL_3V_SET(porin, AIR_NAN, AIR_NAN, AIR_NAN)((porin)[0] = ((airFloatQNaN.f)), (porin)[1] = ((airFloatQNaN
.f)), (porin)[2] = ((airFloatQNaN.f)));
405 }
406
407 for (dirIdx=0; dirIdx<2; dirIdx++) {
408 trceArr[dirIdx] = airArrayNew((void**)(trce + dirIdx), NULL((void*)0),
409 4*sizeof(double), arrIncr);
410 airMopAdd(mop, trceArr[dirIdx], (airMopper)airArrayNuke, airMopAlways);
411 if (recordStrength) {
412 hstrnArr[dirIdx] = airArrayNew((void**)(hstrn + dirIdx), NULL((void*)0),
413 sizeof(double), arrIncr);
414 airMopAdd(mop, hstrnArr[dirIdx], (airMopper)airArrayNuke, airMopAlways);
415 } else {
416 hstrnArr[dirIdx] = NULL((void*)0);
417 hstrn[dirIdx] = NULL((void*)0);
418 }
419 if (pts->fdim && orientTestLen) {
420 horinArr[dirIdx] = airArrayNew((void**)(horin + dirIdx), NULL((void*)0),
421 3*sizeof(double), arrIncr);
422 airMopAdd(mop, horinArr[dirIdx], (airMopper)airArrayNuke, airMopAlways);
423 } else {
424 horinArr[dirIdx] = NULL((void*)0);
425 horin[dirIdx] = NULL((void*)0);
426 }
427 }
428 for (dirIdx=0; dirIdx<2; dirIdx++) {
429 unsigned int step;
430 double dscl;
431 dscl = (!dirIdx ? -1 : +1)*scaleDelta;
432 step = 0;
433 while (1) {
434 if (!step) {
435 /* first step in both directions requires special tricks */
436 if (0 == dirIdx) {
437 /* this is done once, at the very start */
438 /* save constraint sat of seed point */
439 tidx = airArrayLenIncr(trceArr[0], 1);
440 ELL_4V_COPY(trce[0] + 4*tidx, point->pos)((trce[0] + 4*tidx)[0] = (point->pos)[0], (trce[0] + 4*tidx
)[1] = (point->pos)[1], (trce[0] + 4*tidx)[2] = (point->
pos)[2], (trce[0] + 4*tidx)[3] = (point->pos)[3]);
441 if (recordStrength) {
442 airArrayLenIncr(hstrnArr[0], 1);
443 hstrn[0][0] = pullPointScalar(pctx, point, pullInfoStrength,
444 NULL((void*)0), NULL((void*)0));
445 }
446 if (pts->fdim && orientTestLen) {
447 airArrayLenIncr(horinArr[0], 1);
448 ELL_3V_COPY(horin[0] + 3*0, porin)((horin[0] + 3*0)[0] = (porin)[0], (horin[0] + 3*0)[1] = (porin
)[1], (horin[0] + 3*0)[2] = (porin)[2]);
449 }
450 } else {
451 /* re-set position from constraint sat of seed pos */
452 ELL_4V_COPY(point->pos, trce[0] + 4*0)((point->pos)[0] = (trce[0] + 4*0)[0], (point->pos)[1] =
(trce[0] + 4*0)[1], (point->pos)[2] = (trce[0] + 4*0)[2],
(point->pos)[3] = (trce[0] + 4*0)[3]);
453 if (pts->fdim && orientTestLen) {
454 ELL_3V_COPY(porin, horin[0] + 3*0)((porin)[0] = (horin[0] + 3*0)[0], (porin)[1] = (horin[0] + 3
*0)[1], (porin)[2] = (horin[0] + 3*0)[2]);
455 }
456 }
457 }
458 /* nudge position along scale */
459 point->pos[3] += dscl;
460 if (!AIR_IN_OP(ssrange[0], point->pos[3], ssrange[1])((ssrange[0]) < (point->pos[3]) && (point->pos
[3]) < (ssrange[1]))) {
461 /* if we've stepped outside the range of scale for the volume
462 containing the constraint manifold, we're done */
463 pts->whyStop[dirIdx] = pullTraceStopBounds;
464 break;
465 }
466 if (AIR_ABS(point->pos[3] - seedPos[3])((point->pos[3] - seedPos[3]) > 0.0f ? (point->pos[3
] - seedPos[3]) : -(point->pos[3] - seedPos[3])) > halfScaleWin) {
467 /* we've moved along scale as far as allowed */
468 pts->whyStop[dirIdx] = pullTraceStopLength;
469 break;
470 }
471 /* re-assert constraint */
472 /*
473 fprintf(stderr, "%s(%u): pos = %g %g %g %g.... \n", me,
474 point->idtag, point->pos[0], point->pos[1],
475 point->pos[2], point->pos[3]);
476 */
477 if (_pullConstraintSatisfy(pctx->task[0], point,
478 0 /* travmax */, &constrFail)) {
479 biffAddf(PULLpullBiffKey, "%s: dir %u, step %u", me, dirIdx, step);
480 airMopError(mop);
481 return 1;
482 }
483 /*
484 if (pctx->verbose) {
485 fprintf(stderr, "%s(%u): ... %s(%d); pos = %g %g %g %g\n", me,
486 point->idtag,
487 constrFail ? "FAIL" : "(ok)",
488 constrFail, point->pos[0], point->pos[1],
489 point->pos[2], point->pos[3]);
490 }
491 */
492 if (point->status & PULL_STATUS_EDGE_BIT(1<< 4)) {
493 pts->whyStop[dirIdx] = pullTraceStopVolumeEdge;
494 break;
495 }
496 if (constrFail) {
497 /* constraint sat failed; no error, we're just done
498 with stepping for this direction */
499 pts->whyStop[dirIdx] = pullTraceStopConstrFail;
500 pts->whyStopCFail[dirIdx] = constrFail;
501 break;
502 }
503 if (pctx->flag.zeroZ && point->pos[2] != 0) {
504 biffAddf(PULLpullBiffKey, "%s: zeroZ violated (b)", me);
505 airMopError(mop);
506 return 1;
507 }
508 if (pts->fdim && orientTestLen) {
509 if (_pullConstrOrientFind(pctx, point, pts->fdim == 2,
510 polen, porin, porin, &constrFail)) {
511 biffAddf(PULLpullBiffKey, "%s: at dir %u, step %u", me, dirIdx, step);
512 airMopError(mop);
513 return 1;
514 }
515 }
516 if (trceArr[dirIdx]->len >= 2) {
517 /* see if we're moving too fast, by comparing with previous point */
518 /* actually, screw that */
519 }
520 /* else save new point on trace */
521 tidx = airArrayLenIncr(trceArr[dirIdx], 1);
522 ELL_4V_COPY(trce[dirIdx] + 4*tidx, point->pos)((trce[dirIdx] + 4*tidx)[0] = (point->pos)[0], (trce[dirIdx
] + 4*tidx)[1] = (point->pos)[1], (trce[dirIdx] + 4*tidx)[
2] = (point->pos)[2], (trce[dirIdx] + 4*tidx)[3] = (point->
pos)[3]);
523 if (recordStrength) {
524 tidx = airArrayLenIncr(hstrnArr[dirIdx], 1);
525 hstrn[dirIdx][tidx] = pullPointScalar(pctx, point, pullInfoStrength,
526 NULL((void*)0), NULL((void*)0));
527 }
528 if (pts->fdim && orientTestLen) {
529 tidx = airArrayLenIncr(horinArr[dirIdx], 1);
530 ELL_3V_COPY(horin[dirIdx] + 3*tidx, porin)((horin[dirIdx] + 3*tidx)[0] = (porin)[0], (horin[dirIdx] + 3
*tidx)[1] = (porin)[1], (horin[dirIdx] + 3*tidx)[2] = (porin)
[2]);
531 }
532 step++;
533 }
534 }
535
536 /* transfer trace halves to pts->nvert */
537 vertNum = trceArr[0]->len + trceArr[1]->len;
538 if (0 == vertNum || 1 == vertNum || 2 == vertNum) {
539 pts->whyNowhere = pullTraceStopStub;
540 airMopOkay(mop);
541 /* pctx->idtagNext -= 1; / * HACK * / */
542 return 0;
543 }
544
545 if (nrrdMaybeAlloc_va(pts->nvert, nrrdTypeDouble, 2,
546 AIR_CAST(size_t, 4)((size_t)(4)),
547 AIR_CAST(size_t, vertNum)((size_t)(vertNum)))
548 || nrrdMaybeAlloc_va(pts->nstab, nrrdTypeDouble, 2,
549 AIR_CAST(size_t, 2)((size_t)(2)),
550 AIR_CAST(size_t, vertNum)((size_t)(vertNum)))) {
551 biffMovef(PULLpullBiffKey, NRRDnrrdBiffKey, "%s: allocating output", me);
552 airMopError(mop);
553 return 1;
554 }
555 if (recordStrength) {
556 /* doing slicing is a simple form of allocation */
557 if (nrrdSlice(pts->nstrn, pts->nvert, 0 /* axis */, 0 /* pos */)) {
558 biffMovef(PULLpullBiffKey, NRRDnrrdBiffKey, "%s: allocating strength output", me);
559 airMopError(mop);
560 return 1;
561 }
562 }
563 if (pts->fdim && orientTestLen) {
564 /* cropping just to allocate */
565 size_t cmin[2] = {0, 0}, cmax[2] = {2, pts->nvert->axis[1].size-1};
566 if (nrrdCrop(pts->norin, pts->nvert, cmin, cmax)) {
567 biffMovef(PULLpullBiffKey, NRRDnrrdBiffKey, "%s: allocating orientation output", me);
568 airMopError(mop);
569 return 1;
570 }
571 }
572 vert = AIR_CAST(double *, pts->nvert->data)((double *)(pts->nvert->data));
573 if (recordStrength) {
574 strn = AIR_CAST(double *, pts->nstrn->data)((double *)(pts->nstrn->data));
575 } else {
576 strn = NULL((void*)0);
577 }
578 if (pts->fdim && orientTestLen) {
579 orin = AIR_CAST(double *, pts->norin->data)((double *)(pts->norin->data));
580 } else {
581 orin = NULL((void*)0);
582 }
583 stab = AIR_CAST(double *, pts->nstab->data)((double *)(pts->nstab->data));
584 lentmp = trceArr[0]->len;
585 oidx = 0;
586 for (tidx=0; tidx<lentmp; tidx++) {
587 ELL_4V_COPY(vert + 4*oidx, trce[0] + 4*(lentmp - 1 - tidx))((vert + 4*oidx)[0] = (trce[0] + 4*(lentmp - 1 - tidx))[0], (
vert + 4*oidx)[1] = (trce[0] + 4*(lentmp - 1 - tidx))[1], (vert
+ 4*oidx)[2] = (trce[0] + 4*(lentmp - 1 - tidx))[2], (vert +
4*oidx)[3] = (trce[0] + 4*(lentmp - 1 - tidx))[3]);
588 if (strn) {
589 strn[oidx] = hstrn[0][lentmp - 1 - tidx];
590 }
591 if (orin) {
592 ELL_3V_COPY(orin + 3*oidx, horin[0] + 3*(lentmp - 1 - tidx))((orin + 3*oidx)[0] = (horin[0] + 3*(lentmp - 1 - tidx))[0], (
orin + 3*oidx)[1] = (horin[0] + 3*(lentmp - 1 - tidx))[1], (orin
+ 3*oidx)[2] = (horin[0] + 3*(lentmp - 1 - tidx))[2]);
593 }
594 oidx++;
595 }
596 /* the last index written to (before oidx++) was the seed index */
597 pts->seedIdx = oidx-1;
598 lentmp = trceArr[1]->len;
599 for (tidx=0; tidx<lentmp; tidx++) {
600 ELL_4V_COPY(vert + 4*oidx, trce[1] + 4*tidx)((vert + 4*oidx)[0] = (trce[1] + 4*tidx)[0], (vert + 4*oidx)[
1] = (trce[1] + 4*tidx)[1], (vert + 4*oidx)[2] = (trce[1] + 4
*tidx)[2], (vert + 4*oidx)[3] = (trce[1] + 4*tidx)[3]);
601 if (strn) {
602 strn[oidx] = hstrn[1][tidx];
603 }
604 if (orin) {
605 ELL_3V_COPY(orin + 3*oidx, horin[1] + 3*tidx)((orin + 3*oidx)[0] = (horin[1] + 3*tidx)[0], (orin + 3*oidx)
[1] = (horin[1] + 3*tidx)[1], (orin + 3*oidx)[2] = (horin[1] +
3*tidx)[2]);
606 }
607 oidx++;
608 }
609 lentmp = pts->nstab->axis[1].size;
610 if (1 == lentmp) {
611 stab[0 + 2*0] = 0.0;
612 stab[1 + 2*0] = 0.0;
613 } else {
614 for (tidx=0; tidx<lentmp; tidx++) {
615 double *pA, *pB, *p0, *p1, *p2, *rA, *rB, *r0=NULL((void*)0), *r1=NULL((void*)0), *r2=NULL((void*)0);
616 p0 = vert + 4*(tidx-1);
617 p1 = vert + 4*tidx;
618 p2 = vert + 4*(tidx+1);
619 if (orin) {
620 r0 = orin + 3*(tidx-1);
621 r1 = orin + 3*tidx;
622 r2 = orin + 3*(tidx+1);
623 }
624 if (!tidx) {
625 /* first */
626 pA = p1; rA = r1;
627 pB = p2; rB = r2;
628 } else if (tidx < lentmp-1) {
629 /* middle */
630 pA = p0; rA = r0;
631 pB = p2; rB = r2;
632 } else {
633 /* last */
634 pA = p0; rA = r0;
635 pB = p1; rB = r1;
636 }
637 pullTraceStability(stab + 0 + 2*tidx, stab + 1 + 2*tidx,
638 pA, pB, rA, rB, 0.5 /* sigma0 */, pctx);
639 }
640 }
641
642 airMopOkay(mop);
643 /* pctx->idtagNext -= 1; / * HACK * / */
644 return 0;
645}
646
647typedef union {
648 pullTrace ***trace;
649 void **v;
650} blahblahUnion;
651
652pullTraceMulti *
653pullTraceMultiNew(void) {
654 /* static const char me[]="pullTraceMultiNew"; */
655 pullTraceMulti *ret;
656 blahblahUnion bbu;
657
658 ret = AIR_CALLOC(1, pullTraceMulti)(pullTraceMulti*)(calloc((1), sizeof(pullTraceMulti)));
659 if (ret) {
660 ret->trace = NULL((void*)0);
661 ret->traceNum = 0;
662 ret->traceArr = airArrayNew((bbu.trace = &(ret->trace), bbu.v),
663 &(ret->traceNum), sizeof(pullTrace*),
664 _PULL_TRACE_MULTI_INCR1024);
665 airArrayPointerCB(ret->traceArr,
666 NULL((void*)0), /* because we get handed pullTrace structs
667 that have already been allocated
668 (and then we own them) */
669 (void *(*)(void *))pullTraceNix);
670 }
671 return ret;
672}
673
674int
675pullTraceMultiAdd(pullTraceMulti *mtrc, pullTrace *trc, int *addedP) {
676 static const char me[]="pullTraceMultiAdd";
677 unsigned int indx;
678
679 if (!(mtrc && trc && addedP)) {
680 biffAddf(PULLpullBiffKey, "%s: got NULL pointer", me);
681 return 1;
682 }
683 if (!(trc->nvert->data && trc->nvert->axis[1].size >= 3)) {
684 /* for now getting a stub trace is not an error
685 biffAddf(PULL, "%s: got stub trace", me);
686 return 1; */
687 *addedP = AIR_FALSE0;
688 return 0;
689 }
690 if (!(trc->nstab->data
691 && trc->nstab->axis[1].size == trc->nvert->axis[1].size)) {
692 biffAddf(PULLpullBiffKey, "%s: stab data inconsistent", me);
693 return 1;
694 }
695 *addedP = AIR_TRUE1;
696 indx = airArrayLenIncr(mtrc->traceArr, 1);
697 if (!mtrc->trace) {
698 biffAddf(PULLpullBiffKey, "%s: alloc error", me);
699 return 1;
700 }
701 mtrc->trace[indx] = trc;
702 return 0;
703}
704
705int
706pullTraceMultiPlotAdd(Nrrd *nplot, const pullTraceMulti *mtrc,
707 const Nrrd *nfilt, int strengthUse,
708 int smooth, int flatWght,
709 unsigned int trcIdxMin, unsigned int trcNum,
710 Nrrd *nmaskedpos,
711 const Nrrd *nmask) {
712 static const char me[]="pullTraceMultiPlotAdd";
713 double ssRange[2], vRange[2], *plot;
714 unsigned int sizeS, sizeT, trcIdx, trcIdxMax;
715 int *filt;
716 airArray *mop;
717 Nrrd *nsmst; /* smoothed stability */
718 airArray *mposArr;
719 double *mpos;
720 unsigned char *mask;
721
722 if (!(nplot && mtrc)) {
723 biffAddf(PULLpullBiffKey, "%s: got NULL pointer", me);
724 return 1;
725 }
726 if (nrrdCheck(nplot)) {
727 biffMovef(PULLpullBiffKey, NRRDnrrdBiffKey, "%s: trouble with nplot", me);
728 return 1;
729 }
730 if (nfilt) {
731 if (nrrdCheck(nfilt)) {
732 biffMovef(PULLpullBiffKey, NRRDnrrdBiffKey, "%s: trouble with nfilt", me);
733 return 1;
734 }
735 if (!(1 == nfilt->dim && nrrdTypeInt == nfilt->type)) {
736 biffAddf(PULLpullBiffKey, "%s: didn't get 1-D array of %s (got %u-D of %s)", me,
737 airEnumStr(nrrdType, nrrdTypeInt), nfilt->dim,
738 airEnumStr(nrrdType, nfilt->type));
739 return 1;
740 }
741 }
742 if (!(2 == nplot->dim && nrrdTypeDouble == nplot->type)) {
743 biffAddf(PULLpullBiffKey, "%s: didn't get 2-D array of %s (got %u-D of %s)", me,
744 airEnumStr(nrrdType, nrrdTypeDouble), nplot->dim,
745 airEnumStr(nrrdType, nplot->type));
746 return 1;
747 }
748 if (!(trcIdxMin < mtrc->traceNum)) {
749 biffAddf(PULLpullBiffKey, "%s: trcIdxMin %u not < traceNum %u", me,
750 trcIdxMin, mtrc->traceNum);
751 return 1;
752 }
753 if (trcNum) {
754 trcIdxMax = trcIdxMin + trcNum-1;
755 if (!(trcIdxMax < mtrc->traceNum)) {
756 biffAddf(PULLpullBiffKey, "%s: trcIdxMax %u = %u+%u-1 not < traceNum %u", me,
757 trcIdxMax, trcIdxMin, trcNum, mtrc->traceNum);
758 return 1;
759 }
760 } else {
761 trcIdxMax = mtrc->traceNum-1;
762 }
763 if (nmaskedpos || nmask) {
764 if (!( nmaskedpos && nmask )) {
765 biffAddf(PULLpullBiffKey, "%s: need either both or neither of nmaskedpos (%p)"
766 "and nmask (%p)", me, nmaskedpos, nmask);
767 return 1;
768 }
769 if (!( 2 == nmask->dim
770 && nrrdTypeUChar == nmask->type
771 && nplot->axis[0].size == nmask->axis[0].size
772 && nplot->axis[1].size == nmask->axis[1].size )) {
773 biffAddf(PULLpullBiffKey, "%s: got trace mask but wanted "
774 "2-D %s %u-by-%u (not %u-D %s %u-by-%u)\n", me,
775 airEnumStr(nrrdType, nrrdTypeUChar),
776 AIR_CAST(unsigned int, nplot->axis[0].size)((unsigned int)(nplot->axis[0].size)),
777 AIR_CAST(unsigned int, nplot->axis[1].size)((unsigned int)(nplot->axis[1].size)),
778 nmask->dim,
779 airEnumStr(nrrdType, nmask->type),
780 AIR_CAST(unsigned int, nmask->axis[0].size)((unsigned int)(nmask->axis[0].size)),
781 AIR_CAST(unsigned int, nmask->axis[1].size)((unsigned int)(nmask->axis[1].size)));
782 return 1;
783 }
784 mask = AIR_CAST(unsigned char *, nmask->data)((unsigned char *)(nmask->data));
785 } else {
786 mask = NULL((void*)0);
787 }
788
789 ssRange[0] = nplot->axis[0].min;
790 ssRange[1] = nplot->axis[0].max;
791 vRange[0] = nplot->axis[1].min;
792 vRange[1] = nplot->axis[1].max;
793 if (!( AIR_EXISTS(ssRange[0])(((int)(!((ssRange[0]) - (ssRange[0]))))) && AIR_EXISTS(ssRange[1])(((int)(!((ssRange[1]) - (ssRange[1]))))) &&
794 AIR_EXISTS(vRange[0])(((int)(!((vRange[0]) - (vRange[0]))))) && AIR_EXISTS(vRange[1])(((int)(!((vRange[1]) - (vRange[1]))))) )) {
795 biffAddf(PULLpullBiffKey, "%s: need both axis 0 (%g,%g) and 1 (%g,%g) min,max", me,
796 ssRange[0], ssRange[1], vRange[0], vRange[1]);
797 return 1;
798 }
799 if (1 != vRange[0]) {
800 biffAddf(PULLpullBiffKey, "%s: expected vRange[0] == 1 not %g", me, vRange[0]);
801 return 1;
802 }
803 mop = airMopNew();
804
805 mpos = NULL((void*)0);
806 if (nmaskedpos && nmask) {
807 nrrdEmpty(nmaskedpos);
808 mposArr = airArrayNew((void**)(&mpos), NULL((void*)0),
809 4*sizeof(double), 512 /* HEY */);
810 } else {
811 mposArr = NULL((void*)0);
812 }
813
814 nsmst = nrrdNew();
815 airMopAdd(mop, nsmst, (airMopper)nrrdNuke, airMopAlways);
816 plot = AIR_CAST(double *, nplot->data)((double *)(nplot->data));
817 filt = (nfilt
818 ? AIR_CAST(int *, nfilt->data)((int *)(nfilt->data))
819 : NULL((void*)0));
820 sizeS = AIR_CAST(unsigned int, nplot->axis[0].size)((unsigned int)(nplot->axis[0].size));
821 sizeT = AIR_CAST(unsigned int, nplot->axis[1].size)((unsigned int)(nplot->axis[1].size));
822 for (trcIdx=trcIdxMin; trcIdx<=trcIdxMax; trcIdx++) {
823 int pntIdx, pntNum;
824 const pullTrace *trc;
825 const double *vert, *stab, *strn;
826 unsigned int maskInCount;
827 double maskInPos[4];
828
829 if (filt && !filt[trcIdx]) {
830 continue;
831 }
832 trc = mtrc->trace[trcIdx];
833 if (pullTraceStopStub == trc->whyNowhere) {
834 continue;
835 }
836 if (strengthUse && !(trc->nstrn && trc->nstrn->data)) {
837 biffAddf(PULLpullBiffKey, "%s: requesting strength-based weighting, but don't have "
838 "strength info in trace %u", me, trcIdx);
839 airMopError(mop); return 1;
840 }
841 pntNum = AIR_CAST(int, trc->nvert->axis[1].size)((int)(trc->nvert->axis[1].size));
842 vert = AIR_CAST(double *, trc->nvert->data)((double *)(trc->nvert->data));
843 stab = AIR_CAST(double *, trc->nstab->data)((double *)(trc->nstab->data));
844 if (smooth > 0) {
845 double *smst;
846 if (nrrdCopy(nsmst, trc->nstab)) {
847 biffMovef(PULLpullBiffKey, NRRDnrrdBiffKey, "%s: trouble w/ trace %u", me, trcIdx);
848 airMopError(mop); return 1;
849 }
850 smst = AIR_CAST(double *, nsmst->data)((double *)(nsmst->data));
851 for (pntIdx=0; pntIdx<pntNum; pntIdx++) {
852 int ii, jj;
853 double ss, ww, ws;
854 ss = ws = 0;
855 for (jj=-smooth; jj<=smooth; jj++) {
856 ii = pntIdx+jj;
857 ii = AIR_CLAMP(0, ii, pntNum-1)((ii) < (0) ? (0) : ((ii) > (pntNum-1) ? (pntNum-1) : (
ii)));
858 ww = nrrdKernelBSpline3->eval1_d(AIR_AFFINE(-smooth-1, jj,( ((double)(2)-(-2))*((double)(jj)-(-smooth-1)) / ((double)(smooth
+1)-(-smooth-1)) + (-2))
859 smooth+1,-2,2)( ((double)(2)-(-2))*((double)(jj)-(-smooth-1)) / ((double)(smooth
+1)-(-smooth-1)) + (-2)), NULL((void*)0));
860 ws += ww;
861 ss += ww*stab[0 + 2*ii]*stab[1 + 2*ii];
862 }
863 smst[pntIdx] = ss/ws;
864 }
865 /* now redirect stab */
866 stab = smst;
867 }
868 strn = AIR_CAST(double *, (strengthUse && trc->nstrn((double *)((strengthUse && trc->nstrn ? trc->nstrn
->data : ((void*)0))))
869 ? trc->nstrn->data : NULL))((double *)((strengthUse && trc->nstrn ? trc->nstrn
->data : ((void*)0))));
870 /* would be nice to get some graphical indication of this */
871 fprintf(stderr__stderrp, "!%s: trace %u in [%u,%u]: %u points; stops = %s(%s) | %s(%s)\n",
872 me, trcIdx, trcIdxMin, trcIdxMax, pntNum,
873 airEnumStr(pullTraceStop, trc->whyStop[0]),
874 (pullTraceStopConstrFail == trc->whyStop[0]
875 ? airEnumStr(pullConstraintFail, trc->whyStopCFail[0])
876 : ""),
877 airEnumStr(pullTraceStop, trc->whyStop[1]),
878 (pullTraceStopConstrFail == trc->whyStop[1]
879 ? airEnumStr(pullConstraintFail, trc->whyStopCFail[1])
880 : ""));
881 /* */
882
883 if (mask) {
884 maskInCount = 0;
885 ELL_4V_SET(maskInPos, 0, 0, 0, 0)((maskInPos)[0] = (0), (maskInPos)[1] = (0), (maskInPos)[2] =
(0), (maskInPos)[3] = (0));
886 }
887 for (pntIdx=0; pntIdx<pntNum; pntIdx++) {
888 const double *pp;
889 double add, ww;
890 unsigned int sidx, vidx;
891 pp = vert + 4*pntIdx;
892 if (!(AIR_IN_OP(ssRange[0], pp[3], ssRange[1])((ssRange[0]) < (pp[3]) && (pp[3]) < (ssRange[1
])))) {
893 continue;
894 }
895 if (flatWght > 0) {
896 if (!pntIdx || pntIdx == pntNum-1) {
897 continue;
898 }
899 } else if (flatWght < 0) {
900 /* HACK: only show the seed point */
901 if (AIR_CAST(unsigned int, pntIdx)((unsigned int)(pntIdx)) != trc->seedIdx) {
902 continue;
903 }
904 }
905 sidx = airIndex(ssRange[0], pp[3], ssRange[1], sizeS);
906 vidx = airIndexClamp(1, stab[0 + 2*pntIdx]*stab[1 + 2*pntIdx], 0, sizeT);
907 add = strn ? strn[pntIdx] : 1;
908 if (flatWght > 0) {
909 double dx = ( ((vert + 4*(pntIdx+1))[3] - (vert + 4*(pntIdx-1))[3])
910 / (ssRange[1] - ssRange[0]) );
911 /*
912 double dx = ( ((vert + 4*(pntIdx+1))[3] - pp[3])
913 / (ssRange[1] - ssRange[0]) );
914 */
915 double dy = (stab[0 + 2*(pntIdx+1)]*stab[1 + 2*(pntIdx+1)]
916 - stab[0 + 2*(pntIdx-1)]*stab[1 + 2*(pntIdx-1)]);
917 ww = dx/sqrt(dx*dx + dy*dy);
918 } else {
919 ww = 1;
920 }
921 plot[sidx + sizeS*vidx] += AIR_MAX(0, ww*add)((0) > (ww*add) ? (0) : (ww*add));
922 if (mask && mask[sidx + sizeS*vidx] > 200) {
923 ELL_4V_ADD2(maskInPos, maskInPos, pp)((maskInPos)[0] = (maskInPos)[0] + (pp)[0], (maskInPos)[1] = (
maskInPos)[1] + (pp)[1], (maskInPos)[2] = (maskInPos)[2] + (pp
)[2], (maskInPos)[3] = (maskInPos)[3] + (pp)[3]);
924 maskInCount ++;
925 }
926 }
927 if (mask && maskInCount) {
928 unsigned int mpi = airArrayLenIncr(mposArr, 1);
929 ELL_4V_SCALE(mpos + 4*mpi, 1.0/maskInCount, maskInPos)((mpos + 4*mpi)[0] = (maskInPos)[0]*1.0/maskInCount, (mpos + 4
*mpi)[1] = (maskInPos)[1]*1.0/maskInCount, (mpos + 4*mpi)[2] =
(maskInPos)[2]*1.0/maskInCount, (mpos + 4*mpi)[3] = (maskInPos
)[3]*1.0/maskInCount);
930 }
931 }
932 if (mask && mposArr->len) {
933 if (nrrdMaybeAlloc_va(nmaskedpos, nrrdTypeDouble, 2,
934 AIR_CAST(size_t, 4)((size_t)(4)),
935 AIR_CAST(size_t, mposArr->len)((size_t)(mposArr->len)))) {
936 biffAddf(PULLpullBiffKey, "%s: couldn't allocate masked pos", me);
937 airMopError(mop); return 1;
938 }
939 memcpy(nmaskedpos->data, mposArr->data,__builtin___memcpy_chk (nmaskedpos->data, mposArr->data
, 4*(mposArr->len)*sizeof(double), __builtin_object_size (
nmaskedpos->data, 0))
940 4*(mposArr->len)*sizeof(double))__builtin___memcpy_chk (nmaskedpos->data, mposArr->data
, 4*(mposArr->len)*sizeof(double), __builtin_object_size (
nmaskedpos->data, 0));
941 }
942 airMopOkay(mop);
943 return 0;
944}
945
946static size_t
947nsizeof(const Nrrd *nrrd) {
948 return (nrrd
949 ? nrrdElementSize(nrrd)*nrrdElementNumber(nrrd)
950 : 0);
951}
952
953size_t
954pullTraceMultiSizeof(const pullTraceMulti *mtrc) {
955 size_t ret;
956 unsigned int ti;
957
958 if (!mtrc) {
959 return 0;
960 }
961 ret = 0;
962 for (ti=0; ti<mtrc->traceNum; ti++) {
963 ret += sizeof(pullTrace);
964 ret += nsizeof(mtrc->trace[ti]->nvert);
965 ret += nsizeof(mtrc->trace[ti]->nstrn);
966 ret += nsizeof(mtrc->trace[ti]->nstab);
967 }
968 ret += sizeof(pullTrace*)*(mtrc->traceArr->size);
969 return ret;
970}
971
972pullTraceMulti *
973pullTraceMultiNix(pullTraceMulti *mtrc) {
974
975 if (mtrc) {
976 airArrayNuke(mtrc->traceArr);
977 free(mtrc);
978 }
979 return NULL((void*)0);
980}
981
982
983#define PULL_MTRC_MAGIC"PULLMTRC0001" "PULLMTRC0001"
984#define DEMARK_STR"======" "======"
985
986static int
987tracewrite(FILE *file, const pullTrace *trc, unsigned int ti) {
988 static const char me[]="tracewrite";
989
990 /*
991 this was used to get ascii coordinates for a trace,
992 to help isolate (via emacs) one trace from a saved multi-trace
993 NrrdIoState *nio = nrrdIoStateNew();
994 nio->encoding = nrrdEncodingAscii;
995 */
996
997 fprintf(file, "%s %u\n", DEMARK_STR"======", ti);
998 ell_4v_print_d(file, trc->seedPos);
999#define WRITE(FF) \
1000 if (trc->FF && trc->FF->data) { \
1001 if (nrrdWrite(file, trc->FF, NULL((void*)0) /* nio */ )) { \
1002 biffMovef(PULLpullBiffKey, NRRDnrrdBiffKey, "%s: trouble with " #FF , me); \
1003 return 1; \
1004 } \
1005 } else { \
1006 fprintf(file, "NULL"); \
1007 } \
1008 fprintf(file, "\n")
1009
1010 fprintf(file, "nrrds: vert strn stab = %d %d %d\n",
1011 trc->nvert && trc->nvert->data,
1012 trc->nstrn && trc->nstrn->data,
1013 trc->nstab && trc->nstab->data);
1014 WRITE(nvert);
1015 WRITE(nstrn);
1016 WRITE(nstab);
1017 fprintf(file, "%u\n", trc->seedIdx);
1018 fprintf(file, "%s %s %s\n",
1019 airEnumStr(pullTraceStop, trc->whyStop[0]),
1020 airEnumStr(pullTraceStop, trc->whyStop[1]),
1021 airEnumStr(pullTraceStop, trc->whyNowhere));
1022#undef WRITE
1023 return 0;
1024}
1025
1026int
1027pullTraceMultiWrite(FILE *file, const pullTraceMulti *mtrc) {
1028 static const char me[]="pullTraceMultiWrite";
1029 unsigned int ti;
1030
1031 if (!(file && mtrc)) {
1032 biffAddf(PULLpullBiffKey, "%s: got NULL pointer", me);
1033 return 1;
1034 }
1035 fprintf(file, "%s\n", PULL_MTRC_MAGIC"PULLMTRC0001");
1036 fprintf(file, "%u traces\n", mtrc->traceNum);
1037
1038 for (ti=0; ti<mtrc->traceNum; ti++) {
1039 if (tracewrite(file, mtrc->trace[ti], ti)) {
1040 biffAddf(PULLpullBiffKey, "%s: trace %u/%u", me, ti, mtrc->traceNum);
1041 return 1;
1042 }
1043 }
1044 return 0;
1045}
1046
1047static int
1048traceread(pullTrace *trc, FILE *file, unsigned int _ti) {
1049 static const char me[]="traceread";
1050 char line[AIR_STRLEN_MED(256+1)], name[AIR_STRLEN_MED(256+1)];
1051 unsigned int ti, lineLen;
1052 int stops[3], hackhack, vertHN, strnHN, stabHN; /* HN == have nrrd */
1053
1054 sprintf(name, "separator")__builtin___sprintf_chk (name, 0, __builtin_object_size (name
, 2 > 1 ? 1 : 0), "separator");
1055 lineLen = airOneLine(file, line, AIR_STRLEN_MED(256+1));
1056 if (!lineLen) {
1057 biffAddf(PULLpullBiffKey, "%s: didn't get %s line", me, name);
1058 return 1;
1059 }
1060 if (1 != sscanf(line, DEMARK_STR"======" " %u", &ti)) {
1061 biffAddf(PULLpullBiffKey, "%s: \"%s\" doesn't look like %s line", me, line, name);
1062 return 1;
1063 }
1064 if (ti != _ti) {
1065 biffAddf(PULLpullBiffKey, "%s: read trace index %u but wanted %u", me, ti, _ti);
1066 return 1;
1067 }
1068 sprintf(name, "seed pos")__builtin___sprintf_chk (name, 0, __builtin_object_size (name
, 2 > 1 ? 1 : 0), "seed pos");
1069 lineLen = airOneLine(file, line, AIR_STRLEN_MED(256+1));
1070 if (!lineLen) {
1071 biffAddf(PULLpullBiffKey, "%s: didn't get %s line", me, name);
1072 return 1;
1073 }
1074 if (4 != sscanf(line, "%lg %lg %lg %lg", trc->seedPos + 0,
1075 trc->seedPos + 1, trc->seedPos + 2, trc->seedPos + 3)) {
1076 biffAddf(PULLpullBiffKey, "%s: couldn't parse %s line \"%s\" as 4 doubles",
1077 me, name, line);
1078 return 1;
1079 }
1080 sprintf(name, "have nrrds")__builtin___sprintf_chk (name, 0, __builtin_object_size (name
, 2 > 1 ? 1 : 0), "have nrrds");
1081 lineLen = airOneLine(file, line, AIR_STRLEN_MED(256+1));
1082 if (!lineLen) {
1083 biffAddf(PULLpullBiffKey, "%s: didn't get %s line", me, name);
1084 return 1;
1085 }
1086 if (3 != sscanf(line, "nrrds: vert strn stab = %d %d %d",
1087 &vertHN, &strnHN, &stabHN)) {
1088 biffAddf(PULLpullBiffKey, "%s: couldn't parse %s line", me, name);
1089 return 1;
1090 }
1091#define READ(FF)if (FFHN) { if (nrrdRead(trc->nFF, file, ((void*)0))) { biffMovef
(pullBiffKey, nrrdBiffKey, "%s: trouble with " "FF" , me); return
1; } fgetc(file); } else { airOneLine(file, line, (256+1)); } \
1092 if (FF##HN) { \
1093 if (nrrdRead(trc->n##FF, file, NULL((void*)0))) { \
1094 biffMovef(PULLpullBiffKey, NRRDnrrdBiffKey, "%s: trouble with " #FF , me); \
1095 return 1; \
1096 } \
1097 fgetc(file); \
1098 } else { \
1099 airOneLine(file, line, AIR_STRLEN_MED(256+1)); \
1100 }
1101 hackhack = nrrdStateVerboseIO; /* should be fixed in Teem v2 */
1102 nrrdStateVerboseIO = 0;
1103 READ(vert)if (vertHN) { if (nrrdRead(trc->nvert, file, ((void*)0))) {
biffMovef(pullBiffKey, nrrdBiffKey, "%s: trouble with " "vert"
, me); return 1; } fgetc(file); } else { airOneLine(file, line
, (256+1)); };
1104 READ(strn)if (strnHN) { if (nrrdRead(trc->nstrn, file, ((void*)0))) {
biffMovef(pullBiffKey, nrrdBiffKey, "%s: trouble with " "strn"
, me); return 1; } fgetc(file); } else { airOneLine(file, line
, (256+1)); };
1105 READ(stab)if (stabHN) { if (nrrdRead(trc->nstab, file, ((void*)0))) {
biffMovef(pullBiffKey, nrrdBiffKey, "%s: trouble with " "stab"
, me); return 1; } fgetc(file); } else { airOneLine(file, line
, (256+1)); };
1106 nrrdStateVerboseIO = hackhack;
1107
1108 sprintf(name, "seed idx")__builtin___sprintf_chk (name, 0, __builtin_object_size (name
, 2 > 1 ? 1 : 0), "seed idx");
1109 lineLen = airOneLine(file, line, AIR_STRLEN_MED(256+1));
1110 if (!lineLen) {
1111 biffAddf(PULLpullBiffKey, "%s: didn't get %s line", me, name);
1112 return 1;
1113 }
1114 if (1 != sscanf(line, "%u", &(trc->seedIdx))) {
1115 biffAddf(PULLpullBiffKey, "%s: didn't parse uint from %s line \"%s\"",
1116 me, name, line);
1117 return 1;
1118 }
1119 sprintf(name, "stops")__builtin___sprintf_chk (name, 0, __builtin_object_size (name
, 2 > 1 ? 1 : 0), "stops");
1120 lineLen = airOneLine(file, line, AIR_STRLEN_MED(256+1));
1121 if (!lineLen) {
1122 biffAddf(PULLpullBiffKey, "%s: didn't get %s line", me, name);
1123 return 1;
1124 }
1125 if (3 != airParseStrE(stops, line, " ", 3, pullTraceStop)) {
1126 biffAddf(PULLpullBiffKey, "%s: didn't see 3 %s on %s line \"%s\"", me,
1127 pullTraceStop->name, name, line);
1128 return 1;
1129 }
1130
1131 return 0;
1132}
1133int
1134pullTraceMultiRead(pullTraceMulti *mtrc, FILE *file) {
1135 static const char me[]="pullTraceMultiRead";
1136 char line[AIR_STRLEN_MED(256+1)], name[AIR_STRLEN_MED(256+1)];
1137 unsigned int lineLen, ti, tnum;
1138 pullTrace *trc;
1139
1140 if (!(mtrc && file)) {
1
Taking false branch
→
1141 biffAddf(PULLpullBiffKey, "%s: got NULL pointer", me);
1142 return 1;
1143 }
1144 airArrayLenSet(mtrc->traceArr, 0);
1145 sprintf(name, "magic")__builtin___sprintf_chk (name, 0, __builtin_object_size (name
, 2 > 1 ? 1 : 0), "magic");
1146 lineLen = airOneLine(file, line, AIR_STRLEN_MED(256+1));
1147 if (!lineLen) {
2
←
Assuming 'lineLen' is not equal to 0
→
3
←
Taking false branch
→
1148 biffAddf(PULLpullBiffKey, "%s: didn't get %s line", me, name);
1149 return 1;
1150 }
1151 if (strcmp(line, PULL_MTRC_MAGIC"PULLMTRC0001")) {
4
←
Taking false branch
→
1152 biffAddf(PULLpullBiffKey, "%s: %s line \"%s\" not expected \"%s\"",
1153 me, name, line, PULL_MTRC_MAGIC"PULLMTRC0001");
1154 return 1;
1155 }
1156
1157 sprintf(name, "# of traces")__builtin___sprintf_chk (name, 0, __builtin_object_size (name
, 2 > 1 ? 1 : 0), "# of traces");
1158 lineLen = airOneLine(file, line, AIR_STRLEN_MED(256+1));
1159 if (!lineLen) {
5
←
Assuming 'lineLen' is not equal to 0
→
6
←
Taking false branch
→
1160 biffAddf(PULLpullBiffKey, "%s: didn't get %s line", me, name);
1161 return 1;
1162 }
1163 if (1 != sscanf(line, "%u traces", &tnum)) {
7
←
Taking false branch
→
1164 biffAddf(PULLpullBiffKey, "%s: \"%s\" doesn't look like %s line", me, line, name);
1165 return 1;
1166 }
1167 for (ti=0; ti<tnum; ti++) {
8
←
Assuming 'ti' is < 'tnum'
→
9
←
Loop condition is true. Entering loop body
→
1168 int added;
1169 trc = pullTraceNew();
10
←
Calling 'pullTraceNew'
→
14
←
Returned allocated memory
→
1170 if (traceread(trc, file, ti)) {
15
←
Taking true branch
→
1171 biffAddf(PULLpullBiffKey, "%s: on trace %u/%u", me, ti, tnum);
16
←
Potential leak of memory pointed to by 'trc'
1172 return 1;
1173 }
1174 if (pullTraceMultiAdd(mtrc, trc, &added)) {
1175 biffAddf(PULLpullBiffKey, "%s: adding trace %u/%u", me, ti, tnum);
1176 return 1;
1177 }
1178 if (!added) {
1179 trc = pullTraceNix(trc);
1180 }
1181 }
1182
1183 return 0;
1184}