function [J,p] = ell(I, a, b, x0, y0, phi, jcols, jrows, s)
%ELL  Map an image from Cartesian to elliptical coordinates.
% [J,p] = ell(I, a, b, x0, y0, phi, jcols, jrows, s)
%
% Input Arguments:
%   I     - input image
%   a     - length of the major axis
%   b     - length of the minor axis (a > b)
%   x0    - abscissa of the center of the ellipse
%   y0    - ordinate of the center of the ellipse
%   phi   - angle between x-axis and the major axis
%   jcols - number of columns of the output image
%   jrows - number of rows of the output image
%   s     - scale factor of radial axis
%
% Outputs:
%   J     - output image
%   p     - transformation parameters (used for inverse transform)
%

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

if a < b | a < 0 | b < 0,
    error('Invalid arguments.');
end


I = im2double(I);
% signal channel only
if size(I,3) == 3, I = rgb2gray(I);
elseif size(I,3) > 1, I = I(:,:,1); 
end

[irows,icols] = size(I);
J = zeros(jrows, jcols);
[u,v] = meshgrid([0:1/(jcols-1):1] * s, 0:pi*2/(jrows-1):pi*2);

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

% elliptical to Cartesian coordinates mapping
x = c * cosh(u) .* cos(v);
y = c * sinh(u) .* sin(v);

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

% bilinear interpolation
J = interp2(I, X, Y, 'linear', 0);

% save parameters
p.a = a; p.b = b; p.x0 = x0; p.y0 = y0; p.phi = phi;
p.icols = icols; p.irows = irows;
p.s = s;