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 #include "pull.h"
25 #include "privatePull.h"
26
27 /* --------------------------------------------- */
28
29 unsigned int
30 _pullInfoLen[PULL_INFO_MAX+1] = {
31 0, /* pullInfoUnknown */
32 7, /* pullInfoTensor */
33 7, /* pullInfoTensorInverse */
34 9, /* pullInfoHessian */
35 1, /* pullInfoInside */
36 3, /* pullInfoInsideGradient */
37 1, /* pullInfoHeight */
38 3, /* pullInfoHeightGradient */
39 9, /* pullInfoHeightHessian */
40 1, /* pullInfoHeightLaplacian */
41 1, /* pullInfoSeedPreThresh */
42 1, /* pullInfoSeedThresh */
43 1, /* pullInfoLiveThresh */
44 1, /* pullInfoLiveThresh2 */
45 1, /* pullInfoLiveThresh3 */
46 3, /* pullInfoTangent1 */
47 3, /* pullInfoTangent2 */
48 3, /* pullInfoNegativeTangent1 */
49 3, /* pullInfoNegativeTangent2 */
50 1, /* pullInfoIsovalue */
51 3, /* pullInfoIsovalueGradient */
52 9, /* pullInfoIsovalueHessian */
53 1, /* pullInfoStrength */
54 1, /* pullInfoQuality */
55 };
56
57 unsigned int
pullInfoLen(int info)58 pullInfoLen(int info) {
59 unsigned int ret;
60
61 if (!airEnumValCheck(pullInfo, info)) {
62 ret = _pullInfoLen[info];
63 } else {
64 ret = 0;
65 }
66 return ret;
67 }
68
69 unsigned int
pullPropLen(int prop)70 pullPropLen(int prop) {
71 unsigned int ret;
72
73 switch (prop) {
74 case pullPropIdtag:
75 case pullPropIdCC:
76 case pullPropEnergy:
77 case pullPropStepEnergy:
78 case pullPropStepConstr:
79 case pullPropStuck:
80 case pullPropNeighDistMean:
81 case pullPropScale:
82 case pullPropStability:
83 ret = 1;
84 break;
85 case pullPropPosition:
86 case pullPropForce:
87 ret = 4;
88 break;
89 case pullPropNeighCovar:
90 ret = 10;
91 break;
92 case pullPropNeighCovar7Ten:
93 ret = 7;
94 break;
95 case pullPropNeighTanCovar:
96 ret = 6;
97 break;
98 default:
99 ret = 0;
100 break;
101 }
102 return ret;
103 }
104
105 pullInfoSpec *
pullInfoSpecNew(void)106 pullInfoSpecNew(void) {
107 pullInfoSpec *ispec;
108
109 ispec = AIR_CAST(pullInfoSpec *, calloc(1, sizeof(pullInfoSpec)));
110 if (ispec) {
111 ispec->info = pullInfoUnknown;
112 ispec->source = pullSourceUnknown;
113 ispec->volName = NULL;
114 ispec->item = 0; /* should be the unknown item for any kind */
115 ispec->prop = pullPropUnknown;
116 ispec->scale = AIR_NAN;
117 ispec->zero = AIR_NAN;
118 ispec->constraint = AIR_FALSE;
119 ispec->volIdx = UINT_MAX;
120 }
121 return ispec;
122 }
123
124 pullInfoSpec *
pullInfoSpecNix(pullInfoSpec * ispec)125 pullInfoSpecNix(pullInfoSpec *ispec) {
126
127 if (ispec) {
128 airFree(ispec->volName);
129 airFree(ispec);
130 }
131 return NULL;
132 }
133
134 int
pullInfoSpecAdd(pullContext * pctx,pullInfoSpec * ispec)135 pullInfoSpecAdd(pullContext *pctx, pullInfoSpec *ispec) {
136 static const char me[]="pullInfoSpecAdd";
137 unsigned int ii, vi, haveLen, needLen;
138 const gageKind *kind;
139
140 if (!( pctx && ispec )) {
141 biffAddf(PULL, "%s: got NULL pointer", me);
142 return 1;
143 }
144 if (airEnumValCheck(pullInfo, ispec->info)) {
145 biffAddf(PULL, "%s: %d not a valid %s value", me,
146 ispec->info, pullInfo->name);
147 return 1;
148 }
149 if (airEnumValCheck(pullSource, ispec->source)) {
150 biffAddf(PULL, "%s: %d not a valid %s value", me,
151 ispec->source, pullSource->name);
152 return 1;
153 }
154 if (pctx->ispec[ispec->info]) {
155 biffAddf(PULL, "%s: already set info %s (%d)", me,
156 airEnumStr(pullInfo, ispec->info), ispec->info);
157 return 1;
158 }
159 for (ii=0; ii<=PULL_INFO_MAX; ii++) {
160 if (pctx->ispec[ii] == ispec) {
161 biffAddf(PULL, "%s(%s): already got ispec %p as ispec[%u]", me,
162 airEnumStr(pullInfo, ispec->info), AIR_VOIDP(ispec), ii);
163 return 1;
164 }
165 }
166 if (pctx->verbose) {
167 printf("%s: ispec %s from vol %s\n", me,
168 airEnumStr(pullInfo, ispec->info), ispec->volName);
169 }
170 needLen = pullInfoLen(ispec->info);
171 if (pullSourceGage == ispec->source) {
172 vi = _pullVolumeIndex(pctx, ispec->volName);
173 if (UINT_MAX == vi) {
174 biffAddf(PULL, "%s(%s): no volume has name \"%s\"", me,
175 airEnumStr(pullInfo, ispec->info), ispec->volName);
176 return 1;
177 }
178 kind = pctx->vol[vi]->kind;
179 if (airEnumValCheck(kind->enm, ispec->item)) {
180 biffAddf(PULL, "%s(%s): %d not a valid \"%s\" item", me,
181 airEnumStr(pullInfo, ispec->info), ispec->item, kind->name);
182 return 1;
183 }
184 haveLen = kind->table[ispec->item].answerLength;
185 if (needLen != haveLen) {
186 biffAddf(PULL, "%s(%s): need len %u, but \"%s\" item \"%s\" has len %u",
187 me, airEnumStr(pullInfo, ispec->info), needLen,
188 kind->name, airEnumStr(kind->enm, ispec->item), haveLen);
189 return 1;
190 }
191 /* very tricky: seedOnly is initialized to true for everything */
192 if (pullInfoSeedThresh != ispec->info
193 && pullInfoSeedPreThresh != ispec->info) {
194 /* if the info is neither seedthresh nor seedprethresh, then the
195 volume will have to be probed after the first iter, so turn
196 *off* seedOnly */
197 pctx->vol[vi]->seedOnly = AIR_FALSE;
198 }
199 /* less tricky: turn on forSeedPreThresh as needed;
200 its initialized to false */
201 if (pullInfoSeedPreThresh == ispec->info) {
202 pctx->vol[vi]->forSeedPreThresh = AIR_TRUE;
203 if (pctx->verbose) {
204 printf("%s: volume %u %s used for %s\n", me, vi, pctx->vol[vi]->name,
205 airEnumStr(pullInfo, pullInfoSeedPreThresh));
206 }
207 }
208 /* now set item in gage query */
209 if (gageQueryItemOn(pctx->vol[vi]->gctx, pctx->vol[vi]->gpvl,
210 ispec->item)) {
211 biffMovef(PULL, GAGE, "%s: trouble adding item %u to vol %u", me,
212 ispec->item, vi);
213 return 1;
214 }
215 ispec->volIdx = vi;
216 } else if (pullSourceProp == ispec->source) {
217 haveLen = pullPropLen(ispec->prop);
218 if (needLen != haveLen) {
219 biffAddf(PULL, "%s: need len %u, but \"%s\" \"%s\" has len %u",
220 me, needLen, pullProp->name,
221 airEnumStr(pullProp, ispec->prop), haveLen);
222 return 1;
223 }
224
225 } else {
226 biffAddf(PULL, "%s: sorry, source %s unsupported", me,
227 airEnumStr(pullSource, ispec->source));
228 return 1;
229 }
230 if (haveLen > 9) {
231 biffAddf(PULL, "%s: sorry, answer length (%u) > 9 unsupported", me,
232 haveLen);
233 return 1;
234 }
235
236 pctx->ispec[ispec->info] = ispec;
237
238 return 0;
239 }
240
241 /*
242 ** sets:
243 ** pctx->infoIdx[]
244 ** pctx->infoTotalLen
245 ** pctx->constraint
246 ** pctx->constraintDim
247 ** pctx->targetDim (non-trivial logic for scale-space!)
248 */
249 int
_pullInfoSetup(pullContext * pctx)250 _pullInfoSetup(pullContext *pctx) {
251 static const char me[]="_pullInfoSetup";
252 unsigned int ii;
253
254 pctx->infoTotalLen = 0;
255 pctx->constraint = 0;
256 pctx->constraintDim = 0;
257 for (ii=0; ii<=PULL_INFO_MAX; ii++) {
258 if (pctx->ispec[ii]) {
259 pctx->infoIdx[ii] = pctx->infoTotalLen;
260 if (pctx->verbose) {
261 printf("!%s: infoIdx[%u] (%s) = %u\n", me,
262 ii, airEnumStr(pullInfo, ii), pctx->infoIdx[ii]);
263 }
264 pctx->infoTotalLen += pullInfoLen(ii);
265 if (!pullInfoLen(ii)) {
266 biffAddf(PULL, "%s: got zero-length answer for ispec[%u]", me, ii);
267 return 1;
268 }
269 if (pctx->ispec[ii]->constraint) {
270 /* pullVolume *cvol; */
271 pctx->constraint = ii;
272 /* cvol = pctx->vol[pctx->ispec[ii]->volIdx]; */
273 }
274 }
275 }
276 if (pctx->constraint) {
277 pctx->constraintDim = _pullConstraintDim(pctx);
278 if (-1 == pctx->constraintDim) {
279 biffAddf(PULL, "%s: problem learning constraint dimension", me);
280 return 1;
281 }
282 if (!pctx->flag.allowCodimension3Constraints && !pctx->constraintDim) {
283 biffAddf(PULL, "%s: got constr dim 0 but co-dim 3 not allowed", me);
284 return 1;
285 }
286 if (pctx->haveScale) {
287 double *parmS, denS,
288 (*evalS)(double *, double, const double parm[PULL_ENERGY_PARM_NUM]);
289 switch (pctx->interType) {
290 case pullInterTypeUnivariate:
291 pctx->targetDim = 1 + pctx->constraintDim;
292 break;
293 case pullInterTypeSeparable:
294 /* HEY! need to check if this is true given enr and ens! */
295 pctx->targetDim = pctx->constraintDim;
296 break;
297 case pullInterTypeAdditive:
298 parmS = pctx->energySpecS->parm;
299 evalS = pctx->energySpecS->energy->eval;
300 evalS(&denS, _PULL_TARGET_DIM_S_PROBE, parmS);
301 if (denS > 0) {
302 /* at small positive s, derivative was positive ==> attractive */
303 pctx->targetDim = pctx->constraintDim;
304 } else {
305 /* derivative was negative ==> repulsive */
306 pctx->targetDim = 1 + pctx->constraintDim;
307 }
308 break;
309 default:
310 biffAddf(PULL, "%s: sorry, intertype %s not handled here", me,
311 airEnumStr(pullInterType, pctx->interType));
312 break;
313 }
314 } else {
315 pctx->targetDim = pctx->constraintDim;
316 }
317 } else {
318 pctx->constraintDim = 0;
319 pctx->targetDim = 0;
320 }
321 if (pctx->verbose) {
322 printf("!%s: infoTotalLen=%u, constr=%d, constr,targetDim = %d,%d\n",
323 me, pctx->infoTotalLen, pctx->constraint,
324 pctx->constraintDim, pctx->targetDim);
325 }
326 return 0;
327 }
328
329 static void
_infoCopy1(double * dst,const double * src)330 _infoCopy1(double *dst, const double *src) {
331 dst[0] = src[0];
332 }
333
334 static void
_infoCopy2(double * dst,const double * src)335 _infoCopy2(double *dst, const double *src) {
336 dst[0] = src[0]; dst[1] = src[1];
337 }
338
339 static void
_infoCopy3(double * dst,const double * src)340 _infoCopy3(double *dst, const double *src) {
341 dst[0] = src[0]; dst[1] = src[1]; dst[2] = src[2];
342 }
343
344 static void
_infoCopy4(double * dst,const double * src)345 _infoCopy4(double *dst, const double *src) {
346 dst[0] = src[0]; dst[1] = src[1]; dst[2] = src[2]; dst[3] = src[3];
347 }
348
349 static void
_infoCopy5(double * dst,const double * src)350 _infoCopy5(double *dst, const double *src) {
351 dst[0] = src[0]; dst[1] = src[1]; dst[2] = src[2]; dst[3] = src[3];
352 dst[4] = src[4];
353 }
354
355 static void
_infoCopy6(double * dst,const double * src)356 _infoCopy6(double *dst, const double *src) {
357 dst[0] = src[0]; dst[1] = src[1]; dst[2] = src[2]; dst[3] = src[3];
358 dst[4] = src[4]; dst[5] = src[5];
359 }
360
361 static void
_infoCopy7(double * dst,const double * src)362 _infoCopy7(double *dst, const double *src) {
363 dst[0] = src[0]; dst[1] = src[1]; dst[2] = src[2]; dst[3] = src[3];
364 dst[4] = src[4]; dst[5] = src[5]; dst[6] = src[6];
365 }
366
367 static void
_infoCopy8(double * dst,const double * src)368 _infoCopy8(double *dst, const double *src) {
369 dst[0] = src[0]; dst[1] = src[1]; dst[2] = src[2]; dst[3] = src[3];
370 dst[4] = src[4]; dst[5] = src[5]; dst[6] = src[6]; dst[7] = src[7];
371 }
372
373 static void
_infoCopy9(double * dst,const double * src)374 _infoCopy9(double *dst, const double *src) {
375 dst[0] = src[0]; dst[1] = src[1]; dst[2] = src[2]; dst[3] = src[3];
376 dst[4] = src[4]; dst[5] = src[5]; dst[6] = src[6]; dst[7] = src[7];
377 dst[8] = src[8];
378 }
379
380 void (*_pullInfoCopy[10])(double *, const double *) = {
381 NULL,
382 _infoCopy1,
383 _infoCopy2,
384 _infoCopy3,
385 _infoCopy4,
386 _infoCopy5,
387 _infoCopy6,
388 _infoCopy7,
389 _infoCopy8,
390 _infoCopy9
391 };
392
393 int
pullInfoGet(Nrrd * ninfo,int info,pullContext * pctx)394 pullInfoGet(Nrrd *ninfo, int info, pullContext *pctx) {
395 static const char me[]="pullInfoGet";
396 size_t size[2];
397 unsigned int dim, pointNum, pointIdx, binIdx, outIdx, alen, aidx;
398 double *out_d;
399 pullBin *bin;
400 pullPoint *point;
401
402 if (airEnumValCheck(pullInfo, info)) {
403 biffAddf(PULL, "%s: info %d not valid", me, info);
404 return 1;
405 }
406 pointNum = pullPointNumber(pctx);
407 alen = pullInfoLen(info);
408 aidx = pctx->infoIdx[info];
409 if (1 == alen) {
410 dim = 1;
411 size[0] = pointNum;
412 } else {
413 dim = 2;
414 size[0] = alen;
415 size[1] = pointNum;
416 }
417 if (nrrdMaybeAlloc_nva(ninfo, nrrdTypeDouble, dim, size)) {
418 biffMovef(PULL, NRRD, "%s: trouble allocating output", me);
419 return 1;
420 }
421 out_d = AIR_CAST(double *, ninfo->data);
422
423 outIdx = 0;
424 for (binIdx=0; binIdx<pctx->binNum; binIdx++) {
425 bin = pctx->bin + binIdx;
426 for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) {
427 point = bin->point[pointIdx];
428 _pullInfoCopy[alen](out_d + outIdx, point->info + aidx);
429 outIdx += alen;
430 }
431 }
432
433 return 0;
434 }
435
436 /* HEY this was written in a hurry;
437 ** needs to be checked against parsing code */
438 int
pullInfoSpecSprint(char str[AIR_STRLEN_LARGE],const pullContext * pctx,const pullInfoSpec * ispec)439 pullInfoSpecSprint(char str[AIR_STRLEN_LARGE],
440 const pullContext *pctx, const pullInfoSpec *ispec) {
441 static const char me[]="pullInfoSpecSprint";
442 const pullVolume *pvol;
443 char stmp[AIR_STRLEN_LARGE];
444
445 if (!( str && pctx && ispec )) {
446 biffAddf(PULL, "%s: got NULL pointer", me);
447 return 1;
448 }
449 strcpy(str, "");
450 /* HEY: no bounds checking! */
451 strcat(str, airEnumStr(pullInfo, ispec->info));
452 if (ispec->constraint) {
453 strcat(str, "-c");
454 }
455 strcat(str, ":");
456 if (pullSourceGage == ispec->source) {
457 if (UINT_MAX == ispec->volIdx) {
458 biffAddf(PULL, "%s: never learned volIdx for \"%s\"", me,
459 ispec->volName);
460 return 1;
461 }
462 strcat(str, ispec->volName);
463 strcat(str, ":");
464 pvol = pctx->vol[ispec->volIdx];
465 strcat(str, airEnumStr(pvol->kind->enm, ispec->item));
466 } else if (pullSourceProp == ispec->source) {
467 strcat(str, airEnumStr(pullProp, ispec->prop));
468 } else {
469 biffAddf(PULL, "%s: unexplained source %d", me, ispec->source);
470 return 1;
471 }
472 if ( (pullSourceGage == ispec->source
473 && 1 == pullInfoLen(ispec->info))
474 ||
475 (pullSourceProp == ispec->source
476 && 1 == pullPropLen(ispec->prop)) ) {
477 sprintf(stmp, "%g", ispec->zero);
478 strcat(str, stmp);
479 strcat(str, ":");
480 sprintf(stmp, "%g", ispec->scale);
481 strcat(str, stmp);
482 }
483 return 0;
484 }
485
486