Bug Summary

File:src/pull/constraints.c
Location:line 398, column 10
Description:Value stored to '_tmpv' during its initialization 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
25#include "pull.h"
26#include "privatePull.h"
27
28#define DEBUG(0) (0)
29/* #define DEBUG (12 == point->idtag) */
30
31/*
32typedef struct {
33 double val, absval, grad[3];
34} stateIso;
35
36static int
37probeIso(pullTask *task, pullPoint *point, unsigned int iter, int cond,
38 double pos[3],
39 stateIso *state) {
40 static const char me[]="probeIso";
41
42 ELL_3V_COPY(point->pos, pos); / * NB: not touching point->pos[3] * /
43 _pullPointHistAdd(point, cond, AIR_NAN);
44 if (pullProbe(task, point)) {
45 biffAddf(PULL, "%s: on iter %u", me, iter);
46 return 1;
47 }
48 state->val = pullPointScalar(task->pctx, point,
49 pullInfoIsovalue,
50 state->grad, NULL);
51 state->absval = AIR_ABS(state->val);
52 return 0;
53}
54*/
55
56/* NOTE: this assumes variables "iter" (uint) and "me" (char*) */
57#define NORMALIZE_ERR(dir, grad, len)if (task->pctx->flag.zeroZ) grad[2]=0; ((len) = (sqrt((
(((grad)))[0]*(((grad)))[0] + (((grad)))[1]*(((grad)))[1] + (
((grad)))[2]*(((grad)))[2]))), (((dir))[0] = (1.0/(len))*((grad
))[0], ((dir))[1] = (1.0/(len))*((grad))[1], ((dir))[2] = (1.0
/(len))*((grad))[2])); if (!(len)) { biffAddf(pullBiffKey, "%s: got zero grad at (%g,%g,%g,%g) on iter %u\n"
, me, point->pos[0], point->pos[1], point->pos[2], point
->pos[3], iter); return 1; } \
58 if (task->pctx->flag.zeroZ) grad[2]=0; \
59 ELL_3V_NORM((dir), (grad), (len))((len) = (sqrt(((((grad)))[0]*(((grad)))[0] + (((grad)))[1]*(
((grad)))[1] + (((grad)))[2]*(((grad)))[2]))), (((dir))[0] = (
1.0/(len))*((grad))[0], ((dir))[1] = (1.0/(len))*((grad))[1],
((dir))[2] = (1.0/(len))*((grad))[2])); \
60 if (!(len)) { \
61 biffAddf(PULLpullBiffKey, "%s: got zero grad at (%g,%g,%g,%g) on iter %u\n", me,\
62 point->pos[0], point->pos[1], point->pos[2], \
63 point->pos[3], iter); \
64 return 1; \
65 }
66
67#define NORMALIZE(dir, grad, len) \
68 if (task->pctx->flag.zeroZ) grad[2]=0; \
69 ELL_3V_NORM((dir), (grad), (len))((len) = (sqrt(((((grad)))[0]*(((grad)))[0] + (((grad)))[1]*(
((grad)))[1] + (((grad)))[2]*(((grad)))[2]))), (((dir))[0] = (
1.0/(len))*((grad))[0], ((dir))[1] = (1.0/(len))*((grad))[1],
((dir))[2] = (1.0/(len))*((grad))[2])); \
70 if (!(len)) { \
71 ELL_3V_SET((dir), 0, 0, 0)(((dir))[0] = (0), ((dir))[1] = (0), ((dir))[2] = (0)) ; \
72 }
73
74
75/* ------------------------------- isosurface */
76
77
78
79#define PROBE(v, av, g) if (pullProbe(task, point)) { \
80 biffAddf(PULLpullBiffKey, "%s: on iter %u", me, iter); \
81 return 1; \
82 } \
83 (v) = pullPointScalar(task->pctx, point, \
84 pullInfoIsovalue, (g), NULL((void*)0)); \
85 (av) = AIR_ABS(v)((v) > 0.0f ? (v) : -(v))
86#define SAVE(state, aval, val, grad, pos) \
87 state[0] = aval; \
88 state[1] = val; \
89 ELL_3V_COPY(state + 1 + 1, grad)((state + 1 + 1)[0] = (grad)[0], (state + 1 + 1)[1] = (grad)[
1], (state + 1 + 1)[2] = (grad)[2]); \
90 ELL_3V_COPY(state + 1 + 1 + 3, pos)((state + 1 + 1 + 3)[0] = (pos)[0], (state + 1 + 1 + 3)[1] = (
pos)[1], (state + 1 + 1 + 3)[2] = (pos)[2])
91#define RESTORE(aval, val, grad, pos, state) \
92 aval = state[0]; \
93 val = state[1]; \
94 ELL_3V_COPY(grad, state + 1 + 1)((grad)[0] = (state + 1 + 1)[0], (grad)[1] = (state + 1 + 1)[
1], (grad)[2] = (state + 1 + 1)[2]); \
95 ELL_3V_COPY(pos, state + 1 + 1 + 3)((pos)[0] = (state + 1 + 1 + 3)[0], (pos)[1] = (state + 1 + 1
+ 3)[1], (pos)[2] = (state + 1 + 1 + 3)[2])
96
97static int
98constraintSatIso(pullTask *task, pullPoint *point,
99 double stepMax, double constrEps,
100 unsigned int iterMax,
101 /* output */
102 int *constrFailP) {
103 static const char me[]="constraintSatIso";
104 double
105 step, /* current step size */
106 val, aval, /* last and current function values */
107 hack, /* how to control re-tries in the context of a single
108 for-loop, instead of a nested do-while loop */
109 grad[4], dir[3], len, state[1 + 1 + 3 + 3];
110 unsigned int iter = 0; /* 0: initial probe, 1..iterMax: probes in loop */
111
112 PROBE(val, aval, grad);
113 SAVE(state, aval, val, grad, point->pos);
114 hack = 1;
115 for (iter=1; iter<=iterMax; iter++) {
116 /* consider? http://en.wikipedia.org/wiki/Halley%27s_method */
117 NORMALIZE(dir, grad, len);
118 if (!len) {
119 /* no gradient; back off */
120 hack *= task->pctx->sysParm.backStepScale;
121 RESTORE(aval, val, grad, point->pos, state);
122 continue;
123 }
124 step = -val/len; /* the newton-raphson step */
125 step = step > 0 ? AIR_MIN(stepMax, step)((stepMax) < (step) ? (stepMax) : (step)) : AIR_MAX(-stepMax, step)((-stepMax) > (step) ? (-stepMax) : (step));
126 ELL_3V_SCALE_INCR(point->pos, hack*step, dir)((point->pos)[0] += (hack*step)*(dir)[0], (point->pos)[
1] += (hack*step)*(dir)[1], (point->pos)[2] += (hack*step)
*(dir)[2]);
127 PROBE(val, aval, grad);
128 _pullPointHistAdd(point, pullCondConstraintSatA, val);
129 if (aval <= state[0]) { /* we're no further from the root */
130 if (AIR_ABS(step)((step) > 0.0f ? (step) : -(step)) < stepMax*constrEps) { /* HEY stepMax*constrEps vs constrEps */
131 /* we have converged! */
132 break;
133 }
134 SAVE(state, aval, val, grad, point->pos);
135 hack = 1;
136 } else { /* oops, try again, don't update dir or len, reset val */
137 hack *= task->pctx->sysParm.backStepScale;
138 RESTORE(aval, val, grad, point->pos, state);
139 }
140 }
141 if (iter > iterMax) {
142 *constrFailP = pullConstraintFailIterMaxed;
143 } else {
144 *constrFailP = AIR_FALSE0;
145 }
146 return 0;
147}
148
149#undef PROBE
150#undef SAVE
151#undef RESTORE
152
153
154
155/* ------------------------------- laplacian */
156
157
158
159#define PROBE(l) if (pullProbe(task, point)) { \
160 biffAddf(PULLpullBiffKey, "%s: on iter %u", me, iter); \
161 return 1; \
162 } \
163 (l) = pullPointScalar(task->pctx, point, \
164 pullInfoHeightLaplacian, NULL((void*)0), NULL((void*)0));
165#define PROBEG(l, g) \
166 PROBE(l); \
167 pullPointScalar(task->pctx, point, pullInfoHeight, (g), NULL((void*)0));
168
169static int
170constraintSatLapl(pullTask *task, pullPoint *point,
171 double stepMax, double constrEps,
172 unsigned int iterMax,
173 /* output */
174 int *constrFailP) {
175 static const char me[]="constraintSatLapl";
176 double
177 step, /* current step size */
178 valLast, val, /* last and current function values */
179 grad[4], dir[3], len,
180 posOld[3], posNew[3], tmpv[3];
181 double a=0, b=1, s, fa, fb, fs, tmp, diff;
182 int side = 0;
183 unsigned int iter = 0; /* 0: initial probe, 1..iterMax: probes in loop */
184
185 step = stepMax/2;
186 PROBEG(val, grad);
187 if (0 == val) {
188 /* already exactly at the zero, we're done. This actually happens! */
189 /* printf("!%s: a lapl == 0!\n", me); */
190 return 0;
191 }
192 valLast = val;
193 NORMALIZE(dir, grad, len);
194 /* first phase: follow normalized gradient until laplacian sign change */
195 for (iter=1; iter<=iterMax; iter++) {
196 double sgn;
197 ELL_3V_COPY(posOld, point->pos)((posOld)[0] = (point->pos)[0], (posOld)[1] = (point->pos
)[1], (posOld)[2] = (point->pos)[2]);
198 sgn = airSgn(val); /* lapl < 0 => downhill; lapl > 0 => uphill */
199 ELL_3V_SCALE_INCR(point->pos, sgn*step, dir)((point->pos)[0] += (sgn*step)*(dir)[0], (point->pos)[1
] += (sgn*step)*(dir)[1], (point->pos)[2] += (sgn*step)*(dir
)[2]);
200 PROBEG(val, grad);
201 _pullPointHistAdd(point, pullCondConstraintSatA, val);
202 if (val*valLast < 0) {
203 /* laplacian has changed sign; stop looking */
204 break;
205 }
206 valLast = val;
207 NORMALIZE(dir, grad, len);
208 }
209 if (iter > iterMax) {
210 *constrFailP = pullConstraintFailIterMaxed;
211 return 0;
212 }
213 /* second phase: find the zero-crossing, looking between
214 f(posOld)=valLast and f(posNew)=val */
215 ELL_3V_COPY(posNew, point->pos)((posNew)[0] = (point->pos)[0], (posNew)[1] = (point->pos
)[1], (posNew)[2] = (point->pos)[2]);
216 ELL_3V_SUB(tmpv, posNew, posOld)((tmpv)[0] = (posNew)[0] - (posOld)[0], (tmpv)[1] = (posNew)[
1] - (posOld)[1], (tmpv)[2] = (posNew)[2] - (posOld)[2]);
217 len = ELL_3V_LEN(tmpv)(sqrt((((tmpv))[0]*((tmpv))[0] + ((tmpv))[1]*((tmpv))[1] + ((
tmpv))[2]*((tmpv))[2])));
218 fa = valLast;
219 fb = val;
220 if (AIR_ABS(fa)((fa) > 0.0f ? (fa) : -(fa)) < AIR_ABS(fb)((fb) > 0.0f ? (fb) : -(fb))) {
221 ELL_SWAP2(a, b, tmp)((tmp)=(a),(a)=(b),(b)=(tmp)); ELL_SWAP2(fa, fb, tmp)((tmp)=(fa),(fa)=(fb),(fb)=(tmp));
222 }
223 for (iter=1; iter<=iterMax; iter++) {
224 s = AIR_AFFINE(fa, 0, fb, a, b)( ((double)(b)-(a))*((double)(0)-(fa)) / ((double)(fb)-(fa)) +
(a));
225 ELL_3V_LERP(point->pos, s, posOld, posNew)((point->pos)[0] = (((s))*(((posNew)[0]) - ((posOld)[0])) +
((posOld)[0])), (point->pos)[1] = (((s))*(((posNew)[1]) -
((posOld)[1])) + ((posOld)[1])), (point->pos)[2] = (((s))
*(((posNew)[2]) - ((posOld)[2])) + ((posOld)[2])));
226 PROBE(fs);
227 _pullPointHistAdd(point, pullCondConstraintSatB, 0.0);
228 if (0 == fs) {
229 /* exactly nailed the zero, we're done. This actually happens! */
230 printf("!%s: b lapl == 0!\n", me);
231 break;
232 }
233 /* "Illinois" false-position. Dumb, but it works. */
234 if (fs*fb > 0) { /* not between s and b */
235 b = s;
236 fb = fs;
237 if (+1 == side) {
238 fa /= 2;
239 }
240 side = +1;
241 } else { /* not between a and s */
242 a = s;
243 fa = fs;
244 if (-1 == side) {
245 fb /= 2;
246 }
247 side = -1;
248 }
249 diff = (b - a)*len;
250 if (AIR_ABS(diff)((diff) > 0.0f ? (diff) : -(diff)) < stepMax*constrEps) { /* HEY stepMax*constrEps vs constrEps */
251 /* converged! */
252 break;
253 }
254 }
255 if (iter > iterMax) {
256 *constrFailP = pullConstraintFailIterMaxed;
257 } else {
258 *constrFailP = AIR_FALSE0;
259 }
260 return 0;
261}
262#undef PROBE
263#undef PROBEG
264
265
266/* ------------------------------------------- height (line xor surf) */
267
268static int
269probeHeight(pullTask *task, pullPoint *point,
270 /* output */
271 double *heightP, double grad[3], double hess[9]) {
272 static const char me[]="probeHeight";
273
274 if (pullProbe(task, point)) {
275 biffAddf(PULLpullBiffKey, "%s: trouble", me);
276 return 1;
277 }
278 *heightP = pullPointScalar(task->pctx, point, pullInfoHeight, grad, hess);
279 return 0;
280}
281
282/*
283** creaseProj
284**
285** column-space of output posproj spans the directions along which
286** particle is allowed to move *downward* (in height) for constraint
287** satisfaction (according to tangent 1 or tangents 1&2); for seeking
288** minima where 2nd deriv is positive
289**
290** negproj is the same, but for points moving upwards (according to
291** negativetangent1 or negativetangent 1&2); for seeking
292** maxima where 2nd deriv is negative
293*/
294static void
295creaseProj(pullTask *task, pullPoint *point,
296 int tang1Use, int tang2Use,
297 int negtang1Use, int negtang2Use,
298 /* output */
299 double posproj[9], double negproj[9]) {
300 static const char me[]="creaseProj";
301 double pp[9];
302 double *tng;
303
304 ELL_3M_ZERO_SET(posproj)((((posproj)+0)[0] = (0), ((posproj)+0)[1] = (0), ((posproj)+
0)[2] = (0)), (((posproj)+3)[0] = (0), ((posproj)+3)[1] = (0)
, ((posproj)+3)[2] = (0)), (((posproj)+6)[0] = (0), ((posproj
)+6)[1] = (0), ((posproj)+6)[2] = (0)));
305 if (tang1Use) {
306 tng = point->info + task->pctx->infoIdx[pullInfoTangent1];
307 if (DEBUG(0)) {
308 fprintf(stderr__stderrp, "!%s: tng1 = %g %g %g\n", me, tng[0], tng[1], tng[2]);
309 }
310 ELL_3MV_OUTER(pp, tng, tng)((((pp)+0)[0] = ((tng)[0])*((tng))[0], ((pp)+0)[1] = ((tng)[0
])*((tng))[1], ((pp)+0)[2] = ((tng)[0])*((tng))[2]), (((pp)+3
)[0] = ((tng)[1])*((tng))[0], ((pp)+3)[1] = ((tng)[1])*((tng)
)[1], ((pp)+3)[2] = ((tng)[1])*((tng))[2]), (((pp)+6)[0] = ((
tng)[2])*((tng))[0], ((pp)+6)[1] = ((tng)[2])*((tng))[1], ((pp
)+6)[2] = ((tng)[2])*((tng))[2]));
311 ELL_3M_ADD2(posproj, posproj, pp)((posproj)[0] = (posproj)[0] + (pp)[0], (posproj)[1] = (posproj
)[1] + (pp)[1], (posproj)[2] = (posproj)[2] + (pp)[2], (posproj
)[3] = (posproj)[3] + (pp)[3], (posproj)[4] = (posproj)[4] + (
pp)[4], (posproj)[5] = (posproj)[5] + (pp)[5], (posproj)[6] =
(posproj)[6] + (pp)[6], (posproj)[7] = (posproj)[7] + (pp)[7
], (posproj)[8] = (posproj)[8] + (pp)[8]);
312 }
313 if (tang2Use) {
314 tng = point->info + task->pctx->infoIdx[pullInfoTangent2];
315 ELL_3MV_OUTER(pp, tng, tng)((((pp)+0)[0] = ((tng)[0])*((tng))[0], ((pp)+0)[1] = ((tng)[0
])*((tng))[1], ((pp)+0)[2] = ((tng)[0])*((tng))[2]), (((pp)+3
)[0] = ((tng)[1])*((tng))[0], ((pp)+3)[1] = ((tng)[1])*((tng)
)[1], ((pp)+3)[2] = ((tng)[1])*((tng))[2]), (((pp)+6)[0] = ((
tng)[2])*((tng))[0], ((pp)+6)[1] = ((tng)[2])*((tng))[1], ((pp
)+6)[2] = ((tng)[2])*((tng))[2]));
316 ELL_3M_ADD2(posproj, posproj, pp)((posproj)[0] = (posproj)[0] + (pp)[0], (posproj)[1] = (posproj
)[1] + (pp)[1], (posproj)[2] = (posproj)[2] + (pp)[2], (posproj
)[3] = (posproj)[3] + (pp)[3], (posproj)[4] = (posproj)[4] + (
pp)[4], (posproj)[5] = (posproj)[5] + (pp)[5], (posproj)[6] =
(posproj)[6] + (pp)[6], (posproj)[7] = (posproj)[7] + (pp)[7
], (posproj)[8] = (posproj)[8] + (pp)[8]);
317 }
318
319 ELL_3M_ZERO_SET(negproj)((((negproj)+0)[0] = (0), ((negproj)+0)[1] = (0), ((negproj)+
0)[2] = (0)), (((negproj)+3)[0] = (0), ((negproj)+3)[1] = (0)
, ((negproj)+3)[2] = (0)), (((negproj)+6)[0] = (0), ((negproj
)+6)[1] = (0), ((negproj)+6)[2] = (0)));
320 if (negtang1Use) {
321 tng = point->info + task->pctx->infoIdx[pullInfoNegativeTangent1];
322 ELL_3MV_OUTER(pp, tng, tng)((((pp)+0)[0] = ((tng)[0])*((tng))[0], ((pp)+0)[1] = ((tng)[0
])*((tng))[1], ((pp)+0)[2] = ((tng)[0])*((tng))[2]), (((pp)+3
)[0] = ((tng)[1])*((tng))[0], ((pp)+3)[1] = ((tng)[1])*((tng)
)[1], ((pp)+3)[2] = ((tng)[1])*((tng))[2]), (((pp)+6)[0] = ((
tng)[2])*((tng))[0], ((pp)+6)[1] = ((tng)[2])*((tng))[1], ((pp
)+6)[2] = ((tng)[2])*((tng))[2]));
323 ELL_3M_ADD2(negproj, negproj, pp)((negproj)[0] = (negproj)[0] + (pp)[0], (negproj)[1] = (negproj
)[1] + (pp)[1], (negproj)[2] = (negproj)[2] + (pp)[2], (negproj
)[3] = (negproj)[3] + (pp)[3], (negproj)[4] = (negproj)[4] + (
pp)[4], (negproj)[5] = (negproj)[5] + (pp)[5], (negproj)[6] =
(negproj)[6] + (pp)[6], (negproj)[7] = (negproj)[7] + (pp)[7
], (negproj)[8] = (negproj)[8] + (pp)[8]);
324 }
325 if (negtang2Use) {
326 tng = point->info + task->pctx->infoIdx[pullInfoNegativeTangent2];
327 ELL_3MV_OUTER(pp, tng, tng)((((pp)+0)[0] = ((tng)[0])*((tng))[0], ((pp)+0)[1] = ((tng)[0
])*((tng))[1], ((pp)+0)[2] = ((tng)[0])*((tng))[2]), (((pp)+3
)[0] = ((tng)[1])*((tng))[0], ((pp)+3)[1] = ((tng)[1])*((tng)
)[1], ((pp)+3)[2] = ((tng)[1])*((tng))[2]), (((pp)+6)[0] = ((
tng)[2])*((tng))[0], ((pp)+6)[1] = ((tng)[2])*((tng))[1], ((pp
)+6)[2] = ((tng)[2])*((tng))[2]));
328 ELL_3M_ADD2(negproj, negproj, pp)((negproj)[0] = (negproj)[0] + (pp)[0], (negproj)[1] = (negproj
)[1] + (pp)[1], (negproj)[2] = (negproj)[2] + (pp)[2], (negproj
)[3] = (negproj)[3] + (pp)[3], (negproj)[4] = (negproj)[4] + (
pp)[4], (negproj)[5] = (negproj)[5] + (pp)[5], (negproj)[6] =
(negproj)[6] + (pp)[6], (negproj)[7] = (negproj)[7] + (pp)[7
], (negproj)[8] = (negproj)[8] + (pp)[8]);
329 }
330
331 if (!tang1Use && !tang2Use && !negtang1Use && !negtang2Use) {
332 /* we must be after points, and so need freedom to go after them */
333 /* for now we do this via posproj not negproj; see haveNada below */
334 ELL_3M_IDENTITY_SET(posproj)((((posproj)+0)[0] = (1), ((posproj)+0)[1] = (0), ((posproj)+
0)[2] = (0)), (((posproj)+3)[0] = (0), ((posproj)+3)[1] = (1)
, ((posproj)+3)[2] = (0)), (((posproj)+6)[0] = (0), ((posproj
)+6)[1] = (0), ((posproj)+6)[2] = (1)));
335 }
336
337 return;
338}
339
340/* HEY: body of probeHeight could really be expanded in here */
341#define PROBE(height, grad, hess, posproj, negproj) \
342 if (probeHeight(task, point, \
343 &(height), (grad), (hess))) { \
344 biffAddf(PULLpullBiffKey, "%s: trouble on iter %u", me, iter); \
345 return 1; \
346 } \
347 creaseProj(task, point, tang1Use, tang2Use, \
348 negtang1Use, negtang2Use, posproj, negproj)
349#define SAVE(state, height, grad, hess, posproj, negproj, pos) \
350 state[0] = height; \
351 ELL_3V_COPY(state + 1, grad)((state + 1)[0] = (grad)[0], (state + 1)[1] = (grad)[1], (state
+ 1)[2] = (grad)[2]); \
352 ELL_3M_COPY(state + 1 + 3, hess)((((state + 1 + 3)+0)[0] = ((hess)+0)[0], ((state + 1 + 3)+0)
[1] = ((hess)+0)[1], ((state + 1 + 3)+0)[2] = ((hess)+0)[2]),
(((state + 1 + 3)+3)[0] = ((hess)+3)[0], ((state + 1 + 3)+3)
[1] = ((hess)+3)[1], ((state + 1 + 3)+3)[2] = ((hess)+3)[2]),
(((state + 1 + 3)+6)[0] = ((hess)+6)[0], ((state + 1 + 3)+6)
[1] = ((hess)+6)[1], ((state + 1 + 3)+6)[2] = ((hess)+6)[2])); \
353 ELL_3M_COPY(state + 1 + 3 + 9, posproj)((((state + 1 + 3 + 9)+0)[0] = ((posproj)+0)[0], ((state + 1 +
3 + 9)+0)[1] = ((posproj)+0)[1], ((state + 1 + 3 + 9)+0)[2] =
((posproj)+0)[2]), (((state + 1 + 3 + 9)+3)[0] = ((posproj)+
3)[0], ((state + 1 + 3 + 9)+3)[1] = ((posproj)+3)[1], ((state
+ 1 + 3 + 9)+3)[2] = ((posproj)+3)[2]), (((state + 1 + 3 + 9
)+6)[0] = ((posproj)+6)[0], ((state + 1 + 3 + 9)+6)[1] = ((posproj
)+6)[1], ((state + 1 + 3 + 9)+6)[2] = ((posproj)+6)[2])); \
354 ELL_3M_COPY(state + 1 + 3 + 9 + 9, negproj)((((state + 1 + 3 + 9 + 9)+0)[0] = ((negproj)+0)[0], ((state +
1 + 3 + 9 + 9)+0)[1] = ((negproj)+0)[1], ((state + 1 + 3 + 9
+ 9)+0)[2] = ((negproj)+0)[2]), (((state + 1 + 3 + 9 + 9)+3)
[0] = ((negproj)+3)[0], ((state + 1 + 3 + 9 + 9)+3)[1] = ((negproj
)+3)[1], ((state + 1 + 3 + 9 + 9)+3)[2] = ((negproj)+3)[2]), (
((state + 1 + 3 + 9 + 9)+6)[0] = ((negproj)+6)[0], ((state + 1
+ 3 + 9 + 9)+6)[1] = ((negproj)+6)[1], ((state + 1 + 3 + 9 +
9)+6)[2] = ((negproj)+6)[2])); \
355 ELL_3V_COPY(state + 1 + 3 + 9 + 9 + 9, pos)((state + 1 + 3 + 9 + 9 + 9)[0] = (pos)[0], (state + 1 + 3 + 9
+ 9 + 9)[1] = (pos)[1], (state + 1 + 3 + 9 + 9 + 9)[2] = (pos
)[2])
356#define RESTORE(height, grad, hess, posproj, negproj, pos, state) \
357 height = state[0]; \
358 ELL_3V_COPY(grad, state + 1)((grad)[0] = (state + 1)[0], (grad)[1] = (state + 1)[1], (grad
)[2] = (state + 1)[2]); \
359 ELL_3M_COPY(hess, state + 1 + 3)((((hess)+0)[0] = ((state + 1 + 3)+0)[0], ((hess)+0)[1] = ((state
+ 1 + 3)+0)[1], ((hess)+0)[2] = ((state + 1 + 3)+0)[2]), (((
hess)+3)[0] = ((state + 1 + 3)+3)[0], ((hess)+3)[1] = ((state
+ 1 + 3)+3)[1], ((hess)+3)[2] = ((state + 1 + 3)+3)[2]), (((
hess)+6)[0] = ((state + 1 + 3)+6)[0], ((hess)+6)[1] = ((state
+ 1 + 3)+6)[1], ((hess)+6)[2] = ((state + 1 + 3)+6)[2])); \
360 ELL_3M_COPY(posproj, state + 1 + 3 + 9)((((posproj)+0)[0] = ((state + 1 + 3 + 9)+0)[0], ((posproj)+0
)[1] = ((state + 1 + 3 + 9)+0)[1], ((posproj)+0)[2] = ((state
+ 1 + 3 + 9)+0)[2]), (((posproj)+3)[0] = ((state + 1 + 3 + 9
)+3)[0], ((posproj)+3)[1] = ((state + 1 + 3 + 9)+3)[1], ((posproj
)+3)[2] = ((state + 1 + 3 + 9)+3)[2]), (((posproj)+6)[0] = ((
state + 1 + 3 + 9)+6)[0], ((posproj)+6)[1] = ((state + 1 + 3 +
9)+6)[1], ((posproj)+6)[2] = ((state + 1 + 3 + 9)+6)[2])); \
361 ELL_3M_COPY(negproj, state + 1 + 3 + 9 + 9)((((negproj)+0)[0] = ((state + 1 + 3 + 9 + 9)+0)[0], ((negproj
)+0)[1] = ((state + 1 + 3 + 9 + 9)+0)[1], ((negproj)+0)[2] = (
(state + 1 + 3 + 9 + 9)+0)[2]), (((negproj)+3)[0] = ((state +
1 + 3 + 9 + 9)+3)[0], ((negproj)+3)[1] = ((state + 1 + 3 + 9
+ 9)+3)[1], ((negproj)+3)[2] = ((state + 1 + 3 + 9 + 9)+3)[2
]), (((negproj)+6)[0] = ((state + 1 + 3 + 9 + 9)+6)[0], ((negproj
)+6)[1] = ((state + 1 + 3 + 9 + 9)+6)[1], ((negproj)+6)[2] = (
(state + 1 + 3 + 9 + 9)+6)[2])); \
362 ELL_3V_COPY(pos, state + 1 + 3 + 9 + 9 + 9)((pos)[0] = (state + 1 + 3 + 9 + 9 + 9)[0], (pos)[1] = (state
+ 1 + 3 + 9 + 9 + 9)[1], (pos)[2] = (state + 1 + 3 + 9 + 9 +
9)[2])
363#define DNORM(d1, d2, pdir, plen, pgrad, grad, hess, posproj) \
364 ELL_3MV_MUL(pgrad, posproj, grad)((pgrad)[0] = (posproj)[0]*(grad)[0] + (posproj)[1]*(grad)[1]
+ (posproj)[2]*(grad)[2], (pgrad)[1] = (posproj)[3]*(grad)[0
] + (posproj)[4]*(grad)[1] + (posproj)[5]*(grad)[2], (pgrad)[
2] = (posproj)[6]*(grad)[0] + (posproj)[7]*(grad)[1] + (posproj
)[8]*(grad)[2]); \
365 if (task->pctx->flag.zeroZ) pgrad[2]=0; \
366 ELL_3V_NORM(pdir, pgrad, plen)(plen = (sqrt((((pgrad))[0]*((pgrad))[0] + ((pgrad))[1]*((pgrad
))[1] + ((pgrad))[2]*((pgrad))[2]))), ((pdir)[0] = (1.0/plen)
*(pgrad)[0], (pdir)[1] = (1.0/plen)*(pgrad)[1], (pdir)[2] = (
1.0/plen)*(pgrad)[2])); \
367 d1 = ELL_3V_DOT(grad, pdir)((grad)[0]*(pdir)[0] + (grad)[1]*(pdir)[1] + (grad)[2]*(pdir)
[2]); \
368 d2 = ELL_3MV_CONTR(hess, pdir)((hess)[0]*(pdir)[0]*(pdir)[0] + (hess)[1]*(pdir)[1]*(pdir)[0
] + (hess)[2]*(pdir)[2]*(pdir)[0] + (hess)[3]*(pdir)[0]*(pdir
)[1] + (hess)[4]*(pdir)[1]*(pdir)[1] + (hess)[5]*(pdir)[2]*(pdir
)[1] + (hess)[6]*(pdir)[0]*(pdir)[2] + (hess)[7]*(pdir)[1]*(pdir
)[2] + (hess)[8]*(pdir)[2]*(pdir)[2])
369#define PRINT(prefix)fprintf(__stderrp, "-------------- probe results %s (%u @ %g,%g,%g,%g):\n"
, prefix, point->idtag, point->pos[0], point->pos[1]
, point->pos[2], point->pos[3]); fprintf(__stderrp, "-- val = %g\n"
, val); fprintf(__stderrp, "-- grad = %g %g %g\n", grad[0], grad
[1], grad[2]); fprintf(__stderrp,"-- hess = %g %g %g; %g %g %g; %g %g %g\n"
, hess[0], hess[1], hess[2], hess[3], hess[4], hess[5], hess[
6], hess[7], hess[8]); fprintf(__stderrp, "-- posproj = %g %g %g; %g %g %g; %g %g %g\n"
, posproj[0], posproj[1], posproj[2], posproj[3], posproj[4],
posproj[5], posproj[6], posproj[7], posproj[8]); fprintf(__stderrp
, "-- negproj = %g %g %g; %g %g %g; %g %g %g\n", negproj[0]
, negproj[1], negproj[2], negproj[3], negproj[4], negproj[5],
negproj[6], negproj[7], negproj[8]) \
370 fprintf(stderr__stderrp, "-------------- probe results %s (%u @ %g,%g,%g,%g):\n", \
371 prefix, point->idtag, point->pos[0], point->pos[1], \
372 point->pos[2], point->pos[3]); \
373 fprintf(stderr__stderrp, "-- val = %g\n", val); \
374 fprintf(stderr__stderrp, "-- grad = %g %g %g\n", grad[0], grad[1], grad[2]); \
375 fprintf(stderr__stderrp,"-- hess = %g %g %g; %g %g %g; %g %g %g\n", \
376 hess[0], hess[1], hess[2], \
377 hess[3], hess[4], hess[5], \
378 hess[6], hess[7], hess[8]); \
379 fprintf(stderr__stderrp, "-- posproj = %g %g %g; %g %g %g; %g %g %g\n", \
380 posproj[0], posproj[1], posproj[2], \
381 posproj[3], posproj[4], posproj[5], \
382 posproj[6], posproj[7], posproj[8]); \
383 fprintf(stderr__stderrp, "-- negproj = %g %g %g; %g %g %g; %g %g %g\n", \
384 negproj[0], negproj[1], negproj[2], \
385 negproj[3], negproj[4], negproj[5], \
386 negproj[6], negproj[7], negproj[8])
387
388static int
389constraintSatHght(pullTask *task, pullPoint *point,
390 int tang1Use, int tang2Use,
391 int negtang1Use, int negtang2Use,
392 double stepMax, double constrEps, unsigned int iterMax,
393 int *constrFailP) {
394 static const char me[]="constraintSatHght";
395 double val, grad[3], hess[9], posproj[9], negproj[9],
396 state[1+3+9+9+9+3], hack, step,
397 d1, d2, pdir[3], plen, pgrad[3];
398 double _tmpv[3]={0,0,0};
Value stored to '_tmpv' during its initialization is never read
399 /* #endif */
400 int havePos, haveNeg, haveNada, zeroGmagOkay;
401 unsigned int iter = 0; /* 0: initial probe, 1..iterMax: probes in loop */
402 /* http://en.wikipedia.org/wiki/Newton%27s_method_in_optimization */
403
404 zeroGmagOkay = (1 < task->pctx->iter && 0 == task->pctx->constraintDim);
405 havePos = tang1Use || tang2Use;
406 haveNeg = negtang1Use || negtang2Use;
407 haveNada = !havePos && !haveNeg;
408 if (DEBUG(0)) {
409 double stpmin;
410 /* HEY: shouldn't stpmin also be used later in this function? */
411 stpmin = task->pctx->voxelSizeSpace*constrEps;
412 fprintf(stderr__stderrp, "!%s(%u): starting at %g %g %g %g\n", me, point->idtag,
413 point->pos[0], point->pos[1], point->pos[2], point->pos[3]);
414 fprintf(stderr__stderrp, "!%s: pt %d %d nt %d %d (nada %d) "
415 "stepMax %g, iterMax %u\n", me,
416 tang1Use, tang2Use, negtang1Use, negtang2Use, haveNada,
417 stepMax, iterMax);
418 fprintf(stderr__stderrp, "!%s: stpmin = %g = voxsize %g * parm.stepmin %g\n", me,
419 stpmin, task->pctx->voxelSizeSpace, constrEps);
420 }
421 PROBE(val, grad, hess, posproj, negproj);
422 _pullPointHistAdd(point, pullCondOld, val);
423 if (DEBUG(0)) {
424 PRINT("initial probe")fprintf(__stderrp, "-------------- probe results %s (%u @ %g,%g,%g,%g):\n"
, "initial probe", point->idtag, point->pos[0], point->
pos[1], point->pos[2], point->pos[3]); fprintf(__stderrp
, "-- val = %g\n", val); fprintf(__stderrp, "-- grad = %g %g %g\n"
, grad[0], grad[1], grad[2]); fprintf(__stderrp,"-- hess = %g %g %g; %g %g %g; %g %g %g\n"
, hess[0], hess[1], hess[2], hess[3], hess[4], hess[5], hess[
6], hess[7], hess[8]); fprintf(__stderrp, "-- posproj = %g %g %g; %g %g %g; %g %g %g\n"
, posproj[0], posproj[1], posproj[2], posproj[3], posproj[4],
posproj[5], posproj[6], posproj[7], posproj[8]); fprintf(__stderrp
, "-- negproj = %g %g %g; %g %g %g; %g %g %g\n", negproj[0]
, negproj[1], negproj[2], negproj[3], negproj[4], negproj[5],
negproj[6], negproj[7], negproj[8]);
425 }
426 SAVE(state, val, grad, hess, posproj, negproj, point->pos);
427 hack = 1;
428 for (iter=1; iter<=iterMax; iter++) {
429 if (DEBUG(0)) {
430 fprintf(stderr__stderrp, "!%s: =-============= begin iter %u\n", me, iter);
431 }
432 /* HEY: no opportunistic increase of hack? */
433 if (havePos || haveNada) {
434 DNORM(d1, d2, pdir, plen, pgrad, grad, hess, posproj);
435 if (!ELL_3M_FROB(hess)(sqrt((((hess)+0)[0]*((hess)+0)[0] + ((hess)+0)[1]*((hess)+0)
[1] + ((hess)+0)[2]*((hess)+0)[2]) + (((hess)+3)[0]*((hess)+3
)[0] + ((hess)+3)[1]*((hess)+3)[1] + ((hess)+3)[2]*((hess)+3)
[2]) + (((hess)+6)[0]*((hess)+6)[0] + ((hess)+6)[1]*((hess)+6
)[1] + ((hess)+6)[2]*((hess)+6)[2])))) {
436 *constrFailP = pullConstraintFailHessZeroA;
437 return 0;
438 }
439 if (!plen) {
440 if (zeroGmagOkay) {
441 /* getting to actual zero gradient is possible when looking for
442 point extrema (or saddles), and its not a problem, so as a
443 lousy hack we set step=0 and skip to the convergence test */
444 step = 0;
445 goto convtestA;
446 }
447 /* this use to be a biff error, which got to be annoying */
448 *constrFailP = pullConstraintFailProjGradZeroA;
449 return 0;
450 } else if (!AIR_EXISTS(plen)(((int)(!((plen) - (plen)))))) {
451 /* this use to be a biff error, which also got to be annoying */
452 *constrFailP = pullConstraintFailProjLenNonExist;
453 return 0;
454 }
455 step = (d2 > 0 ? -d1/d2 : -plen);
456 if (DEBUG(0)) {
457 fprintf(stderr__stderrp, "!%s: (+) iter %u step = (%g > 0 ? %g=-%g/%g : %g) --> %g\n",
458 me, iter, d2, -d1/d2, d1, d2, -plen, step);
459 }
460 step = step > 0 ? AIR_MIN(stepMax, step)((stepMax) < (step) ? (stepMax) : (step)) : AIR_MAX(-stepMax, step)((-stepMax) > (step) ? (-stepMax) : (step));
461 convtestA:
462 if (d2 > 0 && AIR_ABS(step)((step) > 0.0f ? (step) : -(step)) < constrEps) { /* HEY stepMax*constrEps vs constrEps */
463 /* we're converged because its concave up here
464 and we're close enough to the bottom */
465 if (DEBUG(0)) {
466 fprintf(stderr__stderrp, " |step| %g < %g*%g = %g ==> converged!\n",
467 AIR_ABS(step)((step) > 0.0f ? (step) : -(step)),
468 stepMax, constrEps,
469 stepMax*constrEps);
470 }
471 if (!haveNeg) {
472 break;
473 } else {
474 goto nextstep;
475 }
476 }
477 /* else we have to take a significant step */
478 if (DEBUG(0)) {
479 fprintf(stderr__stderrp, " -> step %g, |pdir| = %g\n",
480 step, ELL_3V_LEN(pdir)(sqrt((((pdir))[0]*((pdir))[0] + ((pdir))[1]*((pdir))[1] + ((
pdir))[2]*((pdir))[2]))));
481 ELL_3V_COPY(_tmpv, point->pos)((_tmpv)[0] = (point->pos)[0], (_tmpv)[1] = (point->pos
)[1], (_tmpv)[2] = (point->pos)[2]);
482 fprintf(stderr__stderrp, " -> pos (%g,%g,%g,%g) += "
483 "%g * %g * (%g,%g,%g)\n",
484 point->pos[0], point->pos[1], point->pos[2], point->pos[3],
485 hack, step, pdir[0], pdir[1], pdir[2]);
486 }
487 ELL_3V_SCALE_INCR(point->pos, hack*step, pdir)((point->pos)[0] += (hack*step)*(pdir)[0], (point->pos)
[1] += (hack*step)*(pdir)[1], (point->pos)[2] += (hack*step
)*(pdir)[2]);
488 if (!ELL_4V_EXISTS(point->pos)((((int)(!(((point->pos)[0]) - ((point->pos)[0]))))) &&
(((int)(!(((point->pos)[1]) - ((point->pos)[1]))))) &&
(((int)(!(((point->pos)[2]) - ((point->pos)[2]))))) &&
(((int)(!(((point->pos)[3]) - ((point->pos)[3]))))))) {
489 biffAddf(PULLpullBiffKey, "%s: pos proj iter %u: pnt %u bad pos (%g,%g,%g,%g); "
490 "hack %g, step %g",
491 me, iter, point->idtag, point->pos[0], point->pos[1],
492 point->pos[2], point->pos[3], hack, step);
493 return 1;
494 }
495 if (DEBUG(0)) {
496 ELL_3V_SUB(_tmpv, _tmpv, point->pos)((_tmpv)[0] = (_tmpv)[0] - (point->pos)[0], (_tmpv)[1] = (
_tmpv)[1] - (point->pos)[1], (_tmpv)[2] = (_tmpv)[2] - (point
->pos)[2]);
497 fprintf(stderr__stderrp, " -> moved to %g %g %g %g\n",
498 point->pos[0], point->pos[1], point->pos[2], point->pos[3]);
499 fprintf(stderr__stderrp, " (moved %g)\n", ELL_3V_LEN(_tmpv)(sqrt((((_tmpv))[0]*((_tmpv))[0] + ((_tmpv))[1]*((_tmpv))[1] +
((_tmpv))[2]*((_tmpv))[2]))));
500 }
501 PROBE(val, grad, hess, posproj, negproj);
502 _pullPointHistAdd(point, pullCondConstraintSatA, val);
503 if (DEBUG(0)) {
504 fprintf(stderr__stderrp, " (+) probed at (%g,%g,%g,%g)\n",
505 point->pos[0], point->pos[1], point->pos[2], point->pos[3]);
506 PRINT("after move")fprintf(__stderrp, "-------------- probe results %s (%u @ %g,%g,%g,%g):\n"
, "after move", point->idtag, point->pos[0], point->
pos[1], point->pos[2], point->pos[3]); fprintf(__stderrp
, "-- val = %g\n", val); fprintf(__stderrp, "-- grad = %g %g %g\n"
, grad[0], grad[1], grad[2]); fprintf(__stderrp,"-- hess = %g %g %g; %g %g %g; %g %g %g\n"
, hess[0], hess[1], hess[2], hess[3], hess[4], hess[5], hess[
6], hess[7], hess[8]); fprintf(__stderrp, "-- posproj = %g %g %g; %g %g %g; %g %g %g\n"
, posproj[0], posproj[1], posproj[2], posproj[3], posproj[4],
posproj[5], posproj[6], posproj[7], posproj[8]); fprintf(__stderrp
, "-- negproj = %g %g %g; %g %g %g; %g %g %g\n", negproj[0]
, negproj[1], negproj[2], negproj[3], negproj[4], negproj[5],
negproj[6], negproj[7], negproj[8]);
507 fprintf(stderr__stderrp, " val(%g,%g,%g,%g)=%g %s state[0]=%g\n",
508 point->pos[0], point->pos[1], point->pos[2], point->pos[3],
509 val, val <= state[0] ? "<=" : ">", state[0]);
510 }
511 if (val <= state[0]) {
512 /* we made progress */
513 if (DEBUG(0)) {
514 fprintf(stderr__stderrp, " (+) progress!\n");
515 }
516 SAVE(state, val, grad, hess, posproj, negproj, point->pos);
517 hack = 1;
518 } else {
519 /* oops, we went uphill instead of down; try again */
520 if (DEBUG(0)) {
521 fprintf(stderr__stderrp, " (+) val *increased* (up from %.17g by %.17g); "
522 "backing hack from %g to %g\n",
523 state[0], val - state[0],
524 hack, hack*task->pctx->sysParm.backStepScale);
525 }
526 hack *= task->pctx->sysParm.backStepScale;
527 RESTORE(val, grad, hess, posproj, negproj, point->pos, state);
528 if (DEBUG(0)) {
529 fprintf(stderr__stderrp, " restored to pos (%g,%g,%g,%g)\n",
530 point->pos[0], point->pos[1], point->pos[2], point->pos[3]);
531 }
532 }
533 }
534 nextstep:
535 if (haveNeg) {
536 /* HEY: copy and paste from above, minus fluff, and with A->B */
537 DNORM(d1, d2, pdir, plen, pgrad, grad, hess, negproj);
538 if (!ELL_3M_FROB(hess)(sqrt((((hess)+0)[0]*((hess)+0)[0] + ((hess)+0)[1]*((hess)+0)
[1] + ((hess)+0)[2]*((hess)+0)[2]) + (((hess)+3)[0]*((hess)+3
)[0] + ((hess)+3)[1]*((hess)+3)[1] + ((hess)+3)[2]*((hess)+3)
[2]) + (((hess)+6)[0]*((hess)+6)[0] + ((hess)+6)[1]*((hess)+6
)[1] + ((hess)+6)[2]*((hess)+6)[2])))) {
539 *constrFailP = pullConstraintFailHessZeroB;
540 return 0;
541 }
542 if (!plen) {
543 if (zeroGmagOkay) {
544 step = 0;
545 goto convtestB;
546 }
547 *constrFailP = pullConstraintFailProjGradZeroB;
548 return 0;
549 }
550 step = (d2 < 0 ? -d1/d2 : plen);
551 if (DEBUG(0)) {
552 fprintf(stderr__stderrp, "!%s: -+) iter %u step = (%g < 0 ? %g : %g) --> %g\n",
553 me, iter, d2, -d1/d2, plen, step);
554 }
555 step = step > 0 ? AIR_MIN(stepMax, step)((stepMax) < (step) ? (stepMax) : (step)) : AIR_MAX(-stepMax, step)((-stepMax) > (step) ? (-stepMax) : (step));
556 convtestB:
557 if (d2 < 0 && AIR_ABS(step)((step) > 0.0f ? (step) : -(step)) < constrEps) { /* HEY stepMax*constrEps vs constrEps */
558 /* we're converged because its concave down here
559 and we're close enough to the top */
560 if (DEBUG(0)) {
561 fprintf(stderr__stderrp, " |step| %g < %g*%g = %g ==> converged!\n",
562 AIR_ABS(step)((step) > 0.0f ? (step) : -(step)),
563 stepMax, constrEps,
564 stepMax*constrEps);
565 }
566 /* no further iteration needed; we're converged */
567 break;
568 }
569 /* else we have to take a significant step */
570 if (DEBUG(0)) {
571 fprintf(stderr__stderrp, " -> step %g, |pdir| = %g\n",
572 step, ELL_3V_LEN(pdir)(sqrt((((pdir))[0]*((pdir))[0] + ((pdir))[1]*((pdir))[1] + ((
pdir))[2]*((pdir))[2]))));
573 ELL_3V_COPY(_tmpv, point->pos)((_tmpv)[0] = (point->pos)[0], (_tmpv)[1] = (point->pos
)[1], (_tmpv)[2] = (point->pos)[2]);
574 fprintf(stderr__stderrp, " -> pos (%g,%g,%g,%g) += "
575 "%g * %g * (%g,%g,%g)\n",
576 point->pos[0], point->pos[1], point->pos[2], point->pos[3],
577 hack, step, pdir[0], pdir[1], pdir[2]);
578 }
579 ELL_3V_SCALE_INCR(point->pos, hack*step, pdir)((point->pos)[0] += (hack*step)*(pdir)[0], (point->pos)
[1] += (hack*step)*(pdir)[1], (point->pos)[2] += (hack*step
)*(pdir)[2]);
580 if (!ELL_4V_EXISTS(point->pos)((((int)(!(((point->pos)[0]) - ((point->pos)[0]))))) &&
(((int)(!(((point->pos)[1]) - ((point->pos)[1]))))) &&
(((int)(!(((point->pos)[2]) - ((point->pos)[2]))))) &&
(((int)(!(((point->pos)[3]) - ((point->pos)[3]))))))) {
581 biffAddf(PULLpullBiffKey, "%s: neg proj iter %u: pnt %u bad pos (%g,%g,%g,%g); "
582 "hack %g, step %g",
583 me, iter, point->idtag, point->pos[0], point->pos[1],
584 point->pos[2], point->pos[3], hack, step);
585 return 1;
586 }
587 if (DEBUG(0)) {
588 ELL_3V_SUB(_tmpv, _tmpv, point->pos)((_tmpv)[0] = (_tmpv)[0] - (point->pos)[0], (_tmpv)[1] = (
_tmpv)[1] - (point->pos)[1], (_tmpv)[2] = (_tmpv)[2] - (point
->pos)[2]);
589 fprintf(stderr__stderrp, " -> moved to %g %g %g %g\n",
590 point->pos[0], point->pos[1], point->pos[2], point->pos[3]);
591 fprintf(stderr__stderrp, " (moved %g)\n", ELL_3V_LEN(_tmpv)(sqrt((((_tmpv))[0]*((_tmpv))[0] + ((_tmpv))[1]*((_tmpv))[1] +
((_tmpv))[2]*((_tmpv))[2]))));
592 }
593 PROBE(val, grad, hess, posproj, negproj);
594 _pullPointHistAdd(point, pullCondConstraintSatB, val);
595 if (DEBUG(0)) {
596 fprintf(stderr__stderrp, " (-) probed at (%g,%g,%g,%g)\n",
597 point->pos[0], point->pos[1], point->pos[2], point->pos[3]);
598 PRINT("after move")fprintf(__stderrp, "-------------- probe results %s (%u @ %g,%g,%g,%g):\n"
, "after move", point->idtag, point->pos[0], point->
pos[1], point->pos[2], point->pos[3]); fprintf(__stderrp
, "-- val = %g\n", val); fprintf(__stderrp, "-- grad = %g %g %g\n"
, grad[0], grad[1], grad[2]); fprintf(__stderrp,"-- hess = %g %g %g; %g %g %g; %g %g %g\n"
, hess[0], hess[1], hess[2], hess[3], hess[4], hess[5], hess[
6], hess[7], hess[8]); fprintf(__stderrp, "-- posproj = %g %g %g; %g %g %g; %g %g %g\n"
, posproj[0], posproj[1], posproj[2], posproj[3], posproj[4],
posproj[5], posproj[6], posproj[7], posproj[8]); fprintf(__stderrp
, "-- negproj = %g %g %g; %g %g %g; %g %g %g\n", negproj[0]
, negproj[1], negproj[2], negproj[3], negproj[4], negproj[5],
negproj[6], negproj[7], negproj[8]);
599 fprintf(stderr__stderrp, " val(%g,%g,%g,%g)=%g %s state[0]=%g\n",
600 point->pos[0], point->pos[1], point->pos[2], point->pos[3],
601 val, val >= state[0] ? ">=" : "<", state[0]);
602 }
603 if (val >= state[0]) {
604 /* we made progress */
605 if (DEBUG(0)) {
606 fprintf(stderr__stderrp, " (-) progress!\n");
607 }
608 SAVE(state, val, grad, hess, posproj, negproj, point->pos);
609 hack = 1;
610 } else {
611 /* oops, we went downhill instead of up; try again */
612 if (DEBUG(0)) {
613 fprintf(stderr__stderrp, " (-) val *decreased* (down from %.17g by %.17g); "
614 "backing hack from %g to %g\n",
615 state[0], state[0] - val,
616 hack, hack*task->pctx->sysParm.backStepScale);
617 }
618 hack *= task->pctx->sysParm.backStepScale;
619 RESTORE(val, grad, hess, posproj, negproj, point->pos, state);
620 if (DEBUG(0)) {
621 fprintf(stderr__stderrp, " restored to pos (%g,%g,%g,%g)\n",
622 point->pos[0], point->pos[1], point->pos[2], point->pos[3]);
623 }
624 }
625 }
626 }
627 if (iter > iterMax) {
628 *constrFailP = pullConstraintFailIterMaxed;
629 } else {
630 *constrFailP = AIR_FALSE0;
631 }
632 _pullPointHistAdd(point, (*constrFailP
633 ? pullCondConstraintFail
634 : pullCondConstraintSuccess), AIR_NAN);
635 if (DEBUG(0)) {
636 printf("!%s: Finished %s: %s (%d)\n", me,
637 *constrFailP ? "with failure" : "OK",
638 *constrFailP ? airEnumStr(pullConstraintFail, *constrFailP) : "",
639 *constrFailP);
640 }
641 return 0;
642}
643#undef PROBE
644#undef DNORM
645#undef SAVE
646#undef RESTORE
647
648double
649_pullSigma(const pullContext *pctx, const double pos[4]) {
650 double ret=0;
651
652 if (pos && pos[3]) {
653 ret = (pctx->flag.scaleIsTau
654 ? gageSigOfTau(pos[3])airSigmaOfTau(pos[3])
655 : pos[3]);
656 }
657 return ret;
658}
659
660/* ------------------------------------------- */
661
662/* have to make sure that scale position point->pos[3]
663** is not modified anywhere in here: constraints are ONLY spatial
664**
665** This uses biff, but only for showstopper problems
666*/
667int
668_pullConstraintSatisfy(pullTask *task, pullPoint *point,
669 double travelMax,
670 /* output */
671 int *constrFailP) {
672 static const char me[]="_pullConstraintSatisfy";
673 double stepMax, constrEps;
674 unsigned int iterMax;
675 double pos3Orig[3], pos3Diff[3], travel;
676
677 ELL_3V_COPY(pos3Orig, point->pos)((pos3Orig)[0] = (point->pos)[0], (pos3Orig)[1] = (point->
pos)[1], (pos3Orig)[2] = (point->pos)[2]);
678 /* HEY the "10*" is based on isolated experiments; what's the principle? */
679 stepMax = 10*task->pctx->voxelSizeSpace;
680 iterMax = task->pctx->iterParm.constraintMax;
681 constrEps = task->pctx->voxelSizeSpace*task->pctx->sysParm.constraintStepMin;
682 /* * (1 + 0.2*_pullSigma(task->pctx, point->pos)); */
683 if (DEBUG(0)) {
684 fprintf(stderr__stderrp, "!%s(%d): hi ==== %g %g %g, stepMax = %g, iterMax = %u\n",
685 me, point->idtag, point->pos[0], point->pos[1], point->pos[2],
686 stepMax, iterMax);
687 }
688 task->pctx->count[pullCountConstraintSatisfy] += 1;
689 switch (task->pctx->constraint) {
690 case pullInfoHeightLaplacian: /* zero-crossing edges */
691 if (constraintSatLapl(task, point, stepMax/4, constrEps,
692 4*iterMax, constrFailP)) {
693 biffAddf(PULLpullBiffKey, "%s: trouble", me);
694 return 1;
695 }
696 break;
697 case pullInfoIsovalue:
698 if (constraintSatIso(task, point, stepMax, constrEps,
699 iterMax, constrFailP)) {
700 biffAddf(PULLpullBiffKey, "%s: trouble", me);
701 return 1;
702 }
703 break;
704 case pullInfoHeight:
705 if (constraintSatHght(task, point,
706 !!task->pctx->ispec[pullInfoTangent1],
707 !!task->pctx->ispec[pullInfoTangent2],
708 !!task->pctx->ispec[pullInfoNegativeTangent1],
709 !!task->pctx->ispec[pullInfoNegativeTangent2],
710 stepMax, constrEps, iterMax, constrFailP)) {
711 biffAddf(PULLpullBiffKey, "%s: trouble", me);
712 return 1;
713 }
714 break;
715 default:
716 fprintf(stderr__stderrp, "%s: constraint on %s (%d) unimplemented!!\n", me,
717 airEnumStr(pullInfo, task->pctx->constraint),
718 task->pctx->constraint);
719 }
720 ELL_3V_SUB(pos3Diff, pos3Orig, point->pos)((pos3Diff)[0] = (pos3Orig)[0] - (point->pos)[0], (pos3Diff
)[1] = (pos3Orig)[1] - (point->pos)[1], (pos3Diff)[2] = (pos3Orig
)[2] - (point->pos)[2]);
721 if (travelMax) {
722 travel = ELL_3V_LEN(pos3Diff)(sqrt((((pos3Diff))[0]*((pos3Diff))[0] + ((pos3Diff))[1]*((pos3Diff
))[1] + ((pos3Diff))[2]*((pos3Diff))[2])))/task->pctx->voxelSizeSpace;
723 if (travel > travelMax) {
724 *constrFailP = pullConstraintFailTravel;
725 if (DEBUG(0)) {
726 fprintf(stderr__stderrp, "!%s: travel %g > travelMax %g\n", me,
727 travel, travelMax);
728 }
729 }
730 }
731 if (DEBUG(0)) {
732 fprintf(stderr__stderrp, "!%s(%u) %s @ (%g,%g,%g) = (%g,%g,%g) + (%g,%g,%g)\n", me,
733 point->idtag,
734 (*constrFailP
735 ? airEnumStr(pullConstraintFail, *constrFailP)
736 : "#GOOD#"),
737 point->pos[0], point->pos[1], point->pos[2],
738 pos3Diff[0], pos3Diff[1], pos3Diff[2],
739 pos3Orig[0], pos3Orig[1], pos3Orig[2]);
740#if PULL_PHIST0
741 if (1) {
742 Nrrd *nhist;
743 FILE *fhist;
744 char fname[AIR_STRLEN_LARGE(512+1)];
745 nhist = nrrdNew();
746 sprintf(fname, "%04u-%04u-phist.nrrd", task->pctx->iter, point->idtag)__builtin___sprintf_chk (fname, 0, __builtin_object_size (fname
, 2 > 1 ? 1 : 0), "%04u-%04u-phist.nrrd", task->pctx->
iter, point->idtag);
747 if (pullPositionHistoryNrrdGet(nhist, task->pctx, point)) {
748 biffAddf(PULLpullBiffKey, "%s: trouble", me);
749 return 1;
750 }
751 if ((fhist = fopen(fname, "w"))) {
752 if (nrrdSave(fname, nhist, NULL((void*)0))) {
753 biffMovef(PULLpullBiffKey, NRRDnrrdBiffKey, "%s: trouble", me);
754 return 1;
755 }
756 fclose(fhist);
757 }
758 nrrdNuke(nhist);
759 }
760#endif
761 }
762 return 0;
763}
764
765#undef NORMALIZE
766
767/*
768** _pullConstraintTangent
769**
770** eigenvectors (with non-zero eigenvalues) of output proj are
771** (hopefully) approximate tangents to the manifold to which particles
772** are constrained. It is *not* the local tangent of the directions
773** along which particles are allowed to move during constraint
774** satisfaction (that is given by creaseProj for creases)
775**
776** this can assume that probe() has just been called
777*/
778void
779_pullConstraintTangent(pullTask *task, pullPoint *point,
780 /* output */
781 double proj[9]) {
782 double vec[4], nvec[3], outer[9], len, posproj[9], negproj[9];
783
784 ELL_3M_IDENTITY_SET(proj)((((proj)+0)[0] = (1), ((proj)+0)[1] = (0), ((proj)+0)[2] = (
0)), (((proj)+3)[0] = (0), ((proj)+3)[1] = (1), ((proj)+3)[2]
= (0)), (((proj)+6)[0] = (0), ((proj)+6)[1] = (0), ((proj)+6
)[2] = (1))); /* NOTE: we are starting with identity . . . */
785 switch (task->pctx->constraint) {
786 case pullInfoHeight:
787 creaseProj(task, point,
788 !!task->pctx->ispec[pullInfoTangent1],
789 !!task->pctx->ispec[pullInfoTangent2],
790 !!task->pctx->ispec[pullInfoNegativeTangent1],
791 !!task->pctx->ispec[pullInfoNegativeTangent2],
792 posproj, negproj);
793 /* .. and subracting out output from creaseProj */
794 ELL_3M_SUB(proj, proj, posproj)((proj)[0] = (proj)[0] - (posproj)[0], (proj)[1] = (proj)[1] -
(posproj)[1], (proj)[2] = (proj)[2] - (posproj)[2], (proj)[3
] = (proj)[3] - (posproj)[3], (proj)[4] = (proj)[4] - (posproj
)[4], (proj)[5] = (proj)[5] - (posproj)[5], (proj)[6] = (proj
)[6] - (posproj)[6], (proj)[7] = (proj)[7] - (posproj)[7], (proj
)[8] = (proj)[8] - (posproj)[8]);
795 ELL_3M_SUB(proj, proj, negproj)((proj)[0] = (proj)[0] - (negproj)[0], (proj)[1] = (proj)[1] -
(negproj)[1], (proj)[2] = (proj)[2] - (negproj)[2], (proj)[3
] = (proj)[3] - (negproj)[3], (proj)[4] = (proj)[4] - (negproj
)[4], (proj)[5] = (proj)[5] - (negproj)[5], (proj)[6] = (proj
)[6] - (negproj)[6], (proj)[7] = (proj)[7] - (negproj)[7], (proj
)[8] = (proj)[8] - (negproj)[8]);
796 break;
797 case pullInfoHeightLaplacian:
798 case pullInfoIsovalue:
799 if (pullInfoHeightLaplacian == task->pctx->constraint) {
800 /* using gradient of height as approx normal to laplacian 0-crossing */
801 pullPointScalar(task->pctx, point, pullInfoHeight, vec, NULL((void*)0));
802 } else {
803 pullPointScalar(task->pctx, point, pullInfoIsovalue, vec, NULL((void*)0));
804 }
805 ELL_3V_NORM(nvec, vec, len)(len = (sqrt((((vec))[0]*((vec))[0] + ((vec))[1]*((vec))[1] +
((vec))[2]*((vec))[2]))), ((nvec)[0] = (1.0/len)*(vec)[0], (
nvec)[1] = (1.0/len)*(vec)[1], (nvec)[2] = (1.0/len)*(vec)[2]
));
806 if (len) {
807 /* .. or and subracting out tensor product of normal with itself */
808 ELL_3MV_OUTER(outer, nvec, nvec)((((outer)+0)[0] = ((nvec)[0])*((nvec))[0], ((outer)+0)[1] = (
(nvec)[0])*((nvec))[1], ((outer)+0)[2] = ((nvec)[0])*((nvec))
[2]), (((outer)+3)[0] = ((nvec)[1])*((nvec))[0], ((outer)+3)[
1] = ((nvec)[1])*((nvec))[1], ((outer)+3)[2] = ((nvec)[1])*((
nvec))[2]), (((outer)+6)[0] = ((nvec)[2])*((nvec))[0], ((outer
)+6)[1] = ((nvec)[2])*((nvec))[1], ((outer)+6)[2] = ((nvec)[2
])*((nvec))[2]));
809 ELL_3M_SUB(proj, proj, outer)((proj)[0] = (proj)[0] - (outer)[0], (proj)[1] = (proj)[1] - (
outer)[1], (proj)[2] = (proj)[2] - (outer)[2], (proj)[3] = (proj
)[3] - (outer)[3], (proj)[4] = (proj)[4] - (outer)[4], (proj)
[5] = (proj)[5] - (outer)[5], (proj)[6] = (proj)[6] - (outer)
[6], (proj)[7] = (proj)[7] - (outer)[7], (proj)[8] = (proj)[8
] - (outer)[8]);
810 }
811 break;
812 }
813 return;
814}
815
816/*
817** returns the *dimension* (not codimension) of the constraint manifold:
818** 0 for points
819** 1 for lines
820** 2 for surfaces
821** This is nontrivial because of the different ways that constraints
822** can be expressed, combined with the possibility of pctx->flag.zeroZ
823**
824** a -1 return value represents a biff-able error
825*/
826int
827_pullConstraintDim(const pullContext *pctx) {
828 static const char me[]="_pullConstraintDim";
829 int ret, t1, t2, nt1, nt2;
830
831 switch (pctx->constraint) {
832 case pullInfoHeightLaplacian: /* zero-crossing edges */
833 ret = (pctx->flag.zeroZ ? 1 : 2);
834 break;
835 case pullInfoIsovalue:
836 ret = (pctx->flag.zeroZ ? 1 : 2);
837 break;
838 case pullInfoHeight:
839 t1 = !!pctx->ispec[pullInfoTangent1];
840 t2 = !!pctx->ispec[pullInfoTangent2];
841 nt1 = !!pctx->ispec[pullInfoNegativeTangent1];
842 nt2 = !!pctx->ispec[pullInfoNegativeTangent2];
843 switch (t1 + t2 + nt1 + nt2) {
844 case 0:
845 ret = 0;
846 break;
847 case 3:
848 if (pctx->flag.zeroZ) {
849 biffAddf(PULLpullBiffKey, "%s: can't have three of (%s,%s,%s,%s) tangents with "
850 "2-D data (pctx->flag.zeroZ)", me,
851 airEnumStr(pullInfo, pullInfoTangent1),
852 airEnumStr(pullInfo, pullInfoTangent2),
853 airEnumStr(pullInfo, pullInfoNegativeTangent1),
854 airEnumStr(pullInfo, pullInfoNegativeTangent2));
855 return -1;
856 } /* else we're in 3D; 3 constraints -> point features */
857 ret = 0;
858 break;
859 case 1:
860 ret = (pctx->flag.zeroZ ? 1 : 2);
861 break;
862 case 2:
863 ret = (pctx->flag.zeroZ ? 0 : 1);
864 break;
865 default:
866 biffAddf(PULLpullBiffKey, "%s: can't simultaneously use all tangents "
867 "(%s,%s,%s,%s) as this implies co-dimension of -1", me,
868 airEnumStr(pullInfo, pullInfoTangent1),
869 airEnumStr(pullInfo, pullInfoTangent2),
870 airEnumStr(pullInfo, pullInfoNegativeTangent1),
871 airEnumStr(pullInfo, pullInfoNegativeTangent2));
872 return -1;
873 }
874 break;
875 default:
876 biffAddf(PULLpullBiffKey, "%s: constraint on %s (%d) unimplemented", me,
877 airEnumStr(pullInfo, pctx->constraint), pctx->constraint);
878 return -1;
879 }
880 return ret;
881}
882