GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: src/mite/txf.c Lines: 0 315 0.0 %
Date: 2017-05-26 Branches: 0 203 0.0 %

Line Branch Exec Source
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
/* learned: don't confuse allocate an array of structs with an array
25
   of pointers to structs.  Don't be surprised when you bus error
26
   because of the difference
27
*/
28
29
#include "mite.h"
30
#include "privateMite.h"
31
32
char
33
miteRangeChar[MITE_RANGE_NUM+1] = "ARGBEadsp";
34
35
const char *
36
_miteStageOpStr[] = {
37
  "(unknown miteStageOp)",
38
  "min",
39
  "max",
40
  "add",
41
  "multiply"
42
};
43
44
const int
45
_miteStageOpVal[] = {
46
  miteStageOpUnknown,
47
  miteStageOpMin,
48
  miteStageOpMax,
49
  miteStageOpAdd,
50
  miteStageOpMultiply
51
};
52
53
const char *
54
_miteStageOpStrEqv[] = {
55
  "min",
56
  "max",
57
  "add", "+",
58
  "multiply", "*", "x",
59
  ""
60
};
61
62
const int
63
_miteStageOpValEqv[] = {
64
  miteStageOpMin,
65
  miteStageOpMax,
66
  miteStageOpAdd, miteStageOpAdd,
67
  miteStageOpMultiply, miteStageOpMultiply, miteStageOpMultiply
68
};
69
70
const airEnum
71
_miteStageOp = {
72
  "miteStageOp",
73
  MITE_STAGE_OP_MAX,
74
  _miteStageOpStr, _miteStageOpVal,
75
  NULL,
76
  _miteStageOpStrEqv, _miteStageOpValEqv,
77
  AIR_FALSE
78
};
79
const airEnum *const
80
miteStageOp = &_miteStageOp;
81
82
/*
83
******** miteVariableParse()
84
**
85
** takes a string (usually the label from a nrrd axis) and parses it
86
** to determine the gageItemSpec from it (which means finding the
87
** kind and item).  The valid formats are:
88
**
89
**   ""                  : NULL kind, 0 item
90
**   <item>              : miteValGageKind (DEPRECATED)
91
**   mite(<item>)        : miteValGageKind
92
**   gage(<item>)        : gageKindScl (DEPRECATED)
93
**   gage(scalar:<item>) : gageKindScl
94
**   gage(vector:<item>) : gageKindVec
95
**   gage(tensor:<item>) : tenGageKind
96
**
97
** Notice that "scalar", "vector", and "tensor" do NOT refer to the type
98
** of the quantity being measured, but rather to the type of volume in
99
** which quantity is measured (i.e., the gageKind used)
100
*/
101
int
102
miteVariableParse(gageItemSpec *isp, const char *label) {
103
  static const char me[]="miteVariableParse";
104
105
  char *buff, *endparen, *kqstr, *col, *kstr, *qstr;
106
  airArray *mop;
107
108
  if (!( isp && label )) {
109
    biffAddf(MITE, "%s: got NULL pointer", me);
110
    return 1;
111
  }
112
  if (0 == strlen(label)) {
113
    /* nothing was specified; we try to indicate that by mimicking
114
       the return of gageItemSpecNew() */
115
    isp->item = 0;
116
    isp->kind = NULL;
117
    return 0;
118
  }
119
  /* else given string was non-empty */
120
  mop = airMopNew();
121
  buff = airStrdup(label);
122
  if (!buff) {
123
    biffAddf(MITE, "%s: couldn't strdup label!", me);
124
    airMopError(mop); return 1;
125
  }
126
  airMopAdd(mop, buff, airFree, airMopAlways);
127
  if (strstr(buff, "gage(") == buff) {
128
    /* txf domain variable is to be measured directly by gage */
129
    if (!(endparen = strstr(buff, ")"))) {
130
      biffAddf(MITE, "%s: didn't see close paren after \"gage(\"", me);
131
      airMopError(mop); return 1;
132
    }
133
    *endparen = 0;
134
    kqstr = buff + strlen("gage(");
135
    /* first see if its a (deprecated) gageKindScl specification */
136
    isp->item = airEnumVal(gageScl, kqstr);
137
    if (0 != isp->item) {
138
      isp->kind = gageKindScl;
139
      fprintf(stderr, "\n%s: WARNING: deprecated use of txf domain "
140
              "\"gage(%s)\" without explicit gage kind specification; "
141
              "should use \"gage(%s:%s)\" instead\n\n",
142
              me, kqstr, gageKindScl->name, kqstr);
143
    } else {
144
      /* should be of form "<kind>:<item>" */
145
      col = strstr(kqstr, ":");
146
      if (!col) {
147
        biffAddf(MITE, "%s: didn't see \":\" separator between gage "
148
                 "kind and item", me);
149
        airMopError(mop); return 1;
150
      }
151
      *col = 0;
152
      kstr = kqstr;
153
      qstr = col+1;
154
      if (!strcmp(gageKindScl->name, kstr)) {
155
        isp->kind = gageKindScl;
156
      } else if (!strcmp(gageKindVec->name, kstr)) {
157
        isp->kind = gageKindVec;
158
      } else if (!strcmp(tenGageKind->name, kstr)) {
159
        isp->kind = tenGageKind;
160
      } else {
161
        biffAddf(MITE, "%s: don't recognized \"%s\" gage kind", me, kstr);
162
        airMopError(mop); return 1;
163
      }
164
      isp->item = airEnumVal(isp->kind->enm, qstr);
165
      if (0 == isp->item) {
166
        biffAddf(MITE, "%s: couldn't parse \"%s\" as a %s variable",
167
                 me, qstr, isp->kind->name);
168
        airMopError(mop); return 1;
169
      }
170
    }
171
  } else if (strstr(buff, "mite(") == buff) {
172
    /* txf domain variable is *not* directly measured by gage */
173
    if (!(endparen = strstr(buff, ")"))) {
174
      biffAddf(MITE, "%s: didn't see close paren after \"mite(\"", me);
175
      airMopError(mop); return 1;
176
    }
177
    *endparen = 0;
178
    qstr = buff + strlen("mite(");
179
    isp->item = airEnumVal(miteVal, qstr);
180
    if (0 == isp->item) {
181
      biffAddf(MITE, "%s: couldn't parse \"%s\" as a miteVal variable",
182
               me, qstr);
183
      airMopError(mop); return 1;
184
    }
185
    isp->kind = miteValGageKind;
186
  } else {
187
    /* didn't start with "gage(" or "mite(" */
188
    isp->item = airEnumVal(miteVal, label);
189
    if (0 != isp->item) {
190
      /* its measured by mite */
191
      isp->kind = miteValGageKind;
192
      fprintf(stderr, "\n%s: WARNING: deprecated use of txf domain "
193
              "\"%s\"; should use \"mite(%s)\" instead\n\n",
194
              me, label, label);
195
    } else {
196
      biffAddf(MITE, "%s: \"%s\" not a recognized variable", me, label);
197
      airMopError(mop); return 1;
198
    }
199
  }
200
  airMopOkay(mop);
201
  return 0;
202
}
203
204
void
205
miteVariablePrint(char *buff, const gageItemSpec *isp) {
206
  static const char me[]="miteVariablePrint";
207
208
  if (!(isp->kind)) {
209
    strcpy(buff, "");
210
  } else if (gageKindScl == isp->kind
211
             || gageKindVec == isp->kind
212
             || tenGageKind == isp->kind) {
213
    sprintf(buff, "gage(%s:%s)", isp->kind->name,
214
            airEnumStr(isp->kind->enm, isp->item));
215
  } else if (miteValGageKind == isp->kind) {
216
    sprintf(buff, "%s(%s)", isp->kind->name,
217
            airEnumStr(isp->kind->enm, isp->item));
218
  } else {
219
    sprintf(buff, "(%s: unknown gageKind!)", me);
220
  }
221
  return;
222
}
223
224
int
225
miteNtxfCheck(const Nrrd *ntxf) {
226
  static const char me[]="miteNtxfCheck";
227
  char *rangeStr, *domStr;
228
  gageItemSpec isp;
229
  unsigned int rii, axi;
230
  int ilog2;
231
232
  if (nrrdCheck(ntxf)) {
233
    biffMovef(MITE, NRRD, "%s: basic nrrd validity check failed", me);
234
    return 1;
235
  }
236
  if (!( nrrdTypeFloat == ntxf->type ||
237
         nrrdTypeDouble == ntxf->type ||
238
         nrrdTypeUChar == ntxf->type )) {
239
    biffAddf(MITE, "%s: need a type %s, %s or %s nrrd (not %s)", me,
240
             airEnumStr(nrrdType, nrrdTypeFloat),
241
             airEnumStr(nrrdType, nrrdTypeDouble),
242
             airEnumStr(nrrdType, nrrdTypeUChar),
243
             airEnumStr(nrrdType, ntxf->type));
244
    return 1;
245
  }
246
  if (!( 2 <= ntxf->dim )) {
247
    biffAddf(MITE, "%s: nrrd dim (%d) isn't at least 2 (for a 1-D txf)",
248
             me, ntxf->dim);
249
    return 1;
250
  }
251
  rangeStr = ntxf->axis[0].label;
252
  if (0 == airStrlen(rangeStr)) {
253
    biffAddf(MITE, "%s: axis[0]'s label doesn't specify txf range", me);
254
    return 1;
255
  }
256
  if (airStrlen(rangeStr) != ntxf->axis[0].size) {
257
    char stmp1[AIR_STRLEN_SMALL], stmp2[AIR_STRLEN_SMALL];
258
    biffAddf(MITE, "%s: axis[0]'s size %s, but label specifies %s values", me,
259
             airSprintSize_t(stmp1, ntxf->axis[0].size),
260
             airSprintSize_t(stmp2, airStrlen(rangeStr)));
261
    return 1;
262
  }
263
  for (rii=0; rii<airStrlen(rangeStr); rii++) {
264
    if (!strchr(miteRangeChar, rangeStr[rii])) {
265
      biffAddf(MITE, "%s: char %d of axis[0]'s label (\"%c\") isn't a valid "
266
               "transfer function range specifier (not in \"%s\")",
267
               me, rii, rangeStr[rii], miteRangeChar);
268
      return 1;
269
    }
270
  }
271
  for (axi=1; axi<ntxf->dim; axi++) {
272
    if (1 == ntxf->axis[axi].size) {
273
      biffAddf(MITE, "%s: # samples on axis %d must be > 1", me, axi);
274
      return 1;
275
    }
276
    domStr = ntxf->axis[axi].label;
277
    if (0 == airStrlen(domStr)) {
278
      biffAddf(MITE, "%s: axis[%d] of txf didn't specify a domain variable",
279
               me, axi);
280
      return 1;
281
    }
282
    if (miteVariableParse(&isp, domStr)) {
283
      biffAddf(MITE, "%s: couldn't parse txf domain \"%s\" for axis %d\n",
284
               me, domStr, axi);
285
      return 1;
286
    }
287
    if (!( 1 == isp.kind->table[isp.item].answerLength ||
288
           3 == isp.kind->table[isp.item].answerLength )) {
289
      biffAddf(MITE, "%s: %s (item %d) not a scalar or vector "
290
               "(answerLength = %d): "
291
               "can't be a txf domain variable", me, domStr, isp.item,
292
               isp.kind->table[isp.item].answerLength);
293
      return 1;
294
    }
295
    if (3 == isp.kind->table[isp.item].answerLength) {
296
      /* has to be right length for one of the quantization schemes */
297
      ilog2 = airLog2(ntxf->axis[axi].size);
298
      if (-1 == ilog2) {
299
        char stmp[AIR_STRLEN_SMALL];
300
        biffAddf(MITE, "%s: txf axis size for %s must be power of 2 (not %s)",
301
                 me, domStr, airSprintSize_t(stmp, ntxf->axis[axi].size));
302
        return 1;
303
      } else {
304
        if (!( AIR_IN_CL(8, ilog2, 16) )) {
305
          biffAddf(MITE, "%s: log_2 of txf axis size for %s should be in "
306
                   "range [8,16] (not %d)", me, domStr, ilog2);
307
          return 1;
308
        }
309
      }
310
    } else {
311
      if (!( AIR_EXISTS(ntxf->axis[axi].min) &&
312
             AIR_EXISTS(ntxf->axis[axi].max) )) {
313
        biffAddf(MITE, "%s: min and max of axis %d aren't both set", me, axi);
314
        return 1;
315
      }
316
      if (!( ntxf->axis[axi].min < ntxf->axis[axi].max )) {
317
        biffAddf(MITE, "%s: min (%g) not less than max (%g) on axis %d",
318
                 me, ntxf->axis[axi].min, ntxf->axis[axi].max, axi);
319
        return 1;
320
      }
321
    }
322
  }
323
324
  return 0;
325
}
326
327
/*
328
******** miteQueryAdd()
329
**
330
** This looks a given gageItemSpec and sets the bits in the
331
** gageKindScl and tenGageKind queries that are required to calculate
332
** the quantity
333
**
334
** NOTE: This does NOT initialize the query{Scl,Vec,Ten}: it
335
** just adds on new items to the existing queries
336
**
337
** HEY: this is really unfortunate: each new gage type use for
338
** volume rendering in mite will have to explicitly added as
339
** arguments here, which is part of the reason that mite may end
340
** up explicitly depending on the libraries supplying those gageKinds
341
** (like how mite now must depend on ten)
342
**
343
** The queryMite argument is a little odd- its not a real gage kind,
344
** but we use it to organize which of the miteVal quantities we take
345
** the time to compute in miteSample().
346
*/
347
void
348
miteQueryAdd(gageQuery queryScl, gageQuery queryVec,
349
             gageQuery queryTen, gageQuery queryMite,
350
             gageItemSpec *isp) {
351
  static const char me[]="miteQueryAdd";
352
353
  if (NULL == isp->kind) {
354
    /* nothing to add */
355
  } else if (gageKindScl == isp->kind) {
356
    GAGE_QUERY_ITEM_ON(queryScl, isp->item);
357
  } else if (gageKindVec == isp->kind) {
358
    GAGE_QUERY_ITEM_ON(queryVec, isp->item);
359
  } else if (tenGageKind == isp->kind) {
360
    GAGE_QUERY_ITEM_ON(queryTen, isp->item);
361
  } else if (miteValGageKind == isp->kind) {
362
    /* regardless of whether the mite query requires scl, vec, or ten
363
       queries, we add it to the quantites that have to be computed
364
       per-sample */
365
    GAGE_QUERY_ITEM_ON(queryMite, isp->item);
366
    /* HEY: some these have useful analogs for tensor data, but I
367
       won't be able to express them.  This means that while Phong
368
       shading of *scalar* volumes can be implemented with transfer
369
       functions, it is currently not possible in *tensor* volumes
370
       (for instance, using the gradient of fractional anisotropy) */
371
    switch(isp->item) {
372
    case miteValVrefN:
373
    case miteValNdotV:
374
    case miteValNdotL:
375
      /* the "N" can be a normalized vector from any of the
376
         available kinds (except miteValGageKind!), and its
377
         associated query will be handled elsewhere, so there's
378
         nothing to add here */
379
      break;
380
    case miteValGTdotV:
381
      GAGE_QUERY_ITEM_ON(queryScl, gageSclGeomTens);
382
      break;
383
    case miteValVdefT:
384
      GAGE_QUERY_ITEM_ON(queryTen, tenGageTensor);
385
    case miteValVdefTdotV:
386
      GAGE_QUERY_ITEM_ON(queryTen, tenGageTensor);
387
      break;
388
    }
389
  } else {
390
    fprintf(stderr, "%s: PANIC: unrecognized non-NULL gageKind\n", me);
391
    exit(1);
392
  }
393
  return;
394
}
395
396
int
397
_miteNtxfCopy(miteRender *mrr, miteUser *muu) {
398
  static const char me[]="_miteNtxfCopy";
399
  int ni, E;
400
401
  mrr->ntxf = AIR_CALLOC(muu->ntxfNum, Nrrd *);
402
  if (!mrr->ntxf) {
403
    biffAddf(MITE, "%s: couldn't calloc %d ntxf pointers", me, muu->ntxfNum);
404
    return 1;
405
  }
406
  mrr->ntxfNum = muu->ntxfNum;
407
  airMopAdd(mrr->rmop, mrr->ntxf, airFree, airMopAlways);
408
  E = 0;
409
  for (ni=0; ni<mrr->ntxfNum; ni++) {
410
    mrr->ntxf[ni] = nrrdNew();
411
    if (!E) airMopAdd(mrr->rmop, mrr->ntxf[ni],
412
                      (airMopper)nrrdNuke, airMopAlways);
413
    if (!( nrrdTypeUChar == muu->ntxf[ni]->type
414
           || nrrdTypeFloat == muu->ntxf[ni]->type
415
           || nrrdTypeDouble == muu->ntxf[ni]->type )) {
416
      biffAddf(MITE,
417
               "%s: sorry, can't handle txf of type %s (only %s, %s, %s)",
418
               me, airEnumStr(nrrdType, muu->ntxf[ni]->type),
419
               airEnumStr(nrrdType, nrrdTypeUChar),
420
               airEnumStr(nrrdType, nrrdTypeFloat),
421
               airEnumStr(nrrdType, nrrdTypeDouble));
422
      return 1;
423
    }
424
    /* note that key/values need to be copied for the sake of
425
       identifying a non-default miteStageOp */
426
    switch(muu->ntxf[ni]->type) {
427
    case nrrdTypeUChar:
428
      if (!E) E |= nrrdUnquantize(mrr->ntxf[ni], muu->ntxf[ni], nrrdTypeUChar);
429
      if (!E) E |= nrrdKeyValueCopy(mrr->ntxf[ni], muu->ntxf[ni]);
430
      break;
431
    case mite_nt:
432
      if (!E) E |= nrrdCopy(mrr->ntxf[ni], muu->ntxf[ni]);
433
      break;
434
    default:  /* will be either float or double (whatever mite_nt isn't) */
435
      if (!E) E |= nrrdConvert(mrr->ntxf[ni], muu->ntxf[ni], mite_nt);
436
      if (!E) E |= nrrdKeyValueCopy(mrr->ntxf[ni], muu->ntxf[ni]);
437
      break;
438
    }
439
  }
440
  if (E) {
441
    biffMovef(MITE, NRRD, "%s: troubling copying/converting all ntxfs", me);
442
    return 1;
443
  }
444
  return 0;
445
}
446
447
int
448
_miteNtxfAlphaAdjust(miteRender *mrr, miteUser *muu) {
449
  static const char me[]="_miteNtxfAlphaAdjust";
450
  int ni, ei, ri, nnum, rnum;
451
  Nrrd *ntxf;
452
  mite_t *data, alpha, frac;
453
454
  if (_miteNtxfCopy(mrr, muu)) {
455
    biffAddf(MITE, "%s: trouble copying/converting transfer functions", me);
456
    return 1;
457
  }
458
  frac = muu->rayStep/muu->refStep;
459
  for (ni=0; ni<mrr->ntxfNum; ni++) {
460
    ntxf = mrr->ntxf[ni];
461
    if (!strchr(ntxf->axis[0].label, miteRangeChar[miteRangeAlpha])) {
462
      continue;
463
    }
464
    /* else this txf sets opacity */
465
    data = (mite_t *)ntxf->data;
466
    rnum = ntxf->axis[0].size;
467
    nnum = nrrdElementNumber(ntxf)/rnum;
468
    for (ei=0; ei<nnum; ei++) {
469
      for (ri=0; ri<rnum; ri++) {
470
        if (ntxf->axis[0].label[ri] == miteRangeChar[miteRangeAlpha]) {
471
          alpha = data[ri + rnum*ei];
472
          data[ri + rnum*ei] = 1 - pow(AIR_MAX(0, 1-alpha), frac);
473
        }
474
      }
475
    }
476
  }
477
  return 0;
478
}
479
480
int
481
_miteStageNum(miteRender *mrr) {
482
  int num, ni;
483
484
  num = 0;
485
  for (ni=0; ni<mrr->ntxfNum; ni++) {
486
    num += mrr->ntxf[ni]->dim - 1;
487
  }
488
  return num;
489
}
490
491
void
492
_miteStageInit(miteStage *stage) {
493
  int rii;
494
495
  stage->val = NULL;
496
  stage->size = -1;
497
  stage->op = miteStageOpUnknown;
498
  stage->qn = NULL;
499
  stage->min = stage->max = AIR_NAN;
500
  stage->data = NULL;
501
  for (rii=0; rii<=MITE_RANGE_NUM-1; rii++) {
502
    stage->rangeIdx[rii] = -1;
503
  }
504
  stage->rangeNum = -1;
505
  stage->label = NULL;
506
  return;
507
}
508
509
double *
510
_miteAnswerPointer(miteThread *mtt, gageItemSpec *isp) {
511
  static const char me[]="_miteAnswerPointer";
512
  double *ret;
513
514
  if (!isp->kind) {
515
    /* we got a NULL kind (as happens with output of
516
       gageItemSpecNew(), or miteVariableParse of an
517
       empty string); only NULL return is sensible */
518
    return NULL;
519
  }
520
521
  if (gageKindScl == isp->kind) {
522
    ret = mtt->ansScl;
523
  } else if (gageKindVec == isp->kind) {
524
    ret = mtt->ansVec;
525
  } else if (tenGageKind == isp->kind) {
526
    ret = mtt->ansTen;
527
  } else if (miteValGageKind == isp->kind) {
528
    ret = mtt->ansMiteVal;
529
  } else {
530
    fprintf(stderr, "\nPANIC: %s: unknown gageKind!\n", me);
531
    exit(1);
532
  }
533
  ret += gageKindAnswerOffset(isp->kind, isp->item);
534
  return ret;
535
}
536
537
/*
538
** _miteStageSet
539
**
540
** ALLOCATES and initializes stage array in a miteThread
541
*/
542
int
543
_miteStageSet(miteThread *mtt, miteRender *mrr) {
544
  static const char me[]="_miteStageSet";
545
  char *value;
546
  int ni, di, stageIdx, rii, stageNum, ilog2;
547
  Nrrd *ntxf;
548
  miteStage *stage;
549
  gageItemSpec isp;
550
  char rc;
551
552
  stageNum = _miteStageNum(mrr);
553
  /* fprintf(stderr, "!%s: stageNum = %d\n", me, stageNum); */
554
  mtt->stage = AIR_CALLOC(stageNum, miteStage);
555
  if (!mtt->stage) {
556
    biffAddf(MITE, "%s: couldn't alloc array of %d stages", me, stageNum);
557
    return 1;
558
  }
559
  airMopAdd(mtt->rmop, mtt->stage, airFree, airMopAlways);
560
  mtt->stageNum = stageNum;
561
  stageIdx = 0;
562
  for (ni=0; ni<mrr->ntxfNum; ni++) {
563
    ntxf = mrr->ntxf[ni];
564
    for (di=ntxf->dim-1; di>=1; di--) {
565
      stage = mtt->stage + stageIdx;
566
      _miteStageInit(stage);
567
      miteVariableParse(&isp, ntxf->axis[di].label);
568
      stage->val = _miteAnswerPointer(mtt, &isp);
569
      stage->label = ntxf->axis[di].label;
570
      /*
571
      fprintf(stderr, "!%s: ans=%p + offset[%d]=%d == %p\n", me,
572
              mtt->ans, dom, kind->ansOffset[dom], stage->val);
573
      */
574
      stage->size = ntxf->axis[di].size;
575
      stage->min =  ntxf->axis[di].min;
576
      stage->max =  ntxf->axis[di].max;
577
      if (di > 1) {
578
        stage->data = NULL;
579
      } else {
580
        stage->data = (mite_t *)ntxf->data;
581
        value = nrrdKeyValueGet(ntxf, "miteStageOp");
582
        if (value) {
583
          stage->op = airEnumVal(miteStageOp, value);
584
          if (miteStageOpUnknown == stage->op) {
585
            stage->op = miteStageOpMultiply;
586
          }
587
        } else {
588
          stage->op = miteStageOpMultiply;
589
        }
590
        if (1 == isp.kind->table[isp.item].answerLength) {
591
          stage->qn = NULL;
592
        } else if (3 == isp.kind->table[isp.item].answerLength) {
593
          char stmp[AIR_STRLEN_SMALL];
594
          ilog2 = airLog2(ntxf->axis[di].size);
595
          switch(ilog2) {
596
          case 8:  stage->qn = limnVtoQN_d[ limnQN8octa]; break;
597
          case 9:  stage->qn = limnVtoQN_d[ limnQN9octa]; break;
598
          case 10: stage->qn = limnVtoQN_d[limnQN10octa]; break;
599
          case 11: stage->qn = limnVtoQN_d[limnQN11octa]; break;
600
          case 12: stage->qn = limnVtoQN_d[limnQN12octa]; break;
601
          case 13: stage->qn = limnVtoQN_d[limnQN13octa]; break;
602
          case 14: stage->qn = limnVtoQN_d[limnQN14octa]; break;
603
          case 15: stage->qn = limnVtoQN_d[limnQN15octa]; break;
604
          case 16: stage->qn = limnVtoQN_d[limnQN16octa]; break;
605
          default:
606
            biffAddf(MITE, "%s: txf axis %d size %s not usable for "
607
                     "vector txf domain variable %s", me, di,
608
                     airSprintSize_t(stmp, ntxf->axis[di].size),
609
                     ntxf->axis[di].label);
610
            return 1;
611
            break;
612
          }
613
        } else {
614
          biffAddf(MITE, "%s: %s not scalar or vector (len = %d): can't be "
615
                   "a txf domain variable", me,
616
                   ntxf->axis[di].label,
617
                   isp.kind->table[isp.item].answerLength);
618
          return 1;
619
        }
620
        stage->rangeNum = ntxf->axis[0].size;
621
        for (rii=0; rii<stage->rangeNum; rii++) {
622
          rc = ntxf->axis[0].label[rii];
623
          stage->rangeIdx[rii] = strchr(miteRangeChar, rc) - miteRangeChar;
624
          /*
625
          fprintf(stderr, "!%s: range: %c -> %d\n", "_miteStageSet",
626
                  ntxf->axis[0].label[rii], stage->rangeIdx[rii]);
627
          */
628
        }
629
      }
630
      stageIdx++;
631
    }
632
  }
633
  return 0;
634
}
635
636
void
637
_miteStageRun(miteThread *mtt, miteUser *muu) {
638
  static const char me[]="_miteStageRun";
639
  int stageIdx, ri, rii;
640
  unsigned int txfIdx, finalIdx;
641
  miteStage *stage;
642
  mite_t *rangeData;
643
  double *dbg=NULL;
644
645
  finalIdx = 0;
646
  if (mtt->verbose) {
647
    dbg = muu->debug + muu->debugIdx;
648
  }
649
  for (stageIdx=0; stageIdx<mtt->stageNum; stageIdx++) {
650
    stage = &(mtt->stage[stageIdx]);
651
    if (stage->qn) {
652
      /* its a vector-valued txf domain variable */
653
      txfIdx = stage->qn(stage->val);
654
      /* right now, we can't store vector-valued txf domain variables */
655
    } else {
656
      /* its a scalar txf domain variable */
657
      txfIdx = airIndexClamp(stage->min, *(stage->val),
658
                             stage->max, stage->size);
659
      if (mtt->verbose) {
660
        fprintf(stderr, "!%s: %s=%g in [%g,%g]/%u -> %u\n", me,
661
                stage->label, *(stage->val),
662
                stage->min, stage->max, stage->size, txfIdx);
663
        dbg[0 + 2*stageIdx] = *(stage->val);
664
      }
665
    }
666
    finalIdx = stage->size*finalIdx + txfIdx;
667
    if (mtt->verbose) {
668
      dbg[1 + 2*stageIdx] = txfIdx;
669
    }
670
    if (stage->data) {
671
      rangeData = stage->data + stage->rangeNum*finalIdx;
672
      for (rii=0; rii<stage->rangeNum; rii++) {
673
        ri = stage->rangeIdx[rii];
674
        switch(stage->op) {
675
        case miteStageOpMin:
676
          mtt->range[ri] = AIR_MIN(mtt->range[ri], rangeData[rii]);
677
          break;
678
        case miteStageOpMax:
679
          mtt->range[ri] = AIR_MAX(mtt->range[ri], rangeData[rii]);
680
          break;
681
        case miteStageOpAdd:
682
          mtt->range[ri] += rangeData[rii];
683
          break;
684
        case miteStageOpMultiply:
685
        default:
686
          mtt->range[ri] *= rangeData[rii];
687
          break;
688
        }
689
      }
690
      finalIdx = 0;
691
    }
692
  }
693
  return;
694
}