Source code

001/*
002 * This file is part of McIDAS-V
003 *
004 * Copyright 2007-2018
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      /**
472       * Apply a median filter to FlatField range If domain dimension == 3, range is filtered as succession of
473       * 2D filtered layers along the slowest varying dimension.
474       * @param fltFld incoming FlatField to be filtered
475       * @param window_lenx filter window (kernel) dimensions; x: fastest varying dimensions
476       * @param window_leny
477       * @return FlatField with filtered (copied) range
478       * @throws VisADException
479       * @throws RemoteException 
480       */
481      public static FlatField medianFilter(FlatField fltFld, int window_lenx, int window_leny) throws VisADException, RemoteException {
482         GriddedSet domSet = (GriddedSet) fltFld.getDomainSet();
483         FlatField filtFld = new FlatField((FunctionType)fltFld.getType(), domSet);
484         
485         int[] lens = domSet.getLengths();
486         int manifoldDimension = domSet.getManifoldDimension();
487         
488         float[][] rngVals = fltFld.getFloats(false);
489         int rngTupleDim = rngVals.length;
490         float[][] filtVals = new float[rngTupleDim][];
491         
492         if (manifoldDimension == 2) {
493            for (int t=0; t<rngTupleDim; t++) {
494               filtVals[t] = medianFilter(rngVals[t], lens[0], lens[1], window_lenx, window_leny);
495            }
496         }
497         else if (manifoldDimension == 3) {
498            int outerDimLen = lens[0];
499            filtVals = new float[rngTupleDim][lens[0]*lens[1]*lens[2]];
500            float[] rngVals2D = new float[lens[1]*lens[2]];
501            
502            for (int k = 0; k<outerDimLen; k++) {
503               int idx = k*lens[1]*lens[2];
504               for (int t=0; t<rngTupleDim; t++) {
505                  System.arraycopy(rngVals[t], idx, rngVals2D, 0, lens[1]*lens[2]);
506                  
507                  float[] fltA = medianFilter(rngVals2D, lens[1], lens[2], window_lenx, window_leny);
508                  
509                  System.arraycopy(fltA, 0, filtVals[t], idx, lens[1]*lens[2]);
510               }
511            }
512         }
513         
514         filtFld.setSamples(filtVals, false);
515         
516         return filtFld;
517      }
518      
519      /**
520       * Apply median filter to 2D array of values
521       * @param A 2D array to filter
522       * @param lenx Dimensions of A, x varies fastest
523       * @param leny
524       * @param window_lenx Dimensions of window (kernel) x fastest varying
525       * @param window_leny
526       * @return
527       * @throws VisADException 
528       */
529      public static float[] medianFilter(float[] A, int lenx, int leny, int window_lenx, int window_leny)
530           throws VisADException {
531        float[] result =  new float[A.length];
532        float[] window =  new float[window_lenx*window_leny];
533        float[] sortedWindow =  new float[window_lenx*window_leny];
534        int[] sort_indexes = new int[window_lenx*window_leny];
535        int[] indexes = new int[window_lenx*window_leny];
536        int[] indexesB = new int[window_lenx*window_leny];
537        
538        int[] numToInsertAt = new int[window_lenx*window_leny];
539        float[][] valsToInsert = new float[window_lenx*window_leny][window_lenx*window_leny];
540        int[][] idxsToInsert = new int[window_lenx*window_leny][window_lenx*window_leny];
541        
542        int[] numBefore = new int[window_lenx*window_leny];
543        float[][] valsBefore = new float[window_lenx*window_leny][window_lenx*window_leny];
544        int[][] idxsBefore = new int[window_lenx*window_leny][window_lenx*window_leny];
545        
546        int[] numAfter = new int[window_lenx*window_leny];
547        float[][] valsAfter = new float[window_lenx*window_leny][window_lenx*window_leny];
548        int[][] idxsAfter = new int[window_lenx*window_leny][window_lenx*window_leny];                
549        
550        float[] sortedArray = new float[window_lenx*window_leny];
551                                                                                                                                    
552        int a_idx;
553        int w_idx;
554                                                                                                                                    
555        int w_lenx = window_lenx/2;
556        int w_leny = window_leny/2;
557                                                                                                                                    
558        int lo;
559        int hi;
560        int ww_jj;
561        int ww_ii;
562        int cnt=0;
563        int ncnt;
564        int midx;
565        float median;
566        
567        int lenA = A.length;
568        
569        for (int i=0; i<lenx; i++) { // zig-zag better? Maybe, but more complicated
570          for (int j=0; j<leny; j++) {             
571            a_idx = j*lenx + i;
572            
573            if (j > 0) {
574              ncnt = 0;
575              for (int t=0; t<cnt; t++) {
576                 // last window index in data coords: A[lenx,leny]
577                 int k = indexes[sort_indexes[t]];
578                 ww_jj = k/lenx;
579                 ww_ii = k % lenx;
580                 
581                 // current window bounds in data coords
582                 int ww_jj_lo = j - w_leny;
583                 int ww_jj_hi = j + w_leny;
584                 int ww_ii_lo = i - w_lenx;
585                 int ww_ii_hi = i + w_lenx;
586                 
587                 if (ww_jj_lo < 0) ww_jj_lo = 0;
588                 if (ww_ii_lo < 0) ww_ii_lo = 0;
589                 if (ww_jj_hi > leny-1) ww_jj_hi = leny-1;
590                 if (ww_ii_hi > lenx-1) ww_ii_hi = lenx-1;
591                 
592                 
593                 // These are the points which overlap between the last and current window
594                 if ((ww_jj >= ww_jj_lo && ww_jj < ww_jj_hi) && (ww_ii >= ww_ii_lo && ww_ii < ww_ii_hi)) {
595                    window[ncnt] = sortedWindow[t];
596                    indexesB[ncnt] = k;
597                    ncnt++;
598                 }
599              }
600              
601              
602              // Add the new points from sliding the window to the overlap points above
603              java.util.Arrays.fill(numToInsertAt, 0);
604              java.util.Arrays.fill(numBefore, 0);
605              java.util.Arrays.fill(numAfter, 0);
606              
607              ww_jj = w_leny-1 + j;
608              for (int w_i=-w_lenx; w_i<w_lenx; w_i++) {
609                 ww_ii = w_i + i;
610                 int k = ww_jj*lenx+ww_ii;
611                  if (k >= 0 && k < lenA) {
612                     float val = A[k];
613                        if (ncnt > 0) {
614                           int t = 0;
615                           if (val < window[t]) {
616                                 valsBefore[0][numBefore[0]] = val;
617                                 idxsBefore[0][numBefore[0]] = k;
618                                 numBefore[0]++;  
619                                 continue;
620                           }                     
621                           t = ncnt-1;
622                           if (val >= window[t]) {
623                                 valsAfter[0][numAfter[0]] = val;
624                                 idxsAfter[0][numAfter[0]] = k;
625                                 numAfter[0]++;  
626                                 continue;
627                           }
628
629                           for (t=0; t<ncnt-1; t++) {
630                              if (val >= window[t] && val < window[t+1]) {
631                                 valsToInsert[t][numToInsertAt[t]] = val;
632                                 idxsToInsert[t][numToInsertAt[t]] = k;
633                                 numToInsertAt[t]++;
634                                 break;
635                              }
636                           } 
637                        }
638                        else if (!Float.isNaN(val)) {
639                                 valsBefore[0][numBefore[0]] = val;
640                                 idxsBefore[0][numBefore[0]] = k;
641                                 numBefore[0]++;  
642                                 continue;                           
643                        }
644                  }
645              }
646
647              // insert new unsorted values into the already sorted overlap window region
648              int tcnt = 0;
649              
650              for (int it=0; it<numBefore[0]; it++) {
651                 sortedArray[tcnt] = valsBefore[0][it];
652                 indexes[tcnt] = idxsBefore[0][it];
653                 tcnt++;
654              }  
655                       
656              for (int t=0; t<ncnt; t++) {
657                 sortedArray[tcnt] = window[t];
658                 indexes[tcnt] = indexesB[t];
659                 tcnt++;
660                 if (numToInsertAt[t] > 0) {
661                    if (numToInsertAt[t] == 2) { // two item sort here to save work for QuickSort
662                       float val0 = valsToInsert[t][0];
663                       float val1 = valsToInsert[t][1];
664                       
665                       if (val0 <= val1) {
666                          sortedArray[tcnt] = val0;
667                          indexes[tcnt] = idxsToInsert[t][0];
668                          tcnt++;
669                          
670                          sortedArray[tcnt] = val1;
671                          indexes[tcnt] = idxsToInsert[t][1];
672                          tcnt++;
673                       }
674                       else {
675                          sortedArray[tcnt] = val1;
676                          indexes[tcnt] = idxsToInsert[t][1];  
677                          tcnt++;
678                          
679                          sortedArray[tcnt] = val0;
680                          indexes[tcnt] = idxsToInsert[t][0];      
681                          tcnt++;
682                       }
683                    }
684                    else if (numToInsertAt[t] == 1) {
685                       sortedArray[tcnt] = valsToInsert[t][0];
686                       indexes[tcnt] = idxsToInsert[t][0];
687                       tcnt++;
688                    }
689                    else {
690                       for (int it=0; it<numToInsertAt[t]; it++) {
691                          sortedArray[tcnt] = valsToInsert[t][it];
692                          indexes[tcnt] = idxsToInsert[t][it];
693                          tcnt++;
694                       }
695                    }
696                 }
697              }
698              
699              for (int it=0; it<numAfter[0]; it++) {
700                 sortedArray[tcnt] = valsAfter[0][it];
701                 indexes[tcnt] = idxsAfter[0][it];
702                 tcnt++;
703              }  
704              
705              // Now sort the new unsorted and overlap sorted points together to get the new window median
706              
707              System.arraycopy(sortedArray, 0, sortedWindow, 0, tcnt);
708              if (tcnt > 0) {
709                 sort_indexes = QuickSort.sort(sortedWindow, 0, tcnt-1);
710                 median = sortedWindow[tcnt/2];
711              }
712              else {
713                 median = Float.NaN;
714              }
715              cnt = tcnt;
716
717            }
718            else { // full sort done once for each row (see note on zigzag above)
719            
720               cnt = 0;
721               for (int w_j=-w_leny; w_j<w_leny; w_j++) {
722                 for (int w_i=-w_lenx; w_i<w_lenx; w_i++) {
723                   ww_jj = w_j + j;
724                   ww_ii = w_i + i;
725                   w_idx = (w_j+w_leny)*window_lenx + (w_i+w_lenx);
726                   if ((ww_jj >= 0) && (ww_ii >=0) && (ww_jj < leny) && (ww_ii < lenx)) {
727                     int k = ww_jj*lenx+ww_ii;
728                     float val = A[k];
729                     if (!Float.isNaN(val)) {
730                       window[cnt] = val;
731                       indexes[cnt] = k;
732                       cnt++;
733                     }
734                   }
735                 }
736               }
737            
738            
739               System.arraycopy(window, 0, sortedWindow, 0, cnt);
740               if (cnt > 0) {
741                  sort_indexes = QuickSort.sort(sortedWindow, 0, cnt-1);
742                  midx = cnt/2;
743                  median = sortedWindow[midx];
744               }
745               else {
746                  median = Float.NaN;
747               }
748               
749            }
750            
751            if (Float.isNaN(A[a_idx])) {
752              result[a_idx] = Float.NaN;
753            }
754            else {
755              result[a_idx] = median;
756            }
757            
758          }
759        }
760        
761        return result;
762      }
763
764}