Finish 2012-11-07 16:00:00 UTC

Black Widow

by Nicholas Howe

Status: Passed
Results: 3005 (cyc: 6, node: 1143)
CPU Time: 52.86
Score: 3074.09
Submitted at: 2012-11-02 15:14:46 UTC
Scored at: 2012-11-02 15:17:53 UTC

Current Rank: 1525th (Highest: 12th )

Comments
Please login or create a profile.
Code
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