% Author:  Nick Howe

function W = pickypicky2(B)

[nrow,ncol] = size(B);
xB(1:2:2*nrow,1:2:2*ncol) = B;
pin = (xB>0);
[Bi,Bj,Bv] = find(xB);
[xnrow,xncol] = size(xB);
s2i = reshape(1:xnrow*xncol,xnrow,xncol);
xC = zeros(xnrow,xncol);
xC(Bi+xnrow*Bj-xnrow) = (1:numel(Bv));
JxC1 = zeros(xnrow,xncol);
JxC2 = zeros(xnrow,xncol);
di = [-1 0 1 0];
dj = [0 -1 0 1];
change = 1;
[V,Vm,Vn] = unique(Bv);
if (numel(V)==1)
    ccn = numel(Bv);
else
    ccn = hist(Bv,V);
end;
rcc(V) = 1:numel(V);
%out = (ccn(Vn)==1);
%Bv(out) = [];
%Bi(out) = [];
%Bj(out) = [];
%Vn(out) = [];
%V = V(ccn>1);
%ccn = ccn(ccn>1);

% evaluate groups
up = zeros(1,numel(ccn));
down = zeros(1,numel(ccn));
for g = 1:numel(ccn)
    ipos = find(Bv==V(g));
    up(g) = ccn(g)*V(g);
    down(g) = (max(Bi(ipos))-min(Bi(ipos)))+(max(Bj(ipos))-min(Bj(ipos)));
end;

% join groups
P1 = zeros(xnrow,xncol);
P2 = zeros(xnrow,xncol);
P1src = zeros(xnrow,xncol);
J1 = zeros(xnrow,xncol);
JJ1 = zeros(xnrow,xncol);
J2 = zeros(xnrow,xncol);
JJ2 = zeros(xnrow,xncol);
JS = zeros(xnrow,xncol);
JJsrc = zeros(xnrow,xncol);
BB1 = zeros(xnrow,xncol);
BB2 = zeros(xnrow,xncol);
BBsrc = zeros(xnrow,xncol);
[sg,og] = sort(up./down);
for g = fliplr(og)
    gBi = Bi(Vn==g);
    gBj = Bj(Vn==g);
    gBk = gBi+gBj*xnrow-xnrow;
    ng = numel(gBk);
    if (ng==1) 
        continue
    end;
    gS = repmat(V(g),1,ng);
    cid = xC(gBk);
    P1(:) = -1;
    P2(:) = -1;
    P1src(:) = 0;
    P1(gBk) = 0;
    P1src(gBk) = 1:ng;
    J1(:) = 0;
    JJ1(:) = 0;
    J2(:) = 0;
    JJ2(:) = 0;
    JS(:) = inf;
    JJsrc(:) = 0;
    BB1(:) = 0;
    BB2(:) = 0;
    BBsrc(:) = 0;
    % find junctions
    for gg = 1:ng
        [P3,JJ3,BB3] = probe(P2,gBi(gg),gBj(gg),0);
        %subplot(2,1,1)
        %imagesc(P1); axis image;
        %subplot(2,1,2)
        %imagesc(P3); axis image;
        Jk = find((P1>-1)&(P3>0));
        Jk = Jk(P1(Jk)+P3(Jk)<JS(Jk));
        J1(Jk) = P1src(Jk);
        JJ1(Jk) = JJsrc(Jk);
        BB1(Jk) = BBsrc(Jk);
        J2(Jk) = gg;
        JJ2(Jk) = JJ3(Jk);
        BB2(Jk) = BB3(Jk);
        JS(Jk) = P1(Jk)+P3(Jk);
        rpl = (P3>=0)&((P1==-1)|(P3<=P1));
        P1(rpl) = P3(rpl);
        P1src(rpl) = gg;
        JJsrc(rpl) = JJ3(rpl);
        BBsrc(rpl) = BB3(rpl);
        %JJ1(rpl) = JJ3(rpl);
        
        %subplot(2,1,1)
        %imagesc(J1);
        %subplot(2,1,2)
        %imagesc(J2);
    end
    [mJS,kJS] = min(JS(:));
    J = find(~isinf(JS));
    while (mJS < V(g))
    %while ~isinf(mJS)
        [rJ,cJ] = ind2sub([xnrow,xncol],kJS);
        j1 = J1(kJS);
        j2 = J2(kJS);
        jj1 = JJ1(kJS);
        jj2 = JJ2(kJS);
        [jj1r,jj1c] = ind2sub([xnrow,xncol],jj1);
        [jj2r,jj2c] = ind2sub([xnrow,xncol],jj2);
        if (BB1(kJS)>0)
            JxC1(BB1(kJS)) = xC(BB1(kJS));
            JxC2(BB1(kJS)) = xC(gBi(j1),gBj(j1));
        end;
        if (BB2(kJS)>0)
            JxC1(BB2(kJS)) = xC(BB2(kJS));
            JxC2(BB2(kJS)) = xC(gBi(j1),gBj(j1));
        end;
        join2(gBi(j1),gBj(j1),rJ,cJ,jj1r,jj1c);
        %imagesc(xB);
        join2(rJ,cJ,gBi(j2),gBj(j2),jj2r,jj2c);
        %imagesc(xB);
        %title(mJS);
        if (BB1(kJS)>0)
            xB(BB1(kJS)) = inf;
        end;
        if (BB2(kJS)>0)
            xB(BB2(kJS)) = inf;
        end;
        jcid = cid(j1);
        cid(cid==cid(j2)) = jcid;
        rmv = (cid(J1(J))==jcid)&(cid(J2(J))==jcid);
        JS(J(rmv)) = inf;
        [mJS,kJS] = min(JS(:));
        %imagesc(xC);
    end;
    pcid = xC(gBk);
    gcid = unique(pcid);
    scid = zeros(1,numel(gcid));
    for i2 = 1:numel(gcid)
        scid(i2) = V(g)*sum(pcid==gcid(i2))-(sum(xC(:)==gcid(i2))-1)/2-25*(sum(JxC1(:)==gcid(i2))+sum(JxC2(:)==gcid(i2)));
    end;
    [gS(g),imscid] = max(scid);
    gID(g) = gcid(imscid);
    for i2 = 1:numel(gcid)
        if (i2 ~= imscid)|(scid(i2)<0)
            % delete these wires
            xB(xC==gcid(i2)) = 0;
            xC(xC==gcid(i2)) = 0;            
            xB(gBk) = V(g);
            xC(gBk) = pcid;
            xB(JxC1==gcid(i2)) = JxC2(JxC1==gcid(i2));
            xB(JxC2==gcid(i2)) = JxC1(JxC2==gcid(i2));
        end;
    end;
