set(0, 'DefaultLineMarkerSize', 16)
N=128; lset=zeros(N,N);
[xxm,yym]=meshgrid(1:N,1:N);
img=(sqrt((xxm-60).^2+(yym-40).^2)<20)+(sqrt((xxm-60).^2+(yym-90).^2)<10);
%img=avsread('huber/phtmt.fld')'; img=img/max(img(:)); 
%img(find(img<0.01))=0.18; img(find(img>0.99))=1.3;
img=img+randn(size(xxm))*0.05;
H=fspecial('gaussian',9,1.2);
img2=filter2(H,img);
img2=img2/max(img2(:));
lset=sqrt((xxm-60).^2*1.5+(yym-57).^2/2)-33;
%lset=sqrt((xxm-70).^2*1.5+(yym-60).^2/2)-13;
nb=(xxm>1 & xxm<128 & yym>1 &yym<128);

speed=zeros(N,N); dlx1=zeros(N,N); dly1=zeros(N,N);
speed2=zeros(N,N); speed3=zeros(N,N);
nbxmbad=zeros(N,N); nbxpbad=zeros(N,N); nbymbad=zeros(N,N); nbypbad=zeros(N,N);

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);

Ki=exp(-sqrt(nbdx(nbx1,nbx2,nbx3,img2).^2+...
	(nbdx(nbx1,nbx2,nbx3,img2')').^2).^2*250);
[dx,dy]=sdxy(Ki);

for i=0:120,
        [level0,level0size]=contour0(lset);
        [closest,closesty,newlset]=exspeedfancy(level0,lset,nb);
        closest=closest+1;
        if mod(i,3)==-1,
                figure(1); image((img-0.18)/1.2*256);
		colormap(gray(256));hold on;
                for j=1:level0size,plot(level0(j,2),level0(j,1),'r.');end
                hold off; axis image; axis off; drawnow;
		eval(sprintf('print -djpeg100 -r60 fancy%d',i/3+1));
        end
        if mod(i,3)==0,
		figure(1); contour(lset(2:127,2:127),(-200:5:200));
                hold on;
                for j=1:level0size,plot(level0(j,2)-1,level0(j,1)-1,'r.');end
                hold off; axis image; axis off; drawnow;
		eval(sprintf('print -djpeg100 -r60 mov4/fancy%d',i/3+1));
		figure(2); image((lset(3:126,3:126)+20)*0.74);
		hold on;
                for j=1:level0size,plot(level0(j,2)-1,level0(j,1)-1,'r.');end
		hold off; axis image; axis off; drawnow;
		eval(sprintf('print -djpeg100 -r60 mov5/fancy%d',i/3+1));
        end
        if mod(i,3)==0,lset=newlset+nbmirror.*lset;end
        mapping=nb(:).*closest(:)+nbmirror(:);
	maps=(round(level0(mapping,2))-1)*N+round(level0(mapping,1));
        speed(:)=Ki(maps).*nb(:);
        speed2(:)=dx(maps).*nb(:);
        speed3(:)=dy(maps).*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=(speed.*ndl.*((dlxx+dlyy)+0.4)+(speed2.*dly+speed3.*dlx))*.6;
        lset=lset+dlst;
end

