Python SimpleITK.sitkNearestNeighbor() Examples

def resample_image(itk_image, out_spacing=[1.0, 1.0, 1.0], is_label=False):
    original_spacing = itk_image.GetSpacing()
    original_size = itk_image.GetSize()

    out_size = [
        int(np.round(original_size[0] * (original_spacing[0] / out_spacing[0]))),
        int(np.round(original_size[1] * (original_spacing[1] / out_spacing[1]))),
        int(np.round(original_size[2] * (original_spacing[2] / out_spacing[2])))
    ]

    resample = sitk.ResampleImageFilter()
    resample.SetOutputSpacing(out_spacing)
    resample.SetSize(out_size)
    resample.SetOutputDirection(itk_image.GetDirection())
    resample.SetOutputOrigin(itk_image.GetOrigin())
    resample.SetTransform(sitk.Transform())
    resample.SetDefaultPixelValue(itk_image.GetPixelIDValue())

    if is_label:
        resample.SetInterpolator(sitk.sitkNearestNeighbor)
    else:
        resample.SetInterpolator(sitk.sitkBSpline)

    return resample.Execute(itk_image) 

def resample_image(itk_image, out_spacing=(1.0, 1.0, 1.0), is_label=False):
    original_spacing = itk_image.GetSpacing()
    original_size = itk_image.GetSize()

    out_size = [int(np.round(original_size[0] * (original_spacing[0] / out_spacing[0]))),
                int(np.round(original_size[1] * (original_spacing[1] / out_spacing[1]))),
                int(np.round(original_size[2] * (original_spacing[2] / out_spacing[2])))]

    resample = sitk.ResampleImageFilter()
    resample.SetOutputSpacing(out_spacing)
    resample.SetSize(out_size)
    resample.SetOutputDirection(itk_image.GetDirection())
    resample.SetOutputOrigin(itk_image.GetOrigin())
    resample.SetTransform(sitk.Transform())
    resample.SetDefaultPixelValue(itk_image.GetPixelIDValue())

    if is_label:
        resample.SetInterpolator(sitk.sitkNearestNeighbor)
    else:
        resample.SetInterpolator(sitk.sitkBSpline)

    return resample.Execute(itk_image) 

def PostPadding(self, seg_post_3d, postProcessList, max_size=(256, 256)):
        """
        Handle : Resizing or post padding operations to get back image to original shape
        """
        # Resize and Pad the output 3d volume to its original dimension 
        # If the ROI extraction resized the image then upsample to the original shape
        boResized = self.input_img['resized']
        if boResized:
            # Upsample the image to max_size = (256, 256)
            seg_resized = np.zeros((seg_post_3d.shape[0], max_size[0], max_size[1]), dtype=np.uint8)
            for slice in range(seg_post_3d.shape[0]):
                seg_resized[slice] = resize_sitk_2D(seg_post_3d[slice], max_size, interpolator=sitk.sitkNearestNeighbor) 
            seg_post_3d = seg_resized

        if 'pad_patch' in postProcessList:
            seg_post_3d = pad_3Dpatch(seg_post_3d, self.input_img['pad_params'])
        # print (seg_post_3d.shape)
        return swapaxes_to_xyz(seg_post_3d) 

def _add_slice_contribution(slice, coverage_sitk):

        #
        slice_sitk = sitk.Image(slice.sitk)
        spacing = np.array(slice_sitk.GetSpacing())
        spacing[-1] = slice.get_slice_thickness()
        slice_sitk.SetSpacing(spacing)

        contrib_nda = sitk.GetArrayFromImage(slice_sitk)
        contrib_nda[:] = 1
        contrib_sitk = sitk.GetImageFromArray(contrib_nda)
        contrib_sitk.CopyInformation(slice_sitk)

        coverage_sitk += sitk.Resample(
            contrib_sitk,
            coverage_sitk,
            sitk.Euler3DTransform(),
            sitk.sitkNearestNeighbor,
            0,
            coverage_sitk.GetPixelIDValue(),
        )

        return coverage_sitk 

