function [F,G]=solver(chart,a,b,fmax)
smin=inf;
p=perilous(chart,fmax);
for t=0:12
[f,g,s]=solve(chart,a,b,fmax,p,t^2);
if s<smin;smin=s;F=f;G=g;end
end
%----------------------------------------------------------------
function [F,G,s]=solve(chart,a,b,fmax,p,tol)
[m,n,two]=size(chart);
xa=1+mod(a-1,m);ya=1+(a-xa)/m;
xb=1+mod(b-1,m);yb=1+(b-xb)/m;
x=xa;y=ya;u=0;v=0;
db=(xb-xa)^2+(yb-ya)^2;
owb=1;omax=500;o=0;fac=1;irs=0;
F=zeros(omax,1);
G=zeros(omax,1);
while 1
wx=chart(x,y,1);
wy=chart(x,y,2);
if owb
dx=xb-x;dy=yb-y;
else
dx=xa-x;dy=ya-y;
end
q=fac*fmax/(abs(dx)+abs(dy));
uu=fix(q*dx);f=uu-u-wx;
vv=fix(q*dy);g=vv-v-wy;
den=abs(f)+abs(g);
if den>0
q=fmax/den;if q<1;f=fix(q*f);g=fix(q*g);end
end
ox=x;ou=u;u=u+wx+f;x=x+u;
oy=y;ov=v;v=v+wy+g;y=y+v;
if x<1
x=1;u=x-ox;f=u-ou-wx;
elseif x>m
x=m;u=x-ox;f=u-ou-wx;
end
if y<1
y=1;v=y-oy;g=v-ov-wy;
elseif y>n
y=n;v=y-oy;g=v-ov-wy;
end
if p(x,y)
for ii=1:4
switch ii
case 1;xx=x-1;yy=y;
case 2;xx=x+1;yy=y;
case 3;xx=x;yy=y-1;
case 4;xx=x;yy=y+1;
end
uu=xx-ox;ff=uu-ou-wx;
vv=yy-oy;gg=vv-ov-wy;
ok=(xx>=1)&&(xx<=m)&&(yy>=1)&&(yy<=n)&&...
(~p(xx,yy))&&(abs(ff)+abs(gg)<=fmax);
if ok
x=xx;u=uu;f=ff;
y=yy;v=vv;g=gg;
break
end
end
end
if((abs(u)+abs(v)==0)||(abs(f)+abs(g)>fmax))&&(irs<20)&&(o>1)
irs=irs+1;fac=0.90*fac;
o=0;x=xa;y=ya;u=0;v=0;owb=1;db=(xb-xa)^2+(yb-ya)^2;
else
o=o+1;F(o)=f;G(o)=g;
if o>=omax
da=(x-xa)^2+(y-ya)^2;s=da+db+sum(abs(F))+sum(abs(G));
return
end
tmp=(x-xb)^2+(y-yb)^2;db=min(db,tmp);
if owb
if db<=tol;owb=0;end
else
da=(x-xa)^2+(y-ya)^2;
if da<=tol
F=F(1:o);G=G(1:o);s=da+db+sum(abs(F))+sum(abs(G));
return
end
end
end
end
%----------------------------------------------------------------
function p=perilous(chart,fmax)
[m,n,two]=size(chart);
wx=chart(:,:,1);
wy=chart(:,:,2);
[y,x]=meshgrid(1:n,1:m);
x=x+wx;
y=y+wy;
p=((x<1)|(x>m)|(y<1)|(y>n))&(abs(wx)+abs(wy)>=fmax);
|
