GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: src/bin/puller.c Lines: 0 250 0.0 %
Date: 2017-05-26 Branches: 0 148 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
/*
25
** NOTE: all the "#ifdef DEFT" directives refer to an FLTK2-based GUI
26
** for some parts of Teem called "Deft".  Unfortunately FLTK2 has been
27
** abandoned, and Deft is not released or supported in any way.  The
28
** Deft-related code is preserved for legacy purposes.
29
*/
30
31
#ifdef DEFT
32
#include "Deft.h"
33
#include "Contour.h"
34
#include "Viewer.h"
35
#include "ViewerUI.h"
36
#include "TensorGlyph.h"
37
#include "TensorGlyphUI.h"
38
#include "Slider.h"
39
#include "TriPlane.h"
40
#include "TriPlaneUI.h"
41
42
#endif
43
44
#include <teem/pull.h>
45
#include <teem/meet.h>
46
47
static const char *info =
48
  ("Command-line interface to the \"pull\" library. "
49
   "Published research using this tool or the \"pull\" library "
50
   "should cite the paper: \n "
51
   "\t\tGordon L. Kindlmann, Ra{\\'u}l San Jos{\\'e} Est{\\'e}par, Stephen M. Smith,\n "
52
   "\t\tCarl-Fredrik Westin. Sampling and Visualizing Creases with Scale-Space\n "
53
   "\t\tParticles. IEEE Trans. on Visualization and Computer Graphics,\n "
54
   "\t\t15(6):1415-1424 (2009).");
