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 (at) gmail.com) 14 15 #include <assert.h> 16 17 #include "./quant_levels_utils.h" 18 19 #define NUM_SYMBOLS 256 20 21 #define MAX_ITER 6 // Maximum number of convergence steps. 22 #define ERROR_THRESHOLD 1e-4 // MSE stopping criterion. 23 24 // ----------------------------------------------------------------------------- 25 // Quantize levels. 26 27 int QuantizeLevels(uint8_t* const data, int width, int height, 28 int num_levels, uint64_t* const sse) { 29 int freq[NUM_SYMBOLS] = { 0 }; 30 int q_level[NUM_SYMBOLS] = { 0 }; 31 double inv_q_level[NUM_SYMBOLS] = { 0 }; 32 int min_s = 255, max_s = 0; 33 const size_t data_size = height * width; 34 int i, num_levels_in, iter; 35 double last_err = 1.e38, err = 0.; 36 const double err_threshold = ERROR_THRESHOLD * data_size; 37 38 if (data == NULL) { 39 return 0; 40 } 41 42 if (width <= 0 || height <= 0) { 43 return 0; 44 } 45 46 if (num_levels < 2 || num_levels > 256) { 47 return 0; 48 } 49 50 { 51 size_t n; 52 num_levels_in = 0; 53 for (n = 0; n < data_size; ++n) { 54 num_levels_in += (freq[data[n]] == 0); 55 if (min_s > data[n]) min_s = data[n]; 56 if (max_s < data[n]) max_s = data[n]; 57 ++freq[data[n]]; 58 } 59 } 60 61 if (num_levels_in <= num_levels) goto End; // nothing to do! 62 63 // Start with uniformly spread centroids. 64 for (i = 0; i < num_levels; ++i) { 65 inv_q_level[i] = min_s + (double)(max_s - min_s) * i / (num_levels - 1); 66 } 67 68 // Fixed values. Won't be changed. 69 q_level[min_s] = 0; 70 q_level[max_s] = num_levels - 1; 71 assert(inv_q_level[0] == min_s); 72 assert(inv_q_level[num_levels - 1] == max_s); 73 74 // k-Means iterations. 75 for (iter = 0; iter < MAX_ITER; ++iter) { 76 double q_sum[NUM_SYMBOLS] = { 0 }; 77 double q_count[NUM_SYMBOLS] = { 0 }; 78 int s, slot = 0; 79 80 // Assign classes to representatives. 81 for (s = min_s; s <= max_s; ++s) { 82 // Keep track of the nearest neighbour 'slot' 83 while (slot < num_levels - 1 && 84 2 * s > inv_q_level[slot] + inv_q_level[slot + 1]) { 85 ++slot; 86 } 87 if (freq[s] > 0) { 88 q_sum[slot] += s * freq[s]; 89 q_count[slot] += freq[s]; 90 } 91 q_level[s] = slot; 92 } 93 94 // Assign new representatives to classes. 95 if (num_levels > 2) { 96 for (slot = 1; slot < num_levels - 1; ++slot) { 97 const double count = q_count[slot]; 98 if (count > 0.) { 99 inv_q_level[slot] = q_sum[slot] / count; 100 } 101 } 102 } 103 104 // Compute convergence error. 105 err = 0.; 106 for (s = min_s; s <= max_s; ++s) { 107 const double error = s - inv_q_level[q_level[s]]; 108 err += freq[s] * error * error; 109 } 110 111 // Check for convergence: we stop as soon as the error is no 112 // longer improving. 113 if (last_err - err < err_threshold) break; 114 last_err = err; 115 } 116 117 // Remap the alpha plane to quantized values. 118 { 119 // double->int rounding operation can be costly, so we do it 120 // once for all before remapping. We also perform the data[] -> slot 121 // mapping, while at it (avoid one indirection in the final loop). 122 uint8_t map[NUM_SYMBOLS]; 123 int s; 124 size_t n; 125 for (s = min_s; s <= max_s; ++s) { 126 const int slot = q_level[s]; 127 map[s] = (uint8_t)(inv_q_level[slot] + .5); 128 } 129 // Final pass. 130 for (n = 0; n < data_size; ++n) { 131 data[n] = map[data[n]]; 132 } 133 } 134 End: 135 // Store sum of squared error if needed. 136 if (sse != NULL) *sse = (uint64_t)err; 137 138 return 1; 139 } 140 141
