This file is indexed.

/usr/include/teem/gage.h is in libteem-dev 1.11.0~svn6057-1.1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

   1
   2
   3
   4
   5
   6
   7
   8
   9
  10
  11
  12
  13
  14
  15
  16
  17
  18
  19
  20
  21
  22
  23
  24
  25
  26
  27
  28
  29
  30
  31
  32
  33
  34
  35
  36
  37
  38
  39
  40
  41
  42
  43
  44
  45
  46
  47
  48
  49
  50
  51
  52
  53
  54
  55
  56
  57
  58
  59
  60
  61
  62
  63
  64
  65
  66
  67
  68
  69
  70
  71
  72
  73
  74
  75
  76
  77
  78
  79
  80
  81
  82
  83
  84
  85
  86
  87
  88
  89
  90
  91
  92
  93
  94
  95
  96
  97
  98
  99
 100
 101
 102
 103
 104
 105
 106
 107
 108
 109
 110
 111
 112
 113
 114
 115
 116
 117
 118
 119
 120
 121
 122
 123
 124
 125
 126
 127
 128
 129
 130
 131
 132
 133
 134
 135
 136
 137
 138
 139
 140
 141
 142
 143
 144
 145
 146
 147
 148
 149
 150
 151
 152
 153
 154
 155
 156
 157
 158
 159
 160
 161
 162
 163
 164
 165
 166
 167
 168
 169
 170
 171
 172
 173
 174
 175
 176
 177
 178
 179
 180
 181
 182
 183
 184
 185
 186
 187
 188
 189
 190
 191
 192
 193
 194
 195
 196
 197
 198
 199
 200
 201
 202
 203
 204
 205
 206
 207
 208
 209
 210
 211
 212
 213
 214
 215
 216
 217
 218
 219
 220
 221
 222
 223
 224
 225
 226
 227
 228
 229
 230
 231
 232
 233
 234
 235
 236
 237
 238
 239
 240
 241
 242
 243
 244
 245
 246
 247
 248
 249
 250
 251
 252
 253
 254
 255
 256
 257
 258
 259
 260
 261
 262
 263
 264
 265
 266
 267
 268
 269
 270
 271
 272
 273
 274
 275
 276
 277
 278
 279
 280
 281
 282
 283
 284
 285
 286
 287
 288
 289
 290
 291
 292
 293
 294
 295
 296
 297
 298
 299
 300
 301
 302
 303
 304
 305
 306
 307
 308
 309
 310
 311
 312
 313
 314
 315
 316
 317
 318
 319
 320
 321
 322
 323
 324
 325
 326
 327
 328
 329
 330
 331
 332
 333
 334
 335
 336
 337
 338
 339
 340
 341
 342
 343
 344
 345
 346
 347
 348
 349
 350
 351
 352
 353
 354
 355
 356
 357
 358
 359
 360
 361
 362
 363
 364
 365
 366
 367
 368
 369
 370
 371
 372
 373
 374
 375
 376
 377
 378
 379
 380
 381
 382
 383
 384
 385
 386
 387
 388
 389
 390
 391
 392
 393
 394
 395
 396
 397
 398
 399
 400
 401
 402
 403
 404
 405
 406
 407
 408
 409
 410
 411
 412
 413
 414
 415
 416
 417
 418
 419
 420
 421
 422
 423
 424
 425
 426
 427
 428
 429
 430
 431
 432
 433
 434
 435
 436
 437
 438
 439
 440
 441
 442
 443
 444
 445
 446
 447
 448
 449
 450
 451
 452
 453
 454
 455
 456
 457
 458
 459
 460
 461
 462
 463
 464
 465
 466
 467
 468
 469
 470
 471
 472
 473
 474
 475
 476
 477
 478
 479
 480
 481
 482
 483
 484
 485
 486
 487
 488
 489
 490
 491
 492
 493
 494
 495
 496
 497
 498
 499
 500
 501
 502
 503
 504
 505
 506
 507
 508
 509
 510
 511
 512
 513
 514
 515
 516
 517
 518
 519
 520
 521
 522
 523
 524
 525
 526
 527
 528
 529
 530
 531
 532
 533
 534
 535
 536
 537
 538
 539
 540
 541
 542
 543
 544
 545
 546
 547
 548
 549
 550
 551
 552
 553
 554
 555
 556
 557
 558
 559
 560
 561
 562
 563
 564
 565
 566
 567
 568
 569
 570
 571
 572
 573
 574
 575
 576
 577
 578
 579
 580
 581
 582
 583
 584
 585
 586
 587
 588
 589
 590
 591
 592
 593
 594
 595
 596
 597
 598
 599
 600
 601
 602
 603
 604
 605
 606
 607
 608
 609
 610
 611
 612
 613
 614
 615
 616
 617
 618
 619
 620
 621
 622
 623
 624
 625
 626
 627
 628
 629
 630
 631
 632
 633
 634
 635
 636
 637
 638
 639
 640
 641
 642
 643
 644
 645
 646
 647
 648
 649
 650
 651
 652
 653
 654
 655
 656
 657
 658
 659
 660
 661
 662
 663
 664
 665
 666
 667
 668
 669
 670
 671
 672
 673
 674
 675
 676
 677
 678
 679
 680
 681
 682
 683
 684
 685
 686
 687
 688
 689
 690
 691
 692
 693
 694
 695
 696
 697
 698
 699
 700
 701
 702
 703
 704
 705
 706
 707
 708
 709
 710
 711
 712
 713
 714
 715
 716
 717
 718
 719
 720
 721
 722
 723
 724
 725
 726
 727
 728
 729
 730
 731
 732
 733
 734
 735
 736
 737
 738
 739
 740
 741
 742
 743
 744
 745
 746
 747
 748
 749
 750
 751
 752
 753
 754
 755
 756
 757
 758
 759
 760
 761
 762
 763
 764
 765
 766
 767
 768
 769
 770
 771
 772
 773
 774
 775
 776
 777
 778
 779
 780
 781
 782
 783
 784
 785
 786
 787
 788
 789
 790
 791
 792
 793
 794
 795
 796
 797
 798
 799
 800
 801
 802
 803
 804
 805
 806
 807
 808
 809
 810
 811
 812
 813
 814
 815
 816
 817
 818
 819
 820
 821
 822
 823
 824
 825
 826
 827
 828
 829
 830
 831
 832
 833
 834
 835
 836
 837
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849
 850
 851
 852
 853
 854
 855
 856
 857
 858
 859
 860
 861
 862
 863
 864
 865
 866
 867
 868
 869
 870
 871
 872
 873
 874
 875
 876
 877
 878
 879
 880
 881
 882
 883
 884
 885
 886
 887
 888
 889
 890
 891
 892
 893
 894
 895
 896
 897
 898
 899
 900
 901
 902
 903
 904
 905
 906
 907
 908
 909
 910
 911
 912
 913
 914
 915
 916
 917
 918
 919
 920
 921
 922
 923
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
/*
  Teem: Tools to process and visualize scientific data and images             .
  Copyright (C) 2013, 2012, 2011, 2010, 2009  University of Chicago
  Copyright (C) 2008, 2007, 2006, 2005  Gordon Kindlmann
  Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998  University of Utah

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public License
  (LGPL) as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  The terms of redistributing and/or modifying this software also
  include exceptions to the LGPL that facilitate static linking.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public License
  along with this library; if not, write to Free Software Foundation, Inc.,
  51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
*/

#ifndef GAGE_HAS_BEEN_INCLUDED
#define GAGE_HAS_BEEN_INCLUDED

#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <math.h>

#include <teem/air.h>
#include <teem/biff.h>
#include <teem/nrrd.h>
#include <teem/ell.h>

#if defined(_WIN32) && !defined(__CYGWIN__) && !defined(TEEM_STATIC)
#  if defined(TEEM_BUILD) || defined(gage_EXPORTS) || defined(teem_EXPORTS)
#    define GAGE_EXPORT extern __declspec(dllexport)
#  else
#    define GAGE_EXPORT extern __declspec(dllimport)
#  endif
#else /* TEEM_STATIC || UNIX */
#  define GAGE_EXPORT extern
#endif

