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