Python SimpleITK.sitkLinear() Examples

def produceRandomlyTranslatedImage(image, label):
    sitkImage = sitk.GetImageFromArray(image, isVector=False)
    sitklabel = sitk.GetImageFromArray(label, isVector=False)

    itemindex = np.where(label > 0)
    randTrans = (0,np.random.randint(-np.min(itemindex[1])/2,(image.shape[1]-np.max(itemindex[1]))/2),np.random.randint(-np.min(itemindex[0])/2,(image.shape[0]-np.max(itemindex[0]))/2))
    translation = sitk.TranslationTransform(3, randTrans)

    resampler = sitk.ResampleImageFilter()
    resampler.SetReferenceImage(sitkImage)
    resampler.SetInterpolator(sitk.sitkLinear)
    resampler.SetDefaultPixelValue(0)
    resampler.SetTransform(translation)

    outimgsitk = resampler.Execute(sitkImage)
    outlabsitk = resampler.Execute(sitklabel)

    outimg = sitk.GetArrayFromImage(outimgsitk)
    outimg = outimg.astype(dtype=float)

    outlbl = sitk.GetArrayFromImage(outlabsitk) > 0
    outlbl = outlbl.astype(dtype=float)

    return outimg, outlbl 

def interpolate_image(image,
                      spacing_new,
                      default_voxel_value=0):
    original_spacing = image.GetSpacing()
    original_size = image.GetSize()
    new_size = [
        round(original_size[0] * (original_spacing[0] / spacing_new[0])),
        round(original_size[1] * (original_spacing[1] / spacing_new[1])),
        round(original_size[2] * (original_spacing[2] / spacing_new[2]))
    ]
    log.debug(f'Got image with spacing {original_spacing} and size ' \
              f'{original_size}. New spacing is {spacing_new}, new size ' \
              f'is {new_size} (before padding).')
    return sitk.Resample(image, new_size, sitk.Transform(),
                         sitk.sitkLinear, image.GetOrigin(), spacing_new,
                         image.GetDirection(), default_voxel_value,
                         image.GetPixelID()) 

def __init__(self):
        self.registration_method = sitk.ImageRegistrationMethod()

        # Similarity metric settings.
        self.registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
        self.registration_method.SetMetricSamplingStrategy(self.registration_method.RANDOM)
        self.registration_method.SetMetricSamplingPercentage(0.01)

        self.registration_method.SetInterpolator(sitk.sitkLinear)

        # Optimizer settings.
        self.registration_method.SetOptimizerAsGradientDescent(learningRate=1.0, 
                                                            numberOfIterations=100, 
                                                            convergenceMinimumValue=1e-6, 
                                                            convergenceWindowSize=10)
        self.registration_method.SetOptimizerScalesFromPhysicalShift()

        # Setup for the multi-resolution framework.            
        self.registration_method.SetShrinkFactorsPerLevel(shrinkFactors = [4,2,1])
        self.registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas=[2,1,0])
        self.registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn() 

def resize_sitk_3D(self, image_array, outputSize=None, interpolator=sitk.sitkLinear):
        """
        Resample 3D images Image:
        For Labels use nearest neighbour
        For image use
        sitkNearestNeighbor = 1,
        sitkLinear = 2,
        sitkBSpline = 3,
        sitkGaussian = 4,
        sitkLabelGaussian = 5,
        """
        image = image_array
        inputSize = image.GetSize()
        inputSpacing = image.GetSpacing()
        outputSpacing = [1.0, 1.0, 1.0]
        if outputSize:
            outputSpacing[0] = inputSpacing[0] * (inputSize[0] /outputSize[0]);
            outputSpacing[1] = inputSpacing[1] * (inputSize[1] / outputSize[1]);
            outputSpacing[2] = inputSpacing[2] * (inputSize[2] / outputSize[2]);
        else:
            # If No outputSize is specified then resample to 1mm spacing
            outputSize = [0.0, 0.0, 0.0]
            outputSize[0] = int(inputSize[0] * inputSpacing[0] / outputSpacing[0] + .5)
            outputSize[1] = int(inputSize[1] * inputSpacing[1] / outputSpacing[1] + .5)
            outputSize[2] = int(inputSize[2] * inputSpacing[2] / outputSpacing[2] + .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)
        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 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_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 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 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_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 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 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_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 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 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 warp_fra(x, warp_field):
    x = apply_warp_fra(x, warp_field,
                       interpolator=sitk.sitkLinear,
                       fill_mode='constant',
                       fill_constant=0,
                       rows_idx=2, cols_idx=3)
    return x 