#ifdef __cplusplus
extern "C" {
#endif

#define GAGE gageBiffKey

/*
** GAGE_DERIV_MAX is the maximum derivative that gage knows how to
** work with. The 0th derivative is the reonstructed value (no
** derivative).  Many short arrays are allocated to 1 + this value
**
** MUST KEEP IN SYNC with GAGE_DV_* macros below
*/
#define GAGE_DERIV_MAX 2

/*
** Convenience macros for arrays declared like blah[GAGE_DERIV_MAX+1]
** (so these need to be kept in sync with GAGE_DERIV_MAX above).  Note
** that "DV" is short for "derivative vector".  With consistent use of
** these, changes to GAGE_DERIV_MAX will either obviously break code
** (in an easy to fix way) or there will be no change needed at all.
**
*/
#define GAGE_DV_SET(v, d0, d1, d2) \
  ((v)[0] = (d0),                  \
   (v)[1] = (d1),                  \
   (v)[2] = (d2))

#define GAGE_DV_EQUAL(a, b)        \
  ((a)[0] == (b)[0] &&             \
   (a)[1] == (b)[1] &&             \
   (a)[2] == (b)[2])

#define GAGE_DV_COPY(a, b)         \
  ((a)[0] = (b)[0],                \
   (a)[1] = (b)[1],                \
   (a)[2] = (b)[2])

/*
** the only extent to which gage treats different axes differently is
** the spacing between samples along the axis.  To have different
** filters for the same function, but along different axes, would be
** too messy.  [Sun Mar 9 13:32:22 EDT 2008: Actually, doing per-axis
** kernels is looking more and more sensible all the time. . .] Thus,
** gage is not very useful as the engine for downsampling: it can't
** tell that along one axis samples should be blurred while they
** should be interpolated along another.  Rather, it assumes that the
** main task of probing is *reconstruction*: of values, of
** derivatives, or lots of different quantities
*/

/*
** terminology: gage is intended to measure multiple things at one
** point in a volume.  The SET of ALL the things that you want
** gage to measure is called the "QUERY".  Each of the many quantities
** comprising the query are called "ITEM"s.
*/

/*
******** gageParm.. enum
**
** these are passed to gageSet.  Look for like-wise named field of
** gageParm for documentation on what these mean.
**
** The following things have to agree:
** - gageParm* enum
** - fields of gageParm struct
** - analogous gageDef* defaults (their declaration and setting)
** - action of gageParmSet (ctx.c)
** - action of gageParmReset (miscGage.c)
*/
enum {
  gageParmUnknown,
  gageParmVerbose,                 /* non-boolean int */
  gageParmRenormalize,             /* int */
  gageParmCheckIntegrals,          /* int */
  gageParmK3Pack,                  /* int */
  gageParmGradMagCurvMin,          /* double */
  gageParmCurvNormalSide,          /* int */
  gageParmKernelIntegralNearZero,  /* double */
  gageParmDefaultCenter,           /* int */
  gageParmStackUse,                /* int */
  gageParmStackNormalizeDeriv,     /* int; does NOT imply gageParmStackUse */
  gageParmStackNormalizeDerivBias, /* double; does NOT imply      "        */
  gageParmStackNormalizeRecon,     /* int; does NOT imply         "        */
  gageParmOrientationFromSpacing,  /* int */
  gageParmGenerateErrStr,          /* int */
  gageParmTwoDimZeroZ,             /* int */
  gageParmLast
};

/* keep in synch with gageErr airEnum */
enum {
  gageErrUnknown,            /* 0: nobody knows */
  gageErrNone,               /* 1: no error, actually, all's well */
  gageErrBoundsSpace,        /* 2: out of 3-D (index-space) bounds */
  gageErrBoundsStack,        /* 3: out of 1-D bounds of stack */
  gageErrStackIntegral,      /* 4: stack recon coeff's sum to 0, so can't
                                   normalize them */
  gageErrStackSearch,        /* 5: for some reason couldn't find the index
                                   position of the probed stack location */
  gageErrStackUnused,        /* 6: can't probe stack without parm.stackUse */
  gageErrLast
};
#define GAGE_ERR_MAX            6

/*
******** gage{Ctx,Pvl}Flag.. enum
**
** organizes all the dependendies within a context and between a
** context and pervolumes.  This logic is to determine the support
** required for a query: different queries need different kernels,
** different kernels have different supports.  The user should not
** have to be concerned about any of this; it should be useful only
** to gageUpdate().
*/
enum {
  gageCtxFlagUnknown,
  gageCtxFlagNeedD,      /*  1: derivatives required for query */
  gageCtxFlagK3Pack,     /*  2: whether to use 3 or 6 kernels */
  gageCtxFlagNeedK,      /*  3: which of the kernels will actually be used */
  gageCtxFlagKernel,     /*  4: any one of the kernels or its parameters */
  gageCtxFlagRadius,     /*  5: radius of support for kernels with query */
  gageCtxFlagShape,      /*  6: a new pervolume shape was set */
  gageCtxFlagLast
};
#define GAGE_CTX_FLAG_MAX    6

enum {
  gagePvlFlagUnknown,
  gagePvlFlagVolume,     /*  1: got a new volume */
  gagePvlFlagQuery,      /*  2: what do you really care about */
  gagePvlFlagNeedD,      /*  3: derivatives required for query */
  gagePvlFlagLast
};
#define GAGE_PVL_FLAG_MAX    3


/*
******** gageKernel.. enum
**
** these are the different kernels that might be used in gage, regardless
** of what kind of volume is being probed.
** (synch with miscGage.c:gageKernel airEnum)
*/
enum {
  gageKernelUnknown,    /*  0: nobody knows */
  gageKernel00,         /*  1: measuring values */
  gageKernel10,         /*  2: reconstructing 1st derivatives */
  gageKernel11,         /*  3: measuring 1st derivatives */
  gageKernel20,         /*  4: reconstructing 1st partials and 2nd deriv.s */
  gageKernel21,         /*  5: measuring 1st partials for a 2nd derivative */
  gageKernel22,         /*  6: measuring 2nd derivatives */
  gageKernelStack,      /*  7: for reconstructing across a stack */
  gageKernelLast
};
#define GAGE_KERNEL_MAX     7

/*
******** GAGE_ITEM_PREREQ_MAXNUM
**
** Max number of prerequisites for any item in *any* kind.
**
** This value has gotten bumped periodically, which used to mean
** that *all* item tables had to be updated, when "-1" was used
** represent the unknown item.  But now that 0 represents the
** unknown item, and because struct elements are implicitly
** initialized to zero, this is no longer the case.
**
** Wed Nov  8 14:12:44 EST 2006 pre-emptively upping this from 6
*/
#define GAGE_ITEM_PREREQ_MAXNUM 8

/*
******** gageItemEntry struct
**
** necessary information about each item supported by the kind, which
** is defined at compile-time in a per-kind table (at least it is for
** the scalar, vector, and tensor kinds)
**
** NOTE!!! YOU CAN NOT re-arrange these variables, because of all the
** compile-time definitions that are done to define the gageKind->table
** for all current kinds.
*/
typedef struct {
  int enumVal;          /* the item's enum value */
  unsigned int
    answerLength;       /* how many double's are needed to store the answer,
                           or (for non-zero items), 0 to represent
                           "the length will be learned later at runtime" */
  int needDeriv,        /* what kind of derivative info is immediately needed
                           for this item (not recursively expanded). This is
                           NO LONGER a bitvector: values are 0, 1, 2, .. */
    prereq[GAGE_ITEM_PREREQ_MAXNUM],
                        /* what are the other items this item depends on
                           (not recusively expanded), you can list up to
                           GAGE_ITEM_PREREQ_MAXNUM of them, and use 0
                           (the unknown item) to fill out the list.
                           _OR_ -1 if this is a dynamic kind and the prereqs
                           will not be known until later in runtime */
    parentItem,         /* the enum value of an item (answerLength > 1)
                           containing the smaller value for which we are
                           merely an alias
                           _OR_ 0 if there's no parent */
    parentIndex,        /* where our value starts in parents value
                           _OR_ 0 if there's no parent */
    needData;           /* whether non-NULL gagePerVolume->data is needed
                           for answering this item */
} gageItemEntry;

/*
** modifying the enums below (scalar, vector, etc query quantities)
** necesitates modifying:
** - the central item table in scl.c, vecGage.c
** - the associated arrays in scl.c, vecGage.c
** - the "answer" method itself in sclanswer.c, vecGage.c
*/

/*
******** gageScl* enum
**
** all the "items" that gage can measure in a scalar volume.
**
** NOTE: although it is currently listed that way, it is not necessary
** that prerequisite measurements are listed before the other measurements
** which need them (that is represented by _gageSclPrereq)
**
** The description for each enum value starts with the numerical value
** followed by a string which identifies the value in the gageScl airEnum.
** The "[N]" indicates how many doubles are used for storing the quantity.
*/
enum {
  gageSclUnknown,      /*  0: nobody knows */
  gageSclValue,        /*  1: "v", data value: [1] */
  gageSclGradVec,      /*  2: "grad", gradient vector, un-normalized: [3] */
  gageSclGradMag,      /*  3: "gm", gradient magnitude: [1] */
  gageSclNormal,       /*  4: "n", gradient vector, normalized: [3] */
  gageSclNProj,        /*  5: "nproj", projection onto normal: [9] */
  gageSclNPerp,        /*  6: "np", projection onto tangent plane: [9] */
  gageSclHessian,      /*  7: "h", Hessian: [9] (column-order) */
  gageSclHessianTen,   /*  8: "ht", Hessian as 7-component tensor: [7]
                           In principle this is a dependency inversion
                           since gage doesn't know about the ten library,
                           where the 7-element tensor is based. */
  gageSclLaplacian,    /*  9: "l", Laplacian: Dxx + Dyy + Dzz: [1] */
  gageSclHessFrob,     /* 10: "hf", Frobenius norm of Hessian: [1] */
  gageSclHessEval,     /* 11: "heval", Hessian's eigenvalues: [3] */
  gageSclHessEval0,    /* 12: "heval0", Hessian's 1st eigenvalue: [1] */
  gageSclHessEval1,    /* 13: "heval1", Hessian's 2nd eigenvalue: [1] */
  gageSclHessEval2,    /* 14: "heval2", Hessian's 3rd eigenvalue: [1] */
  gageSclHessEvec,     /* 15: "hevec", Hessian's eigenvectors: [9] */
  gageSclHessEvec0,    /* 16: "hevec0", Hessian's 1st eigenvector: [3] */
  gageSclHessEvec1,    /* 17: "hevec1", Hessian's 2nd eigenvector: [3] */
  gageSclHessEvec2,    /* 18: "hevec2", Hessian's 3rd eigenvector: [3] */
  gageScl2ndDD,        /* 19: "2d", 2nd dir.deriv. along gradient: [1] */
  gageSclGeomTens,     /* 20: "gten", sym. matx w/ evals {0, K1, K2} and
                              evecs {grad, cdir0, cdir1}: [9] */
  gageSclGeomTensTen,  /* 21: "gtenten", 7-element geometry tensor [7] */
  gageSclK1,           /* 22: "k1", 1st principle curvature: [1] */
  gageSclK2,           /* 23: "k2", 2nd principle curvature (k2 <= k1): [1] */
  gageSclTotalCurv,    /* 24: "tc", L2 norm(K1,K2) (not Koen.'s "C"): [1] */
  gageSclShapeTrace,   /* 25, "st", (K1+K2)/Curvedness: [1] */
  gageSclShapeIndex,   /* 26: "si", Koen.'s shape index, ("S"): [1] */
  gageSclMeanCurv,     /* 27: "mc", mean curvature (K1 + K2)/2: [1] */
  gageSclGaussCurv,    /* 28: "gc", gaussian curvature K1*K2: [1] */
  gageSclCurvDir1,     /* 29: "cdir1", 1st principle curv direction: [3] */
  gageSclCurvDir2,     /* 30: "cdir2", 2nd principle curv direction: [3] */
  gageSclFlowlineCurv, /* 31: "fc", curvature of normal streamline: [1] */
  gageSclMedian,       /* 32: "med", median filter */
  gageSclHessValleyness,  /* 33: "hvalley", vallyness measure: [1] */
  gageSclHessRidgeness,   /* 34: "hridge", ridgeness measure: [1] */
  gageSclHessDotPeakness, /* 35: "hdpeak", positive blobness measure;
                             based on simple dot-product w/ evals : [1] */
  gageSclHessMode,     /* 36: "hmode", Hessian's mode: [1] */
  gageSclLast
};
#define GAGE_SCL_ITEM_MAX 36

/*
******** gageVec* enum
**
** all the "items" that gage knows how to measure in a 3-vector volume
**
** The strings gives one of the gageVec airEnum identifiers, and [x]
** says how many scalars are associated with this value.
*/
enum {
  gageVecUnknown,         /*  0: nobody knows */
  gageVecVector,          /*  1: "v", component-wise-interpolated
                                 (CWI) vec: [3] */
  gageVecVector0,         /*  2: "v0", vector[0]: [1] */
  gageVecVector1,         /*  3: "v1", vector[0]: [1] */
  gageVecVector2,         /*  4: "v2", vector[0]: [1] */
  gageVecLength,          /*  5: "l", length of CWI vector: [1] */
  gageVecNormalized,      /*  6: "n", normalized CWI vector: [3] */
  gageVecJacobian,        /*  7: "j", component-wise Jacobian: [9]
                                0:dv_x/dx  1:dv_x/dy  2:dv_x/dz
                                3:dv_y/dx  4:dv_y/dy  5:dv_y/dz
                                6:dv_z/dx  7:dv_z/dy  8:dv_z/dz */
  gageVecStrain,          /*  8: "S", rate-of-strain tensor: [9] */
  gageVecDivergence,      /*  9: "d", divergence (based on Jacobian): [1] */
  gageVecCurl,            /* 10: "c", curl (based on Jacobian): [3] */
  gageVecCurlNorm,        /* 11: "cm", curl magnitude: [1] */
  gageVecHelicity,        /* 12: "h", helicity: vec . curl: [1] */
  gageVecNormHelicity,    /* 13: "nh", normalized helicity: [1] */
  gageVecSOmega,          /* 14: "somega", S squared + Omega squared: [9] */
  gageVecLambda2,         /* 15: "lambda2", lambda2 criterion: [1] */
  gageVecImaginaryPart,   /* 16: "imag", imag. part of jacobian's
                                  complex-conjugate eigenvalues: [1] */
  gageVecHessian,         /* 17: "vh", second-order derivative: [27]
                                 HEY: indices here need to be double checked
                                 0:d2v_x/dxdx   1:d2v_x/dxdy   2:d2v_x/dxdz
                                 3:d2v_x/dydx   4:d2v_x/dydy   5:d2v_x/dydz
                                 6:d2v_x/dzdx   7:d2v_x/dzdy   8:d2v_x/dzdz
                                 9:d2v_y/dxdx       [..]
                                    [..]
                                24:dv2_z/dzdx  25:d2v_z/dzdy  26:d2v_z/dzdz */
  gageVecDivGradient,     /* 18: "dg", divergence gradient: [3] */
  gageVecCurlGradient,    /* 19: "cg", curl gradient: [9] */
  gageVecCurlNormGrad,    /* 20: "cng", curl norm gradient: [3] */
  gageVecNCurlNormGrad,   /* 21: "ncng", normalized curl norm gradient: [3] */
  gageVecHelGradient,     /* 22: "hg", helicity gradient: [3] */
  gageVecDirHelDeriv,     /* 23: "dhd", directional derivative
                                 of helicity: [1] */
  gageVecProjHelGradient, /* 24: "phg", projected helicity gradient: [3] */
  gageVecGradient0,       /* 25: "g0", gradient of 1st coeff of vector: [3] */
  gageVecGradient1,       /* 26: "g1", gradient of 2nd coeff of vector: [3] */
  gageVecGradient2,       /* 27: "g2", gradient of 3rd coeff of vector: [3] */
  gageVecMultiGrad,       /* 28: "mg", sum of outer products of grads: [9] */
  gageVecMGFrob,          /* 29: "mgfrob", frob norm of multi-gradient: [1] */
  gageVecMGEval,          /* 30: "mgeval", evals of multi-gradient: [3] */
  gageVecMGEvec,          /* 31: "mgevec", evecs of multi-gradient: [9] */
  gageVecLast
};
#define GAGE_VEC_ITEM_MAX     31

struct gageKind_t;       /* dumb forward declaraction, ignore */
struct gagePerVolume_t;  /* dumb forward declaraction, ignore */

/*
******** gageItemPackPart* enum
**
** the different ways that some items can be related for a gagePack
*/
enum {
  gageItemPackPartUnknown,     /*  0 */
  gageItemPackPartScalar,      /*  1 */
  gageItemPackPartGradVec,     /*  2 */
  gageItemPackPartGradMag,     /*  3 */
  gageItemPackPartNormal,      /*  4 */
  gageItemPackPartHessian,     /*  5 */
  gageItemPackPartHessEval0,   /*  6 */
  gageItemPackPartHessEval1,   /*  7 */
  gageItemPackPartHessEval2,   /*  8 */
  gageItemPackPartHessEvec0,   /*  9 */
  gageItemPackPartHessEvec1,   /* 10 */
  gageItemPackPartHessEvec2,   /* 11 */
  gageItemPackPartLast
};
#define GAGE_ITEM_PACK_PART_MAX   11

/*
******** gageShape struct
**
** just a container for all the information related to the "shape"
** of all the volumes associated with a context
**
** Note that the utility of gageShape has extended well beyond doing
** convolution-based measurements in volumes- it has become the
** one-stop place for all of Teem to figure out a reasonable way of
** locating a logically a volume in 3-D space, including using a
** nrrd's full orientation information if it is known.
*/
typedef struct gageShape_t {
  /* ========= INPUT ========= (controls for _gageShapeSet) */
  int defaultCenter,          /* default centering to use when given volume
                                 has no centering set. */
    orientationFromSpacing;   /* only meaningful if nrrd has per-axis spacing,
                                 but not full orientation info. If zero, the
                                 volume is crammed into the bi-unit cube.
                                 If non-zero, gage treats the volume as if it
                                 had axis-aligned spaceDirection vectors, with
                                 the non-zero values determined by the given
                                 per-axis spacing. */
  /* ======== OUTPUT ========= (set by _gageShapeSet) */
  int center,                 /* the sample centering of the volume(s)- this
                                 determines the extent of the locations
                                 that may be probed */
    fromOrientation;          /* non-zero iff the spaceDirections and
                                 spaceOrigin information was used */
  unsigned int size[3];       /* raster dimensions of volume
                                 HEY: shouldn't this be size_t? */
  double spacing[3];          /* spacings for each axis */
  /* fwScale[GAGE_KERNEL_MAX+1][3] has been superceded by the cleaner and
     more general ItoWSubInvTransp and ItoWSubInv matrices below */
  double ItoW[16],            /* homogeneous coord transform from index
                                 to world space (defined either by bi-unit
                                 cube or from full orientation info ) */
    WtoI[16],                 /* inverse of above */
    ItoWSubInvTransp[9],      /* inverse transpose of 3x3 sub-matrix of ItoW,
                                 to transform (covariant) gradients */
    ItoWSubInv[9];            /* tranpose of above, to transform hessians */

} gageShape;

/*
******** gageParm struct
**
** a container for the various switches and knobs which control
** gage, aside from the obvious inputs (kernels, queries, volumes)
*/
typedef struct gageParm_t {
  int renormalize;            /* hack to make sure that sum of
                                 discrete value reconstruction weights
                                 is same as kernel's continuous
                                 integral, and that the 1nd and 2nd
                                 deriv weights really sum to 0.0 */
  int checkIntegrals;         /* call the "integral" method of the
                                 kernel to verify that it is
                                 appropriate for the task for which
                                 the kernel is being set:
                                 reconstruction: 1.0, derivatives: 0.0 */
  int k3pack;                 /* non-zero (true) iff we do not use kernels
                                 gageKernelIJ with I != J. So, we use the
                                 value reconstruction kernel (gageKernel00)
                                 for 1st and 2nd derivative reconstruction,
                                 and so on. This is faster because we can
                                 re-use results from low-order convolutions.
                                 The name "3pack" will likely persist even
                                 with 3rd dervatives, for which "4pack" would
                                 make more sense (for 00, 11, 22, 33) */
  double gradMagCurvMin,      /* punt on computing curvature information if
                                 gradient magnitude is less than this. Yes,
                                 this is scalar-kind-specific, but there's
                                 no other good place for it */
    kernelIntegralNearZero,   /* tolerance with checkIntegrals on derivative
                                 kernels */
    stackNormalizeDerivBias;  /* when using stackNormalizeDeriv, a bias
                                 on the effective scale, used for the
                                 normalization */
  int curvNormalSide,         /* determines direction of gradient that is used
                                 as normal in curvature calculations, exactly
                                 the same as miteUser's normalSide: 1 for
                                 normal pointing to lower values (higher
                                 values are more "inside"); -1 for normal
                                 pointing to higher values (low values more
                                 "inside") */
    defaultCenter,            /* what centering to use when you have to invent
                                 one, because its not set on any axis */
    stackUse,                 /* if non-zero: treat the pvl's (all same kind)
                                 as multiple values of a single logical volume
                                 (e.g. for approximating scale space).
                                 The first pvl is effectively just a buffer;
                                 the N-1 pvls used are at index 1 through N-2.
                                 The query in pvl[0] will determine the
                                 computations done, and answers for the whole
                                 stack will be stored in pvl[0]. */
    stackNormalizeRecon,      /* if non-zero (and if stackUse is non-zero):
                                 the weights used to reconstruct across
                                 stack samples are normalized to unit sum; not
                                 needed if the kernel is accurate enough */
    stackNormalizeDeriv,      /* if non-zero (and if stackUse is non-zero):
                                 derivatives at filter stage (before answer
                                 stage) are renormalized based on the stack
                                 position */
    orientationFromSpacing,   /* only meaningful if nrrd has per-axis spacing,
                                 but not full orientation info. If zero, the
                                 volume is crammed into the bi-unit cube.
                                 If non-zero, gage treats the volume as if it
                                 had axis-aligned spaceDirection vectors, with
                                 the non-zero values determined by the given
                                 per-axis spacing. */
    generateErrStr,           /* when errors happen, biff is never used, but
                                 a descriptive error is sprintf into
                                 gctx->errStr as long as this is non-zero. */
    twoDimZeroZ;              /* a limited way of supporting queries on
                                 two-dimensional images: if this is non-zero
                                 then *some* answers will only involve the 1st
                                 ("X") and 2nd ("Y") coordinates of world
                                 space.  Eigensystems should have only two
                                 two elements, with the 3rd being NaN'd out.
                                 Because gage has always been only about 3D
                                 images; the implementation of this is likely
                                 incomplete, since the responsibility for
                                 correctly handling it ultimately falls to the
                                 "answer" functions of the various
                                 gageKinds */
} gageParm;

/*
******** gagePoint struct
**
** stores *Index Space* location of last query location, which is used
** to determine whether the ctx->fsl, ctx->fw values can be re-used
** (based on the "Frac" values), and, whether all the pvl->iv3 have to
** be refilled (based on the "Idx" values).  The last index (frac[3]
** and idx[3])is for the stack, and can safely stay 0 if the stack
** isn't being used.
**
** with stack usage, stackFwNonZeroNum records how many pvls had
** non-zero stack filter weights, which is used to detect when
** iv3s have to be refilled.  Looking at idx[3] alone is not always
** sufficient for this.
*/
typedef struct gagePoint_t {
  double frac[4];         /* last fractional voxel location */
  unsigned int idx[4],    /* last integral voxel location; actually the
                             *upper* corner of the containing voxel */
    stackFwNonZeroNum;    /* last number of non-zero values of stack filter
                             weights (ctx->stackFw) */
} gagePoint;

/*
******** gageQuery typedef
**
** this short, fixed-length array is used as a bit-vector to store
** all the items that comprise a query.  Its length sets an upper
** bound on the maximum item value that a gageKind may use.
**
** The important thing to keep in mind is that a variable of type
** gageKind ends up being passed by reference, even though the
** syntax of its usage doesn't immediately announce that.
**
** gageKindCheck currently has the role of verifying that a gageKind's
** maximum item value is within the bounds set here. Using
** GAGE_QUERY_BYTES_NUM == 8 gives a max item value of 63, which is
** far above anything being used now.
**
** Sat Jan 21 18:12:01 EST 2006: ha! second derivatives of tensors blew
** past old GAGE_QUERY_BYTES_NUM, now GAGE_QUERY_BYTES_NUM == 16
**
** Tue Nov  7 19:51:05 EST 2006; tenGage items now pushing 127,
** guardedly changing GAGE_QUERY_BYTES_NUM to 24 --> max item 191
**
** Wed Nov 18 17:53:23 CST 2009; yikes, tenGage items now at 190,
** changing GAGE_QUERY_BYTES_NUM to 32 --> max item 255
**
** NOTE: increasing GAGE_QUERY_BYTES_NUM means that the macros below
** have to be redefined as well!
*/
#define GAGE_QUERY_BYTES_NUM 32
#define GAGE_ITEM_MAX ((8*GAGE_QUERY_BYTES_NUM)-1)
typedef unsigned char gageQuery[GAGE_QUERY_BYTES_NUM];

/*
******** GAGE_QUERY_RESET, GAGE_QUERY_TEST
******** GAGE_QUERY_ON, GAGE_QUERY_OFF
**
** Macros for manipulating a gageQuery.
**
** airSanity ensures that an unsigned char is in fact 8 bits
*/
#define GAGE_QUERY_RESET(q) \
  q[ 0] = q[ 1] = q[ 2] = q[ 3] = \
  q[ 4] = q[ 5] = q[ 6] = q[ 7] = \
  q[ 8] = q[ 9] = q[10] = q[11] = \
  q[12] = q[13] = q[14] = q[15] = \
  q[16] = q[17] = q[18] = q[19] = \
  q[20] = q[21] = q[22] = q[23] = \
  q[24] = q[25] = q[26] = q[27] = \
  q[28] = q[29] = q[30] = q[31] = 0

#define GAGE_QUERY_COPY(p, q) \
  p[ 0] = q[ 0]; p[ 1] = q[ 1]; p[ 2] = q[ 2]; p[ 3] = q[ 3]; \
  p[ 4] = q[ 4]; p[ 5] = q[ 5]; p[ 6] = q[ 6]; p[ 7] = q[ 7]; \
  p[ 8] = q[ 8]; p[ 9] = q[ 9]; p[10] = q[10]; p[11] = q[11]; \
  p[12] = q[12]; p[13] = q[13]; p[14] = q[14]; p[15] = q[15]; \
  p[16] = q[16]; p[17] = q[17]; p[18] = q[18]; p[19] = q[19]; \
  p[20] = q[20]; p[21] = q[21]; p[22] = q[22]; p[23] = q[23]; \
  p[24] = q[24]; p[25] = q[25]; p[26] = q[26]; p[27] = q[27]; \
  p[28] = q[28]; p[29] = q[29]; p[30] = q[30]; p[31] = q[31]

#define GAGE_QUERY_ADD(p, q) \
  p[ 0] |= q[ 0]; p[ 1] |= q[ 1]; p[ 2] |= q[ 2]; p[ 3] |= q[ 3]; \
  p[ 4] |= q[ 4]; p[ 5] |= q[ 5]; p[ 6] |= q[ 6]; p[ 7] |= q[ 7]; \
  p[ 8] |= q[ 8]; p[ 9] |= q[ 9]; p[10] |= q[10]; p[11] |= q[11]; \
  p[12] |= q[12]; p[13] |= q[13]; p[14] |= q[14]; p[15] |= q[15]; \
  p[16] |= q[16]; p[17] |= q[17]; p[18] |= q[18]; p[19] |= q[19]; \
  p[20] |= q[20]; p[21] |= q[21]; p[22] |= q[22]; p[23] |= q[23]; \
  p[24] |= q[24]; p[25] |= q[25]; p[26] |= q[26]; p[27] |= q[27]; \
  p[28] |= q[28]; p[29] |= q[29]; p[30] |= q[30]; p[31] |= q[31]

#define GAGE_QUERY_EQUAL(p, q) ( \
  p[ 0] == q[ 0] && p[ 1] == q[ 1] && p[ 2] == q[ 2] && p[ 3] == q[ 3] && \
  p[ 4] == q[ 4] && p[ 5] == q[ 5] && p[ 6] == q[ 6] && p[ 7] == q[ 7] && \
  p[ 8] == q[ 8] && p[ 9] == q[ 9] && p[10] == q[10] && p[11] == q[11] && \
  p[12] == q[12] && p[13] == q[13] && p[14] == q[14] && p[15] == q[15] && \
  p[16] == q[16] && p[17] == q[17] && p[18] == q[18] && p[19] == q[19] && \
  p[20] == q[20] && p[21] == q[21] && p[22] == q[22] && p[23] == q[23] && \
  p[24] == q[24] && p[25] == q[25] && p[26] == q[26] && p[27] == q[27] && \
  p[28] == q[28] && p[29] == q[29] && p[30] == q[30] && p[31] == q[31] )

