1 /*
2 * BilinearInterpolator.java
3 *
4 * Created on June 3, 2008, 4:04 PM
5 *
6 * $Id: BilinearInterpolator.java,v 1.1.1.1 2010/11/30 21:32:00 jeremy Exp $
7 */
8 //
9 package org.lcsim.math.interpolation;
10
11 import static java.lang.Math.abs;
12
13 /**
14 * A class to provide interpolated values for values determined at discrete points
15 * on a 2D grid
16 *
17 * @author Norman Graf
18 */
19 public class BilinearInterpolator
20 {
21 private double[] _x;
22 private double _xmin;
23 private double _xmax;
24 private int _xDim;
25 private double[] _y;
26 private double _ymin;
27 private double _ymax;
28 private int _yDim;
29 private double[][] _val;
30 /**
31 * Creates a new instance of BilinearInterpolator
32 * @param x Array of first independent variable at which values are known
33 * @param y Array of second independent variable at which values are known
34 * @param z Array of values at the (x,y) points
35 */
36 public BilinearInterpolator(double[] x, double[] y, double[][] z)
37 {
38 _x = new double[x.length];;
39 System.arraycopy(x,0,_x,0, x.length);
40 _xDim = _x.length;
41 _xmin = _x[0];
42 _xmax = _x[_xDim-1];
43
44 _y = new double[y.length];
45 System.arraycopy(y,0,_y,0, y.length);
46 _yDim = _y.length;
47 _ymin = _y[0];
48 _ymax = _y[_yDim-1];
49
50
51
52 _val = new double[z.length][z[0].length];
53 for(int i=0; i< z.length; ++i)
54 {
55 System.arraycopy(z[i],0,_val[i],0,z[i].length);
56 }
57 }
58
59 //TODO add protection for values out of range
60 /**
61 * Return the value at an arbitrary x,y point, using bilinear interpolation
62 * @param x the first independent variable
63 * @param y the second independent variable
64 * @return the interpolated value at (x,y)
65 */
66 public double interpolateValueAt(double x, double y)
67 {
68 return trueBilinear( x, y );
69 // return polin2(_x, _y, _val, x, y);
70 }
71
72 /*
73 * Given arrays x1a[1..m] and x2a[1..n] of independent variables, and a submatrix of function
74 * values ya[1..m][1..n], tabulated at the grid points defined by x1a and x2a; and given values
75 * x1 and x2 of the independent variables; this routine returns an interpolated function value y.
76 */
77 double polin2(double[] x1a, double[] x2a, double[][] ya, double x1, double x2)
78 {
79 int m = x1a.length;
80 int n = x2a.length;
81 double[] ymtmp = new double[m];
82 double[] yntmp = new double[n];
83 for (int j=0; j<m; ++j) //Loop over rows.
84 {
85 // copy the row into temporary storage
86 for(int k = 0; k<n; ++k)
87 {
88 yntmp[k] = ya[j][k];
89 }
90 ymtmp[j] = polint(x2a, yntmp, x2 ); //Interpolate answer into temporary storage.
91 }
92 return polint(x1a, ymtmp, x1); //Do the final interpolation.
93 }
94
95 /*
96 * Given arrays xa[1..n] and ya[1..n], and given a value x, this routine returns a value y.
97 * If P(x) is the polynomial of degree N -1 such that P(xai) = ya, i ; i = 1... n, then
98 * the returned value y = P(x).
99 */
100 double polint( double[] xa, double[] ya, double x)
101 {
102 int ns=0;
103 int n = xa.length;
104 double den,dif,dift,ho,hp,w;
105 double dy;
106 double[] c = new double[n];
107 double[] d = new double[n];
108 dif=abs(x-xa[0]);
109
110 for (int i=0;i<n;++i)
111 { //Here we find the index ns of the closest table entry,
112 dift = abs(x-xa[i]);
113 if ( dift < dif)
114 {
115 ns=i;
116 dif=dift;
117 }
118 c[i]=ya[i]; // and initialize the table of c's and d's.
119 d[i]=ya[i];
120 }
121 double y = ya[ns--]; // This is the initial approximation to y.
122 // ns = ns-1;
123 for (int m=1; m<n ; ++m)
124 {
125 //For each column of the table,
126 for (int i=0; i<n-m;++i)
127 { //we loop over the current c's and d's and update them.
128 ho=xa[i]-x;
129 hp=xa[i+m]-x;
130 w=c[i+1]-d[i];
131 den = ho-hp;
132 if (den == 0.0) System.err.println("Error in routine polint");
133 //This error can occur only if two input xa's are (to within roundo ) identical.
134 den=w/den;
135 d[i]=hp*den; //Here the c's and d's are updated.
136 c[i]=ho*den;
137 }
138 if(2*ns < (n-m))
139 {
140 dy = c[ns+1];
141 }
142 else
143 {
144 dy = d[ns];
145 ns = ns-1;
146 }
147 y += dy;
148 //After each column in the tableau is completed, we decide which correction, c or d,
149 //we want to add to our accumulating value of y, i.e., which path to take through the
150 //tableau|forking up or down. We do this in such a way as to take the most \straight
151 //line" route through the tableau to its apex, updating ns accordingly to keep track of
152 //where we are. This route keeps the partial approximations centered (insofar as possible)
153 //on the target x. The last dy added is thus the error indication.
154 }
155 return y;
156 }
157
158 double trueBilinear(double x, double y)
159 {
160 // find bin for x
161 int ixlo = 0;
162 int ixhi = 0;
163 for(int i=0; i<_xDim-1; ++i)
164 {
165 if(x>=_x[i] && x<=_x[i+1])
166 {
167 ixlo = i;
168 ixhi = i+1;
169 break;
170 }
171 }
172
173 // find bin for y
174 int iylo = 0;
175 int iyhi = 0;
176 for(int i=0; i<_yDim-1; ++i)
177 {
178 if(y>=_y[i] && y<=_y[i+1])
179 {
180 iylo = i;
181 iyhi = i+1;
182 break;
183 }
184 }
185 double v1 = _val[ixlo][iylo];
186 double v2 = _val[ixhi][iylo];
187 double v3 = _val[ixhi][iyhi];
188 double v4 = _val[ixlo][iyhi];
189
190 double t = (x - _x[ixlo])/(_x[ixhi] - _x[ixlo]);
191 double u = (y - _y[iylo])/(_y[iyhi] - _y[iylo]);
192
193 return (1-t)*(1-u)*v1 +t*(1-u)*v2+t*u*v3+(1-t)*u*v4;
194 }
195 }