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 |
|
|
|
25 |
|
|
#include "ten.h" |
26 |
|
|
#include "privateTen.h" |
27 |
|
|
|
28 |
|
|
typedef struct { |
29 |
|
|
double weight[3], amount, target; /* tenSizeNormalize */ |
30 |
|
|
/* amount: tenSizeScale */ |
31 |
|
|
double scale; int fixDet; int makePositive; /* tenAnisoScale */ |
32 |
|
|
double min, max; /* tenEigenvalueClamp */ |
33 |
|
|
double expo; /* tenEigenvaluePower */ |
34 |
|
|
double val; /* tenEigenvalueAdd */ |
35 |
|
|
} funcParm; |
36 |
|
|
|
37 |
|
|
enum { |
38 |
|
|
funcUnknown, |
39 |
|
|
funcSizeNormalize, |
40 |
|
|
funcSizeScale, |
41 |
|
|
funcAnisoScale, |
42 |
|
|
funcEigenvalueClamp, |
43 |
|
|
funcEigenvaluePower, |
44 |
|
|
funcEigenvalueAdd, |
45 |
|
|
funcEigenvalueMultiply, |
46 |
|
|
funcLog, |
47 |
|
|
funcExp, |
48 |
|
|
funcLast |
49 |
|
|
}; |
50 |
|
|
|
51 |
|
|
static int |
52 |
|
|
theFunc(Nrrd *nout, const Nrrd *nin, int func, funcParm *parm) { |
53 |
|
|
static const char me[]="theFunc"; |
54 |
|
|
float *tin, *tout, eval[3], evec[9], weight[3], size, mean; |
55 |
|
|
size_t NN, II; |
56 |
|
|
unsigned int ri; |
57 |
|
|
|
58 |
|
|
if (!AIR_IN_OP(funcUnknown, func, funcLast)) { |
59 |
|
|
biffAddf(TEN, "%s: given func %d out of range [%d,%d]", me, func, |
60 |
|
|
funcUnknown+1, funcLast-1); |
61 |
|
|
return 1; |
62 |
|
|
} |
63 |
|
|
if (!(nout && nin && parm)) { |
64 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
65 |
|
|
return 1; |
66 |
|
|
} |
67 |
|
|
if (tenTensorCheck(nin, nrrdTypeFloat, AIR_FALSE, AIR_TRUE)) { |
68 |
|
|
biffAddf(TEN, "%s: didn't get a tensor nrrd", me); |
69 |
|
|
return 1; |
70 |
|
|
} |
71 |
|
|
if (nout != nin) { |
72 |
|
|
if (nrrdCopy(nout, nin)) { |
73 |
|
|
biffMovef(TEN, NRRD, "%s: couldn't allocate output", me); |
74 |
|
|
return 1; |
75 |
|
|
} |
76 |
|
|
} |
77 |
|
|
|
78 |
|
|
tin = (float*)(nin->data); |
79 |
|
|
tout = (float*)(nout->data); |
80 |
|
|
NN = nrrdElementNumber(nin)/7; |
81 |
|
|
switch(func) { |
82 |
|
|
case funcSizeNormalize: |
83 |
|
|
ELL_3V_COPY_TT(weight, float, parm->weight); |
84 |
|
|
size = weight[0] + weight[1] + weight[2]; |
85 |
|
|
if (!size) { |
86 |
|
|
biffAddf(TEN, "%s: some of eigenvalue weights is zero", me); |
87 |
|
|
return 1; |
88 |
|
|
} |
89 |
|
|
weight[0] /= size; |
90 |
|
|
weight[1] /= size; |
91 |
|
|
weight[2] /= size; |
92 |
|
|
for (II=0; II<=NN-1; II++) { |
93 |
|
|
tenEigensolve_f(eval, evec, tin); |
94 |
|
|
size = (weight[0]*AIR_ABS(eval[0]) |
95 |
|
|
+ weight[1]*AIR_ABS(eval[1]) |
96 |
|
|
+ weight[2]*AIR_ABS(eval[2])); |
97 |
|
|
ELL_3V_SET_TT(eval, float, |
98 |
|
|
AIR_AFFINE(0, parm->amount, 1, |
99 |
|
|
eval[0], parm->target*eval[0]/size), |
100 |
|
|
AIR_AFFINE(0, parm->amount, 1, |
101 |
|
|
eval[1], parm->target*eval[1]/size), |
102 |
|
|
AIR_AFFINE(0, parm->amount, 1, |
103 |
|
|
eval[2], parm->target*eval[2]/size)); |
104 |
|
|
tenMakeSingle_f(tout, tin[0], eval, evec); |
105 |
|
|
tin += 7; |
106 |
|
|
tout += 7; |
107 |
|
|
} |
108 |
|
|
break; |
109 |
|
|
case funcSizeScale: |
110 |
|
|
for (II=0; II<=NN-1; II++) { |
111 |
|
|
TEN_T_SET_TT(tout, float, |
112 |
|
|
tin[0], |
113 |
|
|
parm->amount*tin[1], |
114 |
|
|
parm->amount*tin[2], |
115 |
|
|
parm->amount*tin[3], |
116 |
|
|
parm->amount*tin[4], |
117 |
|
|
parm->amount*tin[5], |
118 |
|
|
parm->amount*tin[6]); |
119 |
|
|
tin += 7; |
120 |
|
|
tout += 7; |
121 |
|
|
} |
122 |
|
|
break; |
123 |
|
|
case funcAnisoScale: |
124 |
|
|
for (II=0; II<=NN-1; II++) { |
125 |
|
|
tenEigensolve_f(eval, evec, tin); |
126 |
|
|
if (parm->fixDet) { |
127 |
|
|
eval[0] = AIR_MAX(eval[0], 0.00001f); |
128 |
|
|
eval[1] = AIR_MAX(eval[1], 0.00001f); |
129 |
|
|
eval[2] = AIR_MAX(eval[2], 0.00001f); |
130 |
|
|
ELL_3V_SET_TT(eval, float, log(eval[0]), log(eval[1]), log(eval[2])); |
131 |
|
|
} |
132 |
|
|
mean = (eval[0] + eval[1] + eval[2])/3.0f; |
133 |
|
|
ELL_3V_SET_TT(eval, float, |
134 |
|
|
AIR_LERP(parm->scale, mean, eval[0]), |
135 |
|
|
AIR_LERP(parm->scale, mean, eval[1]), |
136 |
|
|
AIR_LERP(parm->scale, mean, eval[2])); |
137 |
|
|
if (parm->fixDet) { |
138 |
|
|
ELL_3V_SET_TT(eval, float, exp(eval[0]), exp(eval[1]), exp(eval[2])); |
139 |
|
|
} |
140 |
|
|
if (eval[2] < 0 && parm->makePositive) { |
141 |
|
|
eval[0] = AIR_MAX(eval[0], 0.0f); |
142 |
|
|
eval[1] = AIR_MAX(eval[1], 0.0f); |
143 |
|
|
eval[2] = AIR_MAX(eval[2], 0.0f); |
144 |
|
|
} |
145 |
|
|
tenMakeSingle_f(tout, tin[0], eval, evec); |
146 |
|
|
tin += 7; |
147 |
|
|
tout += 7; |
148 |
|
|
} |
149 |
|
|
break; |
150 |
|
|
case funcEigenvalueClamp: |
151 |
|
|
for (II=0; II<=NN-1; II++) { |
152 |
|
|
tenEigensolve_f(eval, evec, tin); |
153 |
|
|
if (AIR_EXISTS(parm->min)) { |
154 |
|
|
ELL_3V_SET_TT(eval, float, |
155 |
|
|
AIR_MAX(eval[0], parm->min), |
156 |
|
|
AIR_MAX(eval[1], parm->min), |
157 |
|
|
AIR_MAX(eval[2], parm->min)); |
158 |
|
|
} |
159 |
|
|
if (AIR_EXISTS(parm->max)) { |
160 |
|
|
ELL_3V_SET_TT(eval, float, |
161 |
|
|
AIR_MIN(eval[0], parm->max), |
162 |
|
|
AIR_MIN(eval[1], parm->max), |
163 |
|
|
AIR_MIN(eval[2], parm->max)); |
164 |
|
|
} |
165 |
|
|
tenMakeSingle_f(tout, tin[0], eval, evec); |
166 |
|
|
tin += 7; |
167 |
|
|
tout += 7; |
168 |
|
|
} |
169 |
|
|
break; |
170 |
|
|
case funcEigenvaluePower: |
171 |
|
|
for (II=0; II<=NN-1; II++) { |
172 |
|
|
tenEigensolve_f(eval, evec, tin); |
173 |
|
|
ELL_3V_SET_TT(eval, float, |
174 |
|
|
pow(eval[0], parm->expo), |
175 |
|
|
pow(eval[1], parm->expo), |
176 |
|
|
pow(eval[2], parm->expo)); |
177 |
|
|
tenMakeSingle_f(tout, tin[0], eval, evec); |
178 |
|
|
tin += 7; |
179 |
|
|
tout += 7; |
180 |
|
|
} |
181 |
|
|
break; |
182 |
|
|
case funcEigenvalueAdd: |
183 |
|
|
for (II=0; II<=NN-1; II++) { |
184 |
|
|
/* HEY: this doesn't require eigensolve */ |
185 |
|
|
tenEigensolve_f(eval, evec, tin); |
186 |
|
|
ELL_3V_SET_TT(eval, float, |
187 |
|
|
eval[0] + parm->val, |
188 |
|
|
eval[1] + parm->val, |
189 |
|
|
eval[2] + parm->val); |
190 |
|
|
tenMakeSingle_f(tout, tin[0], eval, evec); |
191 |
|
|
tin += 7; |
192 |
|
|
tout += 7; |
193 |
|
|
} |
194 |
|
|
break; |
195 |
|
|
case funcEigenvalueMultiply: |
196 |
|
|
for (II=0; II<=NN-1; II++) { |
197 |
|
|
/* HEY: this doesn't require eigensolve */ |
198 |
|
|
tenEigensolve_f(eval, evec, tin); |
199 |
|
|
ELL_3V_SET_TT(eval, float, |
200 |
|
|
eval[0]*parm->val, |
201 |
|
|
eval[1]*parm->val, |
202 |
|
|
eval[2]*parm->val); |
203 |
|
|
tenMakeSingle_f(tout, tin[0], eval, evec); |
204 |
|
|
tin += 7; |
205 |
|
|
tout += 7; |
206 |
|
|
} |
207 |
|
|
break; |
208 |
|
|
case funcLog: |
209 |
|
|
for (II=0; II<=NN-1; II++) { |
210 |
|
|
tenEigensolve_f(eval, evec, tin); |
211 |
|
|
for (ri=0; ri<3; ri++) { |
212 |
|
|
eval[ri] = AIR_CAST(float, log(eval[ri])); |
213 |
|
|
eval[ri] = AIR_EXISTS(eval[ri]) ? eval[ri] : -1000000; |
214 |
|
|
} |
215 |
|
|
tenMakeSingle_f(tout, tin[0], eval, evec); |
216 |
|
|
tin += 7; |
217 |
|
|
tout += 7; |
218 |
|
|
} |
219 |
|
|
break; |
220 |
|
|
case funcExp: |
221 |
|
|
for (II=0; II<=NN-1; II++) { |
222 |
|
|
tenEigensolve_f(eval, evec, tin); |
223 |
|
|
for (ri=0; ri<3; ri++) { |
224 |
|
|
eval[ri] = AIR_CAST(float, exp(eval[ri])); |
225 |
|
|
eval[ri] = AIR_EXISTS(eval[ri]) ? eval[ri] : 0; |
226 |
|
|
} |
227 |
|
|
tenMakeSingle_f(tout, tin[0], eval, evec); |
228 |
|
|
tin += 7; |
229 |
|
|
tout += 7; |
230 |
|
|
} |
231 |
|
|
break; |
232 |
|
|
} |
233 |
|
|
|
234 |
|
|
/* basic and per-axis info handled by nrrdCopy above */ |
235 |
|
|
return 0; |
236 |
|
|
} |
237 |
|
|
|
238 |
|
|
int |
239 |
|
|
tenSizeNormalize(Nrrd *nout, const Nrrd *nin, double _weight[3], |
240 |
|
|
double amount, double target) { |
241 |
|
|
static const char me[]="tenSizeNormalize"; |
242 |
|
|
funcParm parm; |
243 |
|
|
|
244 |
|
|
ELL_3V_COPY(parm.weight, _weight); |
245 |
|
|
parm.amount = amount; |
246 |
|
|
parm.target = target; |
247 |
|
|
if (theFunc(nout, nin, funcSizeNormalize, &parm)) { |
248 |
|
|
biffAddf(TEN, "%s: trouble", me); |
249 |
|
|
return 1; |
250 |
|
|
} |
251 |
|
|
return 0; |
252 |
|
|
} |
253 |
|
|
|
254 |
|
|
int |
255 |
|
|
tenSizeScale(Nrrd *nout, const Nrrd *nin, double amount) { |
256 |
|
|
static const char me[]="tenSizeScale"; |
257 |
|
|
funcParm parm; |
258 |
|
|
|
259 |
|
|
parm.amount = amount; |
260 |
|
|
if (theFunc(nout, nin, funcSizeScale, &parm)) { |
261 |
|
|
biffAddf(TEN, "%s: trouble", me); |
262 |
|
|
return 1; |
263 |
|
|
} |
264 |
|
|
return 0; |
265 |
|
|
} |
266 |
|
|
|
267 |
|
|
|
268 |
|
|
/* |
269 |
|
|
******** tenAnisoScale |
270 |
|
|
** |
271 |
|
|
** scales the "deviatoric" part of a tensor up or down |
272 |
|
|
*/ |
273 |
|
|
int |
274 |
|
|
tenAnisoScale(Nrrd *nout, const Nrrd *nin, double scale, |
275 |
|
|
int fixDet, int makePositive) { |
276 |
|
|
static const char me[]="tenAnisoScale"; |
277 |
|
|
funcParm parm; |
278 |
|
|
|
279 |
|
|
parm.scale = scale; |
280 |
|
|
parm.fixDet = fixDet; |
281 |
|
|
parm.makePositive = makePositive; |
282 |
|
|
if (theFunc(nout, nin, funcAnisoScale, &parm)) { |
283 |
|
|
biffAddf(TEN, "%s: trouble", me); |
284 |
|
|
return 1; |
285 |
|
|
} |
286 |
|
|
return 0; |
287 |
|
|
} |
288 |
|
|
|
289 |
|
|
/* |
290 |
|
|
******** tenEigenvalueClamp |
291 |
|
|
** |
292 |
|
|
** enstates the given value as the lowest eigenvalue |
293 |
|
|
*/ |
294 |
|
|
int |
295 |
|
|
tenEigenvalueClamp(Nrrd *nout, const Nrrd *nin, double min, double max) { |
296 |
|
|
static const char me[]="tenEigenvalueClamp"; |
297 |
|
|
funcParm parm; |
298 |
|
|
|
299 |
|
|
parm.min = min; |
300 |
|
|
parm.max = max; |
301 |
|
|
if (theFunc(nout, nin, funcEigenvalueClamp, &parm)) { |
302 |
|
|
biffAddf(TEN, "%s: trouble", me); |
303 |
|
|
return 1; |
304 |
|
|
} |
305 |
|
|
return 0; |
306 |
|
|
} |
307 |
|
|
|
308 |
|
|
/* |
309 |
|
|
******** tenEigenvaluePower |
310 |
|
|
** |
311 |
|
|
** raises the eigenvalues to some power |
312 |
|
|
*/ |
313 |
|
|
int |
314 |
|
|
tenEigenvaluePower(Nrrd *nout, const Nrrd *nin, double expo) { |
315 |
|
|
static const char me[]="tenEigenvaluePower"; |
316 |
|
|
funcParm parm; |
317 |
|
|
|
318 |
|
|
parm.expo = expo; |
319 |
|
|
if (theFunc(nout, nin, funcEigenvaluePower, &parm)) { |
320 |
|
|
biffAddf(TEN, "%s: trouble", me); |
321 |
|
|
return 1; |
322 |
|
|
} |
323 |
|
|
return 0; |
324 |
|
|
} |
325 |
|
|
|
326 |
|
|
/* |
327 |
|
|
******** tenEigenvalueAdd |
328 |
|
|
** |
329 |
|
|
** adds something to all eigenvalues |
330 |
|
|
*/ |
331 |
|
|
int |
332 |
|
|
tenEigenvalueAdd(Nrrd *nout, const Nrrd *nin, double val) { |
333 |
|
|
static const char me[]="tenEigenvalueAdd"; |
334 |
|
|
funcParm parm; |
335 |
|
|
|
336 |
|
|
parm.val = val; |
337 |
|
|
if (theFunc(nout, nin, funcEigenvalueAdd, &parm)) { |
338 |
|
|
biffAddf(TEN, "%s: trouble", me); |
339 |
|
|
return 1; |
340 |
|
|
} |
341 |
|
|
return 0; |
342 |
|
|
} |
343 |
|
|
|
344 |
|
|
/* |
345 |
|
|
******** tenEigenvalueMultiply |
346 |
|
|
** |
347 |
|
|
** multiplies eigenvalues by something |
348 |
|
|
*/ |
349 |
|
|
int |
350 |
|
|
tenEigenvalueMultiply(Nrrd *nout, const Nrrd *nin, double val) { |
351 |
|
|
static const char me[]="tenEigenvalueMultiply"; |
352 |
|
|
funcParm parm; |
353 |
|
|
|
354 |
|
|
parm.val = val; |
355 |
|
|
if (theFunc(nout, nin, funcEigenvalueMultiply, &parm)) { |
356 |
|
|
biffAddf(TEN, "%s: trouble", me); |
357 |
|
|
return 1; |
358 |
|
|
} |
359 |
|
|
return 0; |
360 |
|
|
} |
361 |
|
|
|
362 |
|
|
/* |
363 |
|
|
******** tenLog |
364 |
|
|
** |
365 |
|
|
** takes the logarithm (by taking the log of the eigenvalues) |
366 |
|
|
*/ |
367 |
|
|
int |
368 |
|
|
tenLog(Nrrd *nout, const Nrrd *nin) { |
369 |
|
|
static const char me[]="tenLog"; |
370 |
|
|
funcParm parm; |
371 |
|
|
|
372 |
|
|
if (theFunc(nout, nin, funcLog, &parm)) { |
373 |
|
|
biffAddf(TEN, "%s: trouble", me); |
374 |
|
|
return 1; |
375 |
|
|
} |
376 |
|
|
return 0; |
377 |
|
|
} |
378 |
|
|
|
379 |
|
|
/* |
380 |
|
|
******** tenExp |
381 |
|
|
** |
382 |
|
|
** takes the exp() (by taking exp() of the eigenvalues) |
383 |
|
|
*/ |
384 |
|
|
int |
385 |
|
|
tenExp(Nrrd *nout, const Nrrd *nin) { |
386 |
|
|
static const char me[]="tenExp"; |
387 |
|
|
funcParm parm; |
388 |
|
|
|
389 |
|
|
if (theFunc(nout, nin, funcExp, &parm)) { |
390 |
|
|
biffAddf(TEN, "%s: trouble", me); |
391 |
|
|
return 1; |
392 |
|
|
} |
393 |
|
|
return 0; |
394 |
|
|
} |