#define GAGE_QUERY_NONZERO(q) ( \
  q[ 0] | q[ 1] | q[ 2] | q[ 3] | \
  q[ 4] | q[ 5] | q[ 6] | q[ 7] | \
  q[ 8] | q[ 9] | q[10] | q[11] | \
  q[12] | q[13] | q[14] | q[15] | \
  q[16] | q[17] | q[18] | q[19] | \
  q[20] | q[21] | q[22] | q[23] | \
  q[24] | q[25] | q[26] | q[27] | \
  q[28] | q[29] | q[30] | q[31] )

#define GAGE_QUERY_ITEM_TEST(q, i) (q[i/8] & (1 << (i % 8)))
#define GAGE_QUERY_ITEM_ON(q, i) (q[i/8] |= (1 << (i % 8)))
#define GAGE_QUERY_ITEM_OFF(q, i) (q[i/8] &= ~(1 << (i % 8)))

/* increment for ctx->pvlArr airArray */
#define GAGE_PERVOLUME_ARR_INCR 32

/* extents of known information about optimal sigma samples for
   Hermite-spline-based scale-space reconstruction */
#define GAGE_OPTIMSIG_SIGMA_MAX 11
#define GAGE_OPTIMSIG_SAMPLES_MAXNUM 11

/*
******** gageContext struct
**
** The information here is specific to the dimensions, scalings, and
** radius of kernel support, but not to kind of volume (all kind-specific
** information is in the gagePerVolume).  One context can be used in
** conjuction with probing multiple volumes.
*/
typedef struct gageContext_t {
  /* INPUT ------------------------- */
  int verbose;                /* verbosity */
  gageParm parm;              /* all parameters */

  /* all the kernels we'll ever need, including the stack kernel */
  NrrdKernelSpec *ksp[GAGE_KERNEL_MAX+1];

  /* all the pervolumes attached to this context.  If using stack,
     the base pvl is the LAST, pvl[pvlNum-1], and the stack samples
     are pvl[0] through pvl[pvlNum-2] */
  struct gagePerVolume_t **pvl;

  /* number of pervolumes currently attached. If using stack,
     this is one more than number of stack samples (because of the
     base volume at the end) */
  unsigned int pvlNum;

  /* airArray for managing pvl and pvlNum */
  airArray *pvlArr;

  /* sizes, spacings, centering, and other geometric aspects of the
     volume */
  gageShape *shape;

  /* if stack is being used, allocated for length pvlNum-1, and
     stackPos[0] through stackPos[pvlNum-2] MUST exist and be
     monotonically increasing stack positions for each volume.
     Otherwise NULL */
  double *stackPos;

  /* INTERNAL ------------------------- */
  /* if using stack: allocated for length pvlNum-1, and filter sample
     locations and weights for reconstruction along the stack.
     Otherwise NULL. */
  double *stackFsl;
  double *stackFw;

  /* all the flags used by gageUpdate() used to describe what changed
     in this context */
  int flag[GAGE_CTX_FLAG_MAX+1];

  /* which value/derivatives need to be calculated for all pervolumes
     (doV, doD1, doD2) */
  int needD[GAGE_DERIV_MAX+1];

  /* which kernels are needed for all pvls.  needK[gageKernelStack]
     is currently not set by the update function that sets needK[] */
  int needK[GAGE_KERNEL_MAX+1];

  /* radius of support of samples needed to satisfy query, given the
     set of kernels.  The "filter diameter" fd == 2*radius.  This is
     incremented by one when filtering across the stack with
     nrrdKernelHermiteScaleSpaceFlag */
  unsigned int radius;

  /* filter sample locations (all axes): logically a fd x 3 array
     (and its 3 because gage works on 3D fields, not because of
     the number of derivatives supported) */
  double *fsl;

  /* filter weights (all axes, all kernels): logically a
     fd x 3 x GAGE_KERNEL_MAX+1 (fast-to-slow) array */
  double *fw;

  /* offsets to other fd^3 samples needed to fill 3D intermediate
     value cache. Allocated size is dependent on kernels, values
     inside are dependent on the dimensions of the volume. It may be
     more correct to be using size_t instead of uint, but the X and Y
     dimensions of the volume would have to be super-outrageous for
     that to be a problem */
  unsigned int *off;

  /* last probe location */
  gagePoint point;

  /* errStr and errNum are for describing errors that happen in gageProbe():
     using biff is too heavy-weight for this, and the idea is that no ill
     should occur if the error is repeatedly ignored. Whether or not
     something is actually sprintf'ed into errStr is controlled by
     parm.generateErrStr.
     NOTE: these variables used to be globals "gageErrStr" and "gageErrNum" */
  char errStr[AIR_STRLEN_LARGE];
  int errNum;                 /* takes values from the gageErr enum */

  /* what fraction of the values in the kernel support had to be invented
     (by bleeding out the margin) in order to satisfy a probe that was near
     the edge (any axis, either high or low) of the volume. However, this
     value is NOT meaningfully set if there is no clamping, and the probe
     location as fallen outside the volume */
  double edgeFrac;
} gageContext;