def reorient_image(image):
    """Reorients an image to standard radiology view."""

    dir = np.array(image.GetDirection()).reshape(len(image.GetSize()), -1)
    ind = np.argmax(np.abs(dir), axis=0)
    new_size = np.array(image.GetSize())[ind]
    new_spacing = np.array(image.GetSpacing())[ind]
    new_extent = new_size * new_spacing
    new_dir = dir[:, ind]

    flip = np.diag(new_dir) < 0
    flip_diag = flip * -1
    flip_diag[flip_diag == 0] = 1
    flip_mat = np.diag(flip_diag)

    new_origin = np.array(image.GetOrigin()) + np.matmul(new_dir, (new_extent * flip))
    new_dir = np.matmul(new_dir, flip_mat)

    resample = sitk.ResampleImageFilter()
    resample.SetOutputSpacing(new_spacing.tolist())
    resample.SetSize(new_size.tolist())
    resample.SetOutputDirection(new_dir.flatten().tolist())
    resample.SetOutputOrigin(new_origin.tolist())
    resample.SetTransform(sitk.Transform())
    resample.SetDefaultPixelValue(image.GetPixelIDValue())
    resample.SetInterpolator(sitk.sitkNearestNeighbor)

    return resample.Execute(image) 

def run(self):

        recon_space = self._target.get_isotropically_resampled_stack(
            extra_frame=self._max_distance,
        )
        mask_sitk = 0 * recon_space.sitk_mask
        dim = mask_sitk.GetDimension()

        for stack in self._stacks:
            stack_mask_sitk = sitk.Resample(
                stack.sitk_mask,
                mask_sitk,
                eval("sitk.Euler%dDTransform()" % dim),
                sitk.sitkNearestNeighbor,
                0,
                mask_sitk.GetPixelIDValue())
            mask_sitk += stack_mask_sitk

        thresholder = sitk.BinaryThresholdImageFilter()
        mask_sitk = thresholder.Execute(mask_sitk, 0, 0.5, 0, 1)

        if self._dilation_radius > 0:
            dilater = sitk.BinaryDilateImageFilter()
            dilater.SetKernelType(eval("sitk.sitk" + self._dilation_kernel))
            dilater.SetKernelRadius(self._dilation_radius)
            mask_sitk = dilater.Execute(mask_sitk)

        self._joint_image_mask = st.Stack.from_sitk_image(
            image_sitk=recon_space.sitk,
            image_sitk_mask=mask_sitk,
            filename=self._target.get_filename(),
            slice_thickness=recon_space.get_slice_thickness(),
        ) 

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 get_sitk_interpolator(interpolator):
    if interpolator == 'nearest':
        return sitk.sitkNearestNeighbor
    elif interpolator == 'linear':
        return sitk.sitkLinear
    elif interpolator == 'cubic':
        return sitk.sitkBSpline
    elif interpolator == 'label_gaussian':
        return sitk.sitkLabelGaussian
    elif interpolator == 'gaussian':
        return sitk.sitkGaussian
    elif interpolator == 'lanczos':
        return sitk.sitkLanczosWindowedSinc
    else:
        raise Exception('invalid interpolator type') 

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 Resampling(Image,Label):
    Size=Image.GetSize()
    Spacing=Image.GetSpacing()
    Origin = Image.GetOrigin()
    Direction = Image.GetDirection()
    ImagePyramid=[]
    LabelPyramid=[]
    for i in range(3):
        NewSpacing = ToSpacing[ResRate[i]]
        NewSize=[int(Size[0]*Spacing[0]/NewSpacing[0]),int(Size[1]*Spacing[1]/NewSpacing[1]),int(Size[2]*Spacing[2]/NewSpacing[2])]       
        Resample = sitk.ResampleImageFilter()
        Resample.SetOutputDirection(Direction)
        Resample.SetOutputOrigin(Origin)
        Resample.SetSize(NewSize)
        Resample.SetInterpolator(sitk.sitkLinear)
        Resample.SetOutputSpacing(NewSpacing)
        NewImage = Resample.Execute(Image)
        ImagePyramid.append(NewImage)
        
        Resample = sitk.ResampleImageFilter()
        Resample.SetOutputDirection(Direction)
        Resample.SetOutputOrigin(Origin)
        Resample.SetSize(NewSize)
        Resample.SetOutputSpacing(NewSpacing)
        Resample.SetInterpolator(sitk.sitkNearestNeighbor)
        NewLabel = Resample.Execute(Label)
        LabelPyramid.append(NewLabel)
    return ImagePyramid,LabelPyramid

