GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: src/pull/energy.c Lines: 0 200 0.0 %
Date: 2017-05-26 Branches: 0 80 0.0 %

Line Branch Exec Source
1
/*
2
  Teem: Tools to process and visualize scientific data and images             .
3
  Copyright (C) 2013, 2012, 2011, 2010, 2009  University of Chicago
4
  Copyright (C) 2008, 2007, 2006, 2005  Gordon Kindlmann
5
  Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998  University of Utah
6
7
  This library is free software; you can redistribute it and/or
8
  modify it under the terms of the GNU Lesser General Public License
9
  (LGPL) as published by the Free Software Foundation; either
10
  version 2.1 of the License, or (at your option) any later version.
11
  The terms of redistributing and/or modifying this software also
12
  include exceptions to the LGPL that facilitate static linking.
13
14
  This library is distributed in the hope that it will be useful,
15
  but WITHOUT ANY WARRANTY; without even the implied warranty of
16
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17
  Lesser General Public License for more details.
18
19
  You should have received a copy of the GNU Lesser General Public License
20
  along with this library; if not, write to Free Software Foundation, Inc.,
21
  51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
22
*/
23
24
#include "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;