GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: src/nrrd/range.c Lines: 36 172 20.9 %
Date: 2017-05-26 Branches: 12 104 11.5 %

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
NrrdRange *
28
nrrdRangeNew(double min, double max) {
29
  NrrdRange *range;
30
31
4
  range = (NrrdRange *)calloc(1, sizeof(NrrdRange));
32
2
  if (range) {
33
2
    range->min = min;
34
2
    range->max = max;
35
2
    range->hasNonExist = nrrdHasNonExistUnknown;
36
2
  }
37
2
  return range;
38
}
39
40
NrrdRange *
41
nrrdRangeCopy(const NrrdRange *rin) {
42
  NrrdRange *rout=NULL;
43
44
  if (rin) {
45
    rout = nrrdRangeNew(rin->min, rin->max);
46
    rout->hasNonExist = rin->hasNonExist;
47
  }
48
  return rout;
49
}
50
51
NrrdRange *
52
nrrdRangeNix(NrrdRange *range) {
53
54
4
  airFree(range);
55
2
  return NULL;
56
}
57
58
void
59
nrrdRangeReset(NrrdRange *range) {
60
61
  if (range) {
62
    range->min = AIR_NAN;
63
    range->max = AIR_NAN;
64
    range->hasNonExist = nrrdHasNonExistUnknown;
65
  }
66
}
67
68
/*
69
** not using biff (obviously)
70
*/
71
void
72
nrrdRangeSet(NrrdRange *range, const Nrrd *nrrd, int blind8BitRange) {
73
4
  NRRD_TYPE_BIGGEST _min, _max;
74
  int blind;
75
76
2
  if (!range) {
77
    return;
78
  }
79
4
  if (nrrd
80
4
      && !airEnumValCheck(nrrdType, nrrd->type)
81
4
      && nrrdTypeBlock != nrrd->type) {
82
2
    blind = (nrrdBlind8BitRangeTrue == blind8BitRange
83
6
             || (nrrdBlind8BitRangeState == blind8BitRange
84
2
                 && nrrdStateBlind8BitRange));
85

4
    if (blind && 1 == nrrdTypeSize[nrrd->type]) {
86
      if (nrrdTypeChar == nrrd->type) {
87
        range->min = SCHAR_MIN;
88
        range->max = SCHAR_MAX;
89
      } else {
90
        range->min = 0;
91
        range->max = UCHAR_MAX;
92
      }
93
      range->hasNonExist = nrrdHasNonExistFalse;
94
    } else {
95
2
      nrrdMinMaxExactFind[nrrd->type](&_min, &_max, &(range->hasNonExist),
96
                                      nrrd);
97
2
      range->min = nrrdDLoad[nrrd->type](&_min);
98
2
      range->max = nrrdDLoad[nrrd->type](&_max);
99
    }
100
  } else {
101
    range->min = range->max = AIR_NAN;
102
    range->hasNonExist = nrrdHasNonExistUnknown;
103
  }
104
2
  return;
105
2
}
106
107
/*
108
******** nrrdRangePercentileSet
109
**
110
** this is called when information about the range of values in the
111
** nrrd is requested; and the learned information is put into "range"
112
** (overwriting whatever is there!)
113
**
114
** uses biff
115
*/
116
int
117
nrrdRangePercentileSet(NrrdRange *range, const Nrrd *nrrd,
118
                       double minPerc, double maxPerc,
119
                       unsigned int hbins, int blind8BitRange) {
120
  static const char me[]="nrrdRangePercentileSet";
121
  airArray *mop;
122
  Nrrd *nhist;
123
  double allmin, allmax, minval, maxval, *hist, sum, total, sumPerc;
124
  unsigned int hi;
125
126
  if (!(range && nrrd)) {
127
    biffAddf(NRRD, "%s: got NULL pointer", me);
128
    return 0;
129
  }
130
  nrrdRangeSet(range, nrrd, blind8BitRange);
131
  /* range->min and range->max
132
     are the full range of nrrd (except maybe for blind8) */
133
  if (!minPerc && !maxPerc) {
134
    /* wanted full range; there is nothing more to do */
135
    return 0;
136
  }
137
  if (!hbins) {
138
    biffAddf(NRRD, "%s: sorry, non-histogram-based percentiles not "
139
             "currently implemented (need hbins > 0)", me);
140
    return 1;
141
  }
142
  if (!(hbins >= 5)) {
143
    biffAddf(NRRD, "%s: # histogram bins %u unreasonably small", me, hbins);
144
    return 1;
145
  }
146
  if (range->hasNonExist) {
147
    biffAddf(NRRD, "%s: sorry, can currently do histogram-based percentiles "
148
             "only in arrays with no non-existent values", me);
149
    return 1;
150
  }
151
152
  mop = airMopNew();
153
  allmin = range->min;
154
  allmax = range->max;
155
156
  nhist = nrrdNew();
157
  airMopAdd(mop, nhist, (airMopper)nrrdNuke, airMopAlways);
158
  /* the histogram is over the entire range of values */
159
  if (nrrdHisto(nhist, nrrd, range, NULL, hbins, nrrdTypeDouble)) {
160
    biffAddf(NRRD, "%s: trouble making histogram", me);
161
    airMopError(mop);
162
    return 1;
163
  }
164
  hist = AIR_CAST(double *, nhist->data);
165
  total = AIR_CAST(double, nrrdElementNumber(nrrd));
166
  if (minPerc) {
167
    minval = AIR_NAN;
168
    sumPerc = AIR_ABS(minPerc)*total/100.0;
169
    sum = hist[0];
170
    for (hi=1; hi<hbins; hi++) {
171
      sum += hist[hi];
172
      if (sum >= sumPerc) {
173
        minval = AIR_AFFINE(0, hi-1, hbins-1,
174
                            nhist->axis[0].min, nhist->axis[0].max);
175
        break;
176
      }
177
    }
178
    if (hi == hbins || !AIR_EXISTS(minval)) {
179
      biffAddf(NRRD, "%s: failed to find lower %g-percentile value",
180
               me, minPerc);
181
      airMopError(mop);
182
      return 1;
183
    }
184
    range->min = (minPerc > 0
185
                  ? minval
186
                  : 2*allmin - minval);
187
  }
188
  if (maxPerc) {
189
    maxval = AIR_NAN;
190
    sumPerc = AIR_ABS(maxPerc)*total/100.0;
191
    sum = hist[hbins-1];
192
    for (hi=hbins-1; hi; hi--) {
193
      sum += hist[hi-1];
194
      if (sum >= sumPerc) {
195
        maxval = AIR_AFFINE(0, hi, hbins-1,
196
                            nhist->axis[0].min, nhist->axis[0].max);
197
        break;
198
      }
199
    }
200
    if (!hi || !AIR_EXISTS(maxval)) {
201
      biffAddf(NRRD, "%s: failed to find upper %g-percentile value", me,
202
               maxPerc);
203
      airMopError(mop);
204
      return 1;
205
    }
206
    range->max = (maxPerc > 0
207
                  ? maxval
208
                  : 2*allmax - maxval);
209
  }
210
211
  airMopOkay(mop);
212
  return 0;
213
}
214
215
/*
216
******** nrrdRangePercentileFromStringSet
217
**
218
** Implements smarts of figuring out a range from no info ("nan"),
219
** a known percentile (e.g. "3%") or a known explicit value (e.g. "3"),
220
** for both min and max.  Used by "unu quantize" and others.
221
*/
222
int
223
nrrdRangePercentileFromStringSet(NrrdRange *range, const Nrrd *nrrd,
224
                                 const char *_minStr, const char *_maxStr,
225
                                 unsigned int hbins, int blind8BitRange) {
226
  static const char me[]="nrrdRangePercentileFromStringSet";
227
  double minVal, maxVal, minPerc, maxPerc;
228
  char *minStr, *maxStr;
229
  unsigned int mmIdx;
230
  airArray *mop;
231
232
  if (!(range && nrrd && _minStr && _maxStr)) {
233
    biffAddf(NRRD, "%s: got NULL pointer", me);
234
    return 1;
235
  }
236
  mop = airMopNew();
237
  minStr = airStrdup(_minStr);
238
  airMopAdd(mop, minStr, airFree, airMopAlways);
239
  maxStr = airStrdup(_maxStr);
240
  airMopAdd(mop, maxStr, airFree, airMopAlways);
241
242
  /* parse min and max */
243
  minVal = maxVal = minPerc = maxPerc = AIR_NAN;
244
  for (mmIdx=0; mmIdx<=1; mmIdx++) {
245
    int percwant;
246
    double val, *mmv, *mmp;
247
    char *mmStr;
248
    if (0 == mmIdx) {
249
      mmv = &minVal;
250
      mmp = &minPerc;
251
      mmStr = minStr;
252
    } else {
253
      mmv = &maxVal;
254
      mmp = &maxPerc;
255
      mmStr = maxStr;
256
    }
257
    if (airEndsWith(mmStr, NRRD_MINMAX_PERC_SUFF)) {
258
      percwant = AIR_TRUE;
259
      mmStr[strlen(mmStr)-strlen(NRRD_MINMAX_PERC_SUFF)] = '\0';
260
    } else {
261
      percwant = AIR_FALSE;
262
    }
263
    if (1 != airSingleSscanf(mmStr, "%lf", &val)) {
264
      biffAddf(NRRD, "%s: couldn't parse \"%s\" for %s", me,
265
               !mmIdx ? _minStr : _maxStr,
266
               !mmIdx ? "minimum" : "maximum");
267
      airMopError(mop);
268
      return 1;
269
    }
270
    if (percwant) {
271
      if (!AIR_EXISTS(val)) {
272
        biffAddf(NRRD, "%s: %s percentile must exist", me,
273
                 !mmIdx ? "minimum" : "maximum");
274
        airMopError(mop);
275
        return 1;
276
      }
277
      /* setting a finite percentile */
278
      *mmp = val;
279
    } else if (!AIR_EXISTS(val)) {
280
      /* don't want a percentile, and given value is nan
281
         => same as full range == zero percentile */
282
      *mmp = 0.0;
283
    } else {
284
      /* value exists: given explicitly */
285
      *mmv = val;
286
    }
287
  }
288
  /* so whenever one end of the range is not given explicitly,
289
     it has been mapped to a statement about the percentile,
290
     which does require learning about the nrrd's values */
291
  if (AIR_EXISTS(minPerc) || AIR_EXISTS(maxPerc)) {
292
    if (nrrdRangePercentileSet(range, nrrd,
293
                               AIR_EXISTS(minPerc) ? minPerc : 0.0,
294
                               AIR_EXISTS(maxPerc) ? maxPerc : 0.0,
295
                               hbins, blind8BitRange)) {
296
      biffAddf(NRRD, "%s: trouble finding percentile range", me);
297
      airMopError(mop);
298
      return 1;
299
    }
300
  }
301
  /* whatever explicit values were given are now stored */
302
  if (AIR_EXISTS(minVal)) {
303
    range->min = minVal;
304
  }
305
  if (AIR_EXISTS(maxVal)) {
306
    range->max = maxVal;
307
  }
308
309
  airMopOkay(mop);
310
  return 0;
311
}
312
313
/*
314
** wrapper around nrrdRangeSet that (effectively) sets range->min
315
** and range->min only if they didn't already exist
316
*/
317
void
318
nrrdRangeSafeSet(NrrdRange *range, const Nrrd *nrrd, int blind8BitRange) {
319
  double minIn, maxIn;
320
321
  if (!range) {
322
    return;
323
  }
324
  minIn = range->min;
325
  maxIn = range->max;
326
  nrrdRangeSet(range, nrrd, blind8BitRange);
327
  if (AIR_EXISTS(minIn)) {
328
    range->min = minIn;
329
  }
330
  if (AIR_EXISTS(maxIn)) {
331
    range->max = maxIn;
332
  }
333
  return;
334
}
335
336
/*
337
** does not use biff
338
*/
339
NrrdRange *
340
nrrdRangeNewSet(const Nrrd *nrrd, int blind8BitRange) {
341
  NrrdRange *range;
342
343
4
  range = nrrdRangeNew(0, 0);  /* doesn't matter what values are used here */
344
2
  nrrdRangeSet(range, nrrd, blind8BitRange);
345
2
  return range;
346
}
347
348
/*
349
******** nrrdHasNonExist
350
**
351
** returns the nrrdHasNonExist* enum value appropriate for a given nrrd.
352
** By cleverness, this value can be used as a regular C boolean, so that
353
** the function will act as you expect.
354
**
355
** (the existence of this function implies that I'll never have an airEnum
356
** of the same name, which would be the usual thing to do with a C enum,
357
** but I don't think an airEnum for this would be useful)
358
*/
359
int
360
nrrdHasNonExist(const Nrrd *nrrd) {
361
2
  NRRD_TYPE_BIGGEST _min, _max;
362
1
  int ret;
363
364
2
  if (nrrd
365
2
      && !airEnumValCheck(nrrdType, nrrd->type)
366
2
      && nrrdTypeBlock != nrrd->type) {
367
1
    if (nrrdTypeIsIntegral[nrrd->type]) {
368
1
      ret = nrrdHasNonExistFalse;
369
1
    } else {
370
      /* HEY: this could be optimized by being more specialized */
371
      nrrdMinMaxExactFind[nrrd->type](&_min, &_max, &ret, nrrd);
372
    }
373
  } else {
374
    ret = nrrdHasNonExistUnknown;
375
  }
376
2
  return ret;
377
1
}