def warp_michal(x, warp_field):
    x = apply_warp(x, warp_field, interpolator=sitk.sitkLinear,
                   fill_mode='constant', fill_constant=0)
    return x 

def apply_warp(x, warp_field, fill_mode='reflect',
               interpolator=sitk.sitkLinear,
               fill_constant=0):
    # Expand deformation field (and later the image), padding for the largest
    # deformation
    warp_field_arr = sitk.GetArrayFromImage(warp_field)
    max_deformation = np.max(np.abs(warp_field_arr))
    pad = np.ceil(max_deformation).astype(np.int32)
    warp_field_padded_arr = pad_image(warp_field_arr, pad_amount=pad,
                                      mode='nearest')
    warp_field_padded = sitk.GetImageFromArray(warp_field_padded_arr,
                                               isVector=True)

    # Warp x, one filter slice at a time
    x_warped = np.zeros(x.shape, dtype=np.float32)
    warp_filter = sitk.WarpImageFilter()
    warp_filter.SetInterpolator(interpolator)
    warp_filter.SetEdgePaddingValue(np.min(x).astype(np.double))
    for i, image in enumerate(x):
        x_tmp = np.zeros(image.shape, dtype=image.dtype)
        for j, channel in enumerate(image):
            image_padded = pad_image(channel, pad_amount=pad, mode=fill_mode,
                                     constant=fill_constant).T
            image_f = sitk.GetImageFromArray(image_padded)
            image_f_warped = warp_filter.Execute(image_f, warp_field_padded)
            image_warped = sitk.GetArrayFromImage(image_f_warped)
            x_tmp[j] = image_warped[pad:-pad, pad:-pad].T
        x_warped[i] = x_tmp
    return x_warped 

def apply_warp(x, warp_field, fill_mode='reflect',
               interpolator=None,
               fill_constant=0, rows_idx=1, cols_idx=2):
    '''Apply an spling warp field on an image'''
    import SimpleITK as sitk
    if interpolator is None:
        interpolator = sitk.sitkLinear
    # Expand deformation field (and later the image), padding for the largest
    # deformation
    warp_field_arr = sitk.GetArrayFromImage(warp_field)
    max_deformation = np.max(np.abs(warp_field_arr))
    pad = np.ceil(max_deformation).astype(np.int32)
    warp_field_padded_arr = pad_image(warp_field_arr, pad_amount=pad,
                                      mode='nearest')
    warp_field_padded = sitk.GetImageFromArray(warp_field_padded_arr,
                                               isVector=True)

    # Warp x, one filter slice at a time
    pattern = [el for el in range(0, x.ndim) if el not in [rows_idx, cols_idx]]
    pattern += [rows_idx, cols_idx]
    inv_pattern = [pattern.index(el) for el in range(x.ndim)]
    x = x.transpose(pattern)  # batch, channel, ...
    x_shape = list(x.shape)
    x = x.reshape([-1] + x_shape[2:])  # *, r, c
    warp_filter = sitk.WarpImageFilter()
    warp_filter.SetInterpolator(interpolator)
    warp_filter.SetEdgePaddingValue(np.min(x).astype(np.double))
    for i in range(x.shape[0]):
        bc_pad = pad_image(x[i], pad_amount=pad, mode=fill_mode,
                           constant=fill_constant).T
        bc_f = sitk.GetImageFromArray(bc_pad)
        bc_f_warped = warp_filter.Execute(bc_f, warp_field_padded)
        bc_warped = sitk.GetArrayFromImage(bc_f_warped)
        x[i] = bc_warped[pad:-pad, pad:-pad].T
    x = x.reshape(x_shape)  # unsquash
    x = x.transpose(inv_pattern)
    return x 

