Source code

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