import SimpleITK as sitk
import numpy as np
def calculate_origin_offset(new_spacing, old_spacing):
return np.subtract(new_spacing, old_spacing)/2
def sitk_resample_to_spacing(image, new_spacing=(1.0, 1.0, 1.0), interpolator=sitk.sitkLinear, default_value=0.):
zoom_factor = np.divide(image.GetSpacing(), new_spacing)
new_size = np.asarray(np.ceil(np.round(np.multiply(zoom_factor, image.GetSize()), decimals=5)), dtype=np.int16)
offset = calculate_origin_offset(new_spacing, image.GetSpacing())
reference_image = sitk_new_blank_image(size=new_size, spacing=new_spacing, direction=image.GetDirection(),
origin=image.GetOrigin() + offset, default_value=default_value)
return sitk_resample_to_image(image, reference_image, interpolator=interpolator, default_value=default_value)
def sitk_resample_to_image(image, reference_image, default_value=0., interpolator=sitk.sitkLinear, transform=None,
output_pixel_type=None):
if transform is None:
transform = sitk.Transform()
transform.SetIdentity()
if output_pixel_type is None:
output_pixel_type = image.GetPixelID()
resample_filter = sitk.ResampleImageFilter()
resample_filter.SetInterpolator(interpolator)
resample_filter.SetTransform(transform)
resample_filter.SetOutputPixelType(output_pixel_type)
resample_filter.SetDefaultPixelValue(default_value)
resample_filter.SetReferenceImage(reference_image)
return resample_filter.Execute(image)
def sitk_new_blank_image(size, spacing, direction, origin, default_value=0.):
image = sitk.GetImageFromArray(np.ones(size, dtype=np.float).T * default_value)
image.SetSpacing(spacing)
image.SetDirection(direction)
image.SetOrigin(origin)
return image
def resample_to_spacing(data, spacing, target_spacing, interpolation="linear", default_value=0.):
image = data_to_sitk_image(data, spacing=spacing)
if interpolation is "linear":
interpolator = sitk.sitkLinear
elif interpolation is "nearest":
interpolator = sitk.sitkNearestNeighbor
else:
raise ValueError("'interpolation' must be either 'linear' or 'nearest'. '{}' is not recognized".format(
interpolation))
resampled_image = sitk_resample_to_spacing(image, new_spacing=target_spacing, interpolator=interpolator,
default_value=default_value)
return sitk_image_to_data(resampled_image)
def data_to_sitk_image(data, spacing=(1., 1., 1.)):
if len(data.shape) == 3:
data = np.rot90(data, 1, axes=(0, 2))
image = sitk.GetImageFromArray(data)
image.SetSpacing(np.asarray(spacing, dtype=np.float))
return image
def sitk_image_to_data(image):
data = sitk.GetArrayFromImage(image)
if len(data.shape) == 3:
data = np.rot90(data, -1, axes=(0, 2))
return data
