third_party/libwebp/utils/quant_levels.c

   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 #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