%% Copyright (c) 2012 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
%%
%%    This program is free software: you can redistribute it and/or modify
%%    it under the terms of the GNU General Public License as published by
%%    the Free Software Foundation, either version 3 of the License, or
%%    any later version.
%%
%%    This program is distributed in the hope that it will be useful,
%%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%%    GNU General Public License for more details.
%%
%%    You should have received a copy of the GNU General Public License
%%    along with this program. If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @deftypefn {Function File} {@var{polyline} = } curve2polyline (@var{curve})
%% @deftypefnx {Function File} {@var{polyline} = } curve2polyline (@dots{},@var{property},@var{value},@dots{})
%% Adaptive sampling of a parametric curve.
%%
%% The @var{curve} is described as a 2-by-N matrix. Rows correspond to the
%% polynomial (compatible with @code{polyval}) describing the respective component
%% of the curve. The curve must be parametrized in the interval [0,1].
%% The vertices of the polyline are accumulated in regions of the curve where
%% the curvature is higher.
%%
%% @strong{Parameters}
%% @table @samp
%% @item 'Nmax'
%% Maximum number of vertices. Not used.
%% @item 'Tol'
%% Tolerance for the error criteria. Default value @code{1e-4}.
%% @item 'MaxIter'
%% Maximum number of iterations. Default value @code{10}.
%% @item 'Method'
%% Not implemented.
%% @end table
%%
%% @seealso{shape2polygon, curveval}
%% @end deftypefn

%% This function is based on the algorithm described in
%% L. H. de Figueiredo (1993). "Adaptive Sampling of Parametric Curves". Graphic Gems III.
%% I had to remove the recursion so this version could be improved.
%% Thursday, April 12 2012 -- JuanPi

function [polyline t bump]= curve2polyline (curve, varargin)
%% TODO make tolerance relative to the "diameter" of the curve.

  # --- Parse arguments --- #
  parser = inputParser ();
  parser.FunctionName = "curve2polyline";
  parser = addParamValue (parser,'Nmax', 32, @(x)x>0);
  parser = addParamValue (parser,'Tol', 1e-4, @(x)x>0);
  parser = addParamValue (parser,'MaxIter', 10, @(x)x>0);
  parser = parse(parser,varargin{:});

  Nmax      = parser.Results.Nmax;
  tol       = parser.Results.Tol;
  MaxIter   = parser.Results.MaxIter;

  clear parser toldef
  # ------ #

  t = [0; 1];
  tf = 1;
  points = 1;
  for iter = 1:MaxIter
    % Add parameter values where error is still bigger than tol.
    t = interleave(t, tf);
    nt = length(t);

    % Update error
    polyline = curveval (curve,t);
    bump = bumpyness(polyline);

    % Check which intervals must be subdivided
    idx = find(bump > tol);
    % The position of the bumps mpas into intervals
    % 1 -> 1 2
    % 2 -> 3 4
    % 3 -> 5 6
    % and so on
    idx = [2*(idx-1)+1; 2*idx](:);
    tf = false (nt-1,1);
    tf(idx) = true;

    if all (!tf)
      break;
    end

  end

endfunction

function f = bumpyness (p)
%% Check for co-linearity
%% TODO implement various method for this
%% -- Area of the triangle close to zero (used currently).
%% -- Angle close to pi.
%% -- abs(p0-pt) + abs(pt-p1) - abs(p0-p1) almost zero.
%% -- Curve's tange at 0,t,1 are almost parallel.
%% -- pt is in chord p0 -> p1.
%% Do this in isParallel.m and remove this function

  PL = p(1:2:end-2,:);
  PC = p(2:2:end-1,:);
  PR = p(3:2:end,:);

  a = PL - PC;
  b = PR - PC;

  f = (a(:,1).*b(:,2) - a(:,2).*b(:,1)).^2;

endfunction

function tt = interleave (t,varargin)

 nt   = length(t);
 ntt = 2 * nt -1;
 tt = zeros(ntt,1);
 tt(1:2:ntt) = t;
 beta = 0.4 + 0.2*rand(nt-1, 1);
 tt(2:2:ntt) = t(1:end-1) + beta.*(t(2:end)-t(1:end-1));

 if nargin > 1
  tf = true (ntt,1);
  tf(2:2:ntt) = varargin{1};
  tt(!tf) = [];
 end

endfunction

%!demo
%! curve    = [0 0 1 0;1 -0.3-1 0.3 0];
%! polyline = curve2polyline(curve,'tol',1e-8);
%!
%! t  = linspace(0,1,100)';
%! pc = curveval(curve,t);
%!
%! plot(polyline(:,1),polyline(:,2),'-o',pc(:,1),pc(:,2),'-r')
