## Gage | |||

teem |

- Terminology
- Self-contained usage example
- The gage_t typedef is no more
- Differentiation
- The steps of gageProbe()
- The Answer array
- Parameter settings: gageParmSet()
- Multi-threaded Gage usage

The Gage library itself can do measurements in scalar and
3-vector volumes, but what makes Gage useful is that it provides
a more general framework for doing measurements
in any other kind of volume. Through this means, the **Ten** library
provides data structures that allow Gage to be used in
diffusion tensor volumes, as well as volumes of diffusion-weighted
images.

- "Item": Any single quantity (scalar or non-scalar) that can be measured.
Users specify items in terms of the
`int`values from the`gageScl*`or`gageVec*`enums (for example). Internally Gage uses the`gageItemEntry`to represent everything about the quantity, such as being scalar or non-scalar, what derivative measurements it requires, what other items are pre-requisites. - "Kind" ==
`gageKind`: The kind of volume dataset that Gage can measure (scalar, vector, tensor, etc.). The`gageKind`is essentially a table of all the items that can be measured in a volume, together with pointers to callbacks which do the work of convolution and calculating the items. Gage supplies`gageKindScl`and`gageKindVec`, for scalar and 3-vector volumes, respectively, but the rest of Gage just takes pointers to`gageKind`s, so they can come from any other Teem library or from the user. - "Item Specification" ==
`gageItemSpec`: The pairing of a Kind and an (integral) Item, which completely specifies a quantity that Gage can measure. - "Query" ==
`gageQuery`: A set of items to be measured at each probe location in a single volume of some Kind. - "Shape" ==
`gageShape`: The dimensions and spacing (or dimensions and orientation) of a volume. If a volume doesn't have full-fledged orientation information (via the`NrrdAxisInfo->spaceDirection`), then Gage invents a notion of axis-aligned world-space in which the volume is isotropically scaled to fit inside [-1,1]^3. By testing whether two volumes have the same shape, Gage can tell whether it makes sense to probe them in tandem. - "PerVolume" ==
`gagePerVolume`: All the parameters and state specific to doing measurements in a single volume of some Kind, including the query, caches for intermediate convolution results, and a space in which all answers are recorded. Gage reads values from the given Nrrd volume, not a copy of it, so changes to the volume are seen immediately. It is the user's responsibility to copy the answers from here after`gageProbe()`has returned. - "Kernel" ==
`NrrdKernel`: Not part of Gage, but fundamental to how it works; see the Nrrd kernel documentation. Gage uses different kernels for reconstructing values (by interpolation, or by blurring), and for measuring first and second derivatives. - "Context" ==
`gageContext`: The top-level data structure that manages all the state relevant to probing one or more volumes, such as which kernels are used, and which derivatives need to be measured. When probing multiple volumes, all volumes must have the same Shape, but they can be of different Kinds.

Note that it is absolutely necessary to callchar me[]="gagedemo"; gageContext *ctx; gagePerVolume *pvl; const double *grad; Nrrd *nin; double kparm[NRRD_KERNEL_PARMS_NUM]; int E; nin = nrrdNew(); if (nrrdLoad(nin, "volume.nrrd", NULL)) { fprintf(stderr, "%s: couldn't load:\n%s\n", me, biffGetDone(NRRD)); return 1; } if (gageKindVolumeCheck(gageKindScl, nin)) { fprintf(stderr, "%s: didn't get a %s volume:\n%s\n", me, gageKindScl->name, biffGetDone(GAGE)); return 1; } kparm[0] = 1.0; /* scale parameter, in units of samples */ kparm[1] = 0.0; /* B */ kparm[2] = 0.5; /* C */ ctx = gageContextNew(); E = 0; if (!E) E |= !(pvl = gagePerVolumeNew(ctx, nin, gageKindScl)); if (!E) E |= gagePerVolumeAttach(ctx, pvl); if (!E) E |= gageKernelSet(ctx, gageKernel00, nrrdKernelBCCubic, kparm); if (!E) E |= gageKernelSet(ctx, gageKernel11, nrrdKernelBCCubicD, kparm); if (!E) E |= gageQueryItemOn(ctx, pvl, gageSclGradVec); if (!E) E |= gageUpdate(ctx); if (E) { fprintf(stderr, "%s: trouble:\n%s\n", me, biffGetDone(GAGE)); return 1; } grad = gageAnswerPointer(ctx, pvl, gageSclGradVec); if (gageProbe(ctx, 1.1, 2.3, 6.8)) { fprintf(stderr, "%s: trouble:\n(%d) %s\n", me, ctx->errNum, ctx->errStr); return 1; } printf("the gradient is (%g,%g,%g)\n", grad[0], grad[1], grad[2]); /* this blows away attached pervolumes, too, but it doesn't do anything to the Nrrds inside the pervolumes */ ctx = gageContextNix(ctx); pvl = NULL; nin = nrrdNuke(nin); /* free the Nrrd and its memory */

