Get Teem at Fast, secure and Free Open Source software downloads
to teem


How teem measures nrrds

Filtered point-sampling in volumes

Gage is a library for probing quantities at arbitrary points in a rectangular lattice volumes, based on separable convolution. It can measure values, and/or take first and second derivatives, and then do whatever math is needed to generate further information from the convolution measurements. For example, once Gage has measured the gradient and Hessian of a scalar field, it can compute intrinsic surface curvature, and all the volume renderings in this paper on curvature-based transfer functions were ultimately based on measurements from Gage. The kernels used to do the measurements are the NrrdKernel defined by Nrrd.

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.


Here basic Gage terminology is described in a bottom-up way, with reference to the C data structures that implement the ideas. It may be helpful to refer to the teem/src/gage/gage.h header file itself.

Self-contained usage example

Here's a code snippet that loads a scalar volume, and measures the gradient at some point, using the Catmull-Rom cubic spline and its first derivative for the kernels. It is very careful about error-checking (but not very careful about memory leaks in case of errors: one should free() the return of biffGetDone(), and gageContextNix(ctx)).
char 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 */
Note that it is absolutely necessary to call gageUpdate() to refresh all internal state in the gageContext prior to calling gageProbe(). Note also that gageProbe() does not use Biff for error handling for speed reasons, but in the case of a problem (such as sampling outside the volume), then error descriptions end up in ctx->errNum and ctx->errStr.

The gage_t typedef is no more

There used to be a gage_t typedef that was used for the return type of all the gage answer computations. This became annoying. After the Teem-1.9 release, double served this role, and gage_t was removed. Apologies for any broken code.


The kernels used in the convolution-based measurements are of two basic types: those for "reconstructing values" and those for "measuring derivatives". Reconstructing values usually means simple interpolation, but not if some smoothing is desired. "Measuring derivatives" means measuring a derivative of the sampled field by convoluting with the derivative of a reconstruction kernel.

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 */
NOTE: This is all well and good, but currently Gage only has support for using gageKernel00, gageKernel11, and gageKernel22. The technical term for this is "3-pack filtering". The ability to use all six kinds of kernels is called "6-pack filtering". The plan is to implement 6-pack filtering for Gage in Teem 2.0.

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 steps of gageProbe()

Here is a rough description of what happens in response to gageProbe.
  1. 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.
  2. 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.
  3. 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.
  4. 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.

The Answer array

After gageProbe() returns, it is the user's responsibility to copy the answer from the PerVolume. The method chosen for storing the results of the filter and answer stages of probing was that of a common scratch-space: a long flat array in which there is space for all answers for all Items in the Kind. This is called the "answer array", stored in gagePerVolume->answer. For convenience, there is also a per-item array of pointers into the answer array, which make it easier to copy a particular item's answer. The user calls gageAnswerPointer() to get this pointer.

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.

Parameter settings: gageParmSet()

There are various ways of modifying the behavior of Gage, with parameters that are set via gageParmSet(). The sensible time to do most of these is shortly after context creation.

Multi-threaded Gage usage

Through the use of the gageContext, Gage was designed to be thread-safe in its operation, so that multiple threads can probe the same volume (without copying the volume), and without clobbering each other's state. The API for this functionality is a little clumsy, but it works for now, and will probably be improved later.

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.