/usr/include/teem/gage.h is in libteem-dev 1.12.0~20160122-2.
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 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 | /*
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
/*
******** gage2Vec* enum
**
** "items" that gage knows how to measure in a 2-vector volume
*/
enum {
gage2VecUnknown, /* 0: nobody knows */
gage2VecVector, /* 1: "v", component-wise-interpolated
(CWI) vec: [2] */
gage2VecVector0, /* 2: "v0", vector[0]: [1] */
gage2VecVector1, /* 3: "v1", vector[0]: [1] */
gage2VecLength, /* 4: "l", length of 2-vector: [1] */
gage2VecLast
};
#define GAGE_2VEC_ITEM_MAX 4
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
/*
******** gageSigmaSampling* enum
**
** The different strategies for locating samples (pre-blurred images) as part
** of scale-space reconstruction. It is partly arbitrary, and partly in the
** interests of tangible familiarity, that of all the possible ways of
** describing scale, we are using sigma (hence the specificity in the name
** gageSigmaSampling)
**
** Adding this enum was motivated by the eventual plan to add a new sampling
** density parameter "rho".
*/
enum {
gageSigmaSamplingUnknown, /* 0 */
gageSigmaSamplingUniformSigma, /* 1 */
gageSigmaSamplingUniformTau, /* 2 */
gageSigmaSamplingOptimal3DL2L2, /* 3 */
gageSigmaSamplingLast
};
#define GAGE_SIGMA_SAMPLING_MAX 3
/*
******** 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"
** for the sake of scale-space analysis.
**
** This was reorganized August 2013 with the new attention to how
** scale-space is handled.
*/
typedef struct {
unsigned int num; /* # of pre-computed blurring scales == allocated
length of the "scale" vector below */
double sigmaRange[2]; /* sigma range for image blurring */
int sigmaSampling; /* from gageSigmaSampling* enum: how to sample the
range from sigmaRange[0] to sigmaRange[1] */
double *sigma; /* when-non-NULL, the sigma parameter for each
blurring level */
NrrdKernelSpec *kspec; /* the kernel with which we do blurring */
int renormalize; /* renormalize kernel weights (associated with the
kernel); passed to nrrdResampleRenormalizeSet */
NrrdBoundarySpec *bspec; /* what do to at image boundaries */
int oneDim, /* for experimental/debugging purposes: blur *only*
along the first (fastest) axis */
needSpatialBlur, /* always do blurring in the spatial domain, even
if frequency space blurring is possible */
verbose; /* verbosity level */
double dgGoodSigmaMax; /* The same info as communicated by
nrrdKernelDiscreteGaussianGoodSigmaMax, but
allowing it to be different. With this limit on
the sigma we pass to nrrdKernelDiscreteGaussian,
instead of doing the blurring in one step (when
doing spatial as opposed to frequency-space
blurring), the diffusion is done iteratively, with
steps in diffusion time of goodSigmaMax^2 */
} 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 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;
/* twovecGage.c */
GAGE_EXPORT const airEnum *const gage2Vec;
GAGE_EXPORT gageKind *const gageKind2Vec;
/* 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 */
/* moved the tau conversion functions to air, so that they can be
called from unu */
#define gageTauOfTee(TEE) airTauOfTime(TEE)
#define gageTeeOfTau(TAU) airTimeOfTau(TAU)
#define gageSigOfTau(TAU) airSigmaOfTau(TAU)
#define gageTauOfSig(SIG) airTauOfSigma(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 const airEnum *const gageSigmaSampling;
GAGE_EXPORT gageStackBlurParm *gageStackBlurParmNew(void);
GAGE_EXPORT int gageStackBlurParmCopy(gageStackBlurParm *sbpDst,
const gageStackBlurParm *sbpSrc);
GAGE_EXPORT void gageStackBlurParmInit(gageStackBlurParm *parm);
GAGE_EXPORT gageStackBlurParm *gageStackBlurParmNix(gageStackBlurParm *sbp);
GAGE_EXPORT int gageStackBlurParmCompare(const gageStackBlurParm *sbpA,
const char *nameA,
const gageStackBlurParm *sbpB,
const char *nameB,
int *differ,
char explain[AIR_STRLEN_LARGE]);
GAGE_EXPORT int gageStackBlurParmScaleSet(gageStackBlurParm *sbp,
unsigned int num,
double sigmaMin,
double sigmaMax,
int uniformSigma, int optimalL2L2);
GAGE_EXPORT int gageStackBlurParmSigmaSet(gageStackBlurParm *sbp,
unsigned int num,
double sigmaMin,
double sigmaMax,
int sigmaSampling);
GAGE_EXPORT int gageStackBlurParmKernelSet(gageStackBlurParm *sbp,
const NrrdKernelSpec *kspec);
GAGE_EXPORT int gageStackBlurParmRenormalizeSet(gageStackBlurParm *sbp,
int renormalize);
GAGE_EXPORT int gageStackBlurParmDgGoodSigmaMaxSet(gageStackBlurParm *sbp,
double dgGoodSigmaMax);
GAGE_EXPORT int gageStackBlurParmBoundarySet(gageStackBlurParm *sbp,
int boundary, double padValue);
GAGE_EXPORT int gageStackBlurParmBoundarySpecSet(gageStackBlurParm *sbp,
const NrrdBoundarySpec *bsp);
GAGE_EXPORT int gageStackBlurParmNeedSpatialBlurSet(gageStackBlurParm *sbp,
int sblur);
GAGE_EXPORT int gageStackBlurParmVerboseSet(gageStackBlurParm *sbp,
int verbose);
GAGE_EXPORT int gageStackBlurParmOneDimSet(gageStackBlurParm *sbp,
int oneDim);
GAGE_EXPORT int gageStackBlurParmCheck(const gageStackBlurParm *sbp);
GAGE_EXPORT int gageStackBlurParmParse(gageStackBlurParm *sbp,
int extraFlags[256],
char **extraParmsP,
const char *str);
GAGE_EXPORT hestCB *gageHestStackBlurParm;
GAGE_EXPORT int gageStackBlurParmSprint(char str[AIR_STRLEN_LARGE],
const gageStackBlurParm *sbp,
int extraFlag[256],
char *extraParm);
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 */
|