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
QuantizeLevels(uint8_t * const data,int width,int height,int num_levels,uint64_t * const sse)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