VTK/Image Interpolators

Image interpolation is done internally by several VTK classes, but there is no easy way for users to interpolate an image directly. Neither is there any straightforward way of adding new interpolation methods to VTK, since each VTK class that uses image interpolation has its own internal code for this purpose. The goal of this project is to add a vtkImageInterpolator class to VTK that provides a clean API for image interpolation.

The interpolation code was merged into VTK on Sept 15, 2011. It is in VTK 5.10 and later releases.

The author of this document is David Gobbi. He can be reached on the vtk-developers mailing list.

Background

Within the imaging pipeline, interpolation is performed by vtkImageReslice, vtkImageResample, and vtkImageResize. These classes use interpolation to resample a whole image, but are not useful when a user wants to interpolate just a few points. Also, the code for vtkImageReslice is quite complex, so adding new interpolation modes to this class is beyond the abilities of most VTK programmers, and any programmer who did so would either have to submit their changes to Kitware, or maintain their own version of vtkImageReslice.

The obvious way to resolve this issue is to move the interpolation code out of vtkImageReslice and into its own interpolator class (or classes), and add a SetInterpolator() method to vtkImageReslice. This is, in fact, something that I have wanted to do for around 10 years, but found very difficult to do because the vtkImageReslice interpolation code was very tightly integrated with the class for performance reasons. It was a requirement of this project that factoring out the interpolation code would not result in any significant loss of performance.

I had also planned to unify the interpolation code between the imaging pipeline and the geometry pipeline, so that these interpolators could be used with e.g. vtkProbeFilter, but decided to delay the unification until a future project. The main difference between the VTK image pipeline and the VTK geometry pipeline is that the former utilizes only scalars and structured coordinates for interpolation, while the latter utilizes all the arrays in the data set's PointData along with the point and connectivity information. It will be possible to provide a unified interface at some point in the future, by providing interpolation methods that operate on all of the PointData instead of just on the scalars.

Hierarchy

• vtkAbstractImageInterpolator - provides abstract interface
  • vtkImageInterpolator - basic linear, cubic, nearest-neighbor interpolation
  • vtkImageSincInterpolator - sinc interpolation
  • vtkImageBSplineInterpolator - b-spline interpolation

Python usage example with vtkImageReslice

interp = vtkImageSincInterpolator()
interp.SetWindowFunctionToLanczos()

reslice = vtkImageReslice()
reslice.SetInterpolator(interp)

Abstract Interface

vtkAbstractImageInterpolator

Basic interface methods

• Initialize(vtkDataObject *data) - set the data to be interpolated (if not used with an image filter)
• SetOutValue(double v) - set value to return for out-of-bounds lookup
• Update() - needs to be called if you call any Set method after having called Initialize()
• double Interpolate(double x, double y, double z, int component) - for Python/Tcl/Java
• bool Interpolate(const double point[3], double *value) - for use from C++
• int GetNumberOfComponents() - get the number of components returned by Interpolate
• ReleaseData() - release the data

Notes:

1. GetNumberOfComponents() gives the number of components that Interpolate(double point[3], double *value) should expect.
2. Interpolate(const double point[3], double *value) returns false if the point was out of bounds, while setting all components to the OutValue.
3. ReleaseData() releases the interpolator's reference to the image scalars.
4. The Initialize() method causes the interpolator to store a pointer to the image scalars array. It does not store a pointer to the whole vtkImageData object. The reason for this is garbage collection efficiency.
5. Calling Update() only updates the interpolator, it does not update the data. You must ensure the data has been updated before calling Initialize().

Advanced methods

• SetComponentOffset(int offset) - set the first component to be interpolated
• SetComponentCount(int count) - set the number of components to interpolate (default: all)
• SetBorderMode(int mode) - control lookups at or beyond the image bounds
• SetTolerance(double t) - set tolerance for considering points to be out-of-bounds

Notes:

