## 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');