## Copyright (C) 2015,2016 Philip Nienhuis
## 
## 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
## (at your option) 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{XYo} =} optimize_branch_cuts (@var{XYi})
## @deftypefnx {Function File} [@var{Xo}, @var{Yo} ] = optimize_branch_cuts (@var{Xi}, @var{Yi})
## @deftypefnx {Function File} [@var{Xo}, @var{Yo}, @var{Zo} ] = optimize_branch_cuts (@var{Xi}, @var{Yi}, @var{Zi})
## Connect outer (filled) polygon and inner (hole) polygons using branch cuts
## such that all polygons are connected as one string of vertices, for
## subsequent plotting polygons with holes using patch().
##
## @var{XYi} can be a 2 or 3 dimensional array; only the X and Y coordinates
## will be optmized and Z-values will be kept with their original vertices.
## Alternatively separate X, Y and optionally Z-vectors can be specified.
## The outer polygon should appear as a first subvector, bounded by NaN values,
## and have a counterclockwise winding direction.  All subsequent inner hole
## polygons should also be bounded by NaN values and have clockwise winding
## directions.
##
## optimize_branch_cuts expects and returns column vectors or matrices where
## each row contains coordinates of a vertex.
##
## @seealso{polyjoin}
## @end deftypefn

## Author: Philip Nienhuis <prnienhuis@users.sf.net>
## Created: 2016-04-30

function [X, Y] = optimize_branch_cuts (XX, YY=[], ZZ=[]);

  ## Check for geometry package
  persistent ofgeom = [];
  if (isempty (ofgeom))
    if (isempty (which ("distancePoints")))
      ofgeom = 0;
    else
      ofgeom = 1;
    endif
  endif
  if (! ofgeom)
    error ("geometry package not loaded (or not even installed ?)");
  endif
  matinp = 0;

  XYsz = size (XX, 2);
  ## Check args
  if (nargin == 1)
    if (ismatrix (XX) && XYsz >= 2)
      ## apparently matrix input rather than two separate X,Y[,Z] vectors
      matinp = 1;
      XY = XX;
    endif
  elseif (nargin >= 2)
    ## Separate vector input
    XY = [XX YY ZZ];
  elseif (nargin < 1)
    error ("optimize_branch_cuts: not enough input arguments");
  endif

  ## Also keep track of Z.  Z isn't (yet) in the branch cut optimization (but
  ## that could be done easily)
  if (isempty (ZZ) || XYsz == 2)
    ## At least provide pointers where Z coordinates have gone in output arrays
    ZZ = [1:size(XY, 1)]';
    XY = [ XY ZZ ];
  endif

  ## Find NAN separators
  ipt = [0; find(isnan (XY(:, 1))); numel(XY(:, 1))+1];
  if (numel (ipt) < 3)
    ## No NaN separators => no subfeatures. Return input args
    if (! matinp)
      X = XX;
      Y = YY;
      if (nargin == 3)
        Z = ZZ;
      endif
    else
      X = XX;
      Y = Z = [];
    endif
    return
  endif
  for ii=1:numel (ipt) - 1
    ## Check for closed polygon
    if (any (abs (XY(ipt(ii)+1, 1:2) - XY(ipt(ii+1)-1, 1:2)) > eps))
      ## Duplicate first vertex as last vertex of subpolygon
      ## First shift all subpolys down
      XY(ipt(ii+1)+1:end+1, :) = XY(ipt(ii+1):end, :);
      ipt(ii+1:end) += 1;
      XY(ipt(ii+1)-1, :) = XY(ipt(ii)+1, :);
    endif
    XY(ipt(ii)+1:ipt(ii+1)-1, 4) = [ ipt(ii)+1:ipt(ii+1)-1 ]';
  endfor

  ## Silence broadcasting warning
  warning ("off",  "Octave:broadcast", "local");
  
  ## Compute all interdistances
  XY(ipt(2:end-1), :) = [];
  dists = distancePoints (XY(:, 1:2), XY(:, 1:2));
  szdst = size (dists);
  dists = dists + tril (Inf (szdst(1)));
  X_Y = XY(1:ipt(2)-1, :);

  ## Keep track of which holes are still unconnected
  processed = [0 (ones (1, numel (ipt) - 2))];
  tt = cumsum (diff (ipt) - 1);
  idx = [tt(1:end-1)+1 tt(2:end)];
  odx = 1 : ipt(2) - 1;
  ody = 1:size (dists, 2);

  ## Find hole polygon with smallest distance to an outer vertex; afterwards
  ## assign that to outer vertex + restart search until all holes are processed.
  ## FIXME Although Octave doesn't draw the branch cuts, it may be better to
  ##       also check that branch cuts do not cross polygons between two
  ##       vertices belonging to other polygons (or self-intersect polygons)
  while (any (processed))
    ## Get slice of dists with outer vertices + vertices connected to it, excl.
    ## columns already processed
    odz = setdiff (ody, odx);
    [~, indx] = min (dists(odx, odz)(:));
    ## Get subscripts into sliced dists matrix
    [r, c] = ind2sub ([numel(odx),numel(odz)], indx);
    ## Recompute subscripts into full dists matrix
    rr = odx(r);             ## Needed to insert new hole into place
    cc = odz(c);             ## To rotate hole such that this vertex has
                             ## smallest distance to outer polygon vertex
    ## Find hole polygon corresponding to these subscripts
    ii = find (cc >= idx(:, 1) & cc < idx(:, 2));
    shft = idx(ii, 1) - cc;
    tmpXY = XY(idx(ii, 1):idx(ii, 2), :);
    if (shft)
      tmpXY(1:end-1, :) = circshift (tmpXY(1:end-1, :), [shft, 0]);
      ## tmpXY is always shifted upward here. Copy toprow coordinates to bottom
      tmpXY(end, :) = tmpXY(1, :);
    endif
    X_Y = [X_Y(1:r, :); tmpXY; X_Y(r:end, :)];
    odx = [odx idx(ii, 1):idx(ii, 2)];
    processed(ii+1) = 0;
    
  endwhile

  if (! matinp)
    X = X_Y(:, 1);
    Y = X_Y(:, 2);
    if (nargin == 3)
      Z = X_Y(:, 3);
    endif
  else
    X = X_Y(:, 1:XYsz);
    Y = Z = [];
  endif

endfunction
