function [cyl_img, mask] = cylinder_projection(img, focal_length, k1, k2)

width = size(img,2);
height = size(img,1);
mask = zeros(height, width);
%img2 = zeros(height, width, 3);
locsx = zeros(height, width);
locsy = zeros(height, width);


for r=1:height
    for c=1:width
        %theta = x_cyl / xc * pi;
        
        %theta = %pi * (1 - c / width);%%(x_cyl - xc) / focal_length;
        theta = (c - width / 2) / focal_length;
        h = (height / 2 - r) / focal_length;
        
        xhat = sin(theta);
        yhat = h;
        zhat = cos(theta);
        
        xn = xhat / zhat;
        yn = yhat / zhat;
        
        r_sqr = xn^2 + yn^2;
        
        xd = xn * (1 + k1 * r_sqr + k2 * r_sqr^2);
        yd = yn * (1 + k1 * r_sqr + k2 * r_sqr^2);
        
        rsrc = floor(height / 2 - (focal_length * yd));
        csrc = floor(width / 2 + (focal_length * xd));
        
        if rsrc > 0 && rsrc <= height && csrc > 0 && csrc <= width
            img2(r, c, 1) = uint8(img(rsrc, csrc, 1));
            img2(r, c, 2) = uint8(img(rsrc, csrc, 2));
            img2(r, c, 3) = uint8(img(rsrc, csrc, 3));
            mask(r, c) = 1;
        else
            img2(r, c, 1) = uint8(0);
            img2(r, c, 2) = uint8(0);
            img2(r, c, 3) = uint8(0);
            mask(r, c) = 0;
        end
    end
end

cyl_img = img2;
        


