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 */
028
029package edu.wisc.ssec.mcidasv.data.hydra;
030
031import java.rmi.RemoteException;
032import visad.Set;
033import visad.Gridded1DSet;
034import visad.CoordinateSystem;
035import visad.RealType;
036import visad.RealTupleType;
037import visad.Linear1DSet;
038import visad.Linear2DSet;
039import visad.Gridded2DSet;
040import visad.SampledSet;
041import visad.Unit;
042import visad.FunctionType;
043import visad.VisADException;
044import visad.QuickSort;
045import visad.FlatField;
046import visad.FieldImpl;
047import java.util.HashMap;
048import java.util.Map;
049import visad.GriddedSet;
050
051
052public abstract class ProfileAlongTrack extends MultiDimensionAdapter {
053
054      private FunctionType mathtype;
055
056      int TrackLen;
057      int VertLen;
058
059      private float[] vertLocs = null;
060      private float[] trackTimes = null;
061      private float[] trackLongitude = null;
062      private float[] trackLatitude = null;
063
064      public static String longitude_name = "Longitude";
065      public static String latitude_name  = "Latitude";
066      public static String trackDim_name  = "TrackDim";
067      public static String vertDim_name  = "VertDim";
068      public static String array_name = "array_name";
069      public static String profileTime_name = "ProfileTime";
070      public static String profileTime_unit = "ProfileTime_unit";
071      public static String altitude_unit = "altitude_unit";
072      public static String sfcElev_name = "SurfaceElev";
073      public static String range_name = "range_name";
074      public static String scale_name = "scale_name";
075      public static String offset_name = "offset_name";
076      public static String fill_value_name = "fill_value_name";
077      public static String valid_range = "valid_range";
078      public static String ancillary_file_name = "ancillary_file";
079      static String product_name = "product_name";
080      
081      String[] rangeName_s  = null;
082      Class[] arrayType_s = null;
083      Unit[] rangeUnit_s  = new Unit[] {null};
084
085      RealType track  = RealType.getRealType(trackDim_name);
086      RealType vert = RealType.getRealType(vertDim_name);
087      RealType[] domainRealTypes = new RealType[2];
088
089      RealType vertLocType;
090      RealType trackTimeType;
091
092      int track_idx      = -1;
093      int vert_idx       = -1;
094      int range_rank     = -1;
095
096      int track_tup_idx;
097      int vert_tup_idx;
098
099      boolean isVertDimAlt = true;
100
101      CoordinateSystem cs = null;
102      
103      int medianFilterTrackWidth = 10;
104      int medianFilterVertWidth = 10;
105
106      public static Map<String, double[]> getEmptySubset() {
107        Map<String, double[]> subset = new HashMap<>();
108        subset.put(trackDim_name, new double[3]);
109        subset.put(vertDim_name, new double[3]);
110        return subset;
111      }
112
113      public static Map<String, Object> getEmptyMetadataTable() {
114        Map<String, Object> metadata = new HashMap<>();
115        metadata.put(array_name, null);
116        metadata.put(trackDim_name, null);
117        metadata.put(vertDim_name, null);
118        metadata.put(longitude_name, null);
119        metadata.put(latitude_name, null);
120        metadata.put(profileTime_name, null);
121        metadata.put(profileTime_unit, null);
122        metadata.put(altitude_unit, null);
123        metadata.put(sfcElev_name, null);
124        metadata.put(scale_name, null);
125        metadata.put(offset_name, null);
126        metadata.put(fill_value_name, null);
127        /*
128        metadata.put(range_name, null);
129        metadata.put(range_unit, null);
130        metadata.put(valid_range, null);
131        */
132        return metadata;
133      }
134
135      public ProfileAlongTrack() {
136      }
137
138      public ProfileAlongTrack(MultiDimensionReader reader, Map<String, Object> metadata) {
139        this(reader, metadata, true);
140      }
141
142      public ProfileAlongTrack(MultiDimensionReader reader, Map<String, Object> metadata, boolean isVertDimAlt) {
143        super(reader, metadata);
144        this.isVertDimAlt = isVertDimAlt;
145        this.init();
146      }
147
148
149      private void init() {
150        for (int k=0; k<array_rank;k++) {
151          if ( ((String)metadata.get(trackDim_name)).equals(array_dim_names[k]) ) {
152            track_idx = k;
153          }
154          if ( ((String)metadata.get(vertDim_name)).equals(array_dim_names[k]) ) {
155            vert_idx = k;
156          }
157        }
158
159        int[] lengths = new int[2];
160
161        if (track_idx < vert_idx) {
162          domainRealTypes[0] = vert;
163          domainRealTypes[1] = track;
164          track_tup_idx = 1;
165          vert_tup_idx = 0;
166          lengths[0] = array_dim_lengths[vert_idx];
167          lengths[1] = array_dim_lengths[track_idx];
168        }
169        else {
170          domainRealTypes[0] = track;
171          domainRealTypes[1] = vert;
172          track_tup_idx = 0;
173          vert_tup_idx = 1;
174          lengths[0] = array_dim_lengths[track_idx];
175          lengths[1] = array_dim_lengths[vert_idx];
176        }
177
178        TrackLen = array_dim_lengths[track_idx];
179        VertLen = array_dim_lengths[vert_idx];
180
181        String rangeName = null;
182        if (metadata.get("range_name") != null) {
183          rangeName = (String)metadata.get("range_name");
184        } 
185        else {
186          rangeName = (String)metadata.get("array_name");
187        }
188        rangeType = RealType.getRealType(rangeName, rangeUnit_s[0]);
189
190        try {
191          rangeProcessor = RangeProcessor.createRangeProcessor(reader, metadata);
192        } 
193        catch (Exception e) {
194          System.out.println("RangeProcessor failed to create");
195          e.printStackTrace();
196        }
197
198        try {
199          if (isVertDimAlt) {
200            vertLocs = getVertBinAltitude();
201          }
202          vertLocType = makeVertLocType();
203          trackTimes = getTrackTimes();
204          trackTimeType = makeTrackTimeType();
205          trackLongitude = getTrackLongitude();
206          trackLatitude = getTrackLatitude();
207        } 
208        catch (Exception e) {
209          System.out.println(e);
210        }
211        
212      }
213
214      public int getTrackLength() {
215        return TrackLen;
216      }
217
218      public int getVertLength() {
219        return VertLen;
220      }
221
222      public int getVertIdx() {
223        return vert_idx;
224      }
225
226      public int getTrackIdx() {
227        return track_idx;
228      }
229
230      public int getVertTupIdx() {
231        return vert_tup_idx;
232      }
233
234      public int getTrackTupIdx() {
235        return track_tup_idx;
236      }
237      
238      public int getMedianFilterWindowWidth() {
239        return medianFilterTrackWidth;
240      }
241      
242      public int getMedianFilterWindowHeight() {
243        return medianFilterVertWidth;
244      }
245                                                                                                                                                     
246      public Set makeDomain(Map<String, double[]> subset) throws Exception {
247        double[] first = new double[2];
248        double[] last = new double[2];
249        int[] length = new int[2];
250
251        Map<String, double[]> domainSubset = new HashMap<>();
252        domainSubset.put(trackDim_name, subset.get(trackDim_name));
253        domainSubset.put(vertDim_name, subset.get(vertDim_name));
254
255        for (int kk=0; kk<2; kk++) {
256          RealType rtype = domainRealTypes[kk];
257          double[] coords = subset.get(rtype.getName());
258          first[kk] = coords[0];
259          last[kk] = coords[1];
260          length[kk] = (int) ((last[kk] - first[kk])/coords[2] + 1);
261          last[kk] = first[kk]+coords[2]*(length[kk]-1);
262        }
263        Linear2DSet domainSet = new Linear2DSet(first[0], last[0], length[0], first[1], last[1], length[1]);
264        final Linear1DSet[] lin1DSet_s = new Linear1DSet[] {domainSet.getLinear1DComponent(0),
265                                           domainSet.getLinear1DComponent(1)};
266
267        float[] new_altitudes = new float[length[vert_tup_idx]];
268        float[] new_times = new float[length[track_tup_idx]];
269
270        int track_start = (int) first[track_tup_idx];
271        int vert_start = (int) first[vert_tup_idx];
272        int vert_skip = (int) (subset.get(vertDim_name))[2];
273        int track_skip = (int) (subset.get(trackDim_name))[2];
274        for (int k=0; k<new_altitudes.length; k++) {
275          new_altitudes[k] = vertLocs[vert_start+(k*vert_skip)];
276        }
277
278        for (int k=0; k<new_times.length; k++) {
279          new_times[k] = trackTimes[track_start+(k*track_skip)];
280        }
281
282        final Gridded1DSet alt_set = new Gridded1DSet(vertLocType, new float[][] {new_altitudes}, new_altitudes.length);
283        final Gridded1DSet time_set = new Gridded1DSet(trackTimeType, new float[][] {new_times}, new_times.length);
284        final float vert_offset = (float) first[vert_tup_idx];
285        final float track_offset = (float) first[track_tup_idx];
286
287        RealTupleType reference = new RealTupleType(vertLocType, trackTimeType);
288        
289        CoordinateSystem cs = null;
290
291        try {
292        cs = new CoordinateSystem(reference, null) {
293          public float[][] toReference(float[][] vals) throws VisADException {
294            int[] indexes = lin1DSet_s[0].valueToIndex(new float[][] {vals[0]});
295            /* ?
296            for (int k=0; k<vals[0].length;k++) {
297               indexes[k] = (int) (vals[vert_tup_idx][k] - vert_offset);
298            }
299            */
300            float[][] alts = alt_set.indexToValue(indexes);
301
302            indexes = lin1DSet_s[1].valueToIndex(new float[][] {vals[1]});
303            /* ?
304            for (int k=0; k<vals[0].length;k++) {
305               indexes[k] = (int) (vals[track_tup_idx][k] - track_offset);
306            }
307            */
308            float[][] times = time_set.indexToValue(indexes);
309
310            return new float[][] {alts[0], times[0]};
311          }
312          public float[][] fromReference(float[][] vals) throws VisADException {
313            int[] indexes = alt_set.valueToIndex(new float[][] {vals[vert_tup_idx]});
314            float[][] vert_coords = lin1DSet_s[vert_tup_idx].indexToValue(indexes);
315            indexes = time_set.valueToIndex(new float[][] {vals[track_tup_idx]});
316            float[][] track_coords = lin1DSet_s[track_tup_idx].indexToValue(indexes);
317            return new float[][] {vert_coords[0], track_coords[0]};
318          }
319          public double[][] toReference(double[][] vals) throws VisADException {
320            return Set.floatToDouble(toReference(Set.doubleToFloat(vals)));
321          }
322          public double[][] fromReference(double[][] vals) throws VisADException {
323            return Set.floatToDouble(fromReference(Set.doubleToFloat(vals)));
324          }
325          public boolean equals(Object obj) {
326            return true;
327          }
328        };
329        }
330        catch (Exception e) {
331        }
332
333        RealTupleType domainTupType = new RealTupleType(domainRealTypes[0], domainRealTypes[1], cs, null);
334        domainSet = new Linear2DSet(domainTupType, first[0], last[0], length[0], first[1], last[1], length[1]);
335
336        return domainSet;
337      }
338
339      public FunctionType getMathType() {
340        return null;
341      }
342
343      public RealType[] getDomainRealTypes() {
344        return domainRealTypes;
345      }
346
347      public Map<String, double[]> getDefaultSubset() {
348        Map<String, double[]> subset = ProfileAlongTrack.getEmptySubset();
349
350        double[] coords = subset.get("TrackDim");
351        coords[0] = 20000.0;
352        coords[1] = (TrackLen - 15000.0) - 1;
353        //-coords[2] = 30.0;
354        coords[2] = 5.0;
355        subset.put("TrackDim", coords);
356
357        coords = subset.get("VertDim");
358        coords[0] = 98.0;
359        coords[1] = (VertLen) - 1;
360        coords[2] = 2.0;
361        subset.put("VertDim", coords);
362        return subset;
363      }
364
365      public int[] getTrackRangeInsideLonLatRect(double minLat, double maxLat, double minLon, double maxLon) {
366        int nn = 100;
367        int skip = TrackLen/nn;
368        double lon;
369        double lat;
370        
371        int idx = 0;
372        while (idx < TrackLen) {
373          lon = (double)trackLongitude[idx];
374          lat = (double)trackLatitude[idx];
375          if (((lon > minLon) && (lon < maxLon)) && ((lat > minLat)&&(lat < maxLat))) break;
376          idx += skip;
377        }
378        if (idx > TrackLen-1) idx = TrackLen-1;
379        if (idx == TrackLen-1) return new int[] {-1,-1};
380
381        int low_idx = idx;
382        while (low_idx > 0) {
383          lon = (double)trackLongitude[low_idx];
384          lat = (double)trackLatitude[low_idx];
385          if (((lon > minLon) && (lon < maxLon)) && ((lat > minLat)&&(lat < maxLat))) {
386            low_idx -= 1;
387            continue;
388          }
389          else {
390            break;
391          }
392        }
393
394        int hi_idx = idx;
395        while (hi_idx < TrackLen-1) {
396          lon = (double)trackLongitude[hi_idx];
397          lat = (double)trackLatitude[hi_idx];
398          if (((lon > minLon) && (lon < maxLon)) && ((lat > minLat)&&(lat < maxLat))) {
399            hi_idx += 1;
400            continue;
401          }
402          else {
403            break;
404          }
405        }
406        return new int[] {low_idx, hi_idx};
407      }
408
409      public Map<String, double[]> getSubsetFromLonLatRect(Map<String, double[]> subset, double minLat, double maxLat, double minLon, double maxLon) {
410        double[] coords = subset.get("TrackDim");
411        int[] idxs = getTrackRangeInsideLonLatRect(minLat, maxLat, minLon, maxLon);
412        coords[0] = (double) idxs[0];
413        coords[1] = (double) idxs[1];
414        if ((coords[0] == -1) || (coords[1] == -1)) return null;
415        return subset;
416      }
417
418      public Map<String, double[]> getSubsetFromLonLatRect(Map<String, double[]> subset, double minLat, double maxLat, double minLon, double maxLon,
419                                             int xStride, int yStride, int zStride) {
420
421        double[] coords = subset.get("TrackDim");
422        int[] idxs = getTrackRangeInsideLonLatRect(minLat, maxLat, minLon, maxLon);
423        coords[0] = (double) idxs[0];
424        coords[1] = (double) idxs[1];
425        if ((coords[0] == -1) || (coords[1] == -1)) return null;
426        if (xStride > 0) {
427          coords[2] = xStride;
428        }
429
430        coords = subset.get("VertDim");
431        if (yStride > 0) {
432          coords[2] = yStride;
433        }
434        return subset;
435      }
436
437      public Map<String, double[]> getSubsetFromLonLatRect(double minLat, double maxLat, double minLon, double maxLon) {
438        return getSubsetFromLonLatRect(getDefaultSubset(), minLat, maxLat, minLon, maxLon);
439      }
440
441      public Map<String, double[]> getSubsetFromLonLatRect(double minLat, double maxLat, double minLon, double maxLon,
442                                             int xStride, int yStride, int zStride) {
443
444        return getSubsetFromLonLatRect(getDefaultSubset(), minLat, maxLat, minLon, maxLon, xStride, yStride, zStride);
445      }
446
447      public abstract float[] getVertBinAltitude() throws Exception;
448      public abstract float[] getTrackTimes() throws Exception;
449      public abstract RealType makeVertLocType() throws Exception;
450      public abstract RealType makeTrackTimeType() throws Exception;
451      public abstract float[] getTrackLongitude() throws Exception;
452      public abstract float[] getTrackLatitude() throws Exception;
453      
454      public static FieldImpl medianFilter(FieldImpl field, int window_lenx, int window_leny) throws VisADException, RemoteException, CloneNotSupportedException  {
455         Set dSet = field.getDomainSet();
456         if (dSet.getManifoldDimension() != 1) {
457            throw new VisADException("medianFilter: outer field domain must have manifoldDimension = 1");
458         }
459         int outerLen = dSet.getLength();
460         
461         FieldImpl filtField = (FieldImpl)field.clone();
462         
463         for (int t=0; t<outerLen; t++) {
464            FlatField ff = (FlatField) filtField.getSample(t, false);
465            medianFilter(ff, window_lenx, window_leny);
466         }
467         
468         return filtField;
469      }
470      
471      public static FlatField medianFilter(FlatField fltFld, int window_lenx, int window_leny) throws VisADException, RemoteException {
472         GriddedSet domSet = (GriddedSet) fltFld.getDomainSet();
473         FlatField filtFld = new FlatField((FunctionType)fltFld.getType(), domSet);
474         
475         int[] lens = domSet.getLengths();
476         int manifoldDimension = domSet.getManifoldDimension();
477         
478         float[][] rngVals = fltFld.getFloats(false);
479         int rngTupleDim = rngVals.length;
480         float[][] filtVals = new float[rngTupleDim][];
481         
482         if (manifoldDimension == 2) {
483            for (int t=0; t<rngTupleDim; t++) {
484               filtVals[t] = medianFilter(rngVals[t], lens[0], lens[1], window_lenx, window_leny);
485            }
486         }
487         else if (manifoldDimension == 3) {
488            int outerDimLen = lens[0];
489            filtVals = new float[rngTupleDim][lens[0]*lens[1]*lens[2]];
490            float[] rngVals2D = new float[lens[1]*lens[2]];
491            
492            for (int k = 0; k<outerDimLen; k++) {
493               int idx = k*lens[1]*lens[2];
494               for (int t=0; t<rngTupleDim; t++) {
495                  System.arraycopy(rngVals[t], idx, rngVals2D, 0, lens[1]*lens[2]);
496                  
497                  float[] fltA = medianFilter(rngVals2D, lens[1], lens[2], window_lenx, window_leny);
498                  
499                  System.arraycopy(fltA, 0, filtVals[t], idx, lens[1]*lens[2]);
500               }
501            }
502         }
503         
504         filtFld.setSamples(filtVals, false);
505         
506         return filtFld;
507      }
508      
509      public static float[] medianFilter(float[] A, int lenx, int leny, int window_lenx, int window_leny)
510           throws VisADException {
511        float[] result =  new float[A.length];
512        float[] window =  new float[window_lenx*window_leny];
513        float[] sortedWindow =  new float[window_lenx*window_leny];
514        int[] sort_indexes = new int[window_lenx*window_leny];
515        int[] indexes = new int[window_lenx*window_leny];
516        int[] indexesB = new int[window_lenx*window_leny];
517        
518        int[] numToInsertAt = new int[window_lenx*window_leny];
519        float[][] valsToInsert = new float[window_lenx*window_leny][window_lenx*window_leny];
520        int[][] idxsToInsert = new int[window_lenx*window_leny][window_lenx*window_leny];
521        
522        int[] numBefore = new int[window_lenx*window_leny];
523        float[][] valsBefore = new float[window_lenx*window_leny][window_lenx*window_leny];
524        int[][] idxsBefore = new int[window_lenx*window_leny][window_lenx*window_leny];
525        
526        int[] numAfter = new int[window_lenx*window_leny];
527        float[][] valsAfter = new float[window_lenx*window_leny][window_lenx*window_leny];
528        int[][] idxsAfter = new int[window_lenx*window_leny][window_lenx*window_leny];                
529        
530        float[] sortedArray = new float[window_lenx*window_leny];
531                                                                                                                                    
532        int a_idx;
533        int w_idx;
534                                                                                                                                    
535        int w_lenx = window_lenx/2;
536        int w_leny = window_leny/2;
537                                                                                                                                    
538        int lo;
539        int hi;
540        int ww_jj;
541        int ww_ii;
542        int cnt=0;
543        int ncnt;
544        int midx;
545        float median;
546        
547        int lenA = A.length;
548        
549        for (int i=0; i<lenx; i++) { // zig-zag better? Maybe, but more complicated
550          for (int j=0; j<leny; j++) {             
551            a_idx = j*lenx + i;
552            
553            if (j > 0) {
554              ncnt = 0;
555              for (int t=0; t<cnt; t++) {
556                 // last window index in data coords: A[lenx,leny]
557                 int k = indexes[sort_indexes[t]];
558                 ww_jj = k/lenx;
559                 ww_ii = k % lenx;
560                 
561                 // current window bounds in data coords
562                 int ww_jj_lo = j - w_leny;
563                 int ww_jj_hi = j + w_leny;
564                 int ww_ii_lo = i - w_lenx;
565                 int ww_ii_hi = i + w_lenx;
566                 
567                 if (ww_jj_lo < 0) ww_jj_lo = 0;
568                 if (ww_ii_lo < 0) ww_ii_lo = 0;
569                 if (ww_jj_hi > leny-1) ww_jj_hi = leny-1;
570                 if (ww_ii_hi > lenx-1) ww_ii_hi = lenx-1;
571                 
572                 
573                 // These are the points which overlap between the last and current window
574                 if ((ww_jj >= ww_jj_lo && ww_jj < ww_jj_hi) && (ww_ii >= ww_ii_lo && ww_ii < ww_ii_hi)) {
575                    window[ncnt] = sortedWindow[t];
576                    indexesB[ncnt] = k;
577                    ncnt++;
578                 }
579              }
580              
581              
582              // Add the new points from sliding the window to the overlap points above
583              java.util.Arrays.fill(numToInsertAt, 0);
584              java.util.Arrays.fill(numBefore, 0);
585              java.util.Arrays.fill(numAfter, 0);
586              
587              ww_jj = w_leny-1 + j;
588              for (int w_i=-w_lenx; w_i<w_lenx; w_i++) {
589                 ww_ii = w_i + i;
590                 int k = ww_jj*lenx+ww_ii;
591                  if (k >= 0 && k < lenA) {
592                     float val = A[k];
593                        if (ncnt > 0) {
594                           int t = 0;
595                           if (val < window[t]) {
596                                 valsBefore[0][numBefore[0]] = val;
597                                 idxsBefore[0][numBefore[0]] = k;
598                                 numBefore[0]++;  
599                                 continue;
600                           }                     
601                           t = ncnt-1;
602                           if (val >= window[t]) {
603                                 valsAfter[0][numAfter[0]] = val;
604                                 idxsAfter[0][numAfter[0]] = k;
605                                 numAfter[0]++;  
606                                 continue;
607                           }
608
609                           for (t=0; t<ncnt-1; t++) {
610                              if (val >= window[t] && val < window[t+1]) {
611                                 valsToInsert[t][numToInsertAt[t]] = val;
612                                 idxsToInsert[t][numToInsertAt[t]] = k;
613                                 numToInsertAt[t]++;
614                                 break;
615                              }
616                           } 
617                        }
618                        else if (!Float.isNaN(val)) {
619                                 valsBefore[0][numBefore[0]] = val;
620                                 idxsBefore[0][numBefore[0]] = k;
621                                 numBefore[0]++;  
622                                 continue;                           
623                        }
624                  }
625              }
626
627              // insert new unsorted values into the already sorted overlap window region
628              int tcnt = 0;
629              
630              for (int it=0; it<numBefore[0]; it++) {
631                 sortedArray[tcnt] = valsBefore[0][it];
632                 indexes[tcnt] = idxsBefore[0][it];
633                 tcnt++;
634              }  
635                       
636              for (int t=0; t<ncnt; t++) {
637                 sortedArray[tcnt] = window[t];
638                 indexes[tcnt] = indexesB[t];
639                 tcnt++;
640                 if (numToInsertAt[t] > 0) {
641                    if (numToInsertAt[t] == 2) { // two item sort here to save work for QuickSort
642                       float val0 = valsToInsert[t][0];
643                       float val1 = valsToInsert[t][1];
644                       
645                       if (val0 <= val1) {
646                          sortedArray[tcnt] = val0;
647                          indexes[tcnt] = idxsToInsert[t][0];
648                          tcnt++;
649                          
650                          sortedArray[tcnt] = val1;
651                          indexes[tcnt] = idxsToInsert[t][1];
652                          tcnt++;
653                       }
654                       else {
655                          sortedArray[tcnt] = val1;
656                          indexes[tcnt] = idxsToInsert[t][1];  
657                          tcnt++;
658                          
659                          sortedArray[tcnt] = val0;
660                          indexes[tcnt] = idxsToInsert[t][0];      
661                          tcnt++;
662                       }
663                    }
664                    else if (numToInsertAt[t] == 1) {
665                       sortedArray[tcnt] = valsToInsert[t][0];
666                       indexes[tcnt] = idxsToInsert[t][0];
667                       tcnt++;
668                    }
669                    else {
670                       for (int it=0; it<numToInsertAt[t]; it++) {
671                          sortedArray[tcnt] = valsToInsert[t][it];
672                          indexes[tcnt] = idxsToInsert[t][it];
673                          tcnt++;
674                       }
675                    }
676                 }
677              }
678              
679              for (int it=0; it<numAfter[0]; it++) {
680                 sortedArray[tcnt] = valsAfter[0][it];
681                 indexes[tcnt] = idxsAfter[0][it];
682                 tcnt++;
683              }  
684              
685              // Now sort the new unsorted and overlap sorted points together to get the new window median
686              
687              System.arraycopy(sortedArray, 0, sortedWindow, 0, tcnt);
688              if (tcnt > 0) {
689                 sort_indexes = QuickSort.sort(sortedWindow, 0, tcnt-1);
690                 median = sortedWindow[tcnt/2];
691              }
692              else {
693                 median = Float.NaN;
694              }
695              cnt = tcnt;
696
697            }
698            else { // full sort done once for each row (see note on zigzag above)
699            
700               cnt = 0;
701               for (int w_j=-w_leny; w_j<w_leny; w_j++) {
702                 for (int w_i=-w_lenx; w_i<w_lenx; w_i++) {
703                   ww_jj = w_j + j;
704                   ww_ii = w_i + i;
705                   w_idx = (w_j+w_leny)*window_lenx + (w_i+w_lenx);
706                   if ((ww_jj >= 0) && (ww_ii >=0) && (ww_jj < leny) && (ww_ii < lenx)) {
707                     int k = ww_jj*lenx+ww_ii;
708                     float val = A[k];
709                     if (!Float.isNaN(val)) {
710                       window[cnt] = val;
711                       indexes[cnt] = k;
712                       cnt++;
713                     }
714                   }
715                 }
716               }
717            
718            
719               System.arraycopy(window, 0, sortedWindow, 0, cnt);
720               if (cnt > 0) {
721                  sort_indexes = QuickSort.sort(sortedWindow, 0, cnt-1);
722                  midx = cnt/2;
723                  median = sortedWindow[midx];
724               }
725               else {
726                  median = Float.NaN;
727               }
728               
729            }
730            
731            if (Float.isNaN(A[a_idx])) {
732              result[a_idx] = Float.NaN;
733            }
734            else {
735              result[a_idx] = median;
736            }
737            
738          }
739        }
740        
741        return result;
742      }
743      
744      public static float[] medianFilterOrg(float[] A, int lenx, int leny, int window_lenx, int window_leny)
745           throws VisADException {
746        float[] result =  new float[A.length];
747        float[] window =  new float[window_lenx*window_leny];
748        float[] new_window =  new float[window_lenx*window_leny];
749        int[] sort_indexes = new int[window_lenx*window_leny];
750                                                                                                                                    
751        int a_idx;
752        int w_idx;
753                                                                                                                                    
754        int w_lenx = window_lenx/2;
755        int w_leny = window_leny/2;
756                                                                                                                                    
757        int lo;
758        int hi;
759        int ww_jj;
760        int ww_ii;
761        int cnt;
762                                                                                                                                    
763        for (int j=0; j<leny; j++) {
764          for (int i=0; i<lenx; i++) {
765            a_idx = j*lenx + i;
766            cnt = 0;
767            for (int w_j=-w_leny; w_j<w_leny; w_j++) {
768              for (int w_i=-w_lenx; w_i<w_lenx; w_i++) {
769                ww_jj = w_j + j;
770                ww_ii = w_i + i;
771                w_idx = (w_j+w_leny)*window_lenx + (w_i+w_lenx);
772                if ((ww_jj >= 0) && (ww_ii >=0) && (ww_jj < leny) && (ww_ii < lenx)) {
773                  window[cnt] = A[ww_jj*lenx+ww_ii];
774                  cnt++;
775                }
776              }
777            }
778            System.arraycopy(window, 0, new_window, 0, cnt);
779            //-sort_indexes = QuickSort.sort(new_window, sort_indexes);
780            sort_indexes = QuickSort.sort(new_window);
781            result[a_idx] = new_window[cnt/2];
782          }
783        }
784        return result;
785      }
786}