Home | History | Annotate | Download | only in utils 
      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 (at) gmail.com)
     14 
     15 #include <assert.h>
     16 
     17 #include "./quant_levels.h"
     18 
     19 #if defined(__cplusplus) || defined(c_plusplus)
     20 extern "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 
     31 int 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