001/*
002 * This file is part of McIDAS-V
003 *
004 * Copyright 2007-2016
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 */
028package edu.wisc.ssec.mcidasv.util;
029
030import java.io.*;
031import java.lang.Math;
032import java.lang.String;
033import java.lang.*;
034import java.rmi.RemoteException;
035
036import visad.*;
037
038public 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}