quant_levels_utils.c 4.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
  1. // Copyright 2011 Google Inc. All Rights Reserved.
  2. //
  3. // Use of this source code is governed by a BSD-style license
  4. // that can be found in the COPYING file in the root of the source
  5. // tree. An additional intellectual property rights grant can be found
  6. // in the file PATENTS. All contributing project authors may
  7. // be found in the AUTHORS file in the root of the source tree.
  8. // -----------------------------------------------------------------------------
  9. //
  10. // Quantize levels for specified number of quantization-levels ([2, 256]).
  11. // Min and max values are preserved (usual 0 and 255 for alpha plane).
  12. //
  13. // Author: Skal (pascal.massimino@gmail.com)
  14. #include <assert.h>
  15. #include "./quant_levels_utils.h"
  16. #define NUM_SYMBOLS 256
  17. #define MAX_ITER 6 // Maximum number of convergence steps.
  18. #define ERROR_THRESHOLD 1e-4 // MSE stopping criterion.
  19. // -----------------------------------------------------------------------------
  20. // Quantize levels.
  21. int QuantizeLevels(uint8_t* const data, int width, int height,
  22. int num_levels, uint64_t* const sse) {
  23. int freq[NUM_SYMBOLS] = { 0 };
  24. int q_level[NUM_SYMBOLS] = { 0 };
  25. double inv_q_level[NUM_SYMBOLS] = { 0 };
  26. int min_s = 255, max_s = 0;
  27. const size_t data_size = height * width;
  28. int i, num_levels_in, iter;
  29. double last_err = 1.e38, err = 0.;
  30. const double err_threshold = ERROR_THRESHOLD * data_size;
  31. if (data == NULL) {
  32. return 0;
  33. }
  34. if (width <= 0 || height <= 0) {
  35. return 0;
  36. }
  37. if (num_levels < 2 || num_levels > 256) {
  38. return 0;
  39. }
  40. {
  41. size_t n;
  42. num_levels_in = 0;
  43. for (n = 0; n < data_size; ++n) {
  44. num_levels_in += (freq[data[n]] == 0);
  45. if (min_s > data[n]) min_s = data[n];
  46. if (max_s < data[n]) max_s = data[n];
  47. ++freq[data[n]];
  48. }
  49. }
  50. if (num_levels_in <= num_levels) goto End; // nothing to do!
  51. // Start with uniformly spread centroids.
  52. for (i = 0; i < num_levels; ++i) {
  53. inv_q_level[i] = min_s + (double)(max_s - min_s) * i / (num_levels - 1);
  54. }
  55. // Fixed values. Won't be changed.
  56. q_level[min_s] = 0;
  57. q_level[max_s] = num_levels - 1;
  58. assert(inv_q_level[0] == min_s);
  59. assert(inv_q_level[num_levels - 1] == max_s);
  60. // k-Means iterations.
  61. for (iter = 0; iter < MAX_ITER; ++iter) {
  62. double q_sum[NUM_SYMBOLS] = { 0 };
  63. double q_count[NUM_SYMBOLS] = { 0 };
  64. int s, slot = 0;
  65. // Assign classes to representatives.
  66. for (s = min_s; s <= max_s; ++s) {
  67. // Keep track of the nearest neighbour 'slot'
  68. while (slot < num_levels - 1 &&
  69. 2 * s > inv_q_level[slot] + inv_q_level[slot + 1]) {
  70. ++slot;
  71. }
  72. if (freq[s] > 0) {
  73. q_sum[slot] += s * freq[s];
  74. q_count[slot] += freq[s];
  75. }
  76. q_level[s] = slot;
  77. }
  78. // Assign new representatives to classes.
  79. if (num_levels > 2) {
  80. for (slot = 1; slot < num_levels - 1; ++slot) {
  81. const double count = q_count[slot];
  82. if (count > 0.) {
  83. inv_q_level[slot] = q_sum[slot] / count;
  84. }
  85. }
  86. }
  87. // Compute convergence error.
  88. err = 0.;
  89. for (s = min_s; s <= max_s; ++s) {
  90. const double error = s - inv_q_level[q_level[s]];
  91. err += freq[s] * error * error;
  92. }
  93. // Check for convergence: we stop as soon as the error is no
  94. // longer improving.
  95. if (last_err - err < err_threshold) break;
  96. last_err = err;
  97. }
  98. // Remap the alpha plane to quantized values.
  99. {
  100. // double->int rounding operation can be costly, so we do it
  101. // once for all before remapping. We also perform the data[] -> slot
  102. // mapping, while at it (avoid one indirection in the final loop).
  103. uint8_t map[NUM_SYMBOLS];
  104. int s;
  105. size_t n;
  106. for (s = min_s; s <= max_s; ++s) {
  107. const int slot = q_level[s];
  108. map[s] = (uint8_t)(inv_q_level[slot] + .5);
  109. }
  110. // Final pass.
  111. for (n = 0; n < data_size; ++n) {
  112. data[n] = map[data[n]];
  113. }
  114. }
  115. End:
  116. // Store sum of squared error if needed.
  117. if (sse != NULL) *sse = (uint64_t)err;
  118. return 1;
  119. }