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