def produceRandomlyDeformedImage(image, label, numcontrolpoints, stdDef):
    sitkImage=sitk.GetImageFromArray(image, isVector=False)
    sitklabel=sitk.GetImageFromArray(label, isVector=False)

    transfromDomainMeshSize=[numcontrolpoints]*sitkImage.GetDimension()

    tx = sitk.BSplineTransformInitializer(sitkImage,transfromDomainMeshSize)


    params = tx.GetParameters()

    paramsNp=np.asarray(params,dtype=float)
    paramsNp = paramsNp + np.random.randn(paramsNp.shape[0])*stdDef

    paramsNp[0:int(len(params)/3)]=0 #remove z deformations! The resolution in z is too bad

    params=tuple(paramsNp)
    tx.SetParameters(params)

    resampler = sitk.ResampleImageFilter()
    resampler.SetReferenceImage(sitkImage)
    resampler.SetInterpolator(sitk.sitkLinear)
    resampler.SetDefaultPixelValue(0)
    resampler.SetTransform(tx)

    resampler.SetDefaultPixelValue(0)
    outimgsitk = resampler.Execute(sitkImage)
    outlabsitk = resampler.Execute(sitklabel)

    outimg = sitk.GetArrayFromImage(outimgsitk)
    outimg = outimg.astype(dtype=np.float32)

    outlbl = sitk.GetArrayFromImage(outlabsitk)
    outlbl = (outlbl>0.5).astype(dtype=np.float32)

    return outimg,outlbl 

def getNumpyGT(self):
        dat = self.getNumpyData(self.sitkGT,sitk.sitkLinear)

        for key in dat:
            dat[key] = (dat[key]>0.5).astype(dtype=np.float32)

        return dat 

def getNumpyImages(self):
        dat = self.getNumpyData(self.sitkImages,sitk.sitkLinear)
        return dat 

def produceRandomlyDeformedImage(image, label, numcontrolpoints, stdDef, seed=1):
        '''
        This function comes from V-net，deform a image by B-spine interpolation
        :param image: images ，numpy array
        :param label: labels，numpy array
        :param numcontrolpoints: control point，B-spine interpolation parameters，take 2 for default
        :param stdDef: Deviation，B-spine interpolation parameters，take 15 for default
        :return: Deformed images and GT in numpy array
        '''
        sitkImage = sitk.GetImageFromArray(image, isVector=False)
        sitklabel = sitk.GetImageFromArray(label, isVector=False)

        transfromDomainMeshSize = [numcontrolpoints] * sitkImage.GetDimension()

        tx = sitk.BSplineTransformInitializer(
            sitkImage, transfromDomainMeshSize)

        params = tx.GetParameters()

        paramsNp = np.asarray(params, dtype=float)
        # 设置种子值，确保多通道时两个通道变换程度一样
        np.random.seed(seed)
        paramsNp = paramsNp + np.random.randn(paramsNp.shape[0]) * stdDef

        # remove z deformations! The resolution in z is too bad
        paramsNp[0:int(len(params) / 3)] = 0

        params = tuple(paramsNp)
        tx.SetParameters(params)

        resampler = sitk.ResampleImageFilter()
        resampler.SetReferenceImage(sitkImage)
        resampler.SetInterpolator(sitk.sitkLinear)
        resampler.SetDefaultPixelValue(0)
        resampler.SetTransform(tx)

        resampler.SetDefaultPixelValue(0)
        outimgsitk = resampler.Execute(sitkImage)
        outlabsitk = resampler.Execute(sitklabel)

        outimg = sitk.GetArrayFromImage(outimgsitk)
        outimg = outimg.astype(dtype=np.float32)

        outlbl = sitk.GetArrayFromImage(outlabsitk)
        # outlbl = (outlbl > 0.5).astype(dtype=np.float32)

        return outimg, outlbl 

