In addition, it would be useful to have a new constant mode.
For comparison, scipy has a zoom function for up-sampling arrays, and has a mode option for selecting the desired mode, with the constant mode used by default.
Motivation:
A recent PyVista discussion revealed that using vtk’s cubic interpolator produces substantially different (and suprising) results when compared to scipy’s zoom interpolator (with order=3 for cubic), and the use (or absence) of a constant mode was found to be the reason why.
In VTK, SetBorderMode() method is a public method, and pyvista should be able to expose it via resample() without any changes to VTK itself. I’m not sure why pyvista does not already do so.
There are reasons why VTK doesn’t implement the constant mode:
it requires an additional piece of information (the constant)
having a default constant value of zero can lead to surprising results
To illustrate (2), consider a typical CT image. Let’s say the raw CT data is in the range [0,4095] and that we also have a copy of the CT image converted to Hounsfield units in the range [-1024,3071] (here, the conversion is done by subtracting 1024 from the raw data). Ideally, we expect to get the same results (to within roundoff error) if we interpolate before converting to Housfield units as compared to interpolating after converting to Hounsfield units.
However if we use the constant border mode with a constant of zero, we would get different results! In order to get the same results, the user would have to be savvy enough to use a constant of zero when working in Hounsfield units, and a constant of 1024 when working with the raw pixel data.
The clamp mode gives the same result in either case, because the clamp mode is invariant under linear transformations of the pixel values. The same is true of the repeat and mirror methods (and the reflect method). The constant mode stands out from the rest because it is not invariant.
Also, consider what most of the border modes are trying to achieve: they are extrapolating the data values past the edge of the image. If the image is smooth (as most medical images are) then we expect that the data values in the border will have similar values to their neighbors in the image. That’s what the clamp mode provides.
I’m not actually opposed to adding a constant mode to vtkImageInterpolator, I’m just trying to illustrate why it isn’t my first choice when interpolating images. Any of the other modes will generally give superior results.
In VTK, SetBorderMode() method is a public method, and pyvista should be able to expose it via resample() without any changes to VTK itself. I’m not sure why pyvista does not already do so.
Thanks for the info! For some reason I couldn’t see it in the docs, I think I was looking at the abstract interpolator docs VTK: vtkImageInterpolator Class Reference. Will be sure to add this option to PyVista’s API.
As for your comments/reasons for not including a constant mode by default, they do make sense, especially for the CT context with Hounsfield units. A similar argument can probably be made for the default value for other image filters like vtkImageConstantPad and vtkImageThreshold where using a default of 0 may not be appropriate for some use cases. (I don’t think these vtk filters actually make this default assumption, but PyVista’s wrappers of these filters do).
When creating filters for PyVista, there is often a tussle between
strictly adhering to VTK’s implementation and defaults
adding quality-of-life changes or defaults
making the API more similar to other pythonic libraries like scipy
In this case, I think we’ll probably stick with the vtk default for the border mode, and add a few notes to the documentation about the effects this has on the output.
Given that a constant mode is not yet available for vtkImageInterpolator, is it sensible to emulate this by simply using vtkImageConstantPad before resampling? And if so, how much padding should be given? E.g. in the linked pyvista discussion, padding the image with zeros seemed to yield different results depending on the number of zeros that the input was padded with… Maybe padding with a single zero at the borders and setting BorderMode to repeat is sufficient to emulate this?
I just want to note that it is not reasonable to pad before resampling long-term, due to the extra memory that would be used by adding this padding and resampling the padding as well.
The number of zeros in the padding depends on the size of the interpolation kernel, i.e. the padding should be half of the kernel width. Yes, you can use repeat to double the padding. Also note that vtkImageResize has an option to crop the image, which can be used to remove the padding without an additional step.
But if the goal is to replicate the scipy zoom function, perhaps you could add a similar zoom method to pyvista? Assuming that scipy zoom uses b-splines (its docs do not say), it could be duplicated almost exactly with VTK, including setting the order of the spline and having the option of turning off the prefilter.
Yes, you can use repeat to double the padding. Also note that vtkImageResize has an option to crop the image, which can be used to remove the padding without an additional step.
Thanks for the feedback! Will probably add a constant mode with this implementation for PyVista’s resample.
But if the goal is to replicate the scipy zoom function, perhaps you could add a similar zoom method to pyvista?
I am personally OK with resample as-is (my original use case for implementing it was to down-sample images), but this is a good suggestion, especially to have a different default border_mode (constant) by default. And I suppose this would also need to use VTK: vtkImageBSplineInterpolator Class Reference ? This is not (yet) an option for resample, so that would also be a key difference. Maybe @MattTheCuber is interested in contributing this kind of filter?
B-spline interpolation is achieved by first applying the prefilter to compute the coefficients, and then applying the interpolator to the output of this prefilter. The prefilter is most likely the reason that it isn’t already available as an interpolator for resample().
And speaking of b-splines, I just remembered another potential issue with adding new border modes to VTK: implementing border modes for the b-spline prefilter is hard. The papers that describe the prefilter are very math-heavy. Though hopefully if scipy is using b-splines, it might already provide an implementation.