How does imwarp interpolate data when using displacement fields?

I want to use imwarp to warp an image according to a displacement field that I specify. I run into two problems/issues: 1) The displacement field required by the MATLAB's imwarp is "inverse" according to what I would expect. 2) The interpolation of missing data using imwarp gives worse results than using a ScatteredInterpolant.

The documentation of imwarp specifies that the displacement field D should be an array of additive displacements in pixel units. It should be the same size as the input image and consist of x displacements and y-displacements. The size of D will therefore be imHeight x imWidth x 2.

The warping I want to do is a rotation of an image, but where different parts of the image is rotated with different angles, almost like Japanese folding fan. The following shows the code and a plot illustrating what I want to do:

% Make 10 lines, and rotate with angles between angle1 and angle2
nLines = 10;
angle1 = 20;
angle2 = 70;

% Create x- and y-coordinates for these lines.
x1 = (-100:100)';
y1 = ones(numel(x1), 1) * linspace(-100,100,nLines);
x1 = repmat(x1, 1, nLines);

% Plot the original lines
figure('Position', [1,1,800,400]);
s1 = subplot(1,2,1);
p1 = plot(x1, y1); 
axis equal; xlim([-110,110]); ylim([-110,110])
title('Original lines')

% Create a vector with angular offsets, different offsets for each line.
thetaOffset = linspace(angle1, angle2, nLines);

% Convert to polar coordinates to do the rotations.
[theta, rho] = cart2pol(x1, y1);
theta = theta + deg2rad(thetaOffset);
[x2, y2] = pol2cart(theta, rho);

% Plot the rotated lines
subplot(1,2,2)
p2 = plot(x2, y2);
c = get(p1, 'Color'); set(p2, {'Color'}, c)
axis equal; xlim([-110,110]); ylim([-110,110])
title('Rotated lines')

Now, I want to do the same warping with an image. I thought the obvious thing was to rotate the pixel coordinates in the same way as rotating the coordinates for the lines above, and set D to the difference between the rotated and the original pixel coordinates:

% Replot the original lines and make and image from the plot
f2=figure('Position', [1, 1, 400, 400]); axes('Position', [0,0,1,1]);
p1 = plot(x1, y1, 'LineWidth', 1); 
axis equal; xlim([-110,110]); ylim([-110,110]); axis off;

frame = getframe(f2);
exampleImage = frame2im(frame);

% Create new angular shifts, one per row in image:
thetaShift = linspace(angle1, angle2, size(exampleImage, 1))';

% Determine size and center of image
imSize = size(exampleImage);
imCenter = mean([1, imSize(1); 1, imSize(2)], 2);

% Create a grid of x and y coordinates with same gridsize as image, but
% centered on zero for before conversion to polar coordinates
[yy1, xx1] = ndgrid((1:imSize(1)) - imCenter(1), (1:imSize(2)) - imCenter(2));

% Find the polar coordinates corresponding to each pixel
[theta, rho] = cart2pol(xx1, yy1);

% Add angular shifts
thetaNew = theta + deg2rad(thetaShift);

% Convert back to cartesian coordinates
[xx2, yy2] = pol2cart(thetaNew, rho);

% Calculate shifts for each image pixel
shiftX = xx2 - xx1;
shiftY = yy2 - yy1;

% Create displacement matrix
D = cat(3, shiftX, shiftY);

% Warp image according to shifts.
[warpedImage, s] = imwarp(exampleImage, D, 'cubic');

f3 = figure('Position', [1,1,800,400]);
subplot(1,2,1); imshow(exampleImage)
subplot(1,2,2); imshow(warpedImage)

The resulting image does not look like expected:

After some digging around within the imwarp function I realized that D is not the displacement which will be added to a pixel coordinate to produce its new position, rather it is a displacement which, when added to a pixel coordinate of the target image describes which pixel coordinate of the original image should be put in that position. This distinction will hopefully be clearer through an example:

% Get grid of original pixel coordinates
[jj1, ii1] = meshgrid(1:imSize(2), 1:imSize(1));

% Find rotated pixel coordinates from the cartesian coordinates.
jj2 = xx2 + imCenter(2); % Shift right
ii2 = imSize(1) - (yy2 + imCenter(1)) + 1; % Shift up and reverse

% Find linear indices of initial (old) pixel coordinates
IND1 = sub2ind(imSize(1:2), ii1(:), jj1(:));

% Ignore pixel indices which have been shifted outside of image boundary.
valid = ii2(:) >= 1 & ii2(:) <= imSize(1) & jj2(:) >= 1 & jj2(:) <= imSize(2);

% For lack of a better idea, assign "invalid" subscripts to the upper left
% corner.
ii2(~valid) = 1;
jj2(~valid) = 1;

% Find linear indices of rotated (new) pixel coordinates
IND2 = sub2ind(imSize(1:2), round(ii2(:)), round(jj2(:)));