def register_patient(self, moving_images, 
                            fixed_image, 
                            save_path,
                            save_transform=True,
                            isotropic=True):
        """
        moving_images : {'key1': path1, 'key2': path2}
        fixed_image :t1c path
        save_path: save path 
        """
        fixed_name = fixed_image.split('/').pop().split('.')[0]
        fixed_image =  sitk.ReadImage(fixed_image, sitk.sitkFloat32)
        coregistration_path = os.path.join(save_path, 'registered')
        isotropic_path = os.path.join(save_path, 'isotropic')
        transform_path = os.path.join(save_path, 'transforms')

        if not os.path.exists(coregistration_path):
            os.makedirs(coregistration_path, exist_ok=True)

        if isotropic:
            if not os.path.exists(isotropic_path):
                os.makedirs(isotropic_path, exist_ok=True)

        if save_transform:
            if not os.path.exists(transform_path):
                os.makedirs(transform_path, exist_ok=True)

        for key in moving_images.keys():
            moving_image = sitk.ReadImage(moving_images[key], sitk.sitkFloat32)
            initial_transform = sitk.CenteredTransformInitializer(fixed_image, 
                                                      moving_image, 
                                                      sitk.VersorRigid3DTransform(), 
                                                      sitk.CenteredTransformInitializerFilter.GEOMETRY)

            
            self.registration_method.SetInitialTransform(initial_transform, inPlace=False)
            final_transform = self.registration_method.Execute(sitk.Cast(fixed_image, sitk.sitkFloat32), 
                                              sitk.Cast(moving_image, sitk.sitkFloat32))
            
            print("[INFO: DeepBrainSeg] (" + strftime("%a, %d %b %Y %H:%M:%S +0000", gmtime()) + ") " +  'Final metric value: {0}'.format(self.registration_method.GetMetricValue()))
            print("[INFO: DeepBrainSeg] (" + strftime("%a, %d %b %Y %H:%M:%S +0000", gmtime()) + ") " +  'Optimizer\'s stopping condition, {0}'.format(self.registration_method.GetOptimizerStopConditionDescription()))
            
            moving_resampled= sitk.Resample(moving_image, 
                                            fixed_image, 
                                            final_transform, 
                                            sitk.sitkLinear, 0.0, 
                                            moving_image.GetPixelID())
            
            sitk.WriteImage(moving_resampled, os.path.join(coregistration_path, key+'.nii.gz'))
            sitk.WriteTransform(final_transform, os.path.join(transform_path, key+'.tfm'))
            # Write Fixed image in nii.gz
            if isotropic:
                print("[INFO: DeepBrainSeg] (" + strftime("%a, %d %b %Y %H:%M:%S +0000", gmtime()) + ") " +  'converting to isotropic volume')
                moving_resized = self.resize_sitk_3D(moving_resampled)
                sitk.WriteImage(moving_resized, os.path.join(isotropic_path, key+'.nii.gz'))

        sitk.WriteImage(fixed_image, os.path.join(coregistration_path, fixed_name+'.nii.gz'))            
        if isotropic:
            fixed_resized = self.resize_sitk_3D(fixed_image)
            sitk.WriteImage(fixed_resized, os.path.join(isotropic_path, fixed_name+'.nii.gz')) 

def _get_best_transform(self, transformations, debug=False):

        if self._refine_pca_initializations:
            transformations = self._run_registrations(transformations)

        warps = []
        for transform_sitk in transformations:
            warped_moving_sitk = sitk.Resample(
                self._moving.sitk,
                self._fixed.sitk,
                transform_sitk,
                sitk.sitkLinear,
            )
            warps.append(
                st.Stack.from_sitk_image(
                    warped_moving_sitk,
                    extract_slices=False,
                    slice_thickness=self._fixed.get_slice_thickness(),
                ))

        image_similarity_evaluator = ise.ImageSimilarityEvaluator(
            stacks=warps,
            reference=self._fixed,
            measures=[self._similarity_measure],
            use_reference_mask=True,
            verbose=False,
        )
        ph.print_info(
            "Find best aligning transform as measured by %s" %
            self._similarity_measure)
        image_similarity_evaluator.compute_similarities()
        similarities = image_similarity_evaluator.get_similarities()

        # get transform which leads to highest similarity
        index = np.argmax(similarities[self._similarity_measure])
        transform_init_sitk = transformations[index]

        if debug:
            labels = ["attempt%d" % (d + 1)
                      for d in range(len(transformations))]
            labels[index] = "best"
            foo = [w.sitk for w in warps]
            foo.insert(0, self._fixed.sitk)
            labels.insert(0, "fixed")
            sitkh.show_sitk_image(
                foo,
                segmentation=self._fixed.sitk_mask,
                label=labels,
            )
            for i in range(len(transformations)):
                print("%s: %.6f" % (
                    labels[1 + i], similarities[self._similarity_measure][i])
                )

        return transform_init_sitk 