/*
******** gagePerVolume
**
** information that is specific to one volume, and to one kind of
** volume.
*/
typedef struct gagePerVolume_t {
  int verbose;                /* verbosity */
  const struct gageKind_t *kind;  /* what kind of volume is this for */
  gageQuery query;            /* the query, recursively expanded */
  int needD[GAGE_DERIV_MAX+1];/* which derivatives need to be calculated for
                                 the query (above) in this volume */
  const Nrrd *nin;            /* the volume to sample within */
  int flag[GAGE_PVL_FLAG_MAX+1];/* for the kind-specific flags .. */
  double *iv3, *iv2, *iv1;    /* 3D, 2D, 1D, value caches.  These are cubical,
                                 square, and linear arrays, all length fd on
                                 each edge.  Currently gageIv3Fill() fills
                                 the iv3 (at a latter point this will be
                                 delegated back to the gageKind to facilitate
                                 bricking), and currently the tuple axis (with
                                 length valLen) always slowest.  However, use
                                 of iv2 and iv1 is entirely up the kind's
                                 filter method. */
  double (*lup)(const void *ptr, size_t I);
                              /* nrrd{F,D}Lookup[] element, according to
                                 nin->type and double */
  double *answer;             /* main buffer to hold all the answers */
  double **directAnswer;      /* array of pointers into answer */
  void *data;                 /* extra data, parameters, buffers, etc.
                                 required for answering some items (as per
                                 the gageItemEntry->needData) managed with
                                 kind->pvlDataNew, kind->pvlDataCopy,
                                 kind->pvlDataUpdate, and kind->pvlDataNix,
                                 so there is no channel for extra info to be
                                 passed into the pvl->data, other that what
                                 was put into kind->data */
} gagePerVolume;

