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) |
29 |
|
|
/* #define DEBUG (12 == point->idtag) */ |
30 |
|
|
|
31 |
|
|
/* |
32 |
|
|
typedef struct { |
33 |
|
|
double val, absval, grad[3]; |
34 |
|
|
} stateIso; |
35 |
|
|
|
36 |
|
|
static int |
37 |
|
|
probeIso(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) \ |
58 |
|
|
if (task->pctx->flag.zeroZ) grad[2]=0; \ |
59 |
|
|
ELL_3V_NORM((dir), (grad), (len)); \ |
60 |
|
|
if (!(len)) { \ |
61 |
|
|
biffAddf(PULL, "%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)); \ |
70 |
|
|
if (!(len)) { \ |
71 |
|
|
ELL_3V_SET((dir), 0, 0, 0) ; \ |
72 |
|
|
} |
73 |
|
|
|
74 |
|
|
|
75 |
|
|
/* ------------------------------- isosurface */ |
76 |
|
|
|
77 |
|
|
|
78 |
|
|
|
79 |
|
|
#define PROBE(v, av, g) if (pullProbe(task, point)) { \ |
80 |
|
|
biffAddf(PULL, "%s: on iter %u", me, iter); \ |
81 |
|
|
return 1; \ |
82 |
|
|
} \ |
83 |
|
|
(v) = pullPointScalar(task->pctx, point, \ |
84 |
|
|
pullInfoIsovalue, (g), NULL); \ |
85 |
|
|
(av) = AIR_ABS(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); \ |
90 |
|
|
ELL_3V_COPY(state + 1 + 1 + 3, pos) |
91 |
|
|
#define RESTORE(aval, val, grad, pos, state) \ |
92 |
|
|
aval = state[0]; \ |
93 |
|
|
val = state[1]; \ |
94 |
|
|
ELL_3V_COPY(grad, state + 1 + 1); \ |
95 |
|
|
ELL_3V_COPY(pos, state + 1 + 1 + 3) |
96 |
|
|
|
97 |
|
|
static int |
98 |
|
|
constraintSatIso(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) : AIR_MAX(-stepMax, step); |
126 |
|
|
ELL_3V_SCALE_INCR(point->pos, hack*step, dir); |
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) < 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_FALSE; |
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(PULL, "%s: on iter %u", me, iter); \ |
161 |
|
|
return 1; \ |
162 |
|
|
} \ |
163 |
|
|
(l) = pullPointScalar(task->pctx, point, \ |
164 |
|
|
pullInfoHeightLaplacian, NULL, NULL); |
165 |
|
|
#define PROBEG(l, g) \ |
166 |
|
|
PROBE(l); \ |
167 |
|
|
pullPointScalar(task->pctx, point, pullInfoHeight, (g), NULL); |
168 |
|
|
|
169 |
|
|
static int |
170 |
|
|
constraintSatLapl(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); |
198 |
|
|
sgn = airSgn(val); /* lapl < 0 => downhill; lapl > 0 => uphill */ |
199 |
|
|
ELL_3V_SCALE_INCR(point->pos, sgn*step, dir); |
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); |
216 |
|
|
ELL_3V_SUB(tmpv, posNew, posOld); |
217 |
|
|
len = ELL_3V_LEN(tmpv); |
218 |
|
|
fa = valLast; |
219 |
|
|
fb = val; |
220 |
|
|
if (AIR_ABS(fa) < AIR_ABS(fb)) { |
221 |
|
|
ELL_SWAP2(a, b, tmp); ELL_SWAP2(fa, fb, tmp); |
222 |
|
|
} |
223 |
|
|
for (iter=1; iter<=iterMax; iter++) { |
224 |
|
|
s = AIR_AFFINE(fa, 0, fb, a, b); |
225 |
|
|
ELL_3V_LERP(point->pos, s, posOld, posNew); |
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) < 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_FALSE; |
259 |
|
|
} |
260 |
|
|
return 0; |
261 |
|
|
} |
262 |
|
|
#undef PROBE |
263 |
|
|
#undef PROBEG |
264 |
|
|
|
265 |
|
|
|
266 |
|
|
/* ------------------------------------------- height (line xor surf) */ |
267 |
|
|
|
268 |
|
|
static int |
269 |
|
|
probeHeight(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(PULL, "%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 |
|
|
*/ |
294 |
|
|
static void |
295 |
|
|
creaseProj(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); |
305 |
|
|
if (tang1Use) { |
306 |
|
|
tng = point->info + task->pctx->infoIdx[pullInfoTangent1]; |
307 |
|
|
if (DEBUG) { |
308 |
|
|
fprintf(stderr, "!%s: tng1 = %g %g %g\n", me, tng[0], tng[1], tng[2]); |
309 |
|
|
} |
310 |
|
|
ELL_3MV_OUTER(pp, tng, tng); |
311 |
|
|
ELL_3M_ADD2(posproj, posproj, pp); |
312 |
|
|
} |
313 |
|
|
if (tang2Use) { |
314 |
|
|
tng = point->info + task->pctx->infoIdx[pullInfoTangent2]; |
315 |
|
|
ELL_3MV_OUTER(pp, tng, tng); |
316 |
|
|
ELL_3M_ADD2(posproj, posproj, pp); |
317 |
|
|
} |
318 |
|
|
|
319 |
|
|
ELL_3M_ZERO_SET(negproj); |
320 |
|
|
if (negtang1Use) { |
321 |
|
|
tng = point->info + task->pctx->infoIdx[pullInfoNegativeTangent1]; |
322 |
|
|
ELL_3MV_OUTER(pp, tng, tng); |
323 |
|
|
ELL_3M_ADD2(negproj, negproj, pp); |
324 |
|
|
} |
325 |
|
|
if (negtang2Use) { |
326 |
|
|
tng = point->info + task->pctx->infoIdx[pullInfoNegativeTangent2]; |
327 |
|
|
ELL_3MV_OUTER(pp, tng, tng); |
328 |
|
|
ELL_3M_ADD2(negproj, negproj, pp); |
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); |
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(PULL, "%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); \ |
352 |
|
|
ELL_3M_COPY(state + 1 + 3, hess); \ |
353 |
|
|
ELL_3M_COPY(state + 1 + 3 + 9, posproj); \ |
354 |
|
|
ELL_3M_COPY(state + 1 + 3 + 9 + 9, negproj); \ |
355 |
|
|
ELL_3V_COPY(state + 1 + 3 + 9 + 9 + 9, pos) |
356 |
|
|
#define RESTORE(height, grad, hess, posproj, negproj, pos, state) \ |
357 |
|
|
height = state[0]; \ |
358 |
|
|
ELL_3V_COPY(grad, state + 1); \ |
359 |
|
|
ELL_3M_COPY(hess, state + 1 + 3); \ |
360 |
|
|
ELL_3M_COPY(posproj, state + 1 + 3 + 9); \ |
361 |
|
|
ELL_3M_COPY(negproj, state + 1 + 3 + 9 + 9); \ |
362 |
|
|
ELL_3V_COPY(pos, state + 1 + 3 + 9 + 9 + 9) |
363 |
|
|
#define DNORM(d1, d2, pdir, plen, pgrad, grad, hess, posproj) \ |
364 |
|
|
ELL_3MV_MUL(pgrad, posproj, grad); \ |
365 |
|
|
if (task->pctx->flag.zeroZ) pgrad[2]=0; \ |
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 (%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, "-- val = %g\n", val); \ |
374 |
|
|
fprintf(stderr, "-- grad = %g %g %g\n", grad[0], grad[1], grad[2]); \ |
375 |
|
|
fprintf(stderr,"-- 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, "-- 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, "-- 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 |
|
|
|
388 |
|
|
static int |
389 |
|
|
constraintSatHght(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}; |
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) { |
409 |
|
|
double stpmin; |
410 |
|
|
/* HEY: shouldn't stpmin also be used later in this function? */ |
411 |
|
|
stpmin = task->pctx->voxelSizeSpace*constrEps; |
412 |
|
|
fprintf(stderr, "!%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, "!%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, "!%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) { |
424 |
|
|
PRINT("initial probe"); |
425 |
|
|
} |
426 |
|
|
SAVE(state, val, grad, hess, posproj, negproj, point->pos); |
427 |
|
|
hack = 1; |
428 |
|
|
for (iter=1; iter<=iterMax; iter++) { |
429 |
|
|
if (DEBUG) { |
430 |
|
|
fprintf(stderr, "!%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)) { |
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)) { |
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) { |
457 |
|
|
fprintf(stderr, "!%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) : AIR_MAX(-stepMax, step); |
461 |
|
|
convtestA: |
462 |
|
|
if (d2 > 0 && AIR_ABS(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) { |
466 |
|
|
fprintf(stderr, " |step| %g < %g*%g = %g ==> converged!\n", |
467 |
|
|
AIR_ABS(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) { |
479 |
|
|
fprintf(stderr, " -> step %g, |pdir| = %g\n", |
480 |
|
|
step, ELL_3V_LEN(pdir)); |
481 |
|
|
ELL_3V_COPY(_tmpv, point->pos); |
482 |
|
|
fprintf(stderr, " -> 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); |
488 |
|
|
if (!ELL_4V_EXISTS(point->pos)) { |
489 |
|
|
biffAddf(PULL, "%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) { |
496 |
|
|
ELL_3V_SUB(_tmpv, _tmpv, point->pos); |
497 |
|
|
fprintf(stderr, " -> moved to %g %g %g %g\n", |
498 |
|
|
point->pos[0], point->pos[1], point->pos[2], point->pos[3]); |
499 |
|
|
fprintf(stderr, " (moved %g)\n", ELL_3V_LEN(_tmpv)); |
500 |
|
|
} |
501 |
|
|
PROBE(val, grad, hess, posproj, negproj); |
502 |
|
|
_pullPointHistAdd(point, pullCondConstraintSatA, val); |
503 |
|
|
if (DEBUG) { |
504 |
|
|
fprintf(stderr, " (+) probed at (%g,%g,%g,%g)\n", |
505 |
|
|
point->pos[0], point->pos[1], point->pos[2], point->pos[3]); |
506 |
|
|
PRINT("after move"); |
507 |
|
|
fprintf(stderr, " 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) { |
514 |
|
|
fprintf(stderr, " (+) 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) { |
521 |
|
|
fprintf(stderr, " (+) 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) { |
529 |
|
|
fprintf(stderr, " 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)) { |
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) { |
552 |
|
|
fprintf(stderr, "!%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) : AIR_MAX(-stepMax, step); |
556 |
|
|
convtestB: |
557 |
|
|
if (d2 < 0 && AIR_ABS(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) { |
561 |
|
|
fprintf(stderr, " |step| %g < %g*%g = %g ==> converged!\n", |
562 |
|
|
AIR_ABS(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) { |
571 |
|
|
fprintf(stderr, " -> step %g, |pdir| = %g\n", |
572 |
|
|
step, ELL_3V_LEN(pdir)); |
573 |
|
|
ELL_3V_COPY(_tmpv, point->pos); |
574 |
|
|
fprintf(stderr, " -> 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); |
580 |
|
|
if (!ELL_4V_EXISTS(point->pos)) { |
581 |
|
|
biffAddf(PULL, "%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) { |
588 |
|
|
ELL_3V_SUB(_tmpv, _tmpv, point->pos); |
589 |
|
|
fprintf(stderr, " -> moved to %g %g %g %g\n", |
590 |
|
|
point->pos[0], point->pos[1], point->pos[2], point->pos[3]); |
591 |
|
|
fprintf(stderr, " (moved %g)\n", ELL_3V_LEN(_tmpv)); |
592 |
|
|
} |
593 |
|
|
PROBE(val, grad, hess, posproj, negproj); |
594 |
|
|
_pullPointHistAdd(point, pullCondConstraintSatB, val); |
595 |
|
|
if (DEBUG) { |
596 |
|
|
fprintf(stderr, " (-) probed at (%g,%g,%g,%g)\n", |
597 |
|
|
point->pos[0], point->pos[1], point->pos[2], point->pos[3]); |
598 |
|
|
PRINT("after move"); |
599 |
|
|
fprintf(stderr, " 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) { |
606 |
|
|
fprintf(stderr, " (-) 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) { |
613 |
|
|
fprintf(stderr, " (-) 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) { |
621 |
|
|
fprintf(stderr, " 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_FALSE; |
631 |
|
|
} |
632 |
|
|
_pullPointHistAdd(point, (*constrFailP |
633 |
|
|
? pullCondConstraintFail |
634 |
|
|
: pullCondConstraintSuccess), AIR_NAN); |
635 |
|
|
if (DEBUG) { |
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 |
|
|
|
648 |
|
|
double |
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]) |
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 |
|
|
*/ |
667 |
|
|
int |
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); |
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) { |
684 |
|
|
fprintf(stderr, "!%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(PULL, "%s: trouble", me); |
694 |
|
|
return 1; |
695 |
|
|
} |
696 |
|
|
break; |
697 |
|
|
case pullInfoIsovalue: |
698 |
|
|
if (constraintSatIso(task, point, stepMax, constrEps, |
699 |
|
|
iterMax, constrFailP)) { |
700 |
|
|
biffAddf(PULL, "%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(PULL, "%s: trouble", me); |
712 |
|
|
return 1; |
713 |
|
|
} |
714 |
|
|
break; |
715 |
|
|
default: |
716 |
|
|
fprintf(stderr, "%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); |
721 |
|
|
if (travelMax) { |
722 |
|
|
travel = ELL_3V_LEN(pos3Diff)/task->pctx->voxelSizeSpace; |
723 |
|
|
if (travel > travelMax) { |
724 |
|
|
*constrFailP = pullConstraintFailTravel; |
725 |
|
|
if (DEBUG) { |
726 |
|
|
fprintf(stderr, "!%s: travel %g > travelMax %g\n", me, |
727 |
|
|
travel, travelMax); |
728 |
|
|
} |
729 |
|
|
} |
730 |
|
|
} |
731 |
|
|
if (DEBUG) { |
732 |
|
|
fprintf(stderr, "!%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_PHIST |
741 |
|
|
if (1) { |
742 |
|
|
Nrrd *nhist; |
743 |
|
|
FILE *fhist; |
744 |
|
|
char fname[AIR_STRLEN_LARGE]; |
745 |
|
|
nhist = nrrdNew(); |
746 |
|
|
sprintf(fname, "%04u-%04u-phist.nrrd", task->pctx->iter, point->idtag); |
747 |
|
|
if (pullPositionHistoryNrrdGet(nhist, task->pctx, point)) { |
748 |
|
|
biffAddf(PULL, "%s: trouble", me); |
749 |
|
|
return 1; |
750 |
|
|
} |
751 |
|
|
if ((fhist = fopen(fname, "w"))) { |
752 |
|
|
if (nrrdSave(fname, nhist, NULL)) { |
753 |
|
|
biffMovef(PULL, NRRD, "%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 |
|
|
*/ |
778 |
|
|
void |
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); /* 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); |
795 |
|
|
ELL_3M_SUB(proj, proj, negproj); |
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); |
802 |
|
|
} else { |
803 |
|
|
pullPointScalar(task->pctx, point, pullInfoIsovalue, vec, NULL); |
804 |
|
|
} |
805 |
|
|
ELL_3V_NORM(nvec, vec, len); |
806 |
|
|
if (len) { |
807 |
|
|
/* .. or and subracting out tensor product of normal with itself */ |
808 |
|
|
ELL_3MV_OUTER(outer, nvec, nvec); |
809 |
|
|
ELL_3M_SUB(proj, proj, outer); |
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 |
|
|
*/ |
826 |
|
|
int |
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(PULL, "%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(PULL, "%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(PULL, "%s: constraint on %s (%d) unimplemented", me, |
877 |
|
|
airEnumStr(pullInfo, pctx->constraint), pctx->constraint); |
878 |
|
|
return -1; |
879 |
|
|
} |
880 |
|
|
return ret; |
881 |
|
|
} |
882 |
|
|
|