function [xs,ys] = polygx(xp,yp,xl,yl)
% XPOLYG   polyline polygon intersection
% [XS,YS] = POLYGX(XP,YP,XL,YL) returns as a NaN separated line the part of the
% NaN separated line XL,YL inside the polygone XP,YP.
%
% [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

%xp = [-1 1 1 0 -1 -1]'; yp = [-1 -1 1 .25 1 -1]';
%xl = [-2 0 -1 -2 2];  yl = [-1 0 0 1 0];
%xl = xl(:).'; yl = yl(:).';

% reshape arguments
xl = xl(:).'; yl = yl(:).'; nl = length(xl);
xp = xp(:)  ; yp = yp(:)  ;
% close polygon
if xp(1) ~= xp(end) || yp(1) ~= yp(end)
 xp = [ xp ; xp(1) ];
 yp = [ yp ; yp(1) ];
end
np = length(xp);

% find line points inside polygon by
% rotating and scaling polygon copies so that each line segment is (0,0)-(1,0)
[sl,xpl,ypl,k23,dxl,dyl] = inpolyg(xl,yl,xp,yp,true);

% line polygon intersection
% line segment may cross polygon side only for jumps Q0<>Q2, Q1<>Q3 and Q0<>Q3
% for those compute and select crossing distance on (0,0)-(1,0) segment
% forget nl-th segment (artificially added and being (-1,0)-(0,0)
l = k23+1;
fl = xpl(k23) - ypl(k23) .* (xpl(l) - xpl(k23)) ./ (ypl(l) - ypl(k23));
[l,kl] = ind2sub([np,nl],k23);
k = find(fl > 0 & fl < 1 & kl < nl);

% build a new line with original and intersection points
kl = [ 1:nl kl(k)' ]; fl = [ zeros(1,nl) fl(k)' ]; sl = [ sl ones(1,length(k)) ];
[l,k] = sort(kl+fl); kl = kl(k); fl = fl(k); sl = sl(k);

% intermediate points on segments defined by two edge points 
ks = find( sl(1:end-1) == 1 & sl(2:end) == 1 );
k = kl(ks); f = ( fl(ks+1) + fl(ks) ) / 2;
xs = xl(k) + dxl(k) .* f;
ys = yl(k) + dyl(k) .* f;
% find and keep those outside polygon
k = find(inpolyg(xs,ys,xp,yp) == 0);
% add them to the line
kl = [ kl kl(ks(k)) ] ; fl = [ fl f(k) ]; sl = [ sl zeros(1,length(k)) ];
[l,k] = sort(kl+fl); kl = kl(k); fl = fl(k); sl = sl(k);
% mark outside points as NaN
fl(sl == 0) = NaN;
xs = xl(kl) + dxl(kl).*fl;
ys = yl(kl) + dyl(kl).*fl;

% make it compact
k = isnan([xs NaN]);
k = find( k(1:end-1) & k(2:end) );
xs(k) = []; ys(k) = [];
if ~isempty(xs) && isnan(xs(1))
 xs(1) = []; ys(1) = [];
end
% if closed move 1st segment at the end
k = min(find(isnan(xs)));
if ~isempty(k) && xs(1) == xs(end) && ys(1) == ys(end)
 l = [k+1:length(xs) 1:k-1];
 xs = xs(l);
 ys = ys(l);
end

function [sl,xpl,ypl,k23,dxl,dyl] = inpolyg(xl,yl,xp,yp,rot)

nl = length(xl);
np = length(xp);

% center a copy of the polygon at each line point
xpl = repmat(xp,1,nl) - repmat(xl,np,1);
ypl = repmat(yp,1,nl) - repmat(yl,np,1);

% rotate and scale polygons so that each line segment is (0,0)-(1,0)
if nargin == 5 && rot
 dxl = [diff(xl) 0]; dxl(end) = dxl(end-1); 
 dyl = [diff(yl) 0]; dyl(end) = dyl(end-1);
 d2 = dxl.^2 + dyl.^2;
 dxd2 = repmat(dxl./d2,np,1);
 dyd2 = repmat(dyl./d2,np,1);
 [xpl,ypl] = deal( dxd2.*xpl + dyd2.*ypl , dxd2.*ypl - dyd2.*xpl );
end

% quandrant jumps; quandrants are labelled ccw Q0 to Q4
xn = xpl < 0;
yn = ypl < 0;
dq = [ diff( (xn & ~yn) + 2*(xn & yn) + 3*(~xn & yn) ) ; zeros(1,nl) ];
% for jumps Q0<>Q2 and Q1<>Q3, derive direction from cross product
k2 = find(abs(dq) == 2); l = k2+1;
dq(k2) = 2*sign(xpl(k2).*ypl(l) - xpl(l).*ypl(k2));
% jump Q0<>Q3 is a -1/+1 jump
k3 = find(abs(dq) == 3); dq(k3) = -sign(dq(k3));
k23 = [k2;k3];
% total jump to marker: +-4 > 2 > inside ; +-2 > 1 > edge ;  0 > 0 > outside
sl = round( abs(sum(dq)) / 2 );




