function [moves, vine] = solver(board, limit)
if limit>15100
[moves,vine] = nick(board,limit);
else
[moves,vine] = alfi(board,limit);
end
end
function [moves, vine] = nick(board, limit)
% Added filter to crosseyed
%
% N. Howe
[vine1,scr1] = trickle2(board);
moves1 = [];
[moves2,vine2,scr2] = pusher(board,limit);
[nrow,ncol] = size(board);
ntot = nrow*ncol;
[moves3,vine3,scr3] = pusher(board',limit);
[moves4,vine4,scr4] = pusher(fliplr(board),limit);
[moves5,vine5,scr5] = pusher(fliplr(board'),limit);
[moves6,vine6,scr6] = dofilter(board,limit);
if (scr1==max([scr1 scr2 scr3 scr4 scr5 scr6]))
moves = moves1;
vine = vine1;
elseif (scr2==max([scr1 scr2 scr3 scr4 scr5 scr6]))
moves = moves2;
vine = vine2;
elseif (scr3==max([scr1 scr2 scr3 scr4 scr5 scr6]))
id = reshape(1:ntot,nrow,ncol);
id3 = id';
moves = id3(moves3);
vine = id3(vine3);
elseif (scr4==max([scr1 scr2 scr3 scr4 scr5 scr6]))
id = reshape(1:ntot,nrow,ncol);
id4 = fliplr(id);
moves = id4(moves4);
vine = id4(vine4);
elseif (scr5==max([scr1 scr2 scr3 scr4 scr5 scr6]))
id = reshape(1:ntot,nrow,ncol);
id5 = fliplr(id');
moves = id5(moves5);
vine = id5(vine5);
else
moves = moves6;
vine = vine6;
end;
end
function [moves, vine, scr] = pusher(board, limit)
moves = zeros(limit,2);
[nrow,ncol] = size(board);
if (limit < nrow)
scr = -inf;
moves = [];
vine = [];
return;
end;
ntot = nrow*ncol;
xg = repmat(1:ncol,nrow,1);
yg = repmat((1:nrow)',1,ncol);
pos = reshape(1:ntot,nrow,ncol);
fitloc = reshape(1:ntot,nrow,ncol);
fitloc(:,2:2:end) = flipud(fitloc(:,2:2:end));
[s,r] = sort(board(:),'descend');
rank = zeros(nrow,ncol);
rank(r) = 1:ntot;
remain = true(1,ntot);
nextfit = 1;
steps = 0;
cboard = board;
for i = 1:ntot
if remain(i)
[cy,cx] = find(rank==i);
tx = xg(fitloc(nextfit));
ty = yg(fitloc(nextfit));
while (cx~=tx)|(cy~=ty)
dx = sign(tx-cx);
dy = sign(ty-cy);
if (dx~=0)&(dy~=0)
if cboard(cy,cx+dx)<=cboard(cy+dy,cx)
dy = 0;
else
dx = 0;
end;
end;
if (rank(cy+dy,cx+dx)>0)
vx = cx+dx;
vy = cy+dy;
if (dx==0)
if (vx>1)&(cboard(vy,vx)>cboard(vy,vx-1))&((vx==ncol)|(cboard(vy,vx+1)>cboard(vy,vx-1)))
% move left
steps = steps+1;
if (steps > limit)
break;
end;
moves(steps,1) = pos(vy,vx);
moves(steps,2) = pos(vy,vx-1);
cboard(vy,vx-1) = cboard(vy,vx);
cboard(vy,vx) = 0;
if (rank(vy,vx-1)>0)
remain(rank(vy,vx-1)) = 0;
end;
rank(vy,vx-1) = rank(vy,vx);
rank(vy,vx) = 0;
%imagesc(cboard);
%axis image;
%drawnow;
elseif (vx < ncol)&(cboard(vy,vx)>cboard(vy,vx+1))
% move right
steps = steps+1;
if (steps > limit)
break;
end;
moves(steps,1) = pos(vy,vx);
moves(steps,2) = pos(vy,vx+1);
cboard(vy,vx+1) = cboard(vy,vx);
cboard(vy,vx) = 0;
if (rank(vy,vx+1)>0)
remain(rank(vy,vx+1)) = 0;
end;
rank(vy,vx+1) = rank(vy,vx);
rank(vy,vx) = 0;
%imagesc(cboard);
%axis image;
%drawnow;
else
remain(rank(cy+dy,cx+dx)) = 0;
end;
else
if (vy>1)&(cboard(vy,vx)>cboard(vy-1,vx))&((vy==nrow)|(cboard(vy+1,vx)>cboard(vy-1,vx)))
% move north
steps = steps+1;
if (steps > limit)
break;
end;
moves(steps,1) = pos(vy,vx);
moves(steps,2) = pos(vy-1,vx);
cboard(vy-1,vx) = cboard(vy,vx);
cboard(vy,vx) = 0;
if (rank(vy-1,vx)>0)
remain(rank(vy-1,vx)) = 0;
end;
rank(vy-1,vx) = rank(vy,vx);
rank(vy,vx) = 0;
%imagesc(cboard);
%axis image;
%drawnow;
elseif (vy < nrow)&(cboard(vy,vx)>cboard(vy+1,vx))
% move right
steps = steps+1;
if (steps > limit)
break;
end;
moves(steps,1) = pos(vy,vx);
moves(steps,2) = pos(vy+1,vx);
cboard(vy+1,vx) = cboard(vy,vx);
cboard(vy,vx) = 0;
if (rank(vy+1,vx)>0)
remain(rank(vy+1,vx)) = 0;
end;
rank(vy+1,vx) = rank(vy,vx);
rank(vy,vx) = 0;
%imagesc(cboard);
%axis image;
%drawnow;
else
remain(rank(cy+dy,cx+dx)) = 0;
end;
end;
end;
steps = steps+1;
if (steps > limit)
break;
end;
moves(steps,1) = pos(cy,cx);
moves(steps,2) = pos(cy+dy,cx+dx);
cboard(cy+dy,cx+dx) = s(i); %cboard(cy,cx);
cboard(cy,cx) = 0;
%rank(cy+dy,cx+dx) = i; %cboard(cy,cx);
rank(cy,cx) = 0;
cx = cx+dx;
cy = cy+dy;
%imagesc(cboard);
%axis image;
%drawnow;
end;
if (steps > limit)
break;
end;
rank(ty,tx) = i;
nextfit = nextfit+1;
end;
end;
moves = moves(1:min(steps,limit),:);
[vine,scr] = trickle2(cboard);
end
function [vine,scr] = trickle2(board)
% Second effort
% Somewhat better handling of ties.
%
%N. Howe
moves = [];
[nrow,ncol] = size(board);
xg = repmat(1:ncol,nrow,1);
yg = repmat((1:nrow)',1,ncol);
limbo = nrow*ncol+1;
vals = [board(:)' -inf];
%up = [zeros(1,ncol); board(1:end-1,:)];
%down = [board(2:end,:);zeros(1,ncol)];
%left = [zeros(nrow,1); board(:,1:end-1)];
%right = [board(:,2:end);zeros(nrow,1)];
ind = reshape(1:nrow*ncol,nrow,ncol);
up = [repmat(limbo,1,ncol); ind(1:end-1,:)];
down = [ind(2:end,:);repmat(limbo,1,ncol)];
left = [repmat(limbo,nrow,1) ind(:,1:end-1)];
right = [ind(:,2:end) repmat(limbo,nrow,1)];
nbr = [[up(:) down(:) left(:) right(:)]; limbo limbo limbo limbo];
%vnbr = vals(nbr);
%big = repmat(board(:),1,4)<=vnbr;
cum = vals;
next = zeros(1,nrow*ncol+1);
next(limbo) = -1;
levels = unique(vals);
done(nrow*ncol+1) = true;
for lvl = levels(end:-1:2)
lvli = find(board==lvl);
lvln = numel(lvli);
for ii = 1:lvln
moves = (vals(nbr(lvli,:))>=lvl).*done(nbr(lvli,:));
[heur,jheur] = max(cum(nbr(lvli,:)).*moves,[],2);
heur(done(lvli)) = 0;
[mheur,mheuri] = max(heur);
if (mheur==0)
done(lvli) = true;
break;
end;
i = lvli(mheuri);
cum(i) = cum(i)+mheur;
next(i) = nbr(i,jheur(mheuri));
done(i) = true;
end;
end;
vine = zeros(1,nrow*ncol);
[scr,vine(1)] = max(cum);
vlen = 1;
while (next(vine(vlen))>0)
vine(vlen+1) = next(vine(vlen));
vlen = vlen+1;
end;
% add on at end
i = vine(vlen);
done(i) = false;
lvl = board(i);
mv = (vals(nbr(i,:))==lvl).*done(nbr(i,:));
while any(mv)
vine(vlen+1) = nbr(i,find(mv,1));
vlen = vlen+1;
i = vine(vlen);
done(i) = false;
lvl = board(i);
mv = (vals(nbr(i,:))==lvl).*done(nbr(i,:));
end;
vine = vine(1:vlen);
% now do augmentation
j = 1;
used = false(nrow,ncol);
used(vine) = true;
while (j < vlen)
i1 = vine(j);
i2 = vine(j+1);
x1 = xg(i1);
x2 = xg(i2);
y1 = yg(i1);
%y2 = yg(i2);
if (x1==x2)
% vertical move
if (x1>1)&&(~used(i1-nrow))&&(~used(i2-nrow))&&(board(i1)<=board(i1-nrow))&&(board(i1-nrow)<=board(i2-nrow))&&(board(i2-nrow)<=board(i2))
% augment to the left
vlen = vlen+2;
vine = [vine(1:j) i1-nrow i2-nrow vine(j+1:end)];
used(i1-nrow) = true;
used(i2-nrow) = true;
elseif (x1<ncol)&&(~used(i1+nrow))&&(~used(i2+nrow))&&(board(i1)<=board(i1+nrow))&&(board(i1+nrow)<=board(i2+nrow))&&(board(i2+nrow)<=board(i2))
% augment to the right
vlen = vlen+2;
vine = [vine(1:j) i1+nrow i2+nrow vine(j+1:end)];
used(i1+nrow) = true;
used(i2+nrow) = true;
else
% no augmentation; move on
j = j+1;
end;
else
% horizontal move
if (y1>1)&&(~used(i1-1))&&(~used(i2-1))&&(board(i1)<=board(i1-1))&&(board(i1-1)<=board(i2-1))&&(board(i2-1)<=board(i2))
% augment to the north
vlen = vlen+2;
vine = [vine(1:j) i1-1 i2-1 vine(j+1:end)];
used(i1-1) = true;
used(i2-1) = true;
elseif (y1<nrow)&&(~used(i1+1))&&(~used(i2+1))&&(board(i1)<=board(i1+1))&&(board(i1+1)<=board(i2+1))&&(board(i2+1)<=board(i2))
% augment to the south
vlen = vlen+2;
vine = [vine(1:j) i1+1 i2+1 vine(j+1:end)];
used(i1+1) = true;
used(i2+1) = true;
else
% no augmentation; move on
j = j+1;
end;
end;
end;
scr = sum(board(vine));
end
function [moves,vine,scr] = dofilter(board,limit)
[nrow,ncol] = size(board);
ntot = nrow*ncol;
moves = zeros(limit,2);
nmove = 0;
pos = reshape(1:ntot,nrow,ncol);
m1 = median(board,1);
m2 = median(board,2);
if (sum(sign(diff(m1)))>0.67*ncol)||(sum(sign(diff(m2)))>0.67*nrow)||(sum(sign(diff(m2)))<-0.67*nrow)||(sum(sign(diff(m2)))<-0.67*nrow)
d = diff(board,1,2);
mb = median(board,2);
md = median(d,2);
rmb = repmat(mb,1,ncol);
rmd = repmat(md,1,ncol);
guess = rmb+kron(md,-(ncol-1)/2:(ncol-1)/2);
bad = (abs(board-guess)>20*abs(rmd));
for i = nrow:-1:1
offset = 0;
for j = ceil(ncol/2):ncol
if bad(i,j)
offset = offset+1;
else
moves(nmove+1:nmove+offset,:) = [pos(i,j:-1:j-offset+1)' pos(i,j-1:-1:j-offset)'];
nmove = nmove+offset;
if (nmove>=limit)
break;
end;
end;
end;
offset = 0;
for j = ceil(ncol/2):-1:1
if bad(i,j)
offset = offset+1;
else
moves(nmove+1:nmove+offset,:) = [pos(i,j:1:j+offset-1)' pos(i,j+1:1:j+offset)'];
nmove = nmove+offset;
if (nmove>=limit)
break;
end;
end;
end;
if (nmove>=limit)
break;
end;
end;
moves = moves(1:min(nmove,limit),:);
for i = 1:size(moves,1)
board(moves(i,2)) = board(moves(i,1));
board(moves(i,1)) = 0;
%imagesc(board);
%axis image;
%drawnow;
end;
[vine,scr] = trickle2(board);
else
scr = -inf;
moves = [];
vine = [];
end;
end
function [moves, vine] = alfi(board, limit)
[moves, vine, gain,newboard] = solver1(board, limit);
% disp(gain)
% imagesc(newboard);hold on;[i1,i2]=ind2sub(size(board),vine);plot(i2,i1,'k.-'); hold off;
end
function [moves, vine, gain, board] = solver1(board, limit)
moves=[];
ab=accumarray(1+board(:),1);
cab=flipud(cumsum(flipud(ab)))-ab(end);
[m,n]=size(board);
board=[zeros(1,n+2);[zeros(m,1),board,zeros(m,1)];zeros(1,n+2)];
[m,n]=size(board);
boardab=cab(1+board);
if max(ab(2:end))>1,
% compute same-valued clusters
[BoardCluster,ClusterValue,IdxList,IdxSegments,Nclusters]=connected(board);
ClusterSize=diff(IdxSegments);
ClusterNeighb=neighb(BoardCluster,board,Nclusters);
% search between-clusters
ClustersOrder = bellman(ClusterNeighb,ClusterValue'.*sqrt(ClusterSize'));
% search within-clusters
iC=bellman_postprocess(ClustersOrder,IdxList,IdxSegments,BoardCluster);
else
% search between-pixels
iC = bellman_pixel(board,limit,boardab);
end
% post-processing moves
if limit>0
[board,iC,moves,limit] = postprocess(board,iC,moves,limit,boardab);
[iC1,iC2]=ind2sub([m,n],moves);
moves=sub2ind([m-2,n-2],iC1-1,iC2-1);
end
[iC1,iC2]=ind2sub([m,n],iC);
vine=sub2ind([m-2,n-2],iC1-1,iC2-1);
vine=vine(end:-1:1);
board=board(2:end-1,2:end-1);
gain=sum(board(vine));
end
function [board,tiC,moves,limit] = postprocess(board,tiC,moves,limit,boardab)
[m,n]=size(board);
nC=numel(tiC);
if nC>1&&limit>0, %%%
% grow laterally
d=abs(diff(tiC))==1;
E=[m*~d+d; m*d+~d];
BorderUp={repmat(1+(0:n-1)*m,[m,1]),repmat((1:m)',[1,n])};
BorderDown={repmat(m+(0:n-1)*m,[m,1]),repmat(m*(n-1)+(1:m)',[1,n])};
tP=board>0;
tP(tiC)=false;
for n1=1:1e3,
[board,moves,tiC,tP,E,limit,ok]=postprocess_lateral(board,moves,tiC,tP,E,limit,BorderDown,BorderUp,boardab);
if ~ok, break; end
end
end
if limit>0
% grow end-points
tP=board>0;
tP(tiC)=false;
for n1=1:1e3,
[board,moves,tiC,tP,limit,ok]=postprocess_endpoint(board,moves,tiC,tP,limit,1);
if ~ok, break; end
end
for n1=1:1e3,
[board,moves,tiC,tP,limit,ok]=postprocess_endpoint(board,moves,tiC,tP,limit,2);
if ~ok, break; end
end
end
end
function [board,moves,tiC,tP,E,limit,ok]=postprocess_lateral(board,moves,tiC,tP,E,limit,BorderDown,BorderUp,boardab)
[m,n]=size(board);
ok=0;
% D=Dijkstra(tP,tiC(1));
% maskout=tP&~isinf(D)&board>board(tiC(1));
% nmaskout=sum(D(maskout));
for ndir=1:4
K=size(tiC,2);
switch(ndir)
case 1, idx1=find(tP(tiC(1:2:K-1)+E(2,1:2:K-1))&tP(tiC(2:2:K)+E(2,1:2:K-1)));
case 2, idx1=find(tP(tiC(1:2:K-1)-E(2,1:2:K-1))&tP(tiC(2:2:K)-E(2,1:2:K-1)));
case 3, idx1=find(tP(tiC(2:2:K-1)+E(2,2:2:K-1))&tP(tiC(3:2:K)+E(2,2:2:K-1)));
case 4, idx1=find(tP(tiC(2:2:K-1)-E(2,2:2:K-1))&tP(tiC(3:2:K)-E(2,2:2:K-1)));
end
%for n2=1:numel(idx1),
for n2=numel(idx1):-1:1, %%%
switch(ndir)
case 1,
tempA=tiC(2*idx1(n2)-1)+E(2,2*idx1(n2)-1):E(2,2*idx1(n2)-1):BorderDown{1+(E(2,2*idx1(n2)-1)==m)}(tiC(2*idx1(n2)-1));
tempB=tiC(2*idx1(n2))+E(2,2*idx1(n2)-1):E(2,2*idx1(n2)-1):BorderDown{1+(E(2,2*idx1(n2)-1)==m)}(tiC(2*idx1(n2)));
case 2,
tempA=tiC(2*idx1(n2)-1)-E(2,2*idx1(n2)-1):-E(2,2*idx1(n2)-1):BorderUp{1+(E(2,2*idx1(n2)-1)==m)}(tiC(2*idx1(n2)-1));
tempB=tiC(2*idx1(n2))-E(2,2*idx1(n2)-1):-E(2,2*idx1(n2)-1):BorderUp{1+(E(2,2*idx1(n2)-1)==m)}(tiC(2*idx1(n2)));
case 3,
tempA=tiC(2*idx1(n2))+E(2,2*idx1(n2)):E(2,2*idx1(n2)):BorderDown{1+(E(2,2*idx1(n2))==m)}(tiC(2*idx1(n2)));
tempB=tiC(2*idx1(n2)+1)+E(2,2*idx1(n2)):E(2,2*idx1(n2)):BorderDown{1+(E(2,2*idx1(n2))==m)}(tiC(2*idx1(n2)+1));
case 4,
tempA=tiC(2*idx1(n2))-E(2,2*idx1(n2)):-E(2,2*idx1(n2)):BorderUp{1+(E(2,2*idx1(n2))==m)}(tiC(2*idx1(n2)));
tempB=tiC(2*idx1(n2)+1)-E(2,2*idx1(n2)):-E(2,2*idx1(n2)):BorderUp{1+(E(2,2*idx1(n2))==m)}(tiC(2*idx1(n2)+1));
end
idxZ=find(~tP(tempA),1,'first');
%idxZ=find(~tP(tempA)|maskout(tempA),1,'first');
switch(ndir)
case {1,2},
%idxA=find(board(tiC(2*idx1(n2)-1))>=board(tempA)&board(tempA)>=board(tiC(2*idx1(n2))),1,'first');
[nill,idxA]=max((1./(abs(board(tiC(2*idx1(n2)-1))-board(tempA))+1)).*(board(tiC(2*idx1(n2)-1))>=board(tempA)&board(tempA)>=board(tiC(2*idx1(n2)))));
if ~nill,idxA=[]; end
case {3,4},
%idxA=find(board(tiC(2*idx1(n2)))>=board(tempA)&board(tempA)>=board(tiC(2*idx1(n2)+1)),1,'first');
[nill,idxA]=max((1./(abs(board(tiC(2*idx1(n2)))-board(tempA))+1)).*(board(tiC(2*idx1(n2)))>=board(tempA)&board(tempA)>=board(tiC(2*idx1(n2)+1))));
if ~nill,idxA=[]; end
end
if ~isempty(idxA)&&idxA<idxZ
idxZ=find(~tP(tempB),1,'first');
%idxZ=find(~tP(tempB)|maskout(tempA),1,'first');
switch(ndir)
case {1,2},
idxB=find(board(tempA(idxA))>=board(tempB)&board(tempB)>=board(tiC(2*idx1(n2))),1,'first');
case {3,4},
idxB=find(board(tempA(idxA))>=board(tempB)&board(tempB)>=board(tiC(2*idx1(n2)+1)),1,'first');
end
if ~isempty(idxB)&&idxB<idxZ&&idxA+idxB-2<=limit,%-nmaskout,%.0*boardab(tiC(1))*(n+m)/2, %%%
ok=1;
switch(ndir)
case 1, tiC=[tiC(1:2*idx1(n2)-1),tiC(2*idx1(n2)-1)+E(2,2*idx1(n2)-1), tiC(2*idx1(n2))+E(2,2*idx1(n2)-1), tiC(2*idx1(n2):end)];
case 2, tiC=[tiC(1:2*idx1(n2)-1),tiC(2*idx1(n2)-1)-E(2,2*idx1(n2)-1), tiC(2*idx1(n2))-E(2,2*idx1(n2)-1), tiC(2*idx1(n2):end)];
case 3, tiC=[tiC(1:2*idx1(n2)),tiC(2*idx1(n2))+E(2,2*idx1(n2)), tiC(2*idx1(n2)+1)+E(2,2*idx1(n2)), tiC(2*idx1(n2)+1:end)];
case 4, tiC=[tiC(1:2*idx1(n2)),tiC(2*idx1(n2))-E(2,2*idx1(n2)), tiC(2*idx1(n2)+1)-E(2,2*idx1(n2)), tiC(2*idx1(n2)+1:end)];
end
E=[E(:,1:2*idx1(n2)-2),E([2,1],2*idx1(n2)-1), E([1,2],2*idx1(n2)-1), E([2,1],2*idx1(n2)-1), E(:,2*idx1(n2):end)];
tP(tiC)=false;
newmoveA=[tempA(2:idxA)',tempA(1:idxA-1)'];
newmoveB=[tempB(2:idxB)',tempB(1:idxB-1)'];
board(newmoveA(:,2))=board(tempA(idxA));
board(newmoveB(:,2))=board(tempB(idxB));
board(newmoveA(:,1))=0;
board(newmoveB(:,1))=0;
moves=[moves;flipud(newmoveA);flipud(newmoveB)];
limit=limit-size(newmoveA,1)-size(newmoveB,1);
%break;
end
end
end
end
end
function [board,moves,tiC,tP,limit,ok]=postprocess_endpoint(board,moves,tiC,tP,limit,type)
[m,n]=size(board);
ok=1;
switch(type)
case 1
mask=tP&board<2*board(tiC(1));
D=Dijkstra(mask,tiC(1));
idx1=find(mask&D-1<=limit&board>=board(tiC(1)));
if isempty(idx1),
mask=tP;
D=Dijkstra(mask,tiC(1));
idx1=find(mask&D-1<=limit&board>=board(tiC(1)));
if isempty(idx1),ok=0;end
end
if ok
[nill,idx2]=min(board(idx1)+1e-3*D(idx1)); %%%t3
idx0=idx1(idx2);
limit=limit-(D(idx0)-1);
for n2=D(idx0)-1:-1:1
if D(idx0+1)==n2, tidx0=idx0+1;
elseif D(idx0-1)==n2, tidx0=idx0-1;
elseif D(idx0+m)==n2, tidx0=idx0+m;
elseif D(idx0-m)==n2, tidx0=idx0-m;
end
moves=[moves; idx0, tidx0];
board(tidx0)=board(idx0);
board(idx0)=0;
idx0=tidx0;
end
tiC=[idx0,tiC];
tP(idx0)=false;
end
case 2
mask=tP;
D=Dijkstra(mask,tiC(end));
idx1=find(mask&D-1<=limit&board<=board(tiC(end))&board>0);
if isempty(idx1),ok=0;end
if ok
[nill,idx2]=min(-board(idx1)+1e-3*D(idx1)); %%%t3
idx0=idx1(idx2);
limit=limit-(D(idx0)-1);
for n2=D(idx0)-1:-1:1
if D(idx0+1)==n2, tidx0=idx0+1;
elseif D(idx0-1)==n2, tidx0=idx0-1;
elseif D(idx0+m)==n2, tidx0=idx0+m;
elseif D(idx0-m)==n2, tidx0=idx0-m;
end
moves=[moves; idx0, tidx0];
board(tidx0)=board(idx0);
board(idx0)=0;
idx0=tidx0;
end
tiC=[tiC,idx0];
tP(idx0)=false;
end
end
end
function iC = bellman_pixel(board,limit,boardab)
[m,n]=size(board);
cD=board;
IDX=zeros(m,n);
touched=board>0;
valid=touched;
while any(touched(:))
touch=find(touched);
touched(touch)=false;
for dx=[1,-1,m,-m],
cDnew=board(touch+dx)+cD(touch)+0*min(limit,boardab(touch+dx)*(n+m)/2).*board(touch+dx)/((n+m)/2); %%%
idx=find(board(touch+dx)>board(touch)&cDnew>cD(touch+dx));
touched(touch(idx)+dx)=valid(touch(idx)+dx);
cD(touch(idx)+dx)=cDnew(idx);
IDX(touch(idx)+dx)=touch(idx);
end
end
[nill,idx]=max(cD(:));
iC=zeros(1,m*n);
for n1=1:m*n,
if idx>0
iC(n1)=idx;
idx=IDX(idx);
else break;
end
end
iC=iC(1:n1-1);
end
function iC = bellman(C,D)
N=size(C,1);
c=cell(1,N); for n1=1:N,c{n1}=find(C(:,n1)); end
%C=full(C);
IDX=zeros(N,1);
cD=D;
touched=true(1,N);
while any(touched)
for touch=find(touched),
touched(touch)=false;
idx=c{touch};%find(C(:,touch));
cDnew=D(idx)+cD(touch);
idx2=find(cD(idx)<cDnew);
if ~isempty(idx2)
cD(idx(idx2))=cDnew(idx2);
touched(idx(idx2))=true;
IDX(idx(idx2))=touch;
end
end
end
[nill,idx]=max(cD);
iC=zeros(N,1);
for n1=1:N,
if idx>0
iC(n1)=idx;
idx=IDX(idx);
else break;
end
end
iC=iC(1:n1-1);
end
function iC=bellman_postprocess(ClustersOrder,IdxList,IdxSegments,BoardCluster)
[m,n]=size(BoardCluster);
iC=[];
for n1=1:numel(ClustersOrder)
idx=IdxList(IdxSegments(ClustersOrder(n1)):IdxSegments(ClustersOrder(n1)+1)-1);
if numel(idx)==1, iC=[iC,idx(1)];
else
% longest shortest-path
ThisCluster=false(m,n);
ThisCluster(idx)=true;
D=Dijkstra(ThisCluster,iC);
if n1==numel(ClustersOrder)
temp=D(ThisCluster);
else
temp=D(ThisCluster).*(BoardCluster([false(1,n);ThisCluster(1:m-1,:)])==ClustersOrder(n1+1)|BoardCluster([ThisCluster(2:m,:);false(1,n)])==ClustersOrder(n1+1)|BoardCluster([false(m,1),ThisCluster(:,1:n-1)])==ClustersOrder(n1+1)|BoardCluster([ThisCluster(:,2:n),false(m,1)])==ClustersOrder(n1+1));
end
[nill,idx0]=max(temp(:));
idx0=idx(idx0);
E=zeros(2,D(idx0)-1);
tiC=zeros(1,D(idx0));
tiC(end)=idx0;
for n2=D(idx0)-1:-1:1
if D(idx0+1)==n2, idx0=idx0+1; E(1,n2)=1;E(2,n2)=m;
elseif D(idx0-1)==n2, idx0=idx0-1; E(1,n2)=1;E(2,n2)=m;
elseif D(idx0+m)==n2, idx0=idx0+m; E(1,n2)=m;E(2,n2)=1;
elseif D(idx0-m)==n2, idx0=idx0-m; E(1,n2)=m;E(2,n2)=1;
end
tiC(n2)=idx0;
end
% now grow it
tP=ThisCluster;
tP(tiC)=false;
ok=1;
while ok
ok=0;
K=size(tiC,2);
idx1=find(tP(tiC(1:2:K-1)+E(2,1:2:K-1))&tP(tiC(2:2:K)+E(2,1:2:K-1)));
for n2=numel(idx1):-1:1
ok=1;
tiC=[tiC(1:2*idx1(n2)-1),tiC(2*idx1(n2)-1)+E(2,2*idx1(n2)-1), tiC(2*idx1(n2))+E(2,2*idx1(n2)-1), tiC(2*idx1(n2):end)];
E=[E(:,1:2*idx1(n2)-2),E([2,1],2*idx1(n2)-1), E([1,2],2*idx1(n2)-1), E([2,1],2*idx1(n2)-1), E(:,2*idx1(n2):end)];
tP(tiC)=false;
end
K=size(tiC,2);
idx1=find(tP(tiC(1:2:K-1)-E(2,1:2:K-1))&tP(tiC(2:2:K)-E(2,1:2:K-1)));
for n2=numel(idx1):-1:1
ok=1;
tiC=[tiC(1:2*idx1(n2)-1),tiC(2*idx1(n2)-1)-E(2,2*idx1(n2)-1), tiC(2*idx1(n2))-E(2,2*idx1(n2)-1), tiC(2*idx1(n2):end)];
E=[E(:,1:2*idx1(n2)-2),E([2,1],2*idx1(n2)-1), E([1,2],2*idx1(n2)-1), E([2,1],2*idx1(n2)-1), E(:,2*idx1(n2):end)];
tP(tiC)=false;
end
K=size(tiC,2);
idx1=find(tP(tiC(2:2:K-1)+E(2,2:2:K-1))&tP(tiC(3:2:K)+E(2,2:2:K-1)));
for n2=numel(idx1):-1:1
ok=1;
tiC=[tiC(1:2*idx1(n2)),tiC(2*idx1(n2))+E(2,2*idx1(n2)), tiC(2*idx1(n2)+1)+E(2,2*idx1(n2)), tiC(2*idx1(n2)+1:end)];
E=[E(:,1:2*idx1(n2)-1),E([2,1],2*idx1(n2)), E([1,2],2*idx1(n2)), E([2,1],2*idx1(n2)), E(:,2*idx1(n2)+1:end)];
tP(tiC)=false;
end
K=size(tiC,2);
idx1=find(tP(tiC(2:2:K-1)-E(2,2:2:K-1))&tP(tiC(3:2:K)-E(2,2:2:K-1)));
for n2=numel(idx1):-1:1
ok=1;
tiC=[tiC(1:2*idx1(n2)),tiC(2*idx1(n2))-E(2,2*idx1(n2)), tiC(2*idx1(n2)+1)-E(2,2*idx1(n2)), tiC(2*idx1(n2)+1:end)];
E=[E(:,1:2*idx1(n2)-1),E([2,1],2*idx1(n2)), E([1,2],2*idx1(n2)), E([2,1],2*idx1(n2)), E(:,2*idx1(n2)+1:end)];
tP(tiC)=false;
end
end
iC=[iC,tiC];
end
end
end
function B = neighb(A,V,nA)
[m,n] = size(A);
B=sparse(A(:,1:n-1),A(:,2:n),double(V(:,1:n-1)>V(:,2:n)),nA,nA)+sparse(A(:,2:n),A(:,1:n-1),double(V(:,2:n)>V(:,1:n-1)),nA,nA)+sparse(A(1:m-1,:),A(2:m,:),double(V(1:m-1,:)>V(2:m,:)),nA,nA)+sparse(A(2:m,:),A(1:m-1,:),double(V(2:m,:)>V(1:m-1,:)),nA,nA);
B=B>0;
end
function [B,C,p,r,nc] = connected(A)
[m,n] = size(A);
N = m*n ;
K = reshape (1:N, m, n) ;
V = A(:,1:n-1) == A(:,2:n);
H = A(1:m-1,:) == A(2:m,:);
G = sparse(K([V,false(m,1)]),K([false(m,1),V]),1,N,N)+sparse(K([H; false(1,n)]),K([false(1,n); H]),1,N,N);
G = G + G' + speye(N);
[p q r] = dmperm(G);
nc = numel(r) - 1;
C = A(p(r(1:nc)));
B = ones(m, n);
for a = 2:nc
B(p(r(a):r(a+1)-1)) = a;
end
end
function D = Dijkstra(A,i1)
[m,n]=size(A);
mn=m*n;
D=inf(m,n);
D(~A)=0;
if isempty(i1),
i1=find(A,1,'first');
mode=0;
else
i1=i1(end);
mode=1;
% if A(i1+1); i1=i1+1;
% elseif A(i1-1); i1=i1-1;
% elseif A(i1+m); i1=i1+m;
% elseif A(i1-m); i1=i1-m;
% end
end
D(i1)=1;
for n1=2:mn,
idx1=D(i1+1)>n1;
idx2=D(i1-1)>n1;
idx3=D(i1+m)>n1;
idx4=D(i1-m)>n1;
%i1=unique([i1(idx1)+1,i1(idx2)-1,i1(idx3)+m,i1(idx4)-m]);
X=false(m,n);X(i1(idx1)+1)=true;X(i1(idx2)-1)=true;X(i1(idx3)+m)=true;X(i1(idx4)-m)=true;
i1=find(X);
if isempty(i1), break; end
D(i1)=n1;
end
if mode, D(D>0)=D(D>0)-1; end
end
|