/*
******** gageKind struct
**
** all the information and functions that are needed to handle one
** kind of volume (such as scalar, vector, etc.)
**
** these are either statically allocated (e.g. gageKindScl, gageKindVec,
** tenGageKind), or dynamically allocated, which case the kind itself
** needs a constructor (e.g. tenDwiGageKindNew()).  The "dynamicAlloc"
** variable indicates this distinction.
**
** Having dynamically allocated kinds raises the possibility of having
** to set and update (or modify and update) their internal state,
** which is currently completely outside the update framework of gage.
** So as far as the core gage functions are concerned, all kinds are
** static, because there is nothing to modify.  It also means that
** those who dynamically create kinds should try to minimize the
** state info that can/must be dynamically modified (i.e. maybe
** the kind constructor should take all the various parameters,
** and set everything in a single shot).
*/
typedef struct gageKind_t {
  int dynamicAlloc;                 /* non-zero if this kind struct was
                                       dynamically allocated */
  char name[AIR_STRLEN_SMALL];      /* short identifying string for kind */
  const airEnum *enm;               /* such as gageScl.  NB: the "unknown"
                                       value in the enum *must* be 0. */
  unsigned int baseDim,             /* dimension that x,y,z axes start on
                                       (e.g. 0 for scalars, 1 for vectors) */
    valLen;                         /* number of scalars per data point,
                                       -or- 0 to represent "this value will
                                       be learned later at runtime" */
  int itemMax;                      /* such as GAGE_SCL_ITEM_MAX */
  gageItemEntry *table;             /* array of gageItemEntry's, indexed
                                       by the item value,
                                       -or- NULL if the table cannot be
                                       statically allocated (not because it
                                       can come in different sizes, but
                                       because it needs to be a modified
                                       version of the compile-time table */
  void (*iv3Print)(FILE *,          /* such as _gageSclIv3Print() */
                   gageContext *,
                   gagePerVolume *),
    (*filter)(gageContext *,        /* such as _gageSclFilter() */
              gagePerVolume *),
    (*answer)(gageContext *,        /* such as _gageSclAnswer() */
              gagePerVolume *),
    /* for allocating, copying, updating, and nixing the pervolume->data */
    /* pvlDataNew and pvlDataCopy can use biff, but:
       --> they must use GAGE key (and not callback's library's key), and
       --> pvlDataNix can not use biff */
    *(*pvlDataNew)(const struct gageKind_t *kind),
    *(*pvlDataCopy)(const struct gageKind_t *kind, const void *pvlDataOld),
    *(*pvlDataNix)(const struct gageKind_t *kind, void *pvlData);
  int (*pvlDataUpdate)(const struct gageKind_t *kind,
                       const gageContext *ctx,
                       const gagePerVolume *pvl,
                       const void *data);
  void *data;                       /* extra information about the kind of
                                       volume that's being probed.  Likely
                                       used by filter, answer, and the
                                       pvlData functions */
} gageKind;

