1 /*
2 Teem: Tools to process and visualize scientific data and images .
3 Copyright (C) 2012, 2011, 2010, 2009 University of Chicago
4 Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
5 Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
6
7 This library is free software; you can redistribute it and/or
8 modify it under the terms of the GNU Lesser General Public License
9 (LGPL) as published by the Free Software Foundation; either
10 version 2.1 of the License, or (at your option) any later version.
11 The terms of redistributing and/or modifying this software also
12 include exceptions to the LGPL that facilitate static linking.
13
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 along with this library; if not, write to Free Software Foundation, Inc.,
21 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22 */
23
24
25 #include "pull.h"
26 #include "privatePull.h"
27
28 #define PRAYING 0
29
30 /*
31 typedef struct {
32 double val, absval, grad[3];
33 } stateIso;
34
35 static int
36 probeIso(pullTask *task, pullPoint *point, unsigned int iter, int cond,
37 double pos[3],
38 stateIso *state) {
39 static const char me[]="probeIso";
40
41 ELL_3V_COPY(point->pos, pos); / * NB: not touching point->pos[3] * /
42 _pullPointHistAdd(point, cond);
43 if (pullProbe(task, point)) {
44 biffAddf(PULL, "%s: on iter %u", me, iter);
45 return 1;
46 }
47 state->val = pullPointScalar(task->pctx, point,
48 pullInfoIsovalue,
49 state->grad, NULL);
50 state->absval = AIR_ABS(state->val);
51 return 0;
52 }
53 */
54
55 /* NOTE: this assumes variables "iter" (uint) and "me" (char*) */
56 #define NORMALIZE_ERR(dir, grad, len) \
57 ELL_3V_NORM((dir), (grad), (len)); \
58 if (!(len)) { \
59 biffAddf(PULL, "%s: got zero grad at (%g,%g,%g,%g) on iter %u\n", me,\
60 point->pos[0], point->pos[1], point->pos[2], \
61 point->pos[3], iter); \
62 return 1; \
63 }
64
65 #define NORMALIZE(dir, grad, len) \
66 ELL_3V_NORM((dir), (grad), (len)); \
67 if (!(len)) { \
68 ELL_3V_SET((dir), 0, 0, 0) ; \
69 }
70
71
72 /* ------------------------------- isosurface */
73
74
75
76 #define PROBE(v, av, g) if (pullProbe(task, point)) { \
77 biffAddf(PULL, "%s: on iter %u", me, iter); \
78 return 1; \
79 } \
80 (v) = pullPointScalar(task->pctx, point, \
81 pullInfoIsovalue, (g), NULL); \
82 (av) = AIR_ABS(v)
83 #define SAVE(state, aval, val, grad, pos) \
84 state[0] = aval; \
85 state[1] = val; \
86 ELL_3V_COPY(state + 1 + 1, grad); \
87 ELL_3V_COPY(state + 1 + 1 + 3, pos)
88 #define RESTORE(aval, val, grad, pos, state) \
89 aval = state[0]; \
90 val = state[1]; \
91 ELL_3V_COPY(grad, state + 1 + 1); \
92 ELL_3V_COPY(pos, state + 1 + 1 + 3)
93
94 static int
constraintSatIso(pullTask * task,pullPoint * point,double stepMax,unsigned int iterMax,int * constrFailP)95 constraintSatIso(pullTask *task, pullPoint *point,
96 double stepMax, unsigned int iterMax,
97 /* output */
98 int *constrFailP) {
99 static const char me[]="constraintSatIso";
100 double
101 step, /* current step size */
102 val, aval, /* last and current function values */
103 hack, /* how to control re-tries in the context of a single
104 for-loop, instead of a nested do-while loop */
105 grad[4], dir[3], len, state[1 + 1 + 3 + 3];
106 unsigned int iter = 0; /* 0: initial probe, 1..iterMax: probes in loop */
107
108 PROBE(val, aval, grad);
109 SAVE(state, aval, val, grad, point->pos);
110 hack = 1;
111 for (iter=1; iter<=iterMax; iter++) {
112 /* consider? http://en.wikipedia.org/wiki/Halley%27s_method */
113 NORMALIZE(dir, grad, len);
114 if (!len) {
115 /* no gradient; back off */
116 hack *= task->pctx->sysParm.backStepScale;
117 RESTORE(aval, val, grad, point->pos, state);
118 continue;
119 }
120 step = -val/len; /* the newton-raphson step */
121 step = step > 0 ? AIR_MIN(stepMax, step) : AIR_MAX(-stepMax, step);
122 ELL_3V_SCALE_INCR(point->pos, hack*step, dir);
123 _pullPointHistAdd(point, pullCondConstraintSatA);
124 PROBE(val, aval, grad);
125 if (aval <= state[0]) { /* we're no further from the root */
126 if (AIR_ABS(step) < stepMax*task->pctx->sysParm.constraintStepMin) {
127 /* we have converged! */
128 break;
129 }
130 SAVE(state, aval, val, grad, point->pos);
131 hack = 1;
132 } else { /* oops, try again, don't update dir or len, reset val */
133 hack *= task->pctx->sysParm.backStepScale;
134 RESTORE(aval, val, grad, point->pos, state);
135 }
136 }
137 if (iter > iterMax) {
138 *constrFailP = pullConstraintFailIterMaxed;
139 } else {
140 *constrFailP = AIR_FALSE;
141 }
142 return 0;
143 }
144
145 #undef PROBE
146 #undef SAVE
147 #undef RESTORE
148
149
150
151 /* ------------------------------- laplacian */
152
153
154
155 #define PROBE(l) if (pullProbe(task, point)) { \
156 biffAddf(PULL, "%s: on iter %u", me, iter); \
157 return 1; \
158 } \
159 (l) = pullPointScalar(task->pctx, point, \
160 pullInfoHeightLaplacian, NULL, NULL);
161 #define PROBEG(l, g) \
162 PROBE(l); \
163 pullPointScalar(task->pctx, point, pullInfoHeight, (g), NULL);
164
165 static int
constraintSatLapl(pullTask * task,pullPoint * point,double stepMax,unsigned int iterMax,int * constrFailP)166 constraintSatLapl(pullTask *task, pullPoint *point,
167 double stepMax, unsigned int iterMax,
168 /* output */
169 int *constrFailP) {
170 static const char me[]="constraintSatLapl";
171 double
172 step, /* current step size */
173 valLast, val, /* last and current function values */
174 grad[4], dir[3], len,
175 posOld[3], posNew[3], tmpv[3];
176 double a=0, b=1, s, fa, fb, fs, tmp, diff;
177 int side = 0;
178 unsigned int iter = 0; /* 0: initial probe, 1..iterMax: probes in loop */
179
180 step = stepMax/2;
181 PROBEG(val, grad);
182 if (0 == val) {
183 /* already exactly at the zero, we're done. This actually happens! */
184 /* printf("!%s: a lapl == 0!\n", me); */
185 return 0;
186 }
187 valLast = val;
188 NORMALIZE(dir, grad, len);
189 /* first phase: follow normalized gradient until laplacian sign change */
190 for (iter=1; iter<=iterMax; iter++) {
191 double sgn;
192 ELL_3V_COPY(posOld, point->pos);
193 sgn = airSgn(val); /* lapl < 0 => downhill; lapl > 0 => uphill */
194 ELL_3V_SCALE_INCR(point->pos, sgn*step, dir);
195 _pullPointHistAdd(point, pullCondConstraintSatA);
196 PROBEG(val, grad);
197 if (val*valLast < 0) {
198 /* laplacian has changed sign; stop looking */
199 break;
200 }
201 valLast = val;
202 NORMALIZE(dir, grad, len);
203 }
204 if (iter > iterMax) {
205 *constrFailP = pullConstraintFailIterMaxed;
206 return 0;
207 }
208 /* second phase: find the zero-crossing, looking between
209 f(posOld)=valLast and f(posNew)=val */
210 ELL_3V_COPY(posNew, point->pos);
211 ELL_3V_SUB(tmpv, posNew, posOld);
212 len = ELL_3V_LEN(tmpv);
213 fa = valLast;
214 fb = val;
215 if (AIR_ABS(fa) < AIR_ABS(fb)) {
216 ELL_SWAP2(a, b, tmp); ELL_SWAP2(fa, fb, tmp);
217 }
218 for (iter=1; iter<=iterMax; iter++) {
219 s = AIR_AFFINE(fa, 0, fb, a, b);
220 ELL_3V_LERP(point->pos, s, posOld, posNew);
221 _pullPointHistAdd(point, pullCondConstraintSatB);
222 PROBE(fs);
223 if (0 == fs) {
224 /* exactly nailed the zero, we're done. This actually happens! */
225 printf("!%s: b lapl == 0!\n", me);
226 break;
227 }
228 /* "Illinois" false-position. Dumb, but it works. */
229 if (fs*fb > 0) { /* not between s and b */
230 b = s;
231 fb = fs;
232 if (+1 == side) {
233 fa /= 2;
234 }
235 side = +1;
236 } else { /* not between a and s */
237 a = s;
238 fa = fs;
239 if (-1 == side) {
240 fb /= 2;
241 }
242 side = -1;
243 }
244 diff = (b - a)*len;
245 if (AIR_ABS(diff) < stepMax*task->pctx->sysParm.constraintStepMin) {
246 /* converged! */
247 break;
248 }
249 }
250 if (iter > iterMax) {
251 *constrFailP = pullConstraintFailIterMaxed;
252 } else {
253 *constrFailP = AIR_FALSE;
254 }
255 return 0;
256 }
257 #undef PROBE
258 #undef PROBEG
259
260
261 /* ------------------------------------------- height (line xor surf) */
262
263 static int
probeHeight(pullTask * task,pullPoint * point,double * heightP,double grad[3],double hess[9])264 probeHeight(pullTask *task, pullPoint *point,
265 /* output */
266 double *heightP, double grad[3], double hess[9]) {
267 static const char me[]="probeHeight";
268
269 if (pullProbe(task, point)) {
270 biffAddf(PULL, "%s: trouble", me);
271 return 1;
272 }
273 *heightP = pullPointScalar(task->pctx, point, pullInfoHeight, grad, hess);
274 return 0;
275 }
276
277 /*
278 ** creaseProj
279 **
280 ** eigenvectors (with non-zero eigenvalues) of output posproj are
281 ** tangents to the directions along which particle is allowed to move
282 ** *downward* (in height) for constraint satisfaction (according to
283 ** tangent 1 or tangents 1&2)
284 **
285 ** negproj is the same, but for points moving upwards (according to
286 ** negativetangent1 or negativetangent 1&2)
287 */
288 static void
creaseProj(pullTask * task,pullPoint * point,int tang1Use,int tang2Use,int negtang1Use,int negtang2Use,double posproj[9],double negproj[9])289 creaseProj(pullTask *task, pullPoint *point,
290 int tang1Use, int tang2Use,
291 int negtang1Use, int negtang2Use,
292 /* output */
293 double posproj[9], double negproj[9]) {
294 #if PRAYING
295 static const char me[]="creaseProj";
296 #endif
297 double pp[9];
298 double *tng;
299
300 ELL_3M_ZERO_SET(posproj);
301 if (tang1Use) {
302 tng = point->info + task->pctx->infoIdx[pullInfoTangent1];
303 #if PRAYING
304 fprintf(stderr, "!%s: tng1 = %g %g %g\n", me, tng[0], tng[1], tng[2]);
305 #endif
306 ELL_3MV_OUTER(pp, tng, tng);
307 ELL_3M_ADD2(posproj, posproj, pp);
308 }
309 if (tang2Use) {
310 tng = point->info + task->pctx->infoIdx[pullInfoTangent2];
311 ELL_3MV_OUTER(pp, tng, tng);
312 ELL_3M_ADD2(posproj, posproj, pp);
313 }
314
315 ELL_3M_ZERO_SET(negproj);
316 if (negtang1Use) {
317 tng = point->info + task->pctx->infoIdx[pullInfoNegativeTangent1];
318 ELL_3MV_OUTER(pp, tng, tng);
319 ELL_3M_ADD2(negproj, negproj, pp);
320 }
321 if (negtang2Use) {
322 tng = point->info + task->pctx->infoIdx[pullInfoNegativeTangent2];
323 ELL_3MV_OUTER(pp, tng, tng);
324 ELL_3M_ADD2(negproj, negproj, pp);
325 }
326
327 if (!tang1Use && !tang2Use && !negtang1Use && !negtang2Use) {
328 /* we must be after points, and so need freedom to go after them */
329 /* for now we do this via posproj not negproj; see haveNada below */
330 ELL_3M_IDENTITY_SET(posproj);
331 }
332
333 return;
334 }
335
336 /* HEY: body of probeHeight could really be expanded in here */
337 #define PROBE(height, grad, hess, posproj, negproj) \
338 if (probeHeight(task, point, \
339 &(height), (grad), (hess))) { \
340 biffAddf(PULL, "%s: trouble on iter %u", me, iter); \
341 return 1; \
342 } \
343 creaseProj(task, point, tang1Use, tang2Use, \
344 negtang1Use, negtang2Use, posproj, negproj)
345 #define SAVE(state, height, grad, hess, posproj, negproj, pos) \
346 state[0] = height; \
347 ELL_3V_COPY(state + 1, grad); \
348 ELL_3M_COPY(state + 1 + 3, hess); \
349 ELL_3M_COPY(state + 1 + 3 + 9, posproj); \
350 ELL_3M_COPY(state + 1 + 3 + 9 + 9, negproj); \
351 ELL_3V_COPY(state + 1 + 3 + 9 + 9 + 9, pos)
352 #define RESTORE(height, grad, hess, posproj, negproj, pos, state) \
353 height = state[0]; \
354 ELL_3V_COPY(grad, state + 1); \
355 ELL_3M_COPY(hess, state + 1 + 3); \
356 ELL_3M_COPY(posproj, state + 1 + 3 + 9); \
357 ELL_3M_COPY(negproj, state + 1 + 3 + 9 + 9); \
358 ELL_3V_COPY(pos, state + 1 + 3 + 9 + 9 + 9)
359 #define POSNORM(d1, d2, pdir, plen, pgrad, grad, hess, posproj) \
360 ELL_3MV_MUL(pgrad, posproj, grad); \
361 ELL_3V_NORM(pdir, pgrad, plen); \
362 d1 = ELL_3V_DOT(grad, pdir); \
363 d2 = ELL_3MV_CONTR(hess, pdir)
364 #define NEGNORM(d1, d2, pdir, plen, pgrad, grad, hess, negproj) \
365 ELL_3MV_MUL(pgrad, negproj, grad); \
366 ELL_3V_NORM(pdir, pgrad, plen); \
367 d1 = -ELL_3V_DOT(grad, pdir); \
368 d2 = -ELL_3MV_CONTR(hess, pdir)
369 #define PRINT(prefix) \
370 fprintf(stderr, "-------------- probe results %s:\n-- val = %g\n", \
371 prefix, val); \
372 fprintf(stderr, "-- grad = %g %g %g\n", grad[0], grad[1], grad[2]); \
373 fprintf(stderr,"-- hess = %g %g %g; %g %g %g; %g %g %g\n", \
374 hess[0], hess[1], hess[2], \
375 hess[3], hess[4], hess[5], \
376 hess[6], hess[7], hess[8]); \
377 fprintf(stderr, "-- posproj = %g %g %g; %g %g %g; %g %g %g\n", \
378 posproj[0], posproj[1], posproj[2], \
379 posproj[3], posproj[4], posproj[5], \
380 posproj[6], posproj[7], posproj[8]); \
381 fprintf(stderr, "-- negproj = %g %g %g; %g %g %g; %g %g %g\n", \
382 negproj[0], negproj[1], negproj[2], \
383 negproj[3], negproj[4], negproj[5], \
384 negproj[6], negproj[7], negproj[8])
385
386 static int
constraintSatHght(pullTask * task,pullPoint * point,int tang1Use,int tang2Use,int negtang1Use,int negtang2Use,double stepMax,unsigned int iterMax,int * constrFailP)387 constraintSatHght(pullTask *task, pullPoint *point,
388 int tang1Use, int tang2Use,
389 int negtang1Use, int negtang2Use,
390 double stepMax, unsigned int iterMax,
391 int *constrFailP) {
392 static const char me[]="constraintSatHght";
393 double val, grad[3], hess[9], posproj[9], negproj[9],
394 state[1+3+9+9+9+3], hack, step,
395 d1, d2, pdir[3], plen, pgrad[3];
396 #if PRAYING
397 double _tmpv[3]={0,0,0};
398 #endif
399 int havePos, haveNeg, haveNada;
400 unsigned int iter = 0; /* 0: initial probe, 1..iterMax: probes in loop */
401 /* http://en.wikipedia.org/wiki/Newton%27s_method_in_optimization */
402
403 havePos = tang1Use || tang2Use;
404 haveNeg = negtang1Use || negtang2Use;
405 haveNada = !havePos && !haveNeg;
406 #if PRAYING
407 {
408 double stpmin;
409 /* HEY: shouldn't stpmin also be used later in this function? */
410 stpmin = task->pctx->voxelSizeSpace*task->pctx->sysParm.constraintStepMin;
411 fprintf(stderr, "!%s(%u): starting at %g %g %g %g\n", me, point->idtag,
412 point->pos[0], point->pos[1], point->pos[2], point->pos[3]);
413 fprintf(stderr, "!%s: pt %d %d nt %d %d (nada %d) "
414 "stepMax %g, iterMax %u\n", me,
415 tang1Use, tang2Use, negtang1Use, negtang2Use, haveNada,
416 stepMax, iterMax);
417 fprintf(stderr, "!%s: stpmin = %g = voxsize %g * parm.stepmin %g\n", me,
418 stpmin, task->pctx->voxelSizeSpace,
419 task->pctx->sysParm.constraintStepMin);
420 }
421 #endif
422 _pullPointHistAdd(point, pullCondOld);
423 PROBE(val, grad, hess, posproj, negproj);
424 #if PRAYING
425 PRINT("initial probe");
426 #endif
427 SAVE(state, val, grad, hess, posproj, negproj, point->pos);
428 hack = 1;
429 for (iter=1; iter<=iterMax; iter++) {
430 #if PRAYING
431 fprintf(stderr, "!%s: =============== begin iter %u\n", me, iter);
432 #endif
433 /* HEY: no opportunistic increase of hack? */
434 if (havePos || haveNada) {
435 POSNORM(d1, d2, pdir, plen, pgrad, grad, hess, posproj);
436 if (!plen) {
437 /* this use to be a biff error, which got to be annoying */
438 *constrFailP = pullConstraintFailProjGradZeroA;
439 return 0;
440 }
441 step = (d2 <= 0 ? -plen : -d1/d2);
442 #if PRAYING
443 fprintf(stderr, "!%s: (+) iter %u step = (%g <= 0 ? %g : %g) --> %g\n",
444 me, iter, d2, -plen, -d1/d2, step);
445 #endif
446 step = step > 0 ? AIR_MIN(stepMax, step) : AIR_MAX(-stepMax, step);
447 if (AIR_ABS(step) < stepMax*task->pctx->sysParm.constraintStepMin) {
448 /* no further iteration needed; we're converged */
449 #if PRAYING
450 fprintf(stderr, " |step| %g < %g*%g = %g ==> converged!\n",
451 AIR_ABS(step),
452 stepMax, task->pctx->sysParm.constraintStepMin,
453 stepMax*task->pctx->sysParm.constraintStepMin);
454 #endif
455 if (!haveNeg) {
456 break;
457 } else {
458 goto nextstep;
459 }
460 }
461 /* else we have to take a significant step */
462 #if PRAYING
463 fprintf(stderr, " -> step %g, |pdir| = %g\n",
464 step, ELL_3V_LEN(pdir));
465 ELL_3V_COPY(_tmpv, point->pos);
466 fprintf(stderr, " -> pos (%g,%g,%g,%g) += %g * %g * (%g,%g,%g)\n",
467 point->pos[0], point->pos[1], point->pos[2], point->pos[3],
468 hack, step, pdir[0], pdir[1], pdir[2]);
469 #endif
470 ELL_3V_SCALE_INCR(point->pos, hack*step, pdir);
471 #if PRAYING
472 ELL_3V_SUB(_tmpv, _tmpv, point->pos);
473 fprintf(stderr, " -> moved to %g %g %g %g\n",
474 point->pos[0], point->pos[1], point->pos[2], point->pos[3]);
475 fprintf(stderr, " (moved %g)\n", ELL_3V_LEN(_tmpv));
476 #endif
477 _pullPointHistAdd(point, pullCondConstraintSatA);
478 PROBE(val, grad, hess, posproj, negproj);
479 #if PRAYING
480 fprintf(stderr, " (+) probed at (%g,%g,%g,%g)\n",
481 point->pos[0], point->pos[1], point->pos[2], point->pos[3]);
482 PRINT("after move");
483 fprintf(stderr, " val(%g,%g,%g,%g)=%g %s state[0]=%g\n",
484 point->pos[0], point->pos[1], point->pos[2], point->pos[3],
485 val, val <= state[0] ? "<=" : ">", state[0]);
486 #endif
487 if (val <= state[0]) {
488 /* we made progress */
489 #if PRAYING
490 fprintf(stderr, " (+) progress!\n");
491 #endif
492 SAVE(state, val, grad, hess, posproj, negproj, point->pos);
493 hack = 1;
494 } else {
495 /* oops, we went uphill instead of down; try again */
496 #if PRAYING
497 fprintf(stderr, " val *increased*; backing hack from %g to %g\n",
498 hack, hack*task->pctx->sysParm.backStepScale);
499 #endif
500 hack *= task->pctx->sysParm.backStepScale;
501 RESTORE(val, grad, hess, posproj, negproj, point->pos, state);
502 #if PRAYING
503 fprintf(stderr, " restored to pos (%g,%g,%g,%g)\n",
504 point->pos[0], point->pos[1], point->pos[2], point->pos[3]);
505 #endif
506 }
507 }
508 nextstep:
509 if (haveNeg) {
510 /* HEY: copy and paste from above, minus fluff */
511 NEGNORM(d1, d2, pdir, plen, pgrad, grad, hess, negproj);
512 if (!plen && !haveNeg) {
513 /* this use to be a biff error, which got to be annoying */
514 *constrFailP = pullConstraintFailProjGradZeroA;
515 return 0;
516 }
517 step = (d2 <= 0 ? -plen : -d1/d2);
518 #if PRAYING
519 fprintf(stderr, "!%s: -+) iter %u step = (%g <= 0 ? %g : %g) --> %g\n",
520 me, iter, d2, -plen, -d1/d2, step);
521 #endif
522 step = step > 0 ? AIR_MIN(stepMax, step) : AIR_MAX(-stepMax, step);
523 if (AIR_ABS(step) < stepMax*task->pctx->sysParm.constraintStepMin) {
524 #if PRAYING
525 fprintf(stderr, " |step| %g < %g*%g = %g ==> converged!\n",
526 AIR_ABS(step),
527 stepMax, task->pctx->sysParm.constraintStepMin,
528 stepMax*task->pctx->sysParm.constraintStepMin);
529 #endif
530 /* no further iteration needed; we're converged */
531 break;
532 }
533 /* else we have to take a significant step */
534 #if PRAYING
535 fprintf(stderr, " -> step %g, |pdir| = %g\n",
536 step, ELL_3V_LEN(pdir));
537 ELL_3V_COPY(_tmpv, point->pos);
538 fprintf(stderr, " -> pos (%g,%g,%g,%g) += %g * %g * (%g,%g,%g)\n",
539 point->pos[0], point->pos[1], point->pos[2], point->pos[3],
540 hack, step, pdir[0], pdir[1], pdir[2]);
541 #endif
542 ELL_3V_SCALE_INCR(point->pos, hack*step, pdir);
543 #if PRAYING
544 ELL_3V_SUB(_tmpv, _tmpv, point->pos);
545 fprintf(stderr, " -> moved to %g %g %g %g\n",
546 point->pos[0], point->pos[1], point->pos[2], point->pos[3]);
547 fprintf(stderr, " (moved %g)\n", ELL_3V_LEN(_tmpv));
548 #endif
549 _pullPointHistAdd(point, pullCondConstraintSatA);
550 PROBE(val, grad, hess, posproj, negproj);
551 #if PRAYING
552 fprintf(stderr, " (-) probed at (%g,%g,%g,%g)\n",
553 point->pos[0], point->pos[1], point->pos[2], point->pos[3]);
554 PRINT("after move");
555 fprintf(stderr, " val(%g,%g,%g,%g)=%g %s state[0]=%g\n",
556 point->pos[0], point->pos[1], point->pos[2], point->pos[3],
557 val, val >= state[0] ? ">=" : "<", state[0]);
558 #endif
559 if (val >= state[0]) {
560 /* we made progress */
561 #if PRAYING
562 fprintf(stderr, " (-) progress!\n");
563 #endif
564 SAVE(state, val, grad, hess, posproj, negproj, point->pos);
565 hack = 1;
566 } else {
567 /* oops, we went uphill instead of down; try again */
568 #if PRAYING
569 fprintf(stderr, " val *increased*; backing hack from %g to %g\n",
570 hack, hack*task->pctx->sysParm.backStepScale);
571 #endif
572 hack *= task->pctx->sysParm.backStepScale;
573 RESTORE(val, grad, hess, posproj, negproj, point->pos, state);
574 #if PRAYING
575 fprintf(stderr, " restored to pos (%g,%g,%g,%g)\n",
576 point->pos[0], point->pos[1], point->pos[2], point->pos[3]);
577 #endif
578 }
579 }
580 }
581 if (iter > iterMax) {
582 *constrFailP = pullConstraintFailIterMaxed;
583 } else {
584 *constrFailP = AIR_FALSE;
585 }
586 /*
587 printf("!%s: %d %s\n", me, *constrFailP,
588 *constrFailP ? "FAILED!" : "ok");
589 */
590 return 0;
591 }
592 #undef PROBE
593 #undef POSNORM
594 #undef NEGNORM
595 #undef SAVE
596 #undef RESTORE
597
598 /* ------------------------------------------- */
599
600 /* HEY: have to make sure that scale position point->pos[3]
601 ** is not modified anywhere in here: constraints are ONLY spatial
602 **
603 ** This uses biff, but only for showstopper problems
604 */
605 int
_pullConstraintSatisfy(pullTask * task,pullPoint * point,double travelMax,int * constrFailP)606 _pullConstraintSatisfy(pullTask *task, pullPoint *point,
607 double travelMax,
608 /* output */
609 int *constrFailP) {
610 static const char me[]="_pullConstraintSatisfy";
611 double stepMax;
612 unsigned int iterMax;
613 double pos3Orig[3], pos3Diff[3], travel;
614
615 ELL_3V_COPY(pos3Orig, point->pos);
616 stepMax = task->pctx->voxelSizeSpace;
617 iterMax = task->pctx->iterParm.constraintMax;
618 /*
619 dlim = _pullDistLimit(task, point);
620 if (iterMax*stepMax > dlim) {
621 stepMax = dlim/iterMax;
622 }
623 */
624 /*
625 fprintf(stderr, "!%s(%d): hi ==== %g %g %g, stepMax = %g, iterMax = %u\n",
626 me, point->idtag, point->pos[0], point->pos[1], point->pos[2],
627 stepMax, iterMax);
628 */
629 task->pctx->count[pullCountConstraintSatisfy] += 1;
630 switch (task->pctx->constraint) {
631 case pullInfoHeightLaplacian: /* zero-crossing edges */
632 if (constraintSatLapl(task, point, stepMax/4, 4*iterMax, constrFailP)) {
633 biffAddf(PULL, "%s: trouble", me);
634 return 1;
635 }
636 break;
637 case pullInfoIsovalue:
638 if (constraintSatIso(task, point, stepMax, iterMax, constrFailP)) {
639 biffAddf(PULL, "%s: trouble", me);
640 return 1;
641 }
642 break;
643 case pullInfoHeight:
644 if (constraintSatHght(task, point,
645 !!task->pctx->ispec[pullInfoTangent1],
646 !!task->pctx->ispec[pullInfoTangent2],
647 !!task->pctx->ispec[pullInfoNegativeTangent1],
648 !!task->pctx->ispec[pullInfoNegativeTangent2],
649 stepMax, iterMax, constrFailP)) {
650 biffAddf(PULL, "%s: trouble", me);
651 return 1;
652 }
653 break;
654 default:
655 fprintf(stderr, "%s: constraint on %s (%d) unimplemented!!\n", me,
656 airEnumStr(pullInfo, task->pctx->constraint),
657 task->pctx->constraint);
658 }
659 ELL_3V_SUB(pos3Diff, pos3Orig, point->pos);
660 travel = ELL_3V_LEN(pos3Diff)/task->pctx->voxelSizeSpace;
661 if (travel > travelMax) {
662 *constrFailP = pullConstraintFailTravel;
663 }
664 /*
665 fprintf(stderr, "!%s(%u) %s @ (%g,%g,%g) = (%g,%g,%g) + (%g,%g,%g)\n", me,
666 point->idtag,
667 (*constrFailP
668 ? airEnumStr(pullConstraintFail, *constrFailP)
669 : "#GOOD#"),
670 point->pos[0], point->pos[1], point->pos[2],
671 pos3Diff[0], pos3Diff[1], pos3Diff[2],
672 pos3Orig[0], pos3Orig[1], pos3Orig[2]);
673 */
674 return 0;
675 }
676
677 #undef NORMALIZE
678
679 /*
680 ** _pullConstraintTangent
681 **
682 ** eigenvectors (with non-zero eigenvalues) of output proj are
683 ** (hopefully) approximate tangents to the manifold to which particles
684 ** are constrained. It is *not* the local tangent of the directions
685 ** along which particles are allowed to move during constraint
686 ** satisfaction (that is given by creaseProj for creases)
687 **
688 ** this can assume that probe() has just been called
689 */
690 void
_pullConstraintTangent(pullTask * task,pullPoint * point,double proj[9])691 _pullConstraintTangent(pullTask *task, pullPoint *point,
692 /* output */
693 double proj[9]) {
694 double vec[4], nvec[3], outer[9], len, posproj[9], negproj[9];
695
696 ELL_3M_IDENTITY_SET(proj); /* NOTE: we are starting with identity . . . */
697 switch (task->pctx->constraint) {
698 case pullInfoHeight:
699 creaseProj(task, point,
700 !!task->pctx->ispec[pullInfoTangent1],
701 !!task->pctx->ispec[pullInfoTangent2],
702 !!task->pctx->ispec[pullInfoNegativeTangent1],
703 !!task->pctx->ispec[pullInfoNegativeTangent2],
704 posproj, negproj);
705 /* .. and subracting out output from creaseProj */
706 ELL_3M_SUB(proj, proj, posproj);
707 ELL_3M_SUB(proj, proj, negproj);
708 break;
709 case pullInfoHeightLaplacian:
710 case pullInfoIsovalue:
711 if (pullInfoHeightLaplacian == task->pctx->constraint) {
712 /* using gradient of height as approx normal to laplacian 0-crossing */
713 pullPointScalar(task->pctx, point, pullInfoHeight, vec, NULL);
714 } else {
715 pullPointScalar(task->pctx, point, pullInfoIsovalue, vec, NULL);
716 }
717 ELL_3V_NORM(nvec, vec, len);
718 if (len) {
719 /* .. or and subracting out tensor product of normal with itself */
720 ELL_3MV_OUTER(outer, nvec, nvec);
721 ELL_3M_SUB(proj, proj, outer);
722 }
723 break;
724 }
725 return;
726 }
727
728 /*
729 ** returns the *dimension* (not codimension) of the constraint manifold:
730 ** 0 for points
731 ** 1 for lines
732 ** 2 for surfaces
733 **
734 ** a -1 return value represents a biff-able error
735 */
736 int
_pullConstraintDim(const pullContext * pctx)737 _pullConstraintDim(const pullContext *pctx) {
738 static const char me[]="_pullConstraintDim";
739 int ret, t1, t2, nt1, nt2;
740
741 switch (pctx->constraint) {
742 case pullInfoHeightLaplacian: /* zero-crossing edges */
743 ret = 2;
744 break;
745 case pullInfoIsovalue:
746 ret = 2;
747 break;
748 case pullInfoHeight:
749 t1 = !!pctx->ispec[pullInfoTangent1];
750 t2 = !!pctx->ispec[pullInfoTangent2];
751 nt1 = !!pctx->ispec[pullInfoNegativeTangent1];
752 nt2 = !!pctx->ispec[pullInfoNegativeTangent2];
753 switch (t1 + t2 + nt1 + nt2) {
754 case 0:
755 case 3:
756 ret = 0;
757 break;
758 case 1:
759 ret = 2;
760 break;
761 case 2:
762 ret = 1;
763 break;
764 default:
765 biffAddf(PULL, "%s: can't simultaneously use all tangents "
766 "(%s,%s,%s,%s) as this implies co-dimension of -1", me,
767 airEnumStr(pullInfo, pullInfoTangent1),
768 airEnumStr(pullInfo, pullInfoTangent2),
769 airEnumStr(pullInfo, pullInfoNegativeTangent1),
770 airEnumStr(pullInfo, pullInfoNegativeTangent2));
771 return -1;
772 }
773 break;
774 default:
775 biffAddf(PULL, "%s: constraint on %s (%d) unimplemented", me,
776 airEnumStr(pullInfo, pctx->constraint), pctx->constraint);
777 return -1;
778 }
779 return ret;
780 }
781
782