Interp2 of image with transformed coordinates

Question

I have 2 greyscale images that i am trying to align using scalar scaling 1 , rotation matrix [2,2] and translation vector [2,1]. I can calculate image1's transformed coordinates as

y = s*R*x + t;

Below the resulting images are shown.

1. The first image is image1 before transformation,
2. the second image is image1 (red) with attempted interpolation using interp2 shown on top of image2 (green)
3. The third image is when i manually insert the pixel values from image1 into an empty array (that has the same size as image2) using the transformed coordinates.

From this we can see that the coordinate transformation must have been successful, as the images are aligned although not perfectly (which is to be expected since only 2 coordinates were used in calculating s, R and t) .

How come interp2 is not producing a result more similar to when i manually insert pixel values? Below the code for doing this is included:

Interpolation code

function [transformed_image] = interpolate_image(im_r,im_t,s,R,t)

[m,n] = size(im_t);

 % doesn't help if i use get_grid that the other function is using here
[~, grid_xr, grid_yr] = get_ipgrid(im_r);
[x_t, grid_xt, grid_yt] = get_ipgrid(im_t); 

y = s*R*x_t + t;
yx = reshape(y(1,:), m,n);
yy = reshape(y(2,:), m,n);

transformed_image = interp2(grid_xr, grid_yr, im_r, yx, yy, 'nearest');
end

function [x, grid_x, grid_y] = get_ipgrid(image)

[m,n] = size(image);
[grid_x,grid_y] = meshgrid(1:n,1:m);
x = [reshape(grid_x, 1, []); reshape(grid_y, 1, [])]; % X is [2xM*N] coordinate pairs
end

The manual code

function [transformed_image] = transform_image(im_r,im_t,s,R,t)

[m,n] = size(im_t);
[x_t, grid_xt, grid_yt] = get_grid(im_t);
y = s*R*x_t + t;
ymat =  reshape(y',m,n,2);
yx = ymat(:,:,1);
yy = ymat(:,:,2);

transformed_image = zeros(m,n);


for i = 1:m
    for j = 1:n
        % make sure coordinates are inside
        if (yx(i,j) < m & yy(i,j) < n & yx(i,j) > 0.5 & yy(i,j) > 0.5)
            transformed_image(round(yx(i,j)),round(yy(i,j))) = im_r(i,j);
        end
    end
end
end

function [x, grid_x, grid_y] = get_grid(image)

[m,n] = size(image);
[grid_y,grid_x] = meshgrid(1:n,1:m);
x = [grid_x(:) grid_y(:)]'; % X is [2xM*N] coordinate pairs
end

Can anyone see what i'm doing wrong with interp2? I feel like i have tried everything

Answer 1

Turns out i got interpolation all wrong.

In my question i calculate the coordinates of im1 in im2. However the way interpolation works is that i need to calculate the coordinates of im2 in im1 such that i can map the image as shown below. This means that i also calculated the wrong s,R and t since they were used to transform im1 -> im2, where as i needed im2 -> im1. (this is also called the inverse transform). Below is the manual code, that is basically the same as interp2 with nearest neighbour interpolation

function [transformed_image] = transform_image(im_r,im_t,s,R,t)
[m,n] = size(im_t);
[x_t, grid_xt, grid_yt] = get_grid(im_t);
y = s*R*x_t + t;
ymat =  reshape(y',m,n,2);
yx = ymat(:,:,1);
yy = ymat(:,:,2);

transformed_image = zeros(m,n);
for i = 1:m
    for j = 1:n
        % make sure coordinates are inside
        if (yx(i,j) < m & yy(i,j) < n & yx(i,j) > 0.5 & yy(i,j) > 0.5)
            transformed_image(i,j) = im_r(round(yx(i,j)),round(yy(i,j)));
        end
    end
end
end