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)

Abstract Interface

vtkAbstractImageInterpolator

The basic set of interface methods are as follows

• Initialize(vtkDataObject *data) - set the data to be interpolated
• SetTolerance(double t) - set tolerance for considering points to be out-of-bounds
• SetOutValue(double v) - set value to return on out-of-bounds lookup
• Update() - needs to be called after any Set method is used
• double Interpolate(double x, double y, double z, int component) - basic interface
• bool Interpolate(const double point[3], double *value) - advanced interface
• int GetNumberOfComponents() - get the number of components
• ReleaseData() - release the data

Notes:

1. The Initialize() method calls Update() automatically. However, if you call any Set method of the interpolator, you must call Update() again.
2. Calling Update() only updates the interpolator, it does not update the data. You must ensure the data has been updated before calling Initialize(). You must call Initialize() again whenever the extent, spacing, or any other image information changes.
3. The Initialize() method causes the interpolator to store a pointer to the image scalars array. It does not store a pointer to the vtkImageData object, because doing so is unnecessary and would, in many use cases, cause a reference loop that VTK's garbage collector would have to untangle.
4. GetNumberOfComponents() gives the number of components that Interpolate(double point[3], double *value) should expect.
5. the Interpolate(const double point[3], double *value) method returns false if the point was out of bounds, while setting all components to the OutValue.
6. SetTolerance() sets the tolerance in fractional units, where 1.0 is one voxel
7. ReleaseData() releases the interpolator's reference to the image scalars

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

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. The tolerance can safely be set to large values, it does not have to be between 0 and 1
6. 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

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.

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.