001 /*
002 * This file is part of McIDAS-V
003 *
004 * Copyright 2007-2013
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 package edu.wisc.ssec.mcidasv.data.hydra;
029
030 import visad.*;
031
032 public 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 }