util1_jan/coord/bint2_inv.m

   1 function [x1,x2]=bint2_inv(a1,a2,v1,v2)\r
   2 % inverse bilinear interpolation\r
   3 % [x1,x2]=bint2_inv(a1,a2,v1,v2)\r
   4 % find x1 x2 such that v1=bint2(a1,x1,x2) and v2=bint2(a2,x1,x2)\r
   5 % it is assumed that the mapping M given by \r
   6 % M: [i1,i2] -> {a1(i1,i2),a2(i1,i2)] is reasonably close to linear\r
   7 % Jan Mandel, February 2006\r
   8 \r
   9 [m,n]=size(a1);\r
  10 \r
  11 % jacobian D at midpoint\r
  12 d11 = (a1(m,1)-a1(1,1)+a1(m,n)-a1(1,n))/(2*(m-1));\r
  13 d21 = (a2(m,1)-a2(1,1)+a2(m,n)-a2(1,n))/(2*(m-1));\r
  14 d12 = (a1(1,n)-a1(1,1)+a1(m,n)-a1(m,1))/(2*(n-1));\r
  15 d22 = (a2(1,n)-a2(1,1)+a2(m,n)-a2(m,1))/(2*(n-1));\r
  16 detj=d11*d22-d12*d21;\r
  17 \r
  18 % start from midpoint: \r
  19 x1=(1+m)/2;\r
  20 x2=(1+n)/2;\r
  21 \r
  22 maxit=20;\r
  23 tol=1e-10;\r
  24 for it=1:maxit\r
  25     % iteration: D*(xnew-x)=v-bint(x);\r
  26     res1=v1-bint2(a1,x1,x2);\r
  27     res2=v2-bint2(a2,x1,x2);\r
  28     err=max(abs(res1),abs(res2));\r
  29     x1=x1+(d22*res1-d12*res2)/detj;\r
  30     x2=x2+(-d21*res1+d11*res2)/detj;\r
  31     it,x1,x2,err\r
  32     if err<tol,\r
  33         break\r
  34     end\r
  35 end\r