Computing coordinates of points of an image after elastic deformation

Question

My task is: given an image and set of points of interest, elastically and randomly deform the image and save it with the modified aforementioned points.

example: (blue points are the points of interest: POIs)

with coordinates:

points = np.array([
    [442, 131],
    [438, 124],
    [438, 136],
    [431, 118],
    [432, 143],
    [413, 113],
    [413, 147],
])

Then i am initializing random deformation field and applying to the image using scipy and numpy packages:

import numpy as np
from scipy.ndimage import map_coordinates, gaussian_filter
from scipy.interpolate import interpn

atlas, points = load_atlas() # local function returns grayscale image and points as two matrices
shape = atlas.shape # returns (512, 256)
rnd_state = np.random.RandomState(3422347) # just some seed
alpha = 512 * 30 # constant
delta = gaussian_filter(input=(rnd_state.rand(*shape) * 2 - 1) * alpha, sigma=42)

# initialize additional elongation along vertical dimension principle is the same as above
elongation = random_elongation(shape)

y, x = np.meshgrid(np.arange(512), np.arange(256))
# apply movement of points and extract them as indices
indices = np.array([
    np.reshape(y + delta + elongation, (-1, 1)),
    np.reshape(x + delta, (-1, 1))
])

# interpolate using scipy.ndimage.map_coordinates
atlas_deformed = map_coordinates(atlas, indices, mode='constant', order=3).reshape(shape)

# 'apply' deformation on points
deformed_y, deformed_x = y - delta - elongation, x - delta

# extract only points of interest from meshgrids
deformed_points = np.zeros(points.shape)
i = 0
for point in points:
    coord_y, coord_x = int(point[0]), int(point[1])
    __y = deformed_y[coord_y, coord_x]
    __x = deformed_x[coord_y, coord_x]
    deformed_points[i] = [__y, __x]
    i += 1

Then i get end-result with various success e.g. in 50% of cases everything works as intended:

but sometimes the points seem to be completely off:

My guess is that I compute points by directly subtracting from initial coordinates while image is interpolated using cubic spline. Is it possible to get the coordinates of the points after spline interpolation?