1 /* 2 * $Id: digit2.i,v 1.2 2010-04-18 10:33:38 thiebaut Exp $ 3 * 2D versions of digitize and interp functions (based on mesh_loc) 4 */ 5 /* Copyright (c) 2005, The Regents of the University of California. 6 * All rights reserved. 7 * This file is part of yorick (http://yorick.sourceforge.net). 8 * Read the accompanying LICENSE file for details. 9 */ 10 11 func digit2(y0,x0, y,x,reg, pt=) 12 /* DOCUMENT digit2(y0,x0, y,x) 13 or digit2(y0,x0, y,x,reg) 14 15 return the index of the zone of the point or points (X0,Y0) 16 in the quadrilateral mesh (X,Y) with the optional region 17 array REG. The result has the same dimensions as the input 18 X0 and Y0 arrays. The result is <=0 at points outside the mesh. 19 20 By default, the zone index is an index into an (M-1)-by-(N-1) 21 array, if X and Y are M-by-N. However, if the keyword pt= is 22 non-nil and non-zero, the return value is the index into an 23 M-by-N array in which the first row and column are non-existent 24 (like the optional REG array). 25 26 SEE ALSO: digitize, interp2, mesh_loc, plm 27 */ 28 { 29 zero= array(0, dimsof(x0,y0)); 30 ndx= mesh_loc(y0+zero,x0+zero, y,x,reg); 31 ndx+= zero; /* works around bug- mesh_loc can return scalar */ 32 if (pt) return ndx; 33 m= dimsof(y)(2); 34 return ndx - m - (ndx-1)/m; 35 } 36 37 func interp2(y0,x0, z,y,x,reg, outside=) 38 /* DOCUMENT z0= interp2(y0,x0, z,y,x) 39 or z0= interp2(y0,x0, z,y,x,reg) 40 41 return the bilinear interpolate of the function Z(X,Y) at the 42 points (X0,Y0). The X, Y, and optional REG arrays specify a 43 quadrilateral mesh as for the plm function. The Z values are 44 specified at the vertices of this mesh, so Z must have the 45 same dimensions as X and Y. 46 47 Points outside the mesh get the value 0.0, unless the outside 48 keyword is non-nil, in which case they get that value. 49 50 SEE ALSO: interp, digit2, mesh_loc, plm 51 */ 52 { 53 scalar= !dimsof(x0,y0)(1); 54 if (scalar) { x0= [x0]; y0= [y0]; } 55 ndx= digit2(y0,x0, y,x,reg, pt=1); 56 mask= ndx>0; 57 58 /* first handle points inside mesh */ 59 list= where(mask); 60 if (numberof(list)) { 61 ndx= ndx(list); 62 zero= array(0., dimsof(x0,y0)); 63 x0= (x0+zero)(list); 64 y0= (y0+zero)(list); 65 /* here are corners of the interpolation tets */ 66 m= dimsof(y)(2); 67 x00= [x(ndx-m-1), y(ndx-m-1), z(ndx-m-1)]; 68 x10= [x(ndx-m), y(ndx-m), z(ndx-m)]; 69 x11= [x(ndx), y(ndx), z(ndx)]; 70 x01= [x(ndx-1), y(ndx-1), z(ndx-1)]; 71 /* form the centroid, make centroid the origin */ 72 xc= 0.25*(x00+x10+x01+x11); 73 x0= xc(,1)-x0; 74 y0= xc(,2)-y0; 75 /* form the median vectors */ 76 xm= 0.5*[x11+x10-x01-x00, x11+x01-x10-x00, x11+x00-x10-x01]; 77 xu= xm(,1,1); yu= xm(,2,1); zu= xm(,3,1); 78 xv= xm(,1,2); yv= xm(,2,2); zv= xm(,3,2); 79 xw= xm(,1,3); yw= xm(,2,3); zw= xm(,3,3); 80 x00= x10= x11= x01= xm= []; 81 /* form various cross products */ 82 cwu= xw*yu - yw*xu; 83 cwv= xw*yv - yw*xv; 84 cuv= xu*yv - yu*xv; 85 cup= xu*y0 - yu*x0; 86 cvp= xv*y0 - yv*x0; 87 cwp= xw*y0 - yw*x0; 88 /* compute the discriminant */ 89 cuv*= 0.5; 90 cwu*= 2.0; 91 cwv*= 2.0; 92 bu= cwp - cuv; 93 bv= cwp + cuv; 94 tmpa= bu*bu; 95 tmpb= bv*bv; 96 use= double(tmpa<=tmpb); 97 dis= sqrt(max(use*(tmpa-cwu*cvp) + (1.-use)*(tmpb-cwv*cup), 0.0)); 98 /* only one solution is the correct one, but there is no way to 99 know which it is without computing both (want to allow bowtied 100 zones and other pathologies) */ 101 mask2= bu>=0.0; 102 list= where(mask2); 103 if (numberof(list)) { 104 tmpa= -bu(list)-dis(list); 105 tmpb= cwu(list); 106 epsa= double(!tmpa)*1.e-99; 107 epsb= double(!tmpb)*1.e-99; 108 u11= tmpa/(tmpb+epsb); 109 u21= cvp(list)/(tmpa+epsa); 110 } 111 list= where(!mask2); 112 if (numberof(list)) { 113 tmpa= -bu(list)+dis(list); 114 tmpb= cwu(list); 115 epsa= double(!tmpa)*1.e-99; 116 epsb= double(!tmpb)*1.e-99; 117 u22= tmpa/(tmpb+epsb); 118 u12= cvp(list)/(tmpa+epsa); 119 } 120 u1= merge(u11,u12,mask2); 121 u2= merge(u21,u22,mask2); 122 u11= u21= u12= u22= []; 123 mask2= bv>=0.0; 124 list= where(mask2); 125 if (numberof(list)) { 126 tmpa= -bv(list)-dis(list); 127 tmpb= cwv(list); 128 epsa= double(!tmpa)*1.e-99; 129 epsb= double(!tmpb)*1.e-99; 130 v21= tmpa/(tmpb+epsb); 131 v11= cup(list)/(tmpa+epsa); 132 } 133 list= where(!mask2); 134 if (numberof(list)) { 135 tmpa= -bv(list)+dis(list); 136 tmpb= cwv(list); 137 epsa= double(!tmpa)*1.e-99; 138 epsb= double(!tmpb)*1.e-99; 139 v12= tmpa/(tmpb+epsb); 140 v22= cup(list)/(tmpa+epsa); 141 } 142 v1= merge(v11,v12,mask2); 143 v2= merge(v21,v22,mask2); 144 v11= v21= v12= v22= []; 145 /* compute the two z values and select the proper one */ 146 xc= xc(,3); 147 z1= zw*2.*u1*v1 + zv*v1 + zu*u1 + xc; 148 z2= zw*2.*u2*v2 + zv*v2 + zu*u2 + xc; 149 mask2= max(abs(v1),abs(u1))<=max(abs(v2),abs(u2)); 150 za= merge(z1(where(mask2)),z2(where(!mask2)),mask2); 151 } 152 153 /* just punt on points outside mesh */ 154 list= where(!mask); 155 if (numberof(list)) { 156 if (is_void(outside)) outside= 0.0; 157 zb= array(outside, numberof(list)); 158 } 159 160 if (scalar) mask= mask(1); 161 return merge(za,zb,mask); 162 } 163