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