1. SetComponentOffset()/Count() allow selection of which components to interpolate. By default, ComponentOffset is 0 and ComponentCount is -1 (all components). These values are clamped to ensure that the interpolator does not try to access components that aren't there.
2. SetBorderModeToClamp() means that any lookups beyond the image bounds but within the Tolerance will be set to the closest point on the image bounds.
3. SetBorderModeToRepeat() will wrap any out-of-bounds lookups to the opposite edge, unless they are beyond the Tolerance.
4. SetBorderModeToMirror() will mirror out-of-bounds lookups, unless they are beyond the Tolerance.
5. SetTolerance() sets the tolerance in fractional units, where 1.0 is one voxel.
6. The tolerance can safely be set to large values, it does not have to be between 0 and 1.
7. Update() must be called if any of these methods are called post-initialization.

Pipeline methods

Pipeline methods are meant to be called by VTK filters that utilize interpolator objects.

• int ComputeNumberOfComponents(int inputComponents) - compute output component count
• void ComputeSupportSize(const double matrix[16], int support[3]) - for update extent computation
• bool CheckBoundsIJK(const double x[3]) - check whether structured coords are within bounds
• void InterpolateIJK(const double x[3], double *value) - interpolate with structured coords
• void GetExtent(int extent[6]) - extent of data being interpolated
• void GetSpacing(double spacing[3]) - spacing of data being interpolated
• void GetOrigin(double origin[3]) - origin of data being interpolated

Notes:

1. The ComputeNumberOfComponents() method is necessary because the user can choose that only some of the components will be utilized.
2. ComputeSupportSize() allows a filter to compute the required support size, which is needed for computing update extents. The matrix is the structured coordinate transformation between output voxels and input voxels (e.g. for vtkImageReslice, it is the IndexMatrix). If the matrix is unknown or unknowable, then it should be set to NULL, and then this method will return the maximum possible support size for the interpolation kernel instead of the minimum possible support size.
3. CheckBounds() uses the tolerance to provide information about which points are out of bounds.
4. InterpolateIJK() ignores the tolerance, it always returns a valid value according to the BorderMode setting.

Separable kernels

Most interpolation methods are separable, meaning that they can be computed separately in the x, y, and z directions. This allows some or all of the interpolation computations to be computed in O(N) time with respect to the kernel size N, instead of O(N^3) time which is required for non-separable kernels.

• bool IsSeparable() - returns true for separable interpolators
• void PrecomputeWeightsForExtent(const double matrix[16], const int extent[6], int checkExtent[6], vtkInterpolationWeights *&weights)
• void FreePrecomputedWeights(vtkInterpolationWeights *&weights)
• void InterpolateRow(vtkInterpolationWeights *&weights, int x, int y, int z, double *value, int rowlen)

These are highly advanced methods, please see the header file for additional documentation.

Concrete Classes

vtkImageInterpolator

Cubic interpolation kernel.

The vtkImageInterpolator class is the default interpolator. It provides nearest-neighbor, linear, and cubic interpolation. I decided not to split each interpolation method into its own class, since doing so would have caused a large amount of code duplication without adding any benefit. One class is easier to maintain than three, in this case.

• SetInterpolationModeToNearest()
• SetInterpolationModeToLinear() - the default
• SetInterpolationModeToCubic() - C(1) continuous

Notes:

1. If the interpolation mode is changed after a call to Initialize(), then you must call Update() before using the interpolator.

vtkImageSincInterpolator

Lanczos interpolation kernel, renormalized. This kernel provides a small amount of edge enhancement.

Blackman interpolation kernel, renormalized. This kernel does no edge enhancement.

Kaiser interpolation kernels. The Blackman, BlackmanHarris3, and BlackmanNuttall3 are nearly identical to Kaiser with alpha=3 and WindowHalfWidth=3. Likewise, BlackmanHarris4, BlackmanNuttall4, and Nuttall are nearly identical to Kaiser with alpha=4 and WindowHalfWidth=4 (not shown).