#We shift the mean value to enhance the darker side 

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 resample_image(image, out_spacing=(1.0, 1.0, 1.0), out_size=None, is_label=False, pad_value=0):
    """Resamples an image to given element spacing and output size."""

    original_spacing = np.array(image.GetSpacing())
    original_size = np.array(image.GetSize())

    if out_size is None:
        out_size = np.round(np.array(original_size * original_spacing / np.array(out_spacing))).astype(int)
    else:
        out_size = np.array(out_size)

    original_direction = np.array(image.GetDirection()).reshape(len(original_spacing),-1)
    original_center = (np.array(original_size, dtype=float) - 1.0) / 2.0 * original_spacing
    out_center = (np.array(out_size, dtype=float) - 1.0) / 2.0 * np.array(out_spacing)

    original_center = np.matmul(original_direction, original_center)
    out_center = np.matmul(original_direction, out_center)
    out_origin = np.array(image.GetOrigin()) + (original_center - out_center)

    resample = sitk.ResampleImageFilter()
    resample.SetOutputSpacing(out_spacing)
    resample.SetSize(out_size.tolist())
    resample.SetOutputDirection(image.GetDirection())
    resample.SetOutputOrigin(out_origin.tolist())
    resample.SetTransform(sitk.Transform())
    resample.SetDefaultPixelValue(pad_value)

    if is_label:
        resample.SetInterpolator(sitk.sitkNearestNeighbor)
    else:
        #resample.SetInterpolator(sitk.sitkBSpline)
        resample.SetInterpolator(sitk.sitkLinear)

    return resample.Execute(sitk.Cast(image, sitk.sitkFloat32)) 

def resample_image_to_ref(image, ref, is_label=False, pad_value=0):
    """Resamples an image to match the resolution and size of a given reference image."""

    resample = sitk.ResampleImageFilter()
    resample.SetReferenceImage(ref)
    resample.SetDefaultPixelValue(pad_value)

    if is_label:
        resample.SetInterpolator(sitk.sitkNearestNeighbor)
    else:
        #resample.SetInterpolator(sitk.sitkBSpline)
        resample.SetInterpolator(sitk.sitkLinear)

    return resample.Execute(image) 

def run_bias_field_correction(self):

        time_start = ph.start_timing()

        bias_field_corrector = sitk.N4BiasFieldCorrectionImageFilter()

        bias_field_corrector.SetBiasFieldFullWidthAtHalfMaximum(
            self._bias_field_fwhm)
        bias_field_corrector.SetConvergenceThreshold(
            self._convergence_threshold)
        bias_field_corrector.SetSplineOrder(self._spline_order)
        bias_field_corrector.SetWienerFilterNoise(self._wiener_filter_noise)

        if self._use_mask:
            image_sitk = bias_field_corrector.Execute(
                self._stack.sitk, self._stack.sitk_mask)
        else:
            image_sitk = bias_field_corrector.Execute(self._stack.sitk)

        # Reading of image might lead to slight differences
        stack_corrected_sitk_mask = sitk.Resample(
            self._stack.sitk_mask,
            image_sitk,
            sitk.Euler3DTransform(),
            sitk.sitkNearestNeighbor,
            0,
            self._stack.sitk_mask.GetPixelIDValue())

        self._stack_corrected = st.Stack.from_sitk_image(
            image_sitk=image_sitk,
            image_sitk_mask=stack_corrected_sitk_mask,
            filename=self._prefix_corrected + self._stack.get_filename(),
            slice_thickness=self._stack.get_slice_thickness(),
        )

        # Get computational time
        self._computational_time = ph.stop_timing(time_start)

        # Debug
        # sitkh.show_stacks([self._stack, self._stack_corrected], label=["orig", "corr"]) 

