MultiMethods

function xy = sleerfid(a, xy, wts)
    xy1 = try09(a, xy, wts);
    n1 = numKnots(xy1,a);
    xy2 = hannes(a, xy, wts);
    n2 = numKnots(xy2,a);
    xy3 = blackwidow(a, xy, wts);
    n3 = numKnots(xy3,a);
    if ((n1<n2)&(n1<n3))
        xy = xy1;
    elseif (n2<n3)
        xy = xy2;
    else
        xy = xy3;
    end;
end

function xyOut = try09(a, xyIn, wts)

N=size(xyIn,1);

D=min(N,mindist(a).^0.25);
mD=mean(D,1);
E=D-repmat(mD,[size(D,1),1])-repmat(mD',[1,size(D,2)])+mean(mD);
[x,nill]=svd(E);
xyOut0=x(:,1:2); 

beta1=10;
beta2=120;
lambda=.1;
thr=.01;
xyOut=N*xyOut0+randn(size(xyOut0))/1e3;
for ncut=1:2
    M=a; M(1:size(M,1)+1:end)=1;
    M=bsxfun(@rdivide,M,sum(M,2));
    Mbeta=(((1-lambda)*eye(N)+lambda*M)^beta1);
    for n1=1:beta2
        xyOut=Mbeta*xyOut;
        mxyOut=mean(xyOut,1);
        cxyOut=cov(xyOut,1);
        [c1,d1,c2]=svd(cxyOut);
        xyOut=(xyOut-mxyOut(ones(N,1),:))*c1*diag(1./sqrt(.1+diag(d1)))*c2'; %sqrtm(pinv(.1*eye(2)+cxyOut));
    end
    [i,j]=find(a>thr); 
    dd=sum((xyOut(i,:)-xyOut(j,:)).^2,2);
    k=find(dd>4*mean(dd));
    if isempty(k), break; end
    a(i(k)+N*(j(k)-1))=thr;
end
xyOut0=xyOut;


xyOut0=sqrt(N)*detrend(xyOut0,'constant')*diag(1./max(eps,std(xyOut0,1,1)));
K=10;
k=1;
[sxyOut0,idxequal]=sortrows(round(K*xyOut0)/K);
idxequal=idxequal(all(~diff(sxyOut0,1,1),2));
xyOut=round(xyOut0*k);
while size(unique(xyOut,'rows'),1)~=N
    k=k*2;
    xyOut=xyOut0*k;
    xyOut(idxequal,:)=(xyOut0(idxequal,:)+randn(numel(idxequal),2)/K)*k;
    xyOut=round(xyOut);
end

if N<30 % small map knot lines
 xyOut=3*xyOut+(randi(3,size(xyOut))-2);
end

dxy = round(wts*(xyOut-xyIn)./sum(wts)); % Howe's
xyOut = bsxfun(@minus,xyOut,dxy);

end

function D=mindist(C)

N=size(C,1);
X=logical(speye(N));
D=inf(N,N);

D(X)=0;
for n=1:N,
    X=(C*X)>0;
    X=X&(D>n);
    if ~any(X(:)),break;end
    D(X)=n;
end

end

function xyOut = hannes(a, xyIn, wts)
sz=size(a,1);
deg=diag(sum(a));
L=diag(sum(a))-a;
[V,D]=eig(L,deg);
xyOut=0;
outdiam=sqrt(sum((max(V(:,2:3))-min(V(:,2:3))).^2));
indiam=sqrt(sum((max(xyIn)-min(xyIn)).^2));
mult=indiam/outdiam;
while size(unique(xyOut,'rows'),1)<sz
xyOut=round(V(:,2:3)*mult);
[u,i]=unique(xyOut,'rows');
m=setdiff(1:sz,i);
xyOut(m,:)=xyOut(m,:)+round(rand(length(m),2))*2-1;
mult=mult*1.5;
end;
end

function xy = blackwidow(a, xy, wts)

nn = length(a);
mid = mean(xy);
dev = std(xy);
[u,s,v] = svd(xy);
if (s(2,2)<0.1)
    xy = xy+bsxfun(@times,2*sqrt(s(1))*randn(nn,1),v(:,2)');
end;

%logical a for outer boundary finding
la = logical(a);
deg = sum(la);
while any(deg==1)
    la(deg==1,:) = 0;
    la(:,deg==1) = 0;
    deg = sum(la);
end;

ia = logical(a*a*a*a+eye(nn));
if ~all(ia(:))
    xy = relax2(la, ia, xy, 10, false);
end;
ia = logical(a*a*a+eye(nn));
if ~all(ia(:))
    xy = relax2(la, ia, xy, 20, false);
end;
ia = logical(a*a+eye(nn));
if ~all(ia(:))
    xy = relax2(la, ia, xy, 20, false);
end;
ia = logical(a+eye(nn));
xy = relax2(la, ia, xy, 10, nn>50);
%wa = edgeWeights(a,xy);
xy = relax2(la, ia, xy, 30, nn>50);

%         gplot(ia,xy,'o-');
%         axis equal;
%         waitforbuttonpress;
% % synthetic spherical deprojection
%         d = sqrt(sum(xy.^2,2));
%         radius = max(d)*1.1;
%         th = acos(d./radius);
%         d2 = (pi/2-th)*radius;
%         xy = xy.*repmat(d2./d,1,2);
%         gplot(ia,xy,'o-');
%         axis equal;
%         waitforbuttonpress;
        
end

function xy = relax2(a, wa, xy, nrep, dohoop)
% Iterative method based on https://akpublic.research.att.com/areas/visualization/papers_videos/pdf/DBLP-journals-camwa-Koren05.pdf

% gplot(a,xy,'o-');
% axis equal;
% waitforbuttonpress;

nn = length(wa);
dev = std(xy);

% xymin = min(xy);
% xymax = max(xy);


% gplot(a,xy,'o-');
% axis equal;
% waitforbuttonpress;

for i = 1:nrep*nn
    j = randi(nn);
    nz = find(wa(j,:));
    w = wa(nz,j);
    wsum = sum(w);
    xy(j,1) = sum(xy(nz,1).*w)./wsum;
    xy(j,2) = sum(xy(nz,2).*w)./wsum;
    
    if (mod(i,nn)==0)
        xy = xy-repmat(mean(xy),nn,1);
        xy = xy.*repmat(dev./std(xy),nn,1);
        
        if dohoop
%             gplot(ia,xy,'o-');
%             axis equal;
%             waitforbuttonpress;
            xy = hoop(xy,a);
%             gplot(ia,xy,'o-');
%             axis equal;
%             waitforbuttonpress;
        end;
        
        [u,s,v] = svd(xy);
        xy  = u(:,1:2)*diag(sqrt([prod(diag(s)) prod(diag(s))]))*v';

%         gplot(ia,xy,'o-');
%         axis equal;
%         waitforbuttonpress;
    end;
end;

% gplot(a,xy,'o-');
% axis equal;
% waitforbuttonpress;

% force to grid
for i = 1:nn
    xy(i,:) = vicinityPoint(xy(i,:),xy(1:i-1,:));
end;
end

function xy = relax(a, ia, xy, nrep, dohoop)
% Iterative method based on https://akpublic.research.att.com/areas/visualization/papers_videos/pdf/DBLP-journals-camwa-Koren05.pdf

nn = length(ia);
dev = std(xy);

for i = 1:nrep*nn
    j = randi(nn);
    xy(j,:) = mean(xy(ia(:,j),:));
    
    if (mod(i,nn)==0)
        xy = xy-repmat(mean(xy),nn,1);
        xy = xy.*repmat(dev./std(xy),nn,1);
        
        if dohoop
            xy = hoop(xy,a);
        end;
        
        [u,s,v] = svd(xy);
        xy  = u(:,1:2)*diag(sqrt([prod(diag(s)) prod(diag(s))]))*v';
    end;
end;

% force to grid
for i = 1:nn
    xy(i,:) = vicinityPoint(xy(i,:),xy(1:i-1,:));
end;
end

function xy = hoop(xy,a)
k = outer(xy,a);
d = sqrt(sum(xy.^2,2));
maxd = max(d);
xy(k,:) = xy(k,:).*repmat(maxd./d(k),1,2);
end

function k = outer(xy,a)
%figure(1);
% clf
% gplot(a,xy,'o-');
% axis equal;
% hold on

nn = size(xy,1);
deg = sum(a);
ink = false(1,nn);
vc = (deg'>0)&(xy(:,1)==min(xy(deg>0,1)));
v = find(vc&(xy(:,2)==min(xy(vc,2))),1);
th = pi/2;
hp = atan2(xy(:,2),xy(:,1));
hp1 = hp(v);
chp = 0;
% plot(xy(v,1),xy(v,2),'r*');
i = 0;
while (chp < 2*pi)&(i<nn)
    i = i+1;
    ink(v) = true;
    hp0 = hp1;
    nbr = find(a(v,:));
    dir0 = atan2(xy(nbr,2)-xy(v,2),xy(nbr,1)-xy(v,1));
    dir = mod(dir0-hp(v),2*pi);
    dir(ink(nbr)) = dir(ink(nbr))-10;
    [~,iv] = max(dir);
    %figure(2);
    %clf; hold on; for i = 1:numel(nbr), plot(xy([v nbr(i)],1),xy([v nbr(i)],2),'o-'); end; axis equal; hold off
    v0 = v;
    v = nbr(iv);
    %figure(1);
%     plot(xy(v,1),xy(v,2),'r*');
%     title(v);
%     xlabel(chp);
    hp1 = hp(v);
    chp = chp+mod(hp0-hp1+pi,2*pi)-pi;
    th = atan2(xy(v0,2)-xy(v,2),xy(v0,1)-xy(v,1));
end;
k = find(ink);
end


function xy = vicinityPoint(xy0,pxy)
r = 0;
xy = [];
rxy0 = round(xy0);
if isempty(pxy)
    xy = rxy0;
end;
while isempty(xy)
    hits = find((pxy(:,1)>=rxy0(1)-r)&(pxy(:,1)<=rxy0(1)+r)&(pxy(:,2)>=rxy0(2)-r)&(pxy(:,2)<=rxy0(2)+r));
    filled = sparse(pxy(hits,2)-rxy0(2)+1+r,pxy(hits,1)-rxy0(1)+1+r,true,2*r+1,2*r+1);
    if ~all(filled(:))
        [xg,yg] = meshgrid(rxy0(1)+(-r:r)-xy0(1),rxy0(2)+(-r:r)-xy0(2));
        d = xg.^2+yg.^2;
        d(filled) = inf;
        [bdy,bdx] = find(d==min(d(:)),1);
        xy = rxy0+[bdx bdy]-r-1;
    else
        r = r+1;
    end
end;
end
function N = numKnots(XY,A)

% This function finds the number of knots in a network of points XY whose
% connectivity is given by A.

% Consider only unique lines (A is symmetric with ones on diagonal)
A = triu(A,1);
nLines = nnz(A);

% Find xy coordinates for first and second point on each line
[p1i,p2i] = find(A);
x1 = repmat(XY(p1i,1),1,nLines);
y1 = repmat(XY(p1i,2),1,nLines);
x2 = repmat(XY(p2i,1),1,nLines);
y2 = repmat(XY(p2i,2),1,nLines);
x3 = repmat(XY(p1i,1)',nLines,1);
y3 = repmat(XY(p1i,2)',nLines,1);
x4 = repmat(XY(p2i,1)',nLines,1);
y4 = repmat(XY(p2i,2)',nLines,1);

pick = tril(true(nLines),-1);
x1 = x1(pick);
y1 = y1(pick);
x2 = x2(pick);
y2 = y2(pick);
x3 = x3(pick);
y3 = y3(pick);
x4 = x4(pick);
y4 = y4(pick);

N = sum(areIntersecting(x1,y1,x2,y2,x3,y3,x4,y4));
end

function bool = areIntersecting(x1,y1,x2,y2,x3,y3,x4,y4)
% Determine if two line segments intersect.
%
% line1 = [x1, y1, x2, y2]
% line2 = [x3, y3, x4, y4]
%
% This function is robust to floating point precision issues assuming all
% input coordinates are integer-valued. If all coordinates have absolute
% values that do not exceed x, then the arithmetic operations used in this
% function could, as an intermediate step, produce a value with magnitude
% as large as 8*x^4. For double-precision floating point, this implies
% that the following must be satisfied:
%
% 8*x^4 < 2^53
%
% or equivalently, x < 5792.
%
% http://en.wikipedia.org/wiki/Line-line_intersection

bool = false(1,numel(x1));
gen = true(numel(x1),1);

% Test to see if the lines share exactly one endpoint. If so, the lines
% intersect if either of the free points lies on the line segment formed by
% the other two points.

s13 = find((x1==x3)&(y1==y3));
%s13 = find((x1==x3)&(y1==y3)&((x2~=x4)|(y2~=y4)));
bool(s13) = (isPointOnSegment(x2(s13), y2(s13), x1(s13), y1(s13), x4(s13), y4(s13))|isPointOnSegment(x4(s13), y4(s13), x1(s13), y1(s13), x2(s13), y2(s13)));
gen(s13) = false;
s14 = find((x1==x4)&(y1==y4));
%s14 = find((x1==x4)&(y1==y4)&((x2~=x3)|(y2~=y3)));
bool(s14) = (isPointOnSegment(x2(s14), y2(s14), x1(s14), y1(s14), x3(s14), y3(s14))|isPointOnSegment(x3(s14), y3(s14), x1(s14), y1(s14), x2(s14), y2(s14)));
gen(s14) = false;
s23 = find((x2==x3)&(y2==y3));
%s23 = find((x2==x3)&(y2==y3)&((x1~=x4)|(y1~=y4)));
bool(s23) = (isPointOnSegment(x1(s23), y1(s23), x2(s23), y2(s23), x4(s23), y4(s23))|isPointOnSegment(x4(s23), y4(s23), x1(s23), y1(s23), x2(s23), y2(s23)));
gen(s23) = false;
s24 = find((x2==x4)&(y2==y4));
%s24 = find((x2==x4)&(y2==y4)&((x1~=x3)|(y1~=y3)));
bool(s24) = (isPointOnSegment(x1(s24), y1(s24), x2(s24), y2(s24), x3(s24), y3(s24))|isPointOnSegment(x3(s24), y3(s24), x1(s24), y1(s24), x2(s24), y2(s24)));
gen(s24) = false;

% Next we check for parallel and coincident lines. Parallel lines
% obviously don't intersect. Coincident lines intersect if an endpoint from
% one of the lines lies on the other segment.

ss = find(gen&haveSameSlope(x1, y1, x2, y2, x3, y3, x4, y4));
gen(ss) = false;
si = ss(haveSameIntercept(x1(ss), y1(ss), x2(ss), y2(ss), x3(ss), y3(ss), x4(ss), y4(ss)));
bool(si) = ((isPointBetween(x3(si), x4(si), x1(si), 1) & isPointBetween(y3(si), y4(si), y1(si), 1)) | ...
               (isPointBetween(x3(si), x4(si), x2(si), 1) & isPointBetween(y3(si), y4(si), y2(si), 1)) | ...
               (isPointBetween(x1(si), x2(si), x3(si), 1) & isPointBetween(y1(si), y2(si), y3(si), 1)) | ...
               (isPointBetween(x1(si), x2(si), x4(si), 1) & isPointBetween(y1(si), y2(si), y4(si), 1)));

% If we get this far, we're in the general case of two well-defined
% non-parallel lines. The lines formed by the two segments must intersect
% somewhere. Find this intersection point and determine if it lies on the
% segments.

gen = find(gen);

Px_n = (x1(gen).*y2(gen) - y1(gen).*x2(gen)).*(x3(gen) - x4(gen)) - (x1(gen) - x2(gen)).*(x3(gen).*y4(gen) - y3(gen).*x4(gen));
Px_d = (x1(gen) - x2(gen)).*(y3(gen) - y4(gen)) - (y1(gen) - y2(gen)).*(x3(gen) - x4(gen));

Py_n = (x1(gen).*y2(gen) - y1(gen).*x2(gen)).*(y3(gen) - y4(gen)) - (y1(gen) - y2(gen)).*(x3(gen).*y4(gen) - y3(gen).*x4(gen));
Py_d = (x1(gen) - x2(gen)).*(y3(gen) - y4(gen)) - (y1(gen) - y2(gen)).*(x3(gen) - x4(gen));

bool(gen) = (isPointBetween(x1(gen), x2(gen), Px_n, Px_d) & ...
    isPointBetween(y1(gen), y2(gen), Py_n, Py_d) & ...
    isPointBetween(x3(gen), x4(gen), Px_n, Px_d) & ...
    isPointBetween(y3(gen), y4(gen), Py_n, Py_d));

end

function bool = haveSameIntercept(x1, y1, x2, y2, x3, y3, x4, y4)
bool = (x4 - x3).*(y1.*(x2 - x1) - x1.*(y2 - y1)) == ...
       (x2 - x1).*(y3.*(x4 - x3) - x3.*(y4 - y3));
end

function bool = haveSameSlope(x1, y1, x2, y2, x3, y3, x4, y4)
bool = (y2 - y1).*(x4 - x3) == (x2 - x1).*(y4 - y3);
end

function bool = isPointOnSegment(x1, y1, x2, y2, x3, y3)
% Determine if point (x1, y1) is on the segment formed by (x2, y2) and (x3, y3).
bool = haveSameSlope(x2, y2, x1, y1, x1, y1, x3, y3) & ...
       isPointBetween(x2, x3, x1, 1) & ...
       isPointBetween(y2, y3, y1, 1);
end

function bool = isPointBetween(pt1, pt2, num, den)
% Determine if either of the following is satisfied:
% pt1 <= num/den <= pt2
% pt1 >= num/den >= pt2

bool = ((pt1.*den <= num) & (num <= pt2.*den)) | ...
       ((pt1.*den >= num) & (num >= pt2.*den));
end