How does imwarp interpolate data when using displacement fields?

10 次查看（过去 30 天）
Eivind Hennestad
Eivind Hennestad 2018-7-30
编辑： Matt J 2023-8-25
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.
  2 个评论
Chris Schierer
Chris Schierer 2022-11-30
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

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Matt J
Matt J 2023-8-25
编辑：Matt J 2023-8-25
Ran in:
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])

类别

在 Help Center 和 File Exchange 中查找有关 Geometric Transformation and Image Registration 的更多信息

标签

产品


版本

R2017b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by