def get_warped_moving(self):

        warped_moving_sitk_mask = sitk.Resample(
            self._moving.sitk_mask,
            self._fixed.sitk,
            self.get_registration_transform_sitk(),
            sitk.sitkNearestNeighbor,
            0,
            self._moving.sitk_mask.GetPixelIDValue(),
        )

        if isinstance(self._moving, st.Stack):
            warped_moving = st.Stack.from_sitk_image(
                image_sitk=self._get_warped_moving_sitk(),
                filename=self._moving.get_filename(),
                image_sitk_mask=warped_moving_sitk_mask,
                slice_thickness=self._fixed.get_slice_thickness(),
            )
        else:
            warped_moving = sl.Slice.from_sitk_image(
                image_sitk=self._get_warped_moving_sitk(),
                filename=self._moving.get_filename(),
                image_sitk_mask=warped_moving_sitk_mask,
                slice_thickness=self._fixed.get_slice_thickness(),
            )
        return warped_moving

    ##
    # Gets the fixed image transformed by the obtained registration transform.
    #
    # The returned image will align the fixed image with the moving image as
    # found during the registration.
    # \date       2017-08-08 16:53:21+0100
    #
    # \param      self  The object
    #
    # \return     The transformed fixed as Stack/Slice object
    # 

def resize_sitk_2D(image_array, outputSize=None, interpolator=sitk.sitkLinear):
    """
    Resample 2D images Image:
    For Labels use nearest neighbour
    For image use 
    sitkNearestNeighbor = 1,
    sitkLinear = 2,
    sitkBSpline = 3,
    sitkGaussian = 4,
    sitkLabelGaussian = 5, 
    """
    image = sitk.GetImageFromArray(image_array) 
    inputSize = image.GetSize()
    inputSpacing = image.GetSpacing()
    outputSpacing = [1.0, 1.0]
    if outputSize:
        outputSpacing[0] = inputSpacing[0] * (inputSize[0] /outputSize[0]);
        outputSpacing[1] = inputSpacing[1] * (inputSize[1] / outputSize[1]);
    else:
        # If No outputSize is specified then resample to 1mm spacing
        outputSize = [0.0, 0.0]
        outputSize[0] = int(inputSize[0] * inputSpacing[0] / outputSpacing[0] + .5)
        outputSize[1] = int(inputSize[1] * inputSpacing[1] / outputSpacing[1] + .5)
    resampler = sitk.ResampleImageFilter()
    resampler.SetSize(outputSize)
    resampler.SetOutputSpacing(outputSpacing)
    resampler.SetOutputOrigin(image.GetOrigin())
    resampler.SetOutputDirection(image.GetDirection())
    resampler.SetInterpolator(interpolator)
    resampler.SetDefaultPixelValue(0)
    image = resampler.Execute(image)
    resampled_arr = sitk.GetArrayFromImage(image)
    return resampled_arr 

def reslice_image(itk_image, itk_ref, is_label=False):
    resample = sitk.ResampleImageFilter()
    resample.SetReferenceImage(itk_ref)

    if is_label:
        resample.SetInterpolator(sitk.sitkNearestNeighbor)
    else:
        resample.SetInterpolator(sitk.sitkBSpline)

    return resample.Execute(itk_image) 