/*
******** gageItemSpec struct
**
** dumb little way to store a kind/item pair.  Formerly known
** as a gageQuerySpec.
*/
typedef struct {
  const gageKind *kind;       /* the kind of the volume to measure */
  int item;                   /* the quantity to measure */
} gageItemSpec;

/*
******** gageItemPack struct
**
** A way of coherently representing a related set of gageItems. This
** kind of inter-item relationship may eventually be part of the
** definition of a gageKind.
**
** These are intended to be set up at compile-time, just like nearly all
** the gageKinds.
**
** This is not to be confused with the gageMultiItem, which is (or
** will be) a list of items, set at run-time, to learn from a given
** volume.
*/
typedef struct {
  const gageKind *kind;
  int item[GAGE_ITEM_PACK_PART_MAX+1];
} gageItemPack;

/*
******** gageStackBlurParm struct
**
** all parameters associated with blurring one volume to form a "stack"
*/
typedef struct {
  unsigned int num;      /* # of blurring scales == # volumes */
  double *scale;         /* scale parameter for each blurred volume */
  double sigmaMax,       /* nrrdKernelDiscreteGaussian is implemented with
                            airBesselInExpScaled, which has numerical issues
                            for very large kernel sizes.  Instead of doing
                            the blurring in one step, the diffusion is done
                            iteratively, with steps in diffusion time
                            of sigmaMax^2 */
    padValue;            /* padding value for nrrdBoundaryPad */
  NrrdKernelSpec *kspec; /* NOTE: parm[0] will get over-written as part
                            of running the resampler at each scale */
  int dataCheck,         /* when checking given stack to see if its the
                            blurring of a volume, actually go in and look at
                            values of first (probably the least blurred)
                            volume */
    boundary,            /* passed to nrrdResampleBoundarySet */
    renormalize,         /* passed to nrrdResampleRenormalizeSet */
    oneDim,              /* for experimental purposes: blur *only* along
                            the first (fastest) axis */
    verbose;
} gageStackBlurParm;

