Home | History | Annotate | Download | only in enc
      1 // Copyright 2012 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 // Author: Jyrki Alakuijala (jyrki (at) google.com)
     11 //
     12 #ifdef HAVE_CONFIG_H
     13 #include "src/webp/config.h"
     14 #endif
     15 
     16 #include <math.h>
     17 
     18 #include "src/enc/backward_references_enc.h"
     19 #include "src/enc/histogram_enc.h"
     20 #include "src/dsp/lossless.h"
     21 #include "src/dsp/lossless_common.h"
     22 #include "src/utils/utils.h"
     23 
     24 #define MAX_COST 1.e38
     25 
     26 // Number of partitions for the three dominant (literal, red and blue) symbol
     27 // costs.
     28 #define NUM_PARTITIONS 4
     29 // The size of the bin-hash corresponding to the three dominant costs.
     30 #define BIN_SIZE (NUM_PARTITIONS * NUM_PARTITIONS * NUM_PARTITIONS)
     31 // Maximum number of histograms allowed in greedy combining algorithm.
     32 #define MAX_HISTO_GREEDY 100
     33 
     34 static void HistogramClear(VP8LHistogram* const p) {
     35   uint32_t* const literal = p->literal_;
     36   const int cache_bits = p->palette_code_bits_;
     37   const int histo_size = VP8LGetHistogramSize(cache_bits);
     38   memset(p, 0, histo_size);
     39   p->palette_code_bits_ = cache_bits;
     40   p->literal_ = literal;
     41 }
     42 
     43 // Swap two histogram pointers.
     44 static void HistogramSwap(VP8LHistogram** const A, VP8LHistogram** const B) {
     45   VP8LHistogram* const tmp = *A;
     46   *A = *B;
     47   *B = tmp;
     48 }
     49 
     50 static void HistogramCopy(const VP8LHistogram* const src,
     51                           VP8LHistogram* const dst) {
     52   uint32_t* const dst_literal = dst->literal_;
     53   const int dst_cache_bits = dst->palette_code_bits_;
     54   const int histo_size = VP8LGetHistogramSize(dst_cache_bits);
     55   assert(src->palette_code_bits_ == dst_cache_bits);
     56   memcpy(dst, src, histo_size);
     57   dst->literal_ = dst_literal;
     58 }
     59 
     60 int VP8LGetHistogramSize(int cache_bits) {
     61   const int literal_size = VP8LHistogramNumCodes(cache_bits);
     62   const size_t total_size = sizeof(VP8LHistogram) + sizeof(int) * literal_size;
     63   assert(total_size <= (size_t)0x7fffffff);
     64   return (int)total_size;
     65 }
     66 
     67 void VP8LFreeHistogram(VP8LHistogram* const histo) {
     68   WebPSafeFree(histo);
     69 }
     70 
     71 void VP8LFreeHistogramSet(VP8LHistogramSet* const histo) {
     72   WebPSafeFree(histo);
     73 }
     74 
     75 void VP8LHistogramStoreRefs(const VP8LBackwardRefs* const refs,
     76                             VP8LHistogram* const histo) {
     77   VP8LRefsCursor c = VP8LRefsCursorInit(refs);
     78   while (VP8LRefsCursorOk(&c)) {
     79     VP8LHistogramAddSinglePixOrCopy(histo, c.cur_pos, NULL, 0);
     80     VP8LRefsCursorNext(&c);
     81   }
     82 }
     83 
     84 void VP8LHistogramCreate(VP8LHistogram* const p,
     85                          const VP8LBackwardRefs* const refs,
     86                          int palette_code_bits) {
     87   if (palette_code_bits >= 0) {
     88     p->palette_code_bits_ = palette_code_bits;
     89   }
     90   HistogramClear(p);
     91   VP8LHistogramStoreRefs(refs, p);
     92 }
     93 
     94 void VP8LHistogramInit(VP8LHistogram* const p, int palette_code_bits) {
     95   p->palette_code_bits_ = palette_code_bits;
     96   HistogramClear(p);
     97 }
     98 
     99 VP8LHistogram* VP8LAllocateHistogram(int cache_bits) {
    100   VP8LHistogram* histo = NULL;
    101   const int total_size = VP8LGetHistogramSize(cache_bits);
    102   uint8_t* const memory = (uint8_t*)WebPSafeMalloc(total_size, sizeof(*memory));
    103   if (memory == NULL) return NULL;
    104   histo = (VP8LHistogram*)memory;
    105   // literal_ won't necessary be aligned.
    106   histo->literal_ = (uint32_t*)(memory + sizeof(VP8LHistogram));
    107   VP8LHistogramInit(histo, cache_bits);
    108   return histo;
    109 }
    110 
    111 VP8LHistogramSet* VP8LAllocateHistogramSet(int size, int cache_bits) {
    112   int i;
    113   VP8LHistogramSet* set;
    114   const int histo_size = VP8LGetHistogramSize(cache_bits);
    115   const size_t total_size =
    116       sizeof(*set) + size * (sizeof(*set->histograms) +
    117       histo_size + WEBP_ALIGN_CST);
    118   uint8_t* memory = (uint8_t*)WebPSafeMalloc(total_size, sizeof(*memory));
    119   if (memory == NULL) return NULL;
    120 
    121   set = (VP8LHistogramSet*)memory;
    122   memory += sizeof(*set);
    123   set->histograms = (VP8LHistogram**)memory;
    124   memory += size * sizeof(*set->histograms);
    125   set->max_size = size;
    126   set->size = size;
    127   for (i = 0; i < size; ++i) {
    128     memory = (uint8_t*)WEBP_ALIGN(memory);
    129     set->histograms[i] = (VP8LHistogram*)memory;
    130     // literal_ won't necessary be aligned.
    131     set->histograms[i]->literal_ = (uint32_t*)(memory + sizeof(VP8LHistogram));
    132     VP8LHistogramInit(set->histograms[i], cache_bits);
    133     memory += histo_size;
    134   }
    135   return set;
    136 }
    137 
    138 // -----------------------------------------------------------------------------
    139 
    140 void VP8LHistogramAddSinglePixOrCopy(VP8LHistogram* const histo,
    141                                      const PixOrCopy* const v,
    142                                      int (*const distance_modifier)(int, int),
    143                                      int distance_modifier_arg0) {
    144   if (PixOrCopyIsLiteral(v)) {
    145     ++histo->alpha_[PixOrCopyLiteral(v, 3)];
    146     ++histo->red_[PixOrCopyLiteral(v, 2)];
    147     ++histo->literal_[PixOrCopyLiteral(v, 1)];
    148     ++histo->blue_[PixOrCopyLiteral(v, 0)];
    149   } else if (PixOrCopyIsCacheIdx(v)) {
    150     const int literal_ix =
    151         NUM_LITERAL_CODES + NUM_LENGTH_CODES + PixOrCopyCacheIdx(v);
    152     ++histo->literal_[literal_ix];
    153   } else {
    154     int code, extra_bits;
    155     VP8LPrefixEncodeBits(PixOrCopyLength(v), &code, &extra_bits);
    156     ++histo->literal_[NUM_LITERAL_CODES + code];
    157     if (distance_modifier == NULL) {
    158       VP8LPrefixEncodeBits(PixOrCopyDistance(v), &code, &extra_bits);
    159     } else {
    160       VP8LPrefixEncodeBits(
    161           distance_modifier(distance_modifier_arg0, PixOrCopyDistance(v)),
    162           &code, &extra_bits);
    163     }
    164     ++histo->distance_[code];
    165   }
    166 }
    167 
    168 // -----------------------------------------------------------------------------
    169 // Entropy-related functions.
    170 
    171 static WEBP_INLINE double BitsEntropyRefine(const VP8LBitEntropy* entropy) {
    172   double mix;
    173   if (entropy->nonzeros < 5) {
    174     if (entropy->nonzeros <= 1) {
    175       return 0;
    176     }
    177     // Two symbols, they will be 0 and 1 in a Huffman code.
    178     // Let's mix in a bit of entropy to favor good clustering when
    179     // distributions of these are combined.
    180     if (entropy->nonzeros == 2) {
    181       return 0.99 * entropy->sum + 0.01 * entropy->entropy;
    182     }
    183     // No matter what the entropy says, we cannot be better than min_limit
    184     // with Huffman coding. I am mixing a bit of entropy into the
    185     // min_limit since it produces much better (~0.5 %) compression results
    186     // perhaps because of better entropy clustering.
    187     if (entropy->nonzeros == 3) {
    188       mix = 0.95;
    189     } else {
    190       mix = 0.7;  // nonzeros == 4.
    191     }
    192   } else {
    193     mix = 0.627;
    194   }
    195 
    196   {
    197     double min_limit = 2 * entropy->sum - entropy->max_val;
    198     min_limit = mix * min_limit + (1.0 - mix) * entropy->entropy;
    199     return (entropy->entropy < min_limit) ? min_limit : entropy->entropy;
    200   }
    201 }
    202 
    203 double VP8LBitsEntropy(const uint32_t* const array, int n,
    204                        uint32_t* const trivial_symbol) {
    205   VP8LBitEntropy entropy;
    206   VP8LBitsEntropyUnrefined(array, n, &entropy);
    207   if (trivial_symbol != NULL) {
    208     *trivial_symbol =
    209         (entropy.nonzeros == 1) ? entropy.nonzero_code : VP8L_NON_TRIVIAL_SYM;
    210   }
    211 
    212   return BitsEntropyRefine(&entropy);
    213 }
    214 
    215 static double InitialHuffmanCost(void) {
    216   // Small bias because Huffman code length is typically not stored in
    217   // full length.
    218   static const int kHuffmanCodeOfHuffmanCodeSize = CODE_LENGTH_CODES * 3;
    219   static const double kSmallBias = 9.1;
    220   return kHuffmanCodeOfHuffmanCodeSize - kSmallBias;
    221 }
    222 
    223 // Finalize the Huffman cost based on streak numbers and length type (<3 or >=3)
    224 static double FinalHuffmanCost(const VP8LStreaks* const stats) {
    225   // The constants in this function are experimental and got rounded from
    226   // their original values in 1/8 when switched to 1/1024.
    227   double retval = InitialHuffmanCost();
    228   // Second coefficient: Many zeros in the histogram are covered efficiently
    229   // by a run-length encode. Originally 2/8.
    230   retval += stats->counts[0] * 1.5625 + 0.234375 * stats->streaks[0][1];
    231   // Second coefficient: Constant values are encoded less efficiently, but still
    232   // RLE'ed. Originally 6/8.
    233   retval += stats->counts[1] * 2.578125 + 0.703125 * stats->streaks[1][1];
    234   // 0s are usually encoded more efficiently than non-0s.
    235   // Originally 15/8.
    236   retval += 1.796875 * stats->streaks[0][0];
    237   // Originally 26/8.
    238   retval += 3.28125 * stats->streaks[1][0];
    239   return retval;
    240 }
    241 
    242 // Get the symbol entropy for the distribution 'population'.
    243 // Set 'trivial_sym', if there's only one symbol present in the distribution.
    244 static double PopulationCost(const uint32_t* const population, int length,
    245                              uint32_t* const trivial_sym) {
    246   VP8LBitEntropy bit_entropy;
    247   VP8LStreaks stats;
    248   VP8LGetEntropyUnrefined(population, length, &bit_entropy, &stats);
    249   if (trivial_sym != NULL) {
    250     *trivial_sym = (bit_entropy.nonzeros == 1) ? bit_entropy.nonzero_code
    251                                                : VP8L_NON_TRIVIAL_SYM;
    252   }
    253 
    254   return BitsEntropyRefine(&bit_entropy) + FinalHuffmanCost(&stats);
    255 }
    256 
    257 // trivial_at_end is 1 if the two histograms only have one element that is
    258 // non-zero: both the zero-th one, or both the last one.
    259 static WEBP_INLINE double GetCombinedEntropy(const uint32_t* const X,
    260                                              const uint32_t* const Y,
    261                                              int length, int trivial_at_end) {
    262   VP8LStreaks stats;
    263   if (trivial_at_end) {
    264     // This configuration is due to palettization that transforms an indexed
    265     // pixel into 0xff000000 | (pixel << 8) in VP8LBundleColorMap.
    266     // BitsEntropyRefine is 0 for histograms with only one non-zero value.
    267     // Only FinalHuffmanCost needs to be evaluated.
    268     memset(&stats, 0, sizeof(stats));
    269     // Deal with the non-zero value at index 0 or length-1.
    270     stats.streaks[1][0] += 1;
    271     // Deal with the following/previous zero streak.
    272     stats.counts[0] += 1;
    273     stats.streaks[0][1] += length - 1;
    274     return FinalHuffmanCost(&stats);
    275   } else {
    276     VP8LBitEntropy bit_entropy;
    277     VP8LGetCombinedEntropyUnrefined(X, Y, length, &bit_entropy, &stats);
    278 
    279     return BitsEntropyRefine(&bit_entropy) + FinalHuffmanCost(&stats);
    280   }
    281 }
    282 
    283 // Estimates the Entropy + Huffman + other block overhead size cost.
    284 double VP8LHistogramEstimateBits(const VP8LHistogram* const p) {
    285   return
    286       PopulationCost(
    287           p->literal_, VP8LHistogramNumCodes(p->palette_code_bits_), NULL)
    288       + PopulationCost(p->red_, NUM_LITERAL_CODES, NULL)
    289       + PopulationCost(p->blue_, NUM_LITERAL_CODES, NULL)
    290       + PopulationCost(p->alpha_, NUM_LITERAL_CODES, NULL)
    291       + PopulationCost(p->distance_, NUM_DISTANCE_CODES, NULL)
    292       + VP8LExtraCost(p->literal_ + NUM_LITERAL_CODES, NUM_LENGTH_CODES)
    293       + VP8LExtraCost(p->distance_, NUM_DISTANCE_CODES);
    294 }
    295 
    296 // -----------------------------------------------------------------------------
    297 // Various histogram combine/cost-eval functions
    298 
    299 static int GetCombinedHistogramEntropy(const VP8LHistogram* const a,
    300                                        const VP8LHistogram* const b,
    301                                        double cost_threshold,
    302                                        double* cost) {
    303   const int palette_code_bits = a->palette_code_bits_;
    304   int trivial_at_end = 0;
    305   assert(a->palette_code_bits_ == b->palette_code_bits_);
    306   *cost += GetCombinedEntropy(a->literal_, b->literal_,
    307                               VP8LHistogramNumCodes(palette_code_bits), 0);
    308   *cost += VP8LExtraCostCombined(a->literal_ + NUM_LITERAL_CODES,
    309                                  b->literal_ + NUM_LITERAL_CODES,
    310                                  NUM_LENGTH_CODES);
    311   if (*cost > cost_threshold) return 0;
    312 
    313   if (a->trivial_symbol_ != VP8L_NON_TRIVIAL_SYM &&
    314       a->trivial_symbol_ == b->trivial_symbol_) {
    315     // A, R and B are all 0 or 0xff.
    316     const uint32_t color_a = (a->trivial_symbol_ >> 24) & 0xff;
    317     const uint32_t color_r = (a->trivial_symbol_ >> 16) & 0xff;
    318     const uint32_t color_b = (a->trivial_symbol_ >> 0) & 0xff;
    319     if ((color_a == 0 || color_a == 0xff) &&
    320         (color_r == 0 || color_r == 0xff) &&
    321         (color_b == 0 || color_b == 0xff)) {
    322       trivial_at_end = 1;
    323     }
    324   }
    325 
    326   *cost +=
    327       GetCombinedEntropy(a->red_, b->red_, NUM_LITERAL_CODES, trivial_at_end);
    328   if (*cost > cost_threshold) return 0;
    329 
    330   *cost +=
    331       GetCombinedEntropy(a->blue_, b->blue_, NUM_LITERAL_CODES, trivial_at_end);
    332   if (*cost > cost_threshold) return 0;
    333 
    334   *cost += GetCombinedEntropy(a->alpha_, b->alpha_, NUM_LITERAL_CODES,
    335                               trivial_at_end);
    336   if (*cost > cost_threshold) return 0;
    337 
    338   *cost +=
    339       GetCombinedEntropy(a->distance_, b->distance_, NUM_DISTANCE_CODES, 0);
    340   *cost +=
    341       VP8LExtraCostCombined(a->distance_, b->distance_, NUM_DISTANCE_CODES);
    342   if (*cost > cost_threshold) return 0;
    343 
    344   return 1;
    345 }
    346 
    347 static WEBP_INLINE void HistogramAdd(const VP8LHistogram* const a,
    348                                      const VP8LHistogram* const b,
    349                                      VP8LHistogram* const out) {
    350   VP8LHistogramAdd(a, b, out);
    351   out->trivial_symbol_ = (a->trivial_symbol_ == b->trivial_symbol_)
    352                        ? a->trivial_symbol_
    353                        : VP8L_NON_TRIVIAL_SYM;
    354 }
    355 
    356 // Performs out = a + b, computing the cost C(a+b) - C(a) - C(b) while comparing
    357 // to the threshold value 'cost_threshold'. The score returned is
    358 //  Score = C(a+b) - C(a) - C(b), where C(a) + C(b) is known and fixed.
    359 // Since the previous score passed is 'cost_threshold', we only need to compare
    360 // the partial cost against 'cost_threshold + C(a) + C(b)' to possibly bail-out
    361 // early.
    362 static double HistogramAddEval(const VP8LHistogram* const a,
    363                                const VP8LHistogram* const b,
    364                                VP8LHistogram* const out,
    365                                double cost_threshold) {
    366   double cost = 0;
    367   const double sum_cost = a->bit_cost_ + b->bit_cost_;
    368   cost_threshold += sum_cost;
    369 
    370   if (GetCombinedHistogramEntropy(a, b, cost_threshold, &cost)) {
    371     HistogramAdd(a, b, out);
    372     out->bit_cost_ = cost;
    373     out->palette_code_bits_ = a->palette_code_bits_;
    374   }
    375 
    376   return cost - sum_cost;
    377 }
    378 
    379 // Same as HistogramAddEval(), except that the resulting histogram
    380 // is not stored. Only the cost C(a+b) - C(a) is evaluated. We omit
    381 // the term C(b) which is constant over all the evaluations.
    382 static double HistogramAddThresh(const VP8LHistogram* const a,
    383                                  const VP8LHistogram* const b,
    384                                  double cost_threshold) {
    385   double cost = -a->bit_cost_;
    386   GetCombinedHistogramEntropy(a, b, cost_threshold, &cost);
    387   return cost;
    388 }
    389 
    390 // -----------------------------------------------------------------------------
    391 
    392 // The structure to keep track of cost range for the three dominant entropy
    393 // symbols.
    394 // TODO(skal): Evaluate if float can be used here instead of double for
    395 // representing the entropy costs.
    396 typedef struct {
    397   double literal_max_;
    398   double literal_min_;
    399   double red_max_;
    400   double red_min_;
    401   double blue_max_;
    402   double blue_min_;
    403 } DominantCostRange;
    404 
    405 static void DominantCostRangeInit(DominantCostRange* const c) {
    406   c->literal_max_ = 0.;
    407   c->literal_min_ = MAX_COST;
    408   c->red_max_ = 0.;
    409   c->red_min_ = MAX_COST;
    410   c->blue_max_ = 0.;
    411   c->blue_min_ = MAX_COST;
    412 }
    413 
    414 static void UpdateDominantCostRange(
    415     const VP8LHistogram* const h, DominantCostRange* const c) {
    416   if (c->literal_max_ < h->literal_cost_) c->literal_max_ = h->literal_cost_;
    417   if (c->literal_min_ > h->literal_cost_) c->literal_min_ = h->literal_cost_;
    418   if (c->red_max_ < h->red_cost_) c->red_max_ = h->red_cost_;
    419   if (c->red_min_ > h->red_cost_) c->red_min_ = h->red_cost_;
    420   if (c->blue_max_ < h->blue_cost_) c->blue_max_ = h->blue_cost_;
    421   if (c->blue_min_ > h->blue_cost_) c->blue_min_ = h->blue_cost_;
    422 }
    423 
    424 static void UpdateHistogramCost(VP8LHistogram* const h) {
    425   uint32_t alpha_sym, red_sym, blue_sym;
    426   const double alpha_cost =
    427       PopulationCost(h->alpha_, NUM_LITERAL_CODES, &alpha_sym);
    428   const double distance_cost =
    429       PopulationCost(h->distance_, NUM_DISTANCE_CODES, NULL) +
    430       VP8LExtraCost(h->distance_, NUM_DISTANCE_CODES);
    431   const int num_codes = VP8LHistogramNumCodes(h->palette_code_bits_);
    432   h->literal_cost_ = PopulationCost(h->literal_, num_codes, NULL) +
    433                      VP8LExtraCost(h->literal_ + NUM_LITERAL_CODES,
    434                                    NUM_LENGTH_CODES);
    435   h->red_cost_ = PopulationCost(h->red_, NUM_LITERAL_CODES, &red_sym);
    436   h->blue_cost_ = PopulationCost(h->blue_, NUM_LITERAL_CODES, &blue_sym);
    437   h->bit_cost_ = h->literal_cost_ + h->red_cost_ + h->blue_cost_ +
    438                  alpha_cost + distance_cost;
    439   if ((alpha_sym | red_sym | blue_sym) == VP8L_NON_TRIVIAL_SYM) {
    440     h->trivial_symbol_ = VP8L_NON_TRIVIAL_SYM;
    441   } else {
    442     h->trivial_symbol_ =
    443         ((uint32_t)alpha_sym << 24) | (red_sym << 16) | (blue_sym << 0);
    444   }
    445 }
    446 
    447 static int GetBinIdForEntropy(double min, double max, double val) {
    448   const double range = max - min;
    449   if (range > 0.) {
    450     const double delta = val - min;
    451     return (int)((NUM_PARTITIONS - 1e-6) * delta / range);
    452   } else {
    453     return 0;
    454   }
    455 }
    456 
    457 static int GetHistoBinIndex(const VP8LHistogram* const h,
    458                             const DominantCostRange* const c, int low_effort) {
    459   int bin_id = GetBinIdForEntropy(c->literal_min_, c->literal_max_,
    460                                   h->literal_cost_);
    461   assert(bin_id < NUM_PARTITIONS);
    462   if (!low_effort) {
    463     bin_id = bin_id * NUM_PARTITIONS
    464            + GetBinIdForEntropy(c->red_min_, c->red_max_, h->red_cost_);
    465     bin_id = bin_id * NUM_PARTITIONS
    466            + GetBinIdForEntropy(c->blue_min_, c->blue_max_, h->blue_cost_);
    467     assert(bin_id < BIN_SIZE);
    468   }
    469   return bin_id;
    470 }
    471 
    472 // Construct the histograms from backward references.
    473 static void HistogramBuild(
    474     int xsize, int histo_bits, const VP8LBackwardRefs* const backward_refs,
    475     VP8LHistogramSet* const image_histo) {
    476   int x = 0, y = 0;
    477   const int histo_xsize = VP8LSubSampleSize(xsize, histo_bits);
    478   VP8LHistogram** const histograms = image_histo->histograms;
    479   VP8LRefsCursor c = VP8LRefsCursorInit(backward_refs);
    480   assert(histo_bits > 0);
    481   while (VP8LRefsCursorOk(&c)) {
    482     const PixOrCopy* const v = c.cur_pos;
    483     const int ix = (y >> histo_bits) * histo_xsize + (x >> histo_bits);
    484     VP8LHistogramAddSinglePixOrCopy(histograms[ix], v, NULL, 0);
    485     x += PixOrCopyLength(v);
    486     while (x >= xsize) {
    487       x -= xsize;
    488       ++y;
    489     }
    490     VP8LRefsCursorNext(&c);
    491   }
    492 }
    493 
    494 // Copies the histograms and computes its bit_cost.
    495 static void HistogramCopyAndAnalyze(
    496     VP8LHistogramSet* const orig_histo, VP8LHistogramSet* const image_histo) {
    497   int i;
    498   const int histo_size = orig_histo->size;
    499   VP8LHistogram** const orig_histograms = orig_histo->histograms;
    500   VP8LHistogram** const histograms = image_histo->histograms;
    501   for (i = 0; i < histo_size; ++i) {
    502     VP8LHistogram* const histo = orig_histograms[i];
    503     UpdateHistogramCost(histo);
    504     // Copy histograms from orig_histo[] to image_histo[].
    505     HistogramCopy(histo, histograms[i]);
    506   }
    507 }
    508 
    509 // Partition histograms to different entropy bins for three dominant (literal,
    510 // red and blue) symbol costs and compute the histogram aggregate bit_cost.
    511 static void HistogramAnalyzeEntropyBin(VP8LHistogramSet* const image_histo,
    512                                        uint16_t* const bin_map,
    513                                        int low_effort) {
    514   int i;
    515   VP8LHistogram** const histograms = image_histo->histograms;
    516   const int histo_size = image_histo->size;
    517   DominantCostRange cost_range;
    518   DominantCostRangeInit(&cost_range);
    519 
    520   // Analyze the dominant (literal, red and blue) entropy costs.
    521   for (i = 0; i < histo_size; ++i) {
    522     UpdateDominantCostRange(histograms[i], &cost_range);
    523   }
    524 
    525   // bin-hash histograms on three of the dominant (literal, red and blue)
    526   // symbol costs and store the resulting bin_id for each histogram.
    527   for (i = 0; i < histo_size; ++i) {
    528     bin_map[i] = GetHistoBinIndex(histograms[i], &cost_range, low_effort);
    529   }
    530 }
    531 
    532 // Compact image_histo[] by merging some histograms with same bin_id together if
    533 // it's advantageous.
    534 static void HistogramCombineEntropyBin(VP8LHistogramSet* const image_histo,
    535                                        VP8LHistogram* cur_combo,
    536                                        const uint16_t* const bin_map,
    537                                        int bin_map_size, int num_bins,
    538                                        double combine_cost_factor,
    539                                        int low_effort) {
    540   VP8LHistogram** const histograms = image_histo->histograms;
    541   int idx;
    542   // Work in-place: processed histograms are put at the beginning of
    543   // image_histo[]. At the end, we just have to truncate the array.
    544   int size = 0;
    545   struct {
    546     int16_t first;    // position of the histogram that accumulates all
    547                       // histograms with the same bin_id
    548     uint16_t num_combine_failures;   // number of combine failures per bin_id
    549   } bin_info[BIN_SIZE];
    550 
    551   assert(num_bins <= BIN_SIZE);
    552   for (idx = 0; idx < num_bins; ++idx) {
    553     bin_info[idx].first = -1;
    554     bin_info[idx].num_combine_failures = 0;
    555   }
    556 
    557   for (idx = 0; idx < bin_map_size; ++idx) {
    558     const int bin_id = bin_map[idx];
    559     const int first = bin_info[bin_id].first;
    560     assert(size <= idx);
    561     if (first == -1) {
    562       // just move histogram #idx to its final position
    563       histograms[size] = histograms[idx];
    564       bin_info[bin_id].first = size++;
    565     } else if (low_effort) {
    566       HistogramAdd(histograms[idx], histograms[first], histograms[first]);
    567     } else {
    568       // try to merge #idx into #first (both share the same bin_id)
    569       const double bit_cost = histograms[idx]->bit_cost_;
    570       const double bit_cost_thresh = -bit_cost * combine_cost_factor;
    571       const double curr_cost_diff =
    572           HistogramAddEval(histograms[first], histograms[idx],
    573                            cur_combo, bit_cost_thresh);
    574       if (curr_cost_diff < bit_cost_thresh) {
    575         // Try to merge two histograms only if the combo is a trivial one or
    576         // the two candidate histograms are already non-trivial.
    577         // For some images, 'try_combine' turns out to be false for a lot of
    578         // histogram pairs. In that case, we fallback to combining
    579         // histograms as usual to avoid increasing the header size.
    580         const int try_combine =
    581             (cur_combo->trivial_symbol_ != VP8L_NON_TRIVIAL_SYM) ||
    582             ((histograms[idx]->trivial_symbol_ == VP8L_NON_TRIVIAL_SYM) &&
    583              (histograms[first]->trivial_symbol_ == VP8L_NON_TRIVIAL_SYM));
    584         const int max_combine_failures = 32;
    585         if (try_combine ||
    586             bin_info[bin_id].num_combine_failures >= max_combine_failures) {
    587           // move the (better) merged histogram to its final slot
    588           HistogramSwap(&cur_combo, &histograms[first]);
    589         } else {
    590           histograms[size++] = histograms[idx];
    591           ++bin_info[bin_id].num_combine_failures;
    592         }
    593       } else {
    594         histograms[size++] = histograms[idx];
    595       }
    596     }
    597   }
    598   image_histo->size = size;
    599   if (low_effort) {
    600     // for low_effort case, update the final cost when everything is merged
    601     for (idx = 0; idx < size; ++idx) {
    602       UpdateHistogramCost(histograms[idx]);
    603     }
    604   }
    605 }
    606 
    607 // Implement a Lehmer random number generator with a multiplicative constant of
    608 // 48271 and a modulo constant of 2^31  1.
    609 static uint32_t MyRand(uint32_t* const seed) {
    610   *seed = (uint32_t)(((uint64_t)(*seed) * 48271u) % 2147483647u);
    611   assert(*seed > 0);
    612   return *seed;
    613 }
    614 
    615 // -----------------------------------------------------------------------------
    616 // Histogram pairs priority queue
    617 
    618 // Pair of histograms. Negative idx1 value means that pair is out-of-date.
    619 typedef struct {
    620   int idx1;
    621   int idx2;
    622   double cost_diff;
    623   double cost_combo;
    624 } HistogramPair;
    625 
    626 typedef struct {
    627   HistogramPair* queue;
    628   int size;
    629   int max_size;
    630 } HistoQueue;
    631 
    632 static int HistoQueueInit(HistoQueue* const histo_queue, const int max_index) {
    633   histo_queue->size = 0;
    634   // max_index^2 for the queue size is safe. If you look at
    635   // HistogramCombineGreedy, and imagine that UpdateQueueFront always pushes
    636   // data to the queue, you insert at most:
    637   // - max_index*(max_index-1)/2 (the first two for loops)
    638   // - max_index - 1 in the last for loop at the first iteration of the while
    639   //   loop, max_index - 2 at the second iteration ... therefore
    640   //   max_index*(max_index-1)/2 overall too
    641   histo_queue->max_size = max_index * max_index;
    642   // We allocate max_size + 1 because the last element at index "size" is
    643   // used as temporary data (and it could be up to max_size).
    644   histo_queue->queue = (HistogramPair*)WebPSafeMalloc(
    645       histo_queue->max_size + 1, sizeof(*histo_queue->queue));
    646   return histo_queue->queue != NULL;
    647 }
    648 
    649 static void HistoQueueClear(HistoQueue* const histo_queue) {
    650   assert(histo_queue != NULL);
    651   WebPSafeFree(histo_queue->queue);
    652   histo_queue->size = 0;
    653   histo_queue->max_size = 0;
    654 }
    655 
    656 // Pop a specific pair in the queue by replacing it with the last one
    657 // and shrinking the queue.
    658 static void HistoQueuePopPair(HistoQueue* const histo_queue,
    659                               HistogramPair* const pair) {
    660   assert(pair >= histo_queue->queue &&
    661          pair < (histo_queue->queue + histo_queue->size));
    662   assert(histo_queue->size > 0);
    663   *pair = histo_queue->queue[histo_queue->size - 1];
    664   --histo_queue->size;
    665 }
    666 
    667 // Check whether a pair in the queue should be updated as head or not.
    668 static void HistoQueueUpdateHead(HistoQueue* const histo_queue,
    669                                  HistogramPair* const pair) {
    670   assert(pair->cost_diff < 0.);
    671   assert(pair >= histo_queue->queue &&
    672          pair < (histo_queue->queue + histo_queue->size));
    673   assert(histo_queue->size > 0);
    674   if (pair->cost_diff < histo_queue->queue[0].cost_diff) {
    675     // Replace the best pair.
    676     const HistogramPair tmp = histo_queue->queue[0];
    677     histo_queue->queue[0] = *pair;
    678     *pair = tmp;
    679   }
    680 }
    681 
    682 // Create a pair from indices "idx1" and "idx2" provided its cost
    683 // is inferior to "threshold", a negative entropy.
    684 // It returns the cost of the pair, or 0. if it superior to threshold.
    685 static double HistoQueuePush(HistoQueue* const histo_queue,
    686                              VP8LHistogram** const histograms, int idx1,
    687                              int idx2, double threshold) {
    688   const VP8LHistogram* h1;
    689   const VP8LHistogram* h2;
    690   HistogramPair pair;
    691   double sum_cost;
    692 
    693   assert(threshold <= 0.);
    694   if (idx1 > idx2) {
    695     const int tmp = idx2;
    696     idx2 = idx1;
    697     idx1 = tmp;
    698   }
    699   pair.idx1 = idx1;
    700   pair.idx2 = idx2;
    701   h1 = histograms[idx1];
    702   h2 = histograms[idx2];
    703   sum_cost = h1->bit_cost_ + h2->bit_cost_;
    704   pair.cost_combo = 0.;
    705   GetCombinedHistogramEntropy(h1, h2, sum_cost + threshold, &pair.cost_combo);
    706   pair.cost_diff = pair.cost_combo - sum_cost;
    707 
    708   // Do not even consider the pair if it does not improve the entropy.
    709   if (pair.cost_diff >= threshold) return 0.;
    710 
    711   // We cannot add more elements than the capacity.
    712   assert(histo_queue->size < histo_queue->max_size);
    713   histo_queue->queue[histo_queue->size++] = pair;
    714   HistoQueueUpdateHead(histo_queue, &histo_queue->queue[histo_queue->size - 1]);
    715 
    716   return pair.cost_diff;
    717 }
    718 
    719 // -----------------------------------------------------------------------------
    720 
    721 // Combines histograms by continuously choosing the one with the highest cost
    722 // reduction.
    723 static int HistogramCombineGreedy(VP8LHistogramSet* const image_histo) {
    724   int ok = 0;
    725   int image_histo_size = image_histo->size;
    726   int i, j;
    727   VP8LHistogram** const histograms = image_histo->histograms;
    728   // Indexes of remaining histograms.
    729   int* const clusters =
    730       (int*)WebPSafeMalloc(image_histo_size, sizeof(*clusters));
    731   // Priority queue of histogram pairs.
    732   HistoQueue histo_queue;
    733 
    734   if (!HistoQueueInit(&histo_queue, image_histo_size) || clusters == NULL) {
    735     goto End;
    736   }
    737 
    738   for (i = 0; i < image_histo_size; ++i) {
    739     // Initialize clusters indexes.
    740     clusters[i] = i;
    741     for (j = i + 1; j < image_histo_size; ++j) {
    742       // Initialize positions array.
    743       HistoQueuePush(&histo_queue, histograms, i, j, 0.);
    744     }
    745   }
    746 
    747   while (image_histo_size > 1 && histo_queue.size > 0) {
    748     const int idx1 = histo_queue.queue[0].idx1;
    749     const int idx2 = histo_queue.queue[0].idx2;
    750     HistogramAdd(histograms[idx2], histograms[idx1], histograms[idx1]);
    751     histograms[idx1]->bit_cost_ = histo_queue.queue[0].cost_combo;
    752     // Remove merged histogram.
    753     for (i = 0; i + 1 < image_histo_size; ++i) {
    754       if (clusters[i] >= idx2) {
    755         clusters[i] = clusters[i + 1];
    756       }
    757     }
    758     --image_histo_size;
    759 
    760     // Remove pairs intersecting the just combined best pair.
    761     for (i = 0; i < histo_queue.size;) {
    762       HistogramPair* const p = histo_queue.queue + i;
    763       if (p->idx1 == idx1 || p->idx2 == idx1 ||
    764           p->idx1 == idx2 || p->idx2 == idx2) {
    765         HistoQueuePopPair(&histo_queue, p);
    766       } else {
    767         HistoQueueUpdateHead(&histo_queue, p);
    768         ++i;
    769       }
    770     }
    771 
    772     // Push new pairs formed with combined histogram to the queue.
    773     for (i = 0; i < image_histo_size; ++i) {
    774       if (clusters[i] != idx1) {
    775         HistoQueuePush(&histo_queue, histograms, idx1, clusters[i], 0.);
    776       }
    777     }
    778   }
    779   // Move remaining histograms to the beginning of the array.
    780   for (i = 0; i < image_histo_size; ++i) {
    781     if (i != clusters[i]) {  // swap the two histograms
    782       HistogramSwap(&histograms[i], &histograms[clusters[i]]);
    783     }
    784   }
    785 
    786   image_histo->size = image_histo_size;
    787   ok = 1;
    788 
    789  End:
    790   WebPSafeFree(clusters);
    791   HistoQueueClear(&histo_queue);
    792   return ok;
    793 }
    794 
    795 // Perform histogram aggregation using a stochastic approach.
    796 // 'do_greedy' is set to 1 if a greedy approach needs to be performed
    797 // afterwards, 0 otherwise.
    798 static int HistogramCombineStochastic(VP8LHistogramSet* const image_histo,
    799                                       int min_cluster_size,
    800                                       int* const do_greedy) {
    801   int iter;
    802   uint32_t seed = 1;
    803   int tries_with_no_success = 0;
    804   int image_histo_size = image_histo->size;
    805   const int outer_iters = image_histo_size;
    806   const int num_tries_no_success = outer_iters / 2;
    807   VP8LHistogram** const histograms = image_histo->histograms;
    808   // Priority queue of histogram pairs. Its size of "kCostHeapSizeSqrt"^2
    809   // impacts the quality of the compression and the speed: the smaller the
    810   // faster but the worse for the compression.
    811   HistoQueue histo_queue;
    812   const int kHistoQueueSizeSqrt = 3;
    813   int ok = 0;
    814 
    815   if (!HistoQueueInit(&histo_queue, kHistoQueueSizeSqrt)) {
    816     goto End;
    817   }
    818   // Collapse similar histograms in 'image_histo'.
    819   ++min_cluster_size;
    820   for (iter = 0; iter < outer_iters && image_histo_size >= min_cluster_size &&
    821                  ++tries_with_no_success < num_tries_no_success;
    822        ++iter) {
    823     double best_cost =
    824         (histo_queue.size == 0) ? 0. : histo_queue.queue[0].cost_diff;
    825     int best_idx1 = -1, best_idx2 = 1;
    826     int j;
    827     const uint32_t rand_range = (image_histo_size - 1) * image_histo_size;
    828     // image_histo_size / 2 was chosen empirically. Less means faster but worse
    829     // compression.
    830     const int num_tries = image_histo_size / 2;
    831 
    832     for (j = 0; j < num_tries; ++j) {
    833       double curr_cost;
    834       // Choose two different histograms at random and try to combine them.
    835       const uint32_t tmp = MyRand(&seed) % rand_range;
    836       const uint32_t idx1 = tmp / (image_histo_size - 1);
    837       uint32_t idx2 = tmp % (image_histo_size - 1);
    838       if (idx2 >= idx1) ++idx2;
    839 
    840       // Calculate cost reduction on combination.
    841       curr_cost =
    842           HistoQueuePush(&histo_queue, histograms, idx1, idx2, best_cost);
    843       if (curr_cost < 0) {  // found a better pair?
    844         best_cost = curr_cost;
    845         // Empty the queue if we reached full capacity.
    846         if (histo_queue.size == histo_queue.max_size) break;
    847       }
    848     }
    849     if (histo_queue.size == 0) continue;
    850 
    851     // Merge the two best histograms.
    852     best_idx1 = histo_queue.queue[0].idx1;
    853     best_idx2 = histo_queue.queue[0].idx2;
    854     assert(best_idx1 < best_idx2);
    855     HistogramAddEval(histograms[best_idx1], histograms[best_idx2],
    856                      histograms[best_idx1], 0);
    857     // Swap the best_idx2 histogram with the last one (which is now unused).
    858     --image_histo_size;
    859     if (best_idx2 != image_histo_size) {
    860       HistogramSwap(&histograms[image_histo_size], &histograms[best_idx2]);
    861     }
    862     histograms[image_histo_size] = NULL;
    863     // Parse the queue and update each pair that deals with best_idx1,
    864     // best_idx2 or image_histo_size.
    865     for (j = 0; j < histo_queue.size;) {
    866       HistogramPair* const p = histo_queue.queue + j;
    867       const int is_idx1_best = p->idx1 == best_idx1 || p->idx1 == best_idx2;
    868       const int is_idx2_best = p->idx2 == best_idx1 || p->idx2 == best_idx2;
    869       int do_eval = 0;
    870       // The front pair could have been duplicated by a random pick so
    871       // check for it all the time nevertheless.
    872       if (is_idx1_best && is_idx2_best) {
    873         HistoQueuePopPair(&histo_queue, p);
    874         continue;
    875       }
    876       // Any pair containing one of the two best indices should only refer to
    877       // best_idx1. Its cost should also be updated.
    878       if (is_idx1_best) {
    879         p->idx1 = best_idx1;
    880         do_eval = 1;
    881       } else if (is_idx2_best) {
    882         p->idx2 = best_idx1;
    883         do_eval = 1;
    884       }
    885       if (p->idx2 == image_histo_size) {
    886         // No need to re-evaluate here as it does not involve a pair
    887         // containing best_idx1 or best_idx2.
    888         p->idx2 = best_idx2;
    889       }
    890       assert(p->idx2 < image_histo_size);
    891       // Make sure the index order is respected.
    892       if (p->idx1 > p->idx2) {
    893         const int tmp = p->idx2;
    894         p->idx2 = p->idx1;
    895         p->idx1 = tmp;
    896       }
    897       if (do_eval) {
    898         // Re-evaluate the cost of an updated pair.
    899         GetCombinedHistogramEntropy(histograms[p->idx1], histograms[p->idx2], 0,
    900                                     &p->cost_diff);
    901         if (p->cost_diff >= 0.) {
    902           HistoQueuePopPair(&histo_queue, p);
    903           continue;
    904         }
    905       }
    906       HistoQueueUpdateHead(&histo_queue, p);
    907       ++j;
    908     }
    909 
    910     tries_with_no_success = 0;
    911   }
    912   image_histo->size = image_histo_size;
    913   *do_greedy = (image_histo->size <= min_cluster_size);
    914   ok = 1;
    915 
    916 End:
    917   HistoQueueClear(&histo_queue);
    918   return ok;
    919 }
    920 
    921 // -----------------------------------------------------------------------------
    922 // Histogram refinement
    923 
    924 // Find the best 'out' histogram for each of the 'in' histograms.
    925 // Note: we assume that out[]->bit_cost_ is already up-to-date.
    926 static void HistogramRemap(const VP8LHistogramSet* const in,
    927                            const VP8LHistogramSet* const out,
    928                            uint16_t* const symbols) {
    929   int i;
    930   VP8LHistogram** const in_histo = in->histograms;
    931   VP8LHistogram** const out_histo = out->histograms;
    932   const int in_size = in->size;
    933   const int out_size = out->size;
    934   if (out_size > 1) {
    935     for (i = 0; i < in_size; ++i) {
    936       int best_out = 0;
    937       double best_bits = MAX_COST;
    938       int k;
    939       for (k = 0; k < out_size; ++k) {
    940         const double cur_bits =
    941             HistogramAddThresh(out_histo[k], in_histo[i], best_bits);
    942         if (k == 0 || cur_bits < best_bits) {
    943           best_bits = cur_bits;
    944           best_out = k;
    945         }
    946       }
    947       symbols[i] = best_out;
    948     }
    949   } else {
    950     assert(out_size == 1);
    951     for (i = 0; i < in_size; ++i) {
    952       symbols[i] = 0;
    953     }
    954   }
    955 
    956   // Recompute each out based on raw and symbols.
    957   for (i = 0; i < out_size; ++i) {
    958     HistogramClear(out_histo[i]);
    959   }
    960 
    961   for (i = 0; i < in_size; ++i) {
    962     const int idx = symbols[i];
    963     HistogramAdd(in_histo[i], out_histo[idx], out_histo[idx]);
    964   }
    965 }
    966 
    967 static double GetCombineCostFactor(int histo_size, int quality) {
    968   double combine_cost_factor = 0.16;
    969   if (quality < 90) {
    970     if (histo_size > 256) combine_cost_factor /= 2.;
    971     if (histo_size > 512) combine_cost_factor /= 2.;
    972     if (histo_size > 1024) combine_cost_factor /= 2.;
    973     if (quality <= 50) combine_cost_factor /= 2.;
    974   }
    975   return combine_cost_factor;
    976 }
    977 
    978 int VP8LGetHistoImageSymbols(int xsize, int ysize,
    979                              const VP8LBackwardRefs* const refs,
    980                              int quality, int low_effort,
    981                              int histo_bits, int cache_bits,
    982                              VP8LHistogramSet* const image_histo,
    983                              VP8LHistogram* const tmp_histo,
    984                              uint16_t* const histogram_symbols) {
    985   int ok = 0;
    986   const int histo_xsize = histo_bits ? VP8LSubSampleSize(xsize, histo_bits) : 1;
    987   const int histo_ysize = histo_bits ? VP8LSubSampleSize(ysize, histo_bits) : 1;
    988   const int image_histo_raw_size = histo_xsize * histo_ysize;
    989   VP8LHistogramSet* const orig_histo =
    990       VP8LAllocateHistogramSet(image_histo_raw_size, cache_bits);
    991   // Don't attempt linear bin-partition heuristic for
    992   // histograms of small sizes (as bin_map will be very sparse) and
    993   // maximum quality q==100 (to preserve the compression gains at that level).
    994   const int entropy_combine_num_bins = low_effort ? NUM_PARTITIONS : BIN_SIZE;
    995   const int entropy_combine =
    996       (orig_histo->size > entropy_combine_num_bins * 2) && (quality < 100);
    997 
    998   if (orig_histo == NULL) goto Error;
    999 
   1000   // Construct the histograms from backward references.
   1001   HistogramBuild(xsize, histo_bits, refs, orig_histo);
   1002   // Copies the histograms and computes its bit_cost.
   1003   HistogramCopyAndAnalyze(orig_histo, image_histo);
   1004 
   1005   if (entropy_combine) {
   1006     const int bin_map_size = orig_histo->size;
   1007     // Reuse histogram_symbols storage. By definition, it's guaranteed to be ok.
   1008     uint16_t* const bin_map = histogram_symbols;
   1009     const double combine_cost_factor =
   1010         GetCombineCostFactor(image_histo_raw_size, quality);
   1011 
   1012     HistogramAnalyzeEntropyBin(orig_histo, bin_map, low_effort);
   1013     // Collapse histograms with similar entropy.
   1014     HistogramCombineEntropyBin(image_histo, tmp_histo, bin_map, bin_map_size,
   1015                                entropy_combine_num_bins, combine_cost_factor,
   1016                                low_effort);
   1017   }
   1018 
   1019   // Don't combine the histograms using stochastic and greedy heuristics for
   1020   // low-effort compression mode.
   1021   if (!low_effort || !entropy_combine) {
   1022     const float x = quality / 100.f;
   1023     // cubic ramp between 1 and MAX_HISTO_GREEDY:
   1024     const int threshold_size = (int)(1 + (x * x * x) * (MAX_HISTO_GREEDY - 1));
   1025     int do_greedy;
   1026     if (!HistogramCombineStochastic(image_histo, threshold_size, &do_greedy)) {
   1027       goto Error;
   1028     }
   1029     if (do_greedy && !HistogramCombineGreedy(image_histo)) {
   1030       goto Error;
   1031     }
   1032   }
   1033 
   1034   // TODO(vikasa): Optimize HistogramRemap for low-effort compression mode also.
   1035   // Find the optimal map from original histograms to the final ones.
   1036   HistogramRemap(orig_histo, image_histo, histogram_symbols);
   1037 
   1038   ok = 1;
   1039 
   1040  Error:
   1041   VP8LFreeHistogramSet(orig_histo);
   1042   return ok;
   1043 }
   1044