GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: src/air/miscAir.c Lines: 123 220 55.9 %
Date: 2017-05-26 Branches: 43 145 29.7 %

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 "air.h"
25
#include "privateAir.h"
26
/* timer functions */
27
#ifdef _WIN32
28
#include <io.h>
29
#include <fcntl.h>
30
#include <time.h>
31
#else
32
#include <sys/time.h>
33
#endif
34
35
/*
36
******** airTeemVersion
37
******** airTeemReleaseDone
38
******** airTeemReleaseDate
39
**
40
** updated with each release to contain a string representation of
41
** the Teem version number and release date.  Originated in version 1.5;
42
** use of TEEM_VERSION #defines started in 1.9
43
*/
44
const char *
45
airTeemVersion = TEEM_VERSION_STRING;
46
const int
47
airTeemReleaseDone = AIR_FALSE;
48
const char *
49
airTeemReleaseDate = "maybe 2014 or 2015";
50
51
/*
52
******** airTeemVersionSprint
53
**
54
** uniform way of printing information about the Teem version
55
*/
56
void
57
airTeemVersionSprint(char buff[AIR_STRLEN_LARGE]) {
58
4
  sprintf(buff, "Teem version %s, %s%s%s",
59
          airTeemVersion,
60
          airTeemReleaseDone ? "released on " : "",
61
          airTeemReleaseDate,
62
          airTeemReleaseDone ? "" : " (not yet released)");
63
2
  return;
64
}
65
66
double
67
_airSanityHelper(double val) {
68
120
  return val*val*val;
69
}
70
71
/*
72
******** airNull()
73
**
74
** returns NULL
75
*/
76
void *
77
airNull(void) {
78
79
298
  return NULL;
80
}
81
82
/*
83
******** airSetNull
84
**
85
** dereferences and sets to NULL, returns NULL
86
*/
87
void *
88
airSetNull(void **ptrP) {
89
90
2260
  if (ptrP) {
91
1130
    *ptrP = NULL;
92
1130
  }
93
1130
  return NULL;
94
}
95
96
/*
97
******** airFree()
98
**
99
** to facilitate setting a newly free()'d pointer; always returns NULL.
100
** also makes sure that NULL is not passed to free().
101
*/
102
void *
103
airFree(void *ptr) {
104
105
257210
  if (ptr) {
106
19438
    free(ptr);
107
19438
  }
108
128605
  return NULL;
109
}
110
111
/*
112
******** airFopen()
113
**
114
** encapsulates that idea that "-" is either standard in or stardard
115
** out, and does McRosopht stuff required to make piping work
116
**
117
** Does not error checking.  If fopen fails, then C' errno and strerror are
118
** left untouched for the caller to access.
119
*/
120
FILE *
121
airFopen(const char *name, FILE *std, const char *mode) {
122
  FILE *ret;
123
124
86
  if (!strcmp(name, "-")) {
125
    ret = std;
126
#ifdef _WIN32
127
    if (strchr(mode, 'b')) {
128
      _setmode(_fileno(ret), _O_BINARY);
129
    }
130
#endif
131
  } else {
132
43
    ret = fopen(name, mode);
133
  }
134
43
  return ret;
135
}
136
137
/*
138
******** airFclose()
139
**
140
** just to facilitate setting a newly fclose()'d file pointer to NULL
141
** also makes sure that NULL is not passed to fclose(), and won't close
142
** stdin, stdout, or stderr (its up to the user to open these correctly)
143
*/
144
FILE *
145
airFclose(FILE *file) {
146
147
110
  if (file) {
148

165
    if (!( stdin == file || stdout == file || stderr == file )) {
149
55
      fclose(file);
150
55
    }
151
  }
152
55
  return NULL;
153
}
154
155
/*
156
******** airSinglePrintf
157
**
158
** a complete stand-in for {f|s}printf(), as long as the given format
159
** string contains exactly one conversion sequence.  The utility of
160
** this is to standardize the printing of IEEE 754 special values:
161
** QNAN, SNAN -> "NaN"
162
** POS_INF -> "+inf"
163
** NEG_INF -> "-inf"
164
** The format string can contain other things besides just the
165
** conversion sequence: airSingleFprintf(f, " (%f)\n", AIR_NAN)
166
** will be the same as fprintf(f, " (%s)\n", "NaN");
167
**
168
** To get fprintf behavior, pass "str" as NULL
169
** to get sprintf bahavior, pass "file" as NULL
170
**
171
** Finding a complete {f|s|}printf replacement is a priority for Teem 2.0
172
*/
173
int
174
airSinglePrintf(FILE *file, char *str, const char *_fmt, ...) {
175
404
  char *fmt, buff[AIR_STRLEN_LARGE];
176
202
  double val=0, gVal, fVal;
177
  int ret, isF, isD, cls;
178
  char *conv=NULL, *p0, *p1, *p2, *p3, *p4, *p5;
179
202
  va_list ap;
180
181
202
  va_start(ap, _fmt);
182
202
  fmt = airStrdup(_fmt);
183
184
  /* this is needlessly complicated; the "l" modifier is a no-op */
185
202
  p0 = strstr(fmt, "%e");
186
202
  p1 = strstr(fmt, "%f");
187
202
  p2 = strstr(fmt, "%g");
188
202
  p3 = strstr(fmt, "%le");
189
202
  p4 = strstr(fmt, "%lf");
190
202
  p5 = strstr(fmt, "%lg");
191
606
  isF = p0 || p1 || p2;
192
606
  isD = p3 || p4 || p5;
193
  /* the code here says "isF" and "isD" as if it means "is float" or
194
     "is double".  It really should be "is2" or "is3", as in,
195
     "is 2-character conv. seq., or "is 3-character conv. seq." */
196
202
  if (isF) {
197
48
    conv = p0 ? p0 : (p1 ? p1 : p2);
198
16
  }
199
202
  if (isD) {
200
    conv = p3 ? p3 : (p4 ? p4 : p5);
201
  }
202
202
  if (isF || isD) {
203
    /* use "double" instead of "float" because var args are _always_
204
       subject to old-style C type promotions: float promotes to double */
205
48
    val = va_arg(ap, double);
206
16
    cls = airFPClass_d(val);
207
16
    switch (cls) {
208
    case airFP_SNAN:
209
    case airFP_QNAN:
210
    case airFP_POS_INF:
211
    case airFP_NEG_INF:
212
      if (isF) {
213
        memcpy(conv, "%s", 2);
214
      } else {
215
        /* this sneakiness allows us to replace a 3-character conversion
216
           sequence for a double (such as %lg) with a 3-character conversion
217
           for a string, which we know has at most 4 characters */
218
        memcpy(conv, "%4s", 3);
219
      }
220
      break;
221
    }
222
#define PRINT(F, S, C, V) ((F) ? fprintf((F),(C),(V)) : sprintf((S),(C),(V)))
223

16
    switch (cls) {
224
    case airFP_SNAN:
225
    case airFP_QNAN:
226
      ret = PRINT(file, str, fmt, "NaN");
227
      break;
228
    case airFP_POS_INF:
229
      ret = PRINT(file, str, fmt, "+inf");
230
      break;
231
    case airFP_NEG_INF:
232
      ret = PRINT(file, str, fmt, "-inf");
233
      break;
234
    default:
235
16
      if (p2 || p5) {
236
        /* got "%g" or "%lg", see if it would be better to use "%f" */
237
16
        sprintf(buff, "%f", val);
238
16
        sscanf(buff, "%lf", &fVal);
239
16
        sprintf(buff, "%g", val);
240
16
        sscanf(buff, "%lf", &gVal);
241
16
        if (fVal != gVal) {
242
          /* using %g (or %lg) lost precision!! Use %f (or %lf) instead */
243
          if (p2) {
244
            memcpy(conv, "%f", 2);
245
          } else {
246
            memcpy(conv, "%lf", 3);
247
          }
248
        }
249
      }
250
48
      ret = PRINT(file, str, fmt, val);
251
16
      break;
252
    }
253
  } else {
254
    /* conversion sequence is neither for float nor double */
255
558
    ret = file ? vfprintf(file, fmt, ap) : vsprintf(str, fmt, ap);
256
  }
257
258
202
  va_end(ap);
259
202
  free(fmt);
260
202
  return ret;
261
202
}
262
263
/*
264
******** airSprintSize_t
265
**
266
** sprints a single size_t to a given string, side-stepping
267
** non-standardized format specifier confusion with printf
268
*/
269
char *
270
airSprintSize_t(char _str[AIR_STRLEN_SMALL], size_t val) {
271
5323474
  char str[AIR_STRLEN_SMALL];
272
  unsigned int si;
273
274
2661737
  if (!_str) {
275
    return NULL;
276
  }
277
  si = AIR_STRLEN_SMALL-1;
278
2661737
  str[si] = '\0';
279
2661737
  do {
280
4324309
    str[--si] = AIR_CAST(char, (val % 10) + '0');
281
4324309
    val /= 10;
282
4324309
  } while (val);
283
2661737
  strcpy(_str, str + si);
284
2661737
  return _str;
285
2661737
}
286
287
/*
288
******** airSprintPtrdiff_t
289
**
290
** sprints a single ptrdiff_t to a given string, side-stepping
291
** non-standardized format specifier confusion with printf
292
*/
293
char *
294
airSprintPtrdiff_t(char _str[AIR_STRLEN_SMALL], ptrdiff_t val) {
295
322
  char str[AIR_STRLEN_SMALL];
296
  unsigned int si;
297
  int sign;
298
299
161
  if (!_str) {
300
    return NULL;
301
  }
302
  si = AIR_STRLEN_SMALL-1;
303
161
  str[si] = '\0';
304
161
  sign = (val < 0 ? -1 : 1);
305
161
  do {
306
    ptrdiff_t dig;
307
217
    dig = val % 10;
308
217
    str[--si] = AIR_CAST(char, dig > 0 ? dig + '0' : -dig + '0');
309
217
    val /= 10;
310
217
  } while (val);
311
161
  if (-1 == sign) {
312
17
    str[--si] = '-';
313
17
  }
314
161
  strcpy(_str, str + si);
315
161
  return _str;
316
161
}
317
318
/* ---- BEGIN non-NrrdIO */
319
320
const int
321
airPresent = 42;
322
323
/*
324
** sprints a length-"len" vector "vec" of size_t values into "dst", which is
325
** simply assumed to be big enough to hold this.  Vector enclosed in "[]" and
326
** values separated by ","
327
*/
328
char *
329
airSprintVecSize_t(char *dst, const size_t *vec, unsigned int len) {
330
1348758
  char stmp[AIR_STRLEN_SMALL];
331
  unsigned int axi;
332
333
  /* if we got NULL to sprint into, there's nothing to do but return it */
334
674379
  if (!dst) {
335
    return dst;
336
  }
337
674379
  strcpy(dst, "[");
338
6671958
  for (axi=0; axi<len; axi++) {
339
2661600
    if (axi) {
340
1987224
      strcat(dst, ",");
341
1987224
    }
342
2661600
    airSprintSize_t(stmp, vec[axi]);
343
2661600
    strcat(dst, stmp);
344
  }
345
674379
  strcat(dst, "]");
346
674379
  return dst;
347
674379
}
348
349
/*
350
******** airPrettySprintSize_t
351
**
352
** sprints a single size_t in a way that is human readable as
353
** bytes, kilobytes (KB), etc
354
*/
355
char *
356
airPrettySprintSize_t(char str[AIR_STRLEN_SMALL], size_t val) {
357
  static const char suff[][7] = {"bytes", "KB", "MB", "GB", "TB", "PB", "EB"};
358
  unsigned int suffIdx, suffNum;
359
  double dval;
360
361
42
  if (!str) {
362
    return NULL;
363
  }
364
  suffIdx = 0;
365
21
  dval = AIR_CAST(double, val);
366
  suffNum = AIR_UINT(sizeof(suff)/sizeof(suff[0]));
367

252
  while (suffIdx < suffNum-1) {  /* while we can go to a bigger suffix */
368
84
    if (dval > 1024) {
369
63
      dval /= 1024;
370
63
      suffIdx++;
371
    } else {
372
      break;
373
    }
374
  }
375
21
  sprintf(str, "%g %s", dval, suff[suffIdx]);
376
21
  return str;
377
21
}
378
379
/*
380
******** airStderr, airStdout, airStdin
381
**
382
** these exist only to give uniform access to FILE *s that might be
383
** annoying to access in higher-level language wrappings around Teem.
384
*/
385
FILE *
386
airStderr(void) {
387
2
  return stderr;
388
}
389
390
FILE *
391
airStdout(void) {
392
2
  return stdout;
393
}
394
395
FILE *
396
airStdin(void) {
397
2
  return stdin;
398
}
399
400
unsigned int
401
airBitsSet(unsigned int vv) {
402
  /* http://graphics.stanford.edu/~seander/bithacks.html#CountBitsSetKernighan */
403
  unsigned int cc;
404
  for (cc=0; vv; cc++) {
405
    /* wherever lowest bit is on in vv; vv-1 will have it off, and leave
406
       unchanged all the higher bits */
407
    vv &= vv - 1;
408
  }
409
  return cc;
410
}
411
412
/*
413
******** AIR_INDEX(i,x,I,L,t)
414
**
415
** READ CAREFULLY!!
416
**
417
** Utility for mapping a floating point x in given range [i,I] to the
418
** index of an array with L elements, AND SAVES THE INDEX INTO GIVEN
419
** VARIABLE t, WHICH MUST BE OF SOME INTEGER TYPE because this relies
420
** on the implicit cast of an assignment to truncate away the
421
** fractional part.  ALSO, t must be of a type large enough to hold
422
** ONE GREATER than L.  So you can't pass a variable of type unsigned
423
** char if L is 256
424
**
425
** DOES NOT DO BOUNDS CHECKING: given an x which is not inside [i,I],
426
** this may produce an index not inside [0,L-1] (but it won't always
427
** do so: the output being outside range [0,L-1] is not a reliable
428
** test of the input being outside range [i, I]).  The mapping is
429
** accomplished by dividing the range from i to I into L intervals,
430
** all but the last of which is half-open; the last one is closed.
431
** For example, the number line from 0 to 3 would be divided as
432
** follows for a call with i = 0, I = 4, L = 4:
433
**
434
** index:       0    1    2    3 = L-1
435
** intervals: [   )[   )[   )[    ]
436
**            |----|----|----|----|
437
** value:     0    1    2    3    4
438
**
439
** The main point of the diagram above is to show how I made the
440
** arbitrary decision to orient the half-open interval, and which
441
** end has the closed interval.
442
**
443
** Note that AIR_INDEX(0,3,4,4,t) and AIR_INDEX(0,4,4,4,t) both set t = 3
444
**
445
** The reason that this macro requires a argument for saving the
446
** result is that this is the easiest way to avoid extra conditionals.
447
** Otherwise, we'd have to do some check to see if x is close enough
448
** to I so that the generated index would be L and not L-1.  "Close
449
** enough" because due to precision problems you can have an x < I
450
** such that (x-i)/(I-i) == 1, which was a bug with the previous version
451
** of this macro.  It is far simpler to just do the index generation
452
** and then do the sneaky check to see if the index is too large by 1.
453
** We are relying on the fact that C _defines_ boolean true to be exactly 1.
454
**
455
** Note also that we are never explicity casting to one kind of int or
456
** another-- the given t can be any integral type, including long long.
457
*/
458
/*
459
#define AIR_INDEX(i,x,I,L,t) ( \
460
(t) = (L) * ((double)(x)-(i)) / ((double)(I)-(i)), \
461
(t) -= ((t) == (L)) )
462
*/
463
/*
464
******* airIndex
465
**
466
** Given a value range [min,max], a single input value val inside the range,
467
** and a number of intervals N spanning and evenly dividing that range,
468
** airIndex returns the index of the interval containing val:
469
**
470
**  input value range: min                     max
471
**   val somewhere in:  |-----|-----|-----|-----|
472
**        N intervals:  [    )[    )[    )[     ]  (for N=4)
473
**       output index:     0     1     2     3
474
**
475
** This can be used (as in nrrdHisto and other histogramming functions) to
476
** represent the range [min,max] with N samples between 0 and N-1.  Those
477
** samples are cell-centered, because "0" is logically located (in the
478
** continuous input range) in the *middle* of the first of the N intervals.
479
** In contrast, the *endpoints* of the N intervals (the 5 "|" in the picture
480
** above) form N+1 *node*-centered samples from min to max.
481
**
482
** If max < min, then "min" and "max" would be swapped in the diagram
483
** above and their roles are switched. This overdue fix was only added
484
** in version 1.11.
485
**
486
** NOTE: This does not do bounds checking; for that use airIndexClamp
487
*/
488
unsigned int
489
airIndex(double min, double val, double max, unsigned int N) {
490
  unsigned int idx;
491
  double mnm;
492
493
85228
  mnm = max - min;
494
42614
  if (mnm > 0) {
495
42614
    idx = AIR_UINT(N*(val - min)/mnm);
496
42614
    idx -= (idx == N);
497
42614
  } else if (mnm < 0) {
498
    idx = AIR_UINT(N*(val - max)/(-mnm));
499
    idx -= (idx == N);
500
    idx = N-1-idx;
501
  } else {
502
    idx = 0;
503
  }
504
42614
  return idx;
505
}
506
507
unsigned int
508
airIndexClamp(double min, double val, double max, unsigned int N) {
509
  unsigned int idx;
510
  double mnm;
511
512
259824
  mnm = max - min;
513
129912
  if (mnm > 0) {
514
129912
    val = AIR_MAX(min, val);
515
129912
    idx = AIR_UINT(N*(val - min)/mnm);
516
129912
    idx = AIR_MIN(idx, N-1);
517
129912
  } else if (mnm < 0) {
518
    val = AIR_MAX(max, val);
519
    idx = AIR_UINT(N*(val - max)/(-mnm));
520
    idx = AIR_MIN(idx, N-1);
521
    idx = N-1-idx;
522
  } else {
523
    idx = 0;
524
  }
525
129912
  return idx;
526
}
527
528
airULLong
529
airIndexULL(double min, double val, double max, airULLong N) {
530
  airULLong idx;
531
  if (min == max) {
532
    return 0;
533
  }
534
  if (min > max) {
535
    return N-1-airIndexULL(max, val, min, N);
536
  }
537
#if defined(_WIN32) && !defined(__CYGWIN__) && !defined(__MINGW32__)
538
  /* compile error on Win32-vs60: "error C2520: conversion from
539
     unsigned __int64 to double not implemented, use signed __int64 */
540
  airLLong sidx;
541
  sidx = AIR_CAST(airLLong, AIR_CAST(double, N)*(val - min)/(max - min));
542
  idx = AIR_CAST(airULLong, sidx);
543
#else
544
  idx = AIR_CAST(airULLong, AIR_CAST(double, N)*(val - min)/(max - min));
545
#endif
546
  idx -= (idx == N);
547
  return idx;
548
}
549
550
airULLong
551
airIndexClampULL(double min, double val, double max, airULLong N) {
552
  airULLong idx;
553
  if (min == max) {
554
    return 0;
555
  }
556
  if (min > max) {
557
    return N-1-airIndexULL(max, val, min, N);
558
  }
559
#if defined(_WIN32) && !defined(__CYGWIN__) && !defined(__MINGW32__)
560
  airLLong sidx;
561
  val = AIR_MAX(min, val); /* see note in airIndexClamp */
562
  sidx = AIR_CAST(airLLong, AIR_CAST(double, N)*(val - min)/(max - min));
563
  idx = AIR_CAST(airULLong, sidx);
564
#else
565
  val = AIR_MAX(min, val); /* see note in airIndexClamp */
566
  idx = AIR_CAST(airULLong, AIR_CAST(double, N)*(val - min)/(max - min));
567
#endif
568
  idx = AIR_MIN(idx, N-1);
569
  return idx;
570
}
571
572
/*
573
******* airDoneStr()
574
**
575
** dinky little utility for generating progress messages of the form
576
** "  1.9%" or " 35.3%" or  "100.0%"
577
**
578
** The message will ALWAYS be six characters, and will ALWAYS be
579
** preceeded by six backspaces.  Thus, you pass in a string to print
580
** into, and it had better be allocated for at least 6+6+1 = 13 chars.
581
*/
582
char *
583
airDoneStr(double start, double here, double end, char *str) {
584
  int perc=0;
585
586
  if (str) {
587
    if (end != start)
588
      perc = (int)(1000*(here-start)/(end-start) + 0.5);
589
    else
590
      perc = 1000;
591
    if (perc < 0) {
592
      sprintf(str, "\b\b\b\b\b\b ---- ");
593
    } else if (perc < 1000) {
594
      sprintf(str, "\b\b\b\b\b\b% 3d.%d%%", perc/10, perc%10);
595
    }
596
    else if (perc == 1000) {
597
      /* the "% 3d" formatting sequence should have taken care
598
         of this, but whatever */
599
      sprintf(str, "\b\b\b\b\b\b100.0%%");
600
    }
601
    else {
602
      sprintf(str, "\b\b\b\b\b\b done.");
603
    }
604
  }
605
606
  return str;
607
}
608
609
/*
610
******** airTime()
611
**
612
** returns current time in seconds (with millisecond resolution only
613
** when not on Windows) as a double.  From "man gettimeofday": The
614
** time is expressed in seconds and microseconds since midnight (0
615
** hour), January 1, 1970.
616
*/
617
double
618
airTime() {
619
#ifdef _WIN32
620
  /* HEY: this has crummy precision */
621
  return (double)clock()/CLOCKS_PER_SEC;
622
#else
623
  double sec, usec;
624
49218
  struct timeval tv;
625
626
24609
  gettimeofday(&tv, NULL);
627
24609
  sec = AIR_CAST(double, tv.tv_sec);
628
24609
  usec = AIR_CAST(double, tv.tv_usec);
629
49218
  return sec + usec*1.0e-6;
630
#endif
631
24609
}
632
633
const char
634
airTypeStr[AIR_TYPE_MAX+1][AIR_STRLEN_SMALL] = {
635
  "(unknown)",
636
  "bool",
637
  "int",
638
  "unsigned int",
639
  "long int",
640
  "unsigned long int",
641
  "size_t",
642
  "float",
643
  "double",
644
  "char",
645
  "string",
646
  "enum",
647
  "other",
648
};
649
650
const size_t
651
airTypeSize[AIR_TYPE_MAX+1] = {
652
  0,
653
  sizeof(int),
654
  sizeof(int),
655
  sizeof(unsigned int),
656
  sizeof(long int),
657
  sizeof(unsigned long int),
658
  sizeof(size_t),
659
  sizeof(float),
660
  sizeof(double),
661
  sizeof(char),
662
  sizeof(char*),
663
  sizeof(int),
664
  0   /* we don't know anything about type "other" */
665
};
666
667
/*
668
******** airEqvSettle()
669
**
670
** takes a mapping map[i], i in [0..len-1], and shifts the range of the
671
** mapping downward so that the range is a contiguous set of uints
672
** starting at 0.
673
**
674
** returns the number of different uints; previous version returned one
675
** less than this (the highest mapping output value, after the settling)
676
**
677
** honestly this doesn't necessarily have anything to do with processing
678
** equivalence classes, but its an operation that is nice to have in
679
** those cases
680
*/
681
unsigned int
682
airEqvSettle(unsigned int *map, unsigned int len) {
683
  unsigned int i, j, count, max, *hit;
684
685
  max = 0;
686
  for (i=0; i<len; i++) {
687
    max = AIR_MAX(max, map[i]);
688
  }
689
  hit = (unsigned int *)calloc(1+max, sizeof(unsigned int));
690
  for (i=0; i<len; i++) {
691
    hit[map[i]] = 1;
692
  }
693
  count = 0;
694
  for (j=0; j<=max; j++) {
695
    if (hit[j]) {
696
      hit[j] = count;
697
      count += 1;
698
    }
699
  }
700
  for (i=0; i<len; i++) {
701
    map[i] = hit[map[i]];
702
  }
703
  free(hit);
704
  return count;
705
}
706
707
/*
708
******** airEqvMap
709
**
710
** takes the equivalence pairs in eqvArr, and an array of uints map of
711
** length len, and puts in map[i] the uint that class i's value should
712
** be changed to.
713
**
714
** based on numerical recipes, C edition, pg. 346
715
** modifications:
716
**  - when resolving ancestors, map to the one with the lower index.
717
**  - applies settling to resulting map
718
**
719
** returns the number of different class IDs
720
*/
721
unsigned int
722
airEqvMap(airArray *eqvArr, unsigned int *map, unsigned int len) {
723
  unsigned int *eqv, j, k, t, eqi;
724
725
  for (j=0; j<len; j++) {
726
    map[j] = j;
727
  }
728
  eqv = (unsigned int*)(eqvArr->data);
729
  for (eqi=0; eqi<eqvArr->len; eqi++) {
730
    j = eqv[0 + 2*eqi];
731
    k = eqv[1 + 2*eqi];
732
    while (map[j] != j) {
733
      j = map[j];
734
    }
735
    while (map[k] != k) {
736
      k = map[k];
737
    }
738
    if (j != k) {
739
      if (j < k) {
740
        t = j; j = k; k = t;
741
      }
742
      map[j] = k;
743
    }
744
  }
745
  for (j=0; j<len; j++) {
746
    while (map[j] != map[map[j]]) {
747
      map[j] = map[map[j]];
748
    }
749
  }
750
  return airEqvSettle(map, len);
751
}
752
753
/*
754
******** airEqvAdd
755
**
756
** adds another equivalence (which may or may not amount to adding
757
** a new class; that will be determined later)
758
*/
759
void
760
airEqvAdd(airArray *eqvArr, unsigned int j, unsigned int k) {
761
  unsigned int *eqv, eqi;
762
763
  /* HEY: would it speed things up at all to enforce j < k? */
764
  if (eqvArr->len) {
765
    /* we have some equivalences, but we're only going to check against
766
       the last one in an effort to reduce duplicate entries */
767
    eqv = AIR_CAST(unsigned int*, eqvArr->data);
768
    eqi = eqvArr->len-1;
769
    if ( (eqv[0 + 2*eqi] == j && eqv[1 + 2*eqi] == k) ||
770
         (eqv[0 + 2*eqi] == k && eqv[1 + 2*eqi] == j) ) {
771
      /* don't add a duplicate */
772
      return;
773
    }
774
  }
775
  eqi = airArrayLenIncr(eqvArr, 1);
776
  eqv = AIR_CAST(unsigned int*, eqvArr->data);
777
  eqv[0 + 2*eqi] = j;
778
  eqv[1 + 2*eqi] = k;
779
  return;
780
}
781
782
/* ---- END non-NrrdIO */