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