nbmirror=1-nb; nbii=find(nb==1);
nbxmbad(2:N,:)=(nb(2:N,:)-nb(1:N-1,:)==1);
nbxpbad(1:N-1,:)=(nb(2:N,:)-nb(1:N-1,:)==-1);
nbxgood=nb.*(1-nbxmbad).*(1-nbxpbad);
nbymbad(:,2:N)=(nb(:,2:N)-nb(:,1:N-1)==1);
nbypbad(:,1:N-1)=(nb(:,2:N)-nb(:,1:N-1)==-1);
nbygood=nb.*(1-nbymbad).*(1-nbypbad);
nbx1=find(nbxgood==1); nbx2=find(nbxmbad==1); nbx3=find(nbxpbad==1);
nby1=find(nbygood'==1); nby2=find(nbymbad'==1); nby3=find(nbypbad'==1);

eps=2; rr=5;

if fuserf==1, i0=-30; else i0=20; end

for i=i0:niter,
	[level0,level0size]=contour0(lset);
	[closest,closesty,newlset]=exspeedfancy(level0,lset,nb);
	closest=closest+1;
	if mod(i,30)==0,
		%figure(1)
		image(img/100+nb*64); colormap(gray(256));hold on;
		for j=1:level0size,plot(level0(j,2),level0(j,1),'r.');end
        	hold off; %axis([20 80 30 90]);
		axis('square');drawnow;%keyboard
	end
	bb(4,i+31,1:level0size,:)=level0;
	%image(img/max(img(:))*256); colormap(gray(256));hold on;
	%for j=1:level0size,plot(level0(j,2),level0(j,1),'r.');end; hold off;
	%axis('square');drawnow; eval(sprintf('print -djpeg iseq%d',i+31));
	if mod(i,3)==0,lset=newlset+nbmirror.*lset;end
        mapping=nb(:).*closest(:)+nbmirror(:);
	if i<20,
		rf=rforce(level0,level0size,(lset<0)+1,mask)/6000;
	else
		rf=hforce(level0,level0size,(lset<0)+1,rr)/10000;
	end
	if i>niter-20, delt=0.2; eps=3200; end
	if i>niter-15, delt=0.1; rr=3; end
	if i>niter-10, delt=0.05; rr=2; end
	if i>niter-5, delt=0.02; end
        speed(:)=rf(mapping).*nb(:);

	dlx=nbdx(nbx1,nbx2,nbx3,lset);
	dly=nbdx(nby1,nby2,nby3,lset')';
        ndl=sqrt(dlx.^2+dly.^2);
        %weirdos=find(ndl+nbmirror==0); ndl(weirdos)=ndl(weirdos-1);
        dlx1(nbii)=dlx(nbii)./ndl(nbii);
        dly1(nbii)=dly(nbii)./ndl(nbii);
	dlxx=nbdx(nbx1,nbx2,nbx3,dlx1);
	dlyy=nbdx(nby1,nby2,nby3,dly1')';
        dlst=ndl.*((dlxx+dlyy)/eps+speed)*delt;
        lset=lset+dlst;
end
