GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: src/nrrd/apply2D.c Lines: 0 151 0.0 %
Date: 2017-05-26 Branches: 0 84 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
#include "nrrd.h"
25
#include "privateNrrd.h"
26
27
enum {
28
  kindLut=0,
29
  kindRmap=1
30
};
31
32
/*
33
** _nrrdApply2DSetUp()
34
**
35
** some error checking and initializing needed for 2D LUTS and regular
36
** maps.  The intent is that if this succeeds, then there is no need
37
** for any further error checking.
38
**
39
** The only thing this function DOES is allocate the output nrrd, and
40
** set meta information.  The rest is just error checking.
41
**
42
** The given NrrdRange has to be fleshed out by the caller: it can't
43
** be NULL, and both range->min and range->max must exist.
44
*/
45
int
46
_nrrdApply2DSetUp(Nrrd *nout, const Nrrd *nin,
47
                  const NrrdRange *range0, const NrrdRange *range1,
48
                  const Nrrd *nmap, int kind, int typeOut,
49
                  int rescale0, int rescale1) {
50
  static const char me[]="_nrrdApply2DSetUp";
51
  char *mapcnt;
52
  char nounStr[][AIR_STRLEN_SMALL]={"2D lut",
53
                                    "2D regular map"};
54
  char verbStr[][AIR_STRLEN_SMALL]={"lut2",
55
                                    "rmap2"};
56
  int mapAxis, copyMapAxis0=AIR_FALSE, axisMap[NRRD_DIM_MAX];
57
  unsigned int dim, entLen;
58
  size_t size[NRRD_DIM_MAX];
59
  double domMin, domMax;
60
61
  if (nout == nin) {
62
    biffAddf(NRRD, "%s: due to laziness, nout==nin always disallowed", me);
63
    return 1;
64
  }
65
  if (airEnumValCheck(nrrdType, typeOut)) {
66
    biffAddf(NRRD, "%s: invalid requested output type %d", me, typeOut);
67
    return 1;
68
  }
69
  if (nrrdTypeBlock == nin->type || nrrdTypeBlock == typeOut) {
70
    biffAddf(NRRD, "%s: input or requested output type is %s, need scalar",
71
             me, airEnumStr(nrrdType, nrrdTypeBlock));
72
    return 1;
73
  }
74
  if (!(2 == nin->axis[0].size)) {
75
    char stmp[AIR_STRLEN_SMALL];
76
    biffAddf(NRRD, "%s: input axis[0] must have size 2 (not %s)", me,
77
             airSprintSize_t(stmp, nin->axis[0].size));
78
    return 1;
79
  }
80
  if (!(nin->dim > 1)) {
81
    biffAddf(NRRD, "%s: input dimension must be > 1 (not %u)", me, nin->dim);
82
    return 1;
83
  }
84
  if (rescale0 && !(range0
85
                    && AIR_EXISTS(range0->min)
86
                    && AIR_EXISTS(range0->max))) {
87
    biffAddf(NRRD, "%s: want axis 0 rescaling but didn't get range, or "
88
             "not both range->{min,max} exist", me);
89
    return 1;
90
  }
91
  if (rescale1 && !(range1
92
                    && AIR_EXISTS(range1->min)
93
                    && AIR_EXISTS(range1->max))) {
94
    biffAddf(NRRD, "%s: want axis 1 rescaling but didn't get range, or "
95
             "not both range->{min,max} exist", me);
96
    return 1;
97
  }
98
  mapAxis = nmap->dim - 2;
99
  if (!(0 == mapAxis || 1 == mapAxis)) {
100
    biffAddf(NRRD, "%s: dimension of %s should be 2 or 3, not %d",
101
             me, nounStr[kind], nmap->dim);
102
    return 1;
103
  }
104
  copyMapAxis0 = (1 == mapAxis);
105
  domMin = _nrrdApplyDomainMin(nmap, AIR_FALSE, mapAxis);
106
  domMax = _nrrdApplyDomainMax(nmap, AIR_FALSE, mapAxis);
107
  if (!( domMin < domMax )) {
108
    biffAddf(NRRD, "%s: (axis %d) domain min (%g) not less than max (%g)", me,
109
             mapAxis, domMin, domMax);
110
    return 1;
111
  }
112
  domMin = _nrrdApplyDomainMin(nmap, AIR_FALSE, mapAxis+1);
113
  domMax = _nrrdApplyDomainMax(nmap, AIR_FALSE, mapAxis+1);
114
  if (!( domMin < domMax )) {
115
    biffAddf(NRRD, "%s: (axis %d) domain min (%g) not less than max (%g)", me,
116
             mapAxis+1, domMin, domMax);
117
    return 1;
118
  }
119
  if (nrrdHasNonExist(nmap)) {
120
    biffAddf(NRRD, "%s: %s nrrd has non-existent values",
121
             me, nounStr[kind]);
122
    return 1;
123
  }
124
  entLen = mapAxis ? AIR_CAST(unsigned int, nmap->axis[0].size) : 1;
125
  if (mapAxis + nin->dim - 1 > NRRD_DIM_MAX) {
126
    biffAddf(NRRD, "%s: input nrrd dim %d through non-scalar %s exceeds "
127
             "NRRD_DIM_MAX %d",
128
             me, nin->dim, nounStr[kind], NRRD_DIM_MAX);
129
    return 1;
130
  }
131
  if (mapAxis) {
132
    size[0] = entLen;
133
    axisMap[0] = -1;
134
  }
135
  for (dim=1; dim<nin->dim; dim++) {
136
    size[dim-1+mapAxis] = nin->axis[dim].size;
137
    axisMap[dim-1+mapAxis] = dim;
138
  }
139
  /*
140
  fprintf(stderr, "##%s: pre maybe alloc: nout->data = %p\n", me, nout->data);
141
  for (dim=0; dim<mapAxis + nin->dim - 1; dim++) {
142
    fprintf(stderr, "    size[%d] = %d\n", dim, (int)size[dim]);
143
  }
144
  fprintf(stderr, "   nout->dim = %u; nout->type = %d = %s; sizes = %u,%u\n",
145
          nout->dim, nout->type,
146
          airEnumStr(nrrdType, nout->type),
147
          AIR_CAST(unsigned int, nout->axis[0].size),
148
          AIR_CAST(unsigned int, nout->axis[1].size));
149
  fprintf(stderr, "   typeOut = %d = %s\n", typeOut,
150
          airEnumStr(nrrdType, typeOut));
151
  */
152
  if (nrrdMaybeAlloc_nva(nout, typeOut, nin->dim - 1 + mapAxis, size)) {
153
    biffAddf(NRRD, "%s: couldn't allocate output nrrd", me);
154
    return 1;
155
  }
156
  /*
157
  fprintf(stderr, "   nout->dim = %d; nout->type = %d = %s\n",
158
          nout->dim, nout->type,
159
          airEnumStr(nrrdType, nout->type));
160
  for (dim=0; dim<nout->dim; dim++) {
161
    fprintf(stderr, "    size[%d] = %d\n", dim, (int)nout->axis[dim].size);
162
  }
163
  fprintf(stderr, "##%s: post maybe alloc: nout->data = %p\n", me, nout->data);
164
  */
165
  if (nrrdAxisInfoCopy(nout, nin, axisMap, NRRD_AXIS_INFO_NONE)) {
166
    biffAddf(NRRD, "%s: trouble copying axis info", me);
167
    return 1;
168
  }
169
  if (copyMapAxis0) {
170
    _nrrdAxisInfoCopy(nout->axis + 0, nmap->axis + 0,
171
                      NRRD_AXIS_INFO_SIZE_BIT);
172
  }
173
174
  mapcnt = _nrrdContentGet(nmap);
175
  if (nrrdContentSet_va(nout, verbStr[kind], nin, "%s", mapcnt)) {
176
    biffAddf(NRRD, "%s:", me);
177
    free(mapcnt); return 1;
178
  }
179
  free(mapcnt);
180
  nrrdBasicInfoInit(nout, (NRRD_BASIC_INFO_DATA_BIT
181
                           | NRRD_BASIC_INFO_TYPE_BIT
182
                           | NRRD_BASIC_INFO_BLOCKSIZE_BIT
183
                           | NRRD_BASIC_INFO_DIMENSION_BIT
184
                           | NRRD_BASIC_INFO_CONTENT_BIT));
185
  return 0;
186
}
187
188
/*
189
** _nrrdApply2DLutOrRegMap()
190
**
191
** the guts of nrrdApply2DLut and nrrdApply2DRegMap
192
**
193
** yikes, does NOT use biff, since we're only supposed to be called
194
** after copious error checking.
195
**
196
** FOR INSTANCE, this allows nout == nin, which could be a big
197
** problem if mapAxis == 1.
198
**
199
** we don't need a typeOut arg because nout has already been allocated
200
** as some specific type; we'll look at that.
201
**
202
** NOTE: non-existent values get passed through regular maps and luts
203
** "unchanged".  However, if the output type is integral, the results
204
** are probaby undefined.  HEY: there is currently no warning message
205
** or error handling based on nrrdStateDisallowIntegerNonExist, but
206
** there really should be.
207
*/
208
int
209
_nrrdApply2DLutOrRegMap(Nrrd *nout, const Nrrd *nin,
210
                        const NrrdRange *range0, const NrrdRange *range1,
211
                        const Nrrd *nmap, int ramps,
212
                        int rescale0, int rescale1) {
213
  static const char me[]="_nrrdApply2DLutOrRegMap";
214
  char *inData, *outData, *mapData, *entData;
215
  size_t N, I;
216
  double (*inLoad)(const void *v), (*mapLup)(const void *v, size_t I),
217
    (*outInsert)(void *v, size_t I, double d),
218
    val0, val1, domMin0, domMax0, domMin1, domMax1;
219
  unsigned int i, mapAxis, mapLen0, mapLen1, mapIdx0, mapIdx1,
220
    entSize, entLen, inSize, outSize;
221
222
  mapAxis = nmap->dim - 2;             /* axis of nmap containing entries */
223
  mapData = (char *)nmap->data;        /* map data, as char* */
224
                                       /* low end of map domain */
225
  domMin0 = _nrrdApplyDomainMin(nmap, ramps, mapAxis + 0);
226
  domMin1 = _nrrdApplyDomainMin(nmap, ramps, mapAxis + 1);
227
                                       /* high end of map domain */
228
  domMax0 = _nrrdApplyDomainMax(nmap, ramps, mapAxis + 0);
229
  domMax1 = _nrrdApplyDomainMax(nmap, ramps, mapAxis + 1);
230
  mapLen0 = AIR_CAST(unsigned int, nmap->axis[mapAxis+0].size);   /* number of entries in map axis 0 */
231
  mapLen1 = AIR_CAST(unsigned int, nmap->axis[mapAxis+1].size);   /* number of entries in map axis 1 */
232
  mapLup = nrrdDLookup[nmap->type];    /* how to get doubles out of map */
233
  inData = (char *)nin->data;          /* input data, as char* */
234
  inLoad = nrrdDLoad[nin->type];       /* how to get doubles out of nin */
235
  inSize = AIR_CAST(unsigned int, nrrdElementSize(nin));       /* size of one input value */
236
  outData = (char *)nout->data;        /* output data, as char* */
237
  outInsert = nrrdDInsert[nout->type]; /* putting doubles into output */
238
  entLen = (mapAxis                    /* number of elements in one entry */
239
            ? AIR_CAST(unsigned int, nmap->axis[0].size)
240
            : 1);
241
  outSize = entLen*AIR_CAST(unsigned int, nrrdElementSize(nout)); /* size of entry in output */
242
  entSize = entLen*AIR_CAST(unsigned int, nrrdElementSize(nmap)); /* size of entry in map */
243
244
  /*
245
  fprintf(stderr, "!%s: entLen = %u, mapLen = %u,%u\n", me,
246
          entLen, mapLen0, mapLen1);
247
  */
248
249
  N = nrrdElementNumber(nin)/2;       /* number of value pairs to be mapped */
250
  /* _VV = 1; */
251
  if (ramps) {
252
    fprintf(stderr, "%s: PANIC: unimplemented\n", me);
253
    exit(1);
254
  } else {
255
    /* lookup table */
256
    for (I=0; I<N; I++) {
257
      val0 = inLoad(inData + 0*inSize);
258
      val1 = inLoad(inData + 1*inSize);
259
      if (rescale0) {
260
        val0 = AIR_AFFINE(range0->min, val0, range0->max, domMin0, domMax0);
261
      }
262
      if (rescale1) {
263
        val1 = AIR_AFFINE(range1->min, val1, range1->max, domMin1, domMax1);
264
      }
265
      if (AIR_EXISTS(val0) && AIR_EXISTS(val1)) {
266
        mapIdx0 = airIndexClamp(domMin0, val0, domMax0, mapLen0);
267
        mapIdx1 = airIndexClamp(domMin1, val1, domMax1, mapLen1);
268
        entData = mapData + entSize*(mapIdx0 + mapLen0*mapIdx1);
269
        for (i=0; i<entLen; i++) {
270
          outInsert(outData, i, mapLup(entData, i));
271
        }
272
      } else {
273
        /* copy non-existent values from input to output */
274
        for (i=0; i<entLen; i++) {
275
          outInsert(outData, i, val0 + val1);  /* HEY this is weird */
276
        }
277
      }
278
      inData += 2*inSize;
279
      outData += outSize;
280
    }
281
  }
282
283
  return 0;
284
}
285
286
/*
287
******** nrrdApply2DLut
288
**
289
** A "lut" is a simple lookup table: the data points are evenly spaced,
290
** with cell-centering assumed, and there is no interpolation except
291
** nearest neighbor.  The axis min and max are used to determine the
292
** range of values that can be mapped with the lut.
293
**
294
** Of the three kinds of 1-D maps, only luts can have output type block.
295
**
296
** If the lut nrrd is 1-D, then the output is a scalar nrrd with the
297
** same dimension as the input.  If the lut nrrd is 2-D, then each
298
** value in the input is mapped to a vector of values from the lut,
299
** which is always a scanline along axis 0.
300
**
301
** This allows lut length to be simply 1.
302
*/
303
int
304
nrrdApply2DLut(Nrrd *nout, const Nrrd *nin, unsigned int domainAxis,
305
               const NrrdRange *_range0, const NrrdRange *_range1,
306
               const Nrrd *nlut,
307
               int typeOut, int rescale0, int rescale1) {
308
  static const char me[]="nrrdApply2DLut";
309
  NrrdRange *range0, *range1;
310
  airArray *mop;
311
  Nrrd *nin0, *nin1;
312
313
  if (!(nout && nlut && nin)) {
314
    biffAddf(NRRD, "%s: got NULL pointer (%p,%p,%p)", me,
315
             AIR_CAST(void*,nout), AIR_CVOIDP(nlut), AIR_CVOIDP(nin));
316
    return 1;
317
  }
318
  if (0 != domainAxis) {
319
    biffAddf(NRRD, "%s: sorry, domainAxis must currently be 0 (not %u)", me,
320
             domainAxis);
321
    return 1;
322
  }
323
  mop = airMopNew();
324
  if (_range0) {
325
    range0 = nrrdRangeCopy(_range0);
326
    airMopAdd(mop, nin0 = nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
327
    if (nrrdSlice(nin0, nin, 0, 0)) {
328
      biffAddf(NRRD, "%s: trouble learning range 0", me);
329
      airMopError(mop); return 1;
330
    }
331
    nrrdRangeSafeSet(range0, nin0, nrrdBlind8BitRangeState);
332
  } else {
333
    range0 = nrrdRangeNewSet(nin, nrrdBlind8BitRangeState);
334
  }
335
  if (_range1) {
336
    range1 = nrrdRangeCopy(_range1);
337
    airMopAdd(mop, nin1 = nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
338
    if (nrrdSlice(nin1, nin, 0, 1)) {
339
      biffAddf(NRRD, "%s: trouble learning range 1", me);
340
      airMopError(mop); return 1;
341
    }
342
    nrrdRangeSafeSet(range1, nin1, nrrdBlind8BitRangeState);
343
  } else {
344
    range1 = nrrdRangeNewSet(nin, nrrdBlind8BitRangeState);
345
  }
346
  airMopAdd(mop, range0, (airMopper)nrrdRangeNix, airMopAlways);
347
  airMopAdd(mop, range1, (airMopper)nrrdRangeNix, airMopAlways);
348
  if (_nrrdApply2DSetUp(nout, nin, range0, range1,
349
                        nlut, kindLut, typeOut, rescale0, rescale1)
350
      || _nrrdApply2DLutOrRegMap(nout, nin, range0, range1,
351
                                 nlut, AIR_FALSE, rescale0, rescale1)) {
352
    biffAddf(NRRD, "%s:", me);
353
    airMopError(mop); return 1;
354
  }
355
  airMopOkay(mop);
356
  return 0;
357
}
358
359
/*
360
******** nrrdApply2DRegMap
361
**
362
363
no, this doesn't actually exist yet....
364
365
int
366
nrrdApply2DRegMap(Nrrd *nout, const Nrrd *nin,
367
                  const NrrdRange *_range0, const NrrdRange *_range1,
368
                  const Nrrd *nmap, int typeOut, int rescale) {
369
  static const char me[]="nrrdApply2DRegMap";
370
  NrrdRange *range;
371
  airArray *mop;
372
373
  if (!(nout && nmap && nin)) {
374
    biffAddf(NRRD, "%s: got NULL pointer", me);
375
    return 1;
376
  }
377
  mop = airMopNew();
378
  if (_range) {
379
    range = nrrdRangeCopy(_range);
380
    nrrdRangeSafeSet(range, nin, nrrdBlind8BitRangeState);
381
  } else {
382
    range = nrrdRangeNewSet(nin, nrrdBlind8BitRangeState);
383
  }
384
  airMopAdd(mop, range, (airMopper)nrrdRangeNix, airMopAlways);
385
  if (_nrrdApply2DSetUp(nout, nin, range, nmap, kindRmap, typeOut,
386
                        rescale, AIR_FALSE)
387
      || _nrrdApply2DLutOrRegMap(nout, nin, range, nmap, AIR_TRUE,
388
                                 rescale, AIR_FALSE)) {
389
    biffAddf(NRRD, "%s:", me);
390
    airMopError(mop); return 1;
391
  }
392
  airMopOkay(mop);
393
  return 0;
394
}
395
396
*/
397