function data = spiralFromEllipse(grid, M, d, offsetR, offsetTh)
% spiralEllipse: Generate an implicit surface function for a wound spiral.
%
%   data = spiralFromEllipse(grid, M, d, offsetR, offsetTh)
%
% Generates an implicit surface function for a wound spiral.  This routine
%   is an attempt to duplicate the wound spiral interface found in Osher &
%   Sethian 88 and many subsequent level set publications discussing
%   curvature dependent flow; for example, O&F figure 4.1.  While it is
%   visually similar, simulations indicate that it probably isn't exactly
%   right.
%
% The spiral is generated by constructing an ellipse in polar coordinates.
%   To get a wound spiral (whose head is more than one revolution removed
%   from its tail), we create duplicated copies of the grid with polar
%   angular coordinates theta running beyond the range [ -pi, +pi ].  The
%   implicit surface function for the ellipsoid is calculated over the
%   entire range of theta, and then we take a minimum over equivalent values
%   of theta.
%
% The implicit surface function for the ellipsoid is
%
%                 f(x) = sqrt(y' * M * y) - d,
%
%   where y = cart2pol(x).
%
% 
% Parameters:
%   grid         Grid on which the spiral will be defined.
%   M            Shape matrix for the ellipsoid.
%   d            Offset for implicit surface function.
%   offsetR      Offset of ellipse's center in the r coordinate.
%   offsetTh     Offset of ellipse's center in the theta coordinate.
%
%   data         Implicit surface function for the wound spiral.
%
%
% Experimentation seems to indicate the following parameters generate a
%   visually reasonable approximation of the spiral from O&F figure 4.1:
%
%                   M = diag([ 0.06, 2 * pi ].^(-2));
%                   d = 1;
%                   offsetR = 0.5;
%                   offsetTh = -pi/3;
%
%                   angleRot = pi / 100;
%                   matrixRot = [ cos(angleRot), -sin(angleRot); ...
%                                 sin(angleRot), cos(angleRot) ];
%                   angleM = matrixRot * M * matrixRot';
%                   data = spiral(g, angleM, d, offsetR, offsetTh);
%
% Note that without the slight rotation by angleRot, a torus is generated
%   instead of a spiral.
  
% Copyright 2004 Ian M. Mitchell (mitchell@cs.ubc.ca).
% This software is used, copied and distributed under the licensing 
%   agreement contained in the file LICENSE in the top directory of 
%   the distribution.
%
% Ian Mitchell 2/19/04
  
  if(grid.dim ~= 2)
    error('Spiral can only be created on two dimensional grids');
  end
  
  % Get polar coordinates.
  [ th1, r1 ] = cart2pol(grid.xs{1}, grid.xs{2});
  
  % The polar coordinates have th \in [ -pi, +pi ].
  %   Make copies with th beyond this range.
  r = { r1; r1; r1; r1; r1 };
  th = { th1 - 4 * pi; th1 - 2 * pi; th1; th1 + 2 * pi; th1 + 4 * pi };
  copyN = length(r);
  
  cellM = num2cell(M);
  data = inf * ones(grid.shape);
  for i = 1 : copyN
    % Calculate the ellipsoidal implicit surface function for each copy.
    state = { r{i} - offsetR; th{i} - offsetTh };
    temp = cellMatrixMultiply(cellM, state);
    temp = cellMatrixMultiply(state', temp);
    implicit = sqrt(temp{1}) - d;
    
    % Union together the copies by taking the minimum.
    data = min(data, implicit);
    
  end