/*
******** gageOptimSigContext struct
**
** The parameters and state needed to optimize (via gageOptimSigCalculate) the
** locations for reconstructing scale-space, in the specific case of using
** nrrdKernelDiscreteGaussian for blurring, and using the
** nrrdKernelHermiteScaleSpaceFlag for reconstructing along scale.  The
** "samples" are the pre-blurred correct blurrings along scale, of which there
** are usually about a dozen or fewer.  The purpose of gageOptimSigCalculate
** is to find where to put those samples so that the error of reconstructing
** at any other scale is in some sense minimized.
**
** This code was re-written August 2013; one of the big changes is that there
** is no longer a need for a big 4-D array of pre-computed blurrings.
*/
typedef struct {
  /* INPUT ------------------------- */
  /* these are set (once) at context creation time */
  unsigned int dim,        /* what dimension of point-spread function to
                              optimize based on; either 1, 2, or 3 */
    sampleNumMax,          /* max number of samples to optimize */
    trueImgNum;            /* how many samples along scale to use for the
                              correct reference blurrings; this discretization
                              determines the accuracy of the error measurement
                              integrated across scales. For allMeasr *other*
                              than nrrdMeasureLinf, this determines the number
                              of error values summarized with allMeasr. With
                              the re-write of this code, setting this is not
                              as consequential as it used to be; now it only
                              controls the allocation of nerr.  One could
                              imagine another tweak that shifted truImgNum
                              to the gageOptimSigCalculate or similar call */
  double sigmaRange[2];    /* range of sigma values that should be studied */
  double cutoff;           /* second parm to nrrdKernelDiscreteGaussian */
  /* NOTE: the image blurring to sample a particular scale is always by
     nrrdKernelDiscreteGaussian; this code is not built for other kinds
     blurring kernels, even though that would be interesting to pursue */

  /* these are set with calls to gageOptimSigCalculate or gageOptimSigPlot */
  NrrdKernelSpec *kssSpec; /* how to interpolate across scale; should be
                              nrrdKernelHermiteScaleSpaceFlag for
                              kspec->kernel; the main purpose of the sigma
                              optimization that this code does */
  unsigned int sampleNum,  /* how many scale samples to optimize */
    maxIter;               /* allowed iterations in optimization */
  int imgMeasr,            /* how to measure error at each reconstructed
                              scale (in the scale-interpolated image) */
    allMeasr;              /* how to summarize errors across all scales */
  double convEps;          /* convergence threshold */

  /* INTERNAL ------------------------- */
  /* NOTE: all internal computations are parameterized by a tau-like
     quantity termed rho, rather than sigma */
  unsigned int sx, sy, sz; /* image size for testing; determined by
                              dim, sigmaRange[1], and cutoff */
  Nrrd *nerr,              /* 1D array of all errors, across scale, with
                              length trueImgNum */
    *ninterp,              /* last scale interpolation result */
    *ndiff;                /* diff between truth and the single recon last
                              computed in ninterp */
  double rhoRange[2],      /* rho(sigmaRange) */
    *kloc,                 /* allocated for length sx to evaluation locations
                              of nrrdKernelDiscreteGaussian */
    *kern, *ktmp1, *ktmp2, /* allocated for length sx to store a high-quality
                              nrrdKernelDiscreteGaussian evaluation, or
                              buffers related to that */
    kone[1];               /* stores 1.0 */
  gageContext *gctx;       /* context around pvlBase, pvlSS, and nsampleImg */
  /* buffers for kernel evaluation */
  /* except for pvlBase, these are allocated for sampleNumMax */
  gagePerVolume *pvlBase,  /* the base pvl for getting answers */
    **pvlSS;               /* for gage; pvlSS is the stack pervolume,
                              for the sampleNum volumes in nsampvol */
  Nrrd **nsampleImg;       /* current set of scale-interpolated samples */
  double *sampleSigma,     /* current locations of nsampleImg in sigma
                              (this array is for the gageStack) */
    *sampleRho,            /* current locations of nsampleImg in rho */
    *sampleTmp,            /* buffer for sample location info */
    *sampleErrMax,         /* for tracking per-sample-pair errors in Linf */
    *step;                 /* per-point stepsize for gradient descent */

  /* OUTPUT ------------------------- */
  double finalErr;         /* error of converged points */
} gageOptimSigContext;

/* defaultsGage.c */
GAGE_EXPORT const char *gageBiffKey;
GAGE_EXPORT int gageDefVerbose;
GAGE_EXPORT double gageDefGradMagCurvMin;
GAGE_EXPORT int gageDefRenormalize;
GAGE_EXPORT int gageDefCheckIntegrals;
GAGE_EXPORT int gageDefK3Pack;
GAGE_EXPORT int gageDefCurvNormalSide;
GAGE_EXPORT double gageDefKernelIntegralNearZero;
GAGE_EXPORT int gageDefDefaultCenter;
GAGE_EXPORT int gageDefStackUse;
GAGE_EXPORT int gageDefStackNormalizeRecon;
GAGE_EXPORT int gageDefStackNormalizeDeriv;
GAGE_EXPORT double gageDefStackNormalizeDerivBias;
GAGE_EXPORT double gageDefStackBlurSigmaMax;
GAGE_EXPORT int gageDefOrientationFromSpacing;
GAGE_EXPORT int gageDefGenerateErrStr;
GAGE_EXPORT int gageDefTwoDimZeroZ;

/* miscGage.c */
GAGE_EXPORT const int gagePresent;
GAGE_EXPORT double gageZeroNormal[3];
GAGE_EXPORT const airEnum *const gageErr;
GAGE_EXPORT const airEnum *const gageKernel;
GAGE_EXPORT const airEnum *const gageItemPackPart;
GAGE_EXPORT void gageParmReset(gageParm *parm);
GAGE_EXPORT void gagePointReset(gagePoint *point);
GAGE_EXPORT gageItemSpec *gageItemSpecNew(void);
GAGE_EXPORT void gageItemSpecInit(gageItemSpec *isp);
GAGE_EXPORT gageItemSpec *gageItemSpecNix(gageItemSpec *isp);

/* kind.c */
GAGE_EXPORT int gageKindCheck(const gageKind *kind);
GAGE_EXPORT unsigned int gageKindTotalAnswerLength(const gageKind *kind);
GAGE_EXPORT unsigned int gageKindAnswerLength(const gageKind *kind, int item);
GAGE_EXPORT int gageKindAnswerOffset(const gageKind *kind, int item);
GAGE_EXPORT int gageKindVolumeCheck(const gageKind *kind, const Nrrd *nrrd);

/* print.c */
GAGE_EXPORT void gageQueryPrint(FILE *file, const gageKind *kind,
                                gageQuery query);

/* sclfilter.c */
typedef void (gageScl3PFilter_t)(gageShape *shape,
                                 double *iv3, double *iv2, double *iv1,
                                 double *fw00, double *fw11, double *fw22,
                                 double *val, double *gvec, double *hess,
                                 const int *needD);
GAGE_EXPORT gageScl3PFilter_t gageScl3PFilter2;
GAGE_EXPORT gageScl3PFilter_t gageScl3PFilter4;
GAGE_EXPORT gageScl3PFilter_t gageScl3PFilter6;
GAGE_EXPORT gageScl3PFilter_t gageScl3PFilter8;
GAGE_EXPORT void gageScl3PFilterN(gageShape *shape, int fd,
                                  double *iv3, double *iv2, double *iv1,
                                  double *fw00, double *fw11, double *fw22,
                                  double *val, double *gvec, double *hess,
                                  const int *needD);

/* scl.c */
GAGE_EXPORT const airEnum *const gageScl;
/* HEY: this should perhaps be a "const gageKind", but pointers for
   general kinds may want to point to this, or to the return of
   a dynamic kind generator like tenDwiGageKindNew(), which is
   most certainly not "const gageKind". */
GAGE_EXPORT gageKind *const gageKindScl;
GAGE_EXPORT const gageItemPack *const gageItemPackSclValue;

/* vecGage.c (together with vecprint.c, these contain everything to
   implement the "vec" kind, and could be used as examples of what it
   takes to create a new gageKind) */
GAGE_EXPORT const airEnum *const gageVec;
GAGE_EXPORT gageKind *const gageKindVec;

/* shape.c */
GAGE_EXPORT void gageShapeReset(gageShape *shp);
GAGE_EXPORT gageShape *gageShapeNew(void);
GAGE_EXPORT gageShape *gageShapeCopy(const gageShape *shp);
GAGE_EXPORT gageShape *gageShapeNix(gageShape *shape);
GAGE_EXPORT int gageShapeSet(gageShape *shp, const Nrrd *nin, int baseDim);
GAGE_EXPORT void gageShapeWtoI(const gageShape *shape,
                               double index[3], const double world[3]);
GAGE_EXPORT void gageShapeItoW(const gageShape *shape,
                               double world[3], const double index[3]);
GAGE_EXPORT int gageShapeEqual(const gageShape *shp1, const char *name1,
                               const gageShape *shp2, const char *name2);