end;
% repick
%imagesc(xB);
for g = fliplr(og)
    gBi = Bi(Vn==g);
    gBj = Bj(Vn==g);
    gBk = gBi+gBj*xnrow-xnrow;
    ng = numel(gBk);
    if (ng==1) 
        continue
    end;
    
    % backup
    xB0 = xB;
    xC0 = xC;
    
    % reset
    xB(xC==gID(g)) = 0;
    xB(pin) = xB0(pin);
    xC(xC==gID(g)) = 0;
    xC(xB==V(g)) = max(xC(:))+(1:sum(xB(:)==V(g)));
    %imagesc(xB);
    
    % refigure
    gS = repmat(V(g),1,ng);
    cid = xC(gBk);
    P1(:) = -1;
    P2(:) = -1;
    P1src(:) = 0;
    P1(gBk) = 0;
    P1src(gBk) = 1:ng;
    J1(:) = 0;
    JJ1(:) = 0;
    J2(:) = 0;
    JJ2(:) = 0;
    JS(:) = inf;
    JJsrc(:) = 0;
    BB1(:) = 0;
    BB2(:) = 0;
    BBsrc(:) = 0;
    % find junctions
    for gg = 1:ng
        [P3,JJ3,BB3] = probe(P2,gBi(gg),gBj(gg),0);
        %subplot(2,1,1)
        %imagesc(P1); axis image;
        %subplot(2,1,2)
        %imagesc(P3); axis image;
        Jk = find((P1>-1)&(P3>0));
        Jk = Jk(P1(Jk)+P3(Jk)<JS(Jk));
        J1(Jk) = P1src(Jk);
        JJ1(Jk) = JJsrc(Jk);
        BB1(Jk) = BBsrc(Jk);
        J2(Jk) = gg;
        JJ2(Jk) = JJ3(Jk);
        BB2(Jk) = BB3(Jk);
        JS(Jk) = P1(Jk)+P3(Jk);
        rpl = (P3>=0)&((P1==-1)|(P3<=P1));
        P1(rpl) = P3(rpl);
        P1src(rpl) = gg;
        JJsrc(rpl) = JJ3(rpl);
        BBsrc(rpl) = BB3(rpl);
        %JJ1(rpl) = JJ3(rpl);
        
        %subplot(2,1,1)
        %imagesc(J1);
        %subplot(2,1,2)
        %imagesc(J2);
    end
    [mJS,kJS] = min(JS(:));
    J = find(~isinf(JS));
    %while (mJS < V(g))
    while ~isinf(mJS)
        [rJ,cJ] = ind2sub([xnrow,xncol],kJS);
        j1 = J1(kJS);
        j2 = J2(kJS);
        jj1 = JJ1(kJS);
        jj2 = JJ2(kJS);
        [jj1r,jj1c] = ind2sub([xnrow,xncol],jj1);
        [jj2r,jj2c] = ind2sub([xnrow,xncol],jj2);
        if (BB1(kJS)>0)
            JxC1(BB1(kJS)) = xC(BB1(kJS));
            JxC2(BB1(kJS)) = xC(gBi(j1),gBj(j1));
        end;
        if (BB2(kJS)>0)
            JxC1(BB2(kJS)) = xC(BB2(kJS));
            JxC2(BB2(kJS)) = xC(gBi(j1),gBj(j1));
        end;
        join2(gBi(j1),gBj(j1),rJ,cJ,jj1r,jj1c);
        %imagesc(xB);
        join2(rJ,cJ,gBi(j2),gBj(j2),jj2r,jj2c);
        %imagesc(xB);
        %title(mJS);
        if (BB1(kJS)>0)
            xB(BB1(kJS)) = inf;
        end;
        if (BB2(kJS)>0)
            xB(BB2(kJS)) = inf;
        end;
        jcid = cid(j1);
        cid(cid==cid(j2)) = jcid;
        rmv = (cid(J1(J))==jcid)&(cid(J2(J))==jcid);
        JS(J(rmv)) = inf;
        [mJS,kJS] = min(JS(:));
        %imagesc(xC);
    end;
    pcid = xC(gBk);
    gcid = unique(pcid);
    scid = zeros(1,numel(gcid));
    for i2 = 1:numel(gcid)
        scid(i2) = V(g)*sum(pcid==gcid(i2))-(sum(xC(:)==gcid(i2))-1)/2-25*(sum(JxC1(:)==gcid(i2))+sum(JxC2(:)==gcid(i2)));
    end;
    [gS2(g),imscid] = max(scid);
    for i2 = 1:numel(gcid)
        if (i2 ~= imscid)|(scid(i2)<0)
            % delete these wires
            xB(xC==gcid(i2)) = 0;
            xC(xC==gcid(i2)) = 0;            
            xB(gBk) = V(g);
            xC(gBk) = pcid;
            xB(JxC1==gcid(i2)) = JxC2(JxC1==gcid(i2));
            xB(JxC2==gcid(i2)) = JxC1(JxC2==gcid(i2));
        end;
    end;
    %if (gS(g)<gS2(g))
    %    % revert
    %    xB = xB0;
    %    xC = xC0;
    %    imagesc(xB);
    %end;
