GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: src/seek/updateSeek.c Lines: 0 533 0.0 %
Date: 2017-05-26 Branches: 0 374 0.0 %

Line Branch Exec Source
1
/*
2
  Teem: Tools to process and visualize scientific data and images             .
3
  Copyright (C) 2009  Thomas Schultz
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 "seek.h"
25
#include "privateSeek.h"
26
27
static int
28
updateNinEtAl(seekContext *sctx) {
29
  static const char me[]="updateNinEtAl";
30
31
  if (sctx->verbose > 5) {
32
    fprintf(stderr, "%s: --------------------\n", me);
33
    fprintf(stderr, "%s: flagData = %d\n", me, sctx->flag[flagData]);
34
  }
35
36
  if (!( sctx->ninscl || sctx->pvl )) {
37
    biffAddf(SEEK, "%s: data never set", me);
38
    return 1;
39
  }
40
41
  if (sctx->flag[flagData]) {
42
    if (sctx->ninscl) {
43
      sctx->nin = sctx->ninscl;
44
      sctx->baseDim = 0;
45
      if (gageShapeSet(sctx->_shape, sctx->ninscl, 0)) {
46
        biffMovef(SEEK, GAGE, "%s: trouble with scalar volume", me);
47
        return 1;
48
      }
49
      sctx->shape = sctx->_shape;
50
    } else {
51
      sctx->nin = sctx->pvl->nin;
52
      sctx->baseDim = sctx->pvl->kind->baseDim;
53
      sctx->shape = sctx->gctx->shape;
54
    }
55
    sctx->flag[flagData] = AIR_FALSE;
56
    sctx->flag[flagNinEtAl] = AIR_TRUE;
57
  }
58
  return 0;
59
}
60
61
static int
62
updateAnswerPointers(seekContext *sctx) {
63
  static const char me[]="updateAnswerPointers";
64
65
  if (sctx->verbose > 5) {
66
    fprintf(stderr, "%s: --------------------\n", me);
67
    fprintf(stderr, "%s: flagItemValue = %d\n", me,
68
            sctx->flag[flagItemValue]);
69
    fprintf(stderr, "%s: flagItemStrength = %d\n", me,
70
            sctx->flag[flagItemStrength]);
71
    fprintf(stderr, "%s: flagItemNormal = %d\n", me,
72
            sctx->flag[flagItemNormal]);
73
    fprintf(stderr, "%s: flagItemGradient = %d\n", me,
74
            sctx->flag[flagItemGradient]);
75
    fprintf(stderr, "%s: flagItemEigensystem = %d\n", me,
76
            sctx->flag[flagItemEigensystem]);
77
    fprintf(stderr, "%s: flagItemHess = %d\n", me,
78
            sctx->flag[flagItemHess]);
79
    fprintf(stderr, "%s: flagNinEtAl = %d\n", me,
80
            sctx->flag[flagNinEtAl]);
81
    fprintf(stderr, "%s: flagNormalsFind = %d\n", me,
82
            sctx->flag[flagNormalsFind]);
83
    fprintf(stderr, "%s: flagStrengthUse = %d\n", me,
84
            sctx->flag[flagStrengthUse]);
85
    fprintf(stderr, "%s: flagType = %d\n", me,
86
            sctx->flag[flagType]);
87
    fprintf(stderr, "%s: flagData = %d\n", me,
88
            sctx->flag[flagData]);
89
  }
90
91
  if (seekTypeUnknown == sctx->type) {
92
    biffAddf(SEEK, "%s: feature type never set", me);
93
    return 1;
94
  }
95
96
  if (sctx->flag[flagItemValue]
97
      || sctx->flag[flagItemStrength]
98
      || sctx->flag[flagItemNormal]
99
      || sctx->flag[flagItemGradient]
100
      || sctx->flag[flagItemEigensystem]
101
      || sctx->flag[flagItemHess]
102
      || sctx->flag[flagNinEtAl]
103
      || sctx->flag[flagNormalsFind]
104
      || sctx->flag[flagStrengthUse]
105
      || sctx->flag[flagType]) {
106
107
    /* this is apt regardless of feature type */
108
    if (sctx->strengthUse) {
109
      if (-1 == sctx->stngItem) {
110
        biffAddf(SEEK, "%s: need to set strength item to use strength", me);
111
        return 1;
112
      }
113
      sctx->stngAns = (gageAnswerPointer(sctx->gctx, sctx->pvl,
114
                                         sctx->stngItem));
115
    } else {
116
      sctx->stngAns = NULL;
117
    }