GAGE_EXPORT void gageShapeBoundingBox(double min[3], double max[3],
                                      const gageShape *shape);

/* the organization of the next two files used to be according to
   what the first argument is, not what appears in the function name,
   but that's just a complete mess now */
/* pvl.c */
GAGE_EXPORT int gageVolumeCheck(const gageContext *ctx, const Nrrd *nin,
                                const gageKind *kind);
GAGE_EXPORT gagePerVolume *gagePerVolumeNew(gageContext *ctx,
                                            const Nrrd *nin,
                                            const gageKind *kind);
GAGE_EXPORT gagePerVolume *gagePerVolumeNix(gagePerVolume *pvl);
GAGE_EXPORT const double *gageAnswerPointer(const gageContext *ctx,
                                            const gagePerVolume *pvl,
                                            int item);
GAGE_EXPORT unsigned int gageAnswerLength(const gageContext *ctx,
                                          const gagePerVolume *pvl,
                                          int item);
GAGE_EXPORT int gageQueryReset(gageContext *ctx, gagePerVolume *pvl);
GAGE_EXPORT int gageQuerySet(gageContext *ctx, gagePerVolume *pvl,
                             gageQuery query);
GAGE_EXPORT int gageQueryAdd(gageContext *ctx, gagePerVolume *pvl,
                             gageQuery query);
GAGE_EXPORT int gageQueryItemOn(gageContext *ctx, gagePerVolume *pvl,
                                int item);

/* optimsig.c */
GAGE_EXPORT int gageOptimSigSet(double *scale, unsigned int num,
                                unsigned int sigmaMax);
GAGE_EXPORT gageOptimSigContext *
  gageOptimSigContextNew(unsigned int dim,
                         unsigned int sampleNumMax,
                         unsigned int trueImgNum,
                         double sigmaMin, double sigmaMax,
                         double cutoff);
GAGE_EXPORT gageOptimSigContext *gageOptimSigContextNix(gageOptimSigContext
                                                        *oscx);
GAGE_EXPORT int gageOptimSigCalculate(gageOptimSigContext *oscx,
                                      /* output */ double *sigma,
                                      unsigned int sigmaNum,
                                      const NrrdKernelSpec *kssSpec,
                                      int imgMeasr, int allMeasr,
                                      unsigned int maxIter, double convEps);
GAGE_EXPORT int gageOptimSigErrorPlot(gageOptimSigContext *oscx,
                                      Nrrd *nout,
                                      const double *sigma,
                                      unsigned int sigmaNum,
                                      const NrrdKernelSpec *kssSpec,
                                      int imgMeasr);
GAGE_EXPORT int gageOptimSigErrorPlotSliding(gageOptimSigContext *oscx,
                                             Nrrd *nout,
                                             double windowRho,
                                             unsigned int sampleNum,
                                             const NrrdKernelSpec *kssSpec,
                                             int imgMeasr);

/* stack.c */
GAGE_EXPORT double gageTauOfTee(double tee);
GAGE_EXPORT double gageTeeOfTau(double tau);
GAGE_EXPORT double gageSigOfTau(double tau);
GAGE_EXPORT double gageTauOfSig(double sig);
GAGE_EXPORT double gageStackWtoI(gageContext *ctx, double swrl,
                                 int *outside);
GAGE_EXPORT double gageStackItoW(gageContext *ctx, double si,
                                 int *outside);
GAGE_EXPORT int gageStackPerVolumeNew(gageContext *ctx,
                                      gagePerVolume **pvlStack,
                                      const Nrrd *const *nblur,
                                      unsigned int blnum,
                                      const gageKind *kind);
GAGE_EXPORT int gageStackPerVolumeAttach(gageContext *ctx,
                                         gagePerVolume *pvlBase,
                                         gagePerVolume **pvlStack,
                                         const double *stackPos,
                                         unsigned int blnum);
GAGE_EXPORT int gageStackProbe(gageContext *ctx,
                               double xi, double yi, double zi, double si);
GAGE_EXPORT int gageStackProbeSpace(gageContext *ctx,
                                    double x, double y, double z, double s,
                                    int indexSpace, int clamp);

/* stackBlur.c */
GAGE_EXPORT gageStackBlurParm *gageStackBlurParmNew(void);
GAGE_EXPORT gageStackBlurParm *gageStackBlurParmNix(gageStackBlurParm *sbp);
GAGE_EXPORT int gageStackBlurParmScaleSet(gageStackBlurParm *sbp,
                                          unsigned int num,
                                          double scaleMin, double scaleMax,
                                          int uniform, int optim);
GAGE_EXPORT int gageStackBlurParmKernelSet(gageStackBlurParm *sbp,
                                           const NrrdKernelSpec *kspec,
                                           int renormalize);
GAGE_EXPORT int gageStackBlurParmSigmaMaxSet(gageStackBlurParm *sbp,
                                             double sigmaMax);
GAGE_EXPORT int gageStackBlurParmBoundarySet(gageStackBlurParm *sbp,
                                             int boundary, double padValue);
GAGE_EXPORT int gageStackBlurParmVerboseSet(gageStackBlurParm *sbp,
                                            int verbose);
GAGE_EXPORT int gageStackBlurParmOneDimSet(gageStackBlurParm *sbp,
                                           int oneDim);
GAGE_EXPORT int gageStackBlurParmCheck(gageStackBlurParm *sbp);
GAGE_EXPORT int gageStackBlur(Nrrd *const nblur[], gageStackBlurParm *sbp,
                              const Nrrd *nin, const gageKind *kind);
GAGE_EXPORT int gageStackBlurCheck(const Nrrd *const nblur[],
                                   gageStackBlurParm *sbp,
                                   const Nrrd *nin, const gageKind *kind);
GAGE_EXPORT int gageStackBlurGet(Nrrd *const nblur[], int *recomputedP,
                                 gageStackBlurParm *sbp,
                                 const char *format,
                                 const Nrrd *nin, const gageKind *kind);
GAGE_EXPORT int gageStackBlurManage(Nrrd ***nblurP, int *recomputedP,
                                    gageStackBlurParm *sbp,
                                    const char *format,
                                    int saveIfComputed, NrrdEncoding *enc,
                                    const Nrrd *nin, const gageKind *kind);

/* ctx.c */
GAGE_EXPORT gageContext *gageContextNew(void);
GAGE_EXPORT gageContext *gageContextCopy(gageContext *ctx);
GAGE_EXPORT gageContext *gageContextNix(gageContext *ctx);
GAGE_EXPORT void gageParmSet(gageContext *ctx, int which, double val);
GAGE_EXPORT int gagePerVolumeIsAttached(const gageContext *ctx,
                                        const gagePerVolume *pvl);
GAGE_EXPORT int gagePerVolumeAttach(gageContext *ctx, gagePerVolume *pvl);
GAGE_EXPORT int gagePerVolumeDetach(gageContext *ctx, gagePerVolume *pvl);
GAGE_EXPORT int gageKernelSet(gageContext *ctx, int which,
                              const NrrdKernel *k, const double *kparm);
GAGE_EXPORT void gageKernelReset(gageContext *ctx);
GAGE_EXPORT int gageProbe(gageContext *ctx, double xi, double yi, double zi);
GAGE_EXPORT int gageProbeSpace(gageContext *ctx, double x, double y, double z,
                               int indexSpace, int clamp);

/* update.c */
GAGE_EXPORT int gageUpdate(gageContext *ctx);

/* st.c */
GAGE_EXPORT int gageStructureTensor(Nrrd *nout, const Nrrd *nin,
                                    int dScale, int iScale, int dsmp);

/* deconvolve.c */
GAGE_EXPORT int gageDeconvolve(Nrrd *nout, double *lastDiffP,
                               const Nrrd *nin, const gageKind *kind,
                               const NrrdKernelSpec *ksp, int typeOut,
                               unsigned int maxIter, int saveAnyway,
                               double step, double epsilon, int verbose);
GAGE_EXPORT int gageDeconvolveSeparableKnown(const NrrdKernelSpec *ksp);
GAGE_EXPORT int gageDeconvolveSeparable(Nrrd *nout, const Nrrd *nin,
                                        const gageKind *kind,
                                        const NrrdKernelSpec *ksp,
                                        int typeOut);

#ifdef __cplusplus
}
#endif

#endif /* GAGE_HAS_BEEN_INCLUDED */