Measuring a derivative in a volume is always accomplished by partial
derivatives, such as reconstruction in Y and Z and a first derivative
in X, to compute a partial derivative with respect to to X. However, the
reconstruction that is done in the context of derivative measurement
(such as along Y and Z for a first partial in X) can be different than
reconstruction of values alone; some smoothing may be desired, for
instance. Gage distinguishes between these two types of
reconstruction. The same applies for second derivatives, but in this
context one has three distinct kernels, for value reconstruction and
for first and second derivative measurement. The way that these
different kernels are identified in gage is with a two-digit scheme,
with the first number giving the context, and the second number giving
the measurement. These two-digits are used in the `gageKernel*`
enum that is used to identify the kind of kernel being specified:

/* ******** gageKernel* enum ** ** these are the different kernels that might be used in gage, regardless ** of what kind of volume is being probed. */ ontinousenum { gageKernelUnknown=-1, /* -1: nobody knows */ gageKernel00, /* 0: measuring values */ gageKernel10, /* 1: reconstructing 1st derivatives */ gageKernel11, /* 2: measuring 1st derivatives */ gageKernel20, /* 3: reconstructing 1st partials and 2nd deriv.s */ gageKernel21, /* 4: measuring 1st partials for a 2nd derivative */ gageKernel22, /* 5: measuring 2nd derivatives */ gageKernelLast };

A basic assumption in Gage is that the interesting quantities (the "items")
are either values in, or derivatives of, the given field, or something
that can be calculated from some combination of these. Because of
this assumption, every item in Gage can be characterized in terms of
the highest-order derivative upon which it depends, together with a
list of its computational prerequisites (and the derivatives upon
which they depend). All this information is stored in the `gageKind`.

- The probe point is located relative to the raster grid of the
volume, and a cube of values (within the support of the necessary
kernels) is copied into a per-volume value cache
(
`gagePerVolume->iv3`). The values are maintained if the last probe location was within the same voxel. - The relationship between the probe location and the image gride
determines the locations at which the c filters are evaluated
to compute the convolution. The filter sample locations are stored
(in
`gageContext->fsl`), and the filters are evaluated to determine the sample weights (in`gageContext->fw`). These values are re-used for all PerVolumes attached to the Context. - The "filter" stage (
`gageKind->filter()`): For each PerVolume, all the various value reconstruction and derivative measurement convolutions are performed, as dot products between the vectors of values in the value cache, and the filter weights. There is where the separability constraint is applied: the weight for each sample is the product of the weights for each axis-aligned filter. Care is taken to re-use intermediate convolution results (stored in`gagePerVolume->iv2`and`gagePerVolume->iv1`), so that no convolution computation is repeated. At the end of the filter stage, quantities like gradients and Hessians are known. - The "answer" stage (
`gageKind->answer()`): For each PerVolume, the results from the filtering stage are mathematically combined into new information. For example, answering a request for`gageSclNormal`requires normalizing the result of`gageSclGradVec`. No convolution happens at this stage.

Some Items are basically just pointers into the answer array of
a larger "parent" item. For example `gageSclHessEvec1`
is the second eigenvector of the Hessian of a scalar field, and
the answer for this is just an offset into the answer for
`gageSclHessEvec`, which is all three eigenvectors.

`gageParmSet(ctx, gageParmVerbose, verb)`

sets a verbosity level`verb`; 0 for silent, and higher values for more feedback. Used for debugging.`gageParmSet(ctx, gageParmRenormalize, AIR_TRUE/AIR_FALSE)`

determines whether to renormalize convolution weights at each probe location (which slows things down slightly). Some kernels are "accurate" in the sense that if you sample them at regular unit spacings, the sum of the samples is constant. This is true of`nrrdKernelTent`, but not of`nrrdKernelGaussian`. The role of renormalization is to make the discrete sum of filter weights equal what is known to be the continuous integral of the filter over its support. This takes on more significance in the context of differentiation: the kernel weights must sum to zero.`gageParmSet(ctx, gageParmCheckIntegrals, AIR_TRUE/AIR_FALSE)`

As a courtesy, verify that kernels being used for value reconstruction have a non-zero continuous integral, and those being used for differentiation have a zero integral.`gageParmSet(ctx, gageParmGradMagCurvMin, gmin)`

When measuring curvature information, if the gradient magnitude is below`gmin`, then don't bother computing curvature information as it will be probably be numerically meaningless.

Calling `gageContextCopy()` immediately after
`gageUpdate()` will result in a new context (complete with
attached PerVolumes), ready to be passed to `gageProbe()`. The
result of `gageContextCopy()` is not functionally the same as
the original: you can't call `gageContextCopy()` on it, and
more importantly, you can't currently modify the parameters of the
copy. Basically, the only thing you should do with a
`gageContext` copy is call `gageProbe()`. If any
parameters need to be changed, delete the copy with
`gageContextNix()`, make necessary `gage*Set()` calls,
then `gageUpdate()`, then `gageContextCopy()` again.
With time this will be fixed.