function [new_image] = moebiusimage(image, inputs, outputs, type, size_final)
% INPUTS
% image(uint8 matrix) = image input;
% inputs(double matrix) = a 3x2 matrix with the initial positions x,y of 3
% points in the region [-pi,pi]x[-pi/2,pi/2], in the form [x_1 y_1; x_2 y_2;
% x_3 y_3];
% outputs(double matrix) = a 3x2 matrix with the final positions x,y of 3
% points without region restriction in R, in the form [nx_1 ny_1; nx_2 ny_2;
% nx_3 ny_3];
% type(string) = type of final projection.
% the types are: equirectangular, perspective and stereographic (complex
% plane);
% size_final(integer vector) = the final wanted resolution MxN pixels in the
% form [M N].
% -----
% OUTPUT
% newimage  = an uint8 matrix that can be read as an image with the
% respective projection trasformation.
% -----
% SYNTAX
% B = imagemoebius(A, inputs, outputs)
% Will do a transformation that sends the pixel points to new positions based 
% on the moebius transformation generated by the cross-ratio of
% <inputs>,<outputs> points and projects the final image on equirectangular
% domain;
% B = imagemoebius(A, inputs, outputs, type)
% Will do a transformation that sends the pixel points to new positions based 
% on the moebius transformation generated by the cross-ratio of
% <inputs>,<outputs> points and projects the final image using a <type>
% projection;
% B = imagemoebius(A, inputs, outputs, type, size_final)
% Will do a transformation that sends the pixel points to new positions based 
% on the moebius transformation generated by the cross-ratio of
% <inputs>,<outputs> points and projects the final image using a <type>
% projection with <size_final> resolution. 
% B = imagemoebius(A, inputs, outputs, type, size_final)
% Will do a transformation that sends the pixel points to new positions based 
% on the moebius transformation generated by the cross-ratio of
% <inputs>,<outputs> points and projects the final image using a projection 
% on equirectangular domain with <size_final> resolution.
% ------------------------------------------------------------------------
% version 2.4
% Leonardo Souto Ferreira, Leonardo Koller Sacht
% Federal University of Santa Catarina
% Mathematics Department
% 2016
% ------------------------------------------------------------------------
%uncatenate the colors of the image
R = double(image(:,:,1));
G = double(image(:,:,2));
B = double(image(:,:,3));
%obtain the size of the image
[nr nc]=size(R);
%Check the parameters
if nargin < 3
    error('Not suficient input values')
end
if nargin > 5
    error('Too many arguments')
end
if nargin == 3
    type = 'equirectangular';
    size_final = [nr nc];
end
if nargin == 4
    size_final = [nr nc];
end
if type == ':';
    type = 'equirectangular';
end
%Define the domain
aa = linspace(-pi,pi,nc);
cc = linspace(pi/2,-pi/2,nr); % had to flip this interval: it seems that meshgrid defines the y axis backwards
[AA,CC]=meshgrid(aa,cc);
%Define the counter-domain
% x = linspace(0,2*pi,size_final(1,2)); % this seems to be wrong
x = linspace(-pi,pi,size_final(1,2));  % this seems to be correct
y = linspace(pi/2,-pi/2,size_final(1,1)); % had to flip this interval: it seems that meshgrid defines the y axis backwards
[X,Y]=meshgrid(x,y);

%concatenate the inputs in a matrix
input_coord = cat(2, inputs, outputs);
if size(input_coord) ~= [3 4]
    error('Wrong entries for input and output points.')
end
N = 3;
%Transform the input points in sphere points
sphere_coord = zeros(3,6);
n = 1;
while n <= N
    sphere_coord(n,1) = cos(input_coord(n,2))*cos(input_coord(n,1));
    sphere_coord(n,2) = cos(input_coord(n,2))*sin(input_coord(n,1));
    sphere_coord(n,3) = sin(input_coord(n,2));
    sphere_coord(n,4) = cos(input_coord(n,4))*cos(input_coord(n,3));
    sphere_coord(n,5) = cos(input_coord(n,4))*sin(input_coord(n,3));
    sphere_coord(n,6) = sin(input_coord(n,4));
    n = n + 1;