118
119
    switch (sctx->type) {
120
    case seekTypeIsocontour:
121
      if (!( sctx->ninscl || -1 != sctx->sclvItem )) {
122
        biffAddf(SEEK,
123
                 "%s: need either scalar volume or value item set for %s",
124
                 me, airEnumStr(seekType, seekTypeIsocontour));
125
        return 1;
126
      }
127
      if (sctx->normalsFind) {
128
        /* NOTE simplifying assumption described in seek.h */
129
        if (!( sctx->ninscl || -1 != sctx->normItem )) {
130
          biffAddf(SEEK, "%s: need either scalar volume or "
131
                   "normal item set for normals for %s",
132
                   me, airEnumStr(seekType, seekTypeIsocontour));
133
          return 1;
134
        }
135
      }
136
      if (sctx->ninscl) {
137
        sctx->sclvAns = NULL;
138
        sctx->normAns = NULL;
139
      } else {
140
        sctx->sclvAns = gageAnswerPointer(sctx->gctx, sctx->pvl,
141
                                          sctx->sclvItem);
142
        sctx->normAns = (sctx->normalsFind
143
                         ? gageAnswerPointer(sctx->gctx, sctx->pvl,
144
                                             sctx->normItem)
145
                         : NULL);
146
      }
147
      if (sctx->flag[flagItemGradient]
148
          || sctx->flag[flagItemEigensystem]
149
          || sctx->flag[flagItemHess]) {
150
        biffAddf(SEEK,
151
                 "%s: can't set gradient, Hessian, or eigensystem for %s",
152
                 me, airEnumStr(seekType, seekTypeIsocontour));
153
        return 1;
154
      }
155
      sctx->gradAns = NULL;
156
      sctx->evalAns = NULL;
157
      sctx->evecAns = NULL;
158
      sctx->hessAns = NULL;
159
      break;
160
    case seekTypeRidgeSurface:
161
    case seekTypeValleySurface:
162
    case seekTypeMaximalSurface:
163
    case seekTypeMinimalSurface:
164
    case seekTypeRidgeSurfaceOP:
165
    case seekTypeRidgeSurfaceT:
166
    case seekTypeValleySurfaceOP:
167
    case seekTypeValleySurfaceT:
168
      if ( !sctx->pvl ) {
169
        biffAddf(SEEK, "%s: can't find %s without a gage context",
170
                 me, airEnumStr(seekType, sctx->type));
171
        return 1;
172
      }
173
      if (!( -1 != sctx->gradItem
174
             && -1 != sctx->evalItem
175
             && -1 != sctx->evecItem )) {
176
        biffAddf(SEEK, "%s: grad, eval, evec items not all set", me);
177
        return 1;
178
      }
179
      if ((sctx->type==seekTypeRidgeSurfaceT ||
180
           sctx->type==seekTypeValleySurfaceT) &&
181
          -1 == sctx->hessItem) {
182
        biffAddf(SEEK, "%s: hess item not set", me);
183
        return 1;
184
      }
185
      if (sctx->normalsFind) {
186
        /* NOTE simplifying assumption described in seek.h */
187
        if (-1 == sctx->normItem) {
188
          biffAddf(SEEK, "%s: need normal item set for normals for %s",
189
                   me, airEnumStr(seekType, sctx->type));
190
          return 1;
191
        }
192
        sctx->normAns = (gageAnswerPointer(sctx->gctx, sctx->pvl,
193
                                           sctx->normItem));
194
      } else {
195
        sctx->normAns = NULL;
196
      }
197
      sctx->sclvAns = NULL;
198
      sctx->gradAns = gageAnswerPointer(sctx->gctx, sctx->pvl, sctx->gradItem);
199
      sctx->evalAns = gageAnswerPointer(sctx->gctx, sctx->pvl, sctx->evalItem);
200
      sctx->evecAns = gageAnswerPointer(sctx->gctx, sctx->pvl, sctx->evecItem);
201
      if (sctx->type==seekTypeRidgeSurfaceT ||
202
          sctx->type==seekTypeValleySurfaceT)
203
        sctx->hessAns = gageAnswerPointer(sctx->gctx, sctx->pvl,
204
                                          sctx->hessItem);
205
      else
206
        sctx->hessAns = NULL;
207
      break;
208
    default:
209
      biffAddf(SEEK, "%s: sorry, %s extraction not implemented", me,
210
               airEnumStr(seekType, sctx->type));
211
      return 1;
212
    }
213
214
    sctx->flag[flagItemValue] = AIR_FALSE;
215
    sctx->flag[flagItemStrength] = AIR_FALSE;
216
    sctx->flag[flagItemNormal] = AIR_FALSE;
217
    sctx->flag[flagItemGradient] = AIR_FALSE;
218
    sctx->flag[flagItemEigensystem] = AIR_FALSE;
219
    sctx->flag[flagItemHess] = AIR_FALSE;
220
    sctx->flag[flagNormalsFind] = AIR_FALSE;
221
    sctx->flag[flagAnswerPointers] = AIR_TRUE;
222
  }
223
224
  return 0;
