After a bit of consideration, here is what I’ve come up with:

the function ''combine_two_images" does the trick…

Also, you need to apply obtained transform to coordinates of the moving image prior to using it.

```
from math import cos, sin
import numpy as np
def matrix_from_euler(x,y,z):
return np.array([[cos(y)*cos(z),-cos(y)*sin(z),sin(y)],
[cos(x)*sin(z) + cos(z)*sin(x)*sin(y), cos(x)*cos(z) - sin(x)*sin(y)*sin(z), -cos(y)*sin(x)],
[sin(x)*sin(z) - cos(x)*cos(z)*sin(y), cos(z)*sin(x) + cos(x)*sin(y)*sin(z), cos(x)*cos(y)]])
def euler_angles_to_matrix(arguments):
params = []
for i in range(len(arguments)):
params.append(float(arguments[i]))
angles = params[:3]
shifts = params[3:]
rotation = matrix_from_euler(*angles)
trans = translation(shifts)
M = np.zeros((4,4))
M[:3,:3] = rotation
#adding translation
for i in range(3):
M[i,3] = shifts[i]
M[3,3] = 1
return M
def apply_transform(points, M):
points = np.hstack([points, np.ones((points.shape[0],1))])
points_transformed = np.dot(M, points.T).T
return points_transformed[:,:3]
def combine_two_images(img1, img2, coords1, coords2, interp_interval):
#combine two coords arrays into one, such that points correspond to values in 'values'-array
coords_comb = np.vstack([coords1,coords2])
#combine two images into one array of values
values_comb = np.hstack([img1.ravel(),img2.ravel()])
#get min and max vals for x, y, and z across two images
x_min, x_max = (np.min(coords_comb[:,0]), np.max(coords_comb[:,0]))
y_min, y_max = (np.min(coords_comb[:,1]), np.max(coords_comb[:,1]))
z_min, z_max = (np.min(coords_comb[:,2]), np.max(coords_comb[:,2]))
n_x_points = int((x_max-x_min)/interp_interval)
n_y_points = int((y_max-y_min)/interp_interval)
n_z_points = int((z_max-z_min)/interp_interval)
grid_x, grid_y, grid_z = np.mgrid[x_min:x_max:eval(str(n_x_points) + 'j'),
y_min:y_max:eval(str(n_y_points) + 'j'),
z_min:z_max:eval(str(n_z_points) + 'j')]
grid = griddata(coords_comb, values_comb, (grid_x, grid_y, grid_z), method='nearest')
return grid
```