I have a 3D map of values that I need to resample, and the map is positive in one region and the rest of the image has zeros. The data format is float 64. However, I only want to resample the non-zero values and avoid including any zeros in the resampling process. My current result is this image:
The image has a standard heatmap scale with red/orange pixels representing higher values, green/blue pixels representing lower values, and the black pixels representing 0. The blue pixels on the border are what concern me. They shouldn’t be that low after resampling but they are because the resampling process includes pixels that are 0. How do I avoid this so that the resampling will ignore any zeros and only resample based on whether the pixel is positive?
The current command I’m doing is the standard SITK resampling method:
Not sure what it is that you were expecting? The linear interpolation does the interpolation across the whole image and those zero values are part of the image content.
Generally speaking, if the most important constraint on the output is not to introduce new values into the resampled image, use nearest neighbor interpolation, per @dzenanz’s advice. Note that this will likely reduce the quality/accuracy of the interpolation.
A compromise of sorts is to combine the two (has its own issues of abrupt discontinuities, but maybe it is good enough for your needs), sample code below:
import SimpleITK as sitk
original_image_size = [128,128]
img = sitk.GaussianSource(sitk.sitkUInt8, original_image_size, [16,16], [64,64])
img_size_scale = 4
ref_img = sitk.Image([img_size_scale*sz for sz in original_image_size], sitk.sitkUInt8)
ref_img.SetSpacing([spc/img_size_scale for spc in img.GetSpacing()])
# Resample using both linear and nearest neighbor, and then combine the two using nearest neighbor near the
# object boundary and linear in its interior.
resampled_img_linear = sitk.Resample(img, ref_img, sitk.Transform(), sitk.sitkLinear, 0, img.GetPixelID())
resampled_img_nn = sitk.Resample(img, ref_img, sitk.Transform(), sitk.sitkNearestNeighbor, 0, img.GetPixelID())
# Create masks for combining the two interpolation results
roi_distance_map = sitk.SignedMaurerDistanceMap(img>0, squaredDistance=False, useImageSpacing=True)
resampled_roi_distance_map = sitk.Resample(roi_distance_map, ref_img, interpolator = sitk.sitkLinear)
offset = 30.0
resampled_img_linear_mask = resampled_roi_distance_map<=-offset
resampled_img_nn_mask = resampled_roi_distance_map>-offset
sitk.Show(img, "original image")
sitk.Show(resampled_img_linear, "linear resampling")
sitk.Show(resampled_img_nn, "nearest neighbor resampling")
sitk.Show(resampled_img_linear*resampled_img_linear_mask + resampled_img_nn*resampled_img_nn_mask, "combination")
Thank you for the comments and ideas @dzenanz and @zivy. My goal is to resample using linear interpolation only within the positive values of the map, so I didn’t want to use nearest neighbor. Could you explain a bit more about what the offset is? What are the units? Is it something like this?
I think I’ve figured out that the signed Maurer distance map is calculating those distances, but I’m wondering if the offset is a threshold number for the distances and I’m not sure why it’s negative. Thanks!
“Inside” can be positive or negative, one is the default. offset in Ziv’s code is distance from the edge. You might want to change it from 30 to 1 (so only the pixels on the very edge of the object get interpolated using nearest neighbor).