end;
W = convertSol(xB);


    %function scr = score(id)


    function join2(r1,c1,r2,c2,r3,c3)
        %v1 = xB(min(r1,r3):max(r1,r3),min(c1,c3):max(c1,c3));
        %v2 = xB(min(r3,r2):max(r3,r2),min(c3,c2):max(c3,c2));
        xB(min(r1,r3):max(r1,r3),min(c1,c3):max(c1,c3)) = xB(r1,c1);
        xB(min(r3,r2):max(r3,r2),min(c3,c2):max(c3,c2)) = xB(r1,c1);
        %if any(v1(2:end-1))|v2(2:end-1)
        %    sv1 = s2i(min(r1,r3):max(r1,r3),min(c1,c3):max(c1,c3));
        %    sv2 = s2i(min(r3,r2):max(r3,r2),min(c3,c2):max(c3,c2));
        %    br = [sv1(find(v1(2:end-1))+1) sv2(find(v2(2:end-1))+1)];
        %    xB(br) = inf;
        %end;
        c = xC(r2,c2);
        xC(min(r1,r3):max(r1,r3),min(c1,c3):max(c1,c3)) = xC(r1,c1);
        xC(min(r3,r2):max(r3,r2),min(c3,c2):max(c3,c2)) = xC(r1,c1);
        if (c > 0) 
            xC(xC==c) = xC(r1,c1);
            JxC1(JxC1==c) = JxC1(r1,c1);
            JxC2(JxC2==c) = JxC2(r1,c1);
        end;
    end


    function W = convertSol(xB)
        v = xB(2:2:end,1:2:end);
        h = xB(1:2:end,2:2:end);
        b = xB(1:2:end,1:2:end);
        [vi,vj] = find(v~=0);
        [hi,hj] = find(h~=0);
        [bi,bj] = find(isinf(b));
        W = [vi vj vi+1 vj; hi hj hi hj+1; bi bj bi bj];
    end


    function [P,JJ,BB] = probe(P,si,sj,base)
        JJ = zeros(xnrow,xncol);
        BB = zeros(xnrow,xncol);
        up = si+1-2*find([1;xB(1:2:si-2,sj)],1,'last'); % up
        left = sj+1-2*find([1 xB(si,1:2:sj-2)],1,'last'); % left
        down = -2+2*find([xB(si+2:2:end,sj);1],1); % down
        right = -2+2*find([xB(si,sj+2:2:end) 1],1); % right
        for ip = 2:2:up
            [P,JJ,BB] = probeH(P,JJ,BB,si-ip,sj,base+ip,sub2ind([xnrow,xncol],si-ip,sj));
        end;
        for ip = 2:2:left
            [P,JJ,BB] = probeV(P,JJ,BB,si,sj-ip,base+ip,sub2ind([xnrow,xncol],si,sj-ip));
        end;
        for ip = 2:2:down
            [P,JJ,BB] = probeH(P,JJ,BB,si+ip,sj,base+ip,sub2ind([xnrow,xncol],si+ip,sj));
        end;
        for ip = 2:2:right
            [P,JJ,BB] = probeV(P,JJ,BB,si,sj+ip,base+ip,sub2ind([xnrow,xncol],si,sj+ip));
        end;
        P(max(1,si-1),sj) = base+1;
        P(min(xnrow,si+1),sj) = base+1;
        P(si,max(1,sj-1)) = base+1;
        P(si,min(xncol,sj+1)) = base+1;
        JJbase = sub2ind([xnrow,xncol],si,sj);
        JJ(max(1,si-1),sj) = JJbase;
        JJ(min(xnrow,si+1),sj) = JJbase;
        JJ(si,max(1,sj-1)) = JJbase;
        JJ(si,min(xncol,sj+1)) = JJbase;
        P(si,sj) = base;
        JJ(si,sj) = JJbase;
        if (si-up-4>=1)&&~pin(si-up-2,sj)&(xB(si-up-4,sj)==0)
            [P,JJ,BB] = probeN(P,JJ,BB,si-up-4,sj,base+up+29,JJbase,s2i(si-up-2,sj));
        end;
        if (si+down+4<=xnrow)&&~pin(si+down+2,sj)&(xB(si+down+4,sj)==0)
            [P,JJ,BB] = probeS(P,JJ,BB,si+down+4,sj,base+down+29,JJbase,s2i(si+down+2,sj));
        end;
        if (sj-left-4>=1)&&~pin(si,sj-left-2)&(xB(si,sj-left-4)==0)
            [P,JJ,BB] = probeW(P,JJ,BB,si,sj-left-4,base+left+29,JJbase,s2i(si,sj-left-2));
        end;
        if (sj+right+4<xncol)&&~pin(si,sj+right+2)&(xB(si,sj+right+4)==0)
            [P,JJ,BB] = probeE(P,JJ,BB,si,sj+right+4,base+right+29,JJbase,s2i(si,sj+right+2));
        end;