55
56
#ifdef DEFT
57
typedef struct {
58
  fltk::FloatInput *scaleVecInput[3];
59
  fltk::ValueInput *glyphScaleRadInput;
60
  Deft::Slider *isoval;
61
  Deft::Slider *strength;
62
  Deft::Slider *quality;
63
  Deft::Slider *alpha, *beta, *cwell, *gamma;
64
  /* Deft::Slider *height; */
65
  Deft::Slider *ccSelect, *rho, *sclMean, *sclWind;
66
  fltk::CheckButton *ccSingle;
67
  Deft::Contour *contour;
68
  Deft::Scene *scene;
69
  Deft::Viewer *viewer;
70
  Deft::TensorGlyph *glyph, *hedge;
71
  fltk::IntInput *iters;
72
  fltk::FloatInput *radius;
73
  fltk::ValueInput *verbose;
74
  pullContext *pctx;
75
  Nrrd *nPosOut, *nTenOut, *nFrcOut, *nten, *ntmp, *nenr, *nscl,
76
    *nidcc, *nstrn, *nqual, *ncovar, *ntcovar, *nstab, *nintern,
77
    *nstuck, *nfrcOld, *nfrcNew, *nposOld, *nposNew, *nrgb, *nccrgb,
78
    *ncval, *ncmap, *ncmapOut, *nblur;
79
  NrrdResampleContext *rsmc;
80
  const Nrrd *norig;
81
  NrrdRange *cvalRange;
82
  limnPolyData *phistLine, *phistTube;
83
  Deft::PolyData *phistSurf;
84
  double icvalr[2], sclMin, sclMax, strnMin, qualMin,
85
    scaleVec[3], glyphScaleRad, energyIncreasePermitFrac;
86
} pullBag;
87
88
void
89
verbose_cb(fltk::Widget *widget, pullBag *bag) {
90
  fltk::ValueInput *val;
91
92
  val = (fltk::ValueInput *)widget;
93
  pullVerboseSet(bag->pctx, (int)val->value());
94
}
95
/* ... DEFT ... */
96
void
97
isovalue_cb(fltk::Widget *widget, pullBag *bag) {
98
  Deft::Slider *slider;
99
100
  slider = (Deft::Slider *)widget;
101
  if (bag->contour) {
102
    bag->contour->extract(slider->value());
103
  }
104
  bag->viewer->redraw();
105
}
106
107
void
108
outputGet(pullBag *bag) {
109
  char me[]="outputGet", *err;
110
  size_t cropMin[2], cropMax[2];
111
112
  if (pullOutputGet(bag->nPosOut, bag->nTenOut,
113
                    bag->nstrn, /* may be NULL */
114
                    bag->scaleVec, bag->glyphScaleRad,
115
                    bag->pctx)
116
      || pullPropGet(bag->nscl, pullPropScale, bag->pctx)
117
      || (bag->pctx->ispec[pullInfoQuality]
118
          ? pullInfoGet(bag->nqual, pullInfoQuality, bag->pctx)
119
          : 0)
120
      || pullPropGet(bag->nenr, pullPropEnergy, bag->pctx)
121
      || pullPropGet(bag->nidcc, pullPropIdCC, bag->pctx)
122
      || pullPropGet(bag->nstuck, pullPropStuck, bag->pctx)
123
      || pullPropGet(bag->ncovar, pullPropNeighCovar7Ten, bag->pctx)
124
#if PULL_TANCOVAR
125
      || pullPropGet(bag->ntcovar, pullPropNeighTanCovar, bag->pctx)
126
#endif
127
      || pullPropGet(bag->nstab, pullPropStability, bag->pctx)
128
      || pullPropGet(bag->nintern, pullPropNeighInterNum, bag->pctx)
129
      || pullPropGet(bag->nFrcOut, pullPropForce, bag->pctx)
130
      || (pullPhistEnabled
131
          ? pullPositionHistoryGet(bag->phistLine, bag->pctx)
132
          : 0)) {
133
    err = biffGetDone(PULL);
134
    fprintf(stderr, "%s: error getting pull output:\n%s\n", me, err);
135
    free(err);
136
    exit(1);
137
  }
138
  /* ... DEFT ... */
139
  cropMin[0] = 0;
140
  cropMin[1] = 0;
141
  cropMax[0] = 2;
142
  cropMax[1] = bag->nPosOut->axis[1].size-1;
143
  if ((!bag->pctx->iter
144
       ? 0
145
       : (nrrdCopy(bag->nfrcOld, bag->nfrcNew)
146
          || nrrdCopy(bag->nposOld, bag->nposNew)))
147
      || nrrdConvert(bag->nten, bag->nTenOut, nrrdTypeFloat)
148
      /* hacks to visualize the (tan) covariance tensors
149
      || nrrdCopy(bag->nten, bag->ncovar)
150
      || nrrdCopy(bag->nten, bag->ntcovar)
151
      */
152
      || nrrdCrop(bag->ntmp, bag->nPosOut, cropMin, cropMax)
153
      || nrrdConvert(bag->nposNew, bag->ntmp, nrrdTypeFloat)
154
      || nrrdCrop(bag->ntmp, bag->nFrcOut, cropMin, cropMax)
155
      || nrrdConvert(bag->nfrcNew, bag->ntmp, nrrdTypeFloat)
156
      || (!bag->pctx->iter
157
          ? (nrrdCopy(bag->nfrcOld, bag->nfrcNew)
158
             || nrrdCopy(bag->nposOld, bag->nposNew))
159
          : 0)) {
160
    err = biffGetDone(NRRD);
161
    fprintf(stderr, "%s: another error 0\n%s\n", me, err);
162
    free(err);
163
    exit(1);
164
  }
165
}
166
/* ... DEFT ... */
167
void
168
outputShow(pullBag *bag) {
169
  char me[]="outputShow", *err;
170
  float *rgb;
171
  unsigned int ii, nn, *idcc;
172
  unsigned char *stuck;
173
  int first;
174
  double *cval, emean, estdv, *pos;
175
176
  /*
177
  if (limnPolyDataSpiralTubeWrap(bag->phistTube, bag->phistLine,
178
                                 (1 << limnPolyDataInfoRGBA)
179
                                 | (1 << limnPolyDataInfoNorm),
180
                                 NULL,
181
                                 8, 8, bag->glyph->glyphScale()/5)) {
182
    err = biffGetDone(LIMN);
183
    fprintf(stderr, "%s: another error 1\n%s\n", me, err);
184
    free(err);
185
    exit(1);
186
  }
187
  */
188
189
  if (pullPhistEnabled) {
190
    bag->phistSurf->changed();
191
  }
192
193
  bag->ncval = bag->nenr;
194
  /* bag->ncval = bag->nstrn; */
195
  /* bag->ncval = bag->nstuck; */
196
  /* bag->ncval = bag->nscl; */
197
198
  /* ... DEFT ... */
199
  if (bag->ncval) {
200
    nrrdRangeSet(bag->cvalRange, bag->ncval, AIR_FALSE);
201
  } else {
202
    bag->cvalRange->min = AIR_NAN;
203
    bag->cvalRange->max = AIR_NAN;
204
  }
205
  if (bag->ncval) {
206
    cval = AIR_CAST(double *, bag->ncval->data);
207
  } else {
208
    cval = NULL;
209
  }
210
  if (cval) {
211
    nn = bag->ncval->axis[0].size;
212
    emean = 0;
213
    for (ii=0; ii<nn; ii++) {
214
      emean += cval[ii];
215
    }
216
    emean /= nn;
217
    estdv = 0;
218
    for (ii=0; ii<nn; ii++) {
219
      estdv += (emean - cval[ii])*(emean - cval[ii]);
220
    }
221
    estdv = sqrt(estdv/nn);
222
    if (bag->cvalRange->hasNonExist) {
223
      fprintf(stderr, "!%s: cval range %g -- %g (%s), mean %g, stdv %g\n", me,
224
              bag->cvalRange->min, bag->cvalRange->max,
225
              bag->cvalRange->hasNonExist ? "HAS non-exist" : "no non-exist",
226
              emean, estdv);
227
    }
228
    bag->cvalRange->min = AIR_LERP(0.7, bag->cvalRange->min, emean - 2*estdv);
229
    bag->cvalRange->max = AIR_LERP(0.7, bag->cvalRange->max, emean + 2*estdv);
230
  }
231
  /* ... DEFT ... */
232
  float *cmapOut;
233
  if (bag->ncmap
234
      && bag->ncval
235
      && AIR_EXISTS(bag->cvalRange->min)
236
      && AIR_EXISTS(bag->cvalRange->max)) {
237
    /* double mmin, mmax; */
238
    fprintf(stderr, "!%s: cval cmap range %g %g ----------- \n", me,
239
            bag->cvalRange->min, bag->cvalRange->max);
240
    /*
241
    mmin = -0.0342937;
242
    mmax = -0.0105725;
243
    */
244
    /* */
245
    /*
246
    bag->cvalRange->min = AIR_LERP(0.05, mmin, mmax);
247
    bag->cvalRange->max = AIR_LERP(0.3, mmin, mmax);
248
    */
249
    /* */
250
    if (nrrdApply1DRegMap(bag->ncmapOut, bag->ncval, bag->cvalRange,
251
                          bag->ncmap, nrrdTypeFloat, AIR_TRUE)) {
252
      err = biffGetDone(NRRD);
253
      fprintf(stderr, "%s: cmap error\n%s\n", me, err);
254
      free(err);  exit(1);
255
    }
256
    cmapOut = AIR_CAST(float *, bag->ncmapOut->data);
257
  } else {
258
    cmapOut = NULL;
259
  }
260
261
  /* ... DEFT ... */
262
263
  idcc = AIR_CAST(unsigned int *, bag->nidcc->data);
264
  stuck = AIR_CAST(unsigned char *, bag->nstuck->data);
265
  nn = bag->nPosOut->axis[1].size;
266
  pos = AIR_CAST(double *, bag->nPosOut->data);
267
268
  first = bag->nrgb->axis[1].size != bag->nPosOut->axis[1].size;
269
  /*
270
  fprintf(stderr, "!%s: %u %u -> %d\n", me,
271
          AIR_UINT(bag->nrgb->axis[1].size),
272
          AIR_UINT(bag->nPosOut->axis[1].size), first);
273
  */
274
  if (first) {
275
    if (nrrdMaybeAlloc_va(bag->nrgb, nrrdTypeFloat, 2,
276
                          AIR_CAST(size_t, 3),
277
                          bag->nPosOut->axis[1].size)) {
278
      err = biffGetDone(NRRD);
279
      fprintf(stderr, "%s: error creating RGB:\n%s\n", me, err);
280
      free(err);
281
      exit(1);
282
    }
283
  }
284
  bag->icvalr[0] = bag->cvalRange->min;
285
  bag->icvalr[1] = bag->cvalRange->max;
286
287
  /* ... DEFT ... */
288
289
  double *strnOut;
290
  strnOut = (bag->nstrn
291
             ? AIR_CAST(double *, bag->nstrn->data)
292
             : NULL);
293
  double *qualOut;
294
  qualOut = (bag->nqual
295
             ? AIR_CAST(double *, bag->nqual->data)
296
             : NULL);
297
  rgb = (float*)bag->nrgb->data;
298
  for (ii=0; ii<nn; ii++) {
299
    float ee, *ccrgb;
300
    ccrgb = (float*)bag->nccrgb->data;
301
    /* ee = bag->cvalRange->min - (bag->icvalr[1] - bag->icvalr[0])/50;*/
302
    ee = bag->cvalRange->min;
303
    if (bag->pctx->CCNum && ccrgb) {
304
      rgb[0 + 3*ii] = ccrgb[0 + 3*idcc[ii]];
305
      rgb[1 + 3*ii] = ccrgb[1 + 3*idcc[ii]];
306
      rgb[2 + 3*ii] = ccrgb[2 + 3*idcc[ii]];
307
    } else if (cmapOut) {
308
      ELL_3V_COPY(rgb + 3*ii, cmapOut + 3*ii);
309
    } else {
310
      ELL_3V_SET(rgb + 3*ii, 0.95, 0.95, 0.95);
311
      /*
312
      if (AIR_EXISTS(cval[ii])) {
313
        rgb[1 + 3*ii] = AIR_AFFINE(ee, cval[ii], bag->cvalRange->max, 0, 1);
314
        rgb[1 + 3*ii] = sqrt(rgb[1 + 3*ii]);
315
        rgb[0 + 3*ii] = rgb[1 + 3*ii];
316
      } else {
317
        rgb[1 + 3*ii] = 0;
318
        rgb[0 + 3*ii] = 1;
319
      }
320
      rgb[2 + 3*ii] = stuck[ii];
321
      */
322
    }
323
  }
324
325
  /* ... DEFT ... */
326
327
  if (1) {
328
    float *ten, *pos;
329
    double *posOut;
330
    ten = AIR_CAST(float *, bag->nten->data);
331
    posOut = AIR_CAST(double *, bag->nPosOut->data);
332
    pos = AIR_CAST(float *, bag->nposNew->data);
333
    for (ii=0; ii<nn; ii++) {
334
      if (!( AIR_IN_CL(bag->sclMin-0.00001, posOut[3],
335
                       bag->sclMax+0.00001) )) {
336
        ten[0] = 0;
337
      } else if (strnOut && strnOut[ii] < bag->strnMin) {
338
        ten[0] = 0;
339
      } else if (qualOut && qualOut[ii] < bag->qualMin-0.000001) {
340
        ten[0] = 0;
341
      } else if (bag->pctx->CCNum
342
                 && (bag->ccSingle->value()
343
                     ? idcc[ii] != bag->ccSelect->value()
344
                     : idcc[ii] > bag->ccSelect->value())) {
345
        ten[0] = 0;
346
      } else {
347
        ten[0] = 1;
348
      }
349
      /*
350
      if (4589 == ii) {
351
        fprintf(stderr, "!%s: point %u/%u at (%g,%g,%g)=(%g,%g,%g,%g) got ten[0] %g\n",
352
                me, ii, nn,
353
                pos[0], pos[1], pos[1],
354
                posOut[0], posOut[1], posOut[2], posOut[3],
355
                ten[0]);
356
      }
357
      */
358
      pos += 3;
359
      posOut += 4;
360
      ten += 7;
361
    }
362
  }
363
364
  /* ... DEFT ... */
365
366
  if (bag->nPosOut->axis[1].size) {
367
    bag->glyph->dataSet(bag->nPosOut->axis[1].size,
368
                        (float*)bag->nten->data, 7,
369
                        (float*)bag->nposNew->data, 3, rgb, 3, NULL);
370
    bag->glyph->update();
371
    /*
372
    bag->hedge->dataSet(bag->nPosOut->axis[1].size,
373
                        (float*)bag->nfrcNew->data, 3,
374
                        (float*)bag->nposOld->data, 3, rgb, 3, NULL);
375
    bag->hedge->update();
376
    */
377
    bag->viewer->redraw();
378
    fltk::flush();
379
  } else {
380
    fprintf(stderr, "!%s: got zero tensors out!\n", me);
381
  }
382
  return;
383
}
384
385
void
386
iter_cb(void *_bag) {
387
  pullBag *bag;
388
  bag = AIR_CAST(pullBag *, _bag);
389
  outputGet(bag);
390
  outputShow(bag);
391
}
392
393
/* ... DEFT ... */
394
395
void
396
step_cb(fltk::Widget *, pullBag *bag) {
397
  /*  static double lastthresh = -42; */
398
  char me[]="step_cb", *err;
399
  static unsigned int itersTotal=0;
400
401
  unsigned int iters = bag->iters->ivalue();
402
  bag->pctx->iterParm.max += iters;
403
  if (pullRun(bag->pctx)) {
404
    err = biffGetDone(PULL);
405
    fprintf(stderr, "%s: error running pull:\n%s\n", me, err);
406
    free(err);
407
    exit(1);
408
  }
409
  itersTotal += iters;
410
  fprintf(stderr, "!%s: enr = %g; time = %g sec; %u iters (%g iters/sec)\n",
411
          me, bag->pctx->energy, bag->pctx->timeRun, itersTotal,
412
          itersTotal/bag->pctx->timeRun);
413
  outputGet(bag);
414
  outputShow(bag);
415
  for (unsigned int ci=pullCountUnknown+1; ci<pullCountLast; ci++) {
416
    if (bag->pctx->count[ci]) {
417
      fprintf(stderr, "  %u: %s\n", bag->pctx->count[ci],
418
              airEnumStr(pullCount, ci));
419
    }
420
  }
421
}
422
423
/* ... DEFT ... */
424
425
void
426
gammaSet_cb(fltk::Widget *, pullBag *bag) {
427
  char me[]="gammaSet_cb";
428
429
  if (pullGammaLearn(bag->pctx)) {
430
    char *err = biffGetDone(PULL);
431
    fprintf(stderr, "%s: problem learning gamma:\n%s", me, err);
432
    free(err);
433
  }
434
  if (bag->pctx->sysParm.gamma > bag->gamma->maximum()) {
435
    bag->gamma->maximum(2*bag->pctx->sysParm.gamma);
436
  }
437
  bag->gamma->value(bag->pctx->sysParm.gamma);
438
}
439
440
/* ... DEFT ... */
441
442
void
443
cc_cb(fltk::Widget *, pullBag *bag) {
444
  char me[]="cc_cb";
445
  unsigned int cc;
446
  float *rgb;
447
448
  if (pullCCFind(bag->pctx)
449
      || pullCCSort(bag->pctx,
450
                    (bag->pctx->ispec[pullInfoQuality]
451
                     ? pullInfoQuality
452
                     : 0), bag->rho->value())) {
453
    char *err = biffGetDone(PULL);
454
    fprintf(stderr, "%s: problem finding/sorting CCs:\n%s", me, err);
455
    free(err);
456
  }
457
  printf("%s: found %u CCs\n", me, bag->pctx->CCNum);
458
  bag->ccSelect->range(0, bag->pctx->CCNum-1);
459
  if (bag->nccrgb->axis[1].size != bag->pctx->CCNum) {
460
    airSrandMT(AIR_UINT(airTime()));
461
    if (nrrdMaybeAlloc_va(bag->nccrgb, nrrdTypeFloat, 2,
462
                          AIR_CAST(size_t, 3),
463
                          AIR_CAST(size_t, bag->pctx->CCNum))) {
464
      char *err = biffGetDone(NRRD);
465
      fprintf(stderr, "%s: problem alloc'ing cc rgb:\n%s", me, err);
466
      free(err);
467
    }
468
    rgb = (float*)bag->nccrgb->data;
469
    ELL_3V_SET(rgb + 0*3, 0.95, 0.95, 0.95);
470
    for (cc=0; cc<bag->pctx->CCNum; cc++) {
471
      rgb[0 + 3*cc] = AIR_AFFINE(0, airDrandMT(), 1, 0.3, 1.0);
472
      rgb[1 + 3*cc] = AIR_AFFINE(0, airDrandMT(), 1, 0.3, 1.0);
473
      rgb[2 + 3*cc] = AIR_AFFINE(0, airDrandMT(), 1, 0.3, 1.0);
474
    }
475
  }
476
  outputGet(bag);
477
  outputShow(bag);
478
}
479
480
/* ... DEFT ... */
481
482
void
483
ccSelect_cb(fltk::Widget *, pullBag *bag) {
484
485
  outputShow(bag);
486
}
487
488
void
489
scaleGlyph_cb(fltk::Widget *, pullBag *bag) {
490
491
  bag->scaleVec[0] = bag->scaleVecInput[0]->fvalue();
492
  bag->scaleVec[1] = bag->scaleVecInput[1]->fvalue();
493
  bag->scaleVec[2] = bag->scaleVecInput[2]->fvalue();
494
  bag->glyphScaleRad = bag->glyphScaleRadInput->value();
495
  outputGet(bag);
496
  outputShow(bag);
497
}
498
499
/* ... DEFT ... */
500
501
void
502
reblur_cb(fltk::Widget *, pullBag *bag) {
503
  static const char me[]="reblur_cb";
504
  double kparm[NRRD_KERNEL_PARMS_NUM], scl;
505
  int E;
506
507
  if (!bag->pctx->haveScale) {
508
    return;
509
  }
510
  scl = bag->sclMean->value();
511
  if (bag->pctx->flag.scaleIsTau) {
512
    kparm[0] = gageSigOfTau(scl);
513
    printf("!%s: tau = %g ---> sigma = %g\n", me, scl, kparm[0]);
514
  } else {
515
    kparm[0] = scl;
516
    printf("!%s: sigma = %g\n", me, kparm[0]);
517
  }
518
  kparm[1] = 3;
519
  E = 0;
520
  for (unsigned int axi=0; axi<3; axi++) {
521
    if (!E) E |= nrrdResampleKernelSet(bag->rsmc, axi,
522
                                       nrrdKernelDiscreteGaussian,
523
                                       kparm);
524
  }
525
  if (!E) E |= nrrdResampleExecute(bag->rsmc, bag->nblur);
526
  if (E) {
527
    char *err = biffGetDone(NRRD);
528
    fprintf(stderr, "%s: problem resampling to scale %g:\n%s",
529
            me, bag->sclMean->value(), err);
530
    free(err);
531
  }
532
  outputShow(bag);
533
  return;
534
}
535
536
/* ... DEFT ... */
537
538
void
539
scale_cb(fltk::Widget *, pullBag *bag) {
540
  double sclMean, sclWind;
541
542
  if (bag->pctx->haveScale) {
543
    sclMean = bag->sclMean->value();
544
    sclWind = bag->sclWind->value();
545
    bag->sclMin = sclMean - sclWind/2;
546
    bag->sclMax = sclMean + sclWind/2;
547
  } else {
548
    bag->sclMin = 0;
549
    bag->sclMax = 0;
550
  }
551
  outputShow(bag);
552
  return;
553
}
554
555
void
556
alpha_cb(fltk::Widget *, pullBag *bag) {
557
558
  pullSysParmSet(bag->pctx, pullSysParmAlpha, bag->alpha->value());
559
}
560
561
void
562
beta_cb(fltk::Widget *, pullBag *bag) {
563
564
  pullSysParmSet(bag->pctx, pullSysParmBeta, bag->beta->value());
565
}
566
567
/* ... DEFT ... */
568
569
void
570
cwell_cb(fltk::Widget *, pullBag *bag) {
571
  double *parm;
572
573
  parm = bag->pctx->energySpecR->parm;
574
  parm[1] = bag->cwell->value();
575
  pullSysParmSet(bag->pctx, pullSysParmEnergyIncreasePermit,
576
                 bag->energyIncreasePermitFrac*bag->cwell->value());
577
  {
578
    unsigned int ii, nn;
579
    double xx, yy, de;
580
    FILE *file;
581
582
    if ((file = fopen("eplot.txt", "w"))) {
583
      nn = 800;
584
      for (ii=0; ii<nn; ii++) {
585
        xx = AIR_AFFINE(0, ii, nn-1, 0.0, 1.0);
586
        yy = bag->pctx->energySpecR->energy->eval(&de, xx, parm);
587
        fprintf(file, "%f %f\n", xx, yy);
588
      }
589
      fclose(file);
590
    }
591
  }
592
}
593
594
void
595
gamma_cb(fltk::Widget *, pullBag *bag) {
596
597
  pullSysParmSet(bag->pctx, pullSysParmGamma, bag->gamma->value());
598
}
599
600
/* ... DEFT ... */
601
602
void
603
strength_cb(fltk::Widget *, pullBag *bag) {
604
605
  bag->strnMin = bag->strength->value();
606
  outputShow(bag);
607
}
608
609
void
610
quality_cb(fltk::Widget *, pullBag *bag) {
611
612
  bag->qualMin = bag->quality->value();
613
  outputShow(bag);
614
}
615
616
void
617
save_cb(fltk::Widget *, pullBag *bag) {
618
  static const char me[]="save_cb";
619
  unsigned int ii, nn, count;
620
  float *ten;
621
  Nrrd *nPosSel, *nStrnSel;
622
  double *posSel, *posAll, *strnSel, *strnAll;
623
624
  if (!( 0.0 == ELL_3V_LEN(bag->scaleVec) )) {
625
    fprintf(stderr, "%s: refusing to save with non-zero scaleVec %g %g %g\n",
626
            me, bag->scaleVec[0], bag->scaleVec[1], bag->scaleVec[2]);
627
    return;
628
  }
629
630
  /* ... DEFT ... */
631
632
  nPosSel = nrrdNew();
633
  if (bag->nstrn) {
634
    nStrnSel = nrrdNew();
635
  } else {
636
    nStrnSel = NULL;
637
  }
638
  outputGet(bag);
639
  outputShow(bag); /* will exploit its masking of ten[0] */
640
  count = 0;
641
  nn = bag->nPosOut->axis[1].size;
642
  ten = AIR_CAST(float *, bag->nten->data);
643
  for (ii=0; ii<nn; ii++) {
644
    count += !!ten[0];
645
    ten += 7;
646
  }
647
  nrrdMaybeAlloc_va(nPosSel, nrrdTypeDouble, 2,
648
                    AIR_CAST(size_t, 4),
649
                    AIR_CAST(size_t, count));
650
  posAll = AIR_CAST(double *, bag->nPosOut->data);
651
  posSel = AIR_CAST(double *, nPosSel->data);
652
  if (bag->nstrn) {
653
    nrrdMaybeAlloc_va(nStrnSel, nrrdTypeDouble, 1,
654
                      AIR_CAST(size_t, count));
655
    strnAll = AIR_CAST(double *, bag->nstrn->data);
656
    strnSel = AIR_CAST(double *, nStrnSel->data);
657
  } else {
658
    strnSel = NULL;
659
    strnAll = NULL;
660
  }
661
  ten = AIR_CAST(float *, bag->nten->data);
662
  count = 0;
663
  for (ii=0; ii<nn; ii++) {
664
    if (ten[0]) {
665
      ELL_4V_COPY(posSel, posAll);
666
      if (strnSel && strnAll) {
667
        strnSel[count] = strnAll[ii];
668
      }
669
      posSel += 4;
670
      count++;
671
    }
672
    posAll += 4;
673
    ten += 7;
674
  }
675
676
  /* ... DEFT ... */
677
678
  nrrdSave("pos-all.nrrd", bag->nPosOut, NULL);
679
  if (0) {
680
    nrrdSave("pos-sel.nrrd", nPosSel, NULL);
681
  } else {
682
    char fname[512];
683
    FILE *ff;
684
    unsigned int ii=0;
685
686
    for (ii=0, ff=NULL; AIR_TRUE; ii++) {
687
      ff = airFclose(ff);
688
      sprintf(fname, "pos-sel-%03u.nrrd", ii);
689
      if (!(ff = fopen(fname, "rb"))) {
690
        break;
691
      }
692
    }
693
    ff = airFclose(ff);
694
    nrrdSave(fname, nPosSel, NULL);
695
  }
696
  nrrdSave("covar-all.nrrd", bag->ncovar, NULL);
697
#if PULL_TANCOVAR
698
  nrrdSave("tcovar-all.nrrd", bag->ntcovar, NULL);
699
#endif
700
  nrrdSave("stab-all.nrrd", bag->nstab, NULL);
701
  nrrdSave("intern-all.nrrd", bag->nintern, NULL);
702
  if (bag->nstrn) {
703
    nrrdSave("strn-all.nrrd", bag->nstrn, NULL);
704
    nrrdSave("strn-sel.nrrd", nStrnSel, NULL);
705
  }
706
  nrrdNuke(nPosSel);
707
  if (bag->nstrn) {
708
    nrrdNuke(nStrnSel);
709
  }
710
  return;
711
}
712
713
#endif  /* DEFT */
714
715
int
716
main(int argc, const char **argv) {
717
  hestOpt *hopt=NULL;
718
  hestParm *hparm;
719
  airArray *mop;
720
  const char *me;
721
722
#ifdef DEFT
723
  float fr[3], at[3], up[3], fovy, neer, faar, dist, bg[3];
724
  int imgSize[2], ortho, rght, atrel, camkeep;
725
726
  float anisoThresh, anisoThreshMin, glyphScale, haloTraceBound,
727
    glyphHaloWidth, glyphNormPow, glyphEvalPow, sqdSharp;
728
  int glyphType, glyphFacetRes, aniso;
729
  double ssrange[2];
730
  pullBag bag;
731
#endif
732
733
  char *err, *outS, *extraOutBaseS, *addLogS, *cachePathSS;
734
  FILE *addLog;
735
  meetPullVol **vspec;
736
  meetPullInfo **idef;
737
  Nrrd *nPosIn=NULL, *nPosOut;
738
  pullEnergySpec *enspR, *enspS, *enspWin;
739
  NrrdKernelSpec *k00, *k11, *k22, *kSSrecon, *kSSblur;
740
  NrrdBoundarySpec *bspec;
741
  pullContext *pctx;
742
  int E=0, ret=0;
743
  unsigned int vspecNum, idefNum;
744
  double scaleVec[3], glyphScaleRad;
745
  /* things that used to be set directly inside pullContext */
746
  int energyFromStrength, nixAtVolumeEdgeSpace, constraintBeforeSeedThresh,
747
    binSingle, liveThresholdOnInit, permuteOnRebin, noPopCntlWithZeroAlpha,
748
    useBetaForGammaLearn, restrictiveAddToBins, noAdd, unequalShapesAllow,
749
    popCntlEnoughTest, convergenceIgnoresPopCntl, zeroZ;
750
  int verbose;
751
  int interType, allowCodimension3Constraints, scaleIsTau, useHalton,
752
    pointPerVoxel;
753
  unsigned int samplesAlongScaleNum, pointNumInitial,
754
    ppvZRange[2], snap, iterMax, stuckIterMax, constraintIterMax,
755
    popCntlPeriod, addDescent, iterCallback, rngSeed, progressBinMod,
756
    threadNum, eipHalfLife, kssOpi, kssFinished, bspOpi, bspFinished;
757
  double jitter, stepInitial, constraintStepMin, radiusSpace, binWidthSpace,
758
    radiusScale, alpha, beta, _gamma, wall, energyIncreasePermit,
759
    backStepScale, opporStepScale, energyDecreaseMin, energyDecreasePopCntlMin,
760
    neighborTrueProb, probeProb, fracNeighNixedMax;
761
762
  mop = airMopNew();
763
  hparm = hestParmNew();
764
  airMopAdd(mop, hparm, (airMopper)hestParmFree, airMopAlways);
765
766
  nPosOut = nrrdNew();
767
  airMopAdd(mop, nPosOut, (airMopper)nrrdNuke, airMopAlways);
768
769
  hparm->respFileEnable = AIR_TRUE;
770
  me = argv[0];
771
772
#ifdef DEFT
773
  unsigned int baryRes;
774
  int saveAndQuit, fog;
775
776
  hestOptAdd(&hopt, "csqvmm", "min max", airTypeDouble, 2, 2,
777
             Deft::colorSclQuantityValueMinMax, "nan nan",
778
             "min/max values for cutting planes of scalar values");
779
  hestOptAdd(&hopt, "saq", "save & quit", airTypeInt, 0, 0, &saveAndQuit,
780
             NULL, "save image and quit, for batch processing");
781
  hestOptAdd(&hopt, "cmap", "nin", airTypeOther, 1, 1, &(bag.ncmap), "",
782
             "colormap for particles", NULL, NULL, nrrdHestNrrd);
783
  hestOptAdd(&hopt, "fr", "from point", airTypeFloat, 3, 3, fr, "3 4 5",
784
             "position of camera, used to determine view vector");
785
  hestOptAdd(&hopt, "at", "at point", airTypeFloat, 3, 3, at, "0 0 0",
786
             "camera look-at point, used to determine view vector");
787
  hestOptAdd(&hopt, "up", "up vector", airTypeFloat, 3, 3, up, "0 0 1",
788
             "camera pseudo-up vector, used to determine view coordinates");
789
  hestOptAdd(&hopt, "rh", NULL, airTypeInt, 0, 0, &rght, NULL,
790
             "normally, use a right-handed UVN frame (V points down), "
791
             "but in Deft this is always true");
792
  hestOptAdd(&hopt, "fv", "fov", airTypeFloat, 1, 1, &fovy, "20",
793
             "vertical field-of-view, in degrees");
794
  hestOptAdd(&hopt, "or", NULL, airTypeInt, 0, 0, &ortho, NULL,
795
             "use orthogonal projection instead of perspective");
796
  hestOptAdd(&hopt, "dn", "near clip", airTypeFloat, 1, 1, &neer, "-2",
797
             "position of near clipping plane, relative to look-at point");
798
  hestOptAdd(&hopt, "di", "image", airTypeFloat, 1, 1, &dist, "0.0",
799
             "normally, distance to image plane, "
800
             "but in Deft this is always 0.0");
801
  hestOptAdd(&hopt, "df", "far clip", airTypeFloat, 1, 1, &faar, "2",
802
             "position of far clipping plane, relative to look-at point");
803
  hestOptAdd(&hopt, "ar", NULL, airTypeInt, 0, 0, &atrel, NULL,
804
             "normally: near, image, and far plane distances are relative to "
805
             "the *at* point, instead of the eye point, "
806
             "but for Deft, this is always true");
807
  hestOptAdd(&hopt, "usecam", NULL, airTypeInt, 0, 0, &camkeep, NULL,
808
             "hack: by default, a camera reset is done to put the volume "
809
             "in view. Use this to say that the camera specified by the "
810
             "flags above should be preserved and used");
811
  hestOptAdd(&hopt, "bg", "R G B", airTypeFloat, 3, 3, bg, "0.2 0.3 0.4",
812
             "background color");
813
  hestOptAdd(&hopt, "fog", NULL, airTypeInt, 0, 0, &fog, NULL,
814
             "hack: turn on fog");
815
  hestOptAdd(&hopt, "is", "su sv", airTypeInt, 2, 2, imgSize, "640 480",
816
             "initial window size");
817
  /* ... DEFT ... */
818
  /* this tensor stuff is here because we're hijacking the tensor glyph
819
     object for doing the particle display ... */
820
  hestOptAdd(&hopt, "a", "aniso", airTypeEnum, 1, 1, &aniso, NULL,
821
             "anisotropy metric to make volume of",
822
             NULL, tenAniso);
823
  hestOptAdd(&hopt, "atr", "aniso thresh", airTypeFloat, 1, 1,
824
             &anisoThresh, "0.85",
825
             "Glyphs will be drawn only for tensors with anisotropy "
826
             "greater than this threshold");
827
  hestOptAdd(&hopt, "atrm", "aniso thresh min", airTypeFloat, 1, 1,
828
             &anisoThreshMin, "0.4",
829
             "lower bound on aniso thresh");
830
  hestOptAdd(&hopt, "g", "glyph shape", airTypeEnum, 1, 1, &glyphType, "sqd",
831
             "shape of glyph to use for display.  Possibilities "
832
             "include \"box\", \"sphere\"=\"sph\", \"cylinder\"=\"cyl\", and "
833
             "\"superquad\"=\"sqd\"", NULL, tenGlyphType);
834
  hestOptAdd(&hopt, "gsc", "scale", airTypeFloat, 1, 1, &glyphScale,
835
             "0.25", "over-all glyph size");
836
  hestOptAdd(&hopt, "htb", "trace", airTypeFloat, 1, 1, &haloTraceBound,
837
             "1.0", "halo trace bound");
838
  hestOptAdd(&hopt, "ghw", "hwidth", airTypeFloat, 1, 1, &glyphHaloWidth,
839
             "0.0", "glyph halo width");
840
  hestOptAdd(&hopt, "gnp", "npow", airTypeFloat, 1, 1, &glyphNormPow,
841
             "1.0", "pow() exponent for compressing range of norms");
842
  hestOptAdd(&hopt, "gep", "epow", airTypeFloat, 1, 1, &glyphEvalPow,
843
             "1.0", "pow() exponent for compressing single eigenvalues");
844
  hestOptAdd(&hopt, "br", "barycentric res", airTypeInt, 1, 1, &baryRes, "50",
845
             "resolution of sampling of tensor shape palette");
846
  hestOptAdd(&hopt, "gr", "glyph res", airTypeInt, 1, 1, &glyphFacetRes, "7",
847
             "resolution of polygonalization of glyphs (other than box)");
848
  hestOptAdd(&hopt, "sh", "sharpness", airTypeFloat, 1, 1, &sqdSharp, "2.5",
849
             "for superquadric glyphs, how much to sharp edges form as a "
850
             "function of differences between eigenvalues.  Higher values "
851
             "mean that edges form more easily");
852
#endif /* DEFT */
853
854
  hestOptAdd(&hopt, "int", "int", airTypeEnum, 1, 1, &interType,
855
             "justr", "inter-particle energy type", NULL, pullInterType);
856
  hestOptAdd(&hopt, "enr", "spec", airTypeOther, 1, 1, &enspR, "cotan",
857
             "inter-particle energy, radial component",
858
             NULL, NULL, pullHestEnergySpec);
859
  hestOptAdd(&hopt, "ens", "spec", airTypeOther, 1, 1, &enspS, "zero",
860
             "inter-particle energy, scale component",
861
             NULL, NULL, pullHestEnergySpec);
862
  hestOptAdd(&hopt, "enw", "spec", airTypeOther, 1, 1, &enspWin,
863
             "butter:16,0.8", "windowing to create locality with additive "
864
             "scale-space interaction (\"-int add\")",
865
             NULL, NULL, pullHestEnergySpec);
866
  hestOptAdd(&hopt, "zz", "bool", airTypeBool, 1, 1, &zeroZ, "false",
867
             "always constrain Z=0, to process 2D images");
868
  hestOptAdd(&hopt, "efs", "bool", airTypeBool, 1, 1,
869
             &energyFromStrength, "false",
870
             "whether or not strength contributes to particle-image energy");
871
  hestOptAdd(&hopt, "nave", "bool", airTypeBool, 1, 1,
872
             &nixAtVolumeEdgeSpace, "false",
873
             "whether or not to nix points at edge of volume, where gage had "
874
             "to invent values for kernel support");
875
  hestOptAdd(&hopt, "cbst", "bool", airTypeBool, 1, 1,
876
             &constraintBeforeSeedThresh, "false",
877
             "during initialization, try constraint satisfaction before "
878
             "testing seedThresh");
879
  hestOptAdd(&hopt, "noadd", NULL, airTypeBool, 0, 0,
880
             &noAdd, NULL, "turn off adding during population control");
881
  hestOptAdd(&hopt, "usa", "bool", airTypeBool, 1, 1,
882
             &unequalShapesAllow, "false",
883
             "allow volumes to have different shapes (false is safe as "
884
             "different volume sizes are often accidental)");
885
  hestOptAdd(&hopt, "pcet", "bool", airTypeBool, 1, 1, &popCntlEnoughTest,
886
             "true", "use neighbor-counting \"enough\" heuristic to "
887
             "bail out of pop cntl");
888
  hestOptAdd(&hopt, "cipc", "bool", airTypeBool, 1, 1, &convergenceIgnoresPopCntl,
889
             "false", "convergence test doesn't care if there has been "
890
             "recent changes due to population control");
891
  hestOptAdd(&hopt, "nobin", NULL, airTypeBool, 0, 0,
892
             &binSingle, NULL,
893
             "turn off spatial binning (which prevents multi-threading "
894
             "from being useful), for debugging or speed-up measurement");
895
  hestOptAdd(&hopt, "lti", "bool", airTypeBool, 1, 1,
896
             &liveThresholdOnInit, "true",
897
             "impose liveThresh on initialization");
898
  hestOptAdd(&hopt, "por", "bool", airTypeBool, 1, 1,
899
             &permuteOnRebin, "true",
900
             "permute points during rebinning");
901
  hestOptAdd(&hopt, "npcwza", "bool", airTypeBool, 1, 1,
902
             &noPopCntlWithZeroAlpha, "false",
903
             "no pop cntl with zero alpha");
904
  hestOptAdd(&hopt, "ubfgl", "bool", airTypeBool, 1, 1,
905
             &useBetaForGammaLearn, "false",
906
             "use beta for gamma learning");
907
  hestOptAdd(&hopt, "ratb", "bool", airTypeBool, 1, 1,
908
             &restrictiveAddToBins, "true",
909
             "be choosy when adding points to bins to avoid overlap");
910
  hestOptAdd(&hopt, "svec", "vec", airTypeDouble, 3, 3, scaleVec, "0 0 0",
911
             "if non-zero (length), vector to use for displaying scale "
912
             "in 3-space");
913
  hestOptAdd(&hopt, "gssr", "rad", airTypeDouble, 1, 1, &glyphScaleRad, "0.0",
914
             "if non-zero (length), scaling of scale to cylindrical tensors");
915
  hestOptAdd(&hopt, "v", "verbosity", airTypeInt, 1, 1, &verbose, "1",
916
             "verbosity level");
917
  hestOptAdd(&hopt, "vol", "vol0 vol1", airTypeOther, 1, -1, &vspec, NULL,
918
             "input volumes, in format <filename>:<kind>:<volname>",
919
             &vspecNum, NULL, meetHestPullVol);
920
  hestOptAdd(&hopt, "info", "info0 info1", airTypeOther, 1, -1, &idef, NULL,
921
             "info definitions, in format "
922
             "<info>[-c]:<volname>:<item>[:<zero>:<scale>]",
923
             &idefNum, NULL, meetHestPullInfo);
924
925
  hestOptAdd(&hopt, "k00", "kern00", airTypeOther, 1, 1, &k00,
926
             "cubic:1,0", "kernel for gageKernel00",
927
             NULL, NULL, nrrdHestKernelSpec);
928
  hestOptAdd(&hopt, "k11", "kern11", airTypeOther, 1, 1, &k11,
929
             "cubicd:1,0", "kernel for gageKernel11",
930
             NULL, NULL, nrrdHestKernelSpec);
931
  hestOptAdd(&hopt, "k22", "kern22", airTypeOther, 1, 1, &k22,
932
             "cubicdd:1,0", "kernel for gageKernel22",
933
             NULL, NULL, nrrdHestKernelSpec);
934
935
  hestOptAdd(&hopt, "sscp", "path", airTypeString, 1, 1, &cachePathSS, "./",
936
             "path (without trailing /) for where to read/write "
937
             "pre-blurred volumes for scale-space");
938
  kssOpi =
939
  hestOptAdd(&hopt, "kssb", "kernel", airTypeOther, 1, 1, &kSSblur,
940
             "dgauss:1,5", "default blurring kernel, to sample scale space",
941
             NULL, NULL, nrrdHestKernelSpec);
942
  bspOpi =
943
  hestOptAdd(&hopt, "bsp", "boundary", airTypeOther, 1, 1, &bspec,
944
             "wrap", "default boundary behavior of scale-space blurring",
945
             NULL, NULL, nrrdHestBoundarySpec);
946
  hestOptAdd(&hopt, "kssr", "kernel", airTypeOther, 1, 1, &kSSrecon,
947
             "hermite", "kernel for reconstructing from scale space samples",
948
             NULL, NULL, nrrdHestKernelSpec);
949
  hestOptAdd(&hopt, "nss", "# scl smpls", airTypeUInt, 1, 1,
950
             &samplesAlongScaleNum, "1",
951
             "if using \"-ppv\", number of samples along scale axis "
952
             "for each spatial position");
953
954
  hestOptAdd(&hopt, "np", "# points", airTypeUInt, 1, 1,
955
             &pointNumInitial, "1000",
956
             "number of points to start in system");
957
  hestOptAdd(&hopt, "halton", NULL, airTypeInt, 0, 0,
958
             &useHalton, NULL,
959
             "use Halton sequence initialization instead of "
960
             "uniform random");
961
  hestOptAdd(&hopt, "ppv", "# pnts/vox", airTypeInt, 1, 1,
962
             &pointPerVoxel, "0",
963
             "number of points per voxel to start in simulation "
964
             "(need to have a seed thresh vol, overrides \"-np\")");
965
  hestOptAdd(&hopt, "ppvzr", "z range", airTypeUInt, 2, 2,
966
             ppvZRange, "1 0",
967
             "range of Z slices (1st num < 2nd num) to do ppv in, or, "
968
             "\"1 0\" for whole volume");
969
  hestOptAdd(&hopt, "jit", "jitter", airTypeDouble, 1, 1,
970
             &jitter, "0",
971
             "amount of jittering to do with ppv");
972
  hestOptAdd(&hopt, "pi", "npos", airTypeOther, 1, 1, &nPosIn, "",
973
             "4-by-N array of positions to start at (overrides \"-np\")",
974
             NULL, NULL, nrrdHestNrrd);
975
  hestOptAdd(&hopt, "step", "step", airTypeDouble, 1, 1,
976
             &stepInitial, "1",
977
             "initial step size for gradient descent");
978
  hestOptAdd(&hopt, "csm", "step", airTypeDouble, 1, 1,
979
             &constraintStepMin, "0.0001",
980
             "convergence criterion for constraint satisfaction");
981
  hestOptAdd(&hopt, "snap", "# iters", airTypeUInt, 1, 1,
982
             &snap, "0",
983
             "if non-zero, # iters between saved snapshots");
984
  hestOptAdd(&hopt, "maxi", "# iters", airTypeUInt, 1, 1,
985
             &iterMax, "0",
986
             "if non-zero, max # iterations to run whole system");
987
  hestOptAdd(&hopt, "stim", "# iters", airTypeUInt, 1, 1,
988
             &stuckIterMax, "5",
989
             "if non-zero, max # iterations to allow a particle "
990
             " to be stuck before nixing");
991
  hestOptAdd(&hopt, "maxci", "# iters", airTypeUInt, 1, 1,
992
             &constraintIterMax, "15",
993
             "if non-zero, max # iterations for contraint enforcement");
994
  hestOptAdd(&hopt, "irad", "scale", airTypeDouble, 1, 1,
995
             &radiusSpace, "1",
996
             "particle radius in spatial domain");
997
  hestOptAdd(&hopt, "srad", "scale", airTypeDouble, 1, 1,
998
             &radiusScale, "1",
999
             "particle radius in scale domain");
1000
  hestOptAdd(&hopt, "bws", "bin width", airTypeDouble, 1, 1,
1001
             &binWidthSpace, "1.001",
1002
             "spatial bin width as multiple of spatial radius");
1003
  hestOptAdd(&hopt, "alpha", "alpha", airTypeDouble, 1, 1,
1004
             &alpha, "0.5",
1005
             "blend between particle-image (alpha=0) and "
1006
             "inter-particle (alpha=1) energies");
1007
  hestOptAdd(&hopt, "beta", "beta", airTypeDouble, 1, 1,
1008
             &beta, "1.0",
1009
             "when using Phi2 energy, blend between pure "
1010
             "space repulsion (beta=0) and "
1011
             "scale attraction (beta=1)");
1012
  hestOptAdd(&hopt, "gamma", "gamma", airTypeDouble, 1, 1,
1013
             &_gamma, "1.0",
1014
             "scaling factor on energy from strength");
1015
  hestOptAdd(&hopt, "wall", "k", airTypeDouble, 1, 1,
1016
             &wall, "0.0",
1017
             "spring constant on walls");
1018
  hestOptAdd(&hopt, "eip", "k", airTypeDouble, 1, 1,
1019
             &energyIncreasePermit, "0.0",
1020
             "amount by which its okay for *per-particle* energy to increase "
1021
             "during gradient descent process");
1022
  hestOptAdd(&hopt, "ess", "scl", airTypeDouble, 1, 1,
1023
             &backStepScale, "0.5",
1024
             "when energy goes up instead of down, scale step "
1025
             "size by this");
1026
  hestOptAdd(&hopt, "oss", "scl", airTypeDouble, 1, 1,
1027
             &opporStepScale, "1.0",
1028
             "opportunistic scaling (hopefully up, >1) of step size "
1029
             "on every iteration");
1030
  hestOptAdd(&hopt, "edmin", "frac", airTypeDouble, 1, 1,
1031
             &energyDecreaseMin, "0.0001",
1032
             "convergence threshold: stop when fractional improvement "
1033
             "(decrease) in energy dips below this");
1034
  hestOptAdd(&hopt, "edpcmin", "frac", airTypeDouble, 1, 1,
1035
             &energyDecreasePopCntlMin, "0.01",
1036
             "population control is triggered when energy improvement "
1037
             "goes below this threshold");
1038
  hestOptAdd(&hopt, "fnnm", "frac", airTypeDouble, 1, 1,
1039
             &fracNeighNixedMax, "0.25",
1040
             "don't nix if this fraction (or more) of neighbors "
1041
             "have been nixed");
1042
  hestOptAdd(&hopt, "pcp", "period", airTypeUInt, 1, 1,
1043
             &popCntlPeriod, "20",
1044
             "# iters to wait between attempts at population control");
1045
  hestOptAdd(&hopt, "iad", "# iters", airTypeUInt, 1, 1,
1046
             &addDescent, "10",
1047
             "# iters to run descent on tentative new points during PC");
1048
  hestOptAdd(&hopt, "icb", "# iters", airTypeUInt, 1, 1,
1049
             &iterCallback, "1",
1050
             "periodicity of calling rendering callback");
1051
1052
  hestOptAdd(&hopt, "ac3c", "ac3c", airTypeBool, 1, 1,
1053
             &allowCodimension3Constraints, "false",
1054
             "allow codimensions 3 constraints");
1055
  hestOptAdd(&hopt, "sit", "sit", airTypeBool, 1, 1, &scaleIsTau, "false",
1056
             "scale is tau");
1057
  hestOptAdd(&hopt, "rng", "seed", airTypeUInt, 1, 1,
1058
             &rngSeed, "42",
1059
             "base seed value for RNGs");
1060
  hestOptAdd(&hopt, "pbm", "mod", airTypeUInt, 1, 1,
1061
             &progressBinMod, "50",
1062
             "progress bin mod");
1063
  hestOptAdd(&hopt, "eiphl", "hl", airTypeUInt, 1, 1, &eipHalfLife, "0",
1064
             "half-life of energyIncreasePermute (\"-eip\")");
1065
  hestOptAdd(&hopt, "nt", "# threads", airTypeInt, 1, 1,
1066
             &threadNum, "1",
1067
             (airThreadCapable
1068
              ? "number of threads hoover should use"
1069
              : "if threads where enabled in this Teem build, this is how "
1070
              "you would control the number of threads to use"));
1071
  hestOptAdd(&hopt, "nprob", "prob", airTypeDouble, 1, 1,
1072
             &neighborTrueProb, "1.0",
1073
             "do full neighbor discovery with this probability");
1074
  hestOptAdd(&hopt, "pprob", "prob", airTypeDouble, 1, 1,
1075
             &probeProb, "1.0",
1076
             "probe local image values with this probability");
1077
1078
  hestOptAdd(&hopt, "addlog", "fname", airTypeString, 1, 1, &addLogS, "",
1079
             "name of file in which to log all particle additions");
1080
  hestOptAdd(&hopt, "o", "nout", airTypeString, 1, 1, &outS, "-",
1081
             "filename for saving computed positions");
1082
  hestOptAdd(&hopt, "eob", "base", airTypeString, 1, 1, &extraOutBaseS, "",
1083
             "save extra info (besides position), and use this string as "
1084
             "the base of the filenames.  Not using this means the extra "
1085
             "info is not saved.");
1086
1087
  hestParseOrDie(hopt, argc-1, argv+1, hparm,
1088
                 me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE);
1089
  airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
1090
  airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
1091
1092
#ifdef DEFT
1093
  if (const char *envS = getenv("DEFT_HOME")) {
1094
    strcpy(Deft::homeDir, envS);
1095
    strcat(Deft::homeDir, "/");
1096
  } else {
1097
    fprintf(stderr, "%s: WARNING: \"DEFT_HOME\" environment variable "
1098
            "not set; assuming \".\"\n", me);
1099
    strcpy(Deft::homeDir, "./");
1100
  }
1101
#endif
1102
  /*
1103
  airEnumPrint(stderr, gageScl);
1104
  exit(0);
1105
  */
1106
  if (airStrlen(addLogS)) {
1107
    if (!(addLog = airFopen(addLogS, stdout, "w"))) {
1108
      fprintf(stderr, "%s: couldn't open %s for writing", me, addLogS);
1109
      airMopError(mop); return 1;
1110
    }
1111
    airMopAdd(mop, addLog, (airMopper)airFclose, airMopAlways);
1112
  } else {
1113
    addLog = NULL;
1114
  }
1115
1116
  pctx = pullContextNew();
1117
  airMopAdd(mop, pctx, (airMopper)pullContextNix, airMopAlways);
1118
  if (pullVerboseSet(pctx, verbose)
1119
      || pullFlagSet(pctx, pullFlagZeroZ, zeroZ)
1120
      || pullFlagSet(pctx, pullFlagEnergyFromStrength, energyFromStrength)
1121
      || pullFlagSet(pctx, pullFlagNixAtVolumeEdgeSpace, nixAtVolumeEdgeSpace)
1122
      || pullFlagSet(pctx, pullFlagConstraintBeforeSeedThresh,
1123
                     constraintBeforeSeedThresh)
1124
      || pullFlagSet(pctx, pullFlagPopCntlEnoughTest, popCntlEnoughTest)
1125
      || pullFlagSet(pctx, pullFlagConvergenceIgnoresPopCntl,
1126
                     convergenceIgnoresPopCntl)
1127
      || pullFlagSet(pctx, pullFlagBinSingle, binSingle)
1128
      || pullFlagSet(pctx, pullFlagNoAdd, noAdd)
1129
      || pullFlagSet(pctx, pullFlagPermuteOnRebin, permuteOnRebin)
1130
      || pullFlagSet(pctx, pullFlagNoPopCntlWithZeroAlpha,
1131
                     noPopCntlWithZeroAlpha)
1132
      || pullFlagSet(pctx, pullFlagUseBetaForGammaLearn,
1133
                     useBetaForGammaLearn)
1134
      || pullFlagSet(pctx, pullFlagRestrictiveAddToBins,
1135
                     restrictiveAddToBins)
1136
      || pullFlagSet(pctx, pullFlagAllowCodimension3Constraints,
1137
                     allowCodimension3Constraints)
1138
      || pullFlagSet(pctx, pullFlagScaleIsTau, scaleIsTau)
1139
      || pullInitUnequalShapesAllowSet(pctx, unequalShapesAllow)
1140
      || pullIterParmSet(pctx, pullIterParmSnap, snap)
1141
      || pullIterParmSet(pctx, pullIterParmMax, iterMax)
1142
      || pullIterParmSet(pctx, pullIterParmStuckMax, stuckIterMax)
1143
      || pullIterParmSet(pctx, pullIterParmConstraintMax, constraintIterMax)
1144
      || pullIterParmSet(pctx, pullIterParmPopCntlPeriod, popCntlPeriod)
1145
      || pullIterParmSet(pctx, pullIterParmAddDescent, addDescent)
1146
      || pullIterParmSet(pctx, pullIterParmCallback, iterCallback)
1147
      || pullIterParmSet(pctx, pullIterParmEnergyIncreasePermitHalfLife,
1148
                         eipHalfLife)
1149
      || pullSysParmSet(pctx, pullSysParmStepInitial, stepInitial)
1150
      || pullSysParmSet(pctx, pullSysParmConstraintStepMin, constraintStepMin)
1151
      || pullSysParmSet(pctx, pullSysParmRadiusSpace, radiusSpace)
1152
      || pullSysParmSet(pctx, pullSysParmRadiusScale, radiusScale)
1153
      || pullSysParmSet(pctx, pullSysParmBinWidthSpace, binWidthSpace)
1154
      || pullSysParmSet(pctx, pullSysParmAlpha, alpha)
1155
      || pullSysParmSet(pctx, pullSysParmBeta, beta)
1156
      || pullSysParmSet(pctx, pullSysParmGamma, _gamma)
1157
      || pullSysParmSet(pctx, pullSysParmWall, wall)
1158
      || pullSysParmSet(pctx, pullSysParmEnergyIncreasePermit,
1159
                        energyIncreasePermit)
1160
      || pullSysParmSet(pctx, pullSysParmEnergyDecreaseMin,
1161
                        energyDecreaseMin)
1162
      || pullSysParmSet(pctx, pullSysParmFracNeighNixedMax,
1163
                        fracNeighNixedMax)
1164
      || pullSysParmSet(pctx, pullSysParmEnergyDecreasePopCntlMin,
1165
                        energyDecreasePopCntlMin)
1166
      || pullSysParmSet(pctx, pullSysParmBackStepScale, backStepScale)
1167
      || pullSysParmSet(pctx, pullSysParmOpporStepScale, opporStepScale)
1168
      || pullSysParmSet(pctx, pullSysParmNeighborTrueProb,
1169
                        neighborTrueProb)
1170
      || pullSysParmSet(pctx, pullSysParmProbeProb, probeProb)
1171
      || pullRngSeedSet(pctx, rngSeed)
1172
      || pullProgressBinModSet(pctx, progressBinMod)
1173
      || pullThreadNumSet(pctx, threadNum)
1174
      || pullInterEnergySet(pctx, interType, enspR, enspS, enspWin)
1175
      || pullInitLiveThreshUseSet(pctx, liveThresholdOnInit)
1176
      || pullLogAddSet(pctx, addLog)) {
1177
    airMopAdd(mop, err = biffGetDone(PULL), airFree, airMopAlways);
1178
    fprintf(stderr, "%s: trouble with flags:\n%s", me, err);
1179
    airMopError(mop); return 1;
1180
  }
1181
1182
  if (nPosIn) {
1183
    E = pullInitGivenPosSet(pctx, nPosIn);
1184
  } else if (pointPerVoxel) {
1185
    E = pullInitPointPerVoxelSet(pctx, pointPerVoxel,
1186
                                 ppvZRange[0], ppvZRange[1],
1187
                                 samplesAlongScaleNum, jitter);
1188
  } else if (useHalton) {
1189
    E = pullInitHaltonSet(pctx, pointNumInitial, 0);
1190
  } else {
1191
    E = pullInitRandomSet(pctx, pointNumInitial);
1192
  }
1193
  if (E) {
1194
    airMopAdd(mop, err = biffGetDone(PULL), airFree, airMopAlways);
1195
    fprintf(stderr, "%s: trouble with flags:\n%s", me, err);
1196
    airMopError(mop); return 1;
1197
  }
1198
  if (meetPullVolStackBlurParmFinishMulti(vspec, vspecNum,
1199
                                          &kssFinished, &bspFinished,
1200
                                          kSSblur, bspec)
1201
      || meetPullVolLoadMulti(vspec, vspecNum, cachePathSS, verbose)
1202
      || meetPullVolAddMulti(pctx, vspec, vspecNum,
1203
                             k00, k11, k22, kSSrecon)
1204
      || meetPullInfoAddMulti(pctx, idef, idefNum)) {
1205
    airMopAdd(mop, err = biffGetDone(MEET), airFree, airMopAlways);
1206
    fprintf(stderr, "%s: trouble with volumes or infos:\n%s", me, err);
1207
    airMopError(mop); return 1;
1208
  }
1209
  if (!kssFinished && hestSourceUser == hopt[kssOpi].source) {
1210
    fprintf(stderr, "\n\n%s: WARNING! Used the -%s flag, but the "
1211
            "meetPullVol specified blurring kernels\n\n\n", me,
1212
            hopt[kssOpi].flag);
1213
  }
1214
  if (!bspFinished && hestSourceUser == hopt[bspOpi].source) {
1215
    fprintf(stderr, "\n\n%s: WARNING! Used the -%s flag, but the "
1216
            "meetPullVol specified boundary specs\n\n\n", me,
1217
            hopt[bspOpi].flag);
1218
  }
1219
  if (pullStart(pctx)) {
1220
    airMopAdd(mop, err = biffGetDone(PULL), airFree, airMopAlways);
1221
    fprintf(stderr, "%s: trouble starting system:\n%s", me, err);
1222
    airMopError(mop); return 1;
1223
  }
1224
1225
  /* -------------------------------------------------- */
1226
1227
  /* not sure when this table was created, don't have heart to nix it
1228
   *
1229
   *                 hght scl   tang1    tang2   mode scl  strength
1230
   *  ridge surface:    -1      evec2      -        -       -eval2
1231
   *     ridge line:    -1      evec2    evec1      -       -eval1
1232
   *     all ridges:    -1      evec2    evec1     +1        ??
1233
   * valley surface:    +1      evec0      -        -        eval0
1234
   *    valley line:    +1      evec0    evec1      -        eval1
1235
   *      all lines:    +1      evec0    evec1     -1
1236
   */
1237
1238
#ifdef DEFT
1239
  ssrange[0] = FLT_MAX;
1240
  ssrange[1] = -FLT_MAX;
1241
  for (vsi=0; vsi<vspecNum; vsi++) {
1242
    meetPullVol *vol;
1243
    vol = vspec[vsi];
1244
    if (vol->numSS) {
1245
      ssrange[0] = AIR_MIN(ssrange[0], vol->rangeSS[0]);
1246
      ssrange[1] = AIR_MAX(ssrange[1], vol->rangeSS[1]);
1247
    }
1248
  }
1249
  if (pctx->flag.scaleIsTau) {
1250
    ssrange[0] = gageTauOfSig(ssrange[0]);
1251
    ssrange[1] = gageTauOfSig(ssrange[1]);
1252
  }
1253
  /* -------------------------------------------------- */
1254
  /* initialize bag and its UI */
1255
  bag.pctx = pctx;
1256
  bag.nPosOut = nrrdNew();
1257
  bag.nTenOut = nrrdNew();
1258
  bag.nFrcOut = nrrdNew();
1259
  bag.nposOld = nrrdNew();
1260
  bag.nposNew = nrrdNew();
1261
  bag.nten = nrrdNew();
1262
  bag.ntmp = nrrdNew();
1263
  bag.nenr = nrrdNew();
1264
  bag.nscl = nrrdNew();
1265
  bag.nidcc = nrrdNew();
1266
  bag.ncovar = nrrdNew();
1267
#if PULL_TANCOVAR
1268
  bag.ntcovar = nrrdNew();
1269
#endif
1270
  bag.nstab = nrrdNew();
1271
  bag.nintern = nrrdNew();
1272
  if (pctx->ispec[pullInfoStrength]) {
1273
    printf("!%s: trouble creating strength nrrd\n", me);
1274
    bag.nstrn = nrrdNew();
1275
  } else {
1276
    bag.nstrn = NULL;
1277
  }
1278
  if (pctx->ispec[pullInfoQuality]) {
1279
    printf("!%s: trouble creating quality nrrd\n", me);
1280
    bag.nqual = nrrdNew();
1281
  } else {
1282
    bag.nqual = NULL;
1283
  }
1284
  bag.nstuck = nrrdNew();
1285
  bag.nfrcOld = nrrdNew();
1286
  bag.nfrcNew = nrrdNew();
1287
  bag.nrgb = nrrdNew();
1288
  bag.nccrgb = nrrdNew();
1289
  bag.ncmapOut = nrrdNew();
1290
  bag.nblur = nrrdNew();
1291
  bag.norig = vspec[0]->nin;
1292
  nrrdCopy(bag.nblur, bag.norig);
1293
  bag.rsmc = nrrdResampleContextNew();
1294
  E = 0;
1295
  if (!E) E |= nrrdResampleDefaultCenterSet(bag.rsmc, nrrdDefaultCenter);
1296
  if (!E) E |= nrrdResampleInputSet(bag.rsmc, bag.norig);
1297
  for (unsigned int axi=0; axi<3; axi++) {
1298
    if (!E) E |= nrrdResampleSamplesSet(bag.rsmc, axi,
1299
                                        bag.norig->axis[axi].size);
1300
    if (!E) E |= nrrdResampleRangeFullSet(bag.rsmc, axi);
1301
  }
1302
  if (!E) E |= nrrdResampleBoundarySet(bag.rsmc, nrrdBoundaryBleed);
1303
  if (!E) E |= nrrdResampleTypeOutSet(bag.rsmc, nrrdTypeDefault);
1304
  if (!E) E |= nrrdResampleRenormalizeSet(bag.rsmc, AIR_TRUE);
1305
  if (E) {
1306
    airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
1307
    fprintf(stderr, "%s: trouble setting up resampler:\n%s", me, err);
1308
    airMopError(mop); return 1;
1309
  }
1310
  /* bag.ncval is just a pointer to other nrrds */
1311
  bag.cvalRange = nrrdRangeNew(AIR_NAN, AIR_NAN);
1312
  ELL_3V_COPY(bag.scaleVec, scaleVec);
1313
  bag.glyphScaleRad = glyphScaleRad;
1314
1315
  bag.scene = new Deft::Scene();
1316
  bag.scene->bgColor(bg[0], bg[1], bg[2]);
1317
1318
  int winy = 10;
1319
  int incy;
1320
  fltk::Window *win = new fltk::Window(400, 600, "pull UI");
1321
  win->begin();
1322
  fltk::Button *stepButton = new fltk::Button(10, winy, 50, incy=20, "step");
1323
  stepButton->callback((fltk::Callback*)step_cb, &bag);
1324
1325
  bag.verbose = new fltk::ValueInput(100, winy, 30, 20, "verb");
1326
  bag.verbose->value(pctx->verbose);
1327
  bag.verbose->callback((fltk::Callback*)verbose_cb, &bag);
1328
1329
  bag.iters = new fltk::IntInput(200, winy, 50, 20, "# iters");
1330
  bag.iters->value(1);
1331
1332
  if (ssrange[1] > ssrange[0]) {
1333
    fltk::Button *gamButton = new fltk::Button(260, winy,
1334
                                               50, 20, "gamma");
1335
    gamButton->callback((fltk::Callback*)gammaSet_cb, &bag);
1336
  }
1337
  fltk::Button *ccButton = new fltk::Button(360, winy,
1338
                                            30, 20, "CC");
1339
  ccButton->callback((fltk::Callback*)cc_cb, &bag);
1340
1341
  winy += incy + 5;
1342
  fltk::Button *saveButton = new fltk::Button(10, winy, 50, 20, "save");
1343
  saveButton->callback((fltk::Callback*)save_cb, &bag);
1344
1345
  bag.scaleVecInput[0] = new fltk::FloatInput(120, winy, 35, 20, "scaleVec");
1346
  bag.scaleVecInput[1] = new fltk::FloatInput(160, winy, 35, 20, "");
1347
  bag.scaleVecInput[2] = new fltk::FloatInput(200, winy, 35, 20, "");
1348
  bag.scaleVecInput[0]->value(scaleVec[0]);
1349
  bag.scaleVecInput[1]->value(scaleVec[1]);
1350
  bag.scaleVecInput[2]->value(scaleVec[2]);
1351
  bag.scaleVecInput[0]->callback((fltk::Callback*)scaleGlyph_cb, &bag);
1352
  bag.scaleVecInput[1]->callback((fltk::Callback*)scaleGlyph_cb, &bag);
1353
  bag.scaleVecInput[2]->callback((fltk::Callback*)scaleGlyph_cb, &bag);
1354
1355
  bag.glyphScaleRadInput = new fltk::ValueInput(300, winy, 45, 20, "gssr");
1356
  bag.glyphScaleRadInput->range(0.0, 100.0);
1357
  bag.glyphScaleRadInput->step(0.1);
1358
  bag.glyphScaleRadInput->linesize(0.1);
1359
  bag.glyphScaleRadInput->value(glyphScaleRad);
1360
  bag.glyphScaleRadInput->callback((fltk::Callback*)scaleGlyph_cb, &bag);
1361
1362
  winy += incy;
1363
  bag.alpha = new Deft::Slider(0, winy, win->w(), incy=55, "alpha");
1364
  bag.alpha->align(fltk::ALIGN_LEFT);
1365
  bag.alpha->range(0, 1);
1366
  bag.alpha->value(pctx->sysParm.alpha);
1367
  bag.alpha->fastUpdate(1);
1368
  bag.alpha->callback((fltk::Callback*)alpha_cb, &bag);
1369
1370
  if (pullInterTypeAdditive == pctx->interType) {
1371
    winy += incy;
1372
    bag.beta = new Deft::Slider(0, winy, win->w(), incy=55, "beta");
1373
    bag.beta->align(fltk::ALIGN_LEFT);
1374
    bag.beta->range(0, 1);
1375
    bag.beta->value(pctx->sysParm.beta);
1376
    bag.beta->fastUpdate(1);
1377
    bag.beta->callback((fltk::Callback*)beta_cb, &bag);
1378
  }
1379
1380
  if (pullEnergyCubicWell == pctx->energySpecR->energy
1381
      || pullEnergyBetterCubicWell == pctx->energySpecR->energy) {
1382
    winy += incy;
1383
    bag.cwell = new Deft::Slider(0, winy, win->w(), incy=55, "well depth");
1384
    bag.cwell->align(fltk::ALIGN_LEFT);
1385
    bag.cwell->range(-0.04, 0);
1386
    bag.cwell->value(bag.pctx->energySpecR->parm[1]);
1387
    bag.cwell->fastUpdate(1);
1388
    bag.cwell->callback((fltk::Callback*)cwell_cb, &bag);
1389
    /* remember eip as fraction of well depth */
1390
    bag.energyIncreasePermitFrac =
1391
      energyIncreasePermit/bag.pctx->energySpecR->parm[1];
1392
  } else {
1393
    bag.energyIncreasePermitFrac = AIR_NAN;
1394
  }
1395
1396
  winy += incy;
1397
  bag.gamma = new Deft::Slider(0, winy, win->w(), incy=55, "gamma");
1398
  bag.gamma->align(fltk::ALIGN_LEFT);
1399
  bag.gamma->range(0, 2*pctx->sysParm.gamma);
1400
  bag.gamma->value(pctx->sysParm.gamma);
1401
  bag.gamma->fastUpdate(1);
1402
  bag.gamma->callback((fltk::Callback*)gamma_cb, &bag);
1403
1404
  winy += incy;
1405
  bag.ccSelect = new Deft::Slider(0, winy, win->w(), incy=55, "CC Select");
1406
  bag.ccSelect->align(fltk::ALIGN_LEFT);
1407
  bag.ccSelect->range(0, 0);
1408
  bag.ccSelect->step(1);
1409
  bag.ccSelect->value(0);
1410
  bag.ccSelect->fastUpdate(1);
1411
  bag.ccSelect->callback((fltk::Callback*)ccSelect_cb, &bag);
1412
1413
  bag.ccSingle = new fltk::CheckButton(130, winy+4, 50, 20, "Single");
1414
  bag.ccSingle->value(0);
1415
  bag.ccSingle->callback((fltk::Callback*)ccSelect_cb, &bag);
1416
1417
  winy += incy;
1418
  bag.rho = new Deft::Slider(0, winy, win->w(), incy=55, "rho");
1419
  bag.rho->align(fltk::ALIGN_LEFT);
1420
  bag.rho->range(0, 1.0);
1421
  bag.rho->value(0.5);
1422
  bag.rho->fastUpdate(1);
1423
  bag.rho->callback((fltk::Callback*)cc_cb, &bag);
1424
1425
  if (ssrange[1] > ssrange[0]) {
1426
    winy += incy;
1427
    bag.sclMean = new Deft::Slider(0, winy, win->w(), incy=55, "scale mean");
1428
    bag.sclMean->align(fltk::ALIGN_LEFT);
1429
    bag.sclMean->range(ssrange[0], ssrange[1]);
1430
    bag.sclMean->value((ssrange[0] + ssrange[1])/2);
1431
    bag.sclMean->fastUpdate(1);
1432
    bag.sclMean->callback((fltk::Callback*)scale_cb, &bag);
1433
1434
    fltk::Button *reblurButton = new fltk::Button(130, winy+4, 50, 20, "reblur");
1435
    reblurButton->callback((fltk::Callback*)reblur_cb, &bag);
1436
1437
    winy += incy;
1438
    bag.sclWind = new Deft::Slider(0, winy, win->w(), incy=55, "scale window");
1439
    bag.sclWind->align(fltk::ALIGN_LEFT);
1440
    bag.sclWind->range(0, ssrange[1] - ssrange[0]);
1441
    bag.sclWind->value(ssrange[1] - ssrange[0]);
1442
    bag.sclWind->fastUpdate(1);
1443
    bag.sclWind->callback((fltk::Callback*)scale_cb, &bag);
1444
    scale_cb(NULL, &bag);
1445
  } else {
1446
    bag.sclMin = bag.sclMax = 0;
1447
  }
1448
1449
  if (pctx->ispec[pullInfoStrength]) {
1450
    winy += incy;
1451
    bag.strength = new Deft::Slider(0, winy, win->w(), incy=55, "strength");
1452
    bag.strength->align(fltk::ALIGN_LEFT);
1453
    bag.strength->range(0, 1);
1454
    bag.strength->callback((fltk::Callback*)strength_cb, &bag);
1455
    bag.strength->fastUpdate(1);
1456
    bag.strength->value(0);
1457
  } else {
1458
    bag.strnMin = 0;
1459
  }
1460
  if (pctx->ispec[pullInfoQuality]) {
1461
    winy += incy;
1462
    bag.quality = new Deft::Slider(0, winy, win->w(), incy=55, "quality");
1463
    bag.quality->align(fltk::ALIGN_LEFT);
1464
    bag.quality->range(-0.1, 1);
1465
    bag.quality->callback((fltk::Callback*)quality_cb, &bag);
1466
    bag.quality->fastUpdate(1);
1467
    bag.quality->value(-0.1);
1468
  } else {
1469
    bag.qualMin = 0;
1470
  }
1471
1472
  /*
1473
  winy += incy;
1474
  bag.height = new Deft::Slider(0, winy, win->w(), incy=55, "height");
1475
  bag.height->align(fltk::ALIGN_LEFT);
1476
  bag.height->range(0, 10);
1477
  bag.height->value(10);
1478
  bag.height->fastUpdate(1);
1479
  bag.height->callback((fltk::Callback*)strength_cb, &bag);
1480
  */
1481
1482
  win->end();
1483
  win->show();
1484
1485
  /* -------------------------------------------------- */
1486
  bag.viewer = new Deft::Viewer(bag.scene, imgSize[0], imgSize[1]);
1487
  if (camkeep) {
1488
    bag.viewer->camera(fr[0], fr[1], fr[2],
1489
                       at[0], at[1], at[2],
1490
                       up[0], up[1], up[2],
1491
                       fovy, neer, faar);
1492
  }
1493
  bag.viewer->resizable(bag.viewer);
1494
  bag.viewer->end();
1495
  const char *fakeArgv[2] = {"Deft_pull", NULL};
1496
  bag.viewer->show(1, (char**)fakeArgv);
1497
  if (ortho) {
1498
    /* total hack */
1499
    bag.viewer->keyboard('p', 0, 0);
1500
  }
1501
  /* bag.viewer->helpPrint(); */
1502
  Deft::ViewerUI *viewerUI = new Deft::ViewerUI(bag.viewer);
1503
  viewerUI->show();
1504
  fltk::flush();
1505
1506
  /* -------------------------------------------------- */
1507
  if (gageKindScl == vspec[0]->kind) {
1508
    bag.contour = new Deft::Contour();
1509
    bag.contour->volumeSet(bag.nblur);
1510
    bag.contour->twoSided(true);
1511
    bag.scene->objectAdd(bag.contour);
1512
  } else {
1513
    bag.contour = NULL;
1514
  }
1515
1516
  /* -------------------------------------------------- */
1517
  Deft::TensorGlyph *glyph = new Deft::TensorGlyph();
1518
  glyph->dynamic(true);
1519
  glyph->twoSided(true);
1520
  glyph->anisoType(aniso);
1521
  glyph->anisoThreshMin(anisoThreshMin);
1522
  glyph->anisoThresh(anisoThresh);
1523
  glyph->glyphType(glyphType);
1524
  glyph->superquadSharpness(sqdSharp);
1525
  glyph->glyphResolution(glyphFacetRes);
1526
  if (tenGlyphTypeBetterquad) {
1527
    glyph->barycentricResolution(baryRes);
1528
  } else {
1529
    glyph->barycentricResolution(20);
1530
  }
1531
  glyph->glyphScale(glyphScale);
1532
  glyph->glyphHaloWidth(glyphHaloWidth);
1533
  glyph->haloTraceBound(haloTraceBound);
1534
  glyph->glyphNormPow(glyphNormPow);
1535
  glyph->glyphEvalPow(glyphEvalPow);
1536
  glyph->rgbEvecParmSet(tenAniso_Cl2, 0, 0.7, 1.0, 2.3, 1.0);
1537
  glyph->rgbEvecParmSet(tenAniso_Cl2, 0, 0, 1.0, 1, 0.0);
1538
  glyph->maskThresh(0.5);
1539
  /*
1540
  void rgbParmSet(int aniso, unsigned int evec,
1541
                  double maxSat, double isoGray,
1542
                  double gamma, double modulate);
1543
  */
1544
1545
  bag.scene->objectAdd(glyph);
1546
  bag.glyph = glyph;
1547
1548
  Deft::TensorGlyphUI *glyphUI = new Deft::TensorGlyphUI(bag.glyph,
1549
                                                         bag.viewer);
1550
  glyphUI->show();
1551
1552
  /* -------------------------------------------------- */
1553
  /*
1554
  Deft::TensorGlyph *hedge = new Deft::TensorGlyph();
1555
  hedge->dynamic(true);
1556
  hedge->anisoType(aniso);
1557
  hedge->anisoThreshMin(anisoThreshMin);
1558
  hedge->anisoThresh(anisoThresh);
1559
  hedge->glyphType(glyphType);
1560
  hedge->superquadSharpness(sqdSharp);
1561
  hedge->glyphResolution(glyphFacetRes);
1562
  hedge->glyphScale(glyphScale/4);
1563
  hedge->rgbParmSet(tenAniso_Cl2, 0, 0.7, 1.0, 2.3, 1.0);
1564
  hedge->rgbParmSet(tenAniso_Cl2, 0, 0, 1.0, 1, 0.0);
1565
  hedge->maskThresh(0.0);
1566
  bag.hedge = hedge;
1567
  bag.scene->objectAdd(hedge);
1568
  Deft::TensorGlyphUI *hedgeUI = new Deft::TensorGlyphUI(bag.hedge,
1569
                                                         bag.viewer);
1570
  hedgeUI->show();
1571
  */
1572
1573
  /* -------------------------------------------------- */
1574
  Deft::Volume *vol = new Deft::Volume(vspec[0]->kind, bag.nblur);
1575
  fprintf(stderr, "!%s: vol = %p *********************\n", me, vol);
1576
  Deft::TriPlane *triplane = new Deft::TriPlane(vol);
1577
  /*
1578
  ** HEY: you can eventually segfault if this isn't set here
1579
  ** shouldn't doing so be optional?
1580
  */
1581
  if (gageKindScl == vspec[0]->kind) {
1582
    triplane->alphaMaskQuantity(Deft::alphaMaskSclQuantityValue);
1583
    triplane->alphaMaskThreshold(0);
1584
    triplane->colorQuantity(Deft::colorSclQuantityValue);
1585
  } else if (tenGageKind == vspec[0]->kind) {
1586
    triplane->alphaMaskQuantity(Deft::alphaMaskTenQuantityConfidence);
1587
    triplane->alphaMaskThreshold(0);
1588
    triplane->colorQuantity(Deft::colorTenQuantityRgbLinear);
1589
  }
1590
1591
  NrrdKernelSpec *ksp = nrrdKernelSpecNew();
1592
  double kparm[10];
1593
1594
  kparm[0] = 1.0;
1595
  nrrdKernelSpecSet(ksp, nrrdKernelTent, kparm);
1596
  triplane->kernel(gageKernel00, ksp);
1597
  ELL_3V_SET(kparm, 1, 1, 0);
1598
  nrrdKernelSpecSet(ksp, nrrdKernelBCCubicD, kparm);
1599
  triplane->kernel(gageKernel11, ksp);
1600
  nrrdKernelSpecSet(ksp, nrrdKernelBCCubicDD, kparm);
1601
  triplane->kernel(gageKernel22, ksp);
1602
  triplane->visible(false);
1603
1604
  bag.scene->groupAdd(triplane);
1605
1606
  Deft::TriPlaneUI *planeUI = new Deft::TriPlaneUI(triplane, bag.viewer);
1607
  planeUI->show();
1608
1609
  /* -------------------------------------------------- */
1610
1611
  if (gageKindScl == vspec[0]->kind) {
1612
    fltk::Window *window = new fltk::Window(400, 80, "isovalue");
1613
    window->begin();
1614
    window->resizable(window);
1615
1616
    winy = 0;
1617
    bag.isoval = new Deft::Slider(0, winy, window->w(), incy=55, "isovalue");
1618
    winy += incy;
1619
    bag.isoval->align(fltk::ALIGN_LEFT);
1620
    bag.isoval->range(bag.contour->minimum(), bag.contour->maximum());
1621
    bag.isoval->value(bag.contour->maximum());
1622
    bag.contour->wireframe(false);
1623
    bag.isoval->fastUpdate(0);
1624
    bag.isoval->callback((fltk::Callback*)isovalue_cb, &bag);
1625
    /* bag.isoval->value(1.0); */
1626
    window->end();
1627
    window->show(argc,(char**)argv);
1628
  }
1629
1630
  if (pullPhistEnabled) {
1631
    bag.phistLine = limnPolyDataNew();
1632
    /* bag.phistTube = limnPolyDataNew(); */
1633
    bag.phistSurf = new Deft::PolyData(bag.phistLine, false);
1634
    bag.phistSurf->wireframe(false);
1635
    bag.phistSurf->normalsUse(true);
1636
    bag.phistSurf->colorUse(true);
1637
    bag.phistSurf->visible(true);
1638
    bag.phistSurf->changed();
1639
    bag.scene->objectAdd(bag.phistSurf);
1640
  }
1641
1642
  fltk::flush();
1643
  fltk::redraw();
1644
1645
  /* this is to get "output" prior to any iterations,
1646
     but after seeding and initialization */
1647
  outputGet(&bag);
1648
  outputShow(&bag);
1649
1650
  if (!camkeep) {
1651
    bag.viewer->cameraReset();
1652
  }
1653
1654
  if (0) {
1655
    Nrrd *nplot;
1656
    nplot = nrrdNew();
1657
    if (pullEnergyPlot(pctx, nplot, 1, 0, 0, 601)) {
1658
      airMopAdd(mop, err = biffGetDone(PULL), airFree, airMopAlways);
1659
      fprintf(stderr, "%s: trouble plotting:\n%s", me, err);
1660
      airMopError(mop); return 1;
1661
    }
1662
    nrrdSave("eplot.nrrd", nplot, NULL);
1663
    nplot = nrrdNuke(nplot);
1664
#if PULL_HINTER
1665
    pctx->nhinter = nrrdNew();
1666
#endif
1667
  }
1668
1669
  if (fog) {
1670
    bag.viewer->fog(true);
1671
    bag.viewer->redraw();
1672
    ret = fltk::wait();
1673
  }
1674
  if (saveAndQuit) {
1675
    bag.scene->draw();
1676
    bag.viewer->screenDump();
1677
  } else {
1678
1679
    /* set up callbacks in pull */
1680
    pullCallbackSet(pctx, iter_cb, &bag);
1681
1682
    /* this returns when the user quits */
1683
    ret = fltk::run();
1684
  }
1685
1686
#else
1687
  /* else not running as Deft, but as command-line */
1688
  if (!E) E |= pullRun(pctx);
1689
  if (E) {
1690
    airMopAdd(mop, err = biffGetDone(PULL), airFree, airMopAlways);
1691
    fprintf(stderr, "%s: trouble 3:\n%s", me, err);
1692
    airMopError(mop); return 1;
1693
  }
1694
  if (pullOutputGet(nPosOut, NULL, NULL, NULL, 0.0, pctx)) {
1695
    airMopAdd(mop, err = biffGetDone(PULL), airFree, airMopAlways);
1696
    fprintf(stderr, "%s: trouble 3.1:\n%s", me, err);
1697
    airMopError(mop); return 1;
1698
  }
1699
  nrrdSave(outS, nPosOut, NULL);
1700
  if (airStrlen(extraOutBaseS)) {
1701
    Nrrd *nstrn, *nstab, *nintern;
1702
    char fname[3][AIR_STRLEN_MED];
1703
    nstrn = nrrdNew();
1704
    nstab = nrrdNew();
1705
    nintern = nrrdNew();
1706
    airMopAdd(mop, nstrn, (airMopper)nrrdNuke, airMopAlways);
1707
    airMopAdd(mop, nstab, (airMopper)nrrdNuke, airMopAlways);
1708
    airMopAdd(mop, nintern, (airMopper)nrrdNuke, airMopAlways);
1709
    if (pullInfoGet(nstrn, pullInfoStrength, pctx)
1710
        || pullPropGet(nstab, pullPropStability, pctx)
1711
        || pullPropGet(nintern, pullPropNeighInterNum, pctx)) {
1712
      airMopAdd(mop, err = biffGetDone(PULL), airFree, airMopAlways);
1713
      fprintf(stderr, "%s: trouble 3.2:\n%s", me, err);
1714
      airMopError(mop); return 1;
1715
    }
1716
    sprintf(fname[0], "%s-strn.nrrd", extraOutBaseS);
1717
    sprintf(fname[1], "%s-stab.nrrd", extraOutBaseS);
1718
    sprintf(fname[2], "%s-intern.nrrd", extraOutBaseS);
1719
    if (nrrdSave(fname[0], nstrn, NULL)
1720
        || nrrdSave(fname[1], nstab, NULL)
1721
        || nrrdSave(fname[2], nintern, NULL)) {
1722
      airMopAdd(mop, err = biffGetDone(PULL), airFree, airMopAlways);
1723
      fprintf(stderr, "%s: trouble 3.3:\n%s", me, err);
1724
      airMopError(mop); return 1;
1725
    }
1726
  }
1727
#endif
1728
1729
  pullFinish(pctx);
1730
  airMopOkay(mop);
1731
  return ret;
1732
}