225
}
226
227
static int
228
updateSxSySz(seekContext *sctx) {
229
  static const char me[]="updateSxSySz";
230
  size_t sizeIn[3], sizeOut[3];
231
  double min, max, scl[3], off[3];
232
  unsigned int axi;
233
234
  if (sctx->verbose > 5) {
235
    fprintf(stderr, "%s: --------------------\n", me);
236
    fprintf(stderr, "%s: flagSamples = %d\n", me, sctx->flag[flagSamples]);
237
    fprintf(stderr, "%s: flagNinEtAl = %d\n", me, sctx->flag[flagNinEtAl]);
238
  }
239
240
  sizeIn[0] = sctx->nin->axis[sctx->baseDim+0].size;
241
  sizeIn[1] = sctx->nin->axis[sctx->baseDim+1].size;
242
  sizeIn[2] = sctx->nin->axis[sctx->baseDim+2].size;
243
244
  if (sctx->flag[flagSamples]
245
      || sctx->flag[flagNinEtAl]) {
246
    if (0 == sctx->samples[0]
247
        || 0 == sctx->samples[1]
248
        || 0 == sctx->samples[2]) {
249
      ELL_3V_COPY(sizeOut, sizeIn);
250
    } else {
251
      if (!sctx->pvl) {
252
        biffAddf(SEEK,
253
                 "%s: can't specify # samples (%u,%u,%u) independent of "
254
                 "volume dimensions (%u,%u,%u) without a gage context", me,
255
                 AIR_CAST(unsigned int, sctx->samples[0]),
256
                 AIR_CAST(unsigned int, sctx->samples[1]),
257
                 AIR_CAST(unsigned int, sctx->samples[2]),
258
                 AIR_CAST(unsigned int, sizeIn[0]),
259
                 AIR_CAST(unsigned int, sizeIn[1]),
260
                 AIR_CAST(unsigned int, sizeIn[2]));
261
        return 1;
262
      }
263
      ELL_3V_COPY(sizeOut, sctx->samples);
264
    }
265
    /* want to make absolutely sure txfIdx was being set ...
266
       if (sctx->sx != sizeOut[0]
267
       || sctx->sy != sizeOut[1]
268
       || sctx->sz != sizeOut[2]) { */
269
    sctx->sx = sizeOut[0];
270
    sctx->sy = sizeOut[1];
271
    sctx->sz = sizeOut[2];
272
    /* there has got to be a better way of doing this... */
273
    /* perhaps refer to the way origin is calculated in new nrrd resampler? */
274
    for (axi=0; axi<3; axi++) {
275
      min = (nrrdCenterCell == sctx->shape->center
276
             ? -0.5
277
             : 0.0);
278
      max = (nrrdCenterCell == sctx->shape->center
279
             ? sizeIn[axi] - 0.5
280
             : sizeIn[axi] - 1.0);
281
      off[axi] = NRRD_POS(sctx->shape->center, min, max, sizeOut[axi], 0);
282
      scl[axi] = (NRRD_POS(sctx->shape->center, min, max, sizeOut[axi], 1)
283
                  - off[axi]);
284
    }
285
    ELL_4V_SET(sctx->txfIdx + 0*4, scl[0],    0.0,    0.0, off[0]);
286
    ELL_4V_SET(sctx->txfIdx + 1*4,    0.0, scl[1],    0.0, off[1]);
287
    ELL_4V_SET(sctx->txfIdx + 2*4,    0.0,    0.0, scl[2], off[2]);
288
    ELL_4V_SET(sctx->txfIdx + 3*4,    0.0,    0.0,    0.0,    1.0);
289
    sctx->flag[flagSxSySz] = AIR_TRUE;
290
    sctx->flag[flagSamples] = AIR_FALSE;
291
  }
292
  return 0;
293
}
294
295
static int
296
updateReverse(seekContext *sctx) {
297
  static const char me[]="updateReverse";
298
299
  if (sctx->verbose > 5) {
300
    fprintf(stderr, "%s: --------------------\n", me);
301
    fprintf(stderr, "%s: flagNinEtAl = %d\n", me,
302
            sctx->flag[flagNinEtAl]);
303
    fprintf(stderr, "%s: flagLowerInside = %d\n", me,
304
            sctx->flag[flagLowerInside]);
305
  }
306
307
  if (sctx->flag[flagNinEtAl]
308
      || sctx->flag[flagLowerInside]) {
309
    double mat[9];
310
    int reverse;
311
312
    ELL_34M_EXTRACT(mat, sctx->shape->ItoW);
313
    reverse = (!!sctx->lowerInside) ^ (ELL_3M_DET(mat) < 0);
314
    if (sctx->reverse != reverse) {
315
      sctx->reverse = reverse;
316
      sctx->flag[flagReverse] = AIR_TRUE;
317
    }
318
  }
319
  return 0;
320
}
321
322
static int
323
updateTxfNormal(seekContext *sctx) {
324
  static const char me[]="updateTxfNormal";
325
326
  if (sctx->verbose > 5) {
327
    fprintf(stderr, "%s: --------------------\n", me);
328
    fprintf(stderr, "%s: flagNinEtAl = %d\n", me,
329
            sctx->flag[flagNinEtAl]);
330
    fprintf(stderr, "%s: flagLowerInside = %d\n", me,
331
            sctx->flag[flagLowerInside]);
332
  }
333
334
  if (sctx->flag[flagNinEtAl]
335
      || sctx->flag[flagLowerInside]) {
336
    double matA[9], matB[9];
337
338
    ELL_34M_EXTRACT(matA, sctx->shape->ItoW);
339
    ell_3m_inv_d(matB, matA);
340
    ELL_3M_TRANSPOSE(sctx->txfNormal, matB);
341
    if (!sctx->lowerInside) {
342
      ELL_3M_SCALE(sctx->txfNormal, -1, sctx->txfNormal);
343
    }
344
    sctx->flag[flagTxfNormal] = AIR_TRUE;
345
  }
346
  sctx->flag[flagLowerInside] = AIR_FALSE;
347
  return 0;
348
}
349
350
static int
351
updateSlabCacheAlloc(seekContext *sctx) {
352
  static const char me[]="updateSlabCacheAlloc";
353
  int E;
354
355
  if (sctx->verbose > 5) {
356
    fprintf(stderr, "%s: --------------------\n", me);
357
    fprintf(stderr, "%s: flagType = %d (type = %s)\n", me,
358
            sctx->flag[flagType], airEnumStr(seekType, sctx->type));
359
    fprintf(stderr, "%s: flagStrengthUse = %d\n", me,
360
            sctx->flag[flagStrengthUse]);
361
    fprintf(stderr, "%s: flagSxSySz = %d\n", me,
362
            sctx->flag[flagSxSySz]);
363
  }
364
365
  E = 0;
366
  if (sctx->flag[flagType]
367
      || sctx->flag[flagStrengthUse]  /* kind of sloppy/overkill */
368
      || sctx->flag[flagSxSySz]) {
369
    if (!E) E |= nrrdMaybeAlloc_va(sctx->nvidx, nrrdTypeInt, 3,
370
                                   AIR_CAST(size_t, 15),
371
                                   sctx->sx,
372
                                   sctx->sy);
373
    if (!E) sctx->vidx = AIR_CAST(int*, sctx->nvidx->data);
374
    if (!E) E |= nrrdMaybeAlloc_va(sctx->ntreated, nrrdTypeChar, 2,
375
                                   sctx->sx,
376
                                   sctx->sy);
377
    if (!E) sctx->treated = AIR_CAST(signed char*, sctx->ntreated->data);
378
    if (sctx->strengthUse) {
379
      if (!E) E |= nrrdMaybeAlloc_va(sctx->nstng, nrrdTypeDouble, 3,
380
                                     AIR_CAST(size_t, 2),
381
                                     sctx->sx,
382
                                     sctx->sy);
383
      if (!E) sctx->stng = AIR_CAST(double*, sctx->nstng->data);
384
    }
385
    if (seekTypeIsocontour == sctx->type) {
386
      if (!E) E |= nrrdMaybeAlloc_va(sctx->nsclv, nrrdTypeDouble, 3,
387
                                     AIR_CAST(size_t, 4),
388
                                     sctx->sx + 2,
389
                                     sctx->sy + 2);
390
      if (!E) sctx->sclv = AIR_CAST(double*, sctx->nsclv->data);
391
    }
392
    if (seekTypeRidgeSurface == sctx->type
393
        || seekTypeValleySurface == sctx->type
394
        || seekTypeMaximalSurface == sctx->type
395
        || seekTypeMinimalSurface == sctx->type
396
        || seekTypeRidgeSurfaceOP == sctx->type
397
        || seekTypeRidgeSurfaceT == sctx->type
398
        || seekTypeValleySurfaceOP == sctx->type
399
        || seekTypeValleySurfaceT == sctx->type) {
400
      if (!E) E |= nrrdMaybeAlloc_va(sctx->ngrad, nrrdTypeDouble, 4,
401
                                     AIR_CAST(size_t, 3),
402
                                     AIR_CAST(size_t, 2),
403
                                     sctx->sx,
404
                                     sctx->sy);
405
      if (!E) sctx->grad = AIR_CAST(double*, sctx->ngrad->data);
406
    } else {
407
      sctx->grad = NULL;
408
    }
409
410
    if (seekTypeRidgeSurface == sctx->type
411
        || seekTypeValleySurface == sctx->type
412
        || seekTypeMaximalSurface == sctx->type
413
        || seekTypeMinimalSurface == sctx->type
414
        || seekTypeRidgeSurfaceOP == sctx->type
415
        || seekTypeValleySurfaceOP == sctx->type) {
416
      if (!E) E |= nrrdMaybeAlloc_va(sctx->neval, nrrdTypeDouble, 4,
417
                                     AIR_CAST(size_t, 3),
418
                                     AIR_CAST(size_t, 2),
419
                                     sctx->sx,
420
                                     sctx->sy);
421
      if (!E) sctx->eval = AIR_CAST(double*, sctx->neval->data);
422
      if (!E) E |= nrrdMaybeAlloc_va(sctx->nevec, nrrdTypeDouble, 4,
423
                                     AIR_CAST(size_t, 9),
424
                                     AIR_CAST(size_t, 2),
425
                                     sctx->sx,
426
                                     sctx->sy);
427
      if (!E) sctx->evec = AIR_CAST(double*, sctx->nevec->data);
428
      if (!E) E |= nrrdMaybeAlloc_va(sctx->nflip, nrrdTypeChar, 3,
429
                                     AIR_CAST(size_t, 5),
430
                                     sctx->sx,
431
                                     sctx->sy);
432
      if (!E) sctx->flip = AIR_CAST(signed char*, sctx->nflip->data);
433
    } else {
434
      sctx->eval = NULL;
435
      sctx->evec = NULL;
436
      sctx->flip = NULL;
437
    }
438
439
    if (seekTypeRidgeSurfaceT == sctx->type ||
440
        seekTypeValleySurfaceT == sctx->type) {
441
      if (!E) E |= nrrdMaybeAlloc_va(sctx->nfacevidx, nrrdTypeInt, 3,
442
                                     AIR_CAST(size_t, 4),
443
                                     sctx->sx,
444
                                     sctx->sy);
445
      if (!E) sctx->facevidx = AIR_CAST(int*, sctx->nfacevidx->data);
446
      if (!E) E |= nrrdMaybeAlloc_va(sctx->nhess, nrrdTypeDouble, 4,
447
                                     AIR_CAST(size_t, 9),
448
                                     AIR_CAST(size_t, 2),
449
                                     sctx->sx,
450
                                     sctx->sy);
451
      if (!E) sctx->hess = AIR_CAST(double*, sctx->nhess->data);
452
      if (!E) E |= nrrdMaybeAlloc_va(sctx->nt, nrrdTypeDouble,  4,
453
                                     AIR_CAST(size_t, 9),
454
                                     AIR_CAST(size_t, 2),
455
                                     sctx->sx,
456
                                     sctx->sy);
457
      if (!E) sctx->t = AIR_CAST(double*, sctx->nt->data);
458
      if (!E) E |= nrrdMaybeAlloc_va(sctx->nedgealpha, nrrdTypeDouble, 4,
459
                                     AIR_CAST(size_t, 5),
460
                                     AIR_CAST(size_t, 3),
461
                                     sctx->sx,
462
                                     sctx->sy);
463
      if (!E) sctx->edgealpha = AIR_CAST(double*, sctx->nedgealpha->data);
464
      if (!E) E |= nrrdMaybeAlloc_va(sctx->nedgenorm, nrrdTypeDouble, 4,
465
                                     AIR_CAST(size_t, 5),
466
                                     AIR_CAST(size_t, 9),
467
                                     sctx->sx,
468
                                     sctx->sy);
469
      if (!E) sctx->edgenorm = AIR_CAST(double*, sctx->nedgenorm->data);
470
      if (!E) E |= nrrdMaybeAlloc_va(sctx->nedgeicoord, nrrdTypeDouble, 4,
471
                                     AIR_CAST(size_t, 5),
472
                                     AIR_CAST(size_t, 9),
473
                                     sctx->sx,
474
                                     sctx->sy);
475
      if (!E) sctx->edgeicoord = AIR_CAST(double*, sctx->nedgeicoord->data);
476
      if (!E) E |= nrrdMaybeAlloc_va(sctx->nfacecoord, nrrdTypeDouble, 4,
477
                                     AIR_CAST(size_t, 4),
478
                                     AIR_CAST(size_t, 2),
479
                                     sctx->sx,
480
                                     sctx->sy);
481
      if (!E) sctx->facecoord = AIR_CAST(double*, sctx->nfacecoord->data);
482
      if (!E) E |= nrrdMaybeAlloc_va(sctx->nfacenorm, nrrdTypeDouble, 4,
483
                                     AIR_CAST(size_t, 4),
484
                                     AIR_CAST(size_t, 3),
485
                                     sctx->sx,
486
                                     sctx->sy);
487
      if (!E) sctx->facenorm = AIR_CAST(double*, sctx->nfacenorm->data);
488
      if (!E) E |= nrrdMaybeAlloc_va(sctx->nfaceicoord, nrrdTypeDouble, 4,
489
                                     AIR_CAST(size_t, 4),
490
                                     AIR_CAST(size_t, 3),
491
                                     sctx->sx,
492
                                     sctx->sy);
493
      if (!E) sctx->faceicoord = AIR_CAST(double*, sctx->nfaceicoord->data);
494
      if (!E) E |= nrrdMaybeAlloc_va(sctx->npairs, nrrdTypeChar, 4,
495
                                     AIR_CAST(size_t, 4),
496
                                     AIR_CAST(size_t, 12),
497
                                     sctx->sx,
498
                                     sctx->sy);
499
      if (!E) sctx->pairs = AIR_CAST(signed char*, sctx->npairs->data);
500
      if (!E) E |= nrrdMaybeAlloc_va(sctx->ngradcontext, nrrdTypeDouble, 4,
501
                                     AIR_CAST(size_t, 3),
502
                                     AIR_CAST(size_t, 2),
503
                                     sctx->sx,
504
                                     sctx->sy);
505
      if (!E) sctx->gradcontext = AIR_CAST(double*, sctx->ngradcontext->data);
506
      if (!E) E |= nrrdMaybeAlloc_va(sctx->nhesscontext, nrrdTypeDouble, 4,
507
                                     AIR_CAST(size_t, 9),
508
                                     AIR_CAST(size_t, 2),
509
                                     sctx->sx,
510
                                     sctx->sy);
511
      if (!E) sctx->hesscontext = AIR_CAST(double*, sctx->nhesscontext->data);
512
      if (!E) E |= nrrdMaybeAlloc_va(sctx->ntcontext, nrrdTypeDouble, 4,
513
                                     AIR_CAST(size_t, 9),
514
                                     AIR_CAST(size_t, 2),
515
                                     sctx->sx,
516
                                     sctx->sy);
517
      if (!E) sctx->tcontext = AIR_CAST(double*, sctx->ntcontext->data);
518
      if (!E) E |= nrrdMaybeAlloc_va(sctx->nstngcontext, nrrdTypeDouble, 2,
519
                                     sctx->sx,
520
                                     sctx->sy);
521
      if (!E) sctx->stngcontext = AIR_CAST(double*, sctx->nstngcontext->data);
522
    } else {
523
      sctx->facevidx = NULL;
524
      sctx->hess = NULL;
525
      sctx->t = NULL;
526
      sctx->edgealpha = NULL;
527
      sctx->edgenorm = NULL;
528
      sctx->facecoord = NULL;
529
      sctx->facenorm = NULL;
530
      sctx->pairs = NULL;
531
      sctx->gradcontext = NULL;
532
      sctx->hesscontext = NULL;
533
      sctx->tcontext = NULL;
534
      sctx->stngcontext = NULL;
535
    }
536
    sctx->flag[flagSlabCacheAlloc] = AIR_TRUE;
537
  }
538
  if (E) {
539
    biffMovef(SEEK, NRRD, "%s: couldn't allocate all slab caches", me);
540
    return 1;
541
  }
542
  sctx->flag[flagStrengthUse] = AIR_FALSE;
543
  sctx->flag[flagSxSySz] = AIR_FALSE;
544
  return 0;
545
}
546
547
static int
548
updateSclDerived(seekContext *sctx) {
549
  static const char me[]="updateSclDerived";
550
  char doneStr[AIR_STRLEN_SMALL];
551
  double *scl, idxIn[4], idxOut[4], val;
552
  unsigned int xi, yi, zi;
553
554
  if (sctx->verbose > 5) {
555
    fprintf(stderr, "%s: --------------------\n", me);
556
    fprintf(stderr, "%s: flagType = %d\n", me,
557
            sctx->flag[flagType]);
558
    fprintf(stderr, "%s: flagNinEtAl = %d\n", me,
559
            sctx->flag[flagNinEtAl]);
560
  }
561
562
  if (sctx->flag[flagType]
563
      || sctx->flag[flagNinEtAl]) {
564
    if (!( seekTypeIsocontour == sctx->type
565
           && sctx->pvl )) {
566
      nrrdEmpty(sctx->nsclDerived);
567
    } else {
568
      if (nrrdMaybeAlloc_va(sctx->nsclDerived, nrrdTypeDouble, 3,
569
                            sctx->sx, sctx->sy, sctx->sz)) {
570
        biffMovef(SEEK, NRRD,
571
                  "%s: couldn't allocated derived scalar volume", me);
572
        return 1;
573
      }
574
      scl = AIR_CAST(double*, sctx->nsclDerived->data);
575
      if (sctx->verbose) {
576
        fprintf(stderr, "%s: pre-computing scalar volume ...       ", me);
577
      }
578
      for (zi=0; zi<sctx->sz; zi++) {
579
        if (sctx->verbose) {
580
          fprintf(stderr, "%s", airDoneStr(0, zi, sctx->sz-1, doneStr));
581
          fflush(stderr);
582
        }
583
        for (yi=0; yi<sctx->sy; yi++) {
584
          for (xi=0; xi<sctx->sx; xi++) {
585
            ELL_4V_SET(idxOut, xi, yi, zi, 1.0);
586
            ELL_4MV_MUL(idxIn, sctx->txfIdx, idxOut);
587
            ELL_34V_HOMOG(idxIn, idxIn);
588
            gageProbe(sctx->gctx, idxIn[0], idxIn[1], idxIn[2]);
589
            val = sctx->sclvAns[0];
590
            if (!AIR_EXISTS(val)) {
591
              biffAddf(SEEK, "%s: probed scalar[%u,%u,%u] %g doesn't exist",
592
                       me, xi, yi, zi, val);
593
              return 1;
594
            }
595
            scl[xi + sctx->sx*(yi + sctx->sy*zi)] = val;
596
          }
597
        }
598
      }
599
      if (sctx->verbose) {
600
        fprintf(stderr, "%s\n", airDoneStr(0, zi, sctx->sz-1, doneStr));
601
      }
602
    }
603
    sctx->flag[flagSclDerived] = AIR_TRUE;
604
  }
605
606
  return 0;
607
}
608
609
static int
610
updateSpanSpaceHist(seekContext *sctx) {
611
  static const char me[]="updateSpanSpaceHist";
612
  unsigned int sx, sy, sz, ss, xi, yi, zi, vi, si, minI, maxI, *spanHist;
613
  double min, max, val;
614
  const void *data;
615
  double (*lup)(const void *, size_t);
616
617
  if (sctx->verbose > 5) {
618
    fprintf(stderr, "%s: --------------------\n", me);
619
    fprintf(stderr, "%s: flagType = %d\n", me,
620
            sctx->flag[flagType]);
621
    fprintf(stderr, "%s: flagSclDerived = %d\n", me,
622
            sctx->flag[flagSclDerived]);
623
    fprintf(stderr, "%s: flagNinEtAl = %d\n", me,
624
            sctx->flag[flagNinEtAl]);
625
  }
626
627
  if (sctx->flag[flagType]
628
      || sctx->flag[flagSclDerived]
629
      || sctx->flag[flagNinEtAl]) {
630
    if (seekTypeIsocontour != sctx->type) {
631
      nrrdEmpty(sctx->nspanHist);
632
      sctx->range->min = AIR_NAN;
633
      sctx->range->max = AIR_NAN;
634
    } else {
635
      nrrdRangeSet(sctx->range,
636
                   (sctx->ninscl ? sctx->ninscl : sctx->nsclDerived),
637
                   nrrdBlind8BitRangeFalse);
638
      if (sctx->range->hasNonExist) {
639
        biffAddf(SEEK, "%s: scalar volume has non-existent values", me);
640
        return 1;
641
      }
642
      sctx->nspanHist->axis[0].min = sctx->range->min;
643
      sctx->nspanHist->axis[1].min = sctx->range->min;
644
      sctx->nspanHist->axis[0].max = sctx->range->max;
645
      sctx->nspanHist->axis[1].max = sctx->range->max;
646
647
      if (sctx->ninscl) {
648
        lup = nrrdDLookup[sctx->ninscl->type];
649
        data = sctx->ninscl->data;
650
      } else {
651
        lup = nrrdDLookup[sctx->nsclDerived->type];
652
        data = sctx->nsclDerived->data;
653
      }
654
655
      /* calculate the span space histogram */
656
      if (nrrdMaybeAlloc_va(sctx->nspanHist, nrrdTypeUInt, 2,
657
                            AIR_CAST(size_t, sctx->spanSize),
658
                            AIR_CAST(size_t, sctx->spanSize))) {
659
        biffMovef(SEEK, NRRD,
660
                  "%s: couldn't allocate space space histogram", me);
661
        return 1;
662
      }
663
      spanHist = AIR_CAST(unsigned int*, sctx->nspanHist->data);
664
      sx = sctx->sx;
665
      sy = sctx->sy;
666
      sz = sctx->sz;
667
      ss = sctx->spanSize;
668
      for (si=0; si<ss*ss; si++) {
669
        spanHist[si] = 0;
670
      }
671
      for (zi=0; zi<sz-1; zi++) {
672
        for (yi=0; yi<sy-1; yi++) {
673
          for (xi=0; xi<sx-1; xi++) {
674
            vi = xi + sx*(yi + sy*zi);
675
            val = lup(data, vi + 0 + 0*sx + 0*sx*sy);
676
            min = max = val;
677
            val = lup(data, vi + 1 + 0*sx + 0*sx*sy);
678
            min = AIR_MIN(min, val);
679
            max = AIR_MAX(max, val);
680
            val = lup(data, vi + 0 + 1*sx + 0*sx*sy);
681
            min = AIR_MIN(min, val);
682
            max = AIR_MAX(max, val);
683
            val = lup(data, vi + 1 + 1*sx + 0*sx*sy);
684
            min = AIR_MIN(min, val);
685
            max = AIR_MAX(max, val);
686
            val = lup(data, vi + 0 + 0*sx + 1*sx*sy);
687
            min = AIR_MIN(min, val);
688
            max = AIR_MAX(max, val);
689
            val = lup(data, vi + 1 + 0*sx + 1*sx*sy);
690
            min = AIR_MIN(min, val);
691
            max = AIR_MAX(max, val);
692
            val = lup(data, vi + 0 + 1*sx + 1*sx*sy);
693
            min = AIR_MIN(min, val);
694
            max = AIR_MAX(max, val);
695
            val = lup(data, vi + 1 + 1*sx + 1*sx*sy);
696
            min = AIR_MIN(min, val);
697
            max = AIR_MAX(max, val);
698
            minI = airIndex(sctx->range->min, min, sctx->range->max, ss);
699
            maxI = airIndex(sctx->range->min, max, sctx->range->max, ss);
700
            spanHist[minI + ss*maxI]++;
701
          }
702
        }
703
      }
704
    }
705
    sctx->flag[flagSclDerived] = AIR_FALSE;
706
    sctx->flag[flagSpanSpaceHist] = AIR_TRUE;
707
  }
708
  return 0;
709
}
710
711
static int
712
updateResult(seekContext *sctx) {
713
  static const char me[]="updateResult";
714
715
  if (sctx->verbose > 5) {
716
    fprintf(stderr, "%s: --------------------\n", me);
717
    fprintf(stderr, "%s: flagIsovalue = %d\n", me,
718
            sctx->flag[flagIsovalue]);
719
    fprintf(stderr, "%s: flagAnswerPointers = %d\n", me,
720
            sctx->flag[flagAnswerPointers]);
721
    fprintf(stderr, "%s: flagType = %d\n", me,
722
            sctx->flag[flagType]);
723
    fprintf(stderr, "%s: flagSlabCacheAlloc = %d\n", me,
724
            sctx->flag[flagSlabCacheAlloc]);
725
    fprintf(stderr, "%s: flagSpanSpaceHist = %d\n", me,
726
            sctx->flag[flagSpanSpaceHist]);
727
    fprintf(stderr, "%s: flagNinEtAl = %d\n", me,
728
            sctx->flag[flagNinEtAl]);
729
    fprintf(stderr, "%s: flagReverse = %d\n", me,
730
            sctx->flag[flagReverse]);
731
    fprintf(stderr, "%s: flagTxfNormal = %d\n", me,
732
            sctx->flag[flagTxfNormal]);
733
  }
734
735
  if (seekTypeIsocontour != sctx->type
736
      && sctx->flag[flagIsovalue]) {
737
    biffAddf(SEEK, "%s: can't set isovalue for %s (only %s)", me,
738
             airEnumStr(seekType, sctx->type),
739
             airEnumStr(seekType, seekTypeIsocontour));
740
    return 1;
741
  }
742
743
  if (sctx->strengthUse && !sctx->stngAns) {
744
    biffAddf(SEEK, "%s: can't use feature strength without a strength item",
745
             me);
746
    return 1;
747
  }
748
749
  /* this seems to be a very pointless exercise */
750
  if (sctx->flag[flagIsovalue]
751
      || sctx->flag[flagEvalDiffThresh]
752
      || sctx->flag[flagAnswerPointers]
753
      || sctx->flag[flagStrengthUse]
754
      || sctx->flag[flagStrength]
755
      || sctx->flag[flagType]
756
      || sctx->flag[flagSlabCacheAlloc]
757
      || sctx->flag[flagSpanSpaceHist]
758
      || sctx->flag[flagNinEtAl]
759
      || sctx->flag[flagReverse]
760
      || sctx->flag[flagTxfNormal]) {
761
762
    sctx->flag[flagResult] = AIR_TRUE;
763
764
    sctx->flag[flagIsovalue] = AIR_FALSE;
765
    sctx->flag[flagEvalDiffThresh] = AIR_FALSE;
766
    sctx->flag[flagAnswerPointers] = AIR_FALSE;
767
    sctx->flag[flagStrengthUse] = AIR_FALSE;
768
    sctx->flag[flagStrength] = AIR_FALSE;
769
    sctx->flag[flagType] = AIR_FALSE;
770
    sctx->flag[flagSlabCacheAlloc] = AIR_FALSE;
771
    sctx->flag[flagSpanSpaceHist] = AIR_FALSE;
772
    sctx->flag[flagNinEtAl] = AIR_FALSE;
773
    sctx->flag[flagReverse] = AIR_FALSE;
774
    sctx->flag[flagTxfNormal] = AIR_FALSE;
775
  }
776
777
  return 0;
778
}
779
780
/*
781
******** seekUpdate
782
**
783
** traverses dependency graph for all of seek stuff, which necessarily
784
** includes nodes specific to, say, isosurfacing, even when isosurfacing
785
** is not being done.  The idea is to use the same graph for all feature
786
** types, for error checking if nothing else, leavings steps as no-ops
787
** as needed.
788
*/
789
int
790
seekUpdate(seekContext *sctx) {
791
  static const char me[]="seekUpdate";
792
  int E;
793
794
  if (!sctx) {
795
    biffAddf(SEEK, "%s: got NULL pointer", me);
796
    return 1;
797
  }
798
  E = 0;
799
  if (!E) E |= updateNinEtAl(sctx);
800
  if (!E) E |= updateAnswerPointers(sctx);
801
  if (!E) E |= updateSxSySz(sctx);
802
  if (!E) E |= updateReverse(sctx);
803
  if (!E) E |= updateTxfNormal(sctx);
804
  if (!E) E |= updateSlabCacheAlloc(sctx);
805
  if (!E) E |= updateSclDerived(sctx);
806
  if (!E) E |= updateSpanSpaceHist(sctx);
807
  if (!E) E |= updateResult(sctx);
808
  if (E) {
809
    biffAddf(SEEK, "%s: trouble updating", me);
810
    return 1;
811
  }
812
813
  return 0;
814
}