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.
Terminology
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.
- "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 gageKinds,
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.
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.
Differentiation
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 */
gageKernelLast
};
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.
- 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.
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.
- 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.
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.
