Strange (wrong?) interpolator behavior

Algorithms
,
1

Hi all,

I was trying various interpolators to use applying a transformation I got from a simple euclid 2d registration. The image I’m applying it to is a CT image, and I’m using a 64bit float image type.

When I resample the image with my transformation, for at least a couple of the interpolators other than linear (Bspline and CoswindowSinc at least), the min and max of the pixel values after re-sampling is less than/greater than what it was prior to.

So the interpolation is resulting in pixel values outside the min/max bounds of the previous image pixels.

IMHO this seems like it shouldn’t be doing that - interpolating a new pixel value shouldn’t give you a value that below/above any of the neighbouring pixels regardless of which algorithm it’s using. …you’re meant to be interpolating between values right???

Is this the correct behavior???

Here is a snippet of the code I’m using where it outputs this behavior, plus an output image from it where it clearly shows the min/max values of the resampled (translated) image have change a lot:

# want to translate the dicom scaled_pixel_array by the transform returned
# want the default pixel value to be take from the average in a corner of
# the array (to match the HU value of the background)

min_value = sitk.GetArrayViewFromImage(dcm_sitk).min()
max_value = sitk.GetArrayViewFromImage(dcm_sitk).max()
print('min/max pre transform', min_value, max_value)
plt.imshow(sitk.GetArrayViewFromImage(dcm_sitk))
plt.show()

# get background values from a corner average
background = np.average(sitk.GetArrayViewFromImage(dcm_sitk)[0:10,0:10])
print('background:', background)

dcm_sitk = sitk.Resample(dcm_sitk, dcm_sitk, transform,
        sitk.sitkBSpline, background)
        #sitk.sitkLinear, background)
        #sitk.sitkCosineWindowedSinc, background)

no_below = (sitk.GetArrayViewFromImage(dcm_sitk) < min_value).sum()
no_above = (sitk.GetArrayViewFromImage(dcm_sitk) > max_value).sum()

new_min_value = sitk.GetArrayViewFromImage(dcm_sitk).min()
new_max_value = sitk.GetArrayViewFromImage(dcm_sitk).max()

print('post transform', new_min_value, new_max_value)
print('number below min: ', no_below)
print('number above max: ', no_above)

plt.imshow(sitk.GetArrayViewFromImage(dcm_sitk))
plt.show()

below_img = (sitk.GetArrayViewFromImage(dcm_sitk) < min_value)
above_img = (sitk.GetArrayViewFromImage(dcm_sitk) > max_value)
plt.imshow(below_img)
plt.show()
plt.imshow(above_img)
plt.show()

Output:

min/max pre transform -3024.0 3071.0
background: -3024.0
post transform -3415.37375494161 3255.1812187911864
number below min:  18477
number above max:  9

Figure_4
Figure_41280×960 9.15 KB
Figure_3
Figure_31280×960 18.5 KB
Figure_2
Figure_21280×960 85 KB
Figure_1
Figure_11280×960 76.4 KB

2

Hi,

I think you are mixing the parameter space and the value space in your mind. Interpolating means finding the value of a data point which is between two other known in parameter space. But depending on the interpolation method you could find values out of the currently observed range. That often the case with spline. See some graphs examples on the matlab doc for spline interpolation.
Regarding your specific case, you interpolate from spatial position (parameter space) to intensity (value space), so your result doesn’t seem wrong to me :slight_smile:

HTH

1 Like
3

When you fit a spline to points, some parts of the spline can, and usually go above the highest point. Second graph from Tim’s link draws it very neatly. Resampled pixels will fall on different points of the spline, and frequently some will fall in the region of the spline which is higher than any of the original points.

This is how spline interpolation works. If you care about your minimum/maximum values, you could calculate them after the resampling. If you want to preserve the min-max range, you could use rescale filter or shift scale filter.

2 Likes