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 | |
15 | #include <assert.h> |
16 | |
17 | #include "src/utils/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 | |