1 /* 2 Teem: Tools to process and visualize scientific data and images . 3 Copyright (C) 2012, 2011, 2010, 2009 University of Chicago 4 Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann 5 Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah 6 7 This library is free software; you can redistribute it and/or 8 modify it under the terms of the GNU Lesser General Public License 9 (LGPL) as published by the Free Software Foundation; either 10 version 2.1 of the License, or (at your option) any later version. 11 The terms of redistributing and/or modifying this software also 12 include exceptions to the LGPL that facilitate static linking. 13 14 This library is distributed in the hope that it will be useful, 15 but WITHOUT ANY WARRANTY; without even the implied warranty of 16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 17 Lesser General Public License for more details. 18 19 You should have received a copy of the GNU Lesser General Public License 20 along with this library; if not, write to Free Software Foundation, Inc., 21 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 22 */ 23 24 #ifndef GAGE_HAS_BEEN_INCLUDED 25 #define GAGE_HAS_BEEN_INCLUDED 26 27 #include <stdio.h> 28 #include <stdlib.h> 29 #include <limits.h> 30 #include <math.h> 31 32 #include <teem/air.h> 33 #include <teem/biff.h> 34 #include <teem/nrrd.h> 35 #include <teem/ell.h> 36 37 #if defined(_WIN32) && !defined(__CYGWIN__) && !defined(TEEM_STATIC) 38 # if defined(TEEM_BUILD) || defined(gage_EXPORTS) || defined(teem_EXPORTS) 39 # define GAGE_EXPORT extern __declspec(dllexport) 40 # else 41 # define GAGE_EXPORT extern __declspec(dllimport) 42 # endif 43 #else /* TEEM_STATIC || UNIX */ 44 # define GAGE_EXPORT extern 45 #endif 46 47 #ifdef __cplusplus 48 extern "C" { 49 #endif 50 51 #define GAGE gageBiffKey 52 53 /* 54 ** GAGE_DERIV_MAX is the maximum derivative that gage knows how to 55 ** work with. The 0th derivative is the reonstructed value (no 56 ** derivative). Many short arrays are allocated to 1 + this value 57 ** 58 ** MUST KEEP IN SYNC with GAGE_DV_* macros below 59 */ 60 #define GAGE_DERIV_MAX 2 61 62 /* 63 ** Convenience macros for arrays declared like blah[GAGE_DERIV_MAX+1] 64 ** (so these need to be kept in sync with GAGE_DERIV_MAX above). Note 65 ** that "DV" is short for "derivative vector". With consistent use of 66 ** these, changes to GAGE_DERIV_MAX will either obviously break code 67 ** (in an easy to fix way) or there will be no change needed at all. 68 ** 69 */ 70 #define GAGE_DV_SET(v, d0, d1, d2) \ 71 ((v)[0] = (d0), \ 72 (v)[1] = (d1), \ 73 (v)[2] = (d2)) 74 75 #define GAGE_DV_EQUAL(a, b) \ 76 ((a)[0] == (b)[0] && \ 77 (a)[1] == (b)[1] && \ 78 (a)[2] == (b)[2]) 79 80 #define GAGE_DV_COPY(a, b) \ 81 ((a)[0] = (b)[0], \ 82 (a)[1] = (b)[1], \ 83 (a)[2] = (b)[2]) 84 85 /* 86 ** the only extent to which gage treats different axes differently is 87 ** the spacing between samples along the axis. To have different 88 ** filters for the same function, but along different axes, would be 89 ** too messy. [Sun Mar 9 13:32:22 EDT 2008: Actually, doing per-axis 90 ** kernels is looking more and more sensible all the time. . .] Thus, 91 ** gage is not very useful as the engine for downsampling: it can't 92 ** tell that along one axis samples should be blurred while they 93 ** should be interpolated along another. Rather, it assumes that the 94 ** main task of probing is *reconstruction*: of values, of 95 ** derivatives, or lots of different quantities 96 */ 97 98 /* 99 ** terminology: gage is intended to measure multiple things at one 100 ** point in a volume. The SET of ALL the things that you want 101 ** gage to measure is called the "QUERY". Each of the many quantities 102 ** comprising the query are called "ITEM"s. 103 */ 104 105 /* 106 ******** gageParm.. enum 107 ** 108 ** these are passed to gageSet. Look for like-wise named field of 109 ** gageParm for documentation on what these mean. 110 ** 111 ** The following things have to agree: 112 ** - gageParm* enum 113 ** - fields of gageParm struct 114 ** - analogous gageDef* defaults (their declaration and setting) 115 ** - action of gageParmSet (ctx.c) 116 ** - action of gageParmReset (miscGage.c) 117 */ 118 enum { 119 gageParmUnknown, 120 gageParmVerbose, /* non-boolean int */ 121 gageParmRenormalize, /* int */ 122 gageParmCheckIntegrals, /* int */ 123 gageParmK3Pack, /* int */ 124 gageParmGradMagCurvMin, /* double */ 125 gageParmCurvNormalSide, /* int */ 126 gageParmKernelIntegralNearZero, /* double */ 127 gageParmDefaultCenter, /* int */ 128 gageParmStackUse, /* int */ 129 gageParmStackNormalizeDeriv, /* int; does NOT imply gageParmStackUse */ 130 gageParmStackNormalizeDerivBias, /* double; does NOT imply " */ 131 gageParmStackNormalizeRecon, /* int; does NOT imply " */ 132 gageParmOrientationFromSpacing, /* int */ 133 gageParmGenerateErrStr, /* int */ 134 gageParmLast 135 }; 136 137 /* keep in synch with gageErr airEnum */ 138 enum { 139 gageErrUnknown, /* 0: nobody knows */ 140 gageErrNone, /* 1: no error, actually, all's well */ 141 gageErrBoundsSpace, /* 2: out of 3-D (index-space) bounds */ 142 gageErrBoundsStack, /* 3: out of 1-D bounds of stack */ 143 gageErrStackIntegral, /* 4: stack recon coeff's sum to 0, so can't 144 normalize them */ 145 gageErrStackSearch, /* 5: for some reason couldn't find the index 146 position of the probed stack location */ 147 gageErrStackUnused, /* 6: can't probe stack without parm.stackUse */ 148 gageErrLast 149 }; 150 #define GAGE_ERR_MAX 6 151 152 /* 153 ******** gage{Ctx,Pvl}Flag.. enum 154 ** 155 ** organizes all the dependendies within a context and between a 156 ** context and pervolumes. This logic is to determine the support 157 ** required for a query: different queries need different kernels, 158 ** different kernels have different supports. The user should not 159 ** have to be concerned about any of this; it should be useful only 160 ** to gageUpdate(). 161 */ 162 enum { 163 gageCtxFlagUnknown, 164 gageCtxFlagNeedD, /* 1: derivatives required for query */ 165 gageCtxFlagK3Pack, /* 2: whether to use 3 or 6 kernels */ 166 gageCtxFlagNeedK, /* 3: which of the kernels will actually be used */ 167 gageCtxFlagKernel, /* 4: any one of the kernels or its parameters */ 168 gageCtxFlagRadius, /* 5: radius of support for kernels with query */ 169 gageCtxFlagShape, /* 6: a new pervolume shape was set */ 170 gageCtxFlagLast 171 }; 172 #define GAGE_CTX_FLAG_MAX 6 173 174 enum { 175 gagePvlFlagUnknown, 176 gagePvlFlagVolume, /* 1: got a new volume */ 177 gagePvlFlagQuery, /* 2: what do you really care about */ 178 gagePvlFlagNeedD, /* 3: derivatives required for query */ 179 gagePvlFlagLast 180 }; 181 #define GAGE_PVL_FLAG_MAX 3 182 183 184 /* 185 ******** gageKernel.. enum 186 ** 187 ** these are the different kernels that might be used in gage, regardless 188 ** of what kind of volume is being probed. 189 ** (synch with miscGage.c:gageKernel airEnum) 190 */ 191 enum { 192 gageKernelUnknown, /* 0: nobody knows */ 193 gageKernel00, /* 1: measuring values */ 194 gageKernel10, /* 2: reconstructing 1st derivatives */ 195 gageKernel11, /* 3: measuring 1st derivatives */ 196 gageKernel20, /* 4: reconstructing 1st partials and 2nd deriv.s */ 197 gageKernel21, /* 5: measuring 1st partials for a 2nd derivative */ 198 gageKernel22, /* 6: measuring 2nd derivatives */ 199 gageKernelStack, /* 7: for reconstructing across a stack */ 200 gageKernelLast 201 }; 202 #define GAGE_KERNEL_MAX 7 203 204 /* 205 ******** GAGE_ITEM_PREREQ_MAXNUM 206 ** 207 ** Max number of prerequisites for any item in *any* kind. 208 ** 209 ** This value has gotten bumped periodically, which used to mean 210 ** that *all* item tables had to be updated, when "-1" was used 211 ** represent the unknown item. But now that 0 represents the 212 ** unknown item, and because struct elements are implicitly 213 ** initialized to zero, this is no longer the case. 214 ** 215 ** Wed Nov 8 14:12:44 EST 2006 pre-emptively upping this from 6 216 */ 217 #define GAGE_ITEM_PREREQ_MAXNUM 8 218 219 /* 220 ******** gageItemEntry struct 221 ** 222 ** necessary information about each item supported by the kind, which 223 ** is defined at compile-time in a per-kind table (at least it is for 224 ** the scalar, vector, and tensor kinds) 225 ** 226 ** NOTE!!! YOU CAN NOT re-arrange these variables, because of all the 227 ** compile-time definitions that are done to define the gageKind->table 228 ** for all current kinds. 229 */ 230 typedef struct { 231 int enumVal; /* the item's enum value */ 232 unsigned int 233 answerLength; /* how many double's are needed to store the answer, 234 or (for non-zero items), 0 to represent 235 "the length will be learned later at runtime" */ 236 int needDeriv, /* what kind of derivative info is immediately needed 237 for this item (not recursively expanded). This is 238 NO LONGER a bitvector: values are 0, 1, 2, .. */ 239 prereq[GAGE_ITEM_PREREQ_MAXNUM], 240 /* what are the other items this item depends on 241 (not recusively expanded), you can list up to 242 GAGE_ITEM_PREREQ_MAXNUM of them, and use 0 243 (the unknown item) to fill out the list. 244 _OR_ -1 if this is a dynamic kind and the prereqs 245 will not be known until later in runtime */ 246 parentItem, /* the enum value of an item (answerLength > 1) 247 containing the smaller value for which we are 248 merely an alias 249 _OR_ 0 if there's no parent */ 250 parentIndex, /* where our value starts in parents value 251 _OR_ 0 if there's no parent */ 252 needData; /* whether non-NULL gagePerVolume->data is needed 253 for answering this item */ 254 } gageItemEntry; 255 256 /* 257 ** modifying the enums below (scalar, vector, etc query quantities) 258 ** necesitates modifying: 259 ** - the central item table in scl.c, vecGage.c 260 ** - the associated arrays in scl.c, vecGage.c 261 ** - the "answer" method itself in sclanswer.c, vecGage.c 262 */ 263 264 /* 265 ******** gageScl* enum 266 ** 267 ** all the "items" that gage can measure in a scalar volume. 268 ** 269 ** NOTE: although it is currently listed that way, it is not necessary 270 ** that prerequisite measurements are listed before the other measurements 271 ** which need them (that is represented by _gageSclPrereq) 272 ** 273 ** The description for each enum value starts with the numerical value 274 ** followed by a string which identifies the value in the gageScl airEnum. 275 ** The "[N]" indicates how many doubles are used for storing the quantity. 276 */ 277 enum { 278 gageSclUnknown, /* 0: nobody knows */ 279 gageSclValue, /* 1: "v", data value: [1] */ 280 gageSclGradVec, /* 2: "grad", gradient vector, un-normalized: [3] */ 281 gageSclGradMag, /* 3: "gm", gradient magnitude: [1] */ 282 gageSclNormal, /* 4: "n", gradient vector, normalized: [3] */ 283 gageSclNProj, /* 5: "nproj", projection onto normal: [9] */ 284 gageSclNPerp, /* 6: "np", projection onto tangent plane: [9] */ 285 gageSclHessian, /* 7: "h", Hessian: [9] (column-order) */ 286 gageSclHessianTen, /* 8: "ht", Hessian as 7-component tensor: [7] 287 In principle this is a dependency inversion 288 since gage doesn't know about the ten library, 289 where the 7-element tensor is based. */ 290 gageSclLaplacian, /* 9: "l", Laplacian: Dxx + Dyy + Dzz: [1] */ 291 gageSclHessFrob, /* 10: "hf", Frobenius norm of Hessian: [1] */ 292 gageSclHessEval, /* 11: "heval", Hessian's eigenvalues: [3] */ 293 gageSclHessEval0, /* 12: "heval0", Hessian's 1st eigenvalue: [1] */ 294 gageSclHessEval1, /* 13: "heval1", Hessian's 2nd eigenvalue: [1] */ 295 gageSclHessEval2, /* 14: "heval2", Hessian's 3rd eigenvalue: [1] */ 296 gageSclHessEvec, /* 15: "hevec", Hessian's eigenvectors: [9] */ 297 gageSclHessEvec0, /* 16: "hevec0", Hessian's 1st eigenvector: [3] */ 298 gageSclHessEvec1, /* 17: "hevec1", Hessian's 2nd eigenvector: [3] */ 299 gageSclHessEvec2, /* 18: "hevec2", Hessian's 3rd eigenvector: [3] */ 300 gageScl2ndDD, /* 19: "2d", 2nd dir.deriv. along gradient: [1] */ 301 gageSclGeomTens, /* 20: "gten", sym. matx w/ evals {0, K1, K2} and 302 evecs {grad, cdir0, cdir1}: [9] */ 303 gageSclGeomTensTen, /* 21: "gtenten", 7-element geometry tensor [7] */ 304 gageSclK1, /* 22: "k1", 1st principle curvature: [1] */ 305 gageSclK2, /* 23: "k2", 2nd principle curvature (k2 <= k1): [1] */ 306 gageSclTotalCurv, /* 24: "tc", L2 norm(K1,K2) (not Koen.'s "C"): [1] */ 307 gageSclShapeTrace, /* 25, "st", (K1+K2)/Curvedness: [1] */ 308 gageSclShapeIndex, /* 26: "si", Koen.'s shape index, ("S"): [1] */ 309 gageSclMeanCurv, /* 27: "mc", mean curvature (K1 + K2)/2: [1] */ 310 gageSclGaussCurv, /* 28: "gc", gaussian curvature K1*K2: [1] */ 311 gageSclCurvDir1, /* 29: "cdir1", 1st principle curv direction: [3] */ 312 gageSclCurvDir2, /* 30: "cdir2", 2nd principle curv direction: [3] */ 313 gageSclFlowlineCurv, /* 31: "fc", curvature of normal streamline: [1] */ 314 gageSclMedian, /* 32: "med", median filter */ 315 gageSclHessValleyness, /* 33: "hvalley", vallyness measure: [1] */ 316 gageSclHessRidgeness, /* 34: "hridge", ridgeness measure: [1] */ 317 gageSclHessMode, /* 35: "hmode", Hessian's mode: [1] */ 318 gageSclLast 319 }; 320 #define GAGE_SCL_ITEM_MAX 35 321 322 /* 323 ******** gageVec* enum 324 ** 325 ** all the "items" that gage knows how to measure in a 3-vector volume 326 ** 327 ** The strings gives one of the gageVec airEnum identifiers, and [x] 328 ** says how many scalars are associated with this value. 329 */ 330 enum { 331 gageVecUnknown, /* 0: nobody knows */ 332 gageVecVector, /* 1: "v", component-wise-interpolated 333 (CWI) vec: [3] */ 334 gageVecVector0, /* 2: "v0", vector[0]: [1] */ 335 gageVecVector1, /* 3: "v1", vector[0]: [1] */ 336 gageVecVector2, /* 4: "v2", vector[0]: [1] */ 337 gageVecLength, /* 5: "l", length of CWI vector: [1] */ 338 gageVecNormalized, /* 6: "n", normalized CWI vector: [3] */ 339 gageVecJacobian, /* 7: "j", component-wise Jacobian: [9] 340 0:dv_x/dx 1:dv_x/dy 2:dv_x/dz 341 3:dv_y/dx 4:dv_y/dy 5:dv_y/dz 342 6:dv_z/dx 7:dv_z/dy 8:dv_z/dz */ 343 gageVecStrain, /* 8: "S", rate-of-strain tensor: [9] */ 344 gageVecDivergence, /* 9: "d", divergence (based on Jacobian): [1] */ 345 gageVecCurl, /* 10: "c", curl (based on Jacobian): [3] */ 346 gageVecCurlNorm, /* 11: "cm", curl magnitude: [1] */ 347 gageVecHelicity, /* 12: "h", helicity: vec . curl: [1] */ 348 gageVecNormHelicity, /* 13: "nh", normalized helicity: [1] */ 349 gageVecSOmega, /* 14: "somega", S squared + Omega squared: [9] */ 350 gageVecLambda2, /* 15: "lambda2", lambda2 criterion: [1] */ 351 gageVecImaginaryPart, /* 16: "imag", imag. part of jacobian's 352 complex-conjugate eigenvalues: [1] */ 353 gageVecHessian, /* 17: "vh", second-order derivative: [27] 354 HEY: indices here need to be double checked 355 0:d2v_x/dxdx 1:d2v_x/dxdy 2:d2v_x/dxdz 356 3:d2v_x/dydx 4:d2v_x/dydy 5:d2v_x/dydz 357 6:d2v_x/dzdx 7:d2v_x/dzdy 8:d2v_x/dzdz 358 9:d2v_y/dxdx [..] 359 [..] 360 24:dv2_z/dzdx 25:d2v_z/dzdy 26:d2v_z/dzdz */ 361 gageVecDivGradient, /* 18: "dg", divergence gradient: [3] */ 362 gageVecCurlGradient, /* 19: "cg", curl gradient: [9] */ 363 gageVecCurlNormGrad, /* 20: "cng", curl norm gradient: [3] */ 364 gageVecNCurlNormGrad, /* 21: "ncng", normalized curl norm gradient: [3] */ 365 gageVecHelGradient, /* 22: "hg", helicity gradient: [3] */ 366 gageVecDirHelDeriv, /* 23: "dhd", directional derivative 367 of helicity: [1] */ 368 gageVecProjHelGradient, /* 24: "phg", projected helicity gradient: [3] */ 369 gageVecGradient0, /* 25: "g0", gradient of 1st coeff of vector: [3] */ 370 gageVecGradient1, /* 26: "g1", gradient of 2nd coeff of vector: [3] */ 371 gageVecGradient2, /* 27: "g2", gradient of 3rd coeff of vector: [3] */ 372 gageVecMultiGrad, /* 28: "mg", sum of outer products of grads: [9] */ 373 gageVecMGFrob, /* 29: "mgfrob", frob norm of multi-gradient: [1] */ 374 gageVecMGEval, /* 30: "mgeval", evals of multi-gradient: [3] */ 375 gageVecMGEvec, /* 31: "mgevec", evecs of multi-gradient: [9] */ 376 gageVecLast 377 }; 378 #define GAGE_VEC_ITEM_MAX 31 379 380 struct gageKind_t; /* dumb forward declaraction, ignore */ 381 struct gagePerVolume_t; /* dumb forward declaraction, ignore */ 382 383 /* 384 ******** gageItemPackPart* enum 385 ** 386 ** the different ways that some items can be related for a gagePack 387 */ 388 enum { 389 gageItemPackPartUnknown, /* 0 */ 390 gageItemPackPartScalar, /* 1 */ 391 gageItemPackPartGradVec, /* 2 */ 392 gageItemPackPartGradMag, /* 3 */ 393 gageItemPackPartNormal, /* 4 */ 394 gageItemPackPartHessian, /* 5 */ 395 gageItemPackPartHessEval0, /* 6 */ 396 gageItemPackPartHessEval1, /* 7 */ 397 gageItemPackPartHessEval2, /* 8 */ 398 gageItemPackPartHessEvec0, /* 9 */ 399 gageItemPackPartHessEvec1, /* 10 */ 400 gageItemPackPartHessEvec2, /* 11 */ 401 gageItemPackPartLast 402 }; 403 #define GAGE_ITEM_PACK_PART_MAX 11 404 405 /* 406 ******** gageShape struct 407 ** 408 ** just a container for all the information related to the "shape" 409 ** of all the volumes associated with a context 410 ** 411 ** Note that the utility of gageShape has extended well beyond doing 412 ** convolution-based measurements in volumes- it has become the 413 ** one-stop place for all of Teem to figure out a reasonable way of 414 ** locating a logically a volume in 3-D space, including using a 415 ** nrrd's full orientation information if it is known. 416 */ 417 typedef struct gageShape_t { 418 /* ========= INPUT ========= (controls for _gageShapeSet) */ 419 int defaultCenter, /* default centering to use when given volume 420 has no centering set. */ 421 orientationFromSpacing; /* only meaningful if nrrd has per-axis spacing, 422 but not full orientation info. If zero, the 423 volume is crammed into the bi-unit cube. 424 If non-zero, gage treats the volume as if it 425 had axis-aligned spaceDirection vectors, with 426 the non-zero values determined by the given 427 per-axis spacing. */ 428 /* ======== OUTPUT ========= (set by _gageShapeSet) */ 429 int center, /* the sample centering of the volume(s)- this 430 determines the extent of the locations 431 that may be probed */ 432 fromOrientation; /* non-zero iff the spaceDirections and 433 spaceOrigin information was used */ 434 unsigned int size[3]; /* raster dimensions of volume 435 HEY: shouldn't this be size_t? */ 436 double spacing[3]; /* spacings for each axis */ 437 /* fwScale[GAGE_KERNEL_MAX+1][3] has been superceded by the cleaner and 438 more general ItoWSubInvTransp and ItoWSubInv matrices below */ 439 double ItoW[16], /* homogeneous coord transform from index 440 to world space (defined either by bi-unit 441 cube or from full orientation info ) */ 442 WtoI[16], /* inverse of above */ 443 ItoWSubInvTransp[9], /* inverse transpose of 3x3 sub-matrix of ItoW, 444 to transform (covariant) gradients */ 445 ItoWSubInv[9]; /* tranpose of above, to transform hessians */ 446 447 } gageShape; 448 449 /* 450 ******** gageParm struct 451 ** 452 ** a container for the various switches and knobs which control 453 ** gage, aside from the obvious inputs (kernels, queries, volumes) 454 */ 455 typedef struct gageParm_t { 456 int renormalize; /* hack to make sure that sum of 457 discrete value reconstruction weights 458 is same as kernel's continuous 459 integral, and that the 1nd and 2nd 460 deriv weights really sum to 0.0 */ 461 int checkIntegrals; /* call the "integral" method of the 462 kernel to verify that it is 463 appropriate for the task for which 464 the kernel is being set: 465 reconstruction: 1.0, derivatives: 0.0 */ 466 int k3pack; /* non-zero (true) iff we do not use kernels 467 gageKernelIJ with I != J. So, we use the 468 value reconstruction kernel (gageKernel00) 469 for 1st and 2nd derivative reconstruction, 470 and so on. This is faster because we can 471 re-use results from low-order convolutions. 472 The name "3pack" will likely persist even 473 with 3rd dervatives, for which "4pack" would 474 make more sense (for 00, 11, 22, 33) */ 475 double gradMagCurvMin, /* punt on computing curvature information if 476 gradient magnitude is less than this. Yes, 477 this is scalar-kind-specific, but there's 478 no other good place for it */ 479 kernelIntegralNearZero, /* tolerance with checkIntegrals on derivative 480 kernels */ 481 stackNormalizeDerivBias; /* when using stackNormalizeDeriv, a bias 482 on the effective scale, used for the 483 normalization */ 484 int curvNormalSide, /* determines direction of gradient that is used 485 as normal in curvature calculations, exactly 486 the same as miteUser's normalSide: 1 for 487 normal pointing to lower values (higher 488 values are more "inside"); -1 for normal 489 pointing to higher values (low values more 490 "inside") */ 491 defaultCenter, /* what centering to use when you have to invent 492 one, because its not set on any axis */ 493 stackUse, /* if non-zero: treat the pvl's (all same kind) 494 as multiple values of a single logical volume 495 (e.g. for approximating scale space). 496 The first pvl is effectively just a buffer; 497 the N-1 pvls used are at index 1 through N-2. 498 The query in pvl[0] will determine the 499 computations done, and answers for the whole 500 stack will be stored in pvl[0]. */ 501 stackNormalizeRecon, /* if non-zero (and if stackUse is non-zero): 502 the weights used to reconstruct across 503 stack samples are normalized to unit sum; not 504 needed if the kernel is accurate enough */ 505 stackNormalizeDeriv, /* if non-zero (and if stackUse is non-zero): 506 derivatives at filter stage (before answer 507 stage) are renormalized based on the stack 508 position */ 509 orientationFromSpacing, /* only meaningful if nrrd has per-axis spacing, 510 but not full orientation info. If zero, the 511 volume is crammed into the bi-unit cube. 512 If non-zero, gage treats the volume as if it 513 had axis-aligned spaceDirection vectors, with 514 the non-zero values determined by the given 515 per-axis spacing. */ 516 generateErrStr; /* when errors happen, biff is never used, but 517 a descriptive error is sprintf into 518 gctx->errStr as long as this is non-zero. */ 519 } gageParm; 520 521 /* 522 ******** gagePoint struct 523 ** 524 ** stores *Index Space* location of last query location, which is used 525 ** to determine whether the ctx->fsl, ctx->fw values can be re-used 526 ** (based on the "Frac" values), and, whether all the pvl->iv3 have to 527 ** be refilled (based on the "Idx" values). The last index (frac[3] 528 ** and idx[3])is for the stack, and can safely stay 0 if the stack 529 ** isn't being used. 530 ** 531 ** with stack usage, stackFwNonZeroNum records how many pvls had 532 ** non-zero stack filter weights, which is used to detect when 533 ** iv3s have to be refilled. Looking at idx[3] alone is not always 534 ** sufficient for this. 535 */ 536 typedef struct gagePoint_t { 537 double frac[4]; /* last fractional voxel location */ 538 unsigned int idx[4], /* last integral voxel location; actually the 539 *upper* corner of the containing voxel */ 540 stackFwNonZeroNum; /* last number of non-zero values of stack filter 541 weights (ctx->stackFw) */ 542 } gagePoint; 543 544 /* 545 ******** gageQuery typedef 546 ** 547 ** this short, fixed-length array is used as a bit-vector to store 548 ** all the items that comprise a query. Its length sets an upper 549 ** bound on the maximum item value that a gageKind may use. 550 ** 551 ** The important thing to keep in mind is that a variable of type 552 ** gageKind ends up being passed by reference, even though the 553 ** syntax of its usage doesn't immediately announce that. 554 ** 555 ** gageKindCheck currently has the role of verifying that a gageKind's 556 ** maximum item value is within the bounds set here. Using 557 ** GAGE_QUERY_BYTES_NUM == 8 gives a max item value of 63, which is 558 ** far above anything being used now. 559 ** 560 ** Sat Jan 21 18:12:01 EST 2006: ha! second derivatives of tensors blew 561 ** past old GAGE_QUERY_BYTES_NUM, now GAGE_QUERY_BYTES_NUM == 16 562 ** 563 ** Tue Nov 7 19:51:05 EST 2006; tenGage items now pushing 127, 564 ** guardedly changing GAGE_QUERY_BYTES_NUM to 24 --> max item 191 565 ** 566 ** Wed Nov 18 17:53:23 CST 2009; yikes, tenGage items now at 190, 567 ** changing GAGE_QUERY_BYTES_NUM to 32 --> max item 255 568 ** 569 ** NOTE: increasing GAGE_QUERY_BYTES_NUM means that the macros below 570 ** have to be redefined as well! 571 */ 572 #define GAGE_QUERY_BYTES_NUM 32 573 #define GAGE_ITEM_MAX ((8*GAGE_QUERY_BYTES_NUM)-1) 574 typedef unsigned char gageQuery[GAGE_QUERY_BYTES_NUM]; 575 576 /* 577 ******** GAGE_QUERY_RESET, GAGE_QUERY_TEST 578 ******** GAGE_QUERY_ON, GAGE_QUERY_OFF 579 ** 580 ** Macros for manipulating a gageQuery. 581 ** 582 ** airSanity ensures that an unsigned char is in fact 8 bits 583 */ 584 #define GAGE_QUERY_RESET(q) \ 585 q[ 0] = q[ 1] = q[ 2] = q[ 3] = \ 586 q[ 4] = q[ 5] = q[ 6] = q[ 7] = \ 587 q[ 8] = q[ 9] = q[10] = q[11] = \ 588 q[12] = q[13] = q[14] = q[15] = \ 589 q[16] = q[17] = q[18] = q[19] = \ 590 q[20] = q[21] = q[22] = q[23] = \ 591 q[24] = q[25] = q[26] = q[27] = \ 592 q[28] = q[29] = q[30] = q[31] = 0 593 594 #define GAGE_QUERY_COPY(p, q) \ 595 p[ 0] = q[ 0]; p[ 1] = q[ 1]; p[ 2] = q[ 2]; p[ 3] = q[ 3]; \ 596 p[ 4] = q[ 4]; p[ 5] = q[ 5]; p[ 6] = q[ 6]; p[ 7] = q[ 7]; \ 597 p[ 8] = q[ 8]; p[ 9] = q[ 9]; p[10] = q[10]; p[11] = q[11]; \ 598 p[12] = q[12]; p[13] = q[13]; p[14] = q[14]; p[15] = q[15]; \ 599 p[16] = q[16]; p[17] = q[17]; p[18] = q[18]; p[19] = q[19]; \ 600 p[20] = q[20]; p[21] = q[21]; p[22] = q[22]; p[23] = q[23]; \ 601 p[24] = q[24]; p[25] = q[25]; p[26] = q[26]; p[27] = q[27]; \ 602 p[28] = q[28]; p[29] = q[29]; p[30] = q[30]; p[31] = q[31] 603 604 #define GAGE_QUERY_ADD(p, q) \ 605 p[ 0] |= q[ 0]; p[ 1] |= q[ 1]; p[ 2] |= q[ 2]; p[ 3] |= q[ 3]; \ 606 p[ 4] |= q[ 4]; p[ 5] |= q[ 5]; p[ 6] |= q[ 6]; p[ 7] |= q[ 7]; \ 607 p[ 8] |= q[ 8]; p[ 9] |= q[ 9]; p[10] |= q[10]; p[11] |= q[11]; \ 608 p[12] |= q[12]; p[13] |= q[13]; p[14] |= q[14]; p[15] |= q[15]; \ 609 p[16] |= q[16]; p[17] |= q[17]; p[18] |= q[18]; p[19] |= q[19]; \ 610 p[20] |= q[20]; p[21] |= q[21]; p[22] |= q[22]; p[23] |= q[23]; \ 611 p[24] |= q[24]; p[25] |= q[25]; p[26] |= q[26]; p[27] |= q[27]; \ 612 p[28] |= q[28]; p[29] |= q[29]; p[30] |= q[30]; p[31] |= q[31] 613 614 #define GAGE_QUERY_EQUAL(p, q) ( \ 615 p[ 0] == q[ 0] && p[ 1] == q[ 1] && p[ 2] == q[ 2] && p[ 3] == q[ 3] && \ 616 p[ 4] == q[ 4] && p[ 5] == q[ 5] && p[ 6] == q[ 6] && p[ 7] == q[ 7] && \ 617 p[ 8] == q[ 8] && p[ 9] == q[ 9] && p[10] == q[10] && p[11] == q[11] && \ 618 p[12] == q[12] && p[13] == q[13] && p[14] == q[14] && p[15] == q[15] && \ 619 p[16] == q[16] && p[17] == q[17] && p[18] == q[18] && p[19] == q[19] && \ 620 p[20] == q[20] && p[21] == q[21] && p[22] == q[22] && p[23] == q[23] && \ 621 p[24] == q[24] && p[25] == q[25] && p[26] == q[26] && p[27] == q[27] && \ 622 p[28] == q[28] && p[29] == q[29] && p[30] == q[30] && p[31] == q[31] ) 623 624 #define GAGE_QUERY_NONZERO(q) ( \ 625 q[ 0] | q[ 1] | q[ 2] | q[ 3] | \ 626 q[ 4] | q[ 5] | q[ 6] | q[ 7] | \ 627 q[ 8] | q[ 9] | q[10] | q[11] | \ 628 q[12] | q[13] | q[14] | q[15] | \ 629 q[16] | q[17] | q[18] | q[19] | \ 630 q[20] | q[21] | q[22] | q[23] | \ 631 q[24] | q[25] | q[26] | q[27] | \ 632 q[28] | q[29] | q[30] | q[31] ) 633 634 #define GAGE_QUERY_ITEM_TEST(q, i) (q[i/8] & (1 << (i % 8))) 635 #define GAGE_QUERY_ITEM_ON(q, i) (q[i/8] |= (1 << (i % 8))) 636 #define GAGE_QUERY_ITEM_OFF(q, i) (q[i/8] &= ~(1 << (i % 8))) 637 638 /* increment for ctx->pvlArr airArray */ 639 #define GAGE_PERVOLUME_ARR_INCR 32 640 641 /* extents of known information about optimal sigma samples for 642 Hermite-spline-based scale-space reconstruction */ 643 #define GAGE_OPTIMSIG_SIGMA_MAX 11 644 #define GAGE_OPTIMSIG_SAMPLES_MAXNUM 11 645 646 /* 647 ******** gageContext struct 648 ** 649 ** The information here is specific to the dimensions, scalings, and 650 ** radius of kernel support, but not to kind of volume (all kind-specific 651 ** information is in the gagePerVolume). One context can be used in 652 ** conjuction with probing multiple volumes. 653 */ 654 typedef struct gageContext_t { 655 /* INPUT ------------------------- */ 656 int verbose; /* verbosity */ 657 gageParm parm; /* all parameters */ 658 659 /* all the kernels we'll ever need, including the stack kernel */ 660 NrrdKernelSpec *ksp[GAGE_KERNEL_MAX+1]; 661 662 /* all the pervolumes attached to this context. If using stack, 663 the base pvl is the LAST, pvl[pvlNum-1], and the stack samples 664 are pvl[0] through pvl[pvlNum-2] */ 665 struct gagePerVolume_t **pvl; 666 667 /* number of pervolumes currently attached. If using stack, 668 this is one more than number of stack samples (because of the 669 base volume at the end) */ 670 unsigned int pvlNum; 671 672 /* airArray for managing pvl and pvlNum */ 673 airArray *pvlArr; 674 675 /* sizes, spacings, centering, and other geometric aspects of the 676 volume */ 677 gageShape *shape; 678 679 /* if stack is being used, allocated for length pvlNum-1, and 680 stackPos[0] through stackPos[pvlNum-2] MUST exist and be 681 monotonically increasing stack positions for each volume. 682 Otherwise NULL */ 683 double *stackPos; 684 685 /* INTERNAL ------------------------- */ 686 /* if using stack: allocated for length pvlNum-1, and filter sample 687 locations and weights for reconstruction along the stack. 688 Otherwise NULL. */ 689 double *stackFsl; 690 double *stackFw; 691 692 /* all the flags used by gageUpdate() used to describe what changed 693 in this context */ 694 int flag[GAGE_CTX_FLAG_MAX+1]; 695 696 /* which value/derivatives need to be calculated for all pervolumes 697 (doV, doD1, doD2) */ 698 int needD[GAGE_DERIV_MAX+1]; 699 700 /* which kernels are needed for all pvls. needK[gageKernelStack] 701 is currently not set by the update function that sets needK[] */ 702 int needK[GAGE_KERNEL_MAX+1]; 703 704 /* radius of support of samples needed to satisfy query, given the 705 set of kernels. The "filter diameter" fd == 2*radius. This is 706 incremented by one when filtering across the stack with 707 nrrdKernelHermiteScaleSpaceFlag */ 708 unsigned int radius; 709 710 /* filter sample locations (all axes): logically a fd x 3 array 711 (and its 3 because gage works on 3D fields, not because of 712 the number of derivatives supported) */ 713 double *fsl; 714 715 /* filter weights (all axes, all kernels): logically a 716 fd x 3 x GAGE_KERNEL_MAX+1 (fast-to-slow) array */ 717 double *fw; 718 719 /* offsets to other fd^3 samples needed to fill 3D intermediate 720 value cache. Allocated size is dependent on kernels, values 721 inside are dependent on the dimensions of the volume. It may be 722 more correct to be using size_t instead of uint, but the X and Y 723 dimensions of the volume would have to be super-outrageous for 724 that to be a problem */ 725 unsigned int *off; 726 727 /* last probe location */ 728 gagePoint point; 729 730 /* errStr and errNum are for describing errors that happen in gageProbe(): 731 using biff is too heavy-weight for this, and the idea is that no ill 732 should occur if the error is repeatedly ignored. Whether or not 733 something is actually sprintf'ed into errStr is controlled by 734 parm.generateErrStr. 735 NOTE: these variables used to be globals "gageErrStr" and "gageErrNum" */ 736 char errStr[AIR_STRLEN_LARGE]; 737 int errNum; /* takes values from the gageErr enum */ 738 739 /* what fraction of the values in the kernel support had to be invented 740 (by bleeding out the margin) in order to satisfy a probe that was near 741 the edge (any axis, either high or low) of the volume */ 742 double edgeFrac; 743 } gageContext; 744 745 /* 746 ******** gagePerVolume 747 ** 748 ** information that is specific to one volume, and to one kind of 749 ** volume. 750 */ 751 typedef struct gagePerVolume_t { 752 int verbose; /* verbosity */ 753 const struct gageKind_t *kind; /* what kind of volume is this for */ 754 gageQuery query; /* the query, recursively expanded */ 755 int needD[GAGE_DERIV_MAX+1];/* which derivatives need to be calculated for 756 the query (above) in this volume */ 757 const Nrrd *nin; /* the volume to sample within */ 758 int flag[GAGE_PVL_FLAG_MAX+1];/* for the kind-specific flags .. */ 759 double *iv3, *iv2, *iv1; /* 3D, 2D, 1D, value caches. These are cubical, 760 square, and linear arrays, all length fd on 761 each edge. Currently gageIv3Fill() fills 762 the iv3 (at a latter point this will be 763 delegated back to the gageKind to facilitate 764 bricking), and currently the tuple axis (with 765 length valLen) always slowest. However, use 766 of iv2 and iv1 is entirely up the kind's 767 filter method. */ 768 double (*lup)(const void *ptr, size_t I); 769 /* nrrd{F,D}Lookup[] element, according to 770 nin->type and double */ 771 double *answer; /* main buffer to hold all the answers */ 772 double **directAnswer; /* array of pointers into answer */ 773 void *data; /* extra data, parameters, buffers, etc. 774 required for answering some items (as per 775 the gageItemEntry->needData) managed with 776 kind->pvlDataNew, kind->pvlDataCopy, 777 kind->pvlDataUpdate, and kind->pvlDataNix, 778 so there is no channel for extra info to be 779 passed into the pvl->data, other that what 780 was put into kind->data */ 781 } gagePerVolume; 782 783 /* 784 ******** gageKind struct 785 ** 786 ** all the information and functions that are needed to handle one 787 ** kind of volume (such as scalar, vector, etc.) 788 ** 789 ** these are either statically allocated (e.g. gageKindScl, gageKindVec, 790 ** tenGageKind), or dynamically allocated, which case the kind itself 791 ** needs a constructor (e.g. tenDwiGageKindNew()). The "dynamicAlloc" 792 ** variable indicates this distinction. 793 ** 794 ** Having dynamically allocated kinds raises the possibility of having 795 ** to set and update (or modify and update) their internal state, 796 ** which is currently completely outside the update framework of gage. 797 ** So as far as the core gage functions are concerned, all kinds are 798 ** static, because there is nothing to modify. It also means that 799 ** those who dynamically create kinds should try to minimize the 800 ** state info that can/must be dynamically modified (i.e. maybe 801 ** the kind constructor should take all the various parameters, 802 ** and set everything in a single shot). 803 */ 804 typedef struct gageKind_t { 805 int dynamicAlloc; /* non-zero if this kind struct was 806 dynamically allocated */ 807 char name[AIR_STRLEN_SMALL]; /* short identifying string for kind */ 808 const airEnum *enm; /* such as gageScl. NB: the "unknown" 809 value in the enum *must* be 0. */ 810 unsigned int baseDim, /* dimension that x,y,z axes start on 811 (e.g. 0 for scalars, 1 for vectors) */ 812 valLen; /* number of scalars per data point, 813 -or- 0 to represent "this value will 814 be learned later at runtime" */ 815 int itemMax; /* such as GAGE_SCL_ITEM_MAX */ 816 gageItemEntry *table; /* array of gageItemEntry's, indexed 817 by the item value, 818 -or- NULL if the table cannot be 819 statically allocated (not because it 820 can come in different sizes, but 821 because it needs to be a modified 822 version of the compile-time table */ 823 void (*iv3Print)(FILE *, /* such as _gageSclIv3Print() */ 824 gageContext *, 825 gagePerVolume *), 826 (*filter)(gageContext *, /* such as _gageSclFilter() */ 827 gagePerVolume *), 828 (*answer)(gageContext *, /* such as _gageSclAnswer() */ 829 gagePerVolume *), 830 /* for allocating, copying, updating, and nixing the pervolume->data */ 831 /* pvlDataNew and pvlDataCopy can use biff, but: 832 --> they must use GAGE key (and not callback's library's key), and 833 --> pvlDataNix can not use biff */ 834 *(*pvlDataNew)(const struct gageKind_t *kind), 835 *(*pvlDataCopy)(const struct gageKind_t *kind, const void *pvlDataOld), 836 *(*pvlDataNix)(const struct gageKind_t *kind, void *pvlData); 837 int (*pvlDataUpdate)(const struct gageKind_t *kind, 838 const gageContext *ctx, 839 const gagePerVolume *pvl, 840 const void *data); 841 void *data; /* extra information about the kind of 842 volume that's being probed. Likely 843 used by filter, answer, and the 844 pvlData functions */ 845 } gageKind; 846 847 /* 848 ******** gageItemSpec struct 849 ** 850 ** dumb little way to store a kind/item pair. Formerly known 851 ** as a gageQuerySpec. 852 */ 853 typedef struct { 854 const gageKind *kind; /* the kind of the volume to measure */ 855 int item; /* the quantity to measure */ 856 } gageItemSpec; 857 858 /* 859 ******** gageItemPack struct 860 ** 861 ** A way of coherently representing a related set of gageItems. This 862 ** kind of inter-item relationship may eventually be part of the 863 ** definition of a gageKind. 864 ** 865 ** These are intended to be set up at compile-time, just like nearly all 866 ** the gageKinds. 867 ** 868 ** This is not to be confused with the gageMultiItem, which is (or 869 ** will be) a list of items, set at run-time, to learn from a given 870 ** volume. 871 */ 872 typedef struct { 873 const gageKind *kind; 874 int item[GAGE_ITEM_PACK_PART_MAX+1]; 875 } gageItemPack; 876 877 /* 878 ******** gageStackBlurParm struct 879 ** 880 ** all parameters associated with blurring one volume to form a "stack" 881 */ 882 typedef struct { 883 unsigned int num; /* # of blurring scales == # volumes */ 884 double *scale; /* scale parameter for each blurred volume */ 885 double sigmaMax, /* nrrdKernelDiscreteGaussian is implemented with 886 airBesselInExpScaled, which has numerical issues 887 for very large kernel sizes. Instead of doing 888 the blurring in one step, the diffusion is done 889 iteratively, with steps in diffusion time 890 of sigmaMax^2 */ 891 padValue; /* padding value for nrrdBoundaryPad */ 892 NrrdKernelSpec *kspec; /* NOTE: parm[0] will get over-written as part 893 of running the resampler at each scale */ 894 int dataCheck, /* when checking given stack to see if its the 895 blurring of a volume, actually go in and look at 896 values of first (probably the least blurred) 897 volume */ 898 boundary, /* passed to nrrdResampleBoundarySet */ 899 renormalize, /* passed to nrrdResampleRenormalizeSet */ 900 verbose; 901 } gageStackBlurParm; 902 903 /* 904 ******** gageOptimSigParm struct 905 ** 906 ** a fairly disorganized mess of parameters. under construction 907 */ 908 typedef struct { 909 /* INPUT ------------------------- */ 910 /* these determine the allocation and (slow) computation of ntruth */ 911 unsigned int dim; /* either 1, 2, or 3 */ 912 double sigmaMax, /* highest sigma in current computation */ 913 cutoff; /* parm[1] for discrete gaussian kernel */ 914 unsigned int measrSampleNum; /* how many samples along sigma to use 915 for measurements of error */ 916 917 /* these can be changed more often */ 918 unsigned int sampleNum; /* how many scale samples to optimize */ 919 int volMeasr, /* how to measure error at each reconstructed 920 scale (interpolated volume) */ 921 lineMeasr, /* how to summarize errors across all scales */ 922 plotting, /* we're plotting, not optimizing */ 923 tentRecon; /* for plotting: use tent instead of hermite */ 924 unsigned int maxIter; /* allowed iterations in optimization */ 925 double convEps; /* convergence threshold */ 926 927 /* INTERNAL ------------------------- */ 928 /* these related to the allocation and (slow) computation of ntruth */ 929 unsigned int sx, sy, sz; /* volume size for testing */ 930 double *sigmatru, /* sigmas for all lines of ntruth, allocated 931 for measrSampleNum */ 932 *truth; /* data pointer of ntruth */ 933 Nrrd *ntruth, /* big array of all truth volumes, logically 934 a sx x sy x sz x measrSampleNum array */ 935 *nerr, /* line of all errors, across scale */ 936 *ntruline, /* *wrapper* around some scanline of ntruth */ 937 *ninterp, /* last recon result */ 938 *ndiff; /* diff between recon and truth, single vol */ 939 /* most of these allocated according sampleNumMax */ 940 unsigned int sampleNumMax; /* largest number of SS samples to look at; 941 this is set at parm creation time and can't 942 be changed safely during parm lifetime */ 943 double *scalePos, /* current SS sample locations, allocated 944 for sampleNumMax */ 945 *step; /* per-point stepsize for descent */ 946 Nrrd **nsampvol; /* current set of SS samples, allocated for 947 sampleNumMax */ 948 gagePerVolume *pvl, **pvlSS; /* for gage; pvlSS allocation different than 949 scalePos or nsampvol */ 950 gageContext *gctx; /* context around nsamplevol */ 951 952 /* OUTPUT ------------------------- */ 953 double finalErr; /* error of converged points */ 954 } gageOptimSigParm; 955 956 /* defaultsGage.c */ 957 GAGE_EXPORT const char *gageBiffKey; 958 GAGE_EXPORT int gageDefVerbose; 959 GAGE_EXPORT double gageDefGradMagCurvMin; 960 GAGE_EXPORT int gageDefRenormalize; 961 GAGE_EXPORT int gageDefCheckIntegrals; 962 GAGE_EXPORT int gageDefK3Pack; 963 GAGE_EXPORT int gageDefCurvNormalSide; 964 GAGE_EXPORT double gageDefKernelIntegralNearZero; 965 GAGE_EXPORT int gageDefDefaultCenter; 966 GAGE_EXPORT int gageDefStackUse; 967 GAGE_EXPORT int gageDefStackNormalizeRecon; 968 GAGE_EXPORT int gageDefStackNormalizeDeriv; 969 GAGE_EXPORT double gageDefStackNormalizeDerivBias; 970 GAGE_EXPORT double gageDefStackBlurSigmaMax; 971 GAGE_EXPORT int gageDefOrientationFromSpacing; 972 GAGE_EXPORT int gageDefGenerateErrStr; 973 974 /* miscGage.c */ 975 GAGE_EXPORT const int gagePresent; 976 GAGE_EXPORT double gageZeroNormal[3]; 977 GAGE_EXPORT const airEnum *const gageErr; 978 GAGE_EXPORT const airEnum *const gageKernel; 979 GAGE_EXPORT const airEnum *const gageItemPackPart; 980 GAGE_EXPORT void gageParmReset(gageParm *parm); 981 GAGE_EXPORT void gagePointReset(gagePoint *point); 982 GAGE_EXPORT gageItemSpec *gageItemSpecNew(void); 983 GAGE_EXPORT void gageItemSpecInit(gageItemSpec *isp); 984 GAGE_EXPORT gageItemSpec *gageItemSpecNix(gageItemSpec *isp); 985 986 /* kind.c */ 987 GAGE_EXPORT int gageKindCheck(const gageKind *kind); 988 GAGE_EXPORT unsigned int gageKindTotalAnswerLength(const gageKind *kind); 989 GAGE_EXPORT unsigned int gageKindAnswerLength(const gageKind *kind, int item); 990 GAGE_EXPORT int gageKindAnswerOffset(const gageKind *kind, int item); 991 GAGE_EXPORT int gageKindVolumeCheck(const gageKind *kind, const Nrrd *nrrd); 992 993 /* print.c */ 994 GAGE_EXPORT void gageQueryPrint(FILE *file, const gageKind *kind, 995 gageQuery query); 996 997 /* sclfilter.c */ 998 typedef void (gageScl3PFilter_t)(gageShape *shape, 999 double *iv3, double *iv2, double *iv1, 1000 double *fw00, double *fw11, double *fw22, 1001 double *val, double *gvec, double *hess, 1002 const int *needD); 1003 GAGE_EXPORT gageScl3PFilter_t gageScl3PFilter2; 1004 GAGE_EXPORT gageScl3PFilter_t gageScl3PFilter4; 1005 GAGE_EXPORT gageScl3PFilter_t gageScl3PFilter6; 1006 GAGE_EXPORT gageScl3PFilter_t gageScl3PFilter8; 1007 GAGE_EXPORT void gageScl3PFilterN(gageShape *shape, int fd, 1008 double *iv3, double *iv2, double *iv1, 1009 double *fw00, double *fw11, double *fw22, 1010 double *val, double *gvec, double *hess, 1011 const int *needD); 1012 1013 /* scl.c */ 1014 GAGE_EXPORT const airEnum *const gageScl; 1015 /* HEY: this should perhaps be a "const gageKind", but pointers for 1016 general kinds may want to point to this, or to the return of 1017 a dynamic kind generator like tenDwiGageKindNew(), which is 1018 most certainly not "const gageKind". */ 1019 GAGE_EXPORT gageKind *const gageKindScl; 1020 GAGE_EXPORT const gageItemPack *const gageItemPackSclValue; 1021 1022 /* vecGage.c (together with vecprint.c, these contain everything to 1023 implement the "vec" kind, and could be used as examples of what it 1024 takes to create a new gageKind) */ 1025 GAGE_EXPORT const airEnum *const gageVec; 1026 GAGE_EXPORT gageKind *const gageKindVec; 1027 1028 /* shape.c */ 1029 GAGE_EXPORT void gageShapeReset(gageShape *shp); 1030 GAGE_EXPORT gageShape *gageShapeNew(void); 1031 GAGE_EXPORT gageShape *gageShapeCopy(const gageShape *shp); 1032 GAGE_EXPORT gageShape *gageShapeNix(gageShape *shape); 1033 GAGE_EXPORT int gageShapeSet(gageShape *shp, const Nrrd *nin, int baseDim); 1034 GAGE_EXPORT void gageShapeWtoI(const gageShape *shape, 1035 double index[3], const double world[3]); 1036 GAGE_EXPORT void gageShapeItoW(const gageShape *shape, 1037 double world[3], const double index[3]); 1038 GAGE_EXPORT int gageShapeEqual(const gageShape *shp1, const char *name1, 1039 const gageShape *shp2, const char *name2); 1040 GAGE_EXPORT void gageShapeBoundingBox(double min[3], double max[3], 1041 const gageShape *shape); 1042 1043 /* the organization of the next two files used to be according to 1044 what the first argument is, not what appears in the function name, 1045 but that's just a complete mess now */ 1046 /* pvl.c */ 1047 GAGE_EXPORT int gageVolumeCheck(const gageContext *ctx, const Nrrd *nin, 1048 const gageKind *kind); 1049 GAGE_EXPORT gagePerVolume *gagePerVolumeNew(gageContext *ctx, 1050 const Nrrd *nin, 1051 const gageKind *kind); 1052 GAGE_EXPORT gagePerVolume *gagePerVolumeNix(gagePerVolume *pvl); 1053 GAGE_EXPORT const double *gageAnswerPointer(const gageContext *ctx, 1054 const gagePerVolume *pvl, 1055 int item); 1056 GAGE_EXPORT unsigned int gageAnswerLength(const gageContext *ctx, 1057 const gagePerVolume *pvl, 1058 int item); 1059 GAGE_EXPORT int gageQueryReset(gageContext *ctx, gagePerVolume *pvl); 1060 GAGE_EXPORT int gageQuerySet(gageContext *ctx, gagePerVolume *pvl, 1061 gageQuery query); 1062 GAGE_EXPORT int gageQueryAdd(gageContext *ctx, gagePerVolume *pvl, 1063 gageQuery query); 1064 GAGE_EXPORT int gageQueryItemOn(gageContext *ctx, gagePerVolume *pvl, 1065 int item); 1066 1067 /* optimsig.c */ 1068 GAGE_EXPORT int gageOptimSigSet(double *scale, unsigned int num, 1069 unsigned int sigmaMax); 1070 GAGE_EXPORT gageOptimSigParm *gageOptimSigParmNew(unsigned int sampleMaxNum); 1071 GAGE_EXPORT gageOptimSigParm *gageOptimSigParmNix(gageOptimSigParm *parm); 1072 GAGE_EXPORT int gageOptimSigTruthSet(gageOptimSigParm *parm, 1073 unsigned int dim, 1074 double sigmaMax, double cutoff, 1075 unsigned int measrSampleNum); 1076 GAGE_EXPORT int gageOptimSigCalculate(gageOptimSigParm *parm, 1077 double *scalePos, unsigned int num, 1078 int volMeasr, int lineMeasr, 1079 double convEps, unsigned int maxIter); 1080 GAGE_EXPORT int gageOptimSigPlot(gageOptimSigParm *parm, Nrrd *nout, 1081 const double *plotpos, 1082 unsigned int plotPosNum, 1083 int volMeasr, int tentRecon); 1084 1085 /* stack.c */ 1086 GAGE_EXPORT double gageTauOfTee(double tee); 1087 GAGE_EXPORT double gageTeeOfTau(double tau); 1088 GAGE_EXPORT double gageSigOfTau(double tau); 1089 GAGE_EXPORT double gageTauOfSig(double sig); 1090 GAGE_EXPORT double gageStackWtoI(gageContext *ctx, double swrl, 1091 int *outside); 1092 GAGE_EXPORT double gageStackItoW(gageContext *ctx, double si, 1093 int *outside); 1094 GAGE_EXPORT int gageStackPerVolumeNew(gageContext *ctx, 1095 gagePerVolume **pvlStack, 1096 const Nrrd *const *nblur, 1097 unsigned int blnum, 1098 const gageKind *kind); 1099 GAGE_EXPORT int gageStackPerVolumeAttach(gageContext *ctx, 1100 gagePerVolume *pvlBase, 1101 gagePerVolume **pvlStack, 1102 const double *stackPos, 1103 unsigned int blnum); 1104 GAGE_EXPORT int gageStackProbe(gageContext *ctx, 1105 double xi, double yi, double zi, double si); 1106 GAGE_EXPORT int gageStackProbeSpace(gageContext *ctx, 1107 double x, double y, double z, double s, 1108 int indexSpace, int clamp); 1109 1110 /* stackBlur.c */ 1111 GAGE_EXPORT gageStackBlurParm *gageStackBlurParmNew(void); 1112 GAGE_EXPORT gageStackBlurParm *gageStackBlurParmNix(gageStackBlurParm *sbp); 1113 GAGE_EXPORT int gageStackBlurParmScaleSet(gageStackBlurParm *sbp, 1114 unsigned int num, 1115 double scaleMin, double scaleMax, 1116 int uniform, int optim); 1117 GAGE_EXPORT int gageStackBlurParmKernelSet(gageStackBlurParm *sbp, 1118 const NrrdKernelSpec *kspec, 1119 int renormalize); 1120 GAGE_EXPORT int gageStackBlurParmBoundarySet(gageStackBlurParm *sbp, 1121 int boundary, double padValue); 1122 GAGE_EXPORT int gageStackBlurParmVerboseSet(gageStackBlurParm *sbp, 1123 int verbose); 1124 GAGE_EXPORT int gageStackBlurParmCheck(gageStackBlurParm *sbp); 1125 GAGE_EXPORT int gageStackBlur(Nrrd *const nblur[], gageStackBlurParm *sbp, 1126 const Nrrd *nin, const gageKind *kind); 1127 GAGE_EXPORT int gageStackBlurCheck(const Nrrd *const nblur[], 1128 gageStackBlurParm *sbp, 1129 const Nrrd *nin, const gageKind *kind); 1130 GAGE_EXPORT int gageStackBlurGet(Nrrd *const nblur[], int *recomputedP, 1131 gageStackBlurParm *sbp, 1132 const char *format, 1133 const Nrrd *nin, const gageKind *kind); 1134 GAGE_EXPORT int gageStackBlurManage(Nrrd ***nblurP, int *recomputedP, 1135 gageStackBlurParm *sbp, 1136 const char *format, 1137 int saveIfComputed, NrrdEncoding *enc, 1138 const Nrrd *nin, const gageKind *kind); 1139 1140 /* ctx.c */ 1141 GAGE_EXPORT gageContext *gageContextNew(void); 1142 GAGE_EXPORT gageContext *gageContextCopy(gageContext *ctx); 1143 GAGE_EXPORT gageContext *gageContextNix(gageContext *ctx); 1144 GAGE_EXPORT void gageParmSet(gageContext *ctx, int which, double val); 1145 GAGE_EXPORT int gagePerVolumeIsAttached(const gageContext *ctx, 1146 const gagePerVolume *pvl); 1147 GAGE_EXPORT int gagePerVolumeAttach(gageContext *ctx, gagePerVolume *pvl); 1148 GAGE_EXPORT int gagePerVolumeDetach(gageContext *ctx, gagePerVolume *pvl); 1149 GAGE_EXPORT int gageKernelSet(gageContext *ctx, int which, 1150 const NrrdKernel *k, const double *kparm); 1151 GAGE_EXPORT void gageKernelReset(gageContext *ctx); 1152 GAGE_EXPORT int gageProbe(gageContext *ctx, double xi, double yi, double zi); 1153 GAGE_EXPORT int gageProbeSpace(gageContext *ctx, double x, double y, double z, 1154 int indexSpace, int clamp); 1155 1156 /* update.c */ 1157 GAGE_EXPORT int gageUpdate(gageContext *ctx); 1158 1159 /* st.c */ 1160 GAGE_EXPORT int gageStructureTensor(Nrrd *nout, const Nrrd *nin, 1161 int dScale, int iScale, int dsmp); 1162 1163 /* deconvolve.c */ 1164 GAGE_EXPORT int gageDeconvolve(Nrrd *nout, double *lastDiffP, 1165 const Nrrd *nin, const gageKind *kind, 1166 const NrrdKernelSpec *ksp, int typeOut, 1167 unsigned int maxIter, int saveAnyway, 1168 double step, double epsilon, int verbose); 1169 GAGE_EXPORT int gageDeconvolveSeparableKnown(const NrrdKernelSpec *ksp); 1170 GAGE_EXPORT int gageDeconvolveSeparable(Nrrd *nout, const Nrrd *nin, 1171 const gageKind *kind, 1172 const NrrdKernelSpec *ksp, 1173 int typeOut); 1174 1175 #ifdef __cplusplus 1176 } 1177 #endif 1178 1179 #endif /* GAGE_HAS_BEEN_INCLUDED */ 1180