end
%Projects the spheric input points onto stereografic plane
stereo_coord = zeros(3,4);
n = 1;
while n <= N
    stereo_coord(n,1) = 2*sphere_coord(n,2)/(sphere_coord(n,1) +1);
    stereo_coord(n,2) = 2*sphere_coord(n,3)/(sphere_coord(n,1) +1);
    stereo_coord(n,3) = 2*sphere_coord(n,5)/(sphere_coord(n,4) +1);
    stereo_coord(n,4) = 2*sphere_coord(n,6)/(sphere_coord(n,4) +1);
    n = n + 1;
end
%transform the stereographic input points onto complex numbers
complex_coord = zeros(3, 2);
n = 1;
while n <= N
    complex_coord(n,1) = stereo_coord(n,1) + stereo_coord(n,2)*1i;
    complex_coord(n,2) = stereo_coord(n,3) + stereo_coord(n,4)*1i;
    n = n + 1;
end
%Calculates the parameters of moebius transformation using a external
%function called crossmoebius
[a b c d] = crossmoebius(complex_coord(1,1), complex_coord(2,1), complex_coord(3,1),...
    complex_coord(1,2), complex_coord(2,2), complex_coord(3,2));
if abs(c) < 1e-10
    c = 1e-10 + 1e-10*1i;
end
%from final equirectangular projection to the sphere
x_sphere = cos(X).*cos(Y);
y_sphere = sin(X).*cos(Y);
z_sphere = sin(Y);
%From sphere to stereographic plane
X_C =(2*y_sphere)./(x_sphere + 1);
Y_C =(2*z_sphere)./(x_sphere + 1);
%From stereographic plane to complex numbers
Complex = X_C + Y_C*1i;
%inverse moebius transformation
Complex = Complex - a/c;
Complex = Complex*(-(c^2/(a*d - b*c)));
Complex = 1./Complex;
Complex = Complex - d/c;
%Returning to stereographic plane
X_C_n = real(Complex);
Y_C_n = imag(Complex);
%Fixing the mirrowing
r_polar = sqrt(X_C_n.^2 + Y_C_n.^2);
theta_polar = atan2(X_C_n,Y_C_n);
theta_polar = theta_polar + pi/2;
X_C_n = r_polar.*cos(theta_polar);
Y_C_n = r_polar.*sin(theta_polar);
%Return to the sphere
% x_sphere2 = (X_C_n.^2 + Y_C_n.^2 - 4)./(X_C_n.^2 + Y_C_n.^2 + 4);
x_sphere2 = -(X_C_n.^2 + Y_C_n.^2 - 4)./(X_C_n.^2 + Y_C_n.^2 + 4); % it needs a minus sign here
% y_sphere2 = (4*X_C_n)./(X_C_n.^2 + Y_C_n.^2 + 4);
y_sphere2 = -(4*X_C_n)./(X_C_n.^2 + Y_C_n.^2 + 4); % it needs a minus sign here
z_sphere2 = (4*Y_C_n)./(X_C_n.^2 + Y_C_n.^2 + 4);

%Return to the original equirectangular domain
U = atan2(y_sphere2, x_sphere2);
V = asin(z_sphere2);

switch type
    case 'equirectangular'
        U_final = U;
        V_final = V;
    case 'perspective'
        U_final = atan(U);
        V_final = atan(V.*cos(U_final));
    case 'complex'
        U_final = atan2((4*U),(4-U.^2-V.^2));
        V_final = asin(4*V./(U.^2+V.^2+4));
end

%Interpolate the points
r=(qinterp2(AA,CC,R,U_final,V_final,0));
g=(qinterp2(AA,CC,G,U_final,V_final,0));
b=(qinterp2(AA,CC,B,U_final,V_final,0));
%final image (concatenating matrix layers)
new_image = (uint8(cat(3,r,g,b)));

end