Cylindrical Image Projection using Matlab's imtransform Function
In this short tutorial we will construct a custom transformation for projecting an image - given it's camera rotation R and calibration K - onto a cylindrical surface.
Contents
Transformation definition
For constructing a tform object we need to define the forward and backward transformation function that take Nx2 input coordinate matrices and return Nx2 output coordinates. First we choose the cylinder to have radius 1 and its principal axis to be parallel to the Z-direction, i.e. x^2 + y^2 = 1. The transformation between an image point x and the corresponding ray direction r then is as follows.
dir2cyl = @(r) [-atan2(r(:, 2), r(:, 1)), ... r(:, 3) ./ sqrt(sum(r(:, 1 : 2).^2))]; cyl2dir = @(xy) [cos(-xy(:, 1)), ... sin(-xy(:, 1)), ... xy(:, 2)];
The final transformation requires the convertion of Euclidean vectors into homogeneous entities and vice versa.
hom2euk = @(x) x(:, 1 : 2) ./ x(:, [3 3]); euk2hom = @(x) [x, ones(size(x, 1), 1)];
Combining these auxiliary functions with the general projection equation x = K * R' * [eye(3), X0] * X we obtain the following functions.
forward = @(U, T) dir2cyl((T.tdata.R * T.tdata.K^-1 * euk2hom(U)')'); inverse = @(X, T) hom2euk((T.tdata.K * T.tdata.R' * cyl2dir(X)')');
Note that we do not specify the camera pose yet, but formulate the transformation depenting on the transformation structure T that is to be specified later.
Demo image
We will now load an example image, i.e. Matlab's "camera man". In the following we will synthetically define a corresponding camera pose and compute the projection onto the unit cylinder.
I = imread('cameraman.tif');
imshow(I);
Camera definition
Before constructing a tform for image transformations, we define the camera pose with rotation matrix R, approximately pointing into X-direction with random pertubations, and a calibration matrix K with focal length c = 100 and principal point at the image center. Note that the calibration matrix has a negative focal length in X-direction in order to take account of the left handed image coordinate system.
rot = @(p, a) circshift(blkdiag([cos(p) -sin(p); sin(p) cos(p)], 1), [a a]);
R = rot(pi / 2, 1) * rot(0.5 * randn, 3) * rot(0.5 * randn, 2) * rot(0.5 * randn, 1);
c = 100;
K = [-c, 0, size(I, 2) / 2;
0, c, size(I, 1) / 2;
0, 0, 1];
Image transformation
Afterwards we can call maketform with both out forward and backward function as well as additional transformation data, i.e. the camera orientation and calibration.
tform = maketform('custom', 2, 2, forward, inverse, struct('K', K, 'R', R));
Using the tform object we can easily transform our image I with Matlab's imtransform function. Since the vertical output coordinates are possibly unbounded, we should specify the output coordinate range. We choose a full revolution in horizontal direction and -1 to 1 for the vertical direction, being equivalent to 90 degrees vertical field of view. Finally, a scale of 1/100 yields an output image of size 628x200 rather than 6x2 pixels.
xdata = [-pi, pi]; ydata = [-1, 1]; scale = 1 / 100; J = imtransform(I, tform, 'XData', xdata, 'YData', ydata, 'XYScale', scale);
To visualize the image with correctly labeled axes, we again need to specify the coordinate range.
imshow(J, 'XData', xdata, 'YData', ydata); axis on; set(gca, 'YDir', 'normal');
Removing the "mirror" image
As we might notice from the result, the image is projected onto the cylinder twice: once in front of the camera and once behind. In order to avoid this mirroring effect - or at least to detect those incorrectly projected pixels - we check the direction of each pixel of the output image.
[y, x] = ndgrid( ... linspace(ydata(1), ydata(2), size(J, 1)), ... linspace(xdata(1), xdata(2), size(J, 2))); r = cyl2dir([x(:), y(:)]);
Pixels are invalid if their direction vector differs from the negative camera-Z-axis by more than 90 degrees - or yields a negative scalar product.
invalid = r * -R(:, 3) < 0;
For this demonstration we dim those pixels to 20 percent brightness.
J(invalid) = 0.2 * J(invalid); imshow(J, 'XData', xdata, 'YData', ydata); axis on; set(gca, 'YDir', 'normal');
