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 |
|
|
******** nrrdArithGamma() |
29 |
|
|
** |
30 |
|
|
** map the values in a nrrd through a power function; essentially: |
31 |
|
|
** val = pow(val, 1/gamma), but this is after the val has been normalized |
32 |
|
|
** to be in the range of 0.0 to 1.0 (assuming that the given min and |
33 |
|
|
** max really are the full range of the values in the nrrd). Thus, |
34 |
|
|
** the given min and max values are fixed points of this |
35 |
|
|
** transformation. Using a negative gamma means that after the pow() |
36 |
|
|
** function has been applied, the value is inverted with respect to |
37 |
|
|
** min and max (like in xv). |
38 |
|
|
*/ |
39 |
|
|
int |
40 |
|
|
nrrdArithGamma(Nrrd *nout, const Nrrd *nin, |
41 |
|
|
const NrrdRange *_range, double Gamma) { |
42 |
|
|
static const char me[]="nrrdArithGamma", func[]="gamma"; |
43 |
|
|
double val, min, max; |
44 |
|
|
size_t I, num; |
45 |
|
|
NrrdRange *range; |
46 |
|
|
airArray *mop; |
47 |
|
|
double (*lup)(const void *, size_t); |
48 |
|
|
double (*ins)(void *, size_t, double); |
49 |
|
|
|
50 |
|
|
if (!(nout && nin)) { |
51 |
|
|
/* _range can be NULL */ |
52 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
53 |
|
|
return 1; |
54 |
|
|
} |
55 |
|
|
if (!( AIR_EXISTS(Gamma) )) { |
56 |
|
|
biffAddf(NRRD, "%s: gamma doesn't exist", me); |
57 |
|
|
return 1; |
58 |
|
|
} |
59 |
|
|
if (!( nrrdTypeBlock != nin->type && nrrdTypeBlock != nout->type )) { |
60 |
|
|
biffAddf(NRRD, "%s: can't deal with %s type", me, |
61 |
|
|
airEnumStr(nrrdType, nrrdTypeBlock)); |
62 |
|
|
return 1; |
63 |
|
|
} |
64 |
|
|
if (nout != nin) { |
65 |
|
|
if (nrrdCopy(nout, nin)) { |
66 |
|
|
biffAddf(NRRD, "%s: couldn't initialize by copy to output", me); |
67 |
|
|
return 1; |
68 |
|
|
} |
69 |
|
|
} |
70 |
|
|
mop = airMopNew(); |
71 |
|
|
if (_range) { |
72 |
|
|
range = nrrdRangeCopy(_range); |
73 |
|
|
nrrdRangeSafeSet(range, nin, nrrdBlind8BitRangeState); |
74 |
|
|
} else { |
75 |
|
|
range = nrrdRangeNewSet(nin, nrrdBlind8BitRangeTrue); |
76 |
|
|
} |
77 |
|
|
airMopAdd(mop, range, (airMopper)nrrdRangeNix, airMopAlways); |
78 |
|
|
min = range->min; |
79 |
|
|
max = range->max; |
80 |
|
|
if (min == max) { |
81 |
|
|
/* this is stupid. We want min < max to avoid making NaNs */ |
82 |
|
|
max += 1; |
83 |
|
|
} |
84 |
|
|
lup = nrrdDLookup[nin->type]; |
85 |
|
|
ins = nrrdDInsert[nout->type]; |
86 |
|
|
Gamma = 1/Gamma; |
87 |
|
|
num = nrrdElementNumber(nin); |
88 |
|
|
if (Gamma < 0.0) { |
89 |
|
|
Gamma = -Gamma; |
90 |
|
|
for (I=0; I<num; I++) { |
91 |
|
|
val = lup(nin->data, I); |
92 |
|
|
val = AIR_AFFINE(min, val, max, 0.0, 1.0); |
93 |
|
|
val = pow(val, Gamma); |
94 |
|
|
val = AIR_AFFINE(1.0, val, 0.0, min, max); |
95 |
|
|
ins(nout->data, I, val); |
96 |
|
|
} |
97 |
|
|
} else { |
98 |
|
|
for (I=0; I<num; I++) { |
99 |
|
|
val = lup(nin->data, I); |
100 |
|
|
val = AIR_AFFINE(min, val, max, 0.0, 1.0); |
101 |
|
|
val = pow(val, Gamma); |
102 |
|
|
val = AIR_AFFINE(0.0, val, 1.0, min, max); |
103 |
|
|
ins(nout->data, I, val); |
104 |
|
|
} |
105 |
|
|
} |
106 |
|
|
if (nrrdContentSet_va(nout, func, nin, "%g,%g,%g", min, max, Gamma)) { |
107 |
|
|
biffAddf(NRRD, "%s:", me); |
108 |
|
|
airMopError(mop); return 1; |
109 |
|
|
} |
110 |
|
|
if (nout != nin) { |
111 |
|
|
nrrdAxisInfoCopy(nout, nin, NULL, NRRD_AXIS_INFO_NONE); |
112 |
|
|
} |
113 |
|
|
/* basic info handled by nrrdCopy above */ |
114 |
|
|
|
115 |
|
|
airMopOkay(mop); |
116 |
|
|
return 0; |
117 |
|
|
} |
118 |
|
|
|
119 |
|
|
/* ---------------------------- unary -------------- */ |
120 |
|
|
|
121 |
|
|
static double _nrrdUnaryOpNegative(double a) {return -a;} |
122 |
|
|
static double _nrrdUnaryOpReciprocal(double a) {return 1.0/a;} |
123 |
|
|
static double _nrrdUnaryOpSin(double a) {return sin(a);} |
124 |
|
|
static double _nrrdUnaryOpCos(double a) {return cos(a);} |
125 |
|
|
static double _nrrdUnaryOpTan(double a) {return tan(a);} |
126 |
|
|
static double _nrrdUnaryOpAsin(double a) {return asin(a);} |
127 |
|
|
static double _nrrdUnaryOpAcos(double a) {return acos(a);} |
128 |
|
|
static double _nrrdUnaryOpAtan(double a) {return atan(a);} |
129 |
|
|
static double _nrrdUnaryOpExp(double a) {return exp(a);} |
130 |
|
|
static double _nrrdUnaryOpLog(double a) {return log(a);} |
131 |
|
|
static double _nrrdUnaryOpLog2(double a) {return log(a)/0.69314718;} |
132 |
|
|
static double _nrrdUnaryOpLog10(double a) {return log10(a);} |
133 |
|
|
/* This code for log1p and expm1 comes from |
134 |
|
|
http://www.plunk.org/~hatch/rightway.php which in turn references |
135 |
|
|
http://www.cs.berkeley.edu/~wkahan/Math128/Sumnfp.pdf from the |
136 |
|
|
great Kahan of IEEE 754 fame, but sadly that URL no longer works |
137 |
|
|
(though the Math128 directory is still there, as are other documents) */ |
138 |
|
|
static double _nrrdUnaryOpLog1p(double a) { |
139 |
|
|
double b; |
140 |
|
|
|
141 |
|
|
b = 1.0 + a; |
142 |
|
|
if (b == 1.0) { |
143 |
|
|
/* "a" was so close to zero that we had underflow when adding to 1, |
144 |
|
|
resulting in something that is exactly equal to 1. So, use the |
145 |
|
|
first term of Taylor expansion of log(x+1) around 0 == x */ |
146 |
|
|
return a; |
147 |
|
|
} |
148 |
|
|
/* else "a" was not so near zero; but GLK doesn't fully grasp |
149 |
|
|
the design of this expression */ |
150 |
|
|
return log(b)*a/(b-1); |
151 |
|
|
} |
152 |
|
|
static double _nrrdUnaryOpExpm1(double x) { |
153 |
|
|
double u; |
154 |
|
|
|
155 |
|
|
u = exp(x); |
156 |
|
|
if (u == 1.0) { |
157 |
|
|
/* "x" was so close to 0.0 that exp return exactly 1; subtracting |
158 |
|
|
1 from this will give a constant for a range of x's. Instead, |
159 |
|
|
use the Taylor expansion of exp(x)-1 around 0 == x */ |
160 |
|
|
return x; |
161 |
|
|
} else if (u-1.0 == -1.0) { |
162 |
|
|
/* "x" was close enough to -inf that exp returned something so close |
163 |
|
|
to 0 that subtracting 1 resulted in exactly -1; return that */ |
164 |
|
|
return -1.0; |
165 |
|
|
} |
166 |
|
|
/* else "x" was neither near 0.0 or -inf, but GLK doesn't fully grasp |
167 |
|
|
the design of this expression */ |
168 |
|
|
return (u-1.0)*x/log(u); |
169 |
|
|
} |
170 |
|
|
static double _nrrdUnaryOpSqrt(double a) {return sqrt(a);} |
171 |
|
|
static double _nrrdUnaryOpCbrt(double a) {return airCbrt(a);} |
172 |
|
|
static double _nrrdUnaryOpErf(double a) {return airErf(a);} |
173 |
|
|
static double _nrrdUnaryOpNerf(double a) {return (1+airErf(a))/2;} |
174 |
|
|
static double _nrrdUnaryOpCeil(double a) {return ceil(a);} |
175 |
|
|
static double _nrrdUnaryOpFloor(double a) {return floor(a);} |
176 |
|
|
static double _nrrdUnaryOpRoundUp(double a) {return AIR_ROUNDUP(a);} |
177 |
|
|
static double _nrrdUnaryOpRoundDown(double a) {return AIR_ROUNDDOWN(a);} |
178 |
|
|
static double _nrrdUnaryOpAbs(double a) {return AIR_ABS(a);} |
179 |
|
|
static double _nrrdUnaryOpSgn(double a) { |
180 |
|
|
return (a < 0.0 ? -1 : (a > 0.0 ? 1 : 0));} |
181 |
|
|
static double _nrrdUnaryOpExists(double a) {return AIR_EXISTS(a);} |
182 |
|
|
static double _nrrdUnaryOpRand(double a) { |
183 |
|
|
AIR_UNUSED(a); |
184 |
|
|
return airDrandMT(); |
185 |
|
|
} |
186 |
|
|
static double _nrrdUnaryOpNormalRand(double a) { |
187 |
|
|
double v; |
188 |
|
|
AIR_UNUSED(a); |
189 |
|
|
airNormalRand(&v, NULL); |
190 |
|
|
return v; |
191 |
|
|
} |
192 |
|
|
static double _nrrdUnaryOpIf(double a) { return (a ? 1 : 0); } |
193 |
|
|
static double _nrrdUnaryOpZero(double a) { |
194 |
|
|
AIR_UNUSED(a); |
195 |
|
|
return 0.0; |
196 |
|
|
} |
197 |
|
|
static double _nrrdUnaryOpOne(double a) { |
198 |
|
|
AIR_UNUSED(a); |
199 |
|
|
return 1.0; |
200 |
|
|
} |
201 |
|
|
static double _nrrdUnaryOpTauOfSigma(double s) { return airTauOfSigma(s); } |
202 |
|
|
static double _nrrdUnaryOpSigmaOfTau(double t) { return airSigmaOfTau(t); } |
203 |
|
|
|
204 |
|
|
double (*_nrrdUnaryOp[NRRD_UNARY_OP_MAX+1])(double) = { |
205 |
|
|
NULL, |
206 |
|
|
_nrrdUnaryOpNegative, |
207 |
|
|
_nrrdUnaryOpReciprocal, |
208 |
|
|
_nrrdUnaryOpSin, |
209 |
|
|
_nrrdUnaryOpCos, |
210 |
|
|
_nrrdUnaryOpTan, |
211 |
|
|
_nrrdUnaryOpAsin, |
212 |
|
|
_nrrdUnaryOpAcos, |
213 |
|
|
_nrrdUnaryOpAtan, |
214 |
|
|
_nrrdUnaryOpExp, |
215 |
|
|
_nrrdUnaryOpLog, |
216 |
|
|
_nrrdUnaryOpLog2, |
217 |
|
|
_nrrdUnaryOpLog10, |
218 |
|
|
_nrrdUnaryOpLog1p, |
219 |
|
|
_nrrdUnaryOpExpm1, |
220 |
|
|
_nrrdUnaryOpSqrt, |
221 |
|
|
_nrrdUnaryOpCbrt, |
222 |
|
|
_nrrdUnaryOpErf, |
223 |
|
|
_nrrdUnaryOpNerf, |
224 |
|
|
_nrrdUnaryOpCeil, |
225 |
|
|
_nrrdUnaryOpFloor, |
226 |
|
|
_nrrdUnaryOpRoundUp, |
227 |
|
|
_nrrdUnaryOpRoundDown, |
228 |
|
|
_nrrdUnaryOpAbs, |
229 |
|
|
_nrrdUnaryOpSgn, |
230 |
|
|
_nrrdUnaryOpExists, |
231 |
|
|
_nrrdUnaryOpRand, |
232 |
|
|
_nrrdUnaryOpNormalRand, |
233 |
|
|
_nrrdUnaryOpIf, |
234 |
|
|
_nrrdUnaryOpZero, |
235 |
|
|
_nrrdUnaryOpOne, |
236 |
|
|
_nrrdUnaryOpTauOfSigma, |
237 |
|
|
_nrrdUnaryOpSigmaOfTau |
238 |
|
|
}; |
239 |
|
|
|
240 |
|
|
int |
241 |
|
|
nrrdArithUnaryOp(Nrrd *nout, int op, const Nrrd *nin) { |
242 |
|
|
static const char me[]="nrrdArithUnaryOp"; |
243 |
|
|
size_t N, I; |
244 |
|
|
int size[NRRD_DIM_MAX]; |
245 |
|
|
double (*insert)(void *v, size_t I, double d), |
246 |
|
|
(*lookup)(const void *v, size_t I), (*uop)(double), val; |
247 |
|
|
|
248 |
|
|
if (!(nout && nin)) { |
249 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
250 |
|
|
return 1; |
251 |
|
|
} |
252 |
|
|
if (nrrdTypeBlock == nin->type) { |
253 |
|
|
biffAddf(NRRD, "%s: can't operate on type %s", me, |
254 |
|
|
airEnumStr(nrrdType, nrrdTypeBlock)); |
255 |
|
|
return 1; |
256 |
|
|
} |
257 |
|
|
if (airEnumValCheck(nrrdUnaryOp, op)) { |
258 |
|
|
biffAddf(NRRD, "%s: unary op %d invalid", me, op); |
259 |
|
|
return 1; |
260 |
|
|
} |
261 |
|
|
if (nout != nin) { |
262 |
|
|
if (nrrdCopy(nout, nin)) { |
263 |
|
|
biffAddf(NRRD, "%s:", me); |
264 |
|
|
return 1; |
265 |
|
|
} |
266 |
|
|
} |
267 |
|
|
nrrdAxisInfoGet_nva(nin, nrrdAxisInfoSize, size); |
268 |
|
|
uop = _nrrdUnaryOp[op]; |
269 |
|
|
|
270 |
|
|
N = nrrdElementNumber(nin); |
271 |
|
|
lookup = nrrdDLookup[nin->type]; |
272 |
|
|
insert = nrrdDInsert[nin->type]; |
273 |
|
|
for (I=0; I<N; I++) { |
274 |
|
|
val = lookup(nin->data, I); |
275 |
|
|
insert(nout->data, I, uop(val)); |
276 |
|
|
} |
277 |
|
|
if (nrrdContentSet_va(nout, airEnumStr(nrrdUnaryOp, op), nin, "")) { |
278 |
|
|
biffAddf(NRRD, "%s:", me); |
279 |
|
|
return 1; |
280 |
|
|
} |
281 |
|
|
nrrdBasicInfoInit(nout, |
282 |
|
|
NRRD_BASIC_INFO_ALL ^ (NRRD_BASIC_INFO_OLDMIN_BIT |
283 |
|
|
| NRRD_BASIC_INFO_OLDMAX_BIT)); |
284 |
|
|
return 0; |
285 |
|
|
} |
286 |
|
|
|
287 |
|
|
/* ---------------------------- binary -------------- */ |
288 |
|
|
|
289 |
|
|
static double _nrrdBinaryOpAdd(double a, double b) {return a + b;} |
290 |
|
256200 |
static double _nrrdBinaryOpSubtract(double a, double b) {return a - b;} |
291 |
|
1556640 |
static double _nrrdBinaryOpMultiply(double a, double b) {return a * b;} |
292 |
|
|
static double _nrrdBinaryOpDivide(double a, double b) {return a / b;} |
293 |
|
|
static double _nrrdBinaryOpPow(double a, double b) {return pow(a,b);} |
294 |
|
|
static double _nrrdBinaryOpSgnPow(double a, double b) {return airSgnPow(a,b);} |
295 |
|
|
static double _nrrdBinaryOpFlippedSgnPow(double a, double b) {return airFlippedSgnPow(a,b);} |
296 |
|
|
static double _nrrdBinaryOpMod(double a, double b) { |
297 |
|
|
return AIR_MOD((int)a,(int)b);} |
298 |
|
|
static double _nrrdBinaryOpFmod(double a, double b) {return fmod(a,b);} |
299 |
|
|
static double _nrrdBinaryOpAtan2(double a, double b) {return atan2(a,b);} |
300 |
|
|
static double _nrrdBinaryOpMin(double a, double b) {return AIR_MIN(a,b);} |
301 |
|
|
static double _nrrdBinaryOpMax(double a, double b) {return AIR_MAX(a,b);} |
302 |
|
|
static double _nrrdBinaryOpLT(double a, double b) {return (a < b);} |
303 |
|
|
static double _nrrdBinaryOpLTE(double a, double b) {return (a <= b);} |
304 |
|
|
static double _nrrdBinaryOpGT(double a, double b) {return (a > b);} |
305 |
|
|
static double _nrrdBinaryOpGTE(double a, double b) {return (a >= b);} |
306 |
|
|
static double _nrrdBinaryOpCompare(double a, double b) { |
307 |
|
|
return (a < b ? -1 : (a > b ? 1 : 0));} |
308 |
|
|
static double _nrrdBinaryOpEqual(double a, double b) {return (a == b);} |
309 |
|
|
static double _nrrdBinaryOpNotEqual(double a, double b) {return (a != b);} |
310 |
|
|
static double _nrrdBinaryOpExists(double a, double b) {return (AIR_EXISTS(a) |
311 |
|
|
? a : b);} |
312 |
|
|
static double _nrrdBinaryOpIf(double a, double b) {return (a ? a : b);} |
313 |
|
|
static double _nrrdBinaryOpNormalRandScaleAdd(double a, double b) { |
314 |
|
10896480 |
double v; |
315 |
|
5448240 |
airNormalRand(&v, NULL); |
316 |
|
10896480 |
return a + b*v; |
317 |
|
5448240 |
} |
318 |
|
|
static double _nrrdBinaryOpRicianRand(double a, double b) { |
319 |
|
|
double vr, vi, rr, ri; |
320 |
|
|
airNormalRand(&rr, &ri); |
321 |
|
|
vr = a + b*rr; |
322 |
|
|
vi = b*ri; |
323 |
|
|
return sqrt(vr*vr + vi*vi); |
324 |
|
|
} |
325 |
|
|
|
326 |
|
|
|
327 |
|
|
double (*_nrrdBinaryOp[NRRD_BINARY_OP_MAX+1])(double, double) = { |
328 |
|
|
NULL, |
329 |
|
|
_nrrdBinaryOpAdd, |
330 |
|
|
_nrrdBinaryOpSubtract, |
331 |
|
|
_nrrdBinaryOpMultiply, |
332 |
|
|
_nrrdBinaryOpDivide, |
333 |
|
|
_nrrdBinaryOpPow, |
334 |
|
|
_nrrdBinaryOpSgnPow, |
335 |
|
|
_nrrdBinaryOpFlippedSgnPow, |
336 |
|
|
_nrrdBinaryOpMod, |
337 |
|
|
_nrrdBinaryOpFmod, |
338 |
|
|
_nrrdBinaryOpAtan2, |
339 |
|
|
_nrrdBinaryOpMin, |
340 |
|
|
_nrrdBinaryOpMax, |
341 |
|
|
_nrrdBinaryOpLT, |
342 |
|
|
_nrrdBinaryOpLTE, |
343 |
|
|
_nrrdBinaryOpGT, |
344 |
|
|
_nrrdBinaryOpGTE, |
345 |
|
|
_nrrdBinaryOpCompare, |
346 |
|
|
_nrrdBinaryOpEqual, |
347 |
|
|
_nrrdBinaryOpNotEqual, |
348 |
|
|
_nrrdBinaryOpExists, |
349 |
|
|
_nrrdBinaryOpIf, |
350 |
|
|
_nrrdBinaryOpNormalRandScaleAdd, |
351 |
|
|
_nrrdBinaryOpRicianRand, |
352 |
|
|
/* for these, the clamping is actually done by the caller */ |
353 |
|
|
_nrrdBinaryOpAdd, /* for nrrdBinaryOpAddClamp */ |
354 |
|
|
_nrrdBinaryOpSubtract, /* for nrrdBinaryOpSubtractClamp */ |
355 |
|
|
_nrrdBinaryOpMultiply, /* for nrrdBinaryOpMultiplyClamp */ |
356 |
|
|
}; |
357 |
|
|
|
358 |
|
|
/* |
359 |
|
|
******** nrrdArithBinaryOp |
360 |
|
|
** |
361 |
|
|
** this is a simplified version of nrrdArithIterBinaryOp, written after |
362 |
|
|
** that, in a hurry, to operate directly on two nrrds, instead with |
363 |
|
|
** the NrrdIter nonsense |
364 |
|
|
*/ |
365 |
|
|
int |
366 |
|
|
nrrdArithBinaryOp(Nrrd *nout, int op, const Nrrd *ninA, const Nrrd *ninB) { |
367 |
|
|
static const char me[]="nrrdArithBinaryOp"; |
368 |
|
|
char *contA, *contB; |
369 |
|
2 |
size_t N, I, size[NRRD_DIM_MAX]; |
370 |
|
|
double (*ins)(void *v, size_t I, double d), (*clmp)(double d), |
371 |
|
|
(*lupA)(const void *v, size_t I), (*lupB)(const void *v, size_t I), |
372 |
|
|
(*bop)(double a, double b), valA, valB; |
373 |
|
|
|
374 |
✓✗✓✗
✗✓ |
3 |
if (!( nout && !nrrdCheck(ninA) && !nrrdCheck(ninB) )) { |
375 |
|
|
biffAddf(NRRD, "%s: NULL pointer or invalid args", me); |
376 |
|
|
return 1; |
377 |
|
|
} |
378 |
✓✗✗✓
|
2 |
if (nrrdTypeBlock == ninA->type || nrrdTypeBlock == ninB->type) { |
379 |
|
|
biffAddf(NRRD, "%s: can't operate on type %s", me, |
380 |
|
|
airEnumStr(nrrdType, nrrdTypeBlock)); |
381 |
|
|
return 1; |
382 |
|
|
} |
383 |
✗✓ |
1 |
if (!nrrdSameSize(ninA, ninB, AIR_TRUE)) { |
384 |
|
|
biffAddf(NRRD, "%s: size mismatch between arguments", me); |
385 |
|
|
return 1; |
386 |
|
|
} |
387 |
✗✓ |
1 |
if (airEnumValCheck(nrrdBinaryOp, op)) { |
388 |
|
|
biffAddf(NRRD, "%s: binary op %d invalid", me, op); |
389 |
|
|
return 1; |
390 |
|
|
} |
391 |
|
|
|
392 |
|
1 |
nrrdAxisInfoGet_nva(ninA, nrrdAxisInfoSize, size); |
393 |
✓✗✓✗
|
2 |
if (!( nout == ninA || nout == ninB)) { |
394 |
✗✓ |
1 |
if (_nrrdMaybeAllocMaybeZero_nva(nout, ninA->type, ninA->dim, size, |
395 |
|
|
AIR_FALSE /* zero when no realloc */)) { |
396 |
|
|
biffAddf(NRRD, "%s: couldn't allocate output nrrd", me); |
397 |
|
|
return 1; |
398 |
|
|
} |
399 |
✗✓ |
1 |
if (nrrdAxisInfoCopy(nout, ninA, NULL, NRRD_AXIS_INFO_NONE)) { |
400 |
|
|
biffAddf(NRRD, "%s:", me); |
401 |
|
|
return 1; |
402 |
|
|
} |
403 |
|
1 |
nrrdBasicInfoCopy(nout, ninA, (NRRD_BASIC_INFO_DATA_BIT |
404 |
|
|
| NRRD_BASIC_INFO_TYPE_BIT |
405 |
|
|
| NRRD_BASIC_INFO_DIMENSION_BIT |
406 |
|
|
| NRRD_BASIC_INFO_CONTENT_BIT |
407 |
|
|
| NRRD_BASIC_INFO_COMMENTS_BIT |
408 |
|
1 |
| (nrrdStateKeyValuePairsPropagate |
409 |
|
|
? 0 |
410 |
|
|
: NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))); |
411 |
|
1 |
} |
412 |
|
1 |
nrrdBasicInfoInit(nout, |
413 |
|
|
NRRD_BASIC_INFO_ALL ^ (NRRD_BASIC_INFO_OLDMIN_BIT |
414 |
|
|
| NRRD_BASIC_INFO_OLDMAX_BIT)); |
415 |
|
1 |
bop = _nrrdBinaryOp[op]; |
416 |
|
|
|
417 |
|
1 |
N = nrrdElementNumber(ninA); |
418 |
|
1 |
lupA = nrrdDLookup[ninA->type]; |
419 |
|
1 |
lupB = nrrdDLookup[ninB->type]; |
420 |
|
1 |
ins = nrrdDInsert[nout->type]; |
421 |
|
1 |
if (nrrdBinaryOpAddClamp == op |
422 |
✗✓✗✓
|
2 |
|| nrrdBinaryOpSubtractClamp == op |
423 |
|
1 |
|| nrrdBinaryOpMultiplyClamp == op) { |
424 |
|
|
clmp = nrrdDClamp[nout->type]; |
425 |
|
|
} else { |
426 |
|
|
clmp = NULL; |
427 |
|
|
} |
428 |
✓✓ |
256202 |
for (I=0; I<N; I++) { |
429 |
|
|
double tmp; |
430 |
|
|
/* HEY: there is a loss of precision issue here with 64-bit ints */ |
431 |
|
128100 |
valA = lupA(ninA->data, I); |
432 |
|
128100 |
valB = lupB(ninB->data, I); |
433 |
|
128100 |
tmp = bop(valA, valB); |
434 |
✗✓ |
128100 |
if (clmp) { |
435 |
|
|
tmp = clmp(tmp); |
436 |
|
|
} |
437 |
|
128100 |
ins(nout->data, I, tmp); |
438 |
|
|
} |
439 |
|
|
|
440 |
|
1 |
contA = _nrrdContentGet(ninA); |
441 |
|
1 |
contB = _nrrdContentGet(ninB); |
442 |
✗✓ |
1 |
if (_nrrdContentSet_va(nout, airEnumStr(nrrdBinaryOp, op), |
443 |
|
|
contA, "%s", contB)) { |
444 |
|
|
biffAddf(NRRD, "%s:", me); |
445 |
|
|
free(contA); free(contB); return 1; |
446 |
|
|
} |
447 |
|
1 |
free(contA); |
448 |
|
1 |
free(contB); |
449 |
|
1 |
return 0; |
450 |
|
1 |
} |
451 |
|
|
|
452 |
|
|
int |
453 |
|
|
nrrdArithIterBinaryOpSelect(Nrrd *nout, int op, |
454 |
|
|
NrrdIter *inA, NrrdIter *inB, |
455 |
|
|
unsigned int which) { |
456 |
|
|
static const char me[]="nrrdArithIterBinaryOpSelect"; |
457 |
|
|
char *contA, *contB; |
458 |
|
32 |
size_t N, I, size[NRRD_DIM_MAX]; |
459 |
|
|
int type; |
460 |
|
|
double (*insert)(void *v, size_t I, double d), (*clmp)(double d), |
461 |
|
|
(*bop)(double a, double b), valA, valB; |
462 |
|
|
const Nrrd *nin; |
463 |
|
|
|
464 |
✗✓ |
16 |
if (!(nout && inA && inB)) { |
465 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
466 |
|
|
return 1; |
467 |
|
|
} |
468 |
✗✓ |
16 |
if (airEnumValCheck(nrrdBinaryOp, op)) { |
469 |
|
|
biffAddf(NRRD, "%s: binary op %d invalid", me, op); |
470 |
|
|
return 1; |
471 |
|
|
} |
472 |
✗✓ |
16 |
if (!( 0 == which || 1 == which )) { |
473 |
|
|
biffAddf(NRRD, "%s: which %u not 0 or 1", me, which); |
474 |
|
|
return 1; |
475 |
|
|
} |
476 |
✓✗ |
32 |
nin = (0 == which |
477 |
✓✗ |
32 |
? _NRRD_ITER_NRRD(inA) |
478 |
|
|
: _NRRD_ITER_NRRD(inB)); |
479 |
✗✓ |
16 |
if (!nin) { |
480 |
|
|
biffAddf(NRRD, "%s: selected input %u is a fixed value", me, which); |
481 |
|
|
return 1; |
482 |
|
|
} |
483 |
|
16 |
type = nin->type; |
484 |
|
16 |
nrrdAxisInfoGet_nva(nin, nrrdAxisInfoSize, size); |
485 |
|
|
|
486 |
✗✓ |
16 |
if (_nrrdMaybeAllocMaybeZero_nva(nout, type, nin->dim, size, |
487 |
|
|
AIR_FALSE /* zero when no realloc */)) { |
488 |
|
|
biffAddf(NRRD, "%s: couldn't allocate output nrrd", me); |
489 |
|
|
return 1; |
490 |
|
|
} |
491 |
|
16 |
nrrdBasicInfoCopy(nout, nin, (NRRD_BASIC_INFO_DATA_BIT |
492 |
|
|
| NRRD_BASIC_INFO_TYPE_BIT |
493 |
|
|
| NRRD_BASIC_INFO_DIMENSION_BIT |
494 |
|
|
| NRRD_BASIC_INFO_CONTENT_BIT |
495 |
|
|
| NRRD_BASIC_INFO_COMMENTS_BIT |
496 |
|
16 |
| (nrrdStateKeyValuePairsPropagate |
497 |
|
|
? 0 |
498 |
|
|
: NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))); |
499 |
|
16 |
nrrdBasicInfoInit(nout, |
500 |
|
|
NRRD_BASIC_INFO_ALL ^ (NRRD_BASIC_INFO_OLDMIN_BIT |
501 |
|
|
| NRRD_BASIC_INFO_OLDMAX_BIT)); |
502 |
|
16 |
bop = _nrrdBinaryOp[op]; |
503 |
|
|
|
504 |
|
|
/* |
505 |
|
|
fprintf(stderr, "%s: inA->left = %d, inB->left = %d\n", me, |
506 |
|
|
(int)(inA->left), (int)(inB->left)); |
507 |
|
|
*/ |
508 |
|
16 |
N = nrrdElementNumber(nin); |
509 |
|
16 |
insert = nrrdDInsert[type]; |
510 |
|
16 |
if (nrrdBinaryOpAddClamp == op |
511 |
✗✓✗✓
|
32 |
|| nrrdBinaryOpSubtractClamp == op |
512 |
|
16 |
|| nrrdBinaryOpMultiplyClamp == op) { |
513 |
|
|
clmp = nrrdDClamp[nout->type]; |
514 |
|
|
} else { |
515 |
|
|
clmp = NULL; |
516 |
|
|
} |
517 |
✓✓ |
12453152 |
for (I=0; I<N; I++) { |
518 |
|
|
double tmp; |
519 |
|
|
/* HEY: there is a loss of precision issue here with 64-bit ints */ |
520 |
|
6226560 |
valA = nrrdIterValue(inA); |
521 |
|
6226560 |
valB = nrrdIterValue(inB); |
522 |
|
6226560 |
tmp = bop(valA, valB); |
523 |
✗✓ |
6226560 |
if (clmp) { |
524 |
|
|
tmp = clmp(tmp); |
525 |
|
|
} |
526 |
|
6226560 |
insert(nout->data, I, tmp); |
527 |
|
|
} |
528 |
|
16 |
contA = nrrdIterContent(inA); |
529 |
|
16 |
contB = nrrdIterContent(inB); |
530 |
✗✓ |
16 |
if (_nrrdContentSet_va(nout, airEnumStr(nrrdBinaryOp, op), |
531 |
|
|
contA, "%s", contB)) { |
532 |
|
|
biffAddf(NRRD, "%s:", me); |
533 |
|
|
free(contA); free(contB); return 1; |
534 |
|
|
} |
535 |
✓✓ |
16 |
if (nout != nin) { |
536 |
|
8 |
nrrdAxisInfoCopy(nout, nin, NULL, NRRD_AXIS_INFO_NONE); |
537 |
|
8 |
} |
538 |
|
16 |
free(contA); |
539 |
|
16 |
free(contB); |
540 |
|
16 |
return 0; |
541 |
|
16 |
} |
542 |
|
|
|
543 |
|
|
int |
544 |
|
|
nrrdArithIterBinaryOp(Nrrd *nout, int op, NrrdIter *inA, NrrdIter *inB) { |
545 |
|
|
static const char me[]="nrrdArithIterBinaryOp"; |
546 |
|
|
unsigned int which; |
547 |
|
|
|
548 |
✗✓ |
32 |
if (!(nout && inA && inB)) { |
549 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
550 |
|
|
return 1; |
551 |
|
|
} |
552 |
✓✗✗✓
✗✗ |
48 |
which = (_NRRD_ITER_NRRD(inA) |
553 |
|
|
? 0 |
554 |
|
|
: (_NRRD_ITER_NRRD(inB) |
555 |
|
|
? 1 |
556 |
|
|
: 2)); |
557 |
✗✓ |
16 |
if (2 == which) { |
558 |
|
|
biffAddf(NRRD, "%s: can't operate on two fixed values", me); |
559 |
|
|
return 1; |
560 |
|
|
} |
561 |
✗✓ |
16 |
if (nrrdArithIterBinaryOpSelect(nout, op, inA, inB, which)) { |
562 |
|
|
biffAddf(NRRD, "%s: trouble", me); |
563 |
|
|
return 1; |
564 |
|
|
} |
565 |
|
16 |
return 0; |
566 |
|
16 |
} |
567 |
|
|
|
568 |
|
|
/* ---------------------------- ternary -------------- */ |
569 |
|
|
|
570 |
|
|
static double _nrrdTernaryOpAdd(double a, double b, double c) { |
571 |
|
|
return a + b + c; |
572 |
|
|
} |
573 |
|
|
static double _nrrdTernaryOpMultiply(double a, double b, double c) { |
574 |
|
|
return a * b * c; |
575 |
|
|
} |
576 |
|
|
static double _nrrdTernaryOpMin(double a, double b, double c) { |
577 |
|
|
b = AIR_MIN(b, c); |
578 |
|
|
return AIR_MIN(a, b); |
579 |
|
|
} |
580 |
|
|
/* |
581 |
|
|
** minsmooth(x, w, M) is like min(x,M), but starting at value M-w, values |
582 |
|
|
** are lowered (via erf), so that the output is asymptotic to M |
583 |
|
|
*/ |
584 |
|
|
static double _nrrdTernaryOpMinSmooth(double x, double width, double max) { |
585 |
|
|
double tran; |
586 |
|
|
tran = max - width; |
587 |
|
|
return (tran < max /* using the function as intended */ |
588 |
|
|
? (x < tran |
589 |
|
|
? x |
590 |
|
|
: airErf((x-tran)*0.886226925452758/(max - tran))*(max - tran) + tran) |
591 |
|
|
: AIR_MIN(x, max)); /* transition in wrong place; revert to simple max() */ |
592 |
|
|
} |
593 |
|
|
static double _nrrdTernaryOpMax(double a, double b, double c) { |
594 |
|
|
b = AIR_MAX(b, c); |
595 |
|
|
return AIR_MAX(a, b); |
596 |
|
|
} |
597 |
|
|
/* |
598 |
|
|
** maxsmooth(m, w, x) is like max(m,x), but starting at value m+w, values |
599 |
|
|
** are raised (via erf), so that the output is asymptotic to m |
600 |
|
|
*/ |
601 |
|
|
static double _nrrdTernaryOpMaxSmooth(double min, double width, double x) { |
602 |
|
|
double tran; |
603 |
|
|
tran = min + width; |
604 |
|
|
return (min < tran /* using the function as intended */ |
605 |
|
|
? (tran < x |
606 |
|
|
? x |
607 |
|
|
: airErf((x-tran)*0.886226925452758/(min - tran))*(min - tran) + tran) |
608 |
|
|
: AIR_MAX(x, min)); /* transition in wrong place; revert to simple max() */ |
609 |
|
|
} |
610 |
|
|
static double _nrrdTernaryOpLTSmooth(double a, double w, double b) { |
611 |
|
|
return AIR_AFFINE(-1.0, airErf((b-a)/w), 1.0, 0.0, 1.0); |
612 |
|
|
} |
613 |
|
|
static double _nrrdTernaryOpGTSmooth(double a, double w, double b) { |
614 |
|
|
return AIR_AFFINE(-1.0, airErf((a-b)/w), 1.0, 0.0, 1.0); |
615 |
|
|
} |
616 |
|
|
static double _nrrdTernaryOpClamp(double a, double b, double c) { |
617 |
|
|
return AIR_CLAMP(a, b, c); |
618 |
|
|
} |
619 |
|
|
static double _nrrdTernaryOpIfElse(double a, double b, double c) { |
620 |
|
|
return (a ? b : c); |
621 |
|
|
} |
622 |
|
|
static double _nrrdTernaryOpLerp(double a, double b, double c) { |
623 |
|
|
/* we do something more than the simple lerp here because |
624 |
|
|
we want to facilitate usage as something which can get around |
625 |
|
|
non-existent values (b and c as NaN or Inf) without |
626 |
|
|
getting polluted by them. */ |
627 |
|
|
|
628 |
|
|
if (0.0 == a) { |
629 |
|
|
return b; |
630 |
|
|
} else if (1.0 == a) { |
631 |
|
|
return c; |
632 |
|
|
} else { |
633 |
|
|
return AIR_LERP(a, b, c); |
634 |
|
|
} |
635 |
|
|
} |
636 |
|
|
static double _nrrdTernaryOpExists(double a, double b, double c) { |
637 |
|
|
return (AIR_EXISTS(a) ? b : c); |
638 |
|
|
} |
639 |
|
|
static double _nrrdTernaryOpInOpen(double a, double b, double c) { |
640 |
|
|
return (AIR_IN_OP(a, b, c)); |
641 |
|
|
} |
642 |
|
|
static double _nrrdTernaryOpInClosed(double a, double b, double c) { |
643 |
|
|
return (AIR_IN_CL(a, b, c)); |
644 |
|
|
} |
645 |
|
|
static double _nrrdTernaryOpGaussian(double x, double mu, double sig) { |
646 |
|
|
return airGaussian(x, mu, sig); |
647 |
|
|
} |
648 |
|
|
static double _nrrdTernaryOpRician(double x, double mu, double sig) { |
649 |
|
|
return airRician(x, mu, sig); |
650 |
|
|
} |
651 |
|
|
double (*_nrrdTernaryOp[NRRD_TERNARY_OP_MAX+1])(double, double, double) = { |
652 |
|
|
NULL, |
653 |
|
|
_nrrdTernaryOpAdd, |
654 |
|
|
_nrrdTernaryOpMultiply, |
655 |
|
|
_nrrdTernaryOpMin, |
656 |
|
|
_nrrdTernaryOpMinSmooth, |
657 |
|
|
_nrrdTernaryOpMax, |
658 |
|
|
_nrrdTernaryOpMaxSmooth, |
659 |
|
|
_nrrdTernaryOpLTSmooth, |
660 |
|
|
_nrrdTernaryOpGTSmooth, |
661 |
|
|
_nrrdTernaryOpClamp, |
662 |
|
|
_nrrdTernaryOpIfElse, |
663 |
|
|
_nrrdTernaryOpLerp, |
664 |
|
|
_nrrdTernaryOpExists, |
665 |
|
|
_nrrdTernaryOpInOpen, |
666 |
|
|
_nrrdTernaryOpInClosed, |
667 |
|
|
_nrrdTernaryOpGaussian, |
668 |
|
|
_nrrdTernaryOpRician |
669 |
|
|
}; |
670 |
|
|
|
671 |
|
|
/* |
672 |
|
|
******** nrrdArithTerneryOp |
673 |
|
|
** |
674 |
|
|
** HEY: UNTESTED UNTESTED UNTESTED UNTESTED UNTESTED UNTESTED UNTESTED |
675 |
|
|
** |
676 |
|
|
** this is a simplified version of nrrdArithIterTernaryOp, written after |
677 |
|
|
** that, in a hurry, to operate directly on three nrrds, instead with |
678 |
|
|
** the NrrdIter nonsense |
679 |
|
|
*/ |
680 |
|
|
int |
681 |
|
|
nrrdArithTernaryOp(Nrrd *nout, int op, const Nrrd *ninA, |
682 |
|
|
const Nrrd *ninB, const Nrrd *ninC) { |
683 |
|
|
static const char me[]="nrrdArithTernaryOp"; |
684 |
|
|
char *contA, *contB, *contC; |
685 |
|
|
size_t N, I, size[NRRD_DIM_MAX]; |
686 |
|
|
double (*ins)(void *v, size_t I, double d), |
687 |
|
|
(*lupA)(const void *v, size_t I), (*lupB)(const void *v, size_t I), |
688 |
|
|
(*lupC)(const void *v, size_t I), |
689 |
|
|
(*top)(double a, double b, double c), valA, valB, valC; |
690 |
|
|
|
691 |
|
|
if (!( nout && !nrrdCheck(ninA) && !nrrdCheck(ninB) && !nrrdCheck(ninC) )) { |
692 |
|
|
biffAddf(NRRD, "%s: NULL pointer or invalid args", me); |
693 |
|
|
return 1; |
694 |
|
|
} |
695 |
|
|
if (!( nrrdSameSize(ninA, ninB, AIR_TRUE) && |
696 |
|
|
nrrdSameSize(ninA, ninC, AIR_TRUE) )) { |
697 |
|
|
biffAddf(NRRD, "%s: size mismatch between arguments", me); |
698 |
|
|
return 1; |
699 |
|
|
} |
700 |
|
|
if (airEnumValCheck(nrrdTernaryOp, op)) { |
701 |
|
|
biffAddf(NRRD, "%s: ternary op %d invalid", me, op); |
702 |
|
|
return 1; |
703 |
|
|
} |
704 |
|
|
|
705 |
|
|
nrrdAxisInfoGet_nva(ninA, nrrdAxisInfoSize, size); |
706 |
|
|
if (!( nout == ninA || nout == ninB || nout == ninC)) { |
707 |
|
|
if (_nrrdMaybeAllocMaybeZero_nva(nout, ninA->type, ninA->dim, size, |
708 |
|
|
AIR_FALSE /* zero when no realloc */)) { |
709 |
|
|
biffAddf(NRRD, "%s: couldn't allocate output nrrd", me); |
710 |
|
|
return 1; |
711 |
|
|
} |
712 |
|
|
if (nrrdAxisInfoCopy(nout, ninA, NULL, NRRD_AXIS_INFO_NONE)) { |
713 |
|
|
biffAddf(NRRD, "%s:", me); |
714 |
|
|
return 1; |
715 |
|
|
} |
716 |
|
|
nrrdBasicInfoCopy(nout, ninA, (NRRD_BASIC_INFO_DATA_BIT |
717 |
|
|
| NRRD_BASIC_INFO_TYPE_BIT |
718 |
|
|
| NRRD_BASIC_INFO_DIMENSION_BIT |
719 |
|
|
| NRRD_BASIC_INFO_CONTENT_BIT |
720 |
|
|
| NRRD_BASIC_INFO_COMMENTS_BIT |
721 |
|
|
| (nrrdStateKeyValuePairsPropagate |
722 |
|
|
? 0 |
723 |
|
|
: NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))); |
724 |
|
|
} |
725 |
|
|
nrrdBasicInfoInit(nout, |
726 |
|
|
NRRD_BASIC_INFO_ALL ^ (NRRD_BASIC_INFO_OLDMIN_BIT |
727 |
|
|
| NRRD_BASIC_INFO_OLDMAX_BIT)); |
728 |
|
|
top = _nrrdTernaryOp[op]; |
729 |
|
|
|
730 |
|
|
N = nrrdElementNumber(ninA); |
731 |
|
|
lupA = nrrdDLookup[ninA->type]; |
732 |
|
|
lupB = nrrdDLookup[ninB->type]; |
733 |
|
|
lupC = nrrdDLookup[ninC->type]; |
734 |
|
|
ins = nrrdDInsert[nout->type]; |
735 |
|
|
for (I=0; I<N; I++) { |
736 |
|
|
/* HEY: there is a loss of precision issue here with 64-bit ints */ |
737 |
|
|
valA = lupA(ninA->data, I); |
738 |
|
|
valB = lupB(ninB->data, I); |
739 |
|
|
valC = lupC(ninC->data, I); |
740 |
|
|
ins(nout->data, I, top(valA, valB, valC)); |
741 |
|
|
} |
742 |
|
|
|
743 |
|
|
contA = _nrrdContentGet(ninA); |
744 |
|
|
contB = _nrrdContentGet(ninB); |
745 |
|
|
contC = _nrrdContentGet(ninC); |
746 |
|
|
if (_nrrdContentSet_va(nout, airEnumStr(nrrdTernaryOp, op), |
747 |
|
|
contA, "%s,%s", contB, contC)) { |
748 |
|
|
biffAddf(NRRD, "%s:", me); |
749 |
|
|
free(contA); free(contB); free(contC); return 1; |
750 |
|
|
} |
751 |
|
|
free(contA); |
752 |
|
|
free(contB); |
753 |
|
|
free(contC); |
754 |
|
|
|
755 |
|
|
return 0; |
756 |
|
|
} |
757 |
|
|
|
758 |
|
|
int |
759 |
|
|
nrrdArithIterTernaryOpSelect(Nrrd *nout, int op, |
760 |
|
|
NrrdIter *inA, NrrdIter *inB, NrrdIter *inC, |
761 |
|
|
unsigned int which) { |
762 |
|
|
static const char me[]="nrrdArithIterTernaryOpSelect"; |
763 |
|
|
char *contA, *contB, *contC; |
764 |
|
|
size_t N, I, size[NRRD_DIM_MAX]; |
765 |
|
|
int type; |
766 |
|
|
double (*insert)(void *v, size_t I, double d), |
767 |
|
|
(*top)(double a, double b, double c), valA, valB, valC; |
768 |
|
|
const Nrrd *nin; |
769 |
|
|
|
770 |
|
|
if (!(nout && inA && inB && inC)) { |
771 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
772 |
|
|
return 1; |
773 |
|
|
} |
774 |
|
|
if (airEnumValCheck(nrrdTernaryOp, op)) { |
775 |
|
|
biffAddf(NRRD, "%s: ternary op %d invalid", me, op); |
776 |
|
|
return 1; |
777 |
|
|
} |
778 |
|
|
if (!( 0 == which || 1 == which || 2 == which )) { |
779 |
|
|
biffAddf(NRRD, "%s: which %u not valid, want 0, 1, or 2", me, which); |
780 |
|
|
return 1; |
781 |
|
|
} |
782 |
|
|
nin = (0 == which |
783 |
|
|
? _NRRD_ITER_NRRD(inA) |
784 |
|
|
: (1 == which |
785 |
|
|
? _NRRD_ITER_NRRD(inB) |
786 |
|
|
: _NRRD_ITER_NRRD(inC))); |
787 |
|
|
if (!nin) { |
788 |
|
|
biffAddf(NRRD, "%s: selected input %u is a fixed value", me, which); |
789 |
|
|
return 1; |
790 |
|
|
} |
791 |
|
|
type = nin->type; |
792 |
|
|
nrrdAxisInfoGet_nva(nin, nrrdAxisInfoSize, size); |
793 |
|
|
if (_nrrdMaybeAllocMaybeZero_nva(nout, type, nin->dim, size, |
794 |
|
|
AIR_FALSE /* zero when no realloc */)) { |
795 |
|
|
biffAddf(NRRD, "%s: couldn't allocate output nrrd", me); |
796 |
|
|
return 1; |
797 |
|
|
} |
798 |
|
|
nrrdBasicInfoCopy(nout, nin, (NRRD_BASIC_INFO_DATA_BIT |
799 |
|
|
| NRRD_BASIC_INFO_TYPE_BIT |
800 |
|
|
| NRRD_BASIC_INFO_DIMENSION_BIT |
801 |
|
|
| NRRD_BASIC_INFO_CONTENT_BIT |
802 |
|
|
| NRRD_BASIC_INFO_COMMENTS_BIT |
803 |
|
|
| (nrrdStateKeyValuePairsPropagate |
804 |
|
|
? 0 |
805 |
|
|
: NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))); |
806 |
|
|
nrrdBasicInfoInit(nout, |
807 |
|
|
NRRD_BASIC_INFO_ALL ^ (NRRD_BASIC_INFO_OLDMIN_BIT |
808 |
|
|
| NRRD_BASIC_INFO_OLDMAX_BIT)); |
809 |
|
|
top = _nrrdTernaryOp[op]; |
810 |
|
|
|
811 |
|
|
/* |
812 |
|
|
fprintf(stderr, "%!s: inA->left = %d, inB->left = %d\n", me, |
813 |
|
|
(int)(inA->left), (int)(inB->left)); |
814 |
|
|
*/ |
815 |
|
|
N = nrrdElementNumber(nin); |
816 |
|
|
insert = nrrdDInsert[type]; |
817 |
|
|
for (I=0; I<N; I++) { |
818 |
|
|
/* HEY: there is a loss of precision issue here with 64-bit ints */ |
819 |
|
|
valA = nrrdIterValue(inA); |
820 |
|
|
valB = nrrdIterValue(inB); |
821 |
|
|
valC = nrrdIterValue(inC); |
822 |
|
|
/* |
823 |
|
|
if (!(I % 1000)) { |
824 |
|
|
fprintf(stderr, "!%s: %d: top(%g,%g,%g) = %g\n", me, (int)I, |
825 |
|
|
valA, valB, valC, |
826 |
|
|
top(valA, valB, valC)); |
827 |
|
|
} |
828 |
|
|
*/ |
829 |
|
|
insert(nout->data, I, top(valA, valB, valC)); |
830 |
|
|
} |
831 |
|
|
contA = nrrdIterContent(inA); |
832 |
|
|
contB = nrrdIterContent(inB); |
833 |
|
|
contC = nrrdIterContent(inC); |
834 |
|
|
if (_nrrdContentSet_va(nout, airEnumStr(nrrdTernaryOp, op), |
835 |
|
|
contA, "%s,%s", contB, contC)) { |
836 |
|
|
biffAddf(NRRD, "%s:", me); |
837 |
|
|
free(contA); free(contB); free(contC); return 1; |
838 |
|
|
} |
839 |
|
|
if (nout != nin) { |
840 |
|
|
nrrdAxisInfoCopy(nout, nin, NULL, NRRD_AXIS_INFO_NONE); |
841 |
|
|
} |
842 |
|
|
free(contA); |
843 |
|
|
free(contB); |
844 |
|
|
free(contC); |
845 |
|
|
return 0; |
846 |
|
|
} |
847 |
|
|
|
848 |
|
|
int |
849 |
|
|
nrrdArithIterTernaryOp(Nrrd *nout, int op, |
850 |
|
|
NrrdIter *inA, NrrdIter *inB, NrrdIter *inC) { |
851 |
|
|
static const char me[]="nrrdArithIterTernaryOp"; |
852 |
|
|
unsigned int which; |
853 |
|
|
|
854 |
|
|
if (!(nout && inA && inB && inC)) { |
855 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
856 |
|
|
return 1; |
857 |
|
|
} |
858 |
|
|
which = (_NRRD_ITER_NRRD(inA) |
859 |
|
|
? 0 |
860 |
|
|
: (_NRRD_ITER_NRRD(inB) |
861 |
|
|
? 1 |
862 |
|
|
: (_NRRD_ITER_NRRD(inC) |
863 |
|
|
? 2 |
864 |
|
|
: 3 ))); |
865 |
|
|
if (3 == which) { |
866 |
|
|
biffAddf(NRRD, "%s: can't operate on 3 fixed values", me); |
867 |
|
|
return 1; |
868 |
|
|
} |
869 |
|
|
if (nrrdArithIterTernaryOpSelect(nout, op, inA, inB, inC, which)) { |
870 |
|
|
biffAddf(NRRD, "%s: trouble", me); |
871 |
|
|
return 1; |
872 |
|
|
} |
873 |
|
|
return 0; |
874 |
|
|
} |
875 |
|
|
|
876 |
|
|
int |
877 |
|
|
nrrdArithAffine(Nrrd *nout, double minIn, |
878 |
|
|
const Nrrd *nin, double maxIn, |
879 |
|
|
double minOut, double maxOut, int clamp) { |
880 |
|
|
static const char me[]="nrrdArithAffine"; |
881 |
|
|
size_t I, N; |
882 |
|
|
double (*ins)(void *v, size_t I, double d), |
883 |
|
|
(*lup)(const void *v, size_t I), mmin, mmax; |
884 |
|
|
|
885 |
|
|
if ( !nout || nrrdCheck(nin) ) { |
886 |
|
|
biffAddf(NRRD, "%s: got NULL pointer or invalid input", me); |
887 |
|
|
return 1; |
888 |
|
|
} |
889 |
|
|
if (nout != nin) { |
890 |
|
|
if (nrrdCopy(nout, nin)) { |
891 |
|
|
biffAddf(NRRD, "%s: couldn't initialize output", me); |
892 |
|
|
return 1; |
893 |
|
|
} |
894 |
|
|
} |
895 |
|
|
N = nrrdElementNumber(nin); |
896 |
|
|
ins = nrrdDInsert[nout->type]; |
897 |
|
|
lup = nrrdDLookup[nin->type]; |
898 |
|
|
mmin = AIR_MIN(minOut, maxOut); |
899 |
|
|
mmax = AIR_MAX(minOut, maxOut); |
900 |
|
|
for (I=0; I<N; I++) { |
901 |
|
|
double val; |
902 |
|
|
val = lup(nin->data, I); |
903 |
|
|
val = AIR_AFFINE(minIn, val, maxIn, minOut, maxOut); |
904 |
|
|
if (clamp) { |
905 |
|
|
val = AIR_CLAMP(mmin, val, mmax); |
906 |
|
|
} |
907 |
|
|
ins(nout->data, I, val); |
908 |
|
|
} |
909 |
|
|
/* HEY: it would be much better if the ordering here was the same as in |
910 |
|
|
AIR_AFFINE, but that's not easy with the way the content functions are |
911 |
|
|
now set up */ |
912 |
|
|
if (nrrdContentSet_va(nout, "affine", nin, |
913 |
|
|
"%g,%g,%g,%g", minIn, maxIn, |
914 |
|
|
minOut, maxOut)) { |
915 |
|
|
biffAddf(NRRD, "%s:", me); |
916 |
|
|
} |
917 |
|
|
return 0; |
918 |
|
|
} |
919 |
|
|
|
920 |
|
|
int |
921 |
|
|
nrrdArithIterAffine(Nrrd *nout, NrrdIter *minIn, |
922 |
|
|
NrrdIter *in, NrrdIter *maxIn, |
923 |
|
|
NrrdIter *minOut, NrrdIter *maxOut, int clamp) { |
924 |
|
|
static const char me[]="nrrdArithInterAffine"; |
925 |
|
|
double (*ins)(void *v, size_t I, double d), |
926 |
|
|
mini, vin, maxi, mino, maxo, vout; |
927 |
|
|
const Nrrd *nin; |
928 |
|
|
char *contA, *contB, *contC, *contD, *contE; |
929 |
|
|
size_t I, N; |
930 |
|
|
|
931 |
|
|
if (!(nout && minIn && in && maxIn && minOut && maxOut)) { |
932 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
933 |
|
|
return 1; |
934 |
|
|
} |
935 |
|
|
nin = (_NRRD_ITER_NRRD(in) |
936 |
|
|
? _NRRD_ITER_NRRD(in) |
937 |
|
|
: (_NRRD_ITER_NRRD(minIn) |
938 |
|
|
? _NRRD_ITER_NRRD(minIn) |
939 |
|
|
: (_NRRD_ITER_NRRD(maxIn) |
940 |
|
|
? _NRRD_ITER_NRRD(maxIn) |
941 |
|
|
: (_NRRD_ITER_NRRD(minOut) |
942 |
|
|
? _NRRD_ITER_NRRD(minOut) |
943 |
|
|
: _NRRD_ITER_NRRD(maxOut))))); |
944 |
|
|
if (!nin) { |
945 |
|
|
biffAddf(NRRD, "%s: can't operate solely on fixed values", me); |
946 |
|
|
return 1; |
947 |
|
|
} |
948 |
|
|
if (nrrdCopy(nout, nin)) { |
949 |
|
|
biffAddf(NRRD, "%s: couldn't initialize output", me); |
950 |
|
|
return 1; |
951 |
|
|
} |
952 |
|
|
N = nrrdElementNumber(nin); |
953 |
|
|
ins = nrrdDInsert[nout->type]; |
954 |
|
|
for (I=0; I<N; I++) { |
955 |
|
|
mini = nrrdIterValue(minIn); |
956 |
|
|
vin = nrrdIterValue(in); |
957 |
|
|
maxi = nrrdIterValue(maxIn); |
958 |
|
|
mino = nrrdIterValue(minOut); |
959 |
|
|
maxo = nrrdIterValue(maxOut); |
960 |
|
|
vout = AIR_AFFINE(mini, vin, maxi, mino, maxo); |
961 |
|
|
if (clamp) { |
962 |
|
|
double mmin = AIR_MIN(mino, maxo); |
963 |
|
|
double mmax = AIR_MAX(mino, maxo); |
964 |
|
|
vout = AIR_CLAMP(mmin, vout, mmax); |
965 |
|
|
} |
966 |
|
|
ins(nout->data, I, vout); |
967 |
|
|
} |
968 |
|
|
contA = nrrdIterContent(in); |
969 |
|
|
contB = nrrdIterContent(minIn); |
970 |
|
|
contC = nrrdIterContent(maxIn); |
971 |
|
|
contD = nrrdIterContent(maxOut); |
972 |
|
|
contE = nrrdIterContent(maxOut); |
973 |
|
|
/* HEY: same annoyance about order of arguments as in function above */ |
974 |
|
|
if (_nrrdContentSet_va(nout, "affine", contA, "%s,%s,%s,%s", |
975 |
|
|
contB, contC, contD, contE)) { |
976 |
|
|
biffAddf(NRRD, "%s:", me); |
977 |
|
|
free(contA); free(contB); free(contC); free(contD); free(contE); return 1; |
978 |
|
|
} |
979 |
|
|
free(contA); free(contB); free(contC); free(contD); free(contE); |
980 |
|
|
|
981 |
|
|
return 0; |
982 |
|
|
} |
983 |
|
|
|
984 |
|
|
unsigned int |
985 |
|
|
nrrdCRC32(const Nrrd *nin, int endian) { |
986 |
|
|
size_t nn; |
987 |
|
|
|
988 |
|
|
/* NULL nrrd or data */ |
989 |
✗✓ |
36 |
if (!nin |
990 |
✓✗ |
24 |
|| !(nin->data) |
991 |
✓✗ |
24 |
|| !(nn = nrrdElementSize(nin)*nrrdElementNumber(nin)) |
992 |
✓✗ |
24 |
|| airEnumValCheck(airEndian, endian)) { |
993 |
|
|
return 0; |
994 |
|
|
} |
995 |
|
|
|
996 |
|
24 |
return airCRC32(AIR_CAST(const unsigned char *, nin->data), |
997 |
|
12 |
nn, nrrdElementSize(nin), |
998 |
|
12 |
endian == airMyEndian() ? AIR_FALSE : AIR_TRUE); |
999 |
|
12 |
} |
1000 |
|
|
|