001 /*
002 * This file is part of McIDAS-V
003 *
004 * Copyright 2007-2013
005 * Space Science and Engineering Center (SSEC)
006 * University of Wisconsin - Madison
007 * 1225 W. Dayton Street, Madison, WI 53706, USA
008 * https://www.ssec.wisc.edu/mcidas
009 *
010 * All Rights Reserved
011 *
012 * McIDAS-V is built on Unidata's IDV and SSEC's VisAD libraries, and
013 * some McIDAS-V source code is based on IDV and VisAD source code.
014 *
015 * McIDAS-V is free software; you can redistribute it and/or modify
016 * it under the terms of the GNU Lesser Public License as published by
017 * the Free Software Foundation; either version 3 of the License, or
018 * (at your option) any later version.
019 *
020 * McIDAS-V is distributed in the hope that it will be useful,
021 * but WITHOUT ANY WARRANTY; without even the implied warranty of
022 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
023 * GNU Lesser Public License for more details.
024 *
025 * You should have received a copy of the GNU Lesser Public License
026 * along with this program. If not, see http://www.gnu.org/licenses.
027 */
028 package edu.wisc.ssec.mcidasv.util;
029
030 import java.io.*;
031 import java.lang.Math;
032 import java.lang.String;
033 import java.lang.*;
034 import java.rmi.RemoteException;
035
036 import visad.*;
037
038 public class Interpolation {
039
040 static int icnt,jcnt;
041
042 static double ang = Math.PI/2.0;
043
044 static double[][] array = new double[20][20];
045 static double fi;
046 static double fj;
047 static double aa;
048 static double aaa;
049 static double aer;
050 static double ae1 = 0.0;
051 static double ae2 = 0.0;
052
053 static int ijij = 0;
054
055 public static double biquad(double[][] ain, int nx, int ny, double x, double y) {
056
057 double MISSING = 9.9E31;
058 double BADVAL = 1.0E20;
059 double OutQuadInterp;
060
061 double[] dm = new double[4];
062 double[] a = new double[nx*ny];
063 double dst;
064 double dst2;
065 double dst3;
066 double y0 = MISSING;
067 double y1 = MISSING;
068 double y2 = MISSING;
069 double y3 = MISSING;
070 double y12 = MISSING;
071 double y03 = MISSING;
072 double y123 = MISSING;
073
074 int inx;
075 int iny;
076 int nxny;
077 int iz;
078 int imx;
079 int ibad;
080 int ix,iy;
081 int derive03;
082
083 for (ix=0;ix<nx;ix++) {
084 for (iy=0;iy<ny;iy++) {
085 a[(ix*nx)+iy] = ain[ix][iy];
086 }
087 }
088
089 if (x<0.0) x = x + 0.0000001;
090 if (y<0.0) y = y + 0.0000001;
091 if (x>(float)(nx-1)) x = x - 0.0000001;
092 if (y>(float)(ny-1)) y = y - 0.0000001;
093 inx = (int)(Math.floor(x));
094 iny = (int)(Math.floor(y));
095 // nxny = nx*(iny-2) + inx - 1;
096 nxny = nx*(iny-1) + inx -1;
097 dst = (x+1000.0)%1.0;
098 dst2 = dst*dst;
099 dst3 = dst*dst2;
100 iny--;
101 for (iz=0;iz<4;iz++) {
102 y0 = MISSING;
103 y1 = MISSING;
104 y2 = MISSING;
105 y3 = MISSING;
106 // System.out.println("nxny="+nxny+" inx="+inx+" iny="+iny);
107 ibad = 1;
108 // System.out.println("iz="+iz);
109 if ((iny<0)||(iny>=ny)) {
110 // System.out.println("bad");
111 ibad = -1;
112 }
113 else {
114 // System.out.println("here");
115 imx = inx;
116 // System.out.println("imx="+imx);
117 if ((imx>=0)&&(imx<nx)) y1 = a[nxny+1];
118 imx++;
119 // System.out.println("imx="+imx);
120 if ((imx>=0)&&(imx<nx)) y2 = a[nxny+2];
121 // System.out.println("y1="+y1+" y2="+y2);
122 if ((y1>1.0e29)&&(y2>1.0e29)) {
123 ibad = -1;
124 }
125 else {
126 imx = inx - 1;
127 // System.out.println("imx="+imx);
128 if ((imx>=0)&&(imx<nx)&&(a[nxny]<BADVAL)) y0 = a[nxny];
129 imx = imx + 3;
130 // System.out.println("imx="+imx);
131 if ((imx>=0)&&(imx<nx)&&(a[nxny+3]<BADVAL)) y3 = a[nxny+3];
132 // System.out.println("y0="+y0+" y3="+y3);
133 if (y1>BADVAL) {
134 if ((y2>BADVAL)||(y3>BADVAL)) {
135 ibad = -1;
136 }
137 else {
138 y1 = y2 + y2 -y3;
139 }
140 }
141 if ((y2>BADVAL)&&(ibad==1)) {
142 if ((y0>BADVAL)||(y1>BADVAL)) {
143 ibad = -1;
144 }
145 else {
146 y2 = y1 + y1 -y0;
147 }
148 }
149 if (ibad==1) {
150 if (y0>BADVAL) y0 = y1 + y1 - y2;
151 if (y3>BADVAL) y3 = y2 + y2 - y1;
152 y12 = y1 - y2;
153 y03 = y0 - y3;
154 y123 = y12 + y12 + y12;
155 }
156 }
157 }
158 if (ibad==-1) {
159 dm[iz] = MISSING;
160 }
161 else {
162 // System.out.println("dst="+dst+" dst2="+dst2+" dst3="+dst3);
163 // System.out.println("y0="+y0+" y1="+y1+" y2="+y2+" y3="+y3);
164 // System.out.println("y03="+y03+" y12="+y12+" y123="+y123);
165 dm[iz] = (dst3*(y123-y03)) + (dst2*(y0+y03-y1-y123-y12)) + (dst*(y2-y0)) + y1 + y1;
166 }
167 iny++;
168 nxny = nxny+nx;
169 }
170 ibad = 1;
171
172 // System.out.println("dm0= "+dm[0]+" dm1= "+dm[1]+" dm2= "+dm[2]+" dm3= "+dm[3]);
173 if ((dm[1]>BADVAL)&&(dm[2]>BADVAL)) {
174 // bad 1 and 2 values -- cannot interpolate
175 derive03 = -1;
176 }
177 else {
178 if (dm[1]<BADVAL) {
179 // 1 is good, check 2
180 if (dm[2]<BADVAL) {
181 // 2 is good, can derive 0 and 3, if necessary
182 derive03 = 1;
183 }
184 else {
185 // 2 is bad, check 0 to see if 2 can be estimated
186 if (dm[0]<BADVAL) {
187 // 0 is good, estimate 2 and derive 3, if necessary
188 dm[2] = dm[1] + dm[1] - dm[0];
189 derive03 = 1;
190 }
191 else {
192 // 0 is bad, cannot estimate 2 -- cannot interpolate
193 derive03 = -1;
194 }
195 }
196 }
197 else {
198 // 1 is bad, but 2 is good, check 3 to see if 1 can be estimated
199 if (dm[3]<BADVAL) {
200 // 3 is good, estimate 1 and derive 0, if necessary
201 dm[1] = dm[2] + dm[2] - dm[3];
202 derive03 = 1;
203 }
204 else {
205 // 3 is bad, cannot estimate 1 -- cannot interpolate
206 derive03 = -1;
207 }
208 }
209 }
210 if (derive03==1) {
211 // values 1 and 2 are good, will derive 0 and 3, if necessary
212 if(dm[0]>BADVAL) dm[0] = dm[1] + dm[1] - dm[2];
213 if(dm[3]>BADVAL) dm[3] = dm[2] + dm[2] - dm[1];
214 dst = (y+1000.0)%1.0;
215 dst2 = dst*dst;
216 dst3 = dst*dst2;
217 y12 = dm[1] - dm[2];
218 y03 = dm[0] - dm[3];
219 y123 = y12 + y12 + y12;
220 OutQuadInterp = 0.25*(dst3*(y123-y03) + (dst2*(dm[0]+y03-dm[1]-y123-y12)) + (dst*(dm[2]-dm[0])) + dm[1] + dm[1]);
221 }
222 else {
223 // cannot interpolate, return missing value
224 OutQuadInterp = MISSING;
225 }
226
227 return OutQuadInterp;
228 }
229
230 public static FlatField biquad(FlatField fromFld, FlatField toFld) throws VisADException, RemoteException {
231 Gridded2DSet fromSet = (Gridded2DSet) fromFld.getDomainSet();
232 int[] lens = fromSet.getLengths();
233 double[][] fromVals = fromFld.getValues(false);
234 Set toSet = toFld.getDomainSet();
235
236
237 float[][] coords = transform(toSet, fromSet);
238 coords = ((GriddedSet)toSet).valueToGrid(coords);
239
240 double[] newVals = biquad(fromVals, lens[1], lens[0], Set.floatToDouble(coords));
241 FlatField newFld = new FlatField((FunctionType) toFld.getType(), toSet);
242 newFld.setSamples(new double[][] {newVals});
243
244 return newFld;
245 }
246
247 public static double[] biquad(double[][] a, int nx, int ny, double[][] xy) {
248 int numPts = xy[0].length;
249 double[] interpVals = new double[numPts];
250 for (int k=0; k<numPts; k++) {
251 interpVals[k] = biquad(a, nx, ny, xy[0][k], xy[1][k]);
252 }
253 return interpVals;
254 }
255
256 public static double bilinear(double[][] a, int nx, int ny, double x, double y) {
257
258 double MISSING = 9.9E31;
259 double OutBiInterp;
260 double dx;
261 double dy;
262 double vy1;
263 double vy2;
264
265 int nx1;
266 int nx2;
267 int ny1;
268 int ny2;
269
270 if (((x<=0.0)||(x>=(float)(nx)))||((y<=0.0)||(y>=(float)(ny)))) {
271 OutBiInterp = MISSING;
272 }
273 else {
274 // nx1 = Math.max((int)(Math.floor(x-1)),0);
275 // nx2 = nx1 + 1;
276 // ny1 = Math.max((int)(Math.floor(y-1)),0);
277 // ny2 = ny1 + 1;
278 // dx = x - (float)(nx1) - 1; // added -1
279 // dy = y - (float)(ny1) - 1; // added -1
280 nx1 = Math.max((int)(Math.floor(x)),0);
281 nx2 = nx1 + 1;
282 ny1 = Math.max((int)(Math.floor(y)),0);
283 ny2 = ny1 + 1;
284 dx = x - (float)(nx1); // added -1
285 dy = y - (float)(ny1); // added -1
286 // System.out.println("a11="+a[nx1][ny1]+" a21="+a[nx2][ny1]);
287 // System.out.println("a22="+a[nx2][ny2]+" a12="+a[nx2][ny1]);
288 vy1 = a[nx1][ny1] + dx*(a[nx2][ny1] - a[nx1][ny1]);
289 vy2 = a[nx1][ny2] + dx*(a[nx2][ny2] - a[nx1][ny2]);
290 // System.out.println("vy1="+vy1+" vy2="+vy2);
291 OutBiInterp = vy1 + dy*(vy2 - vy1);
292 }
293
294 return OutBiInterp;
295 }
296
297
298 public static float[][] transform(Set fromSet, Set toSet) throws VisADException, RemoteException {
299 CoordinateSystem fromCS = fromSet.getCoordinateSystem();
300 CoordinateSystem toCS = toSet.getCoordinateSystem();
301 Unit[] fromUnits = fromSet.getSetUnits();
302 Unit[] toUnits = toSet.getSetUnits();
303
304 int[] wedge = fromSet.getWedge();
305 float[][] values = fromSet.indexToValue(wedge);
306
307 values = CoordinateSystem.transformCoordinates(
308 ((SetType) fromSet.getType()).getDomain(), fromCS, fromUnits, null,
309 ((SetType) toSet.getType()).getDomain(), toCS, toUnits, null,
310 values);
311
312 return values;
313
314 //values = ((GriddedSet)toSet).valueToGrid(values);
315 }
316
317 public static void main(String[] args) throws IOException {
318
319 for (icnt=0;icnt<20;icnt++) {
320 for (jcnt=0;jcnt<20;jcnt++) {
321 array[icnt][jcnt] = Math.sin(ang*((float)(icnt)/19.0)) + Math.sin(ang*((float)(jcnt)/19.0));
322 // System.out.println("icnt="+icnt+" jcnt="+jcnt+" Array: "+array[icnt][jcnt]);
323 }
324 }
325
326 for (icnt=5;icnt<=189;icnt=icnt+19) {
327 for (jcnt=9;jcnt<=189;jcnt=jcnt+18) {
328 fi = (double)(icnt) * 0.1;
329 fj = (double)(jcnt) * 0.1;
330 aa = (fi * fi) + (fj + fj);
331 aa = Math.sin(ang*((fi)/19.0)) + Math.sin(ang*((fj)/19.0));
332 aaa = bilinear(array,20,20,fi,fj);
333 aer = (aa - aaa)/aa;
334 ae1 = ae1 + Math.abs(aer);
335 System.out.println("Line: "+icnt+" "+jcnt+" "+aa+" "+aaa+" "+aer);
336 aaa = biquad(array,20,20,fi,fj);
337 aer = (aa - aaa)/aa;
338 ae2 = ae2 + Math.abs(aer);
339 System.out.println("Quad: "+icnt+" "+jcnt+" "+aa+" "+aaa+" "+aer);
340 ijij++;
341 }
342 }
343 ae1 = ae1/(float)(ijij);
344 ae2 = ae2/(float)(ijij);
345 System.out.println("ae1,ae2="+ae1+" "+ae2);
346
347 }
348 }