def est_nl_transf(im_ref, im_mov, fixed_image_mask='None'): # TODO: provide input parameters
    im_ref = sitk.GetImageFromArray(im_ref)
    im_mov = sitk.GetImageFromArray(im_mov)

    im_ref =  sitk.Cast(im_ref, sitk.sitkFloat32)
    im_mov = sitk.Cast(im_mov, sitk.sitkFloat32)

    # TODO: initialize the registration
    reg_method = sitk.ImageRegistrationMethod()

    # TODO: create initial identity transformation.
    transform_to_displacement_field_filter = sitk.TransformToDisplacementFieldFilter()
    transform_to_displacement_field_filter.SetReferenceImage(im_ref)
    initial_transform = sitk.DisplacementFieldTransform(transform_to_displacement_field_filter.Execute(sitk.AffineTransform(im_ref.GetDimension())))

    # TODO: regularization. The update field refers to fluid regularization; the total field to elastic regularization.
    initial_transform.SetSmoothingGaussianOnUpdate(varianceForUpdateField=2, varianceForTotalField=2)

    # TODO: set the initial transformation
    reg_method.SetInitialTransform(initial_transform)

    # TODO: use the function 'SetMetricAsDemons' to be able to perform Demons registration.
    # Be aware that you will need to provide a parameter (the intensity difference threshold) as input:
    # during the registration, intensities are considered to be equal if their difference is less than the given threshold.
    reg_method.SetMetricAsDemons(1)

    # TODO: evaluate the metrics only in the mask, if provided as an input
    if fixed_image_mask != 'None':
        reg_method.SetMetricFixedMask(fixed_image_mask)

    # TODO: set a linear interpolator
    reg_method.SetInterpolator(sitk.sitkLinear)

    # TODO: set a gradient descent optimizer
    reg_method.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=200, convergenceMinimumValue=1e-6, convergenceWindowSize=10)
    reg_method.SetOptimizerScalesFromPhysicalShift()
    print('Optimizer stop condition: {0}'.format(reg_method.GetOptimizerStopConditionDescription()))
    print('Number of iterations: {0}'.format(reg_method.GetOptimizerIteration()))
    return reg_method.Execute(im_ref, im_mov)

def apply_nl_transf(im_ref, im_mov, nl_xfm):
    """
    Apply given non-linear transform 'nl_xfm' to 'im_mov' and return the transformed image.
    """
    im_ref = sitk.GetImageFromArray(im_ref)
    im_mov = sitk.GetImageFromArray(im_mov)

    resampler = sitk.ResampleImageFilter()
    # Set the reference image
    resampler.SetReferenceImage(im_ref)
    # Use a non linear interpolator
    resampler.SetInterpolator(sitk.sitkLinear)
    # Set the desired transformationr
    resampler.SetTransform(nl_xfm)
    im_mov_resampled = resampler.Execute(im_mov)
    trans_im = sitk.GetArrayFromImage(im_mov_resampled)

    return trans_im