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