function result=fastint(I,s1,s2)
%interpolates (bilinearly) matrix I given at equally spaced grid points
%to points displaced by vector field s with components s1 and s2
[l1,l2]=size(I);
[Y,X]=meshgrid(1:l2,(1:l1)');
XI=X+s1;M=XI>l1.*ones(l1,l2);XI=(1-M).*XI+M.*l1.*ones(l1,l2);
M=XI<ones(l1,l2);XI=(1-M).*XI+M.*ones(l1,l2);
YI=Y+s2;M=YI>l2.*ones(l1,l2);YI=(1-M).*YI+M.*l2.*ones(l1,l2);
M=YI<ones(l1,l2);YI=(1-M).*YI+M.*ones(l1,l2);


result=interp2(Y,X,I,YI,XI);