1 package nom.tam.fits.compression.algorithm.quant;
2
3 /*
4 * #%L
5 * nom.tam FITS library
6 * %%
7 * Copyright (C) 1996 - 2024 nom-tam-fits
8 * %%
9 * This is free and unencumbered software released into the public domain.
10 *
11 * Anyone is free to copy, modify, publish, use, compile, sell, or
12 * distribute this software, either in source code form or as a compiled
13 * binary, for any purpose, commercial or non-commercial, and by any
14 * means.
15 *
16 * In jurisdictions that recognize copyright laws, the author or authors
17 * of this software dedicate any and all copyright interest in the
18 * software to the public domain. We make this dedication for the benefit
19 * of the public at large and to the detriment of our heirs and
20 * successors. We intend this dedication to be an overt act of
21 * relinquishment in perpetuity of all present and future rights to this
22 * software under copyright law.
23 *
24 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
25 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
26 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
27 * IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
28 * OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
29 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
30 * OTHER DEALINGS IN THE SOFTWARE.
31 * #L%
32 */
33
34 import java.util.Arrays;
35
36 /**
37 * (<i>for internal use</i>) Determines the optimal quantization to use for floating-point data. It estimates the noise
38 * level in the data to determine qhat quantization should be use to lose no information above the noise level.
39 *
40 * @deprecated (<i>for internal use</i>) This class sohuld have visibility reduced to the package level
41 */
42 @Deprecated
43 @SuppressWarnings("javadoc")
44 public class Quantize {
45
46 @Deprecated
47 class DoubleArrayPointer {
48
49 private final double[] array;
50
51 private int startIndex;
52
53 @Deprecated
54 DoubleArrayPointer(double[] arrayIn) {
55 array = arrayIn;
56 }
57
58 @Deprecated
59 public DoubleArrayPointer copy(long l) {
60 DoubleArrayPointer result = new DoubleArrayPointer(array);
61 result.startIndex = (int) l;
62 return result;
63 }
64
65 @Deprecated
66 public double get(int ii) {
67 return array[ii + startIndex];
68 }
69 }
70
71 private static final double DEFAULT_QUANT_LEVEL = 4.;
72
73 private static final double MAX_INT_AS_DOUBLE = Integer.MAX_VALUE;
74
75 private static final int MINIMUM_PIXEL_WIDTH = 9;
76
77 /**
78 * number of reserved values, starting with
79 */
80 private static final long N_RESERVED_VALUES = 10;
81
82 private static final int N4 = 4;
83
84 private static final int N6 = 6;
85
86 private static final double NOISE_2_MULTIPLICATOR = 1.0483579;
87
88 private static final double NOISE_3_MULTIPLICATOR = 0.6052697;
89
90 private static final double NOISE_5_MULTIPLICATOR = 0.1772048;
91
92 private final QuantizeOption parameter;
93
94 /**
95 * maximum non-null value
96 */
97 private double maxValue;
98
99 /**
100 * minimum non-null value
101 */
102 private double minValue;
103
104 /**
105 * number of good, non-null pixels?
106 */
107 private long ngood;
108
109 /**
110 * returned 2nd order MAD of all non-null pixels
111 */
112 private double noise2;
113
114 /**
115 * returned 3rd order MAD of all non-null pixels
116 */
117 private double noise3;
118
119 /* returned 5th order MAD of all non-null pixels */
120 private double noise5;
121
122 private double xmaxval;
123
124 private double xminval;
125
126 private double xnoise2;
127
128 private double xnoise3;
129
130 private double xnoise5;
131
132 @Deprecated
133 public Quantize(QuantizeOption quantizeOption) {
134 parameter = quantizeOption;
135 }
136
137 /**
138 * Estimate the median and background noise in the input image using 2nd, 3rd and 5th order Median Absolute
139 * Differences. The noise in the background of the image is calculated using the MAD algorithms developed for
140 * deriving the signal to noise ratio in spectra (see issue #42 of the ST-ECF newsletter,
141 * http://www.stecf.org/documents/newsletter/) 3rd order: noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) -
142 * flux(i-2) - flux(i+2))) The returned estimates are the median of the values that are computed for each row of the
143 * image.
144 *
145 * @param arrayIn 2 dimensional tiledImageOperation of image pixels
146 * @param nx number of pixels in each row of the image
147 * @param ny number of rows in the image
148 * @param nullcheck check for null values, if true
149 * @param nullvalue value of null pixels, if nullcheck is true
150 */
151 private void calculateNoise(double[] arrayIn, int nx, int ny) {
152 DoubleArrayPointer array = new DoubleArrayPointer(arrayIn);
153 initializeNoise();
154 if (nx < MINIMUM_PIXEL_WIDTH) {
155 // treat entire tiledImageOperation as an image with a single row
156 nx = nx * ny;
157 ny = 1;
158 }
159 if (calculateNoiseShortRow(array, nx, ny)) {
160 return;
161 }
162 DoubleArrayPointer rowpix;
163 int nrows = 0, nrows2 = 0;
164 long ngoodpix = 0;
165 /* allocate arrays used to compute the median and noise estimates */
166 double[] differences2 = new double[nx];
167 double[] differences3 = new double[nx];
168 double[] differences5 = new double[nx];
169 double[] diffs2 = new double[ny];
170 double[] diffs3 = new double[ny];
171 double[] diffs5 = new double[ny];
172 /* loop over each row of the image */
173 for (int jj = 0; jj < ny; jj++) {
174 rowpix = array.copy(jj * nx); /* point to first pixel in the row */
175 int ii = 0;
176 ii = findNextValidPixelWithNullCheck(nx, rowpix, ii);
177 if (ii == nx) {
178 continue; /* hit end of row */
179 }
180 double v1 = getNextPixelAndCheckMinMax(rowpix, ii);
181 ngoodpix++;
182 ii = findNextValidPixelWithNullCheck(nx, rowpix, ++ii);
183 if (ii == nx) {
184 continue; /* hit end of row */
185 }
186 double v2 = getNextPixelAndCheckMinMax(rowpix, ii);
187 ngoodpix++;
188 ii = findNextValidPixelWithNullCheck(nx, rowpix, ++ii);
189 if (ii == nx) {
190 continue; /* hit end of row */
191 }
192 double v3 = getNextPixelAndCheckMinMax(rowpix, ii);
193 ngoodpix++;
194 ii = findNextValidPixelWithNullCheck(nx, rowpix, ++ii);
195 if (ii == nx) {
196 continue; /* hit end of row */
197 }
198 double v4 = getNextPixelAndCheckMinMax(rowpix, ii);
199 ngoodpix++;
200 ii = findNextValidPixelWithNullCheck(nx, rowpix, ++ii);
201 if (ii == nx) {
202 continue; /* hit end of row */
203 }
204 double v5 = getNextPixelAndCheckMinMax(rowpix, ii);
205 ngoodpix++;
206 ii = findNextValidPixelWithNullCheck(nx, rowpix, ++ii);
207 if (ii == nx) {
208 continue; /* hit end of row */
209 }
210 double v6 = getNextPixelAndCheckMinMax(rowpix, ii);
211 ngoodpix++;
212 ii = findNextValidPixelWithNullCheck(nx, rowpix, ++ii);
213 if (ii == nx) {
214 continue; /* hit end of row */
215 }
216 double v7 = getNextPixelAndCheckMinMax(rowpix, ii);
217 ngoodpix++;
218 ii = findNextValidPixelWithNullCheck(nx, rowpix, ++ii);
219 if (ii == nx) {
220 continue; /* hit end of row */
221 }
222 double v8 = getNextPixelAndCheckMinMax(rowpix, ii);
223 ngoodpix++;
224 // now populate the differences arrays for the remaining pixels in
225 // the row */
226 int nvals = 0;
227 int nvals2 = 0;
228 for (ii++; ii < nx; ii++) {
229 ii = findNextValidPixelWithNullCheck(nx, rowpix, ii);
230 if (ii == nx) {
231 continue; /* hit end of row */
232 }
233 double v9 = getNextPixelAndCheckMinMax(rowpix, ii);
234 /* construct tiledImageOperation of absolute differences */
235 if (!(v5 == v6 && v6 == v7)) {
236 differences2[nvals2] = Math.abs(v5 - v7);
237 nvals2++;
238 }
239 if (!(v3 == v4 && v4 == v5 && v5 == v6 && v6 == v7)) {
240 differences3[nvals] = Math.abs(2 * v5 - v3 - v7);
241 differences5[nvals] = Math.abs(N6 * v5 - N4 * v3 - N4 * v7 + v1 + v9);
242 nvals++;
243 } else {
244 /* ignore constant background regions */
245 ngoodpix++;
246 }
247 /* shift over 1 pixel */
248 v1 = v2;
249 v2 = v3;
250 v3 = v4;
251 v4 = v5;
252 v5 = v6;
253 v6 = v7;
254 v7 = v8;
255 v8 = v9;
256 } /* end of loop over pixels in the row */
257 // compute the median diffs Note that there are 8 more pixel values
258 // than there are diffs values.
259 ngoodpix += nvals;
260 if (nvals == 0) {
261 continue; /* cannot compute medians on this row */
262 }
263 if (nvals == 1) {
264 if (nvals2 == 1) {
265 diffs2[nrows2] = differences2[0];
266 nrows2++;
267 }
268 diffs3[nrows] = differences3[0];
269 diffs5[nrows] = differences5[0];
270 } else {
271 /* quick_select returns the median MUCH faster than using qsort */
272 if (nvals2 > 1) {
273 diffs2[nrows2] = quickSelect(differences2, nvals);
274 nrows2++;
275 }
276 diffs3[nrows] = quickSelect(differences3, nvals);
277 diffs5[nrows] = quickSelect(differences5, nvals);
278 }
279 nrows++;
280 } /* end of loop over rows */
281 computeMedianOfValuesEachRow(nrows, nrows2, diffs2, diffs3, diffs5);
282 setNoiseResult(ngoodpix);
283 }
284
285 private boolean calculateNoiseShortRow(DoubleArrayPointer array, int nx, int ny) {
286 /* rows must have at least 9 pixels */
287 if (nx < MINIMUM_PIXEL_WIDTH) {
288 int ngoodpix = 0;
289 for (int index = 0; index < nx; index++) {
290 if (isNull(array.get(index))) {
291 continue;
292 }
293 if (array.get(index) < xminval) {
294 xminval = array.get(index);
295 }
296 if (array.get(index) > xmaxval) {
297 xmaxval = array.get(index);
298 }
299 ngoodpix++;
300 }
301 setNoiseResult(ngoodpix);
302 return true;
303 }
304 return false;
305 }
306
307 @Deprecated
308 protected void computeMedianOfValuesEachRow(int nrows, int nrows2, double[] diffs2, double[] diffs3, double[] diffs5) {
309 // compute median of the values for each row.
310 if (nrows == 0) {
311 xnoise3 = 0;
312 xnoise5 = 0;
313 } else if (nrows == 1) {
314 xnoise3 = diffs3[0];
315 xnoise5 = diffs5[0];
316 } else {
317 Arrays.sort(diffs3, 0, nrows);
318 Arrays.sort(diffs5, 0, nrows);
319 xnoise3 = (diffs3[(nrows - 1) / 2] + diffs3[nrows / 2]) / 2.;
320 xnoise5 = (diffs5[(nrows - 1) / 2] + diffs5[nrows / 2]) / 2.;
321 }
322 if (nrows2 == 0) {
323 xnoise2 = 0;
324 } else if (nrows2 == 1) {
325 xnoise2 = diffs2[0];
326 } else {
327 Arrays.sort(diffs2, 0, nrows2);
328 xnoise2 = (diffs2[(nrows2 - 1) / 2] + diffs2[nrows2 / 2]) / 2.;
329 }
330 }
331
332 @Deprecated
333 protected int findNextValidPixelWithNullCheck(int nx, DoubleArrayPointer rowpix, int ii) {
334 return ii;
335 }
336
337 private double getNextPixelAndCheckMinMax(DoubleArrayPointer rowpix, int ii) {
338 double pixelValue = rowpix.get(ii); /* store the good pixel value */
339 if (pixelValue < xminval) {
340 xminval = pixelValue;
341 }
342 if (pixelValue > xmaxval) {
343 xmaxval = pixelValue;
344 }
345 return pixelValue;
346 }
347
348 @Deprecated
349 protected double getNoise2() {
350 return noise2;
351 }
352
353 @Deprecated
354 protected double getNoise3() {
355 return noise3;
356 }
357
358 @Deprecated
359 protected double getNoise5() {
360 return noise5;
361 }
362
363 private void initializeNoise() {
364 xnoise2 = 0;
365 xnoise3 = 0;
366 xnoise5 = 0;
367 xminval = Double.MAX_VALUE;
368 xmaxval = Double.MIN_VALUE;
369 }
370
371 @Deprecated
372 protected boolean isNull(double d) {
373 return false;
374 }
375
376 /**
377 * arguments: long row i: tile number = row number in the binary table double fdata[] i: tiledImageOperation of
378 * image pixels to be compressed long nxpix i: number of pixels in each row of fdata long nypix i: number of rows in
379 * fdata nullcheck i: check for nullvalues in fdata? double in_null_value i: value used to represent undefined
380 * pixels in fdata float qlevel i: quantization level int dither_method i; which dithering method to use int idata[]
381 * o: values of fdata after applying bzero and bscale double bscale o: scale factor double bzero o: zero offset int
382 * iminval o: minimum quantized value that is returned int imaxval o: maximum quantized value that is returned The
383 * function value will be one if the input fdata were copied to idata; in this case the parameters bscale and bzero
384 * can be used to convert back to nearly the original floating point values: fdata ~= idata * bscale + bzero. If the
385 * function value is zero, the data were not copied to idata.
386 * <p>
387 * In earlier implementations of the compression code, we only used the noise3 value as the most reliable estimate
388 * of the background noise in an image. If it is not possible to compute a noise3 value, then this serves as a red
389 * flag to indicate that quantizing the image could cause a loss of significant information in the image.
390 * </p>
391 * <p>
392 * At some later date, we decided to take the more conservative approach of using the minimum of all three of the
393 * noise values (while still requiring that noise3 has a defined value) as the best estimate of the noise. Note that
394 * if an image contains pure Gaussian distributed noise, then noise2, noise3, and noise5 will have exactly the same
395 * value (within statistical measurement errors).
396 * </p>
397 *
398 * @param fdata the data to quantinize
399 * @param nxpix the image width
400 * @param nypix the image hight
401 *
402 * @return true if the quantification was possible
403 */
404 @Deprecated
405 public boolean quantize(double[] fdata, int nxpix, int nypix) {
406 // MAD 2nd, 3rd, and 5th order noise values
407 double stdev;
408 double bScale; /* bscale, 1 in intdata = delta in fdata */
409
410 long nx = (long) nxpix * (long) nypix;
411 if (nx <= 1L) {
412 parameter.setBScale(1.);
413 parameter.setBZero(0.);
414 return false;
415 }
416 if (parameter.getQLevel() >= 0.) {
417 /* estimate background noise using MAD pixel differences */
418 calculateNoise(fdata, nxpix, nypix);
419 // special case of an image filled with Nulls
420 if (parameter.isCheckNull() && ngood == 0) {
421 /* set parameters to dummy values, which are not used */
422 minValue = 0.;
423 maxValue = 1.;
424 stdev = 1;
425 } else {
426 // use the minimum of noise2, noise3, and noise5 as the best
427 // noise value
428 stdev = noise3;
429 if (noise2 != 0. && noise2 < stdev) {
430 stdev = noise2;
431 }
432 if (noise5 != 0. && noise5 < stdev) {
433 stdev = noise5;
434 }
435 }
436 if (parameter.getQLevel() == 0.) {
437 bScale = stdev / DEFAULT_QUANT_LEVEL; /* default quantization */
438 } else {
439 bScale = stdev / parameter.getQLevel();
440 }
441 if (bScale == 0.) {
442 return false; /* don't quantize */
443 }
444 } else {
445 /* negative value represents the absolute quantization level */
446 bScale = -parameter.getQLevel();
447 /* only nned to calculate the min and max values */
448 calculateNoise(fdata, nxpix, nypix);
449 }
450 /* check that the range of quantized levels is not > range of int */
451 if ((maxValue - minValue) / bScale > 2. * MAX_INT_AS_DOUBLE - N_RESERVED_VALUES) {
452 return false; /* don't quantize */
453 }
454
455 parameter.setBScale(bScale);
456 parameter.setMinValue(minValue);
457 parameter.setMaxValue(maxValue);
458 parameter.setCheckNull(parameter.isCheckNull() && ngood != nx);
459 return true; /* yes, data have been quantized */
460 }
461
462 private double quickSelect(double[] arr, int n) {
463 int low, high;
464 int median;
465 int middle, ll, hh;
466
467 low = 0;
468 high = n - 1;
469 median = low + high >>> 1; // was (low + high) / 2;
470 for (;;) {
471 if (high <= low) {
472 return arr[median];
473 }
474
475 if (high == low + 1) { /* Two elements only */
476 if (arr[low] > arr[high]) {
477 swapElements(arr, low, high);
478 }
479 return arr[median];
480 }
481
482 /* Find median of low, middle and high items; swap into position low */
483 middle = low + high >>> 1; // was (low + high) / 2;
484 if (arr[middle] > arr[high]) {
485 swapElements(arr, middle, high);
486 }
487 if (arr[low] > arr[high]) {
488 swapElements(arr, low, high);
489 }
490 if (arr[middle] > arr[low]) {
491 swapElements(arr, middle, low);
492 }
493
494 /* Swap low item (now in position middle) into position (low+1) */
495 swapElements(arr, middle, low + 1);
496
497 /* Nibble from each end towards middle, swapping items when stuck */
498 ll = low + 1;
499 hh = high;
500 for (;;) {
501 do {
502 ll++;
503 } while (arr[low] > arr[ll]);
504 do {
505 hh--;
506 } while (arr[hh] > arr[low]);
507
508 if (hh < ll) {
509 break;
510 }
511
512 swapElements(arr, ll, hh);
513 }
514
515 /* Swap middle item (in position low) back into correct position */
516 swapElements(arr, low, hh);
517
518 /* Re-set active partition */
519 if (hh <= median) {
520 low = ll;
521 }
522 if (hh >= median) {
523 high = hh - 1;
524 }
525 }
526 }
527
528 private void setNoiseResult(long ngoodpix) {
529 minValue = xminval;
530 maxValue = xmaxval;
531 ngood = ngoodpix;
532 noise2 = NOISE_2_MULTIPLICATOR * xnoise2;
533 noise3 = NOISE_3_MULTIPLICATOR * xnoise3;
534 noise5 = NOISE_5_MULTIPLICATOR * xnoise5;
535 }
536
537 private void swapElements(double[] array, int one, int second) {
538 double value = array[one];
539 array[one] = array[second];
540 array[second] = value;
541 }
542
543 }