function I = iell(J, p)
%IELL  Map an image from elliptical to Cartesian coordinates.
% I = iell(J, p)

%
% Copyright (c) 2007-2008 Kang Li (kangl@cmu.edu).
% All rights reserved.
%

[jrows,jcols] = size(J);
I = zeros(p.irows, p.icols);
[X,Y] = meshgrid(1:p.icols, 1:p.irows);

% apply rotation and translation
X = X - p.x0;
Y = Y - p.y0;
x = cos(p.phi) * X + sin(p.phi) * Y;
y = cos(p.phi) * Y - sin(p.phi) * X;

x(abs(x) < 1) = 1; % FIXME: eliminate numerical errors
y(abs(y) < 1) = 1; %

% compute focus
c = sqrt(p.a*p.a - p.b*p.b);

x = x / c;
y = y / c;

% map coordinates from elliptical to Cartesian
u1 = real(acosh(-sqrt(2)*x./sqrt(1+x.^2+y.^2-sqrt(4*y.^2+(-1+x.^2+y.^2).^2))));
v1 = real(acos (-sqrt(1+x.^2+y.^2 - sqrt(4*y.^2+(-1+x.^2+y.^2).^2)) /sqrt(2)));
u2 = real(acosh( sqrt(2)*x./sqrt(1+x.^2+y.^2-sqrt(4*y.^2+(-1+x.^2+y.^2).^2))));
v2 = real(acos ( sqrt(1+x.^2+y.^2 - sqrt(4*y.^2+(-1+x.^2+y.^2).^2)) /sqrt(2)));

u = zeros(size(X));
v = zeros(size(Y));

% process four quadrants separately
u(x< 0&y>=0) = u1(x< 0&y>=0);
v(x< 0&y>=0) = v1(x< 0&y>=0);
u(x< 0&y< 0) = u2(x< 0&y< 0);
v(x< 0&y< 0) = v2(x< 0&y< 0) + pi;

u(x>=0&y>=0) = u2(x>=0&y>=0);
v(x>=0&y>=0) = v2(x>=0&y>=0);
u(x>=0&y< 0) = u1(x>=0&y< 0);
v(x>=0&y< 0) = v1(x>=0&y< 0) + pi;

u = u / p.s * (jcols - 1);
v = v / (pi * 2) * (jrows - 1);
I = interp2(J, u + 1, v + 1, 'linear', 0);