Bug Summary

File:src/nrrd/apply1D.c
Location:line 846, column 27
Description:Assigned value is garbage or undefined

Annotated Source Code

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#include "nrrd.h"
25#include "privateNrrd.h"
26
27/*
28** learned: even when using doubles, because of limited floating point
29** precision, you can get different results between quantizing
30** unrescaled (value directly from nrrd, map domain set to nrrd range,
31** as with early behavior of unu rmap) and rescaled (value from nrrd
32** scaled to fit in existing map domain, as with unu imap -r) value,
33** to the exact same index space.
34*/
35
36/*
37** I won't try to support mapping individual values through a
38** colormap, as with a function evaluation on a single passed value.
39** That will be handled in an upcoming library...
40*/
41
42/*
43** this identifies the different kinds of 1D maps, useful for the
44** functions in this file only
45*/
46enum {
47 kindLut=0,
48 kindRmap=1,
49 kindImap=2
50};
51
52double
53_nrrdApplyDomainMin(const Nrrd *nmap, int ramps, int mapAxis) {
54 double ret;
55
56 AIR_UNUSED(ramps)(void)(ramps);
57 ret = nmap->axis[mapAxis].min;
58 ret = AIR_EXISTS(ret)(((int)(!((ret) - (ret))))) ? ret : 0;
59 return ret;
60}
61
62double
63_nrrdApplyDomainMax(const Nrrd *nmap, int ramps, int mapAxis) {
64 double ret;
65
66 ret = nmap->axis[mapAxis].max;
67 if (!AIR_EXISTS(ret)(((int)(!((ret) - (ret)))))) {
68 ret = AIR_CAST(double, nmap->axis[mapAxis].size)((double)(nmap->axis[mapAxis].size));
69 ret = ramps ? ret-1 : ret;
70 }
71 return ret;
72}
73
74/*
75** _nrrdApply1DSetUp()
76**
77** some error checking and initializing needed for 1D LUTS, regular,
78** and irregular maps. The intent is that if this succeeds, then
79** there is no need for any further error checking.
80**
81** The only thing this function DOES is allocate the output nrrd, and
82** set meta information. The rest is just error checking.
83**
84** The given NrrdRange has to be fleshed out by the caller: it can't
85** be NULL, and both range->min and range->max must exist.
86*/
87int
88_nrrdApply1DSetUp(Nrrd *nout, const Nrrd *nin, const NrrdRange *range,
89 const Nrrd *nmap, int kind, int typeOut,
90 int rescale, int multi) {
91 static const char me[]="_nrrdApply1DSetUp";
92 char *mapcnt;
93 char nounStr[][AIR_STRLEN_SMALL(128+1)]={"lut",
94 "regular map",
95 "irregular map"};
96 char mnounStr[][AIR_STRLEN_SMALL(128+1)]={"multi lut",
97 "multi regular map",
98 "multi irregular map"};
99 /* wishful thinking */
100 char verbStr[][AIR_STRLEN_SMALL(128+1)]={"lut",
101 "rmap",
102 "imap"};
103 char mverbStr[][AIR_STRLEN_SMALL(128+1)]={"mlut",
104 "mrmap",
105 "mimap"}; /* wishful thinking */
106 int mapAxis, copyMapAxis0=AIR_FALSE0, axisMap[NRRD_DIM_MAX16];
107 unsigned int ax, dim, entLen;
108 size_t size[NRRD_DIM_MAX16];
109 double domMin, domMax;
110
111 if (nout == nin) {
112 biffAddf(NRRDnrrdBiffKey, "%s: due to laziness, nout==nin always disallowed", me);
113 return 1;
114 }
115 if (airEnumValCheck(nrrdType, typeOut)) {
116 biffAddf(NRRDnrrdBiffKey, "%s: invalid requested output type %d", me, typeOut);
117 return 1;
118 }
119 if (nrrdTypeBlock == nin->type || nrrdTypeBlock == typeOut) {
120 biffAddf(NRRDnrrdBiffKey, "%s: input or requested output type is %s, need scalar",
121 me, airEnumStr(nrrdType, nrrdTypeBlock));
122 return 1;
123 }
124 if (rescale) {
125 if (!range) {
126 biffAddf(NRRDnrrdBiffKey, "%s: want rescaling but didn't get a range", me);
127 return 1;
128 }
129 if (!( AIR_EXISTS(range->min)(((int)(!((range->min) - (range->min))))) && AIR_EXISTS(range->max)(((int)(!((range->max) - (range->max))))) )) {
130 biffAddf(NRRDnrrdBiffKey, "%s: want rescaling but not both "
131 "range->{min,max} %g %g exist", me, range->min, range->max);
132 return 1;
133 }
134 /* HEY: consider adding an error check for range->min == range->max
135 here; the code below now guards
136 AIR_AFFINE(range->min, val, range->max, domMin, domMax)
137 against producing NaNs, but maybe that's too clever */
138 }
139 if (kindLut == kind || kindRmap == kind) {
140 if (!multi) {
141 mapAxis = nmap->dim - 1;
142 if (!(0 == mapAxis || 1 == mapAxis)) {
143 biffAddf(NRRDnrrdBiffKey, "%s: dimension of %s should be 1 or 2, not %d",
144 me, nounStr[kind], nmap->dim);
145 return 1;
146 }
147 copyMapAxis0 = (1 == mapAxis);
148 } else {
149 mapAxis = nmap->dim - nin->dim - 1;
150 if (!(0 == mapAxis || 1 == mapAxis)) {
151 biffAddf(NRRDnrrdBiffKey, "%s: dimension of %s should be %d or %d, not %d",
152 me, mnounStr[kind],
153 nin->dim + 1, nin->dim + 2, nmap->dim);
154 return 1;
155 }
156 copyMapAxis0 = (1 == mapAxis);
157 /* need to make sure the relevant sizes match */
158 for (ax=0; ax<nin->dim; ax++) {
159 unsigned int taxi = mapAxis + 1 + ax;
160 if (taxi > NRRD_DIM_MAX16-1) {
161 biffAddf(NRRDnrrdBiffKey, "%s: test axis index %u exceeds NRRD_DIM_MAX-1 %u",
162 me, taxi, NRRD_DIM_MAX16-1);
163 return 1;
164 }
165 if (nin->axis[ax].size != nmap->axis[taxi].size) {
166 char stmp1[AIR_STRLEN_SMALL(128+1)], stmp2[AIR_STRLEN_SMALL(128+1)];
167 biffAddf(NRRDnrrdBiffKey, "%s: input and mmap don't have compatible sizes: "
168 "nin->axis[%d].size (%s) "
169 "!= nmap->axis[%d].size (%s): ", me,
170 ax, airSprintSize_t(stmp1, nin->axis[ax].size),
171 mapAxis + 1 + ax,
172 airSprintSize_t(stmp2, nmap->axis[taxi].size));
173 return 1;
174 }
175 }
176 }
177 domMin = _nrrdApplyDomainMin(nmap, AIR_FALSE0, mapAxis);
178 domMax = _nrrdApplyDomainMax(nmap, AIR_FALSE0, mapAxis);
179 if (!( domMin < domMax )) {
180 biffAddf(NRRDnrrdBiffKey, "%s: (axis %d) domain min (%g) not less than max (%g)",
181 me, mapAxis, domMin, domMax);
182 return 1;
183 }
184 if (nrrdHasNonExist(nmap)) {
185 biffAddf(NRRDnrrdBiffKey, "%s: %s nrrd has non-existent values",
186 me, multi ? mnounStr[kind] : nounStr[kind]);
187 return 1;
188 }
189 entLen = mapAxis ? AIR_CAST(unsigned int, nmap->axis[0].size)((unsigned int)(nmap->axis[0].size)) : 1;
190 } else {
191 if (multi) {
192 biffAddf(NRRDnrrdBiffKey, "%s: sorry, multi irregular maps not implemented", me);
193 return 1;
194 }
195 /* its an irregular map */
196 if (nrrd1DIrregMapCheck(nmap)) {
197 biffAddf(NRRDnrrdBiffKey, "%s: problem with irregular map", me);
198 return 1;
199 }
200 /* mapAxis has no meaning for irregular maps, but we'll pretend ... */
201 mapAxis = nmap->axis[0].size == 2 ? 0 : 1;
202 copyMapAxis0 = AIR_TRUE1;
203 entLen = AIR_CAST(unsigned int, nmap->axis[0].size-1)((unsigned int)(nmap->axis[0].size-1));
204 }
205 if (mapAxis + nin->dim > NRRD_DIM_MAX16) {
206 biffAddf(NRRDnrrdBiffKey, "%s: input nrrd dim %d through non-scalar %s exceeds "
207 "NRRD_DIM_MAX %d",
208 me, nin->dim,
209 multi ? mnounStr[kind] : nounStr[kind], NRRD_DIM_MAX16);
210 return 1;
211 }
212 nrrdAxisInfoGet_nva(nin, nrrdAxisInfoSize, size+mapAxis);
213 if (mapAxis) {
214 size[0] = entLen;
215 axisMap[0] = -1;
216 }
217 for (dim=0; dim<nin->dim; dim++) {
218 axisMap[dim+mapAxis] = dim;
219 }
220 /*
221 fprintf(stderr, "##%s: pre maybe alloc: nout->data = %p\n", me, nout->data);
222 for (dim=0; dim<mapAxis + nin->dim; dim++) {
223 fprintf(stderr, " size[%d] = %d\n", d, (int)size[d]);
224 }
225 fprintf(stderr, " nout->dim = %d; nout->type = %d = %s; sizes = %d,%d\n",
226 nout->dim, nout->type,
227 airEnumStr(nrrdType, nout->type));
228 fprintf(stderr, " typeOut = %d = %s\n", typeOut,
229 airEnumStr(nrrdType, typeOut));
230 */
231 if (nrrdMaybeAlloc_nva(nout, typeOut, mapAxis + nin->dim, size)) {
232 biffAddf(NRRDnrrdBiffKey, "%s: couldn't allocate output nrrd", me);
233 return 1;
234 }
235 /*
236 fprintf(stderr, " nout->dim = %d; nout->type = %d = %s\n",
237 nout->dim, nout->type,
238 airEnumStr(nrrdType, nout->type),
239 nout->axis[0].size, nout->axis[1].size);
240 for (d=0; d<nout->dim; d++) {
241 fprintf(stderr, " size[%d] = %d\n", d, (int)nout->axis[d].size);
242 }
243 fprintf(stderr, "##%s: post maybe alloc: nout->data = %p\n", me, nout->data);
244 */
245 if (nrrdAxisInfoCopy(nout, nin, axisMap, NRRD_AXIS_INFO_NONE0)) {
246 biffAddf(NRRDnrrdBiffKey, "%s: trouble copying axis info", me);
247 return 1;
248 }
249 if (copyMapAxis0) {
250 _nrrdAxisInfoCopy(nout->axis + 0, nmap->axis + 0,
251 NRRD_AXIS_INFO_SIZE_BIT(1<< 1));
252 }
253
254 mapcnt = _nrrdContentGet(nmap);
255 if (nrrdContentSet_va(nout, multi ? mverbStr[kind] : verbStr[kind],
256 nin, "%s", mapcnt)) {
257 biffAddf(NRRDnrrdBiffKey, "%s:", me);
258 free(mapcnt); return 1;
259 }
260 free(mapcnt);
261 if (nrrdBasicInfoCopy(nout, nin,
262 NRRD_BASIC_INFO_DATA_BIT(1<< 1)
263 | NRRD_BASIC_INFO_TYPE_BIT(1<< 2)
264 | NRRD_BASIC_INFO_BLOCKSIZE_BIT(1<< 3)
265 | NRRD_BASIC_INFO_DIMENSION_BIT(1<< 4)
266 | NRRD_BASIC_INFO_CONTENT_BIT(1<< 5)
267 | (nrrdStateKeyValuePairsPropagate
268 ? 0
269 : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT(1<<15)))) {
270 biffAddf(NRRDnrrdBiffKey, "%s:", me);
271 return 1;
272 }
273 return 0;
274}
275
276/*
277** _nrrdApply1DLutOrRegMap()
278**
279** the guts of nrrdApply1DLut and nrrdApply1DRegMap
280**
281** yikes, does NOT use biff, since we're only supposed to be called
282** after copious error checking.
283**
284** FOR INSTANCE, this allows nout == nin, which could be a big
285** problem if mapAxis == 1.
286**
287** we don't need a typeOut arg because nout has already been allocated
288** as some specific type; we'll look at that.
289**
290** NOTE: non-existent values get passed through regular maps and luts
291** "unchanged". However, if the output type is integral, the results
292** are probaby undefined. HEY: there is currently no warning message
293** or error handling based on nrrdStateDisallowIntegerNonExist, but
294** there really should be.
295*/
296int
297_nrrdApply1DLutOrRegMap(Nrrd *nout, const Nrrd *nin, const NrrdRange *range,
298 const Nrrd *nmap, int ramps, int rescale, int multi) {
299 /* static const char me[]="_nrrdApply1DLutOrRegMap"; */
300 char *inData, *outData, *mapData, *entData0, *entData1;
301 size_t N, I;
302 double (*inLoad)(const void *v), (*mapLup)(const void *v, size_t I),
303 (*outInsert)(void *v, size_t I, double d),
304 val, mapIdxFrac, domMin, domMax;
305 unsigned int i, mapAxis, mapLen, mapIdx, entSize, entLen, inSize, outSize;
306
307 if (!multi) {
308 mapAxis = nmap->dim - 1; /* axis of nmap containing entries */
309 } else {
310 mapAxis = nmap->dim - nin->dim - 1;
311 }
312 mapData = (char *)nmap->data; /* map data, as char* */
313 /* low end of map domain */
314 domMin = _nrrdApplyDomainMin(nmap, ramps, mapAxis);
315 /* high end of map domain */
316 domMax = _nrrdApplyDomainMax(nmap, ramps, mapAxis);
317 mapLen = AIR_CAST(unsigned int, nmap->axis[mapAxis].size)((unsigned int)(nmap->axis[mapAxis].size)); /* number of entries in map */
318 mapLup = nrrdDLookup[nmap->type]; /* how to get doubles out of map */
319 inData = (char *)nin->data; /* input data, as char* */
320 inLoad = nrrdDLoad[nin->type]; /* how to get doubles out of nin */
321 inSize = AIR_CAST(unsigned int, nrrdElementSize(nin))((unsigned int)(nrrdElementSize(nin))); /* size of one input value */
322 outData = (char *)nout->data; /* output data, as char* */
323 outInsert = nrrdDInsert[nout->type]; /* putting doubles into output */
324 entLen = (mapAxis /* number of elements in one entry */
325 ? AIR_CAST(unsigned int, nmap->axis[0].size)((unsigned int)(nmap->axis[0].size))
326 : 1);
327 outSize = entLen*AIR_CAST(unsigned int, nrrdElementSize(nout))((unsigned int)(nrrdElementSize(nout))); /* size of entry in output */
328 entSize = entLen*AIR_CAST(unsigned int, nrrdElementSize(nmap))((unsigned int)(nrrdElementSize(nmap))); /* size of entry in map */
329
330 N = nrrdElementNumber(nin); /* the number of values to be mapped */
331 if (ramps) {
332 /* regular map */
333 for (I=0; I<N; I++) {
334 /* ...
335 if (!(I % 100)) fprintf(stderr, "I = %d\n", (int)I);
336 ... */
337 val = inLoad(inData);
338 /* ...
339 fprintf(stderr, "##%s: val = \na% 31.15f --> ", me, val);
340 ... */
341 if (rescale) {
342 val = (range->min != range->max
343 ? AIR_AFFINE(range->min, val, range->max, domMin, domMax)( ((double)(domMax)-(domMin))*((double)(val)-(range->min))
/ ((double)(range->max)-(range->min)) + (domMin))
344 : domMin);
345 /* ...
346 fprintf(stderr, "\nb% 31.15f (min,max = %g,%g)--> ", val,
347 range->min, range->max);
348 ... */
349 }
350 /* ...
351 fprintf(stderr, "\nc% 31.15f --> clamp(%g,%g), %d --> ",
352 val, domMin, domMax, mapLen);
353 ... */
354 if (AIR_EXISTS(val)(((int)(!((val) - (val)))))) {
355 val = AIR_CLAMP(domMin, val, domMax)((val) < (domMin) ? (domMin) : ((val) > (domMax) ? (domMax
) : (val)));
356 mapIdxFrac = AIR_AFFINE(domMin, val, domMax, 0, mapLen-1)( ((double)(mapLen-1)-(0))*((double)(val)-(domMin)) / ((double
)(domMax)-(domMin)) + (0));
357 /* ...
358 fprintf(stderr, "mapIdxFrac = \nd% 31.15f --> ",
359 mapIdxFrac);
360 ... */
361 mapIdx = (unsigned int)mapIdxFrac;
362 mapIdx -= mapIdx == mapLen-1;
363 mapIdxFrac -= mapIdx;
364 /* ...
365 fprintf(stderr, "%s: (%d,\ne% 31.15f) --> \n",
366 me, mapIdx, mapIdxFrac);
367 ... */
368 entData0 = mapData + mapIdx*entSize;
369 entData1 = mapData + (mapIdx+1)*entSize;
370 for (i=0; i<entLen; i++) {
371 val = ((1-mapIdxFrac)*mapLup(entData0, i) +
372 mapIdxFrac*mapLup(entData1, i));
373 outInsert(outData, i, val);
374 /* ...
375 fprintf(stderr, "f% 31.15f\n", val);
376 ... */
377 }
378 } else {
379 /* copy non-existent values from input to output */
380 for (i=0; i<entLen; i++) {
381 outInsert(outData, i, val);
382 }
383 }
384 inData += inSize;
385 outData += outSize;
386 if (multi) {
387 mapData += mapLen*entSize;
388 }
389 }
390 } else {
391 /* lookup table */
392 for (I=0; I<N; I++) {
393 val = inLoad(inData);
394 if (rescale) {
395 val = (range->min != range->max
396 ? AIR_AFFINE(range->min, val, range->max, domMin, domMax)( ((double)(domMax)-(domMin))*((double)(val)-(range->min))
/ ((double)(range->max)-(range->min)) + (domMin))
397 : domMin);
398 }
399 if (AIR_EXISTS(val)(((int)(!((val) - (val)))))) {
400 mapIdx = airIndexClamp(domMin, val, domMax, mapLen);
401 entData0 = mapData + mapIdx*entSize;
402 for (i=0; i<entLen; i++) {
403 outInsert(outData, i, mapLup(entData0, i));
404 }
405 } else {
406 /* copy non-existent values from input to output */
407 for (i=0; i<entLen; i++) {
408 outInsert(outData, i, val);
409 }
410 }
411 inData += inSize;
412 outData += outSize;
413 if (multi) {
414 mapData += mapLen*entSize;
415 }
416 }
417 }
418
419 return 0;
420}
421
422/*
423******** nrrdApply1DLut
424**
425** A "lut" is a simple lookup table: the data points are evenly spaced,
426** with cell-centering assumed, and there is no interpolation except
427** nearest neighbor. The axis min and max are used to determine the
428** range of values that can be mapped with the lut.
429**
430** Of the three kinds of 1-D maps, only luts can have output type block.
431**
432** If the lut nrrd is 1-D, then the output is a scalar nrrd with the
433** same dimension as the input. If the lut nrrd is 2-D, then each
434** value in the input is mapped to a vector of values from the lut,
435** which is always a scanline along axis 0.
436**
437** This allows lut length to be simply 1.
438*/
439int
440nrrdApply1DLut(Nrrd *nout, const Nrrd *nin,
441 const NrrdRange *_range, const Nrrd *nlut,
442 int typeOut, int rescale) {
443 static const char me[]="nrrdApply1DLut";
444 NrrdRange *range;
445 airArray *mop;
446
447 if (!(nout && nlut && nin)) {
448 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
449 return 1;
450 }
451 mop = airMopNew();
452 if (_range) {
453 range = nrrdRangeCopy(_range);
454 nrrdRangeSafeSet(range, nin, nrrdBlind8BitRangeState);
455 } else {
456 range = nrrdRangeNewSet(nin, nrrdBlind8BitRangeState);
457 }
458 airMopAdd(mop, range, (airMopper)nrrdRangeNix, airMopAlways);
459 if (_nrrdApply1DSetUp(nout, nin, range, nlut, kindLut, typeOut,
460 rescale, AIR_FALSE0 /* multi */)
461 || _nrrdApply1DLutOrRegMap(nout, nin, range, nlut, AIR_FALSE0 /* ramps */,
462 rescale, AIR_FALSE0 /* multi */)) {
463 biffAddf(NRRDnrrdBiffKey, "%s:", me);
464 airMopError(mop); return 1;
465 }
466 airMopOkay(mop);
467 return 0;
468}
469
470int
471nrrdApplyMulti1DLut(Nrrd *nout, const Nrrd *nin,
472 const NrrdRange *_range, const Nrrd *nmlut,
473 int typeOut, int rescale) {
474 static const char me[]="nrrdApplyMulti1DLut";
475 NrrdRange *range;
476 airArray *mop;
477
478 if (!(nout && nmlut && nin)) {
479 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
480 return 1;
481 }
482 mop = airMopNew();
483 if (_range) {
484 range = nrrdRangeCopy(_range);
485 nrrdRangeSafeSet(range, nin, nrrdBlind8BitRangeState);
486 } else {
487 range = nrrdRangeNewSet(nin, nrrdBlind8BitRangeState);
488 }
489 airMopAdd(mop, range, (airMopper)nrrdRangeNix, airMopAlways);
490 if (_nrrdApply1DSetUp(nout, nin, range, nmlut, kindLut, typeOut,
491 rescale, AIR_TRUE1 /* multi */)
492 || _nrrdApply1DLutOrRegMap(nout, nin, range, nmlut,
493 AIR_FALSE0 /* ramps */,
494 rescale, AIR_TRUE1 /* multi */)) {
495 biffAddf(NRRDnrrdBiffKey, "%s:", me);
496 airMopError(mop); return 1;
497 }
498 airMopOkay(mop);
499 return 0;
500}
501
502/*
503******** nrrdApply1DRegMap
504**
505** A regular map has data points evenly spaced, with node-centering and
506** and linear interpolation assumed. As with luts, the axis min and
507** max determine the range of values that are mapped. This function
508** is used in nrrdHistoEq() and is the basis of the very popular "unu rmap".
509**
510** If the lut nrrd is 1-D, then the output is a scalar nrrd with the
511** same dimension as the input. If the lut nrrd is 2-D, then each
512** value in the input is mapped to a linear weighting of vectors
513** from the map; the vectors are the scanlines along axis 0.
514**
515** NB: this function makes NO provisions for non-existent input values.
516** There won't be any memory errors, but the results are undefined.
517**
518** This allows map length to be simply 1.
519*/
520int
521nrrdApply1DRegMap(Nrrd *nout, const Nrrd *nin,
522 const NrrdRange *_range, const Nrrd *nmap,
523 int typeOut, int rescale) {
524 static const char me[]="nrrdApply1DRegMap";
525 NrrdRange *range;
526 airArray *mop;
527
528 if (!(nout && nmap && nin)) {
529 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
530 return 1;
531 }
532 mop = airMopNew();
533 if (_range) {
534 range = nrrdRangeCopy(_range);
535 nrrdRangeSafeSet(range, nin, nrrdBlind8BitRangeState);
536 } else {
537 range = nrrdRangeNewSet(nin, nrrdBlind8BitRangeState);
538 }
539 airMopAdd(mop, range, (airMopper)nrrdRangeNix, airMopAlways);
540 if (_nrrdApply1DSetUp(nout, nin, range, nmap, kindRmap, typeOut,
541 rescale, AIR_FALSE0 /* multi */)
542 || _nrrdApply1DLutOrRegMap(nout, nin, range, nmap, AIR_TRUE1 /* ramps */,
543 rescale, AIR_FALSE0 /* multi */)) {
544 biffAddf(NRRDnrrdBiffKey, "%s:", me);
545 airMopError(mop); return 1;
546 }
547 airMopOkay(mop);
548 return 0;
549}
550
551int
552nrrdApplyMulti1DRegMap(Nrrd *nout, const Nrrd *nin,
553 const NrrdRange *_range, const Nrrd *nmmap,
554 int typeOut, int rescale) {
555 static const char me[]="nrrdApplyMulti1DRegMap";
556 NrrdRange *range;
557 airArray *mop;
558
559 if (!(nout && nmmap && nin)) {
560 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
561 return 1;
562 }
563 mop = airMopNew();
564 if (_range) {
565 range = nrrdRangeCopy(_range);
566 nrrdRangeSafeSet(range, nin, nrrdBlind8BitRangeState);
567 } else {
568 range = nrrdRangeNewSet(nin, nrrdBlind8BitRangeState);
569 }
570 airMopAdd(mop, range, (airMopper)nrrdRangeNix, airMopAlways);
571 if (_nrrdApply1DSetUp(nout, nin, range, nmmap, kindRmap, typeOut,
572 rescale, AIR_TRUE1 /* multi */)
573 || _nrrdApply1DLutOrRegMap(nout, nin, range, nmmap, AIR_TRUE1 /* ramps */,
574 rescale, AIR_TRUE1 /* multi */)) {
575 biffAddf(NRRDnrrdBiffKey, "%s:", me);
576 airMopError(mop); return 1;
577 }
578 airMopOkay(mop);
579 return 0;
580}
581
582/*
583******** nrrd1DIrregMapCheck()
584**
585** return zero only for the valid forms of 1D irregular map.
586** imap must be 2D, both sizes >= 2, non-block-type, no non-existent
587** values in range. If the first point's position is non-existent,
588** than the first three points positions must be -inf, NaN, and +inf,
589** and none of the other points locations can be non-existent, and
590** they must increase monotonically. There must be at least two
591** points with existent positions.
592*/
593int
594nrrd1DIrregMapCheck(const Nrrd *nmap) {
595 static const char me[]="nrrd1DIrregMapCheck";
596 double (*mapLup)(const void *v, size_t I);
597 int i, entLen, mapLen, baseI;
598 size_t min[2], max[2];
599 Nrrd *nrange;
600
601 if (!nmap) {
602 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
603 return 1;
604 }
605 if (nrrdCheck(nmap)) {
606 biffAddf(NRRDnrrdBiffKey, "%s: ", me);
607 return 1;
608 }
609 if (nrrdTypeBlock == nmap->type) {
610 biffAddf(NRRDnrrdBiffKey, "%s: map is %s type, need scalar",
611 me, airEnumStr(nrrdType, nrrdTypeBlock));
612 return 1;
613 }
614 if (2 != nmap->dim) {
615 biffAddf(NRRDnrrdBiffKey, "%s: map needs to have dimension 2, not %d",
616 me, nmap->dim);
617 return 1;
618 }
619 entLen = AIR_CAST(unsigned int,nmap->axis[0].size)((unsigned int)(nmap->axis[0].size));
620 mapLen = AIR_CAST(unsigned int,nmap->axis[1].size)((unsigned int)(nmap->axis[1].size));
621 if (!( entLen >= 2 && mapLen >= 2 )) {
622 biffAddf(NRRDnrrdBiffKey, "%s: both map's axes sizes should be >= 2 (not %d,%d)",
623 me, entLen, mapLen);
624 return 1;
625 }
626 min[0] = 1; max[0] = nmap->axis[0].size-1;
627 min[1] = 0; max[1] = nmap->axis[1].size-1;
628 if (nrrdCrop(nrange=nrrdNew(), nmap, min, max)) {
629 biffAddf(NRRDnrrdBiffKey, "%s: couldn't crop to isolate range of map", me);
630 nrrdNuke(nrange); return 1;
631 }
632 if (nrrdHasNonExist(nrange)) {
633 biffAddf(NRRDnrrdBiffKey, "%s: map has non-existent values in its range", me);
634 nrrdNuke(nrange); return 1;
635 }
636 nrrdNuke(nrange);
637 mapLup = nrrdDLookup[nmap->type];
638 if (AIR_EXISTS(mapLup(nmap->data, 0))(((int)(!((mapLup(nmap->data, 0)) - (mapLup(nmap->data,
0))))))) {
639 baseI = 0;
640 } else {
641 baseI = 3;
642 if (!( mapLen >= 5 )) {
643 biffAddf(NRRDnrrdBiffKey, "%s: length of map w/ non-existent locations must "
644 "be >= 5 (not %d)", me, mapLen);
645 return 1;
646 }
647 if (!( airFP_NEG_INF == airFPClass_d(mapLup(nmap->data, 0*entLen)) &&
648 airFP_QNAN == airFPClass_d(mapLup(nmap->data, 1*entLen)) &&
649 airFP_POS_INF == airFPClass_d(mapLup(nmap->data, 2*entLen)) )) {
650 biffAddf(NRRDnrrdBiffKey, "%s: 1st entry's position non-existent, but position "
651 "of 1st three entries (%g:%d,%g:%d,%g:%d) not "
652 "-inf, NaN, and +inf", me,
653 mapLup(nmap->data, 0*entLen),
654 airFPClass_d(mapLup(nmap->data, 0*entLen)),
655 mapLup(nmap->data, 1*entLen),
656 airFPClass_d(mapLup(nmap->data, 1*entLen)),
657 mapLup(nmap->data, 2*entLen),
658 airFPClass_d(mapLup(nmap->data, 2*entLen)));
659 return 1;
660 }
661 }
662 for (i=baseI; i<mapLen; i++) {
663 if (!AIR_EXISTS(mapLup(nmap->data, i*entLen))(((int)(!((mapLup(nmap->data, i*entLen)) - (mapLup(nmap->
data, i*entLen))))))) {
664 biffAddf(NRRDnrrdBiffKey, "%s: entry %d has non-existent position", me, i);
665 return 1;
666 }
667 }
668 for (i=baseI; i<mapLen-1; i++) {
669 if (!( mapLup(nmap->data, i*entLen) < mapLup(nmap->data, (i+1)*entLen) )) {
670 biffAddf(NRRDnrrdBiffKey, "%s: map entry %d pos (%g) not < entry %d pos (%g)",
671 me, i, mapLup(nmap->data, i*entLen),
672 i+1, mapLup(nmap->data, (i+1)*entLen));
673 return 1;
674 }
675 }
676 return 0;
677}
678
679/*
680******** nrrd1DIrregAclCheck()
681**
682** returns zero only on valid accelerators for 1D irregular mappings
683*/
684int
685nrrd1DIrregAclCheck(const Nrrd *nacl) {
686 static const char me[]="nrrd1DIrregAclCheck";
687
688 if (!nacl) {
689 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
690 return 1;
691 }
692 if (nrrdCheck(nacl)) {
693 biffAddf(NRRDnrrdBiffKey, "%s: ", me);
694 return 1;
695 }
696 if (nrrdTypeUShort != nacl->type) {
697 biffAddf(NRRDnrrdBiffKey, "%s: type should be %s, not %s", me,
698 airEnumStr(nrrdType, nrrdTypeUShort),
699 airEnumStr(nrrdType, nacl->type));
700 return 1;
701 }
702 if (2 != nacl->dim) {
703 biffAddf(NRRDnrrdBiffKey, "%s: dimension should be 2, not %d", me, nacl->dim);
704 return 1;
705 }
706 if (!( nacl->axis[0].size == 2 && nacl->axis[1].size >= 2 )) {
707 char stmp1[AIR_STRLEN_SMALL(128+1)], stmp2[AIR_STRLEN_SMALL(128+1)];
708 biffAddf(NRRDnrrdBiffKey, "%s: sizes (%s,%s) not (2,>=2)", me,
709 airSprintSize_t(stmp1, nacl->axis[0].size),
710 airSprintSize_t(stmp2, nacl->axis[1].size));
711 return 1;
712 }
713
714 return 0;
715}
716
717/*
718** _nrrd1DIrregMapDomain()
719**
720** ALLOCATES an array of doubles storing the existent control point
721** locations, and sets its length in *poslenP. If there are the three
722** points with non-existent locations, these are ignored.
723**
724** Assumes that nrrd1DIrregMapCheck has been called on "nmap".
725*/
726double *
727_nrrd1DIrregMapDomain(int *posLenP, int *baseIP, const Nrrd *nmap) {
728 static const char me[]="_nrrd1DIrregMapDomain";
729 int i, entLen, baseI, posLen;
730 double *pos, (*mapLup)(const void *v, size_t I);
731
732 mapLup = nrrdDLookup[nmap->type];
733 baseI = AIR_EXISTS(mapLup(nmap->data, 0))(((int)(!((mapLup(nmap->data, 0)) - (mapLup(nmap->data,
0)))))) ? 0 : 3;
734 if (baseIP) {
735 *baseIP = baseI;
736 }
737 entLen = AIR_CAST(unsigned int,nmap->axis[0].size)((unsigned int)(nmap->axis[0].size));
738 posLen = AIR_CAST(unsigned int,nmap->axis[1].size)((unsigned int)(nmap->axis[1].size)) - baseI;
739 if (posLenP) {
740 *posLenP = posLen;
741 }
742 pos = (double*)malloc(posLen * sizeof(double));
743 if (!pos) {
744 biffAddf(NRRDnrrdBiffKey, "%s: couldn't allocate %d doubles\n", me, posLen);
745 return NULL((void*)0);
746 }
747 for (i=0; i<posLen; i++) {
748 pos[i] = mapLup(nmap->data, (baseI+i)*entLen);
749 }
750 return pos;
751}
752
753/*
754** _nrrd1DIrregFindInterval()
755**
756** The hard part of doing 1D irregular mapping: given an array of
757** control point locations, and a value, find which interval the value
758** lies in. The lowest and highest possible indices are given in
759** "loI" and "hiI". Results are undefined if these do not in fact
760** bound the location of correct interval, or if loI > hiI, or if the
761** query positon "p" is not in the domain vector "pos". Intervals are
762** identified by the integral index of the LOWER of the two control
763** points spanning the interval.
764**
765** This imposes the same structure of half-open intervals that
766** is done by airIndex(). That is, a value p is in interval i
767** if pos[i] <= p < pos[i+1] for all but the last interval, and
768** pos[i] <= p <= pos[i+1] for the last interval (in which case
769** i == hiI)
770*/
771int
772_nrrd1DIrregFindInterval(const double *pos, double p, int loI, int hiI) {
773 int midI;
774
775 /*
776 fprintf(stderr, "##%s(%g): hi: %d/%g-%g | %d/%g-%g\n",
777 "_nrrd1DIrregFindInterval", p,
778 loI, pos[loI], pos[loI+1],
779 hiI, pos[hiI], pos[hiI+1]);
780 */
781 while (loI < hiI) {
782 midI = (loI + hiI)/2;
783 if ( pos[midI] <= p && ((midI < hiI && p < pos[midI+1]) ||
784 (midI == hiI && p <= pos[midI+1])) ) {
785 /* p is between (pos[midI],pos[midI+1]): we're done */
786 loI = hiI = midI;
787 } else if (pos[midI] > p) {
788 /* p is below interval midI: midI-1 is valid upper bound */
789 hiI = midI-1;
790 } else {
791 /* p is above interval midI: midI+1 is valid lower bound */
792 loI = midI+1;
793 }
794 /*
795 fprintf(stderr, "##%s(%g): %d/%g-%g | %d/%g-%g | %d/%g-%g\n",
796 "_nrrd1DIrregFindInterval", p,
797 loI, pos[loI], pos[loI+1],
798 midI, pos[midI], pos[midI+1],
799 hiI, pos[hiI], pos[hiI+1]);
800 */
801 }
802 return loI;
803}
804
805/*
806******** nrrd1DIrregAclGenerate()
807**
808** Generates the "acl" that is used to speed up the action of
809** nrrdApply1DIrregMap(). Basically, the domain of the map
810** is quantized into "acllen" bins, and for each bin, the
811** lowest and highest possible map interval is stored. This
812** either obviates or speeds up the task of finding which
813** interval contains a given value.
814**
815** Assumes that nrrd1DIrregMapCheck has been called on "nmap".
816*/
817int
818nrrd1DIrregAclGenerate(Nrrd *nacl, const Nrrd *nmap, size_t aclLen) {
819 static const char me[]="nrrd1DIrregAclGenerate";
820 int posLen;
821 unsigned int ii;
822 unsigned short *acl;
823 double lo, hi, min, max, *pos;
824
825 if (!(nacl && nmap)) {
1
Taking false branch
→
826 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
827 return 1;
828 }
829 if (!(aclLen >= 2)) {
2
←
Assuming 'aclLen' is >= 2
→
3
←
Taking false branch
→
830 char stmp[AIR_STRLEN_SMALL(128+1)];
831 biffAddf(NRRDnrrdBiffKey, "%s: given acl length (%s) is too small", me,
832 airSprintSize_t(stmp, aclLen));
833 return 1;
834 }
835 if (nrrdMaybeAlloc_va(nacl, nrrdTypeUShort, 2,
4
←
Taking false branch
→
836 AIR_CAST(size_t, 2)((size_t)(2)), AIR_CAST(size_t, aclLen)((size_t)(aclLen)))) {
837 biffAddf(NRRDnrrdBiffKey, "%s: ", me);
838 return 1;
839 }
840 acl = (unsigned short *)nacl->data;
841 pos = _nrrd1DIrregMapDomain(&posLen, NULL((void*)0), nmap);
842 if (!pos) {
5
←
Taking false branch
→
843 biffAddf(NRRDnrrdBiffKey, "%s: couldn't determine domain", me);
844 return 1;
845 }
846 nacl->axis[1].min = min = pos[0];
6
←
Assigned value is garbage or undefined
847 nacl->axis[1].max = max = pos[posLen-1];
848 for (ii=0; ii<=aclLen-1; ii++) {
849 lo = AIR_AFFINE(0, ii, aclLen, min, max)( ((double)(max)-(min))*((double)(ii)-(0)) / ((double)(aclLen
)-(0)) + (min));
850 hi = AIR_AFFINE(0, ii+1, aclLen, min, max)( ((double)(max)-(min))*((double)(ii+1)-(0)) / ((double)(aclLen
)-(0)) + (min));
851 acl[0 + 2*ii] = _nrrd1DIrregFindInterval(pos, lo, 0, posLen-2);
852 acl[1 + 2*ii] = _nrrd1DIrregFindInterval(pos, hi, 0, posLen-2);
853 }
854 free(pos);
855
856 return 0;
857}
858
859/*
860******** nrrdApply1DIrregMap()
861**
862** Linear interpolation between irregularly spaced control points.
863** Obviously, the location of the control point has to be given
864** explicitly. The map nrrd must have dimension 2, and each
865** control point is represented by a scanline along axis 0. The
866** first value is the position of the control point, and the remaining
867** value(s) are linearly weighted according to the position of the
868** input value among the control point locations.
869**
870** To allow "coloring" of non-existent values -inf, NaN, and +inf, if
871** the very first value of the map (the location of the first control
872** point) is non-existent, then the first three control point locations
873** must be -inf, NaN, and +inf, in that order, and the information
874** about these points will be used for corresponding input values.
875** Doing this makes everything slower, however, because airFPClass_f()
876** is called on every single value.
877**
878** This assumes that nrrd1DIrregMapCheck has been called on "nmap",
879** and that nrrd1DIrregAclCheck has been called on "nacl" (if it is
880** non-NULL).
881*/
882int
883nrrdApply1DIrregMap(Nrrd *nout, const Nrrd *nin, const NrrdRange *_range,
884 const Nrrd *nmap, const Nrrd *nacl,
885 int typeOut, int rescale) {
886 static const char me[]="nrrdApply1DIrregMap";
887 size_t N, I;
888 int i, *acl, entLen, posLen, aclLen, mapIdx, aclIdx,
889 entSize, colSize, inSize, lo, hi, baseI;
890 double val, *pos, domMin, domMax, mapIdxFrac,
891 (*mapLup)(const void *v, size_t I),
892 (*inLoad)(const void *v), (*outInsert)(void *v, size_t I, double d);
893 char *inData, *outData, *entData0, *entData1;
894 NrrdRange *range;
895 airArray *mop;
896
897 /*
898 fprintf(stderr, "!%s: rescale = %d\n", me, rescale);
899 */
900 if (!(nout && nmap && nin)) {
901 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
902 return 1;
903 }
904 mop = airMopNew();
905 if (_range) {
906 range = nrrdRangeCopy(_range);
907 nrrdRangeSafeSet(range, nin, nrrdBlind8BitRangeState);
908 } else {
909 range = nrrdRangeNewSet(nin, nrrdBlind8BitRangeState);
910 }
911 airMopAdd(mop, range, (airMopper)nrrdRangeNix, airMopAlways);
912 if (_nrrdApply1DSetUp(nout, nin, range, nmap,
913 kindImap, typeOut, rescale, AIR_FALSE0 /* multi */)) {
914 biffAddf(NRRDnrrdBiffKey, "%s:", me);
915 airMopError(mop); return 1;
916 }
917 if (nacl && nrrd1DIrregAclCheck(nacl)) {
918 biffAddf(NRRDnrrdBiffKey, "%s: given acl isn't valid", me);
919 airMopError(mop); return 1;
920 }
921
922 if (nacl) {
923 acl = (int *)nacl->data;
924 aclLen = AIR_CAST(unsigned int,nacl->axis[1].size)((unsigned int)(nacl->axis[1].size));
925 } else {
926 acl = NULL((void*)0);
927 aclLen = 0;
928 }
929 pos = _nrrd1DIrregMapDomain(&posLen, &baseI, nmap);
930 if (!pos) {
931 biffAddf(NRRDnrrdBiffKey, "%s: couldn't determine domain", me);
932 airMopError(mop); return 1;
933 }
934 airMopAdd(mop, pos, airFree, airMopAlways);
935 mapLup = nrrdDLookup[nmap->type];
936
937 inData = (char *)nin->data;
938 inLoad = nrrdDLoad[nin->type];
939 inSize = AIR_CAST(unsigned int,nrrdElementSize(nin))((unsigned int)(nrrdElementSize(nin)));
940 mapLup = nrrdDLookup[nmap->type];
941 entLen = AIR_CAST(unsigned int,nmap->axis[0].size)((unsigned int)(nmap->axis[0].size)); /* entLen is really 1 + entry length */
942 entSize = entLen*AIR_CAST(unsigned int,nrrdElementSize(nmap))((unsigned int)(nrrdElementSize(nmap)));
943 colSize = (entLen-1)*AIR_CAST(unsigned int,nrrdTypeSize[typeOut])((unsigned int)(nrrdTypeSize[typeOut]));
944 outData = (char *)nout->data;
945 outInsert = nrrdDInsert[nout->type];
946 domMin = pos[0];
947 domMax = pos[posLen-1];
948 /*
949 fprintf(stderr, "!%s: domMin, domMax = %g, %g\n", me, domMin, domMax);
950 */
951
952 N = nrrdElementNumber(nin);
953 for (I=0;
954 I<N;
955 I++, inData += inSize, outData += colSize) {
956 val = inLoad(inData);
957 /*
958 fprintf(stderr, "!%s: (%d) val = % 31.15f\n", me, (int)I, val);
959 */
960 if (!AIR_EXISTS(val)(((int)(!((val) - (val)))))) {
961 /* got a non-existent value */
962 if (baseI) {
963 /* and we know how to deal with them */
964 switch (airFPClass_d(val)) {
965 case airFP_NEG_INF:
966 mapIdx = 0;
967 break;
968 case airFP_SNAN:
969 case airFP_QNAN:
970 mapIdx = 1;
971 break;
972 case airFP_POS_INF:
973 mapIdx = 2;
974 break;
975 default:
976 mapIdx = 0;
977 fprintf(stderr__stderrp, "%s: PANIC: non-existent value/class %g/%d "
978 "not handled\n",
979 me, val, airFPClass_d(val));
980 exit(1);
981 }
982 entData0 = (char*)(nmap->data) + mapIdx*entSize;
983 for (i=1; i<entLen; i++) {
984 outInsert(outData, i-1, mapLup(entData0, i));
985 }
986 continue; /* we're done! (with this value) */
987 } else {
988 /* we don't know how to properly deal with this non-existent value:
989 we use the first entry, and then fall through to code below */
990 mapIdx = 0;
991 mapIdxFrac = 0.0;
992 }
993 } else {
994 /* we have an existent value */
995 if (rescale) {
996 val = (range->min != range->max
997 ? AIR_AFFINE(range->min, val, range->max, domMin, domMax)( ((double)(domMax)-(domMin))*((double)(val)-(range->min))
/ ((double)(range->max)-(range->min)) + (domMin))
998 : domMin);
999 }
1000 val = AIR_CLAMP(domMin, val, domMax)((val) < (domMin) ? (domMin) : ((val) > (domMax) ? (domMax
) : (val)));
1001 if (acl) {
1002 aclIdx = airIndex(domMin, val, domMax, aclLen);
1003 lo = acl[0 + 2*aclIdx];
1004 hi = acl[1 + 2*aclIdx];
1005 } else {
1006 lo = 0;
1007 hi = posLen-2;
1008 }
1009 if (lo < hi) {
1010 mapIdx = _nrrd1DIrregFindInterval(pos, val, lo, hi);
1011 } else {
1012 /* acl did its job ==> lo == hi */
1013 mapIdx = lo;
1014 }
1015 /*
1016 fprintf(stderr, "!%s: --> val = %g; lo,hi = %d,%d, mapIdx = %d\n",
1017 me, val, lo, hi, mapIdx);
1018 */
1019 }
1020 mapIdxFrac = AIR_AFFINE(pos[mapIdx], val, pos[mapIdx+1], 0.0, 1.0)( ((double)(1.0)-(0.0))*((double)(val)-(pos[mapIdx])) / ((double
)(pos[mapIdx+1])-(pos[mapIdx])) + (0.0));
1021 /*
1022 fprintf(stderr, "!%s: mapIdxFrac = %g\n", me, mapIdxFrac);
1023 */
1024 entData0 = (char*)(nmap->data) + (baseI+mapIdx)*entSize;
1025 entData1 = (char*)(nmap->data) + (baseI+mapIdx+1)*entSize;
1026 for (i=1; i<entLen; i++) {
1027 val = ((1-mapIdxFrac)*mapLup(entData0, i) +
1028 mapIdxFrac*mapLup(entData1, i));
1029 /*
1030 fprintf(stderr, "!%s: %g * %g + %g * %g = %g\n", me,
1031 1-mapIdxFrac, mapLup(entData0, i),
1032 mapIdxFrac, mapLup(entData1, i), val);
1033 */
1034 outInsert(outData, i-1, val);
1035 }
1036 }
1037 airMopOkay(mop);
1038 return 0;
1039}
1040
1041/*
1042******** nrrdApply1DSubstitution
1043**
1044** A "subst" is a substitution table, i.e. a list of pairs that
1045** describes what values should be substituted with what (substitution
1046** rules). So, nsubst must be a scalar 2xN array. The output is a
1047** copy of the input with values substituted using this table.
1048**
1049** Unlike with lookup tables and maps (irregular and regular), we
1050** aren't at liberty to change the dimensionality of the data (can't
1051** substitute a grayscale with a color). The ability to apply
1052** substitutions to non-scalar data will be in a different function.
1053** Also unlike lookup tables and maps, the output type is the SAME as
1054** the input type; the output does NOT inherit the type of the
1055** substitution
1056*/
1057int
1058nrrdApply1DSubstitution(Nrrd *nout, const Nrrd *nin, const Nrrd *_nsubst) {
1059 static const char me[]="nrrdApply1DSubstitution";
1060 double (*lup)(const void *, size_t);
1061 double (*ins)(void *, size_t, double);
1062 Nrrd *nsubst;
1063 double val, *subst;
1064 size_t ii, jj, num, asize0, asize1;
1065 int changed;
1066 airArray *mop;
1067
1068 if (!(nout && _nsubst && nin)) {
1069 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
1070 return 1;
1071 }
1072 if (nrrdTypeBlock == nin->type || nrrdTypeBlock == _nsubst->type) {
1073 biffAddf(NRRDnrrdBiffKey, "%s: input or substitution type is %s, need scalar",
1074 me, airEnumStr(nrrdType, nrrdTypeBlock));
1075 return 1;
1076 }
1077 if (2 != _nsubst->dim) {
1078 biffAddf(NRRDnrrdBiffKey, "%s: substitution table has to be 2-D, not %d-D",
1079 me, _nsubst->dim);
1080 return 1;
1081 }
1082 /* nrrdAxisInfoGet_va(_nsubst, nrrdAxisInfoSize, &asize0, &asize1); */
1083 asize0 = _nsubst->axis[0].size;
1084 asize1 = _nsubst->axis[1].size;
1085 if (2 != asize0) {
1086 biffAddf(NRRDnrrdBiffKey, "%s: substitution table has to be 2xN, not %uxN",
1087 me, AIR_CAST(unsigned int, asize0)((unsigned int)(asize0)));
1088 return 1;
1089 }
1090 if (nout != nin) {
1091 if (nrrdCopy(nout, nin)) {
1092 biffAddf(NRRDnrrdBiffKey, "%s: couldn't initialize by copy to output", me);
1093 return 1;
1094 }
1095 }
1096
1097 mop = airMopNew();
1098 nsubst = nrrdNew();
1099 airMopAdd(mop, nsubst, (airMopper)nrrdNuke, airMopAlways);
1100 if (nrrdConvert(nsubst, _nsubst, nrrdTypeDouble)) {
1101 biffAddf(NRRDnrrdBiffKey, "%s: couldn't create double copy of substitution table",
1102 me);
1103 airMopError(mop); return 1;
1104 }
1105 lup = nrrdDLookup[nout->type];
1106 ins = nrrdDInsert[nout->type];
1107 subst = (double *)nsubst->data;
1108 num = nrrdElementNumber(nout);
1109 for (ii=0; ii<num; ii++) {
1110 val = lup(nout->data, ii);
1111 changed = AIR_FALSE0;
1112 for (jj=0; jj<asize1; jj++) {
1113 if (val == subst[jj*2+0]) {
1114 val = subst[jj*2+1];
1115 changed = AIR_TRUE1;
1116 }
1117 }
1118 if (changed) {
1119 ins(nout->data, ii, val);
1120 }
1121 }
1122
1123 airMopOkay(mop);
1124 return 0;
1125}