end


    function [P,JJ,BB] = probeV(P,JJ,BB,si,sj,base,jj)
        vup = si+1-2*find([1;xB(1:2:si-2,sj)],1,'last'); % up
        vdown = -2+2*find([xB(si+2:2:end,sj);1],1); % down
        P((si-vup):2:si,sj) = [(base+vup):-2:base]';
        P(si:2:(si+vdown),sj) = [base:2:(base+vdown)]';
        JJ((si-vup):2:si,sj) = jj;
        JJ(si:2:(si+vdown),sj) = jj;
        if (si-vup-4>=1)&&~pin(si-vup-2,sj)&(xB(si-vup-4,sj)==0)
            [P,JJ,BB] = probeN(P,JJ,BB,si-vup-4,sj,base+vup+29,jj,s2i(si-vup-2,sj));
        end;
        if (si+vdown+4<=xnrow)&&~pin(si+vdown+2,sj)&(xB(si+vdown+4,sj)==0)
            [P,JJ,BB] = probeS(P,JJ,BB,si+vdown+4,sj,base+vdown+29,jj,s2i(si+vdown+2,sj));
        end;
    end


    function [P,JJ,BB] = probeH(P,JJ,BB,si,sj,base,jj)
        hleft = sj+1-2*find([1 xB(si,1:2:sj-2)],1,'last'); % left
        hright = -2+2*find([xB(si,sj+2:2:end) 1],1); % right
        P(si,(sj-hleft):2:sj) = [(base+hleft):-2:base];
        P(si,sj:2:(sj+hright)) = base:2:(base+hright);
        JJ(si,(sj-hleft):2:sj) = jj;
        JJ(si,sj:2:(sj+hright)) = jj;
        if (sj-hleft-4>=1)&&~pin(si,sj-hleft-2)&(xB(si,sj-hleft-4)==0)
            [P,JJ,BB] = probeW(P,JJ,BB,si,sj-hleft-4,base+hleft+29,jj,s2i(si,sj-hleft-2));
        end;
        if (sj+hright+4<xncol)&&~pin(si,sj+hright+2)&(xB(si,sj+hright+4)==0)
            [P,JJ,BB] = probeE(P,JJ,BB,si,sj+hright+4,base+hright+29,jj,s2i(si,sj+hright+2));
        end;
    end


    function [P,JJ,BB] = probeN(P,JJ,BB,si,sj,base,jj,bb)
        nup = si+1-2*find([1;xB(1:2:si-2,sj)],1,'last'); % up
        P((si-nup):2:si,sj) = [(base+nup):-2:base]';
        JJ((si-nup):2:si,sj) = jj;
        BB((si-nup):2:si,sj) = bb;
    end


    function [P,JJ,BB] = probeS(P,JJ,BB,si,sj,base,jj,bb)
        sdown = -2+2*find([xB(si+2:2:end,sj);1],1); % down
        P(si:2:(si+sdown),sj) = [base:2:(base+sdown)]';
        JJ(si:2:(si+sdown),sj) = jj;
        BB(si:2:(si+sdown),sj) = bb;
    end


    function [P,JJ,BB] = probeW(P,JJ,BB,si,sj,base,jj,bb)
        wleft = sj+1-2*find([1 xB(si,1:2:sj-2)],1,'last'); % left
        P(si,(sj-wleft):2:sj) = [(base+wleft):-2:base];
        JJ(si,(sj-wleft):2:sj) = jj;
        BB(si,(sj-wleft):2:sj) = bb;
    end


    function [P,JJ,BB] = probeE(P,JJ,BB,si,sj,base,jj,bb)
        eright = -2+2*find([xB(si,sj+2:2:end) 1],1); % right
        P(si,sj:2:(sj+eright)) = base:2:(base+eright);
        JJ(si,sj:2:(sj+eright)) = jj;
        BB(si,sj:2:(sj+eright)) = bb;
    end
end