Sinc interpolation allows perfect reconstruction of a band-limited signal from discretely sampled data. The vtkImageSincInterpolator provides a very close approximation to sinc interpolation, where the approximation improves as the kernel size increases. Since these are only approximations of sinc interpolation, there is a tradeoff between image blurring and ringing, where any decrease in blurring results in a corresponding increase in ringing. The Lanczos, Cosine, Hann, and Hamming window functions provide good sharpness but at the cost of a small amount of ringing. The Blackman family of windows are closer to the ideal, they do not provide quite as much sharpness but have almost no ringing. The Kaiser window has a parameter to allow you to control the tradeoff between blurring and ringing. The rationale behind these windows is described in Nuttall AH, "Some Windows with Very Good Sidelobe Behavior," IEEE Transactions on Acoustics, Speech, and Signal Processing, Vol. ASSP-29(1), 1981.

• SetWindowFunctionToLanczos() - the default, it does some edge enhancement
• SetWindowFunctionToKaiser() - call UseWindowParameterOn() and SetWindowParameter(3)
• SetWindowFunctionToCosine()
• SetWindowFunctionToHann()
• SetWindowFunctionToHamming()
• SetWindowFunctionToBlackman() - use WindowHalfWidth of 3 or higher
• SetWindowFunctionToBlackmanHarris3() - use WindowHalfWidth of 3 or higher
• SetWindowFunctionToBlackmanHarris4() - use WindowHalfWidth of 4 or higher
• SetWindowFunctionToNuttall() - use WindowHalfWidth of 4 or higher
• SetWindowFunctionToBlackmanNuttall3() - use WindowHalfWidth of 3 or higher
• SetWindowFunctionToBlackmanNuttall4() - use WindowHalfWidth of 4 or higher
• SetWindowHalfWidth(int width) - default is 3, max is 16
• UseWindowParameterOn() - allow tuning of the window
• SetWindowParameter(double alpha) - tune the window, requires UseWindowParameterOn()
• SetBlurFactors(double f[3]) - blur the image while interpolating
• AntialiasingOn() - automatically set the BlurFactors so that the cut-off frequency is the Nyquist frequency
• RenormalizationOff() - do not renormalize the kernel before using it

Notes:

1. The kernel size is twice the half-width, so the default kernel size is 6.
2. The number of sinc lobes is one less than the kernel size, so the default is 5 sinc lobes.
3. All window functions use precomputed lookup tables for efficiency.
4. The Kaiser window can be tuned by setting the WindowParameter.
5. The blur factors allow bandlimiting in the x, y, and z directions (blur factor must be greater than or equal to 1). Note that the size of the kernel is automatically increased by the blur factor, so a large blur factor increases the computation time.
6. Antialiasing automatically sets the blur factors so that Nyquist is satisfied for the output sampling rate. When antialiasing is on, the blur factors are computed when ComputeSupportSize() is called (this method is automatically called by any image filters that use the interpolator).
7. The equations for most of the Sinc kernels result in interpolation coefficients that do not sum to unity, in fact the sum of the coefficients is a function of the fractional offset and results in small ripples in the output image. By default, vtkImageSincInterpolator renormalizes the coefficients to eliminate these ripples, but by calling RenormalizationOff() you can skip the renormalization step.

vtkImageBSplineInterpolator

BSpline interpolation kernels of varying degree. Only the inner lobes are shown here.

The b-spline interpolator provides interpolation of images with b-splines of degree 2 through 9, with a default degree of 3, which provides a cubic b-spline.

• SetSplineDegree(int d)

Unlike the other interpolators, b-splines cannot operate on the image directly, instead the image must first be pre-filtered with vtkImageBSplineCoefficients to compute the control points for the b-spline (there will be exactly one control point per image voxel, so filtering with vtkImageBSplineCoefficients results in an image that is the same size as the input image). The b-spline kernels are approximately Gaussian in shape, and applying a b-spline without pre-filtering the image will result in a blurring of the image. The vtkImageBSplineCoefficients class is described in Uniform B-Splines for the VTK Imaging Pipeline (VTK Journal publication).

