github.com/cellofellow/gopkg@v0.0.0-20140722061823-eec0544a62ad/image/webp/libwebp/src/utils/quant_levels.c (about)

     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