001/*
002 * This file is part of McIDAS-V
003 *
004 * Copyright 2007-2016
005 * Space Science and Engineering Center (SSEC)
006 * University of Wisconsin - Madison
007 * 1225 W. Dayton Street, Madison, WI 53706, USA
008 * https://www.ssec.wisc.edu/mcidas
009 * 
010 * All Rights Reserved
011 * 
012 * McIDAS-V is built on Unidata's IDV and SSEC's VisAD libraries, and
013 * some McIDAS-V source code is based on IDV and VisAD source code.  
014 * 
015 * McIDAS-V is free software; you can redistribute it and/or modify
016 * it under the terms of the GNU Lesser Public License as published by
017 * the Free Software Foundation; either version 3 of the License, or
018 * (at your option) any later version.
019 * 
020 * McIDAS-V is distributed in the hope that it will be useful,
021 * but WITHOUT ANY WARRANTY; without even the implied warranty of
022 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
023 * GNU Lesser Public License for more details.
024 * 
025 * You should have received a copy of the GNU Lesser Public License
026 * along with this program.  If not, see http://www.gnu.org/licenses.
027 */
028package edu.wisc.ssec.mcidasv.data.hydra;
029
030import visad.*;
031
032public class HistogramField {
033
034    Linear2DSet histSet;
035    Linear1DSet set0;
036    Linear1DSet set1;
037    int len0;
038    int len1;
039    int[] count;
040    int[][] indexes;
041    FlatField field_0;
042    FlatField field_1;
043    FlatField mask_field;
044    float[][] maskRange;
045    Class rangeType;
046    byte[][] mask = new byte[3][];
047    byte[] order = new byte[3];
048
049    public FlatField scatterDensityField;
050
051    public HistogramField(FlatField field_0, FlatField field_1,
052            FlatField mask_field,
053            int n_bins, int bin_size)
054            throws Exception {
055        this.field_0 = field_0;
056        this.field_1 = field_1;
057        this.mask_field = mask_field;
058        maskRange = mask_field.getFloats(false);
059
060        mask[0] = new byte[maskRange[0].length];
061        mask[1] = new byte[maskRange[0].length];
062        mask[2] = new byte[maskRange[0].length];
063        java.util.Arrays.fill(mask[0], Byte.MAX_VALUE);
064        java.util.Arrays.fill(mask[1], Byte.MAX_VALUE);
065        java.util.Arrays.fill(mask[2], Byte.MAX_VALUE);
066        java.util.Arrays.fill(order, Byte.MAX_VALUE);
067
068        Set[] rangeSets = field_0.getRangeSets();
069        Set rngSet = rangeSets[0];
070
071        if (rngSet instanceof FloatSet) {
072            rangeType = Float.TYPE;
073        } else if (rngSet instanceof DoubleSet) {
074            rangeType = Double.TYPE;
075        } else if (rngSet instanceof IntegerSet) {
076            rangeType = Integer.TYPE;
077        }
078
079        double[] minmax_0 = {Double.MAX_VALUE, -Double.MAX_VALUE};
080        double[] minmax_1 = {Double.MAX_VALUE, -Double.MAX_VALUE};
081
082
083        if (rangeType == Integer.TYPE) {
084          //Ghansham: Dont do any allocation here. Do based on the individual ranges of fieldX and fieldY respectively
085        } else {
086            indexes = new int[n_bins * n_bins][1];
087            count = new int[n_bins * n_bins];
088        }
089
090
091        double[][] val = new double[2][1];
092        int[] histIdx = null;
093
094        if (rangeType == Double.TYPE) {
095            double[][] vals_0 = field_0.getValues(false);
096            double[][] vals_1 = field_1.getValues(false);
097            int n_samples = vals_0[0].length;
098            for (int k = 0; k < n_samples; k++) {
099                double v0 = vals_0[0][k];
100                if (v0 < minmax_0[0]) {
101                    minmax_0[0] = v0;
102                } else if (v0 > minmax_0[1]) {
103                    minmax_0[1] = v0;
104                }
105
106                double v1 = vals_1[0][k];
107                if (v1 < minmax_1[0]) {
108                    minmax_1[0] = v1;
109                } else if (v1 > minmax_1[1]) {
110                    minmax_1[1] = v1;
111                }
112            }
113
114            histSet = new Linear2DSet(minmax_0[0], minmax_0[1], n_bins,
115                    minmax_1[0], minmax_1[1], n_bins);
116
117            for (int k = 0; k < n_samples; k++) {
118                val[0][0] = vals_0[0][k];
119                val[1][0] = vals_1[0][k];
120                histIdx = histSet.doubleToIndex(val);
121                if (histIdx[0] >= 0) {
122                    int len = indexes[histIdx[0]].length;
123                    if (count[histIdx[0]] > len - 1) { //-grow array
124                        int[] tmp = new int[len + bin_size];
125                        System.arraycopy(indexes[histIdx[0]], 0, tmp, 0, len);
126                        indexes[histIdx[0]] = tmp;
127                    }
128                    indexes[histIdx[0]][count[histIdx[0]]++] = k;
129                }
130            }
131        } else if (rangeType == Float.TYPE) {
132            float[][] vals_0 = field_0.getFloats(false);
133            float[][] vals_1 = field_1.getFloats(false);
134            int n_samples = vals_0[0].length;
135            for (int k = 0; k < n_samples; k++) {
136                double v0 = vals_0[0][k];
137                if (v0 < minmax_0[0]) {
138                    minmax_0[0] = v0;
139                } else if (v0 > minmax_0[1]) {
140                    minmax_0[1] = v0;
141                }
142                double v1 = vals_1[0][k];
143                if (v1 < minmax_1[0]) {
144                    minmax_1[0] = v1;
145                } else if (v1 > minmax_1[1]) {
146                    minmax_1[1] = v1;
147                }
148
149            }
150            histSet = new Linear2DSet(minmax_0[0], minmax_0[1], n_bins, minmax_1[0], minmax_1[1], n_bins);
151            for (int k = 0; k < n_samples; k++) {
152                val[0][0] = vals_0[0][k];
153                val[1][0] = vals_1[0][k];
154                histIdx = histSet.doubleToIndex(val);
155                if (histIdx[0] >= 0) {
156                    int len = indexes[histIdx[0]].length;
157                    if (count[histIdx[0]] > len - 1) { //-grow array
158
159                        int[] tmp = new int[len + bin_size];
160                        System.arraycopy(indexes[histIdx[0]], 0, tmp, 0, len);
161                        indexes[histIdx[0]] = tmp;
162                    }
163                    indexes[histIdx[0]][count[histIdx[0]]++] = k;
164                }
165            }
166        } else if (rangeType == Integer.TYPE) {
167            float[][] vals_0 = field_0.getFloats(false);
168            float[][] vals_1 = field_1.getFloats(false);
169            int n_samples = vals_0[0].length;
170            for (int k = 0; k < n_samples; k++) {
171                double v0 = vals_0[0][k];
172                if (v0 < minmax_0[0]) {
173                    minmax_0[0] = v0;
174                } else if (v0 > minmax_0[1]) {
175                    minmax_0[1] = v0;
176                }
177                double v1 = vals_1[0][k];
178                if (v1 < minmax_1[0]) {
179                    minmax_1[0] = v1;
180                } else if (v1 > minmax_1[1]) {
181                    minmax_1[1] = v1;
182                }
183            }
184
185
186            int startX = (int) minmax_0[0];
187            int endX = (int) minmax_0[1];
188            int startY = (int) minmax_1[0];
189            int endY = (int) minmax_1[1];
190            int lenX = endX - startX + 1;
191            int lenY = endY - startY + 1;
192            histSet = new Linear2DSet(startX, endX, lenX, startY, endY, lenY);
193
194
195
196            //Ghansham:Allocate here based on length of lenghts of XField and YField
197            indexes = new int[lenY * lenX][]; //Ghansham:Dont allocate here if not required.
198            count = new int[lenY * lenX];
199
200            //First calculate frequency of each grey count.
201            for (int k = 0; k < n_samples; k++) {
202                val[0][0] = vals_0[0][k];
203                val[1][0] = vals_1[0][k];
204                histIdx = histSet.doubleToIndex(val);
205                if (histIdx[0] >= 0) {
206                    count[histIdx[0]]++;
207                }
208            }
209
210            for (int k = 0; k < n_samples; k++) {
211                val[0][0] = vals_0[0][k];
212                val[1][0] = vals_1[0][k];
213                histIdx = histSet.doubleToIndex(val);
214                if (histIdx[0] >= 0) {
215                    if (indexes[histIdx[0]] == null) { //Tricky stuff is here: encountering a particular grey count first time.
216                        indexes[histIdx[0]] = new int[count[histIdx[0]]]; //Allocating the values based on the frequency of this grey count (calculate earlier). No extra allocation at all
217                        count[histIdx[0]] = 0;  //Resetting the frequency to 0.
218                    }
219                    indexes[histIdx[0]][count[histIdx[0]]++] = k;
220                }
221            }
222        }
223
224
225        set0 = histSet.getLinear1DComponent(0);
226        set1 = histSet.getLinear1DComponent(1);
227        len0 = set0.getLength();
228        len1 = set1.getLength();
229
230
231        Linear2DSet dSet = (Linear2DSet) histSet.changeMathType(new RealTupleType(RealType.XAxis, RealType.YAxis));
232        scatterDensityField = new FlatField(
233            new FunctionType(((SetType)dSet.getType()).getDomain(), RealType.getRealType("ScatterDensity")), dSet);
234        float[][] fltCount = new float[1][count.length];
235        for (int i=0; i<count.length; i++) { 
236            fltCount[0][i] = (float) count[i];
237            if (count[i] == 0) {
238               fltCount[0][i] = Float.NaN;
239            }
240            else {
241               fltCount[0][i] = (float) java.lang.Math.log((double)fltCount[0][i]);
242            }
243        }
244        scatterDensityField.setSamples(fltCount);
245    }
246
247    public FlatField getScatterDensityField() {
248        return scatterDensityField;
249    }
250
251    public void markMaskFieldByRange(double[] lowhi_0, double[] lowhi_1, float maskVal)
252            throws Exception {
253        reorder((byte)maskVal);
254
255        int[] hist0 = set0.doubleToIndex(new double[][] {{lowhi_0[0], lowhi_0[1]}});
256        int[] hist1 = set1.doubleToIndex(new double[][] {{lowhi_1[0], lowhi_1[1]}});
257
258        if (hist0[0] < 0) {
259            if (lowhi_0[0] < lowhi_0[1]) {
260                hist0[0] = 0;
261            } else {
262                hist0[0] = len0 - 1;
263            }
264        }
265        if (hist0[1] < 0) {
266            if (lowhi_0[0] < lowhi_0[1]) {
267                hist0[1] = len0 - 1;
268            } else {
269                hist0[1] = 0;
270            }
271        }
272
273        if (hist1[0] < 0) {
274            if (lowhi_1[0] < lowhi_1[1]) {
275                hist1[0] = 0;
276            } else {
277                hist1[0] = len1 - 1;
278            }
279        }
280        if (hist1[1] < 0) {
281            if (lowhi_1[0] < lowhi_1[1]) {
282                hist1[1] = len1 - 1;
283            } else {
284                hist1[1] = 0;
285            }
286        }
287
288        int h00, h01, h10, h11;
289
290
291        h10 = hist1[1];
292        h11 = hist1[0];
293        if (hist1[0] < hist1[1]) {
294            h10 = hist1[0];
295            h11 = hist1[1];
296        }
297
298        h00 = hist0[1];
299        h01 = hist0[0];
300        if (hist0[0] < hist0[1]) {
301            h00 = hist0[0];
302            h01 = hist0[1];
303        }
304
305        for (int k = 0; k < maskRange[0].length; k++) {
306            if (maskRange[0][k] == maskVal) {
307                maskRange[0][k] = Float.NaN;
308                mask[(byte)maskVal][k] = Byte.MAX_VALUE;
309            }
310        }
311
312        int lenX = set0.getLengthX();
313
314        for (int j = h10; j <= h11; j++) {
315            int col_factor = j * lenX;
316            for (int i = h00; i <= h01; i++) {
317                int idx = col_factor + i;
318                for (int k = 0; k < count[idx]; k++) {
319                    maskRange[0][indexes[idx][k]] = maskVal;
320                    mask[(byte)maskVal][indexes[idx][k]] = (byte)maskVal;
321                }
322            }
323        }
324
325        mask_field.setSamples(maskRange, false);
326    }
327
328    public void markMaskFieldByCurve(float[][] curve, float maskVal) throws Exception {
329        reorder((byte)maskVal);
330        float[][] samples0 = set0.getSamples();
331        float[][] samples1 = set1.getSamples();
332
333        boolean[][] checked = null;
334        boolean[][] inside = null;
335
336
337        if (rangeType == Integer.TYPE) { //Dealing with rangeSet constructed fields separately.
338            float[][] vals_0 = field_0.getFloats(false);
339            float[][] vals_1 = field_1.getFloats(false);
340            int lenX = set0.getLength();
341            int lenY = set1.getLength();
342            checked = new boolean[lenX][lenY];
343            inside = new boolean[lenX][lenY];
344            for (int jj = 0; jj < lenX; jj++) {
345                java.util.Arrays.fill(checked[jj], false);
346                java.util.Arrays.fill(inside[jj], false);
347            }
348            for (int jj = 0; jj < lenY - 1; jj++) {
349                for (int ii = 0; ii < lenX - 1; ii++) {
350                    int idx = jj * lenX + ii; //Calclualting the index value in the start only.
351                    if (count[idx] > 0) { //No need to do go further if the frequency of particular grey count is zero.
352                        int inside_cnt = 0;
353                        if (!checked[ii][jj]) {
354                            float x = samples0[0][ii];
355                            float y = samples1[0][jj];
356                            if (DelaunayCustom.inside(curve, x, y)) {
357                                inside_cnt++;
358                                inside[ii][jj] = true;
359                            }
360                            checked[ii][jj] = true;
361                        } else if (inside[ii][jj]) {
362                            inside_cnt++;
363                        }
364
365                        if (!checked[ii + 1][jj]) {
366                            float x = samples0[0][ii + 1];
367                            float y = samples1[0][jj];
368                            if (DelaunayCustom.inside(curve, x, y)) {
369                                inside_cnt++;
370                                inside[ii + 1][jj] = true;
371                            }
372                            checked[ii + 1][jj] = true;
373                        } else if (inside[ii + 1][jj]) {
374                            inside_cnt++;
375                        }
376
377                        if (!checked[ii][jj + 1]) {
378                            float x = samples0[0][ii];
379                            float y = samples1[0][jj + 1];
380                            if (DelaunayCustom.inside(curve, x, y)) {
381                                inside_cnt++;
382                                inside[ii][jj + 1] = true;
383                            }
384                            checked[ii][jj + 1] = true;
385                        } else if (inside[ii][jj + 1]) {
386                            inside_cnt++;
387                        }
388
389                        if (!checked[ii + 1][jj + 1]) {
390                            float x = samples0[0][ii + 1];
391                            float y = samples1[0][jj + 1];
392                            if (DelaunayCustom.inside(curve, x, y)) {
393                                inside_cnt++;
394                                inside[ii + 1][jj + 1] = true;
395                            }
396                            checked[ii + 1][jj + 1] = true;
397                        } else if (inside[ii + 1][jj + 1]) {
398                            inside_cnt++;
399                        }
400
401                        if (inside_cnt == 0) {
402                            continue;
403                        } else if (inside_cnt == 4) {
404                            for (int k = 0; k < count[idx]; k++) {
405                                maskRange[0][indexes[idx][k]] = maskVal;
406                            }
407                        } else if (inside_cnt > 0 && inside_cnt < 4) {
408                            for (int k = 0; k < count[idx]; k++) {
409                                float xx = vals_0[0][indexes[idx][k]];
410                                float yy = vals_1[0][indexes[idx][k]];
411                                if (DelaunayCustom.inside(curve, xx, yy)) {
412                                    maskRange[0][indexes[idx][k]] = maskVal;
413                                }
414                            }
415                        }
416                    }
417                }
418            }
419        } else {
420            int len = set0.getLength();
421            checked = new boolean[len][len];
422            inside = new boolean[len][len];
423            for (int jj = 0; jj < len; jj++) {
424                java.util.Arrays.fill(checked[jj], false);
425                java.util.Arrays.fill(inside[jj], false);
426            }
427            for (int jj = 0; jj < len - 1; jj++) {
428                for (int ii = 0; ii < len - 1; ii++) {
429                    int idx = jj * set0.getLengthX() + ii; //Calclualting the index value in the start only.
430                    if (count[idx] > 0) { //No need to do go further if the frequency of particular value is zero.
431                        int inside_cnt = 0;
432                        if (!checked[ii][jj]) {
433                            float x = samples0[0][ii];
434                            float y = samples1[0][jj];
435                            if (DelaunayCustom.inside(curve, x, y)) {
436                                inside_cnt++;
437                                inside[ii][jj] = true;
438                            }
439                            checked[ii][jj] = true;
440                        } else if (inside[ii][jj]) {
441                            inside_cnt++;
442                        }
443
444                        if (!checked[ii + 1][jj]) {
445                            float x = samples0[0][ii + 1];
446                            float y = samples1[0][jj];
447                            if (DelaunayCustom.inside(curve, x, y)) {
448                                inside_cnt++;
449                                inside[ii + 1][jj] = true;
450                            }
451                            checked[ii + 1][jj] = true;
452                        } else if (inside[ii + 1][jj]) {
453                            inside_cnt++;
454                        }
455
456                        if (!checked[ii][jj + 1]) {
457                            float x = samples0[0][ii];
458                            float y = samples1[0][jj + 1];
459                            if (DelaunayCustom.inside(curve, x, y)) {
460                                inside_cnt++;
461                                inside[ii][jj + 1] = true;
462                            }
463                            checked[ii][jj + 1] = true;
464                        } else if (inside[ii][jj + 1]) {
465                            inside_cnt++;
466                        }
467
468                        if (!checked[ii + 1][jj + 1]) {
469                            float x = samples0[0][ii + 1];
470                            float y = samples1[0][jj + 1];
471                            if (DelaunayCustom.inside(curve, x, y)) {
472                                inside_cnt++;
473                                inside[ii + 1][jj + 1] = true;
474                            }
475                            checked[ii + 1][jj + 1] = true;
476                        } else if (inside[ii + 1][jj + 1]) {
477                            inside_cnt++;
478                        }
479                        if (inside_cnt == 0) {
480                            continue;
481                        }
482
483                        if (rangeType == Float.TYPE) {
484                            float[][] vals_0 = field_0.getFloats(false);
485                            float[][] vals_1 = field_1.getFloats(false);
486                            if (inside_cnt == 4) {
487
488                                for (int k = 0; k < count[idx]; k++) {
489                                    maskRange[0][indexes[idx][k]] = maskVal;
490                                    mask[(byte)maskVal][indexes[idx][k]] = (byte)maskVal;
491                                }
492                            }
493                            if (inside_cnt > 0 && inside_cnt < 4) {
494
495                                for (int k = 0; k < count[idx]; k++) {
496                                    float xx = vals_0[0][indexes[idx][k]];
497                                    float yy = vals_1[0][indexes[idx][k]];
498                                    if (DelaunayCustom.inside(curve, xx, yy)) {
499                                        maskRange[0][indexes[idx][k]] = maskVal;
500                                        mask[(byte)maskVal][indexes[idx][k]] = (byte)maskVal;
501                                    }
502                                }
503                            }
504                        } else if (rangeType == Double.TYPE) {
505                            double[][] vals_0 = field_0.getValues(false);
506                            double[][] vals_1 = field_1.getValues(false);
507                            if (inside_cnt == 4) {
508                                for (int k = 0; k < count[idx]; k++) {
509                                    maskRange[0][indexes[idx][k]] = maskVal;
510                                    mask[(byte)maskVal][indexes[idx][k]] = (byte)maskVal;
511                                }
512                            }
513                            if (inside_cnt > 0 && inside_cnt < 4) {
514
515                                for (int k = 0; k < count[idx]; k++) {
516                                    double xx = vals_0[0][indexes[idx][k]];
517                                    double yy = vals_1[0][indexes[idx][k]];
518                                    if (DelaunayCustom.inside(curve, (float) xx, (float) yy)) {
519                                        maskRange[0][indexes[idx][k]] = maskVal;
520                                        mask[(byte)maskVal][indexes[idx][k]] = (byte)maskVal;
521                                    }
522                                }
523                            }
524                        }
525                    }
526                }
527            }
528        }
529
530        mask_field.setSamples(maskRange, false);
531    }
532
533    private void reorder(byte maskVal) {
534       order[2] = order[1];
535       order[1] = order[0];
536       order[0] = maskVal;
537    }
538
539    public void clearMaskField(float maskVal) {
540        for (int k = 0; k < maskRange[0].length; k++) {
541            maskRange[0][k] = Float.NaN;
542            mask[(byte)maskVal][k] = Byte.MAX_VALUE;
543        }
544
545        for (int t=0; t<order.length; t++) {
546           if (order[t] == (byte)maskVal) {
547               order[t] = Byte.MAX_VALUE;
548           }
549        }
550
551        for (int t=order.length-1; t >=0; t--) {
552            if (order[t] != Byte.MAX_VALUE) {
553               for (int k=0; k<maskRange[0].length; k++) {
554                   if (mask[order[t]][k] != Byte.MAX_VALUE) {
555                      maskRange[0][k] = (float) order[t];
556                   }
557               }
558            }
559        }
560    }
561
562    public void resetMaskField(float maskVal) throws Exception {
563        clearMaskField(maskVal);
564        mask_field.setSamples(maskRange, false);
565    }
566}