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 */ |