def test_ground_truth_affine_transforms_no_motion(self):

        filename_HR_volume = "HR_volume_postmortem"
        HR_volume = st.Stack.from_filename(
            os.path.join(self.dir_test_data, filename_HR_volume + ".nii.gz"),
            os.path.join(self.dir_test_data,
                         filename_HR_volume + "_mask.nii.gz")
        )

        # Run simulation for Nearest Neighbor interpolation (shouldn't not
        # matter anyway and is quicker)
        slice_acquistion = sa.SliceAcqusition(HR_volume)
        slice_acquistion.set_interpolator_type("NearestNeighbor")
        slice_acquistion.set_motion_type("NoMotion")

        slice_acquistion.run_simulation_view_1()
        slice_acquistion.run_simulation_view_2()
        slice_acquistion.run_simulation_view_3()

        # Get simulated stack of slices + corresponding ground truth affine
        #  transforms indicating the correct acquisition of the slice
        #  within the volume
        stacks_simulated = slice_acquistion.get_simulated_stacks()
        affine_transforms_ground_truth, rigid_motion_transforms_ground_truth = slice_acquistion.get_ground_truth_data()

        for i in range(0, len(stacks_simulated)):
            stack = st.Stack.from_stack(stacks_simulated[i])

            slices = stack.get_slices()
            N_slices = stack.get_number_of_slices()

            for j in range(0, N_slices):
                # print("Stack %s: Slice %s/%s" %(i,j,N_slices-1))
                slice = slices[j]
                # slice.update_motion_correction(rigid_motion_transforms_ground_truth[i][j])

                HR_volume_resampled_slice_sitk = sitk.Resample(
                    HR_volume.sitk, slice.sitk, sitk.Euler3DTransform(
                    ), sitk.sitkNearestNeighbor, 0.0, slice.sitk.GetPixelIDValue()
                )

                self.assertEqual(np.round(
                    np.linalg.norm(sitk.GetArrayFromImage(slice.sitk - HR_volume_resampled_slice_sitk)), decimals=self.accuracy), 0)

    # Test whether the ground truth affine transforms set during the
    #  simulation correspond to the actually acquired positions within the
    #  sliced volume whereby motion is applied to the HR volume 

def test_ground_truth_affine_transforms_with_motion_NearestNeighbor(self):

        filename_HR_volume = "HR_volume_postmortem"
        HR_volume = st.Stack.from_filename(
            os.path.join(self.dir_test_data, filename_HR_volume + ".nii.gz"),
            os.path.join(self.dir_test_data,
                         filename_HR_volume + "_mask.nii.gz")
        )

        # Run simulation for Nearest Neighbor interpolation (shouldn't not
        # matter anyway and is quicker)
        slice_acquistion = sa.SliceAcqusition(HR_volume)
        slice_acquistion.set_interpolator_type("NearestNeighbor")
        slice_acquistion.set_motion_type("Random")

        slice_acquistion.run_simulation_view_1()
        slice_acquistion.run_simulation_view_2()
        slice_acquistion.run_simulation_view_3()

        # Get simulated stack of slices + corresponding ground truth affine
        #  transforms indicating the correct acquisition of the slice
        #  within the volume
        stacks_simulated = slice_acquistion.get_simulated_stacks()
        affine_transforms_ground_truth, rigid_motion_transforms_ground_truth = slice_acquistion.get_ground_truth_data()

        for i in range(0, len(stacks_simulated)):
            stack = stacks_simulated[i]

            slices = stack.get_slices()
            N_slices = stack.get_number_of_slices()

            for j in range(0, N_slices):
                # print("Stack %s: Slice %s/%s" %(i,j,N_slices-1))
                slice = slices[j]
                slice.update_motion_correction(
                    rigid_motion_transforms_ground_truth[i][j])

                HR_volume_resampled_slice_sitk = sitk.Resample(
                    HR_volume.sitk, slice.sitk, sitk.Euler3DTransform(
                    ), sitk.sitkNearestNeighbor, 0.0, slice.sitk.GetPixelIDValue()
                )

                self.assertEqual(np.round(
                    np.linalg.norm(sitk.GetArrayFromImage(slice.sitk - HR_volume_resampled_slice_sitk)), decimals=self.accuracy), 0)

    # Test whether the ground truth affine transforms set during the
    #  simulation correspond to the actually acquired positions within the
    #  sliced volume whereby motion is applied to the HR volume 

