How to estimate a sub-pixel shift between images with non uniform translation/rotation

I have several images which are fully overlapping the same scene. But there is a small shift between all the images, something like 1px or less, so a sub-pixel shift. Let's say it's the problem (1): how can I estimate this sub-pixel shift between 2 images (actually, I know how and I write the code about this below). I used python here.

In addition to problem (1), there is the problem (2) which is about a non uniform shift on the full image. Let's give image A and image B, on the top left, image A is shifted about 1px from image B in x and y axis but, on the center, image A is shifted in 0.5px from image B also in x and y axis. The shift between image A and B is not uniform among the total surface of the image. The problem is how can I estimate this non-uniform shift, let's nammed it a shifted surface, for all the pixels for all images (taking one as reference) (I have a solution also for this question, I will explain it below).

Finally, problem (3) is about to shift the image with the estimated shift surface (calculated on (2)). I know how to shift an image to 0.5 px on X-axis and 1.2 px on Y-axis for example. But I don't know how to shift an array with a specific shift for each pixel.

My solutions:

Problem (1): This problem can be solved using a cross-correlation in fourrier space. A function already existe in the scipy library : register_translation reference here, I just need to give two images as parameters and the float precision I want.

Problem (2): Remember, the shift is not uniform among all the surface of the image. What I did, is basically, on window of 500x500 px, the shift is uniform and can be easily estimated from problem (1). So, I calculated the shift among all the surface of the image with window of 500x500px and a step of 100px. I have thus now, the non uniform shift estimated as shown below . I can then, interpolate a surface from this ponctual estimated shift, which will give me an estimated shift for each pixel of the image. To do that, I have to interpolate a surface with the same resolution of the image. I did it, using numpy.griddata. Here is the result for both components (x and y) . I have thus, estimated the non uniform shift among all the surface of the image.

Problem (3): I now want to apply this shift to all the image. I don't know how to do that. To shift an image at sub-pixel, you can use a function from scipy.ndimage nammed fourier_shift that you can find here but you can only give a single shift for all the image. Here, I want to give a shift for each pixel of the image.

Do you guys have any ideas to solve the problem (3) ? Also, if you think that there is a simpliest way to solve Problem 1 and 3, it can still be usefull ! For information, I have 7 images of 16000x26000px, so it take some time to solve Problem (2) as I do.

You now need to interpolate the original image at locations (x + x_shift(x,y), y + y_shift(x,y)). Likely scipy.interpolate.interpn is the most efficient way to do this.

I think your code would look something like this (not tested):

import numpy as np
import scipy

# ... (load data, find shifts, etc.)

input_coords = (np.arange(x_size), np.arange(y_size))
output_coords = np.column_stack((
    ( x_shift + input_coords[0] ).ravel(),
    ( y_shift + input_coords[1][None,:] ).ravel(),
))
output_image = scipy.interpolate.interpn(
    input_coords,
    original_image,
    output_coords,
    method='linear',
    bounds_error=False
)