## 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 .
## -*- 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
## 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