% Initialize Dx and Dy for the Displacement Matrix D.
Dx = zeros(imSize(1:2));
Dy = zeros(imSize(1:2));

% This is the key, shifts have to be negated and placed according 
% to rotated pixel coordinates, not original ones.
Dx(IND2) = -shiftX(IND1); % Change sign of shifts.

% Create a mask representing area outside the image borders of warped image.
imRegionMask = Dx ~= 0;
% Use imdilate followed by imerode to remove missing values within the
% displacement field.
imRegionMask = imdilate(imRegionMask, ones(3,3));
imRegionMask = imerode(imRegionMask, ones(3,3));

% Set all missing values to nan and use fillmissing to determine missing 
% shift values.
Dx(1,1) = nan;
Dx(Dx == 0) = nan;
Dx = fillmissing(Dx, 'linear'); %, 'movmean', 3);

% Do the same for Y. No need to negate shifts because image coordinates are
% reverse of cartesian coordinates.
Dy(IND2) = shiftY(IND1);
Dy(1,1) = nan;
Dy(Dy == 0) = nan;
Dy = fillmissing(Dy, 'linear'); %, 'movmean', 3);

% Put Dx and Dy into displacement matrix:
D = cat(3, Dx, Dy);

% Warp image and make sure area outside warped image is 0.
[warpedImage, s] = imwarp(exampleImage, D, 'cubic');
imRegionMask3 = repmat(imRegionMask, 1, 1, 3);
warpedImage(~imRegionMask3) = 0;

f4 = figure('Position', [1,1,800,400]);
subplot(1,2,1); imshow(exampleImage); title('Original Image')
subplot(1,2,2); imshow(warpedImage); title('Warped Image v2')

My first question is this. Is there a simpler/more elegant way for me to estimate the displacement matrix? Also, does anyone agree with me that the documentation of imwarp is easy to misunderstand?

As I was struggling to find out why imwarp did not return what I expected, I stumbled upon another advice using a scatteredInterpolant that also gives a working solution:

% Use scattered interpolant to transform image into new pixelcoordinates
% and interpolate missing values.

% I havent figured out why, but I need to flip angles upside down for this
% example.
thetaShift = flipud(thetaShift);
thetaShifted2 = theta + deg2rad(thetaShift);
[xx3, yy3] = pol2cart(thetaShifted2, rho);

% Calculate shifts for each image pixel
shiftX2 = xx3 - xx1;
shiftY2 = yy3 - yy1;

warpedImage2 = zeros(size(exampleImage), 'like', exampleImage);
ii3 = ii1 + shiftY2; jj3 = jj1 + shiftX2;
for i = 1:3
    exampleImageTmp = exampleImage(:, :, i);
    F = scatteredInterpolant(ii3(:), jj3(:), double(exampleImageTmp(:)), 'linear', 'none');
    warpedImage2(:, :, i) = uint8(F(jj1, ii1));
end

f4 = figure('Position', [1,1,800,400]);
subplot(1,2,1); imshow(exampleImage); title('Original Image')
subplot(1,2,2); imshow(warpedImage2); title('Warped Image v3')

My second question: It might not be very obvious from these images, but the last version using the scatteredInterpolant gives a smoother resulting image than the imwarp version. Does anyone know why, and is there a way to improve the interpolation used in imwarp? The scatteredInterpolant is much slower than imwarp, and also does not mask the image borders properly.

[Eivind Hennestad]

bump

[Chris Schierer]

I don't have an answer, but the coordinates in imwarp seem backwards to me too. A simple displacement field of a positive constant value for all elements translates the resulting image to the upper-left, even though the pixel indices increase to the lower-right.

Further, a displacement field of the form 1:N for all rows of D(:,:,1), and 1:M for all columns of D(:,:,2), would seem to compress the image to a single pixel (using the coordinate directions described above), but instead reduces the image by half into the upper left corner. This is not intuitive to me.

Your version/path may vary, but the underlying functions are:

C:\Program Files\MATLAB\R2019a\toolbox\images\images\+images\+geotrans\+internal\applyDisplacementField.m

which calls to do the interpolation.

C:\Program Files\MATLAB\R2019a\toolbox\images\images\+images\+internal\interp2d.m

Answer 1

[Matt J]

No, in terms of what I would expect at least, a simple left shift is giving me the right thing (see below).

I = imread('cameraman.tif');

D=ones(size(I))*50;
D=cat(3,D,0*D);  %left shift intended

Iw=imwarp(I,D);

immontage({I,Iw},'Border',[0,5] ,'Back','w')

It should be expected that the image will shift to the left when the pixel coordinates are incremented to the right, similar to the way that, given a 1D function f(t), the transformation f(t+50) will shift the graph of the function to the left.

figure;
f=@(t) t.^2;
fplot(f);hold on
fplot(@(t)f(t+50),'--');hold off
legend  Original Transformed
xlim([-200,200])