def produceRandomlyDeformedImage(image, label, numcontrolpoints=2, stdDef=15):
    sitkImage=sitk.GetImageFromArray(image, isVector=False)
    sitklabel=sitk.GetImageFromArray(label, isVector=False)

    transfromDomainMeshSize=[numcontrolpoints]*sitkImage.GetDimension()

    tx = sitk.BSplineTransformInitializer(sitkImage,transfromDomainMeshSize)


    params = tx.GetParameters()

    paramsNp=np.asarray(params,dtype=float)
    paramsNp = paramsNp + np.random.randn(paramsNp.shape[0])*stdDef
    #remove z deformations! The resolution in z is too bad in case of 3D or its channels in 2D
    paramsNp[0:int(len(params)/3)]=0 #remove z deformations! The resolution in z is too bad in case of 3D or its channels

    params=tuple(paramsNp)
    tx.SetParameters(params)
    # print (sitkImage.GetSize(), sitklabel.GetSize(), transfromDomainMeshSize, paramsNp.shape)

    resampler = sitk.ResampleImageFilter()
    resampler.SetReferenceImage(sitkImage)
    resampler.SetInterpolator(sitk.sitkLinear)
    resampler.SetDefaultPixelValue(0)
    resampler.SetTransform(tx)

    outimgsitk = resampler.Execute(sitkImage)

    # For Label use nearest neighbour
    resampler.SetReferenceImage(sitklabel)
    resampler.SetInterpolator(sitk.sitkLabelGaussian)
    resampler.SetDefaultPixelValue(0)
    outlabsitk = resampler.Execute(sitklabel)

    outimg = sitk.GetArrayFromImage(outimgsitk)
    outimg = outimg.astype(dtype=np.float32)

    outlbl = sitk.GetArrayFromImage(outlabsitk)
    outlbl = (outlbl).astype(dtype=np.uint8)
    return outimg, outlbl

# ********************************Augmentation Transforms**************************# 

def resample(input_image,
             transform,
             output_size,
             output_spacing=None,
             output_origin=None,
             output_direction=None,
             interpolator=None,
             output_pixel_type=None,
             default_pixel_value=None):
    """
    Resample a given input image according to a transform.

    :param input_image: ITK image
        the input image
    :param transform: SimpleITK.Transform
        the (composite) transform to be applied
    :param output_size: list of int
        default is same as input image
    :param output_spacing: list of float
        default is input spacing from input_image
    :param output_direction: list of float
        default is input direction from input_image
    :param default_pixel_value:
        default is zero
    :param output_origin: list of int
        Default is zero-origin for each dimension
    :param interpolator: SimpleITK.InterpolatorEnum
        Default is SimpleITK.sitkLinear.
    :return: the resampled image
    """
    image_dim = input_image.GetDimension()
    transform_dim = transform.GetDimension()
    assert image_dim == transform_dim, 'image and transform dim must be equal, are ' + str(image_dim) + ' and ' + str(transform_dim)
    output_spacing = output_spacing or [1] * image_dim
    output_origin = output_origin or [0] * image_dim
    output_direction = output_direction or np.eye(image_dim).flatten().tolist()
    interpolator = interpolator or 'linear'

    sitk_interpolator = get_sitk_interpolator(interpolator)

    # resample the image
    resample_filter = sitk.ResampleImageFilter()
    resample_filter.SetSize(output_size)
    resample_filter.SetInterpolator(sitk_interpolator)
    resample_filter.SetOutputSpacing(output_spacing)
    resample_filter.SetOutputOrigin(output_origin)
    resample_filter.SetOutputDirection(output_direction)
    resample_filter.SetTransform(transform)
    if default_pixel_value is not None:
        resample_filter.SetDefaultPixelValue(default_pixel_value)
    if output_pixel_type is None:
        resample_filter.SetOutputPixelType(input_image.GetPixelID())
    else:
        resample_filter.SetOutputPixelType(output_pixel_type)

    # perform resampling
    output_image = resample_filter.Execute(input_image)

    return output_image