def generate_mask_from_stack_mask_intersections(self,
                                                    mask_dilation_radius=0,
                                                    mask_dilation_kernel="Ball",
                                                    ):

        # Define helpers to obtain averaged stack
        shape = sitk.GetArrayFromImage(self._HR_volume.sitk).shape
        array_mask = np.ones(shape, dtype=np.uint8)

        # Average over domain specified by the joint mask ("union mask")
        for i in range(0, self._N_stacks):

            # Resample warped stack masks
            stack_sitk_mask = sitk.Resample(
                self._stacks[i].sitk_mask,
                self._HR_volume.sitk_mask,
                sitk.Euler3DTransform(),
                sitk.sitkNearestNeighbor,
                0,
                self._HR_volume.sitk_mask.GetPixelIDValue())

            # Get arrays of resampled warped stack and mask
            array_mask_tmp = sitk.GetArrayFromImage(
                stack_sitk_mask).astype(np.uint8)

            # Sum intensities of stack and mask
            array_mask *= array_mask_tmp

        # Create (joint) binary mask. Mask represents union of all masks
        array_mask[array_mask > 0] = 1

        HR_volume_mask_sitk = sitk.GetImageFromArray(array_mask)
        HR_volume_mask_sitk.CopyInformation(self._HR_volume.sitk)

        if mask_dilation_radius > 0:
            dilater = sitk.BinaryDilateImageFilter()
            dilater.SetKernelType(eval("sitk.sitk" + mask_dilation_kernel))
            dilater.SetKernelRadius(mask_dilation_radius)
            HR_volume_mask_sitk = dilater.Execute(HR_volume_mask_sitk)

        self._HR_volume = st.Stack.from_sitk_image(
            image_sitk=self._HR_volume.sitk,
            filename=self._HR_volume.get_filename(),
            image_sitk_mask=HR_volume_mask_sitk,
            slice_thickness=self._HR_volume.get_slice_thickness(),
        ) 

def _get_jacobian_slice_in_slice_neighbours_fit(self,
                                                    slice,
                                                    transform_sitk,
                                                    transform_itk):

        # Get slice(T(theta, x))
        slice_sitk = sitk.Resample(
            slice.sitk,
            self._slice_grid_2D_sitk,
            transform_sitk,
            self._interpolator_sitk)

        # Get d[slice(T(theta, x))]/dx as (Ny x Nx x dim)-array
        dslice_nda = self._get_gradient_image_nda_from_sitk_image(slice_sitk)

        # Get slice data array (used for intensity correction parameter
        # gradient)
        slice_nda = sitk.GetArrayFromImage(slice_sitk)

        if self._use_stack_mask_neighbour_fit_term:
            slice_sitk_mask = sitk.Resample(
                slice.sitk_mask,
                self._slice_grid_2D_sitk,
                transform_sitk,
                sitk.sitkNearestNeighbor)
            slice_nda_mask = sitk.GetArrayFromImage(slice_sitk_mask)

            # slice_nda *= slice_nda_mask[:,:,np.newaxis]
            slice_nda *= slice_nda_mask

        # Get Jacobian of slice w.r.t to transform parameters
        jacobian_slice_nda = \
            self._get_gradient_with_respect_to_transform_parameters(
                dslice_nda, transform_itk, slice_sitk)

        # Add Jacobian w.r.t. to intensity correction parameters
        jacobian_slice_nda = \
            self._add_gradient_with_respect_to_intensity_correction_parameters[
                self._intensity_correction_type_slice_neighbour_fit](
                    jacobian_slice_nda, slice_nda)

        return jacobian_slice_nda

    ##
    # Gets the gradient with respect to transform parameters of all voxels
    # within a slice.
    # \date       2017-07-15 23:03:10+0100
    #
    # \param      self           The object
    # \param      dslice_nda     The dslice nda
    # \param      transform_itk  The transform itk
    # \param      slice_sitk     The slice sitk
    #
    # \return     The gradient with respect to transform parameters;
    #             (N_slice_voxels x transform_type_dofs) numpy array
    # 

