Source code

001/*
002 * This file is part of McIDAS-V
003 *
004 * Copyright 2007-2017
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.data.hydra;
029
030import visad.Data;
031import visad.FlatField;
032import visad.VisADException;
033import visad.FunctionType;
034import visad.RealType;
035import visad.RealTupleType;
036import visad.Linear1DSet;
037import visad.Linear2DSet;
038import visad.Gridded2DSet;
039import visad.CoordinateSystem;
040import visad.CommonUnit;
041import visad.SetType;
042import visad.georef.MapProjection;
043import visad.util.ThreadManager;
044
045import java.util.Arrays;
046import java.rmi.RemoteException;
047
048public class ReprojectSwath {
049  private static int count = 0;
050
051  Linear2DSet grid;
052  Linear2DSet swathDomain;
053  FunctionType ftype;
054  float[][] swathRange;
055  CoordinateSystem swathCoordSys;
056  CoordinateSystem gridCoordSys;
057
058  float[][] allSwathGridCoords;
059  int[] allSwathGridIndexs;
060  float[][] swathGridCoord;
061  int[] swathIndexAtGrid;
062
063  int trackLen;
064  int xtrackLen;
065
066  int gridXLen;
067  int gridYLen;
068  int gridLen;
069
070  float[][] gridRange;
071  int rngTupDim;
072  FlatField grdFF;
073
074  int[][][] quads;
075  int mode;
076  
077  public static final int NEAREST = 1;
078  public static final int BILINEAR_VISAD = 0;
079  
080  int numProc = Runtime.getRuntime().availableProcessors();
081  private static boolean doParallel = false;
082
083  private static ReprojectSwath lastReproject = null;
084
085  public static void setDoParallel(boolean enable) {
086     doParallel = enable;
087  }
088  
089  public static FlatField swathToGrid(Linear2DSet grid, FlatField[] swaths, int mode) throws Exception {
090     return swathToGrid(grid, swaths, mode, true);
091  }
092  
093  public static FlatField swathToGrid(Linear2DSet grid, FlatField[] swaths, int mode, boolean filter) throws Exception {
094    FunctionType ftype = (FunctionType) swaths[0].getType();
095    visad.Set domSet = swaths[0].getDomainSet();
096
097    FlatField swath = new FlatField(new FunctionType(ftype.getDomain(),
098        new RealTupleType(new RealType[] 
099           {RealType.getRealType("redimage_"+count), RealType.getRealType("greenimage_"+count), RealType.getRealType("blueimage_"+count)})), domSet);
100
101    swath.setSamples(new float[][]
102        {swaths[0].getFloats(false)[0], swaths[1].getFloats(false)[0], swaths[2].getFloats(false)[0]});
103
104    count++;
105
106    return swathToGrid(grid, swath, mode);
107  }
108  
109  public static FlatField swathToGrid(Linear2DSet grid, FlatField swath, int mode) throws Exception {
110      return swathToGrid(grid, swath, mode, true);
111  }
112  
113  public static FlatField swathToGrid(Linear2DSet grid, FlatField swath, int mode, boolean filter) throws Exception {
114    ReprojectSwath obj = null;
115    FlatField ff = null;
116
117    if (lastReproject != null) {
118       if (grid.equals(lastReproject.grid) && (swath.getDomainSet()).equals(lastReproject.swathDomain)) {
119         obj = lastReproject;
120         ff = obj.reproject(swath, mode, filter);
121       }
122       else {
123         obj = new ReprojectSwath(grid, swath);
124         ff = obj.reproject(mode, filter);
125       }
126    }
127    else {
128      obj = new ReprojectSwath(grid, swath);
129      ff = obj.reproject(mode, filter);
130    }
131    lastReproject = obj;
132
133    return ff;
134  }
135  
136  public ReprojectSwath() {
137  }
138  
139  public ReprojectSwath(Linear2DSet grid, FlatField swath) throws Exception {
140      
141    init(grid, swath);
142    
143    if (trackLen < 200 || doParallel == false) {
144        numProc = 1;
145    }
146
147    projectSwathToGrid();
148  }
149  
150  private void init(Linear2DSet grid, FlatField swath) throws VisADException {
151     this.grid = grid;
152     gridLen = grid.getLength();
153     int[] lens = grid.getLengths();
154     gridXLen = lens[0];
155     gridYLen = lens[1];
156     gridCoordSys = grid.getCoordinateSystem();
157     
158     swathDomain = (Linear2DSet) swath.getDomainSet();
159     lens = swathDomain.getLengths();
160     trackLen = lens[1];
161     xtrackLen = lens[0];
162     int swathLen = trackLen*xtrackLen;
163     swathCoordSys = swathDomain.getCoordinateSystem();
164     swathRange = swath.getFloats(false);
165     ftype = (FunctionType) swath.getType();
166    
167     allSwathGridCoords = new float[2][swathLen];
168     allSwathGridIndexs = new int[swathLen];
169     swathGridCoord = new float[2][gridLen];
170     swathIndexAtGrid = new int[gridLen];
171    
172     quads = new int[gridYLen][gridXLen][4];
173   }
174
175  /*
176  public FlatField reproject(Linear2DSet grid, FlatField swath, int mode, boolean filter) throws Exception {
177    init();
178    
179    if (trackLen < 200 || doParallel == false) {
180        numProc = 1;
181    }
182  
183    return reproject(mode, filter);
184  }
185  */
186    
187  public FlatField reproject(int mode, boolean filter) throws Exception {
188    this.mode = mode;
189    
190    initGrid();
191
192    getBoundingQuadAtGridPts();
193
194    interpolateToGrid();
195    
196    if (filter) {
197       grdFF.setSamples(filter(), false);
198    }
199    else {
200       grdFF.setSamples(gridRange, false);
201    }
202    
203    return grdFF;
204  }
205
206  private FlatField reproject(FlatField swath, int mode, boolean filter) throws Exception {
207     this.mode = mode;
208     swathRange = swath.getFloats(false);
209     ftype = (FunctionType) swath.getType();
210     
211     initGrid();
212     
213     getBoundingQuadAtGridPts();
214     
215     interpolateToGrid();
216     
217     if (filter) {
218       grdFF.setSamples(filter(), false);
219     }
220     else {
221       grdFF.setSamples(gridRange, false);
222     }
223     
224     return grdFF;
225  }
226  
227   private void getBoundingQuadAtGridPts() throws VisADException, RemoteException {
228    int ystart = 3;
229    int ystop = gridYLen-4;
230    int subLen = ((ystop - ystart)+1)/numProc;
231    int rem = ((ystop - ystart)+1) % numProc;
232    
233    
234    ThreadManager threadManager = new ThreadManager("getBoundingQuadAtGridPts");
235    for (int i=0; i<numProc; i++) {
236        final int start = i*subLen + ystart;
237        final int stop = (i != numProc-1 ) ? (start + subLen - 1): (start + subLen + rem - 1);
238          threadManager.addRunnable(new ThreadManager.MyRunnable() {
239                  public void run()  throws Exception {
240                     getBoundingQuadAtGridPts(start, stop);
241                  }
242              });
243    }
244    
245    if (numProc == 1 || !doParallel) {
246       threadManager.runSequentially();
247    }
248    else {
249       threadManager.runAllParallel();
250    }
251  }
252
253  // start to stop inclusive
254  private void getBoundingQuadAtGridPts(int grdYstart, int grdYstop) {
255
256    for (int j=grdYstart; j<=grdYstop; j++) {
257       for (int i=3; i<gridXLen-3; i++) {
258          int grdIdx = i + j*gridXLen;
259
260          int ll = findSwathGridLoc(grdIdx, swathGridCoord, gridYLen, gridXLen, "LL");
261          quads[j][i][0] = ll;
262
263          int lr = findSwathGridLoc(grdIdx, swathGridCoord, gridYLen, gridXLen, "LR");
264          quads[j][i][1] = lr;
265
266          int ul = findSwathGridLoc(grdIdx, swathGridCoord, gridYLen, gridXLen, "UL");
267          quads[j][i][2] = ul;
268
269          int ur = findSwathGridLoc(grdIdx, swathGridCoord, gridYLen, gridXLen, "UR");
270          quads[j][i][3] = ur;
271       }
272    }
273  }
274  
275  public void interpolateToGrid() throws VisADException, RemoteException {
276    int ystart = 3;
277    int ystop = gridYLen-4;
278    int subLen = ((ystop - ystart)+1)/numProc;
279    int rem = ((ystop - ystart)+1) % numProc;
280    
281    ThreadManager threadManager = new ThreadManager("interpolateToGrid");
282    for (int i=0; i<numProc; i++) {
283        final int start = i*subLen + ystart;
284        final int stop = (i != numProc-1 ) ? (start + subLen - 1): (start + subLen + rem - 1);
285          threadManager.addRunnable(new ThreadManager.MyRunnable() {
286                  public void run()  throws Exception {
287                     interpolateToGrid(start, stop);
288                  }
289              });
290    }
291    
292    if (numProc == 1 || !doParallel) {
293       threadManager.runSequentially();
294    }
295    else {
296       threadManager.runAllParallel();
297    }
298  }
299
300  // start to stop inclusive
301  public void interpolateToGrid(int grdYstart, int grdYstop) throws VisADException, RemoteException {
302
303    float[][] corners = new float[2][4];
304    float[][] rngVals = new float[rngTupDim][4];
305    float[] values = new float[4];
306    float gx;
307    float gy;
308
309    for (int j=grdYstart; j<=grdYstop; j++) {
310       for (int i=3; i<gridXLen-3; i++) {
311          int grdIdx = i + j*gridXLen;
312          gx = (float) (grdIdx % gridXLen);
313          gy = (float) (grdIdx / gridXLen);
314
315          java.util.Arrays.fill(corners[0], Float.NaN);
316          java.util.Arrays.fill(corners[1], Float.NaN);
317        
318          int ll = quads[j][i][0];
319          int lr = quads[j][i][1];
320          int ul = quads[j][i][2];
321          int ur = quads[j][i][3];
322
323          if (ll >= 0) {
324             corners[0][0] = swathGridCoord[0][ll] - gx;
325             corners[1][0] = swathGridCoord[1][ll] - gy;
326          }
327          if (lr >= 0) {
328             corners[0][1] = swathGridCoord[0][lr] - gx;
329             corners[1][1] = swathGridCoord[1][lr] - gy;
330          }
331          if (ul >= 0) {
332             corners[0][2] = swathGridCoord[0][ul] - gx;
333             corners[1][2] = swathGridCoord[1][ul] - gy;
334          }
335          if (ur >= 0) {
336             corners[0][3] = swathGridCoord[0][ur] - gx;
337             corners[1][3] = swathGridCoord[1][ur] - gy;
338          }
339
340          if (mode == NEAREST) { // Nearest neighbor
341             for (int t=0; t<rngTupDim; t++) {
342                java.util.Arrays.fill(values, Float.NaN);
343                if (ll >=0) {
344                   values[0] = swathRange[t][swathIndexAtGrid[ll]];
345                }
346                if (lr >= 0) {
347                   values[1] = swathRange[t][swathIndexAtGrid[lr]];
348                }
349                if (ul >= 0) {
350                   values[2] = swathRange[t][swathIndexAtGrid[ul]];
351                }
352                if (ur >= 0) {
353                   values[3] = swathRange[t][swathIndexAtGrid[ur]];
354                }
355                float val = nearest(0f, 0f, corners, values);
356                gridRange[t][grdIdx] = val;
357             }
358          }
359          else if (mode == BILINEAR_VISAD) {  //from VisAD
360             if (!(ll >= 0 && lr >= 0 && ul >= 0 && ur >= 0)) {
361                continue;
362             }
363             corners[0][0] = swathGridCoord[0][ll];
364             corners[1][0] = swathGridCoord[1][ll];
365             corners[0][1] = swathGridCoord[0][lr];
366             corners[1][1] = swathGridCoord[1][lr];
367             corners[0][2] = swathGridCoord[0][ul];
368             corners[1][2] = swathGridCoord[1][ul];
369             corners[0][3] = swathGridCoord[0][ur];
370             corners[1][3] = swathGridCoord[1][ur];
371             for (int t=0; t<rngTupDim; t++) {
372                java.util.Arrays.fill(values, Float.NaN);
373                values[0] = swathRange[t][swathIndexAtGrid[ll]];
374                values[1] = swathRange[t][swathIndexAtGrid[lr]];
375                values[2] = swathRange[t][swathIndexAtGrid[ul]];
376                values[3] = swathRange[t][swathIndexAtGrid[ur]];
377                float val = visad2D(gy, gx, corners, values);
378                gridRange[t][grdIdx] = val;
379             }
380          }
381          else if (mode == 1) { //TODO: not working yet
382             if (!(ll >= 0 && lr >= 0 && ul >= 0 && ur >= 0)) {
383                continue;
384             }
385             gx -= swathGridCoord[0][ll];
386             gy -= swathGridCoord[1][ll];
387             corners[0][0] = 0f;
388             corners[1][0] = 0f;
389             corners[0][1] = swathGridCoord[0][ul] - swathGridCoord[0][ll];
390             corners[1][1] = swathGridCoord[1][ul] - swathGridCoord[1][ll];
391             corners[0][2] = swathGridCoord[0][ur] - swathGridCoord[0][ll];
392             corners[1][2] = swathGridCoord[1][ur] - swathGridCoord[1][ll];
393             corners[0][3] = swathGridCoord[0][lr] - swathGridCoord[0][ll];
394             corners[1][3] = swathGridCoord[1][lr] - swathGridCoord[1][ll];
395             for (int t=0; t<rngTupDim; t++) {
396                java.util.Arrays.fill(values, Float.NaN);
397                values[0] = swathRange[t][swathIndexAtGrid[ll]];
398                values[1] = swathRange[t][swathIndexAtGrid[ul]];
399                values[2] = swathRange[t][swathIndexAtGrid[ur]];
400                values[3] = swathRange[t][swathIndexAtGrid[lr]];
401                float val = biLinearIntrp(gy, gx, corners, values);
402                gridRange[t][grdIdx] = val;
403             }
404
405          }
406       }
407    }
408
409  }
410
411 public void projectSwathToGrid() throws VisADException, RemoteException {
412    int subLen = trackLen/numProc;
413    int rem = trackLen % numProc;
414    
415    ThreadManager threadManager = new ThreadManager("projectSwathToGrid");
416    for (int i=0; i<numProc; i++) {
417        final int start = i*subLen;
418        final int stop = (i != numProc-1 ) ? (start + subLen - 1): (start + subLen + rem - 1);
419          threadManager.addRunnable(new ThreadManager.MyRunnable() {
420                  public void run()  throws Exception {
421                     projectSwathToGrid(start, stop);
422                  }
423              });
424    }
425    
426    if (numProc == 1 || !doParallel) {
427       threadManager.runSequentially();
428    }
429    else {
430       threadManager.runAllParallel();
431    }
432 }
433 
434 public void projectSwathToGrid(int trackStart, int trackStop) throws VisADException, RemoteException {
435
436    for (int j=trackStart; j <= trackStop; j++) {
437       for (int i=0; i < xtrackLen; i++) {
438         int swathIdx = j*xtrackLen + i;
439
440         float[][] swathCoord = swathDomain.indexToValue(new int[] {swathIdx});
441         float[][] swathEarthCoord = swathCoordSys.toReference(swathCoord);
442
443         float[][] gridValue = gridCoordSys.fromReference(swathEarthCoord);
444         float[][] gridCoord = grid.valueToGrid(gridValue);
445         float g0 = gridCoord[0][0];
446         float g1 = gridCoord[1][0];
447         int grdIdx = (g0 != g0 || g1 != g1) ? -1 : ((int) g0) + gridXLen * ((int) g1);
448         int m=0;
449         int n=0;
450         int k = grdIdx + (m + n*gridXLen);
451         
452         allSwathGridCoords[0][swathIdx] = g0;
453         allSwathGridCoords[1][swathIdx] = g1;
454         allSwathGridIndexs[swathIdx] = k;
455             
456       }
457    }
458 }
459 
460 public void initGrid() throws VisADException {
461    Arrays.fill(swathGridCoord[0], -999.9f);
462    Arrays.fill(swathGridCoord[1], -999.9f);
463    Arrays.fill(swathIndexAtGrid, -1);
464
465    for (int j=0; j < trackLen; j++) {
466       for (int i=0; i < xtrackLen; i++) {
467         int swathIdx = j*xtrackLen + i;
468         float val = swathRange[0][swathIdx];
469         int k = allSwathGridIndexs[swathIdx];
470         
471         if ( !(Float.isNaN(val)) && ((k >=0) && (k < gridLen)) ) { // val or val[rngTupDim] ?
472            if (swathIndexAtGrid[k] == -1) {
473               swathGridCoord[0][k] = allSwathGridCoords[0][swathIdx];
474               swathGridCoord[1][k] = allSwathGridCoords[1][swathIdx];
475               swathIndexAtGrid[k] = swathIdx;
476            }
477         }
478       }
479    }
480    
481    for (int j=0; j<gridYLen; j++) {
482       for (int i=0; i<gridXLen; i++) {
483          java.util.Arrays.fill(quads[j][i], -1);
484       }
485    }
486    
487    RealTupleType rtt = ((SetType)grid.getType()).getDomain();
488    grdFF = new FlatField(new FunctionType(rtt, ftype.getRange()), grid);
489    gridRange = grdFF.getFloats(false);
490    rngTupDim = gridRange.length;
491    for (int t=0; t<rngTupDim; t++) {
492       java.util.Arrays.fill(gridRange[t], Float.NaN);
493    }
494 }
495
496 private float[][] filter() throws VisADException, RemoteException {
497
498    double mag = 3.0;
499    double sigma = 0.4;
500
501    float[][] weights = new float[3][3];
502
503    float sumWeights = 0f;
504    for (int n=-1; n<=1; n++) {
505       for (int m=-1; m<=1; m++) {
506          double del_0 = m;
507          double del_1 = n;
508          double dst = Math.sqrt(del_0*del_0 + del_1*del_1);
509
510          weights[n+1][m+1] = (float) (mag/Math.exp(dst/sigma));
511
512          sumWeights += weights[n+1][m+1];
513       }
514    }
515
516    for (int n=-1; n<=1; n++) {
517       for (int m=-1; m<=1; m++) {
518          weights[n+1][m+1] /= sumWeights;
519        }
520    }
521
522    float[][] newRange = new float[rngTupDim][gridLen];
523    for (int t=0; t<rngTupDim; t++) {
524       java.util.Arrays.fill(newRange[t], Float.NaN);
525    }
526    float[] sum = new float[rngTupDim];
527
528    for (int j=2; j<gridYLen-2; j++) {
529       for (int i=2; i<gridXLen-2; i++) {
530         int grdIdx = i + j*gridXLen;
531
532         java.util.Arrays.fill(sum, 0f);
533         for (int n=-1; n<=1; n++) {
534            for (int m=-1; m<=1; m++) {
535               int k = grdIdx + (m + n*gridXLen);
536
537               for (int t=0; t<rngTupDim; t++) {
538                  sum[t] += weights[n+1][m+1]*gridRange[t][k];
539               }
540            }
541         }
542
543         for (int t=0; t<rngTupDim; t++) {
544            newRange[t][grdIdx] = sum[t];
545         }
546       }
547    }
548
549    return newRange;
550 }
551 private static int findSwathGridLoc(int grdIdx, float[][] swathGridCoord, int gridYLen, int gridXLen, String which) {
552  
553    int idx = -1;
554
555    int gy = grdIdx/gridXLen;
556    int gx = grdIdx % gridXLen;
557    
558    int yu1 = (gy+1)*gridXLen;
559    int yu2 = (gy+2)*gridXLen;
560    int yu3 = (gy+3)*gridXLen;    
561    
562    int yd1 = (gy-1)*gridXLen;
563    int yd2 = (gy-2)*gridXLen;
564    int yd3 = (gy-3)*gridXLen;
565
566    switch (which) {
567       case "LL":
568
569          idx = yd1 + (gx-1);
570          if (swathGridCoord[0][idx] != -999.9f) {
571             break;
572          }
573          idx = yd2 + (gx-1);
574          if (swathGridCoord[0][idx] != -999.9f) {
575             break;
576          }
577          idx = yd1 + (gx-2);
578          if (swathGridCoord[0][idx] != -999.9f) {
579             break;
580          }
581          idx = yd2 + (gx-2);
582          if (swathGridCoord[0][idx] != -999.9f) {
583             break;
584          }
585
586
587          idx = yd2 + (gx-3);
588          if (swathGridCoord[0][idx] != -999.9f) {
589             break;
590          }
591          idx = yd3 + (gx-2);
592          if (swathGridCoord[0][idx] != -999.9f) {
593             break;
594          }
595          idx = yd3 + (gx-1);
596          if (swathGridCoord[0][idx] != -999.9f) {
597             break;
598          }
599          idx = yd1 + (gx-3);
600          if (swathGridCoord[0][idx] != -999.9f) {
601             break;
602          }
603          idx = yd3 + (gx-3);
604          if (swathGridCoord[0][idx] != -999.9f) {
605             break;
606          }
607
608          idx = -1;
609          break;
610       case "UL":
611          idx = (gy)*gridXLen + (gx-1);
612          if (swathGridCoord[0][idx] != -999.9f) {
613             break;
614          }
615          idx = (gy)*gridXLen + (gx-2);
616          if (swathGridCoord[0][idx] != -999.9f) {
617             break;
618          }
619          idx = yu1 + (gx-1);
620          if (swathGridCoord[0][idx] != -999.9f) {
621             break;
622          }
623          idx = yu1 + (gx-2);
624          if (swathGridCoord[0][idx] != -999.9f) {
625             break;
626          }
627
628          idx = yu1 + (gx-3);
629          if (swathGridCoord[0][idx] != -999.9f) {
630             break;
631          }
632          idx = yu2 + (gx-3);
633          if (swathGridCoord[0][idx] != -999.9f) {
634             break;
635          }
636          idx = yu2 + (gx-2);
637          if (swathGridCoord[0][idx] != -999.9f) {
638             break;
639          }
640          idx = (gy)*gridXLen + (gx-3);
641          if (swathGridCoord[0][idx] != -999.9f) {
642             break;
643          }
644          idx = yd1 + (gx-3);
645          if (swathGridCoord[0][idx] != -999.9f) {
646             break;
647          }
648
649          idx = -1;
650          break;
651       case "UR":
652          idx = (gy)*gridXLen + (gx);
653          if (swathGridCoord[0][idx] != -999.9f) {
654             break;
655          }
656          idx = yu1 + (gx);
657          if (swathGridCoord[0][idx] != -999.9f) {
658             break;
659          }
660          idx = (gy)*gridXLen + (gx+1);
661          if (swathGridCoord[0][idx] != -999.9f) {
662             break;
663          }
664          idx = yu1 + (gx+1);
665          if (swathGridCoord[0][idx] != -999.9f) {
666             break;
667          }
668
669          idx = yu2 + (gx+1);
670          if (swathGridCoord[0][idx] != -999.9f) {
671             break;
672          }
673          idx = yu1 + (gx+2);
674          if (swathGridCoord[0][idx] != -999.9f) {
675             break;
676          }
677          idx = yu2 + (gx+2);
678          if (swathGridCoord[0][idx] != -999.9f) {
679             break;
680          }
681          idx = (gy)*gridXLen + (gx+2);
682          if (swathGridCoord[0][idx] != -999.9f) {
683             break;
684          }
685          idx = yd1 + (gx+2);
686          if (swathGridCoord[0][idx] != -999.9f) {
687             break;
688          }
689
690          idx = -1;
691          break;
692       case "LR":
693          idx = yd1 + (gx);
694          if (swathGridCoord[0][idx] != -999.9f) {
695             break;
696          }
697          idx = yd1 + (gx+1);
698          if (swathGridCoord[0][idx] != -999.9f) {
699             break;
700          }
701          idx = yd2 + (gx);
702          if (swathGridCoord[0][idx] != -999.9f) {
703             break;
704          }
705          idx = yd2 + (gx+1);
706          if (swathGridCoord[0][idx] != -999.9f) {
707             break;
708          }
709
710          idx = yd1 + (gx+2);
711          if (swathGridCoord[0][idx] != -999.9f) {
712             break;
713          }
714          idx = yd2 + (gx+2);
715          if (swathGridCoord[0][idx] != -999.9f) {
716             break;
717          }
718          idx = yd3 + (gx);
719          if (swathGridCoord[0][idx] != -999.9f) {
720             break;
721          }
722          idx = yd3 + (gx+1);
723          if (swathGridCoord[0][idx] != -999.9f) {
724             break;
725          }
726          idx = yd3 + (gx+2);
727          if (swathGridCoord[0][idx] != -999.9f) {
728             break;
729          }
730 
731          idx = -1;
732          break;
733    }
734
735    return idx;
736 }
737
738 /* Reference: David W. Zingg, University of Toronto, Downsview, Ontario, Canada
739               Maurice Yarrow, Sterling Software, Arnes Research Center, Moffett Field, California. 
740               NASA Technical Memorandum 102213
741
742     y -> q, x -> p; first point (x=0, y=0) and clockwise around
743  */
744 public static float biLinearIntrp(float gy, float gx, float[][] corners, float[] values) {
745    // transform physical space (gy, gx) to unit square (q, p)
746    // bilinear mapping coefficients
747/*
748    float a = corners[0][3];
749    float b = corners[0][1];
750    float c = (corners[0][2] - corners[0][3] - corners[0][1]);
751
752    float d = corners[1][3];
753    float e = corners[1][1];
754    float f = (corners[1][2] - corners[1][3] - corners[1][1]);
755*/
756
757    float a = corners[0][1];
758    float b = corners[0][3];
759    float c = (corners[0][2] - corners[0][1] - corners[0][3]);
760
761    float d = corners[1][1];
762    float e = corners[1][3];
763    float f = (corners[1][2] - corners[1][1] - corners[1][3]);
764
765
766    // quadratic terms to determine p
767    // p = (-coef1 +/- sqrt(coef1^2 - 4*coef2*coef0))/2*coef2
768
769    float coef2 = (c*d - a*f);  // A 
770    float coef1 = (-c*gy + b*d + gx*f - a*e);  // B
771    float coef0 = (-gy*b + gx*e);  // C
772
773    // two p vals from quadratic:
774    float p0 =  (-coef1 + ((float)Math.sqrt(coef1*coef1 - 4*coef2*coef0)) )/2f*coef2;
775    float p1 =  (-coef1 - ((float)Math.sqrt(coef1*coef1 - 4*coef2*coef0)) )/2f*coef2;
776
777    // the corresponding q values for both p solutions:
778    float q0 = (gx - a*p0)/(b + c*p0);
779    float q1 = (gx - a*p1)/(b + c*p1);
780
781    // test  which point to use. One will be outside unit square:
782    float p = Float.NaN;
783    float q = Float.NaN;
784    if ((p0 >= 0f && p0 <= 1f) && (q0 >= 0f && q0 <= 1f)) {
785       p = p0;
786       q = q0;
787    }
788    else if ((p1 >= 0f && p1 <= 1f) && (q1 >= 0f && q1 <= 1f)) {
789       p = p1;
790       q = q1;
791    }
792
793    // bilinear interpolation within the unit square:
794
795    float intrp = values[0]*(1f-p)*(1f-q) + values[1]*(1f-p)*q + values[2]*p*q + values[3]*p*(1f-q);
796
797    return intrp;
798 }
799
800 private static float nearest(float gy, float gx, float[][] corners, float[] values) {
801   float minDist = Float.MAX_VALUE;
802
803   float delx;
804   float dely;
805   int closest = 0;
806   for (int k=0; k<4; k++) {
807      delx = corners[0][k] - gx;
808      dely = corners[1][k] - gy;
809      float dst_sqrd = delx*delx + dely*dely;
810
811      if (dst_sqrd < minDist) {
812         minDist = dst_sqrd;
813         closest = k;
814      }
815   }
816
817   return values[closest];
818 }
819 public static float visad2D(float gy, float gx, float[][] corners, float[] values) {
820    boolean Pos = true;
821
822    // A:0, B:1, C:2, D:3
823    float v0x = corners[0][0];
824    float v0y = corners[1][0];
825    float v1x = corners[0][1];
826    float v1y = corners[1][1];
827    float v2x = corners[0][2];
828    float v2y = corners[1][2];
829    float v3x = corners[0][3];
830    float v3y = corners[1][3];
831
832    float bdx = v2x-v1x;
833    float bdy = v2y-v1y;
834
835    float bpx = gx-v1x;
836    float bpy = gy-v1y;
837
838    float dpx = gx-v2x;
839    float dpy = gy-v2y;
840
841    // lower triangle first
842
843    boolean insideTri = true;
844
845    float abx = v1x-v0x;
846    float aby = v1y-v0y;
847    float dax = v0x-v2x;
848    float day = v0y-v2y;
849    float apx = gx-v0x;
850    float apy = gy-v0y;
851    
852    float tval1 = abx*apy-aby*apx;
853    float tval2 = bdx*bpy-bdy*bpx;
854    float tval3 = dax*dpy-day*dpx;
855    boolean test1 = (tval1 == 0) || ((tval1 > 0) == Pos);
856    boolean test2 = (tval2 == 0) || ((tval2 > 0) == Pos);
857    boolean test3 = (tval3 == 0) || ((tval3 > 0) == Pos);
858
859    if (!test1 && !test2) {      // Go UP & RIGHT
860       insideTri = false;
861    }
862    else if (!test2 && !test3) { // Go DOWN & LEFT
863       insideTri = false;
864    }
865    else if (!test1 && !test3) { // Go UP & LEFT
866       insideTri = false;
867    }
868    else if (!test1) {           // Go UP
869       insideTri = false;
870    }
871    else if (!test3) {           // Go LEFT
872       insideTri = false;
873    }
874
875    insideTri = (insideTri && test2);
876 
877    float gxx = Float.NaN;
878    float gyy = Float.NaN;
879
880    if (insideTri) {
881       gxx = ((gx-v0x)*(v2y-v0y)
882                        + (v0y-gy)*(v2x-v0x))
883                       / ((v1x-v0x)*(v2y-v0y)
884                        + (v0y-v1y)*(v2x-v0x));
885
886       gyy =  ((gx-v0x)*(v1y-v0y)
887                        + (v0y-gy)*(v1x-v0x))
888                       / ((v2x-v0x)*(v1y-v0y)
889                        + (v0y-v2y)*(v1x-v0x)); 
890    }
891    else {
892       insideTri = true;
893
894       float bcx = v3x-v1x;
895       float bcy = v3y-v1y;
896       float cdx = v2x-v3x;
897       float cdy = v2y-v3y;
898       float cpx = gx-v3x;
899       float cpy = gy-v3y;
900
901       tval1 = bcx*bpy-bcy*bpx;
902       tval2 = cdx*cpy-cdy*cpx;
903       tval3 = bdx*dpy-bdy*dpx;
904       test1 = (tval1 == 0) || ((tval1 > 0) == Pos);
905       test2 = (tval2 == 0) || ((tval2 > 0) == Pos);
906       test3 = (tval3 == 0) || ((tval3 < 0) == Pos);
907       if (!test1 && !test3) {      // Go UP & RIGHT
908          insideTri = false;
909       }
910       else if (!test2 && !test3) { // Go DOWN & LEFT
911          insideTri = false;
912       }
913       else if (!test1 && !test2) { // Go DOWN & RIGHT
914          insideTri = false;
915       }
916       else if (!test1) {           // Go RIGHT
917          insideTri = false;
918       }
919       else if (!test2) {           // Go DOWN
920          insideTri = false;
921       }
922       insideTri = (insideTri && test3);
923
924       if (insideTri) {
925            // Found correct grid triangle
926            // Solve the point with the reverse interpolation
927            gxx = ((v3x-gx)*(v1y-v3y)
928                        + (gy-v3y)*(v1x-v3x))
929                       / ((v2x-v3x)*(v1y-v3y)
930                        - (v2y-v3y)*(v1x-v3x)) + 1;
931            gyy = ((v2y-v3y)*(v3x-gx)
932                        + (v2x-v3x)*(gy-v3y))
933                       / ((v1x-v3x)*(v2y-v3y)
934                        - (v2x-v3x)*(v1y-v3y)) + 1;
935       }
936
937    }
938
939    // bilinear interpolation within the unit square:
940
941    float intrp = values[0]*(1f-gxx)*(1f-gyy) + values[2]*(1f-gxx)*gyy + values[3]*gxx*gyy + values[1]*gxx*(1f-gyy);
942
943
944    return intrp;
945 }
946
947}