Why imwarp isn't working for cylindrical projection?

I am trying to find the cylindrical projection of an image using the cylindrical projection equations. After finding the warpped values for the X and Y corodinates, I am trying to use the inbuilt MATLAB imwarp function to project the image onto the new coordinates, so that I can use the interpolation for smooth image. However, I get a different warpped image. Here is the output what I am getting:

What I expect the image to be:

Below is the cylindrical projection code:

clc; close all; clear;

% Inputs
fileName = 'checker.jpg';

% Focal lengths
fx = 50;
fy = 50;

% Read image
image = (imread(fileName));

% Get image size
[ydim, xdim, bypixs] = size(image);

% Camera intrinsics
K = [fx, 0, xdim/2; 0, fy, ydim/2; 0, 0, 1];

% Distortion coefficients [k1, k2, k3, p1, p2]
DC = [0, 0, 0, 0, 0];

% Get distrotion coefficients
fx = K(1,1);
fy = K(2,2);
k1 = DC(1);
k2 = DC(2);
k3 = DC(3);
p1 = DC(4);
p2 = DC(5);

% Get image size
[ydim, xdim, bypixs] = size(image);

% Get the center of image
xc = xdim/2;
yc = ydim/2;

% Create X and Y coordinates grid
[X,Y] = meshgrid(1:xdim, 1:ydim);

% Perform the cylindrical projection
theta = (X - xc) / fx;
h = (Y - yc) / fy;

% Cylindrical coordinates to Cartesian
xcap = sin(theta);
ycap = h;
zcap = cos(theta);

xyz_cap = cat(3, xcap, ycap, zcap);
xyz_cap = reshape(xyz_cap,[],3);

% Normalized coords
xyz_cap_norm = (K * xyz_cap')';
xn = xyz_cap_norm(:,1) ./ xyz_cap_norm(:,3);
yn = xyz_cap_norm(:,2) ./ xyz_cap_norm(:,3);

% Radial and tangential distortion 
r = xn.^2 + yn.^2;
xd_r = xn .* (1 + k1 * r.^2 + k2 * r.^4 + k3 * r.^6);
yd_r = yn .* (1 + k1 * r.^2 + k2 * r.^4 + k3 * r.^6);

xd_t = 2 * p1 * xn .* yn + p2 * (r.^2 + 2 * xn.^2);
yd_t = p1 * (r.^2 + 2 * yn.^2) + 2 * p2 * xn .* yn;

xd = xd_r + xd_t;
yd = yd_r + yd_t;

% Convert to 2D array
xd = reshape(ceil(xd),[ydim, xdim]);
yd = reshape(ceil(yd),[ydim, xdim]);

% Get projections
imageCylindrical = imwarp(image, cat(3,xd, yd), 'cubic');
imshow(imageCylindrical)

Please let me know what I am doing wrong using the imwarp function. Any help is appreciated!

Answer 1

[Avni Agrawal]

Hi @Preetham Manjunatha,

I understand that you are trying to plot cylindrical projection of checkers image using imwarp.

The issue you're encountering with imwarp not working as expected for cylindrical projection stems from a misunderstanding of how imwarp functions and the way you're attempting to use it. imwarp is designed to warp an image based on a geometric transformation object (such as affine2d, projective2d, etc.) or a spatial referencing object, rather than directly mapping pixel coordinates.

For cylindrical projection or any complex geometric transformation that cannot be represented by a simple affine or projective transformation, you typically need to perform the mapping manually or use a different approach. Your current method attempts to directly pass the warped pixel coordinates (xd, yd) as if they were a transformation object, which is not how imwarp operates.

Here's an approach to achieve cylindrical projection by manually interpolating the pixel values based on the computed (xd, yd) coordinates. This involves using interpolation functions such as interp2 to map the original image pixels to their new locations:

clc; close all; clear;

% Inputs
fileName = 'checker.jpg';

% Focal lengths
fx = 300;
fy = 300;

% Read image
image = (imread(fileName));

% Get image size
[ydim, xdim, bypixs] = size(image);

% Camera intrinsics
K = [fx, 0, xdim/2; 0, fy, ydim/2; 0, 0, 1];

% Distortion coefficients [k1, k2, k3, p1, p2]
DC = [0, 0, 0, 0, 0];

% Get distrotion coefficients
fx = K(1,1);
fy = K(2,2);
k1 = DC(1);
k2 = DC(2);
k3 = DC(3);
p1 = DC(4);
p2 = DC(5);

% Get image size
[ydim, xdim, bypixs] = size(image);

% Get the center of image
xc = xdim/2;
yc = ydim/2;

% Create X and Y coordinates grid
[X,Y] = meshgrid(1:xdim, 1:ydim);

% Perform the cylindrical projection
theta = (X - xc) / fx;
h = (Y - yc) / fy;

% Cylindrical coordinates to Cartesian
xcap = sin(theta);
ycap = h;
zcap = cos(theta);

xyz_cap = cat(3, xcap, ycap, zcap);
xyz_cap = reshape(xyz_cap,[],3);

% Normalized coords
xyz_cap_norm = (K * xyz_cap')';
xn = xyz_cap_norm(:,1) ./ xyz_cap_norm(:,3);
yn = xyz_cap_norm(:,2) ./ xyz_cap_norm(:,3);

% Radial and tangential distortion 
r = xn.^2 + yn.^2;
xd_r = xn .* (1 + k1 * r.^2 + k2 * r.^4 + k3 * r.^6);
yd_r = yn .* (1 + k1 * r.^2 + k2 * r.^4 + k3 * r.^6);

xd_t = 2 * p1 * xn .* yn + p2 * (r.^2 + 2 * xn.^2);
yd_t = p1 * (r.^2 + 2 * yn.^2) + 2 * p2 * xn .* yn;

xd = xd_r + xd_t;
yd = yd_r + yd_t;

% Convert to 2D array
xd = reshape(ceil(xd),[ydim, xdim]);
yd = reshape(ceil(yd),[ydim, xdim]);

% Get projections
% Assuming xd and yd are the new positions for each pixel
% Convert them to the coordinate system used by interp2
% Interpolate for each color channel
for k = 1:size(image, 3)
    imageCylindrical(:,:,k) = interp2(X, Y, double(image(:,:,k)), xd, yd, 'cubic', 0);
end

% Display the result
imageCylindrical = uint8(imageCylindrical);
imshow(imageCylindrical);

This method creates a more complex cylindrical projection that includes perspective distortion, simulating a cylinder that tilts away at the edges. Adjust the fx and fy values to control the extent of the distortion and the appearance of the cylinder.

I hope this helps!

[Preetham Manjunatha]

Thank you, @Avni Agrawal! Later, I also implemented interp2 solution. Out of curiosity, how is the image produced by you is so crisp along the edges? I have some aliasing issues and spur pixels at the edges using your code.