def get_increased_stack(self, extra_slices_z=0):

        interpolator = sitk.sitkNearestNeighbor

        # Read original spacing (voxel dimension) and size of target stack:
        spacing = np.array(self.sitk.GetSpacing())
        size = np.array(self.sitk.GetSize()).astype("int")
        origin = np.array(self.sitk.GetOrigin())
        direction = self.sitk.GetDirection()

        # Update information according to isotropic resolution
        size[2] += extra_slices_z

        # Resample image and its mask to isotropic grid
        default_pixel_value = 0.0

        isotropic_resampled_stack_sitk = sitk.Resample(
            self.sitk,
            size,
            sitk.Euler3DTransform(),
            interpolator,
            origin,
            spacing,
            direction,
            default_pixel_value,
            self.sitk.GetPixelIDValue())

        isotropic_resampled_stack_sitk_mask = sitk.Resample(
            self.sitk_mask,
            size,
            sitk.Euler3DTransform(),
            sitk.sitkNearestNeighbor,
            origin,
            spacing,
            direction,
            0,
            self.sitk_mask.GetPixelIDValue())

        # Create Stack instance
        stack = self.from_sitk_image(
            isotropic_resampled_stack_sitk, "zincreased_" + self._filename, isotropic_resampled_stack_sitk_mask)

        return stack 

def get_isotropically_resampled_stack(self, resolution=None, interpolator="Linear", extra_frame=0, filename=None, mask_dilation_radius=0, mask_dilation_kernel="Ball"):

        # Choose interpolator
        try:
            interpolator_str = interpolator
            interpolator = eval("sitk.sitk" + interpolator_str)
        except:
            raise ValueError("Error: interpolator is not known")

        if resolution is None:
            spacing = self.sitk.GetSpacing()[0]
        else:
            spacing = resolution

        # Resample image and its mask to isotropic grid
        resampler = simplereg.resampler.Resampler
        isotropic_resampled_stack_sitk = resampler.get_resampled_image_sitk(
            image_sitk=self.sitk,
            spacing=spacing,
            interpolator=interpolator,
            padding=0.0,
            add_to_grid=extra_frame,
            add_to_grid_unit="mm",
        )
        isotropic_resampled_stack_sitk_mask = resampler.get_resampled_image_sitk(
            image_sitk=self.sitk_mask,
            spacing=spacing,
            interpolator=sitk.sitkNearestNeighbor,
            padding=0,
            add_to_grid=extra_frame,
            add_to_grid_unit="mm",
        )

        if mask_dilation_radius > 0:
            dilater = sitk.BinaryDilateImageFilter()
            dilater.SetKernelType(eval("sitk.sitk" + mask_dilation_kernel))
            dilater.SetKernelRadius(mask_dilation_radius)
            isotropic_resampled_stack_sitk_mask = dilater.Execute(
                isotropic_resampled_stack_sitk_mask)

        # Create Stack instance
        if filename is None:
            filename = self._filename + "_" + interpolator_str + "Iso"
        stack = self.from_sitk_image(
            image_sitk=isotropic_resampled_stack_sitk,
            filename=filename,
            slice_thickness=isotropic_resampled_stack_sitk.GetSpacing()[-1],
            image_sitk_mask=isotropic_resampled_stack_sitk_mask,
        )

        return stack

    # Increase stack by adding zero voxels in respective directions
    #  \remark Used for MS project to add empty slices on top of (chopped) brain
    #  \param[in] resolution length of voxel side, scalar
    #  \param[in] interpolator choose type of interpolator for resampling
    #  \param[in] extra_frame additional extra frame of zero intensities surrounding the stack in mm
    #  \return isotropically, resampled stack as Stack object 

def reslice_image(itk_image, itk_ref, is_label=False):
    resample = sitk.ResampleImageFilter()
    resample.SetReferenceImage(itk_ref)

    if is_label:
        resample.SetInterpolator(sitk.sitkNearestNeighbor)
    else:
        resample.SetInterpolator(sitk.sitkBSpline)

    return resample.Execute(itk_image)