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 "pull.h" |
25 |
|
|
#include "privatePull.h" |
26 |
|
|
|
27 |
|
|
const char * |
28 |
|
|
_pullInterTypeStr[PULL_INTER_TYPE_MAX+1] = { |
29 |
|
|
"(unknown_inter)", |
30 |
|
|
"justR", |
31 |
|
|
"univariate", |
32 |
|
|
"separable", |
33 |
|
|
"additive" |
34 |
|
|
}; |
35 |
|
|
|
36 |
|
|
const char * |
37 |
|
|
_pullInterTypeStrEqv[] = { |
38 |
|
|
"r", "justr", |
39 |
|
|
"univariate", "univar", "uni", |
40 |
|
|
"separable", "separ", "sep", |
41 |
|
|
"additive", "add", |
42 |
|
|
"" |
43 |
|
|
}; |
44 |
|
|
|
45 |
|
|
const int |
46 |
|
|
_pullInterTypeValEqv[] = { |
47 |
|
|
pullInterTypeJustR, pullInterTypeJustR, |
48 |
|
|
pullInterTypeUnivariate, pullInterTypeUnivariate, pullInterTypeUnivariate, |
49 |
|
|
pullInterTypeSeparable, pullInterTypeSeparable, pullInterTypeSeparable, |
50 |
|
|
pullInterTypeAdditive, pullInterTypeAdditive |
51 |
|
|
}; |
52 |
|
|
|
53 |
|
|
const airEnum |
54 |
|
|
_pullInterType = { |
55 |
|
|
"interaction type", |
56 |
|
|
PULL_INTER_TYPE_MAX, |
57 |
|
|
_pullInterTypeStr, NULL, |
58 |
|
|
NULL, |
59 |
|
|
_pullInterTypeStrEqv, _pullInterTypeValEqv, |
60 |
|
|
AIR_FALSE |
61 |
|
|
}; |
62 |
|
|
const airEnum *const |
63 |
|
|
pullInterType = &_pullInterType; |
64 |
|
|
|
65 |
|
|
#define SPRING "spring" |
66 |
|
|
#define GAUSS "gauss" |
67 |
|
|
#define BSPLN "bspln" |
68 |
|
|
#define BUTTER "butter" |
69 |
|
|
#define COTAN "cotan" |
70 |
|
|
#define CUBIC "cubic" |
71 |
|
|
#define QUARTIC "quartic" |
72 |
|
|
#define CWELL "cwell" |
73 |
|
|
#define BWELL "bwell" |
74 |
|
|
#define QWELL "qwell" |
75 |
|
|
#define HWELL "hwell" |
76 |
|
|
#define ZERO "zero" |
77 |
|
|
#define BPARAB "bparab" |
78 |
|
|
|
79 |
|
|
const char * |
80 |
|
|
_pullEnergyTypeStr[PULL_ENERGY_TYPE_MAX+1] = { |
81 |
|
|
"(unknown_energy)", |
82 |
|
|
SPRING, |
83 |
|
|
GAUSS, |
84 |
|
|
BSPLN, |
85 |
|
|
BUTTER, |
86 |
|
|
COTAN, |
87 |
|
|
CUBIC, |
88 |
|
|
QUARTIC, |
89 |
|
|
CWELL, |
90 |
|
|
BWELL, |
91 |
|
|
QWELL, |
92 |
|
|
HWELL, |
93 |
|
|
ZERO, |
94 |
|
|
BPARAB |
95 |
|
|
}; |
96 |
|
|
|
97 |
|
|
const char * |
98 |
|
|
_pullEnergyTypeDesc[PULL_ENERGY_TYPE_MAX+1] = { |
99 |
|
|
"unknown_energy", |
100 |
|
|
"Hooke's law-based potential, with a tunable region of attraction", |
101 |
|
|
"Gaussian potential", |
102 |
|
|
"uniform cubic B-spline basis function", |
103 |
|
|
"like a Gaussian, but a lot wider", |
104 |
|
|
"Cotangent-based potential (from Meyer et al. SMI '05)", |
105 |
|
|
"Cubic thing", |
106 |
|
|
"Quartic thing", |
107 |
|
|
"Piecewice cubic with tunable well location and depth", |
108 |
|
|
"Better piecewice cubic with tunable well location and depth", |
109 |
|
|
"Single quartic with tunable well location", |
110 |
|
|
"Single heptic with tunable well location", |
111 |
|
|
"no energy", |
112 |
|
|
"butterworth-windowed spatial repel and scale attract" |
113 |
|
|
}; |
114 |
|
|
|
115 |
|
|
const airEnum |
116 |
|
|
_pullEnergyType = { |
117 |
|
|
"energy", |
118 |
|
|
PULL_ENERGY_TYPE_MAX, |
119 |
|
|
_pullEnergyTypeStr, NULL, |
120 |
|
|
_pullEnergyTypeDesc, |
121 |
|
|
NULL, NULL, |
122 |
|
|
AIR_FALSE |
123 |
|
|
}; |
124 |
|
|
const airEnum *const |
125 |
|
|
pullEnergyType = &_pullEnergyType; |
126 |
|
|
|
127 |
|
|
double |
128 |
|
|
_pullEnergyNoWell(double *wx, const double *parm) { |
129 |
|
|
|
130 |
|
|
AIR_UNUSED(wx); |
131 |
|
|
AIR_UNUSED(parm); |
132 |
|
|
return 0.0; |
133 |
|
|
} |
134 |
|
|
|
135 |
|
|
/* ---------------------------------------------------------------- |
136 |
|
|
** ------------------------------ UNKNOWN ------------------------- |
137 |
|
|
** ---------------------------------------------------------------- |
138 |
|
|
*/ |
139 |
|
|
double |
140 |
|
|
_pullEnergyUnknownEval(double *denr, double dist, const double *parm) { |
141 |
|
|
static const char me[]="_pullEnergyUnknownEval"; |
142 |
|
|
|
143 |
|
|
AIR_UNUSED(dist); |
144 |
|
|
AIR_UNUSED(parm); |
145 |
|
|
*denr = AIR_NAN; |
146 |
|
|
fprintf(stderr, "%s: ERROR- using unknown energy.\n", me); |
147 |
|
|
return AIR_NAN; |
148 |
|
|
} |
149 |
|
|
|
150 |
|
|
pullEnergy |
151 |
|
|
_pullEnergyUnknown = { |
152 |
|
|
"unknown", |
153 |
|
|
0, |
154 |
|
|
_pullEnergyNoWell, |
155 |
|
|
_pullEnergyUnknownEval |
156 |
|
|
}; |
157 |
|
|
const pullEnergy *const |
158 |
|
|
pullEnergyUnknown = &_pullEnergyUnknown; |
159 |
|
|
|
160 |
|
|
/* ---------------------------------------------------------------- |
161 |
|
|
** ------------------------------ SPRING -------------------------- |
162 |
|
|
** ---------------------------------------------------------------- |
163 |
|
|
** 1 parms: |
164 |
|
|
** parm[0]: width of pull region. Used to be width beyond 1.0, but |
165 |
|
|
** now things are scrunched to fit both repelling and attractive |
166 |
|
|
** region inside [0,1] |
167 |
|
|
** |
168 |
|
|
** learned: "1/2" is not 0.5 !!!!! |
169 |
|
|
*/ |
170 |
|
|
double |
171 |
|
|
_pullEnergySpringEval(double *denr, double dist, const double *parm) { |
172 |
|
|
/* static const char me[]="_pullEnergySpringEval"; */ |
173 |
|
|
double enr, xx, pull; |
174 |
|
|
|
175 |
|
|
pull = parm[0]; |
176 |
|
|
/* support used to be [0,1 + pull], but now is scrunched to [0,1], |
177 |
|
|
so hack "dist" to match old parameterization */ |
178 |
|
|
dist = AIR_AFFINE(0, dist, 1, 0, 1+pull); |
179 |
|
|
xx = dist - 1.0; |
180 |
|
|
if (xx > pull) { |
181 |
|
|
enr = 0; |
182 |
|
|
*denr = 0; |
183 |
|
|
} else if (xx > 0) { |
184 |
|
|
enr = xx*xx*(xx*xx/(4*pull*pull) - 2*xx/(3*pull) + 1.0/2.0); |
185 |
|
|
*denr = xx*(xx*xx/(pull*pull) - 2*xx/pull + 1); |
186 |
|
|
} else { |
187 |
|
|
enr = xx*xx/2; |
188 |
|
|
*denr = xx; |
189 |
|
|
} |
190 |
|
|
/* |
191 |
|
|
if (!AIR_EXISTS(ret)) { |
192 |
|
|
printf("!%s: dist=%g, pull=%g, blah=%d --> ret=%g\n", |
193 |
|
|
me, dist, pull, blah, ret); |
194 |
|
|
} |
195 |
|
|
*/ |
196 |
|
|
return enr; |
197 |
|
|
} |
198 |
|
|
|
199 |
|
|
int |
200 |
|
|
_pullEnergySpringWell(const double *parm) { |
201 |
|
|
|
202 |
|
|
return (parm[0] > 0.0); |
203 |
|
|
} |
204 |
|
|
|
205 |
|
|
const pullEnergy |
206 |
|
|
_pullEnergySpring = { |
207 |
|
|
SPRING, |
208 |
|
|
1, |
209 |
|
|
_pullEnergyNoWell, /* HEY: is this true? */ |
210 |
|
|
_pullEnergySpringEval |
211 |
|
|
}; |
212 |
|
|
const pullEnergy *const |
213 |
|
|
pullEnergySpring = &_pullEnergySpring; |
214 |
|
|
|
215 |
|
|
/* ---------------------------------------------------------------- |
216 |
|
|
** ------------------------------ GAUSS -------------------------- |
217 |
|
|
** ---------------------------------------------------------------- |
218 |
|
|
** 0 parms: for simplicity we're now always cutting off at 4 sigmas |
219 |
|
|
*/ |
220 |
|
|
/* HEY: copied from teem/src/nrrd/kernel.c */ |
221 |
|
|
#define _GAUSS(x, sig, cut) ( \ |
222 |
|
|
x >= sig*cut ? 0 \ |
223 |
|
|
: exp(-x*x/(2.0*sig*sig))) |
224 |
|
|
|
225 |
|
|
#define _DGAUSS(x, sig, cut) ( \ |
226 |
|
|
x >= sig*cut ? 0 \ |
227 |
|
|
: -exp(-x*x/(2.0*sig*sig))*(x/(sig*sig))) |
228 |
|
|
|
229 |
|
|
double |
230 |
|
|
_pullEnergyGaussEval(double *denr, double dist, const double *parm) { |
231 |
|
|
|
232 |
|
|
AIR_UNUSED(parm); |
233 |
|
|
*denr = _DGAUSS(dist, 0.25, 4); |
234 |
|
|
return _GAUSS(dist, 0.25, 4); |
235 |
|
|
} |
236 |
|
|
|
237 |
|
|
const pullEnergy |
238 |
|
|
_pullEnergyGauss = { |
239 |
|
|
GAUSS, |
240 |
|
|
0, |
241 |
|
|
_pullEnergyNoWell, |
242 |
|
|
_pullEnergyGaussEval |
243 |
|
|
}; |
244 |
|
|
const pullEnergy *const |
245 |
|
|
pullEnergyGauss = &_pullEnergyGauss; |
246 |
|
|
|
247 |
|
|
/* ---------------------------------------------------------------- |
248 |
|
|
** ------------------------------ BSPLN -------------------------- |
249 |
|
|
** ---------------------------------------------------------------- |
250 |
|
|
** 0 parms |
251 |
|
|
*/ |
252 |
|
|
/* HEY: copied from teem/src/nrrd/bsplKernel.c */ |
253 |
|
|
#define BSPL(ret, t, x) \ |
254 |
|
|
if (x < 1) { \ |
255 |
|
|
ret = (4 + 3*(-2 + x)*x*x)/6; \ |
256 |
|
|
} else if (x < 2) { \ |
257 |
|
|
t = (-2 + x); \ |
258 |
|
|
ret = -t*t*t/6; \ |
259 |
|
|
} else { \ |
260 |
|
|
ret = 0; \ |
261 |
|
|
} |
262 |
|
|
|
263 |
|
|
#define DBSPL(ret, t, x) \ |
264 |
|
|
if (x < 1) { \ |
265 |
|
|
ret = (-4 + 3*x)*x/2; \ |
266 |
|
|
} else if (x < 2) { \ |
267 |
|
|
t = (-2 + x); \ |
268 |
|
|
ret = -t*t/2; \ |
269 |
|
|
} else { \ |
270 |
|
|
ret = 0; \ |
271 |
|
|
} |
272 |
|
|
|
273 |
|
|
double |
274 |
|
|
_pullEnergyBsplnEval(double *denr, double dist, const double *parm) { |
275 |
|
|
double tmp, ret; |
276 |
|
|
|
277 |
|
|
AIR_UNUSED(parm); |
278 |
|
|
dist *= 2; |
279 |
|
|
DBSPL(*denr, tmp, dist); |
280 |
|
|
*denr *= 2; |
281 |
|
|
BSPL(ret, tmp, dist); |
282 |
|
|
return ret; |
283 |
|
|
} |
284 |
|
|
|
285 |
|
|
const pullEnergy |
286 |
|
|
_pullEnergyBspln = { |
287 |
|
|
BSPLN, |
288 |
|
|
0, |
289 |
|
|
_pullEnergyNoWell, |
290 |
|
|
_pullEnergyBsplnEval |
291 |
|
|
}; |
292 |
|
|
const pullEnergy *const |
293 |
|
|
pullEnergyBspln = &_pullEnergyBspln; |
294 |
|
|
|
295 |
|
|
|
296 |
|
|
/* ---------------------------------------------------------------- |
297 |
|
|
** ------------------------------ BUTTER -------------------------- |
298 |
|
|
** ---------------------------------------------------------------- |
299 |
|
|
** 2 parms: order (an integer) and "cut-ff" (where height==0.5) |
300 |
|
|
*/ |
301 |
|
|
|
302 |
|
|
double |
303 |
|
|
_pullEnergyButterworthEval(double *denr, double x, const double *parm) { |
304 |
|
|
int n; |
305 |
|
|
double cut, denom, enr; |
306 |
|
|
|
307 |
|
|
n = AIR_CAST(int, parm[0]); |
308 |
|
|
cut = parm[1]; |
309 |
|
|
denom = 1 + airIntPow(x/cut, 2*n); |
310 |
|
|
enr = 1/denom; |
311 |
|
|
*denr = -2*n*airIntPow(x/cut, 2*n - 1)*enr*enr/cut; |
312 |
|
|
return enr; |
313 |
|
|
} |
314 |
|
|
|
315 |
|
|
const pullEnergy |
316 |
|
|
_pullEnergyButterworth= { |
317 |
|
|
BUTTER, |
318 |
|
|
2, |
319 |
|
|
_pullEnergyNoWell, |
320 |
|
|
_pullEnergyButterworthEval |
321 |
|
|
}; |
322 |
|
|
const pullEnergy *const |
323 |
|
|
pullEnergyButterworth = &_pullEnergyButterworth; |
324 |
|
|
|
325 |
|
|
/* ---------------------------------------------------------------- |
326 |
|
|
** ------------------------------ COTAN --------------------------- |
327 |
|
|
** ---------------------------------------------------------------- |
328 |
|
|
** 0 parms! |
329 |
|
|
*/ |
330 |
|
|
double |
331 |
|
|
_pullEnergyCotanEval(double *denr, double dist, const double *parm) { |
332 |
|
|
double pot, cc, enr; |
333 |
|
|
|
334 |
|
|
AIR_UNUSED(parm); |
335 |
|
|
pot = AIR_PI/2.0; |
336 |
|
|
cc = 1.0/(FLT_MIN + tan(dist*pot)); |
337 |
|
|
enr = dist > 1 ? 0 : cc + dist*pot - pot; |
338 |
|
|
*denr = dist > 1 ? 0 : -cc*cc*pot; |
339 |
|
|
return enr; |
340 |
|
|
} |
341 |
|
|
|
342 |
|
|
const pullEnergy |
343 |
|
|
_pullEnergyCotan = { |
344 |
|
|
COTAN, |
345 |
|
|
0, |
346 |
|
|
_pullEnergyNoWell, |
347 |
|
|
_pullEnergyCotanEval |
348 |
|
|
}; |
349 |
|
|
const pullEnergy *const |
350 |
|
|
pullEnergyCotan = &_pullEnergyCotan; |
351 |
|
|
|
352 |
|
|
/* ---------------------------------------------------------------- |
353 |
|
|
** ------------------------------ CUBIC --------------------------- |
354 |
|
|
** ---------------------------------------------------------------- |
355 |
|
|
** 0 parms! |
356 |
|
|
*/ |
357 |
|
|
double |
358 |
|
|
_pullEnergyCubicEval(double *denr, double dist, const double *parm) { |
359 |
|
|
double omr, enr; |
360 |
|
|
|
361 |
|
|
AIR_UNUSED(parm); |
362 |
|
|
if (dist <= 1) { |
363 |
|
|
omr = 1 - dist; |
364 |
|
|
enr = omr*omr*omr; |
365 |
|
|
*denr = -3*omr*omr; |
366 |
|
|
} else { |
367 |
|
|
enr = *denr = 0; |
368 |
|
|
} |
369 |
|
|
return enr; |
370 |
|
|
} |
371 |
|
|
|
372 |
|
|
const pullEnergy |
373 |
|
|
_pullEnergyCubic = { |
374 |
|
|
CUBIC, |
375 |
|
|
0, |
376 |
|
|
_pullEnergyNoWell, |
377 |
|
|
_pullEnergyCubicEval |
378 |
|
|
}; |
379 |
|
|
const pullEnergy *const |
380 |
|
|
pullEnergyCubic = &_pullEnergyCubic; |
381 |
|
|
|
382 |
|
|
/* ---------------------------------------------------------------- |
383 |
|
|
** ----------------------------- QUARTIC -------------------------- |
384 |
|
|
** ---------------------------------------------------------------- |
385 |
|
|
** 0 parms! |
386 |
|
|
*/ |
387 |
|
|
double |
388 |
|
|
_pullEnergyQuarticEval(double *denr, double dist, const double *parm) { |
389 |
|
|
double omr, enr; |
390 |
|
|
|
391 |
|
|
AIR_UNUSED(parm); |
392 |
|
|
if (dist <= 1) { |
393 |
|
|
omr = 1 - dist; |
394 |
|
|
enr = 2.132*omr*omr*omr*omr; |
395 |
|
|
*denr = -4*2.132*omr*omr*omr; |
396 |
|
|
} else { |
397 |
|
|
enr = *denr = 0; |
398 |
|
|
} |
399 |
|
|
return enr; |
400 |
|
|
} |
401 |
|
|
|
402 |
|
|
const pullEnergy |
403 |
|
|
_pullEnergyQuartic = { |
404 |
|
|
QUARTIC, |
405 |
|
|
0, |
406 |
|
|
_pullEnergyNoWell, |
407 |
|
|
_pullEnergyQuarticEval |
408 |
|
|
}; |
409 |
|
|
const pullEnergy *const |
410 |
|
|
pullEnergyQuartic = &_pullEnergyQuartic; |
411 |
|
|
|
412 |
|
|
/* ---------------------------------------------------------------- |
413 |
|
|
** ------------------ tunable piece-wise cubic -------------------- |
414 |
|
|
** ---------------------------------------------------------------- |
415 |
|
|
** 2 parm: wellX, wellY |
416 |
|
|
*/ |
417 |
|
|
double |
418 |
|
|
_pullEnergyCubicWellEval(double *denr, double x, const double *parm) { |
419 |
|
|
double a, b, c, d, e, wx, wy, enr; |
420 |
|
|
|
421 |
|
|
wx = parm[0]; |
422 |
|
|
wy = parm[1]; |
423 |
|
|
a = (3*(-1 + wy))/wx; |
424 |
|
|
b = (-3*(-1 + wy))/(wx*wx); |
425 |
|
|
c = -(1 - wy)/(wx*wx*wx); |
426 |
|
|
d = (-3*wy)/((wx-1)*(wx-1)); |
427 |
|
|
e = (-2*wy)/((wx-1)*(wx-1)*(wx-1)); |
428 |
|
|
if (x < wx) { |
429 |
|
|
enr = 1 + x*(a + x*(b + c*x)); |
430 |
|
|
*denr = a + x*(2*b + 3*c*x); |
431 |
|
|
} else if (x < 1) { |
432 |
|
|
double _x; |
433 |
|
|
_x = x - wx; |
434 |
|
|
enr = wy + _x*_x*(d + e*_x); |
435 |
|
|
*denr = _x*(2*d + 3*e*_x); |
436 |
|
|
} else { |
437 |
|
|
enr = 0; |
438 |
|
|
*denr = 0; |
439 |
|
|
} |
440 |
|
|
return enr; |
441 |
|
|
} |
442 |
|
|
|
443 |
|
|
double |
444 |
|
|
_pullEnergyCubicWellWell(double *wx, const double *parm) { |
445 |
|
|
|
446 |
|
|
*wx = parm[0]; |
447 |
|
|
return AIR_MIN(0.0, parm[1]); |
448 |
|
|
} |
449 |
|
|
|
450 |
|
|
const pullEnergy |
451 |
|
|
_pullEnergyCubicWell = { |
452 |
|
|
CWELL, |
453 |
|
|
2, |
454 |
|
|
_pullEnergyCubicWellWell, |
455 |
|
|
_pullEnergyCubicWellEval |
456 |
|
|
}; |
457 |
|
|
const pullEnergy *const |
458 |
|
|
pullEnergyCubicWell = &_pullEnergyCubicWell; |
459 |
|
|
|
460 |
|
|
/* ---------------------------------------------------------------- |
461 |
|
|
** --------------- better tunable piece-wise cubic ---------------- |
462 |
|
|
** ---------------------------------------------------------------- |
463 |
|
|
** 2 parm: wellX, wellY |
464 |
|
|
*/ |
465 |
|
|
double |
466 |
|
|
_pullEnergyBetterCubicWellEval(double *denr, double x, const double *parm) { |
467 |
|
|
double a, b, c, d, e, f, g, wx, wy, xmo, xmoo, xmooo, enr; |
468 |
|
|
|
469 |
|
|
/* HEY: it would be good if there was a place to store these |
470 |
|
|
intermediate values, so that they don't need to be computed |
471 |
|
|
for *every*single* inter-particle interaction */ |
472 |
|
|
wx = parm[0]; |
473 |
|
|
wy = parm[1]; |
474 |
|
|
xmo = wx-1; |
475 |
|
|
xmoo = xmo*xmo; |
476 |
|
|
xmooo = xmoo*xmo; |
477 |
|
|
a = -3*(xmoo + (-1 + 2*wx)*wy)/(xmoo*wx); |
478 |
|
|
b = 3*(xmoo + (-1 + wx*(2 + wx))*wy)/(xmoo*wx*wx); |
479 |
|
|
c = (-1 + wy - wx*(-2 + wx + 2*(1 + wx)*wy))/(xmoo*wx*wx*wx); |
480 |
|
|
d = ((-1 + 3*wx)*wy)/xmooo; |
481 |
|
|
e = -(6*wx*wy)/xmooo; |
482 |
|
|
f = (3*(1 + wx)*wy)/xmooo; |
483 |
|
|
g = -(2*wy)/xmooo; |
484 |
|
|
if (x < wx) { |
485 |
|
|
enr = 1 + x*(a + x*(b + x*c)); |
486 |
|
|
*denr = a + x*(2*b + x*3*c); |
487 |
|
|
} else if (x < 1) { |
488 |
|
|
enr = d + x*(e + x*(f + x*g)); |
489 |
|
|
*denr = e + x*(2*f + x*3*g); |
490 |
|
|
} else { |
491 |
|
|
enr = 0; |
492 |
|
|
*denr = 0; |
493 |
|
|
} |
494 |
|
|
return enr; |
495 |
|
|
} |
496 |
|
|
|
497 |
|
|
double |
498 |
|
|
_pullEnergyBetterCubicWellWell(double *wx, const double *parm) { |
499 |
|
|
|
500 |
|
|
*wx = parm[0]; |
501 |
|
|
return AIR_MIN(0.0, parm[1]); |
502 |
|
|
} |
503 |
|
|
|
504 |
|
|
const pullEnergy |
505 |
|
|
_pullEnergyBetterCubicWell = { |
506 |
|
|
BWELL, |
507 |
|
|
2, |
508 |
|
|
_pullEnergyBetterCubicWellWell, |
509 |
|
|
_pullEnergyBetterCubicWellEval |
510 |
|
|
}; |
511 |
|
|
const pullEnergy *const |
512 |
|
|
pullEnergyBetterCubicWell = &_pullEnergyBetterCubicWell; |
513 |
|
|
|
514 |
|
|
/* ---------------------------------------------------------------- |
515 |
|
|
** ----------------- tunable single quartic well ------------------ |
516 |
|
|
** ---------------------------------------------------------------- |
517 |
|
|
** 1 parm: well radius |
518 |
|
|
*/ |
519 |
|
|
double |
520 |
|
|
_pullEnergyQuarticWellEval(double *denr, double x, const double *parm) { |
521 |
|
|
double a, b, c, d, w, enr; |
522 |
|
|
|
523 |
|
|
w = parm[0]; |
524 |
|
|
a = (12*w)/(1 - 4*w); |
525 |
|
|
b = 3 + 9/(-1 + 4*w); |
526 |
|
|
c = (8 + 4*w)/(1 - 4*w); |
527 |
|
|
d = 3/(-1 + 4*w); |
528 |
|
|
enr = 1 + x*(a + x*(b + x*(c + x*d))); |
529 |
|
|
*denr = a + x*(2*b + x*(3*c + x*4*d)); |
530 |
|
|
return enr; |
531 |
|
|
} |
532 |
|
|
|
533 |
|
|
double |
534 |
|
|
_pullEnergyQuarticWellWell(double *wx, const double *parm) { |
535 |
|
|
double t; |
536 |
|
|
|
537 |
|
|
*wx = parm[0]; |
538 |
|
|
t = *wx - 1; |
539 |
|
|
return t*t*t*t/(4*(*wx) - 1); |
540 |
|
|
} |
541 |
|
|
|
542 |
|
|
const pullEnergy |
543 |
|
|
_pullEnergyQuarticWell = { |
544 |
|
|
QWELL, |
545 |
|
|
1, |
546 |
|
|
_pullEnergyQuarticWellWell, |
547 |
|
|
_pullEnergyQuarticWellEval |
548 |
|
|
}; |
549 |
|
|
const pullEnergy *const |
550 |
|
|
pullEnergyQuarticWell = &_pullEnergyQuarticWell; |
551 |
|
|
|
552 |
|
|
/* ---------------------------------------------------------------- |
553 |
|
|
** ------------------ tunable single heptic well ------------------ |
554 |
|
|
** ---------------------------------------------------------------- |
555 |
|
|
** 1 parm: well radius |
556 |
|
|
*/ |
557 |
|
|
double |
558 |
|
|
_pullEnergyHepticWellEval(double *denr, double x, const double *parm) { |
559 |
|
|
double a, b, c, d, e, f, g, w, enr; |
560 |
|
|
|
561 |
|
|
w = parm[0]; |
562 |
|
|
a = (42*w)/(1 - 7*w); |
563 |
|
|
b = 15 + 36/(-1 + 7*w); |
564 |
|
|
c = -20 + 90/(1 - 7*w); |
565 |
|
|
d = (105*(1 + w))/(-1 + 7*w); |
566 |
|
|
e = -((42*(2 + w))/(-1 + 7*w)); |
567 |
|
|
f = (7*(5 + w))/(-1 + 7*w); |
568 |
|
|
g = 6/(1 - 7*w); |
569 |
|
|
enr = 1 + x*(a + x*(b + x*(c + x*(d + x*(e + x*(f + g*x)))))); |
570 |
|
|
*denr = a + x*(2*b + x*(3*c + x*(4*d + x*(5*e + x*(6*f + x*7*g))))); |
571 |
|
|
return enr; |
572 |
|
|
} |
573 |
|
|
|
574 |
|
|
double |
575 |
|
|
_pullEnergyHepticWellWell(double *wx, const double *parm) { |
576 |
|
|
double t; |
577 |
|
|
|
578 |
|
|
*wx = parm[0]; |
579 |
|
|
t = *wx - 1; |
580 |
|
|
return t*t*t*t*t*t*t/(7*(*wx) - 1); |
581 |
|
|
} |
582 |
|
|
|
583 |
|
|
const pullEnergy |
584 |
|
|
_pullEnergyHepticWell = { |
585 |
|
|
HWELL, |
586 |
|
|
1, |
587 |
|
|
_pullEnergyHepticWellWell, |
588 |
|
|
_pullEnergyHepticWellEval |
589 |
|
|
}; |
590 |
|
|
const pullEnergy *const |
591 |
|
|
pullEnergyHepticWell = &_pullEnergyHepticWell; |
592 |
|
|
|
593 |
|
|
/* ---------------------------------------------------------------- |
594 |
|
|
** ------------------------------- ZERO --------------------------- |
595 |
|
|
** ---------------------------------------------------------------- |
596 |
|
|
** 0 parms: |
597 |
|
|
*/ |
598 |
|
|
double |
599 |
|
|
_pullEnergyZeroEval(double *denr, double dist, const double *parm) { |
600 |
|
|
|
601 |
|
|
AIR_UNUSED(dist); |
602 |
|
|
AIR_UNUSED(parm); |
603 |
|
|
*denr = 0; |
604 |
|
|
return 0; |
605 |
|
|
} |
606 |
|
|
|
607 |
|
|
const pullEnergy |
608 |
|
|
_pullEnergyZero = { |
609 |
|
|
ZERO, |
610 |
|
|
0, |
611 |
|
|
_pullEnergyNoWell, |
612 |
|
|
_pullEnergyZeroEval |
613 |
|
|
}; |
614 |
|
|
const pullEnergy *const |
615 |
|
|
pullEnergyZero = &_pullEnergyZero; |
616 |
|
|
|
617 |
|
|
/* ---------------------------------------------------------------- |
618 |
|
|
** ------------------------------- BPARAB ------------------------- |
619 |
|
|
** ---------------------------------------------------------------- |
620 |
|
|
** 3 parms, the first two are for butterworth, |
621 |
|
|
** parm[2] is a shift (probably negative) on the parabola |
622 |
|
|
*/ |
623 |
|
|
double |
624 |
|
|
_pullEnergyBParabEval(double *denr, double x, const double *parm) { |
625 |
|
|
double ben, dben; |
626 |
|
|
|
627 |
|
|
ben = _pullEnergyButterworthEval(&dben, x, parm); |
628 |
|
|
*denr = 2*x*ben + x*x*dben; |
629 |
|
|
return (x*x + parm[2])*ben; |
630 |
|
|
} |
631 |
|
|
|
632 |
|
|
const pullEnergy |
633 |
|
|
_pullEnergyButterworthParabola = { |
634 |
|
|
BPARAB, |
635 |
|
|
3, |
636 |
|
|
_pullEnergyNoWell, |
637 |
|
|
_pullEnergyBParabEval |
638 |
|
|
}; |
639 |
|
|
const pullEnergy *const |
640 |
|
|
pullEnergyButterworthParabola = &_pullEnergyButterworthParabola; |
641 |
|
|
|
642 |
|
|
/* ---------------------------------------------------------------- |
643 |
|
|
** ---------------------------------------------------------------- |
644 |
|
|
** ---------------------------------------------------------------- |
645 |
|
|
*/ |
646 |
|
|
|
647 |
|
|
const pullEnergy *const pullEnergyAll[PULL_ENERGY_TYPE_MAX+1] = { |
648 |
|
|
&_pullEnergyUnknown, /* 0 */ |
649 |
|
|
&_pullEnergySpring, /* 1 */ |
650 |
|
|
&_pullEnergyGauss, /* 2 */ |
651 |
|
|
&_pullEnergyBspln, /* 3 */ |
652 |
|
|
&_pullEnergyButterworth, /* 4 */ |
653 |
|
|
&_pullEnergyCotan, /* 5 */ |
654 |
|
|
&_pullEnergyCubic, /* 6 */ |
655 |
|
|
&_pullEnergyQuartic, /* 7 */ |
656 |
|
|
&_pullEnergyCubicWell, /* 8 */ |
657 |
|
|
&_pullEnergyBetterCubicWell, /* 9 */ |
658 |
|
|
&_pullEnergyQuarticWell, /* 10 */ |
659 |
|
|
&_pullEnergyHepticWell, /* 11 */ |
660 |
|
|
&_pullEnergyZero, /* 12 */ |
661 |
|
|
&_pullEnergyButterworthParabola /* 13 */ |
662 |
|
|
}; |
663 |
|
|
|
664 |
|
|
pullEnergySpec * |
665 |
|
|
pullEnergySpecNew() { |
666 |
|
|
pullEnergySpec *ensp; |
667 |
|
|
int pi; |
668 |
|
|
|
669 |
|
|
ensp = (pullEnergySpec *)calloc(1, sizeof(pullEnergySpec)); |
670 |
|
|
if (ensp) { |
671 |
|
|
ensp->energy = pullEnergyUnknown; |
672 |
|
|
for (pi=0; pi<PULL_ENERGY_PARM_NUM; pi++) { |
673 |
|
|
ensp->parm[pi] = AIR_NAN; |
674 |
|
|
} |
675 |
|
|
} |
676 |
|
|
return ensp; |
677 |
|
|
} |
678 |
|
|
|
679 |
|
|
void |
680 |
|
|
pullEnergySpecSet(pullEnergySpec *ensp, const pullEnergy *energy, |
681 |
|
|
const double parm[PULL_ENERGY_PARM_NUM]) { |
682 |
|
|
unsigned int pi; |
683 |
|
|
|
684 |
|
|
if (ensp && energy && parm) { |
685 |
|
|
ensp->energy = energy; |
686 |
|
|
for (pi=0; pi<PULL_ENERGY_PARM_NUM; pi++) { |
687 |
|
|
ensp->parm[pi] = parm[pi]; |
688 |
|
|
} |
689 |
|
|
} |
690 |
|
|
return; |
691 |
|
|
} |
692 |
|
|
|
693 |
|
|
void |
694 |
|
|
pullEnergySpecCopy(pullEnergySpec *esDst, const pullEnergySpec *esSrc) { |
695 |
|
|
|
696 |
|
|
if (esDst && esSrc) { |
697 |
|
|
pullEnergySpecSet(esDst, esSrc->energy, esSrc->parm); |
698 |
|
|
} |
699 |
|
|
return; |
700 |
|
|
} |
701 |
|
|
|
702 |
|
|
pullEnergySpec * |
703 |
|
|
pullEnergySpecNix(pullEnergySpec *ensp) { |
704 |
|
|
|
705 |
|
|
airFree(ensp); |
706 |
|
|
return NULL; |
707 |
|
|
} |
708 |
|
|
|
709 |
|
|
int |
710 |
|
|
pullEnergySpecParse(pullEnergySpec *ensp, const char *_str) { |
711 |
|
|
static const char me[]="pullEnergySpecParse"; |
712 |
|
|
char *str, *col, *_pstr, *pstr; |
713 |
|
|
int etype; |
714 |
|
|
unsigned int pi, haveParm; |
715 |
|
|
airArray *mop; |
716 |
|
|
double pval; |
717 |
|
|
|
718 |
|
|
if (!( ensp && _str )) { |
719 |
|
|
biffAddf(PULL, "%s: got NULL pointer", me); |
720 |
|
|
return 1; |
721 |
|
|
} |
722 |
|
|
|
723 |
|
|
/* see if its the name of something that needs no parameters */ |
724 |
|
|
etype = airEnumVal(pullEnergyType, _str); |
725 |
|
|
if (pullEnergyTypeUnknown != etype) { |
726 |
|
|
/* the string is the name of some energy */ |
727 |
|
|
ensp->energy = pullEnergyAll[etype]; |
728 |
|
|
if (0 != ensp->energy->parmNum) { |
729 |
|
|
biffAddf(PULL, "%s: need %u parm%s for %s energy, but got none", me, |
730 |
|
|
ensp->energy->parmNum, |
731 |
|
|
(1 == ensp->energy->parmNum ? "" : "s"), |
732 |
|
|
ensp->energy->name); |
733 |
|
|
return 1; |
734 |
|
|
} |
735 |
|
|
/* the energy needs 0 parameters */ |
736 |
|
|
for (pi=0; pi<PULL_ENERGY_PARM_NUM; pi++) { |
737 |
|
|
ensp->parm[pi] = AIR_NAN; |
738 |
|
|
} |
739 |
|
|
return 0; |
740 |
|
|
} |
741 |
|
|
|
742 |
|
|
/* start parsing parms after ':' */ |
743 |
|
|
mop = airMopNew(); |
744 |
|
|
str = airStrdup(_str); |
745 |
|
|
airMopAdd(mop, str, (airMopper)airFree, airMopAlways); |
746 |
|
|
col = strchr(str, ':'); |
747 |
|
|
if (!col) { |
748 |
|
|
biffAddf(PULL, "%s: either \"%s\" is not a recognized energy, " |
749 |
|
|
"or it is an energy with parameters, and there's no " |
750 |
|
|
"\":\" separator to indicate parameters", me, str); |
751 |
|
|
airMopError(mop); return 1; |
752 |
|
|
} |
753 |
|
|
*col = '\0'; |
754 |
|
|
etype = airEnumVal(pullEnergyType, str); |
755 |
|
|
if (pullEnergyTypeUnknown == etype) { |
756 |
|
|
biffAddf(PULL, "%s: didn't recognize \"%s\" as a %s", me, |
757 |
|
|
str, pullEnergyType->name); |
758 |
|
|
airMopError(mop); return 1; |
759 |
|
|
} |
760 |
|
|
|
761 |
|
|
ensp->energy = pullEnergyAll[etype]; |
762 |
|
|
if (0 == ensp->energy->parmNum) { |
763 |
|
|
biffAddf(PULL, "%s: \"%s\" energy has no parms, but got something", me, |
764 |
|
|
ensp->energy->name); |
765 |
|
|
return 1; |
766 |
|
|
} |
767 |
|
|
|
768 |
|
|
_pstr = pstr = col+1; |
769 |
|
|
/* code lifted from teem/src/nrrd/kernel.c, should probably refactor. . . */ |
770 |
|
|
for (haveParm=0; haveParm<ensp->energy->parmNum; haveParm++) { |
771 |
|
|
if (!pstr) { |
772 |
|
|
break; |
773 |
|
|
} |
774 |
|
|
if (1 != sscanf(pstr, "%lg", &pval)) { |
775 |
|
|
biffAddf(PULL, "%s: trouble parsing \"%s\" as double (in \"%s\")", |
776 |
|
|
me, _pstr, _str); |
777 |
|
|
airMopError(mop); return 1; |
778 |
|
|
} |
779 |
|
|
ensp->parm[haveParm] = pval; |
780 |
|
|
if ((pstr = strchr(pstr, ','))) { |
781 |
|
|
pstr++; |
782 |
|
|
if (!*pstr) { |
783 |
|
|
biffAddf(PULL, "%s: nothing after last comma in \"%s\" (in \"%s\")", |
784 |
|
|
me, _pstr, _str); |
785 |
|
|
airMopError(mop); return 1; |
786 |
|
|
} |
787 |
|
|
} |
788 |
|
|
} |
789 |
|
|
/* haveParm is now the number of parameters that were parsed. */ |
790 |
|
|
if (haveParm < ensp->energy->parmNum) { |
791 |
|
|
biffAddf(PULL, "%s: parsed only %u of %u required parms (for %s energy)" |
792 |
|
|
"from \"%s\" (in \"%s\")", |
793 |
|
|
me, haveParm, ensp->energy->parmNum, |
794 |
|
|
ensp->energy->name, _pstr, _str); |
795 |
|
|
airMopError(mop); return 1; |
796 |
|
|
} else { |
797 |
|
|
if (pstr) { |
798 |
|
|
biffAddf(PULL, "%s: \"%s\" (in \"%s\") has more than %u doubles", |
799 |
|
|
me, _pstr, _str, ensp->energy->parmNum); |
800 |
|
|
airMopError(mop); return 1; |
801 |
|
|
} |
802 |
|
|
} |
803 |
|
|
|
804 |
|
|
airMopOkay(mop); |
805 |
|
|
return 0; |
806 |
|
|
} |
807 |
|
|
|
808 |
|
|
int |
809 |
|
|
_pullHestEnergyParse(void *ptr, char *str, char err[AIR_STRLEN_HUGE]) { |
810 |
|
|
static const char me[]="_pullHestForceParse"; |
811 |
|
|
pullEnergySpec **enspP; |
812 |
|
|
char *perr; |
813 |
|
|
|
814 |
|
|
if (!(ptr && str)) { |
815 |
|
|
sprintf(err, "%s: got NULL pointer", me); |
816 |
|
|
return 1; |
817 |
|
|
} |
818 |
|
|
enspP = (pullEnergySpec **)ptr; |
819 |
|
|
*enspP = pullEnergySpecNew(); |
820 |
|
|
if (pullEnergySpecParse(*enspP, str)) { |
821 |
|
|
perr = biffGetDone(PULL); |
822 |
|
|
airStrcpy(err, AIR_STRLEN_HUGE, perr); |
823 |
|
|
free(perr); |
824 |
|
|
return 1; |
825 |
|
|
} |
826 |
|
|
return 0; |
827 |
|
|
} |
828 |
|
|
|
829 |
|
|
hestCB |
830 |
|
|
_pullHestEnergySpec = { |
831 |
|
|
sizeof(pullEnergySpec*), |
832 |
|
|
"energy specification", |
833 |
|
|
_pullHestEnergyParse, |
834 |
|
|
(airMopper)pullEnergySpecNix |
835 |
|
|
}; |
836 |
|
|
|
837 |
|
|
hestCB * |
838 |
|
|
pullHestEnergySpec = &_pullHestEnergySpec; |