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