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