function W=solver(B);

W=zeros(0,4);
sB=size(B);
B0=zeros(sB);
idx0=find(B(:));
[b,I,J]=unique(B(idx0));

L0=zeros(size(b));
for n1=1:length(b),
	idx10{n1}=find(J==n1);
	L0(n1)=length(idx10{n1});
end
L0(L0==1)=0;
[nill,sortN]=sort(L0.*b);

for nsortN=length(sortN):-1:1+sum(L0==0),
	N=sortN(nsortN);
	L1=L0(N);
	idx1=idx10{N};

	if L1>1,
		R=cell(2*L1-1,2*L1-1);C=R;Q=R; % connecting paths
		D=zeros(2*L1-1,2*L1-1); D(:)=inf; % distances
		E=zeros(2*L1-1,2*L1-1); % gain
		IDX=zeros(2*L1-1,1); IDX(1:L1)=1; % valid clusters
		[r,c]=ind2sub(sB,idx0(idx1));
		for n1=1:L1, RR{n1}=r(n1); CC{n1}=c(n1); QQ{n1}=[]; end % clusters
		EE=zeros(2*L1-1,1); EE(1:L1)=b(N);
		WW=cell(2*L1-1,1);
		
		% first pass
		LL=2/L1/(L1-1);
		e=0; idx=[];
		for n1=1:L1, 
			for n2=n1+1:L1, 
				[R{n1,n2},C{n1,n2},Q{n1,n2},D(n1,n2)]=findpath(B,r(n1),c(n1),r(n2),c(n2)); 
				E(n1,n2)=EE(n1)+EE(n2)-D(n1,n2)+1;
				if E(n1,n2)>e || (E(n1,n2)==e&&e>0&&rand<LL), e=E(n1,n2); idx=[n1,n2]; end
			end
		end
		K=L1; 
		if e>0, % something
			while ~isempty(idx),
				% new cluster
				K=K+1; 
				Rnow=R{idx(1),idx(2)}; Cnow=C{idx(1),idx(2)}; Qnow=Q{idx(1),idx(2)}; 
				RR{K}=cat(1,RR{idx(1)},Rnow((2:end-1)'),RR{idx(2)});
				CC{K}=cat(1,CC{idx(1)},Cnow((2:end-1)'),CC{idx(2)});
				QQ{K}=cat(1,QQ{idx(1)},Qnow,QQ{idx(2)});
				EE(K)=EE(idx(1))+EE(idx(2))-D(idx(1),idx(2))+1;
				WW{K}=cat(1,WW{idx(1)},WW{idx(2)},[Rnow(1:end-1),Cnow(1:end-1),Rnow(2:end),Cnow(2:end)]);
				IDX(idx)=0;
				n1=K; N2=find(IDX);
				IDX(K)=1;
				if ~isempty(N2), 
					for n2=N2(:)', 
						if idx(1)<n2, mi1=idx(1);ma1=n2; else, mi1=n2;ma1=idx(1); end
						if idx(2)<n2, mi2=idx(2);ma2=n2; else, mi2=n2;ma2=idx(2); end
						if D(mi1,ma1)<D(mi2,ma2)
							R{n2,n1}=R{mi1,ma1};
							C{n2,n1}=C{mi1,ma1};
							D(n2,n1)=D(mi1,ma1);
							E(n2,n1)=EE(n1)+EE(n2)-D(n2,n1)+1;
						else, 
							R{n2,n1}=R{mi2,ma2};
							C{n2,n1}=C{mi2,ma2};
							D(n2,n1)=D(mi2,ma2);
							E(n2,n1)=EE(n1)+EE(n2)-D(n2,n1)+1;
						end
						if D(idx(1),idx(2))>2,
							% recompute distances
							LL=1/length(N2);
							[Rt,Ct,Qt,Dt]=findpath(B,Rnow((2:end-1)'),Cnow((2:end-1)'),RR{n2},CC{n2},1); 
							if Dt<D(n2,n1), R{n2,n1}=Rt;C{n2,n1}=Ct;Q{n2,n1}=Qt;D(n2,n1)=Dt; end
							%[R{n2,n1},C{n2,n1},Q{n2,n1},D(n2,n1)]=findpath(B,RR{n1},CC{n1},RR{n2},CC{n2},1); 
							E(n2,n1)=EE(n1)+EE(n2)-D(n2,n1)+1;
						end
					end
					
					idx=[];
					IDX1=find(IDX);N1=length(IDX1);d=D(IDX1,IDX1);
					[sortd,idxd]=sort(d(:)); % closest path
					[r1,c1]=ind2sub([N1,N1],idxd(1:N1*(N1-1)/2));
					idxfirst=find(E(IDX1(r1)+size(E,1)*(IDX1(c1)-1))>max(EE(IDX1(r1)),EE(IDX1(c1))));  % first increasing gain
					if ~isempty(idxfirst), idx=[IDX1(r1(idxfirst(1))),IDX1(c1(idxfirst(1)))]; end
				else, idx=[]; end
			end
		end
		if K>L1,
			for n1=1:K,EE(n1)=EE(n1)-25*sum(B0(RR{n1}+sB(1)*(CC{n1}-1))>0);end
			[maxEE,K]=max(EE);
			if maxEE>0,
				W=cat(1,W,WW{K});
				%for n3=1:size(QQ{K},1), 
				%	B(QQ{K}(n3,1),QQ{K}(n3,2))=-1; 
				%end
				for n3=1:length(RR{K}), 
					B(RR{K}(n3),CC{K}(n3))=-1; 
					if B0(RR{K}(n3),CC{K}(n3))>0, W=cat(1,W,[RR{K}(n3),CC{K}(n3),RR{K}(n3),CC{K}(n3)]);
					else, B0(RR{K}(n3),CC{K}(n3))=1; end
				end
			end
		end
	end
end

function [R,C,Q,L]=findpath(B,r1,c1,r2,c2,xtra);
% finds straight path from r1 c1 to r2 c2 in map B
if nargin<6, xtra=0; end
Q=[];
L1=length(r1); L2=length(r2);
A=abs(r1(:,ones(1,L2))-r2(:,ones(1,L1))')+abs(c1(:,ones(1,L2))-c2(:,ones(1,L1))');
[a,idx]=sort(A(:));
for n1=1:length(a),
	i2=ceil(idx(n1)/L1);
	i1=idx(n1)-(i2-1)*L1;
	% From B(c1(i1),r1(i1)) to B(c2(i2),r2(i2))
	if c1(i1)==c2(i2), % vertical
		C=c1(i1);
		if r2(i2)<r1(i1), R=(r1(i1):-1:r2(i2))';
		else, R=(r1(i1):r2(i2))'; end
		L=length(R);
		if L<3||~any(B(R(2:end-1),C)), % yup
			C=C(ones(L,1));
			return;
		end
	end
	if r1(i1)==r2(i2), % horizontal
		R=r1(i1);
		if c2(i2)<c1(i1), C=(c1(i1):-1:c2(i2))';
		else, C=(c1(i1):c2(i2))'; end
		L=length(C);
		if L<3||~any(B(R,C(2:end-1))), % yup
			R=R(ones(L,1));
			return;
		end
	end
	if ~B(r2(i2),c1(i1)),
		C1=c1(i1); % vertical portion
		if r2(i2)<r1(i1), R1=(r1(i1):-1:r2(i2))';
		else, R1=(r1(i1):r2(i2))'; end
		L=length(R1);
		if ~any(B(R1(2:end-1),C1)), % maybe
			C1=C1(ones(L,1));
			R2=r2(i2); % horizontal portion
			if c2(i2)<c1(i1), C2=(c1(i1):-1:c2(i2))';
			else, C2=(c1(i1):c2(i2))'; end
			L=length(C2);
			if L<3||~any(B(R2,C2(2:end-1))), % yup
				R2=R2(ones(L,1));
				R=cat(1,R1,R2(2:end));
				C=cat(1,C1,C2(2:end));
				Q=[r2(i2),c1(i1)];
				L=length(R);
				return;
			end
		end
	end
	if ~B(r1(i1),c2(i2)),
		R1=r1(i1); % horizontal portion
		if c2(i2)<c1(i1), C1=(c1(i1):-1:c2(i2))';
		else, C1=(c1(i1):c2(i2))'; end
		L=length(C1);
		if ~any(B(R1,C1(2:end-1))), % maybe
			R1=R1(ones(L,1));
			C2=c2(i2); % vertical portion
			if r2(i2)<r1(i1), R2=(r1(i1):-1:r2(i2))';
			else, R2=(r1(i1):r2(i2))'; end
			L=length(R2);
			if L<3||~any(B(R2(2:end-1),C2)), % yup
				C2=C2(ones(L,1));
				R=cat(1,R1,R2(2:end));
				C=cat(1,C1,C2(2:end));
				Q=[r1(i1),c2(i2)];
				L=length(R);
				return;
			end
		end
	end
end
if ~xtra,
	for n1=1:length(a),
		i2=ceil(idx(n1)/L1);
		i1=idx(n1)-(i2-1)*L1;
		if r1(i1)>1&&~B(r1(i1)-1,c1(i1)),
			[R,C,Q,L]=findpath(B,r1(i1)-1,c1(i1),r2(i2),c2(i2),1);
			if ~isempty(R), 
				R=cat(1,r1(i1),R);
				C=cat(1,c1(i1),C);
				Q=cat(1,[R(1),C(1)],Q);
				L=L+1;
				return;
			end
		end
		if r1(i1)<size(B,1)&&~B(r1(i1)+1,c1(i1)),
			[R,C,Q,L]=findpath(B,r1(i1)+1,c1(i1),r2(i2),c2(i2),1);
			if ~isempty(R), 
				R=cat(1,r1(i1),R);
				C=cat(1,c1(i1),C);
				Q=cat(1,[R(1),C(1)],Q);
				L=L+1;
				return;
			end
		end
		if c1(i1)>1&&~B(r1(i1),c1(i1)-1),
			[R,C,Q,L]=findpath(B,r1(i1),c1(i1)-1,r2(i2),c2(i2),1);
			if ~isempty(R), 
				R=cat(1,r1(i1),R);
				C=cat(1,c1(i1),C);
				Q=cat(1,[R(1),C(1)],Q);
				L=L+1;
				return;
			end
		end
		if c1(i1)<size(B,2)&&~B(r1(i1),c1(i1)+1),
			[R,C,Q,L]=findpath(B,r1(i1),c1(i1)+1,r2(i2),c2(i2),1);
			if ~isempty(R), 
				R=cat(1,r1(i1),R);
				C=cat(1,c1(i1),C);
				Q=cat(1,[R(1),C(1)],Q);
				L=L+1;
				return;
			end
		end
		if r2(i2)>1&&~B(r2(i2)-1,c2(i2)),
			[R,C,Q,L]=findpath(B,r1(i1),c1(i1),r2(i2)-1,c2(i2),1);
			if ~isempty(R), 
				R=cat(1,R,r2(i2));
				C=cat(1,C,c2(i2));
				Q=cat(1,[R(1),C(1)],Q);
				L=L+1;
				return;
			end
		end
		if r2(i2)<size(B,1)&&~B(r2(i2)+1,c2(i2)),
			[R,C,Q,L]=findpath(B,r1(i1),c1(i1),r2(i2)+1,c2(i2),1);
			if ~isempty(R), 
				R=cat(1,R,r2(i2));
				C=cat(1,C,c2(i2));
				Q=cat(1,[R(1),C(1)],Q);
				L=L+1;
				return;
			end
		end
		if c2(i2)>1&&~B(r2(i2),c2(i2)-1),
			[R,C,Q,L]=findpath(B,r1(i1),c1(i1),r2(i2),c2(i2)-1,1);
			if ~isempty(R), 
				R=cat(1,R,r2(i2));
				C=cat(1,C,c2(i2));
				Q=cat(1,[R(1),C(1)],Q);
				L=L+1;
				return;
			end
		end
		if c2(i2)<size(B,2)&&~B(r2(i2),c2(i2)+1),
			[R,C,Q,L]=findpath(B,r1(i1),c1(i1),r2(i2),c2(i2)+1,1);
			if ~isempty(R), 
				R=cat(1,R,r2(i2));
				C=cat(1,C,c2(i2));
				Q=cat(1,[R(1),C(1)],Q);
				L=L+1;
				return;
			end
		end
	end
end
C=[];R=[];L=inf;