Implementation Details

The interpolator object has a pair of function pointer members called InterpolationFuncFloat and InterpolationFuncDouble that point to the functions that perform the interpolation. These function pointers are set by the Update() method, which chooses which function to use based on the data type (short, float, etc) and the interpolation mode.

The InterpolateIJK() method is inlined, it simply calls the function pointer and imposes no overhead:

inline void vtkAbstractImageInterpolator::InterpolateIJK(const double point[3], double *value)
{
  this->InterpolationFuncDouble(this->InterpolationInfo, point, value);
}

The InterpolationInfo member is a pointer to a vtkInterpolationInfo helper struct, which contains information about the image data such as the scalar pointer and the increments. I chose to use function pointers instead of virtual methods because, with a function pointer, I can easily predict what code the compiler will produce for the function call. This was important to me because the Interpolate methods are usually called in very tight, performance-critical loops. It also allows the interpolation functions to be written in pure C code, though I am not yet certain what the value of that is.

Future Possibilities

An efficient image resize filter

As mentioned in the section above on separable interpolators, for an interpolation kernel of size N it is possible to resample an image in O(N) time instead of O(N^3). Currently, vtkImageReslice and vtkImageResample use a hybrid method where, when resizing an image, the interpolation coefficients are computed in O(N) but the convolution is performed in O(N^2) or O(N^3) time depending on the image dimensionality. I have written a vtkImageResize filter that performs the entire operation in O(N) time with minimal memory overhead through the use of a sliding window memory buffer. I can add this filter to VTK as soon as the image interpolators have been added. Added on Dec 29, 2011.

Gradients

It would be possible to add methods to the interpolator classes so that they could return the gradient of the image at any (X,Y,Z) point. This would allow a future vtkVolumeMapper class to be written that used these interpolators, enabling VTK volume rendering to performed with higher-order interpolation kernels. It would also allow the interpolators to be used with vtkImplicitVolume, which requires gradient computation. There is no set timeline for this feature.

GPU Interpolators

For the GPU to use an interpolator, the interpolation kernel would have to be converted into an array that could then be loaded onto the GPU as a texture. In the future a method could be added to the interpolator classes that would compute such an array. Rendering classes could then load this array into OpenGL and use it via a fragment shader. There is no timeline set for this, in fact it is unlikely that I will add such a feature myself.

Interpolators for the geometry pipeline

VTK's abstract data interface provides methods for interpolating any cell type, but there is no way for people to apply their own interpolation strategies. I began work on an abstract vtkDataInterpolator, but had to stop because it was taking too much time. I hope to continue the work eventually, but for now all of the information provided here is purely hypothetical and is subject to review. An abstract data interpolator would have a method like this:

Interpolate(vtkPointData *inPD, vtkPointData *outPD, vtkIdType outPointID, const double point[3])

The interpolator would have to internally store info about the points and cells (and it should probably store a locator as well). It would interpolate the input PointData at position (x,y,z) and then store the results in the output PointData at location outPointID.

Since it might be necessary to store some state information for efficiency, e.g. for cell-walking, the interpolation method should probably take a threadID so that concrete interpolation strategies can store per-thread state information:

Interpolate(vtkPointData *inPD, vtkPointData *outPD, vtkIdType outPointID, const double point[3], int threadID)

After the abstract data interpolator is complete, the interpolator hierarchy would look like this:

• vtkAbstractDataSetInterpolator - abstract interface
  • vtkDataSetInterpolator - provides existing vtkCell interpolation methods
  • vtkMyNamedInterpolator - various custom interpolators
  • vtkAbstractImageInterpolator - additional interfaces for image data
    • vtkImageInterpolator - basic image interpolation
    • vtkImageSincInterpolator
    • vtkImageBSplineInterpolator
    • vtkImageMyNamedInterpolator - other image interpolators