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.

Background

Within the imaging pipeline, interpolation is performed by vtkImageReslice and its subclass vtkImageResample. These classes are use interpolation to resample a whole image, they are not useful when a user wants to interpolate just a single point. 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.

The interpolator code is currently under review: [1]

Hierarchy

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

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 using with an image filter like vtkImageReslice)
• 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 will values will be silently 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

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

Sinc interpolation allows perfect reconstruction of a band-limited signal from discretely sampled data. The vtkImageSincInterpolator provides a very close approximation of sinc interpolation , where the approximation improves and the kernel size increases.

• SetWindowFunctionToLanczos() - the default
• SetWindowFunctionToKaiser()
• SetWindowHalfWidth(int width) - default is 3, max is 7

Notes:

1. The kernel size is twice the half-width, so the default kernel size is 6.
2. All window functions use precomputed lookup tables for efficiency.
3. The Hamming and Blackman windows are likely to be added soon.
4. Tunable parameters for Kaiser are also likely to be added.
5. Optional bandlimiting to the Nyquist frequency for the output sampling will be added soon. By default, the sinc width is determined by the Nyquist frequency for the input sampling. Satisfying Nyquist for the output sampling will provide antialiasing in situations where the image is subsampled.

vtkImageBSplineInterpolator

I have a version of vtkImageReslice posted in the VTK Journal that performs b-spline interpolation. At some point in the future (no date has been set), this code will be refactored to create a vtkImageBSplineInterpolator that will become part of the main VTK code base.

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:

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 mainly 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.

= An efficient image resize filter

As mentioned in the section above on separable interpolators, interpolation can be done in O(N) time instead of O(N^3), where N is the kernel size. Currently, vtkImageReslice and vtkImageResample use a hybrid method that computes the interpolation coefficients in O(N) time, but perform the convolution in O(N^2) or O(N^3) time depending on the image dimensionality. I have written a vtkImageResize filter that performs the entire computation 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.

B-Splines

Also as mentioned above, a B-Spline interpolator will be added to VTK in the future. No timeline has been set, but the amount of work required is fairly small, as a version of vtkImageReslice with B-Spline interpolation code already exists.

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.

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 finish the work eventually, but do not know when that would be. An abstract data interpolator would have a method like this:

Interpolate(vtkPointData *inPD, vtkPointData *outPD, int 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.

The plan would be to have a hierarchy like this:

• vtkAbstractDataSetInterpolator - abstract interface
  • vtkDataSetInterpolator - provides existing vtkCell interpolation methods
  • vtkAbstractImageInterpolator - additional interfaces for image data
    • vtkImageInterpolator - basic image interpolation