Source code

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.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 = 0;
078  public static final int BILINEAR_VISAD = 2;
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
189    projectSwathToGrid();
190    
191    initGrid();
192
193    getBoundingQuadAtGridPts();
194
195    interpolateToGrid();
196    
197    if (filter) {
198       grdFF.setSamples(filter(), false);
199    }
200    else {
201       grdFF.setSamples(gridRange, false);
202    }
203    
204    return grdFF;
205  }
206
207  private FlatField reproject(FlatField swath, int mode, boolean filter) throws Exception {
208     this.mode = mode;
209     swathRange = swath.getFloats(false);
210     ftype = (FunctionType) swath.getType();
211     
212     initGrid();
213     
214     getBoundingQuadAtGridPts();
215     
216     interpolateToGrid();
217     
218     if (filter) {
219       grdFF.setSamples(filter(), false);
220     }
221     else {
222       grdFF.setSamples(gridRange, false);
223     }
224     
225     return grdFF;
226  }
227  
228   private void getBoundingQuadAtGridPts() throws VisADException, RemoteException {
229    int ystart = 3;
230    int ystop = gridYLen-4;
231    int subLen = ((ystop - ystart)+1)/numProc;
232    int rem = ((ystop - ystart)+1) % numProc;
233    
234    
235    ThreadManager threadManager = new ThreadManager("getBoundingQuadAtGridPts");
236    for (int i=0; i<numProc; i++) {
237        final int start = i*subLen + ystart;
238        final int stop = (i != numProc-1 ) ? (start + subLen - 1): (start + subLen + rem - 1);
239          threadManager.addRunnable(new ThreadManager.MyRunnable() {
240                  public void run()  throws Exception {
241                     getBoundingQuadAtGridPts(start, stop);
242                  }
243              });
244    }
245    
246    if (numProc == 1 || !doParallel) {
247       threadManager.runSequentially();
248    }
249    else {
250       threadManager.runAllParallel();
251    }
252  }
253
254  // start to stop inclusive
255  private void getBoundingQuadAtGridPts(int grdYstart, int grdYstop) {
256
257    for (int j=grdYstart; j<=grdYstop; j++) {
258       for (int i=3; i<gridXLen-3; i++) {
259          int grdIdx = i + j*gridXLen;
260
261          int ll = findSwathGridLoc(grdIdx, swathGridCoord, gridYLen, gridXLen, "LL");
262          quads[j][i][0] = ll;
263
264          int lr = findSwathGridLoc(grdIdx, swathGridCoord, gridYLen, gridXLen, "LR");
265          quads[j][i][1] = lr;
266
267          int ul = findSwathGridLoc(grdIdx, swathGridCoord, gridYLen, gridXLen, "UL");
268          quads[j][i][2] = ul;
269
270          int ur = findSwathGridLoc(grdIdx, swathGridCoord, gridYLen, gridXLen, "UR");
271          quads[j][i][3] = ur;
272       }
273    }
274  }
275  
276  public void interpolateToGrid() throws VisADException, RemoteException {
277    int ystart = 3;
278    int ystop = gridYLen-4;
279    int subLen = ((ystop - ystart)+1)/numProc;
280    int rem = ((ystop - ystart)+1) % numProc;
281    
282    ThreadManager threadManager = new ThreadManager("interpolateToGrid");
283    for (int i=0; i<numProc; i++) {
284        final int start = i*subLen + ystart;
285        final int stop = (i != numProc-1 ) ? (start + subLen - 1): (start + subLen + rem - 1);
286          threadManager.addRunnable(new ThreadManager.MyRunnable() {
287                  public void run()  throws Exception {
288                     interpolateToGrid(start, stop);
289                  }
290              });
291    }
292    
293    if (numProc == 1 || !doParallel) {
294       threadManager.runSequentially();
295    }
296    else {
297       threadManager.runAllParallel();
298    }
299  }
300
301  // start to stop inclusive
302  public void interpolateToGrid(int grdYstart, int grdYstop) throws VisADException, RemoteException {
303
304    float[][] corners = new float[2][4];
305    float[][] rngVals = new float[rngTupDim][4];
306    float[] values = new float[4];
307    float gx;
308    float gy;
309
310    for (int j=grdYstart; j<=grdYstop; j++) {
311       for (int i=3; i<gridXLen-3; i++) {
312          int grdIdx = i + j*gridXLen;
313          gx = (float) (grdIdx % gridXLen);
314          gy = (float) (grdIdx / gridXLen);
315
316          java.util.Arrays.fill(corners[0], Float.NaN);
317          java.util.Arrays.fill(corners[1], Float.NaN);
318        
319          int ll = quads[j][i][0];
320          int lr = quads[j][i][1];
321          int ul = quads[j][i][2];
322          int ur = quads[j][i][3];
323
324          if (ll >= 0) {
325             corners[0][0] = swathGridCoord[0][ll] - gx;
326             corners[1][0] = swathGridCoord[1][ll] - gy;
327          }
328          if (lr >= 0) {
329             corners[0][1] = swathGridCoord[0][lr] - gx;
330             corners[1][1] = swathGridCoord[1][lr] - gy;
331          }
332          if (ul >= 0) {
333             corners[0][2] = swathGridCoord[0][ul] - gx;
334             corners[1][2] = swathGridCoord[1][ul] - gy;
335          }
336          if (ur >= 0) {
337             corners[0][3] = swathGridCoord[0][ur] - gx;
338             corners[1][3] = swathGridCoord[1][ur] - gy;
339          }
340
341          if (mode == NEAREST) { // Nearest neighbor
342             for (int t=0; t<rngTupDim; t++) {
343                java.util.Arrays.fill(values, Float.NaN);
344                if (ll >=0) {
345                   values[0] = swathRange[t][swathIndexAtGrid[ll]];
346                }
347                if (lr >= 0) {
348                   values[1] = swathRange[t][swathIndexAtGrid[lr]];
349                }
350                if (ul >= 0) {
351                   values[2] = swathRange[t][swathIndexAtGrid[ul]];
352                }
353                if (ur >= 0) {
354                   values[3] = swathRange[t][swathIndexAtGrid[ur]];
355                }
356                float val = nearest(0f, 0f, corners, values);
357                gridRange[t][grdIdx] = val;
358             }
359          }
360          else if (mode == BILINEAR_VISAD) {  //from VisAD
361             if (!(ll >= 0 && lr >= 0 && ul >= 0 && ur >= 0)) {
362                continue;
363             }
364             corners[0][0] = swathGridCoord[0][ll];
365             corners[1][0] = swathGridCoord[1][ll];
366             corners[0][1] = swathGridCoord[0][lr];
367             corners[1][1] = swathGridCoord[1][lr];
368             corners[0][2] = swathGridCoord[0][ul];
369             corners[1][2] = swathGridCoord[1][ul];
370             corners[0][3] = swathGridCoord[0][ur];
371             corners[1][3] = swathGridCoord[1][ur];
372             for (int t=0; t<rngTupDim; t++) {
373                java.util.Arrays.fill(values, Float.NaN);
374                values[0] = swathRange[t][swathIndexAtGrid[ll]];
375                values[1] = swathRange[t][swathIndexAtGrid[lr]];
376                values[2] = swathRange[t][swathIndexAtGrid[ul]];
377                values[3] = swathRange[t][swathIndexAtGrid[ur]];
378                float val = visad2D(gy, gx, corners, values);
379                gridRange[t][grdIdx] = val;
380             }
381          }
382          else if (mode == 1) { //TODO: not working yet
383             if (!(ll >= 0 && lr >= 0 && ul >= 0 && ur >= 0)) {
384                continue;
385             }
386             gx -= swathGridCoord[0][ll];
387             gy -= swathGridCoord[1][ll];
388             corners[0][0] = 0f;
389             corners[1][0] = 0f;
390             corners[0][1] = swathGridCoord[0][ul] - swathGridCoord[0][ll];
391             corners[1][1] = swathGridCoord[1][ul] - swathGridCoord[1][ll];
392             corners[0][2] = swathGridCoord[0][ur] - swathGridCoord[0][ll];
393             corners[1][2] = swathGridCoord[1][ur] - swathGridCoord[1][ll];
394             corners[0][3] = swathGridCoord[0][lr] - swathGridCoord[0][ll];
395             corners[1][3] = swathGridCoord[1][lr] - swathGridCoord[1][ll];
396             for (int t=0; t<rngTupDim; t++) {
397                java.util.Arrays.fill(values, Float.NaN);
398                values[0] = swathRange[t][swathIndexAtGrid[ll]];
399                values[1] = swathRange[t][swathIndexAtGrid[ul]];
400                values[2] = swathRange[t][swathIndexAtGrid[ur]];
401                values[3] = swathRange[t][swathIndexAtGrid[lr]];
402                float val = biLinearIntrp(gy, gx, corners, values);
403                gridRange[t][grdIdx] = val;
404             }
405
406          }
407       }
408    }
409
410  }
411
412 public void projectSwathToGrid() throws VisADException, RemoteException {
413    int subLen = trackLen/numProc;
414    int rem = trackLen % numProc;
415    
416    ThreadManager threadManager = new ThreadManager("projectSwathToGrid");
417    for (int i=0; i<numProc; i++) {
418        final int start = i*subLen;
419        final int stop = (i != numProc-1 ) ? (start + subLen - 1): (start + subLen + rem - 1);
420          threadManager.addRunnable(new ThreadManager.MyRunnable() {
421                  public void run()  throws Exception {
422                     projectSwathToGrid(start, stop);
423                  }
424              });
425    }
426    
427    if (numProc == 1 || !doParallel) {
428       threadManager.runSequentially();
429    }
430    else {
431       threadManager.runAllParallel();
432    }
433 }
434 
435 public void projectSwathToGrid(int trackStart, int trackStop) throws VisADException, RemoteException {
436
437    for (int j=trackStart; j <= trackStop; j++) {
438       for (int i=0; i < xtrackLen; i++) {
439         int swathIdx = j*xtrackLen + i;
440
441         float[][] swathCoord = swathDomain.indexToValue(new int[] {swathIdx});
442         float[][] swathEarthCoord = swathCoordSys.toReference(swathCoord);
443
444         float[][] gridValue = gridCoordSys.fromReference(swathEarthCoord);
445         float[][] gridCoord = grid.valueToGrid(gridValue);
446         float g0 = gridCoord[0][0];
447         float g1 = gridCoord[1][0];
448         int grdIdx = (g0 != g0 || g1 != g1) ? -1 : ((int) g0) + gridXLen * ((int) g1);
449         int m=0;
450         int n=0;
451         int k = grdIdx + (m + n*gridXLen);
452         
453         allSwathGridCoords[0][swathIdx] = g0;
454         allSwathGridCoords[1][swathIdx] = g1;
455         allSwathGridIndexs[swathIdx] = k;
456             
457       }
458    }
459 }
460 
461 public void initGrid() throws VisADException {
462    Arrays.fill(swathGridCoord[0], Float.NaN);
463    Arrays.fill(swathGridCoord[1], Float.NaN);
464    Arrays.fill(swathIndexAtGrid, -1);
465
466    for (int j=0; j < trackLen; j++) {
467       for (int i=0; i < xtrackLen; i++) {
468         int swathIdx = j*xtrackLen + i;
469         float val = swathRange[0][swathIdx];
470         int k = allSwathGridIndexs[swathIdx];
471         
472         if ( !(Float.isNaN(val)) && ((k >=0) && (k < gridLen)) ) { // val or val[rngTupDim] ?
473            if (swathIndexAtGrid[k] == -1) {
474               swathGridCoord[0][k] = allSwathGridCoords[0][swathIdx];
475               swathGridCoord[1][k] = allSwathGridCoords[1][swathIdx];
476               swathIndexAtGrid[k] = swathIdx;
477            }
478         }
479       }
480    }
481    
482    for (int j=0; j<gridYLen; j++) {
483       for (int i=0; i<gridXLen; i++) {
484          java.util.Arrays.fill(quads[j][i], -1);
485       }
486    }
487    
488    RealTupleType rtt = ((SetType)grid.getType()).getDomain();
489    grdFF = new FlatField(new FunctionType(rtt, ftype.getRange()), grid);
490    gridRange = grdFF.getFloats(false);
491    rngTupDim = gridRange.length;
492    for (int t=0; t<rngTupDim; t++) {
493       java.util.Arrays.fill(gridRange[t], Float.NaN);
494    }
495 }
496
497 private float[][] filter() throws VisADException, RemoteException {
498
499    double mag = 3.0;
500    double sigma = 0.4;
501
502    float[][] weights = new float[3][3];
503
504    float sumWeights = 0f;
505    for (int n=-1; n<=1; n++) {
506       for (int m=-1; m<=1; m++) {
507          double del_0 = m;
508          double del_1 = n;
509          double dst = Math.sqrt(del_0*del_0 + del_1*del_1);
510
511          weights[n+1][m+1] = (float) (mag/Math.exp(dst/sigma));
512
513          sumWeights += weights[n+1][m+1];
514       }
515    }
516
517    for (int n=-1; n<=1; n++) {
518       for (int m=-1; m<=1; m++) {
519          weights[n+1][m+1] /= sumWeights;
520        }
521    }
522
523    float[][] newRange = new float[rngTupDim][gridLen];
524    for (int t=0; t<rngTupDim; t++) {
525       java.util.Arrays.fill(newRange[t], Float.NaN);
526    }
527    float[] sum = new float[rngTupDim];
528
529    for (int j=2; j<gridYLen-2; j++) {
530       for (int i=2; i<gridXLen-2; i++) {
531         int grdIdx = i + j*gridXLen;
532
533         java.util.Arrays.fill(sum, 0f);
534         for (int n=-1; n<=1; n++) {
535            for (int m=-1; m<=1; m++) {
536               int k = grdIdx + (m + n*gridXLen);
537
538               for (int t=0; t<rngTupDim; t++) {
539                  sum[t] += weights[n+1][m+1]*gridRange[t][k];
540               }
541            }
542         }
543
544         for (int t=0; t<rngTupDim; t++) {
545            newRange[t][grdIdx] = sum[t];
546         }
547       }
548    }
549
550    return newRange;
551 }
552
553 private static int findSwathGridLoc(int grdIdx, float[][] swathGridCoord, int gridYLen, int gridXLen, String which) {
554  
555    int idx = -1;
556
557    int gy = grdIdx/gridXLen;
558    int gx = grdIdx % gridXLen;
559
560    switch (which) {
561       case "LL":
562
563          idx = (gy-1)*gridXLen + (gx-1);
564          if (!Float.isNaN(swathGridCoord[0][idx])) {
565             break;
566          }
567          idx = (gy-2)*gridXLen + (gx-1);
568          if (!Float.isNaN(swathGridCoord[0][idx])) {
569             break;
570          }
571          idx = (gy-1)*gridXLen + (gx-2);
572          if (!Float.isNaN(swathGridCoord[0][idx])) {
573             break;
574          }
575          idx = (gy-2)*gridXLen + (gx-2);
576          if (!Float.isNaN(swathGridCoord[0][idx])) {
577             break;
578          }
579
580
581          idx = (gy-2)*gridXLen + (gx-3);
582          if (!Float.isNaN(swathGridCoord[0][idx])) {
583             break;
584          }
585          idx = (gy-3)*gridXLen + (gx-2);
586          if (!Float.isNaN(swathGridCoord[0][idx])) {
587             break;
588          }
589          idx = (gy-3)*gridXLen + (gx-1);
590          if (!Float.isNaN(swathGridCoord[0][idx])) {
591             break;
592          }
593          idx = (gy-1)*gridXLen + (gx-3);
594          if (!Float.isNaN(swathGridCoord[0][idx])) {
595             break;
596          }
597          idx = (gy-3)*gridXLen + (gx-3);
598          if (!Float.isNaN(swathGridCoord[0][idx])) {
599             break;
600          }
601
602          idx = -1;
603          break;
604       case "UL":
605          idx = (gy)*gridXLen + (gx-1);
606          if (!Float.isNaN(swathGridCoord[0][idx])) {
607             break;
608          }
609          idx = (gy)*gridXLen + (gx-2);
610          if (!Float.isNaN(swathGridCoord[0][idx])) {
611             break;
612          }
613          idx = (gy+1)*gridXLen + (gx-1);
614          if (!Float.isNaN(swathGridCoord[0][idx])) {
615             break;
616          }
617          idx = (gy+1)*gridXLen + (gx-2);
618          if (!Float.isNaN(swathGridCoord[0][idx])) {
619             break;
620          }
621
622          idx = (gy+1)*gridXLen + (gx-3);
623          if (!Float.isNaN(swathGridCoord[0][idx])) {
624             break;
625          }
626          idx = (gy+2)*gridXLen + (gx-3);
627          if (!Float.isNaN(swathGridCoord[0][idx])) {
628             break;
629          }
630          idx = (gy+2)*gridXLen + (gx-2);
631          if (!Float.isNaN(swathGridCoord[0][idx])) {
632             break;
633          }
634          idx = (gy)*gridXLen + (gx-3);
635          if (!Float.isNaN(swathGridCoord[0][idx])) {
636             break;
637          }
638          idx = (gy-1)*gridXLen + (gx-3);
639          if (!Float.isNaN(swathGridCoord[0][idx])) {
640             break;
641          }
642
643          idx = -1;
644          break;
645       case "UR":
646          idx = (gy)*gridXLen + (gx);
647          if (!Float.isNaN(swathGridCoord[0][idx])) {
648             break;
649          }
650          idx = (gy+1)*gridXLen + (gx);
651          if (!Float.isNaN(swathGridCoord[0][idx])) {
652             break;
653          }
654          idx = (gy)*gridXLen + (gx+1);
655          if (!Float.isNaN(swathGridCoord[0][idx])) {
656             break;
657          }
658          idx = (gy+1)*gridXLen + (gx+1);
659          if (!Float.isNaN(swathGridCoord[0][idx])) {
660             break;
661          }
662
663          idx = (gy+2)*gridXLen + (gx+1);
664          if (!Float.isNaN(swathGridCoord[0][idx])) {
665             break;
666          }
667          idx = (gy+1)*gridXLen + (gx+2);
668          if (!Float.isNaN(swathGridCoord[0][idx])) {
669             break;
670          }
671          idx = (gy+2)*gridXLen + (gx+2);
672          if (!Float.isNaN(swathGridCoord[0][idx])) {
673             break;
674          }
675          idx = (gy)*gridXLen + (gx+2);
676          if (!Float.isNaN(swathGridCoord[0][idx])) {
677             break;
678          }
679          idx = (gy-1)*gridXLen + (gx+2);
680          if (!Float.isNaN(swathGridCoord[0][idx])) {
681             break;
682          }
683
684          idx = -1;
685          break;
686       case "LR":
687          idx = (gy-1)*gridXLen + (gx);
688          if (!Float.isNaN(swathGridCoord[0][idx])) {
689             break;
690          }
691          idx = (gy-1)*gridXLen + (gx+1);
692          if (!Float.isNaN(swathGridCoord[0][idx])) {
693             break;
694          }
695          idx = (gy-2)*gridXLen + (gx);
696          if (!Float.isNaN(swathGridCoord[0][idx])) {
697             break;
698          }
699          idx = (gy-2)*gridXLen + (gx+1);
700          if (!Float.isNaN(swathGridCoord[0][idx])) {
701             break;
702          }
703
704          idx = (gy-1)*gridXLen + (gx+2);
705          if (!Float.isNaN(swathGridCoord[0][idx])) {
706             break;
707          }
708          idx = (gy-2)*gridXLen + (gx+2);
709          if (!Float.isNaN(swathGridCoord[0][idx])) {
710             break;
711          }
712          idx = (gy-3)*gridXLen + (gx);
713          if (!Float.isNaN(swathGridCoord[0][idx])) {
714             break;
715          }
716          idx = (gy-3)*gridXLen + (gx+1);
717          if (!Float.isNaN(swathGridCoord[0][idx])) {
718             break;
719          }
720          idx = (gy-3)*gridXLen + (gx+2);
721          if (!Float.isNaN(swathGridCoord[0][idx])) {
722             break;
723          }
724 
725          idx = -1;
726          break;
727    }
728
729    return idx;
730 }
731
732 /* Reference: David W. Zingg, University of Toronto, Downsview, Ontario, Canada
733               Maurice Yarrow, Sterling Software, Arnes Research Center, Moffett Field, California. 
734               NASA Technical Memorandum 102213
735
736     y -> q, x -> p; first point (x=0, y=0) and clockwise around
737  */
738 public static float biLinearIntrp(float gy, float gx, float[][] corners, float[] values) {
739    // transform physical space (gy, gx) to unit square (q, p)
740    // bilinear mapping coefficients
741/*
742    float a = corners[0][3];
743    float b = corners[0][1];
744    float c = (corners[0][2] - corners[0][3] - corners[0][1]);
745
746    float d = corners[1][3];
747    float e = corners[1][1];
748    float f = (corners[1][2] - corners[1][3] - corners[1][1]);
749*/
750
751    float a = corners[0][1];
752    float b = corners[0][3];
753    float c = (corners[0][2] - corners[0][1] - corners[0][3]);
754
755    float d = corners[1][1];
756    float e = corners[1][3];
757    float f = (corners[1][2] - corners[1][1] - corners[1][3]);
758
759
760    // quadratic terms to determine p
761    // p = (-coef1 +/- sqrt(coef1^2 - 4*coef2*coef0))/2*coef2
762
763    float coef2 = (c*d - a*f);  // A 
764    float coef1 = (-c*gy + b*d + gx*f - a*e);  // B
765    float coef0 = (-gy*b + gx*e);  // C
766
767    // two p vals from quadratic:
768    float p0 =  (-coef1 + ((float)Math.sqrt(coef1*coef1 - 4*coef2*coef0)) )/2f*coef2;
769    float p1 =  (-coef1 - ((float)Math.sqrt(coef1*coef1 - 4*coef2*coef0)) )/2f*coef2;
770
771    // the corresponding q values for both p solutions:
772    float q0 = (gx - a*p0)/(b + c*p0);
773    float q1 = (gx - a*p1)/(b + c*p1);
774
775    // test  which point to use. One will be outside unit square:
776    float p = Float.NaN;
777    float q = Float.NaN;
778    if ((p0 >= 0f && p0 <= 1f) && (q0 >= 0f && q0 <= 1f)) {
779       p = p0;
780       q = q0;
781    }
782    else if ((p1 >= 0f && p1 <= 1f) && (q1 >= 0f && q1 <= 1f)) {
783       p = p1;
784       q = q1;
785    }
786
787    // bilinear interpolation within the unit square:
788
789    float intrp = values[0]*(1f-p)*(1f-q) + values[1]*(1f-p)*q + values[2]*p*q + values[3]*p*(1f-q);
790
791    return intrp;
792 }
793
794 private static float nearest(float gy, float gx, float[][] corners, float[] values) {
795   float minDist = Float.MAX_VALUE;
796
797   float delx;
798   float dely;
799   int closest = 0;
800   for (int k=0; k<4; k++) {
801      delx = corners[0][k] - gx;
802      dely = corners[1][k] - gy;
803      float dst_sqrd = delx*delx + dely*dely;
804
805      if (dst_sqrd < minDist) {
806         minDist = dst_sqrd;
807         closest = k;
808      }
809   }
810
811   return values[closest];
812 }
813
814 private static float[] v0 = new float[2];  // A
815 private static float[] v1 = new float[2];  // B
816 private static float[] v2 = new float[2];  // D
817 private static float[] v3 = new float[2];  // C
818
819 private static float[] bd = new float[2];
820 private static float[] bp = new float[2];
821 private static float[] dp = new float[2];
822
823 private static float[] ab = new float[2];
824 private static float[] da = new float[2];
825 private static float[] ap = new float[2];
826
827 private static float[] bc = new float[2];
828 private static float[] cd = new float[2];
829 private static float[] cp = new float[2];
830
831
832 public static float visad2D(float gy, float gx, float[][] corners, float[] values) {
833
834    boolean Pos = true;
835
836    v0[0] = corners[0][0];
837    v0[1] = corners[1][0];
838    v1[0] = corners[0][1];
839    v1[1] = corners[1][1];
840    v2[0] = corners[0][2];
841    v2[1] = corners[1][2];
842    v3[0] = corners[0][3];
843    v3[1] = corners[1][3];
844
845    bd[0] = v2[0]-v1[0];
846    bd[1] = v2[1]-v1[1];
847
848    bp[0] = gx-v1[0];
849    bp[1] = gy-v1[1];
850
851    dp[0] = gx-v2[0];
852    dp[1] = gy-v2[1];
853
854    // lower triangle first
855
856    boolean insideTri = true;
857
858    ab[0] = v1[0]-v0[0];
859    ab[1] = v1[1]-v0[1];
860    da[0] = v0[0]-v2[0];
861    da[1] = v0[1]-v2[1];
862    ap[0] = gx-v0[0];
863    ap[1] = gy-v0[1];
864    float tval1 = ab[0]*ap[1]-ab[1]*ap[0];
865    float tval2 = bd[0]*bp[1]-bd[1]*bp[0];
866    float tval3 = da[0]*dp[1]-da[1]*dp[0];
867    boolean test1 = (tval1 == 0) || ((tval1 > 0) == Pos);
868    boolean test2 = (tval2 == 0) || ((tval2 > 0) == Pos);
869    boolean test3 = (tval3 == 0) || ((tval3 > 0) == Pos);
870
871    if (!test1 && !test2) {      // Go UP & RIGHT
872       insideTri = false;
873    }
874    else if (!test2 && !test3) { // Go DOWN & LEFT
875       insideTri = false;
876    }
877    else if (!test1 && !test3) { // Go UP & LEFT
878       insideTri = false;
879    }
880    else if (!test1) {           // Go UP
881       insideTri = false;
882    }
883    else if (!test3) {           // Go LEFT
884       insideTri = false;
885    }
886
887    insideTri = (insideTri && test2);
888 
889    float gxx = Float.NaN;
890    float gyy = Float.NaN;
891
892    if (insideTri) {
893       gxx = ((gx-v0[0])*(v2[1]-v0[1])
894                        + (v0[1]-gy)*(v2[0]-v0[0]))
895                       / ((v1[0]-v0[0])*(v2[1]-v0[1])
896                        + (v0[1]-v1[1])*(v2[0]-v0[0]));
897
898       gyy =  ((gx-v0[0])*(v1[1]-v0[1])
899                        + (v0[1]-gy)*(v1[0]-v0[0]))
900                       / ((v2[0]-v0[0])*(v1[1]-v0[1])
901                        + (v0[1]-v2[1])*(v1[0]-v0[0])); 
902    }
903    else {
904       insideTri = true;
905
906       bc[0] = v3[0]-v1[0];
907       bc[1] = v3[1]-v1[1];
908       cd[0] = v2[0]-v3[0];
909       cd[1] = v2[1]-v3[1];
910       cp[0] = gx-v3[0];
911       cp[1] = gy-v3[1];
912
913       tval1 = bc[0]*bp[1]-bc[1]*bp[0];
914       tval2 = cd[0]*cp[1]-cd[1]*cp[0];
915       tval3 = bd[0]*dp[1]-bd[1]*dp[0];
916       test1 = (tval1 == 0) || ((tval1 > 0) == Pos);
917       test2 = (tval2 == 0) || ((tval2 > 0) == Pos);
918       test3 = (tval3 == 0) || ((tval3 < 0) == Pos);
919       if (!test1 && !test3) {      // Go UP & RIGHT
920          insideTri = false;
921       }
922       else if (!test2 && !test3) { // Go DOWN & LEFT
923          insideTri = false;
924       }
925       else if (!test1 && !test2) { // Go DOWN & RIGHT
926          insideTri = false;
927       }
928       else if (!test1) {           // Go RIGHT
929          insideTri = false;
930       }
931       else if (!test2) {           // Go DOWN
932          insideTri = false;
933       }
934       insideTri = (insideTri && test3);
935
936       if (insideTri) {
937            // Found correct grid triangle
938            // Solve the point with the reverse interpolation
939            gxx = ((v3[0]-gx)*(v1[1]-v3[1])
940                        + (gy-v3[1])*(v1[0]-v3[0]))
941                       / ((v2[0]-v3[0])*(v1[1]-v3[1])
942                        - (v2[1]-v3[1])*(v1[0]-v3[0])) + 1;
943            gyy = ((v2[1]-v3[1])*(v3[0]-gx)
944                        + (v2[0]-v3[0])*(gy-v3[1]))
945                       / ((v1[0]-v3[0])*(v2[1]-v3[1])
946                        - (v2[0]-v3[0])*(v1[1]-v3[1])) + 1;
947       }
948
949    }
950
951    // bilinear interpolation within the unit square:
952
953    float intrp = values[0]*(1f-gxx)*(1f-gyy) + values[2]*(1f-gxx)*gyy + values[3]*gxx*gyy + values[1]*gxx*(1f-gyy);
954
955
956    return intrp;
957 }
958
959
960}