tor-browser

histogram_enc.c (44593B)

---

// Copyright 2012 Google Inc. All Rights Reserved.
//
// Use of this source code is governed by a BSD-style license
// that can be found in the COPYING file in the root of the source
// tree. An additional intellectual property rights grant can be found
// in the file PATENTS. All contributing project authors may
// be found in the AUTHORS file in the root of the source tree.
// -----------------------------------------------------------------------------
//
// Author: Jyrki Alakuijala (jyrki@google.com)
//
#ifdef HAVE_CONFIG_H
#include "src/webp/config.h"
#endif

#include <assert.h>
#include <stdlib.h>
#include <string.h>

#include "src/dsp/lossless.h"
#include "src/dsp/lossless_common.h"
#include "src/enc/backward_references_enc.h"
#include "src/enc/histogram_enc.h"
#include "src/enc/vp8i_enc.h"
#include "src/utils/utils.h"
#include "src/webp/encode.h"
#include "src/webp/format_constants.h"
#include "src/webp/types.h"

// Number of partitions for the three dominant (literal, red and blue) symbol
// costs.
#define NUM_PARTITIONS 4
// The size of the bin-hash corresponding to the three dominant costs.
#define BIN_SIZE (NUM_PARTITIONS * NUM_PARTITIONS * NUM_PARTITIONS)
// Maximum number of histograms allowed in greedy combining algorithm.
#define MAX_HISTO_GREEDY 100

// Enum to meaningfully access the elements of the Histogram arrays.
typedef enum {
 LITERAL = 0,
 RED,
 BLUE,
 ALPHA,
 DISTANCE
} HistogramIndex;

// Return the size of the histogram for a given cache_bits.
static int GetHistogramSize(int cache_bits) {
 const int literal_size = VP8LHistogramNumCodes(cache_bits);
 const size_t total_size = sizeof(VP8LHistogram) + sizeof(int) * literal_size;
 assert(total_size <= (size_t)0x7fffffff);
 return (int)total_size;
}

static void HistogramStatsClear(VP8LHistogram* const h) {
 int i;
 for (i = 0; i < 5; ++i) {
   h->trivial_symbol[i] = VP8L_NON_TRIVIAL_SYM;
   // By default, the histogram is assumed to be used.
   h->is_used[i] = 1;
 }
 h->bit_cost = 0;
 memset(h->costs, 0, sizeof(h->costs));
}

static void HistogramClear(VP8LHistogram* const h) {
 uint32_t* const literal = h->literal;
 const int cache_bits = h->palette_code_bits;
 const int histo_size = GetHistogramSize(cache_bits);
 memset(h, 0, histo_size);
 h->palette_code_bits = cache_bits;
 h->literal = literal;
 HistogramStatsClear(h);
}

// Swap two histogram pointers.
static void HistogramSwap(VP8LHistogram** const h1, VP8LHistogram** const h2) {
 VP8LHistogram* const tmp = *h1;
 *h1 = *h2;
 *h2 = tmp;
}

static void HistogramCopy(const VP8LHistogram* const src,
                         VP8LHistogram* const dst) {
 uint32_t* const dst_literal = dst->literal;
 const int dst_cache_bits = dst->palette_code_bits;
 const int literal_size = VP8LHistogramNumCodes(dst_cache_bits);
 const int histo_size = GetHistogramSize(dst_cache_bits);
 assert(src->palette_code_bits == dst_cache_bits);
 memcpy(dst, src, histo_size);
 dst->literal = dst_literal;
 memcpy(dst->literal, src->literal, literal_size * sizeof(*dst->literal));
}

void VP8LFreeHistogram(VP8LHistogram* const h) { WebPSafeFree(h); }

void VP8LFreeHistogramSet(VP8LHistogramSet* const histograms) {
 WebPSafeFree(histograms);
}

void VP8LHistogramCreate(VP8LHistogram* const h,
                        const VP8LBackwardRefs* const refs,
                        int palette_code_bits) {
 if (palette_code_bits >= 0) {
   h->palette_code_bits = palette_code_bits;
 }
 HistogramClear(h);
 VP8LHistogramStoreRefs(refs, /*distance_modifier=*/NULL,
                        /*distance_modifier_arg0=*/0, h);
}

void VP8LHistogramInit(VP8LHistogram* const h, int palette_code_bits,
                      int init_arrays) {
 h->palette_code_bits = palette_code_bits;
 if (init_arrays) {
   HistogramClear(h);
 } else {
   HistogramStatsClear(h);
 }
}

VP8LHistogram* VP8LAllocateHistogram(int cache_bits) {
 VP8LHistogram* histo = NULL;
 const int total_size = GetHistogramSize(cache_bits);
 uint8_t* const memory = (uint8_t*)WebPSafeMalloc(total_size, sizeof(*memory));
 if (memory == NULL) return NULL;
 histo = (VP8LHistogram*)memory;
 // 'literal' won't necessary be aligned.
 histo->literal = (uint32_t*)(memory + sizeof(VP8LHistogram));
 VP8LHistogramInit(histo, cache_bits, /*init_arrays=*/ 0);
 return histo;
}

// Resets the pointers of the histograms to point to the bit buffer in the set.
static void HistogramSetResetPointers(VP8LHistogramSet* const set,
                                     int cache_bits) {
 int i;
 const int histo_size = GetHistogramSize(cache_bits);
 uint8_t* memory = (uint8_t*) (set->histograms);
 memory += set->max_size * sizeof(*set->histograms);
 for (i = 0; i < set->max_size; ++i) {
   memory = (uint8_t*) WEBP_ALIGN(memory);
   set->histograms[i] = (VP8LHistogram*) memory;
   // 'literal' won't necessary be aligned.
   set->histograms[i]->literal = (uint32_t*)(memory + sizeof(VP8LHistogram));
   memory += histo_size;
 }
}

// Returns the total size of the VP8LHistogramSet.
static size_t HistogramSetTotalSize(int size, int cache_bits) {
 const int histo_size = GetHistogramSize(cache_bits);
 return (sizeof(VP8LHistogramSet) + size * (sizeof(VP8LHistogram*) +
         histo_size + WEBP_ALIGN_CST));
}

VP8LHistogramSet* VP8LAllocateHistogramSet(int size, int cache_bits) {
 int i;
 VP8LHistogramSet* set;
 const size_t total_size = HistogramSetTotalSize(size, cache_bits);
 uint8_t* memory = (uint8_t*)WebPSafeMalloc(total_size, sizeof(*memory));
 if (memory == NULL) return NULL;

 set = (VP8LHistogramSet*)memory;
 memory += sizeof(*set);
 set->histograms = (VP8LHistogram**)memory;
 set->max_size = size;
 set->size = size;
 HistogramSetResetPointers(set, cache_bits);
 for (i = 0; i < size; ++i) {
   VP8LHistogramInit(set->histograms[i], cache_bits, /*init_arrays=*/ 0);
 }
 return set;
}

void VP8LHistogramSetClear(VP8LHistogramSet* const set) {
 int i;
 const int cache_bits = set->histograms[0]->palette_code_bits;
 const int size = set->max_size;
 const size_t total_size = HistogramSetTotalSize(size, cache_bits);
 uint8_t* memory = (uint8_t*)set;

 memset(memory, 0, total_size);
 memory += sizeof(*set);
 set->histograms = (VP8LHistogram**)memory;
 set->max_size = size;
 set->size = size;
 HistogramSetResetPointers(set, cache_bits);
 for (i = 0; i < size; ++i) {
   set->histograms[i]->palette_code_bits = cache_bits;
 }
}

// Removes the histogram 'i' from 'set'.
static void HistogramSetRemoveHistogram(VP8LHistogramSet* const set, int i) {
 set->histograms[i] = set->histograms[set->size - 1];
 --set->size;
 assert(set->size > 0);
}

// -----------------------------------------------------------------------------

static void HistogramAddSinglePixOrCopy(
   VP8LHistogram* const histo, const PixOrCopy* const v,
   int (*const distance_modifier)(int, int), int distance_modifier_arg0) {
 if (PixOrCopyIsLiteral(v)) {
   ++histo->alpha[PixOrCopyLiteral(v, 3)];
   ++histo->red[PixOrCopyLiteral(v, 2)];
   ++histo->literal[PixOrCopyLiteral(v, 1)];
   ++histo->blue[PixOrCopyLiteral(v, 0)];
 } else if (PixOrCopyIsCacheIdx(v)) {
   const int literal_ix =
       NUM_LITERAL_CODES + NUM_LENGTH_CODES + PixOrCopyCacheIdx(v);
   assert(histo->palette_code_bits != 0);
   ++histo->literal[literal_ix];
 } else {
   int code, extra_bits;
   VP8LPrefixEncodeBits(PixOrCopyLength(v), &code, &extra_bits);
   ++histo->literal[NUM_LITERAL_CODES + code];
   if (distance_modifier == NULL) {
     VP8LPrefixEncodeBits(PixOrCopyDistance(v), &code, &extra_bits);
   } else {
     VP8LPrefixEncodeBits(
         distance_modifier(distance_modifier_arg0, PixOrCopyDistance(v)),
         &code, &extra_bits);
   }
   ++histo->distance[code];
 }
}

void VP8LHistogramStoreRefs(const VP8LBackwardRefs* const refs,
                           int (*const distance_modifier)(int, int),
                           int distance_modifier_arg0,
                           VP8LHistogram* const histo) {
 VP8LRefsCursor c = VP8LRefsCursorInit(refs);
 while (VP8LRefsCursorOk(&c)) {
   HistogramAddSinglePixOrCopy(histo, c.cur_pos, distance_modifier,
                               distance_modifier_arg0);
   VP8LRefsCursorNext(&c);
 }
}

// -----------------------------------------------------------------------------
// Entropy-related functions.

static WEBP_INLINE uint64_t BitsEntropyRefine(const VP8LBitEntropy* entropy) {
 uint64_t mix;
 if (entropy->nonzeros < 5) {
   if (entropy->nonzeros <= 1) {
     return 0;
   }
   // Two symbols, they will be 0 and 1 in a Huffman code.
   // Let's mix in a bit of entropy to favor good clustering when
   // distributions of these are combined.
   if (entropy->nonzeros == 2) {
     return DivRound(99 * ((uint64_t)entropy->sum << LOG_2_PRECISION_BITS) +
                         entropy->entropy,
                     100);
   }
   // No matter what the entropy says, we cannot be better than min_limit
   // with Huffman coding. I am mixing a bit of entropy into the
   // min_limit since it produces much better (~0.5 %) compression results
   // perhaps because of better entropy clustering.
   if (entropy->nonzeros == 3) {
     mix = 950;
   } else {
     mix = 700;  // nonzeros == 4.
   }
 } else {
   mix = 627;
 }

 {
   uint64_t min_limit = (uint64_t)(2 * entropy->sum - entropy->max_val)
                        << LOG_2_PRECISION_BITS;
   min_limit =
       DivRound(mix * min_limit + (1000 - mix) * entropy->entropy, 1000);
   return (entropy->entropy < min_limit) ? min_limit : entropy->entropy;
 }
}

uint64_t VP8LBitsEntropy(const uint32_t* const array, int n) {
 VP8LBitEntropy entropy;
 VP8LBitsEntropyUnrefined(array, n, &entropy);

 return BitsEntropyRefine(&entropy);
}

static uint64_t InitialHuffmanCost(void) {
 // Small bias because Huffman code length is typically not stored in
 // full length.
 static const uint64_t kHuffmanCodeOfHuffmanCodeSize = CODE_LENGTH_CODES * 3;
 // Subtract a bias of 9.1.
 return (kHuffmanCodeOfHuffmanCodeSize << LOG_2_PRECISION_BITS) -
        DivRound(91ll << LOG_2_PRECISION_BITS, 10);
}

// Finalize the Huffman cost based on streak numbers and length type (<3 or >=3)
static uint64_t FinalHuffmanCost(const VP8LStreaks* const stats) {
 // The constants in this function are empirical and got rounded from
 // their original values in 1/8 when switched to 1/1024.
 uint64_t retval = InitialHuffmanCost();
 // Second coefficient: Many zeros in the histogram are covered efficiently
 // by a run-length encode. Originally 2/8.
 uint32_t retval_extra = stats->counts[0] * 1600 + 240 * stats->streaks[0][1];
 // Second coefficient: Constant values are encoded less efficiently, but still
 // RLE'ed. Originally 6/8.
 retval_extra += stats->counts[1] * 2640 + 720 * stats->streaks[1][1];
 // 0s are usually encoded more efficiently than non-0s.
 // Originally 15/8.
 retval_extra += 1840 * stats->streaks[0][0];
 // Originally 26/8.
 retval_extra += 3360 * stats->streaks[1][0];
 return retval + ((uint64_t)retval_extra << (LOG_2_PRECISION_BITS - 10));
}

// Get the symbol entropy for the distribution 'population'.
// Set 'trivial_sym', if there's only one symbol present in the distribution.
static uint64_t PopulationCost(const uint32_t* const population, int length,
                              uint16_t* const trivial_sym,
                              uint8_t* const is_used) {
 VP8LBitEntropy bit_entropy;
 VP8LStreaks stats;
 VP8LGetEntropyUnrefined(population, length, &bit_entropy, &stats);
 if (trivial_sym != NULL) {
   *trivial_sym = (bit_entropy.nonzeros == 1) ? bit_entropy.nonzero_code
                                              : VP8L_NON_TRIVIAL_SYM;
 }
 if (is_used != NULL) {
   // The histogram is used if there is at least one non-zero streak.
   *is_used = (stats.streaks[1][0] != 0 || stats.streaks[1][1] != 0);
 }

 return BitsEntropyRefine(&bit_entropy) + FinalHuffmanCost(&stats);
}

static WEBP_INLINE void GetPopulationInfo(const VP8LHistogram* const histo,
                                         HistogramIndex index,
                                         const uint32_t** population,
                                         int* length) {
 switch (index) {
   case LITERAL:
     *population = histo->literal;
     *length = VP8LHistogramNumCodes(histo->palette_code_bits);
     break;
   case RED:
     *population = histo->red;
     *length = NUM_LITERAL_CODES;
     break;
   case BLUE:
     *population = histo->blue;
     *length = NUM_LITERAL_CODES;
     break;
   case ALPHA:
     *population = histo->alpha;
     *length = NUM_LITERAL_CODES;
     break;
   case DISTANCE:
     *population = histo->distance;
     *length = NUM_DISTANCE_CODES;
     break;
 }
}

// trivial_at_end is 1 if the two histograms only have one element that is
// non-zero: both the zero-th one, or both the last one.
// 'index' is the index of the symbol in the histogram (literal, red, blue,
// alpha, distance).
static WEBP_INLINE uint64_t GetCombinedEntropy(const VP8LHistogram* const h1,
                                              const VP8LHistogram* const h2,
                                              HistogramIndex index) {
 const uint32_t* X;
 const uint32_t* Y;
 int length;
 VP8LStreaks stats;
 VP8LBitEntropy bit_entropy;
 const int is_h1_used = h1->is_used[index];
 const int is_h2_used = h2->is_used[index];
 const int is_trivial = h1->trivial_symbol[index] != VP8L_NON_TRIVIAL_SYM &&
                        h1->trivial_symbol[index] == h2->trivial_symbol[index];

 if (is_trivial || !is_h1_used || !is_h2_used) {
   if (is_h1_used) return h1->costs[index];
   return h2->costs[index];
 }
 assert(is_h1_used && is_h2_used);

 GetPopulationInfo(h1, index, &X, &length);
 GetPopulationInfo(h2, index, &Y, &length);
 VP8LGetCombinedEntropyUnrefined(X, Y, length, &bit_entropy, &stats);
 return BitsEntropyRefine(&bit_entropy) + FinalHuffmanCost(&stats);
}

// Estimates the Entropy + Huffman + other block overhead size cost.
uint64_t VP8LHistogramEstimateBits(const VP8LHistogram* const h) {
 int i;
 uint64_t cost = 0;
 for (i = 0; i < 5; ++i) {
   int length;
   const uint32_t* population;
   GetPopulationInfo(h, (HistogramIndex)i, &population, &length);
   cost += PopulationCost(population, length, /*trivial_sym=*/NULL,
                          /*is_used=*/NULL);
 }
 cost += ((uint64_t)(VP8LExtraCost(h->literal + NUM_LITERAL_CODES,
                                   NUM_LENGTH_CODES) +
                     VP8LExtraCost(h->distance, NUM_DISTANCE_CODES))
          << LOG_2_PRECISION_BITS);
 return cost;
}

// -----------------------------------------------------------------------------
// Various histogram combine/cost-eval functions

// Set a + b in b, saturating at WEBP_INT64_MAX.
static WEBP_INLINE void SaturateAdd(uint64_t a, int64_t* b) {
 if (*b < 0 || (int64_t)a <= WEBP_INT64_MAX - *b) {
   *b += (int64_t)a;
 } else {
   *b = WEBP_INT64_MAX;
 }
}

// Returns 1 if the cost of the combined histogram is less than the threshold.
// Otherwise returns 0 and the cost is invalid due to early bail-out.
WEBP_NODISCARD static int GetCombinedHistogramEntropy(
   const VP8LHistogram* const a, const VP8LHistogram* const b,
   int64_t cost_threshold_in, uint64_t* cost, uint64_t costs[5]) {
 int i;
 const uint64_t cost_threshold = (uint64_t)cost_threshold_in;
 assert(a->palette_code_bits == b->palette_code_bits);
 if (cost_threshold_in <= 0) return 0;
 *cost = 0;

 // No need to add the extra cost for length and distance as it is a constant
 // that does not influence the histograms.
 for (i = 0; i < 5; ++i) {
   costs[i] = GetCombinedEntropy(a, b, (HistogramIndex)i);
   *cost += costs[i];
   if (*cost >= cost_threshold) return 0;
 }

 return 1;
}

static WEBP_INLINE void HistogramAdd(const VP8LHistogram* const h1,
                                    const VP8LHistogram* const h2,
                                    VP8LHistogram* const hout) {
 int i;
 assert(h1->palette_code_bits == h2->palette_code_bits);

 for (i = 0; i < 5; ++i) {
   int length;
   const uint32_t *p1, *p2, *pout_const;
   uint32_t* pout;
   GetPopulationInfo(h1, (HistogramIndex)i, &p1, &length);
   GetPopulationInfo(h2, (HistogramIndex)i, &p2, &length);
   GetPopulationInfo(hout, (HistogramIndex)i, &pout_const, &length);
   pout = (uint32_t*)pout_const;
   if (h2 == hout) {
     if (h1->is_used[i]) {
       if (hout->is_used[i]) {
         VP8LAddVectorEq(p1, pout, length);
       } else {
         memcpy(pout, p1, length * sizeof(pout[0]));
       }
     }
   } else {
     if (h1->is_used[i]) {
       if (h2->is_used[i]) {
         VP8LAddVector(p1, p2, pout, length);
       } else {
         memcpy(pout, p1, length * sizeof(pout[0]));
       }
     } else if (h2->is_used[i]) {
       memcpy(pout, p2, length * sizeof(pout[0]));
     } else {
       memset(pout, 0, length * sizeof(pout[0]));
     }
   }
 }

 for (i = 0; i < 5; ++i) {
   hout->trivial_symbol[i] = h1->trivial_symbol[i] == h2->trivial_symbol[i]
                                 ? h1->trivial_symbol[i]
                                 : VP8L_NON_TRIVIAL_SYM;
   hout->is_used[i] = h1->is_used[i] || h2->is_used[i];
 }
}

static void UpdateHistogramCost(uint64_t bit_cost, uint64_t costs[5],
                               VP8LHistogram* const h) {
 int i;
 h->bit_cost = bit_cost;
 for (i = 0; i < 5; ++i) {
   h->costs[i] = costs[i];
 }
}

// Performs out = a + b, computing the cost C(a+b) - C(a) - C(b) while comparing
// to the threshold value 'cost_threshold'. The score returned is
//  Score = C(a+b) - C(a) - C(b), where C(a) + C(b) is known and fixed.
// Since the previous score passed is 'cost_threshold', we only need to compare
// the partial cost against 'cost_threshold + C(a) + C(b)' to possibly bail-out
// early.
// Returns 1 if the cost is less than the threshold.
// Otherwise returns 0 and the cost is invalid due to early bail-out.
WEBP_NODISCARD static int HistogramAddEval(const VP8LHistogram* const a,
                                          const VP8LHistogram* const b,
                                          VP8LHistogram* const out,
                                          int64_t cost_threshold) {
 const uint64_t sum_cost = a->bit_cost + b->bit_cost;
 uint64_t bit_cost, costs[5];
 SaturateAdd(sum_cost, &cost_threshold);
 if (!GetCombinedHistogramEntropy(a, b, cost_threshold, &bit_cost, costs)) {
   return 0;
 }

 HistogramAdd(a, b, out);
 UpdateHistogramCost(bit_cost, costs, out);
 return 1;
}

// Same as HistogramAddEval(), except that the resulting histogram
// is not stored. Only the cost C(a+b) - C(a) is evaluated. We omit
// the term C(b) which is constant over all the evaluations.
// Returns 1 if the cost is less than the threshold.
// Otherwise returns 0 and the cost is invalid due to early bail-out.
WEBP_NODISCARD static int HistogramAddThresh(const VP8LHistogram* const a,
                                            const VP8LHistogram* const b,
                                            int64_t cost_threshold,
                                            int64_t* cost_out) {
 uint64_t cost, costs[5];
 assert(a != NULL && b != NULL);
 SaturateAdd(a->bit_cost, &cost_threshold);
 if (!GetCombinedHistogramEntropy(a, b, cost_threshold, &cost, costs)) {
   return 0;
 }

 *cost_out = (int64_t)cost - (int64_t)a->bit_cost;
 return 1;
}

// -----------------------------------------------------------------------------

// The structure to keep track of cost range for the three dominant entropy
// symbols.
typedef struct {
 uint64_t literal_max;
 uint64_t literal_min;
 uint64_t red_max;
 uint64_t red_min;
 uint64_t blue_max;
 uint64_t blue_min;
} DominantCostRange;

static void DominantCostRangeInit(DominantCostRange* const c) {
 c->literal_max = 0;
 c->literal_min = WEBP_UINT64_MAX;
 c->red_max = 0;
 c->red_min = WEBP_UINT64_MAX;
 c->blue_max = 0;
 c->blue_min = WEBP_UINT64_MAX;
}

static void UpdateDominantCostRange(
   const VP8LHistogram* const h, DominantCostRange* const c) {
 if (c->literal_max < h->costs[LITERAL]) c->literal_max = h->costs[LITERAL];
 if (c->literal_min > h->costs[LITERAL]) c->literal_min = h->costs[LITERAL];
 if (c->red_max < h->costs[RED]) c->red_max = h->costs[RED];
 if (c->red_min > h->costs[RED]) c->red_min = h->costs[RED];
 if (c->blue_max < h->costs[BLUE]) c->blue_max = h->costs[BLUE];
 if (c->blue_min > h->costs[BLUE]) c->blue_min = h->costs[BLUE];
}

static void ComputeHistogramCost(VP8LHistogram* const h) {
 int i;
 // No need to add the extra cost for length and distance as it is a constant
 // that does not influence the histograms.
 for (i = 0; i < 5; ++i) {
   const uint32_t* population;
   int length;
   GetPopulationInfo(h, i, &population, &length);
   h->costs[i] = PopulationCost(population, length, &h->trivial_symbol[i],
                                &h->is_used[i]);
 }
 h->bit_cost = h->costs[LITERAL] + h->costs[RED] + h->costs[BLUE] +
               h->costs[ALPHA] + h->costs[DISTANCE];
}

static int GetBinIdForEntropy(uint64_t min, uint64_t max, uint64_t val) {
 const uint64_t range = max - min;
 if (range > 0) {
   const uint64_t delta = val - min;
   return (int)((NUM_PARTITIONS - 1e-6) * delta / range);
 } else {
   return 0;
 }
}

static int GetHistoBinIndex(const VP8LHistogram* const h,
                           const DominantCostRange* const c, int low_effort) {
 int bin_id =
     GetBinIdForEntropy(c->literal_min, c->literal_max, h->costs[LITERAL]);
 assert(bin_id < NUM_PARTITIONS);
 if (!low_effort) {
   bin_id = bin_id * NUM_PARTITIONS +
            GetBinIdForEntropy(c->red_min, c->red_max, h->costs[RED]);
   bin_id = bin_id * NUM_PARTITIONS +
            GetBinIdForEntropy(c->blue_min, c->blue_max, h->costs[BLUE]);
   assert(bin_id < BIN_SIZE);
 }
 return bin_id;
}

// Construct the histograms from backward references.
static void HistogramBuild(
   int xsize, int histo_bits, const VP8LBackwardRefs* const backward_refs,
   VP8LHistogramSet* const image_histo) {
 int x = 0, y = 0;
 const int histo_xsize = VP8LSubSampleSize(xsize, histo_bits);
 VP8LHistogram** const histograms = image_histo->histograms;
 VP8LRefsCursor c = VP8LRefsCursorInit(backward_refs);
 assert(histo_bits > 0);
 VP8LHistogramSetClear(image_histo);
 while (VP8LRefsCursorOk(&c)) {
   const PixOrCopy* const v = c.cur_pos;
   const int ix = (y >> histo_bits) * histo_xsize + (x >> histo_bits);
   HistogramAddSinglePixOrCopy(histograms[ix], v, NULL, 0);
   x += PixOrCopyLength(v);
   while (x >= xsize) {
     x -= xsize;
     ++y;
   }
   VP8LRefsCursorNext(&c);
 }
}

// Copies the histograms and computes its bit_cost.
static void HistogramCopyAndAnalyze(VP8LHistogramSet* const orig_histo,
                                   VP8LHistogramSet* const image_histo) {
 int i;
 VP8LHistogram** const orig_histograms = orig_histo->histograms;
 VP8LHistogram** const histograms = image_histo->histograms;
 assert(image_histo->max_size == orig_histo->max_size);
 image_histo->size = 0;
 for (i = 0; i < orig_histo->max_size; ++i) {
   VP8LHistogram* const histo = orig_histograms[i];
   ComputeHistogramCost(histo);

   // Skip the histogram if it is completely empty, which can happen for tiles
   // with no information (when they are skipped because of LZ77).
   if (!histo->is_used[LITERAL] && !histo->is_used[RED] &&
       !histo->is_used[BLUE] && !histo->is_used[ALPHA] &&
       !histo->is_used[DISTANCE]) {
     // The first histogram is always used.
     assert(i > 0);
     orig_histograms[i] = NULL;
   } else {
     // Copy histograms from orig_histo[] to image_histo[].
     HistogramCopy(histo, histograms[image_histo->size]);
     ++image_histo->size;
   }
 }
}

// Partition histograms to different entropy bins for three dominant (literal,
// red and blue) symbol costs and compute the histogram aggregate bit_cost.
static void HistogramAnalyzeEntropyBin(VP8LHistogramSet* const image_histo,
                                      int low_effort) {
 int i;
 VP8LHistogram** const histograms = image_histo->histograms;
 const int histo_size = image_histo->size;
 DominantCostRange cost_range;
 DominantCostRangeInit(&cost_range);

 // Analyze the dominant (literal, red and blue) entropy costs.
 for (i = 0; i < histo_size; ++i) {
   UpdateDominantCostRange(histograms[i], &cost_range);
 }

 // bin-hash histograms on three of the dominant (literal, red and blue)
 // symbol costs and store the resulting bin_id for each histogram.
 for (i = 0; i < histo_size; ++i) {
   histograms[i]->bin_id =
       GetHistoBinIndex(histograms[i], &cost_range, low_effort);
 }
}

// Merges some histograms with same bin_id together if it's advantageous.
// Sets the remaining histograms to NULL.
// 'combine_cost_factor' has to be divided by 100.
static void HistogramCombineEntropyBin(VP8LHistogramSet* const image_histo,
                                      VP8LHistogram* cur_combo, int num_bins,
                                      int32_t combine_cost_factor,
                                      int low_effort) {
 VP8LHistogram** const histograms = image_histo->histograms;
 int idx;
 struct {
   int16_t first;    // position of the histogram that accumulates all
                     // histograms with the same bin_id
   uint16_t num_combine_failures;   // number of combine failures per bin_id
 } bin_info[BIN_SIZE];

 assert(num_bins <= BIN_SIZE);
 for (idx = 0; idx < num_bins; ++idx) {
   bin_info[idx].first = -1;
   bin_info[idx].num_combine_failures = 0;
 }

 for (idx = 0; idx < image_histo->size;) {
   const int bin_id = histograms[idx]->bin_id;
   const int first = bin_info[bin_id].first;
   if (first == -1) {
     bin_info[bin_id].first = idx;
     ++idx;
   } else if (low_effort) {
     HistogramAdd(histograms[idx], histograms[first], histograms[first]);
     HistogramSetRemoveHistogram(image_histo, idx);
   } else {
     // try to merge #idx into #first (both share the same bin_id)
     const uint64_t bit_cost = histograms[idx]->bit_cost;
     const int64_t bit_cost_thresh =
         -DivRound((int64_t)bit_cost * combine_cost_factor, 100);
     if (HistogramAddEval(histograms[first], histograms[idx], cur_combo,
                          bit_cost_thresh)) {
       const int max_combine_failures = 32;
       // Try to merge two histograms only if the combo is a trivial one or
       // the two candidate histograms are already non-trivial.
       // For some images, 'try_combine' turns out to be false for a lot of
       // histogram pairs. In that case, we fallback to combining
       // histograms as usual to avoid increasing the header size.
       int try_combine =
           cur_combo->trivial_symbol[RED] != VP8L_NON_TRIVIAL_SYM &&
           cur_combo->trivial_symbol[BLUE] != VP8L_NON_TRIVIAL_SYM &&
           cur_combo->trivial_symbol[ALPHA] != VP8L_NON_TRIVIAL_SYM;
       if (!try_combine) {
         try_combine =
             histograms[idx]->trivial_symbol[RED] == VP8L_NON_TRIVIAL_SYM ||
             histograms[idx]->trivial_symbol[BLUE] == VP8L_NON_TRIVIAL_SYM ||
             histograms[idx]->trivial_symbol[ALPHA] == VP8L_NON_TRIVIAL_SYM;
         try_combine &=
             histograms[first]->trivial_symbol[RED] == VP8L_NON_TRIVIAL_SYM ||
             histograms[first]->trivial_symbol[BLUE] == VP8L_NON_TRIVIAL_SYM ||
             histograms[first]->trivial_symbol[ALPHA] == VP8L_NON_TRIVIAL_SYM;
       }
       if (try_combine ||
           bin_info[bin_id].num_combine_failures >= max_combine_failures) {
         // move the (better) merged histogram to its final slot
         HistogramSwap(&cur_combo, &histograms[first]);
         HistogramSetRemoveHistogram(image_histo, idx);
       } else {
         ++bin_info[bin_id].num_combine_failures;
         ++idx;
       }
     } else {
       ++idx;
     }
   }
 }
 if (low_effort) {
   // for low_effort case, update the final cost when everything is merged
   for (idx = 0; idx < image_histo->size; ++idx) {
     ComputeHistogramCost(histograms[idx]);
   }
 }
}

// Implement a Lehmer random number generator with a multiplicative constant of
// 48271 and a modulo constant of 2^31 - 1.
static uint32_t MyRand(uint32_t* const seed) {
 *seed = (uint32_t)(((uint64_t)(*seed) * 48271u) % 2147483647u);
 assert(*seed > 0);
 return *seed;
}

// -----------------------------------------------------------------------------
// Histogram pairs priority queue

// Pair of histograms. Negative idx1 value means that pair is out-of-date.
typedef struct {
 int idx1;
 int idx2;
 int64_t cost_diff;
 uint64_t cost_combo;
 uint64_t costs[5];
} HistogramPair;

typedef struct {
 HistogramPair* queue;
 int size;
 int max_size;
} HistoQueue;

static int HistoQueueInit(HistoQueue* const histo_queue, const int max_size) {
 histo_queue->size = 0;
 histo_queue->max_size = max_size;
 // We allocate max_size + 1 because the last element at index "size" is
 // used as temporary data (and it could be up to max_size).
 histo_queue->queue = (HistogramPair*)WebPSafeMalloc(
     histo_queue->max_size + 1, sizeof(*histo_queue->queue));
 return histo_queue->queue != NULL;
}

static void HistoQueueClear(HistoQueue* const histo_queue) {
 assert(histo_queue != NULL);
 WebPSafeFree(histo_queue->queue);
 histo_queue->size = 0;
 histo_queue->max_size = 0;
}

// Pop a specific pair in the queue by replacing it with the last one
// and shrinking the queue.
static void HistoQueuePopPair(HistoQueue* const histo_queue,
                             HistogramPair* const pair) {
 assert(pair >= histo_queue->queue &&
        pair < (histo_queue->queue + histo_queue->size));
 assert(histo_queue->size > 0);
 *pair = histo_queue->queue[histo_queue->size - 1];
 --histo_queue->size;
}

// Check whether a pair in the queue should be updated as head or not.
static void HistoQueueUpdateHead(HistoQueue* const histo_queue,
                                HistogramPair* const pair) {
 assert(pair->cost_diff < 0);
 assert(pair >= histo_queue->queue &&
        pair < (histo_queue->queue + histo_queue->size));
 assert(histo_queue->size > 0);
 if (pair->cost_diff < histo_queue->queue[0].cost_diff) {
   // Replace the best pair.
   const HistogramPair tmp = histo_queue->queue[0];
   histo_queue->queue[0] = *pair;
   *pair = tmp;
 }
}

// Replaces the bad_id with good_id in the pair.
static void HistoQueueFixPair(int bad_id, int good_id,
                             HistogramPair* const pair) {
 if (pair->idx1 == bad_id) pair->idx1 = good_id;
 if (pair->idx2 == bad_id) pair->idx2 = good_id;
 if (pair->idx1 > pair->idx2) {
   const int tmp = pair->idx1;
   pair->idx1 = pair->idx2;
   pair->idx2 = tmp;
 }
}

// Update the cost diff and combo of a pair of histograms. This needs to be
// called when the histograms have been merged with a third one.
// Returns 1 if the cost diff is less than the threshold.
// Otherwise returns 0 and the cost is invalid due to early bail-out.
WEBP_NODISCARD static int HistoQueueUpdatePair(const VP8LHistogram* const h1,
                                              const VP8LHistogram* const h2,
                                              int64_t cost_threshold,
                                              HistogramPair* const pair) {
 const int64_t sum_cost = h1->bit_cost + h2->bit_cost;
 SaturateAdd(sum_cost, &cost_threshold);
 if (!GetCombinedHistogramEntropy(h1, h2, cost_threshold, &pair->cost_combo,
                                  pair->costs)) {
   return 0;
 }
 pair->cost_diff = (int64_t)pair->cost_combo - sum_cost;
 return 1;
}

// Create a pair from indices "idx1" and "idx2" provided its cost
// is inferior to "threshold", a negative entropy.
// It returns the cost of the pair, or 0 if it superior to threshold.
static int64_t HistoQueuePush(HistoQueue* const histo_queue,
                             VP8LHistogram** const histograms, int idx1,
                             int idx2, int64_t threshold) {
 const VP8LHistogram* h1;
 const VP8LHistogram* h2;
 HistogramPair pair;

 // Stop here if the queue is full.
 if (histo_queue->size == histo_queue->max_size) return 0;
 assert(threshold <= 0);
 if (idx1 > idx2) {
   const int tmp = idx2;
   idx2 = idx1;
   idx1 = tmp;
 }
 pair.idx1 = idx1;
 pair.idx2 = idx2;
 h1 = histograms[idx1];
 h2 = histograms[idx2];

 // Do not even consider the pair if it does not improve the entropy.
 if (!HistoQueueUpdatePair(h1, h2, threshold, &pair)) return 0;

 histo_queue->queue[histo_queue->size++] = pair;
 HistoQueueUpdateHead(histo_queue, &histo_queue->queue[histo_queue->size - 1]);

 return pair.cost_diff;
}

// -----------------------------------------------------------------------------

// Combines histograms by continuously choosing the one with the highest cost
// reduction.
static int HistogramCombineGreedy(VP8LHistogramSet* const image_histo) {
 int ok = 0;
 const int image_histo_size = image_histo->size;
 int i, j;
 VP8LHistogram** const histograms = image_histo->histograms;
 // Priority queue of histogram pairs.
 HistoQueue histo_queue;

 // image_histo_size^2 for the queue size is safe. If you look at
 // HistogramCombineGreedy, and imagine that UpdateQueueFront always pushes
 // data to the queue, you insert at most:
 // - image_histo_size*(image_histo_size-1)/2 (the first two for loops)
 // - image_histo_size - 1 in the last for loop at the first iteration of
 //   the while loop, image_histo_size - 2 at the second iteration ...
 //   therefore image_histo_size*(image_histo_size-1)/2 overall too
 if (!HistoQueueInit(&histo_queue, image_histo_size * image_histo_size)) {
   goto End;
 }

 // Initialize the queue.
 for (i = 0; i < image_histo_size; ++i) {
   for (j = i + 1; j < image_histo_size; ++j) {
     HistoQueuePush(&histo_queue, histograms, i, j, 0);
   }
 }

 while (histo_queue.size > 0) {
   const int idx1 = histo_queue.queue[0].idx1;
   const int idx2 = histo_queue.queue[0].idx2;
   HistogramAdd(histograms[idx2], histograms[idx1], histograms[idx1]);
   UpdateHistogramCost(histo_queue.queue[0].cost_combo,
                       histo_queue.queue[0].costs, histograms[idx1]);

   // Remove merged histogram.
   HistogramSetRemoveHistogram(image_histo, idx2);

   // Remove pairs intersecting the just combined best pair.
   for (i = 0; i < histo_queue.size;) {
     HistogramPair* const p = histo_queue.queue + i;
     if (p->idx1 == idx1 || p->idx2 == idx1 ||
         p->idx1 == idx2 || p->idx2 == idx2) {
       HistoQueuePopPair(&histo_queue, p);
     } else {
       HistoQueueFixPair(image_histo->size, idx2, p);
       HistoQueueUpdateHead(&histo_queue, p);
       ++i;
     }
   }

   // Push new pairs formed with combined histogram to the queue.
   for (i = 0; i < image_histo->size; ++i) {
     if (i == idx1) continue;
     HistoQueuePush(&histo_queue, image_histo->histograms, idx1, i, 0);
   }
 }

 ok = 1;

End:
 HistoQueueClear(&histo_queue);
 return ok;
}

// Perform histogram aggregation using a stochastic approach.
// 'do_greedy' is set to 1 if a greedy approach needs to be performed
// afterwards, 0 otherwise.
static int HistogramCombineStochastic(VP8LHistogramSet* const image_histo,
                                     int min_cluster_size,
                                     int* const do_greedy) {
 int j, iter;
 uint32_t seed = 1;
 int tries_with_no_success = 0;
 const int outer_iters = image_histo->size;
 const int num_tries_no_success = outer_iters / 2;
 VP8LHistogram** const histograms = image_histo->histograms;
 // Priority queue of histogram pairs. Its size of 'kHistoQueueSize'
 // impacts the quality of the compression and the speed: the smaller the
 // faster but the worse for the compression.
 HistoQueue histo_queue;
 const int kHistoQueueSize = 9;
 int ok = 0;

 if (image_histo->size < min_cluster_size) {
   *do_greedy = 1;
   return 1;
 }

 if (!HistoQueueInit(&histo_queue, kHistoQueueSize)) goto End;

 // Collapse similar histograms in 'image_histo'.
 for (iter = 0; iter < outer_iters && image_histo->size >= min_cluster_size &&
               ++tries_with_no_success < num_tries_no_success;
     ++iter) {
   int64_t best_cost =
       (histo_queue.size == 0) ? 0 : histo_queue.queue[0].cost_diff;
   int best_idx1 = -1, best_idx2 = 1;
   const uint32_t rand_range = (image_histo->size - 1) * (image_histo->size);
   // (image_histo->size) / 2 was chosen empirically. Less means faster but
   // worse compression.
   const int num_tries = (image_histo->size) / 2;

   // Pick random samples.
   for (j = 0; image_histo->size >= 2 && j < num_tries; ++j) {
     int64_t curr_cost;
     // Choose two different histograms at random and try to combine them.
     const uint32_t tmp = MyRand(&seed) % rand_range;
     uint32_t idx1 = tmp / (image_histo->size - 1);
     uint32_t idx2 = tmp % (image_histo->size - 1);
     if (idx2 >= idx1) ++idx2;

     // Calculate cost reduction on combination.
     curr_cost =
         HistoQueuePush(&histo_queue, histograms, idx1, idx2, best_cost);
     if (curr_cost < 0) {  // found a better pair?
       best_cost = curr_cost;
       // Empty the queue if we reached full capacity.
       if (histo_queue.size == histo_queue.max_size) break;
     }
   }
   if (histo_queue.size == 0) continue;

   // Get the best histograms.
   best_idx1 = histo_queue.queue[0].idx1;
   best_idx2 = histo_queue.queue[0].idx2;
   assert(best_idx1 < best_idx2);
   // Merge the histograms and remove best_idx2 from the queue.
   HistogramAdd(histograms[best_idx2], histograms[best_idx1],
                histograms[best_idx1]);
   UpdateHistogramCost(histo_queue.queue[0].cost_combo,
                       histo_queue.queue[0].costs, histograms[best_idx1]);
   HistogramSetRemoveHistogram(image_histo, best_idx2);
   // Parse the queue and update each pair that deals with best_idx1,
   // best_idx2 or image_histo_size.
   for (j = 0; j < histo_queue.size;) {
     HistogramPair* const p = histo_queue.queue + j;
     const int is_idx1_best = p->idx1 == best_idx1 || p->idx1 == best_idx2;
     const int is_idx2_best = p->idx2 == best_idx1 || p->idx2 == best_idx2;
     // The front pair could have been duplicated by a random pick so
     // check for it all the time nevertheless.
     if (is_idx1_best && is_idx2_best) {
       HistoQueuePopPair(&histo_queue, p);
       continue;
     }
     // Any pair containing one of the two best indices should only refer to
     // best_idx1. Its cost should also be updated.
     if (is_idx1_best || is_idx2_best) {
       HistoQueueFixPair(best_idx2, best_idx1, p);
       // Re-evaluate the cost of an updated pair.
       if (!HistoQueueUpdatePair(histograms[p->idx1], histograms[p->idx2], 0,
                                 p)) {
         HistoQueuePopPair(&histo_queue, p);
         continue;
       }
     }
     HistoQueueFixPair(image_histo->size, best_idx2, p);
     HistoQueueUpdateHead(&histo_queue, p);
     ++j;
   }
   tries_with_no_success = 0;
 }
 *do_greedy = (image_histo->size <= min_cluster_size);
 ok = 1;

End:
 HistoQueueClear(&histo_queue);
 return ok;
}

// -----------------------------------------------------------------------------
// Histogram refinement

// Find the best 'out' histogram for each of the 'in' histograms.
// At call-time, 'out' contains the histograms of the clusters.
// Note: we assume that out[]->bit_cost is already up-to-date.
static void HistogramRemap(const VP8LHistogramSet* const in,
                          VP8LHistogramSet* const out,
                          uint32_t* const symbols) {
 int i;
 VP8LHistogram** const in_histo = in->histograms;
 VP8LHistogram** const out_histo = out->histograms;
 const int in_size = out->max_size;
 const int out_size = out->size;
 if (out_size > 1) {
   for (i = 0; i < in_size; ++i) {
     int best_out = 0;
     int64_t best_bits = WEBP_INT64_MAX;
     int k;
     if (in_histo[i] == NULL) {
       // Arbitrarily set to the previous value if unused to help future LZ77.
       symbols[i] = symbols[i - 1];
       continue;
     }
     for (k = 0; k < out_size; ++k) {
       int64_t cur_bits;
       if (HistogramAddThresh(out_histo[k], in_histo[i], best_bits,
                              &cur_bits)) {
         best_bits = cur_bits;
         best_out = k;
       }
     }
     symbols[i] = best_out;
   }
 } else {
   assert(out_size == 1);
   for (i = 0; i < in_size; ++i) {
     symbols[i] = 0;
   }
 }

 // Recompute each out based on raw and symbols.
 VP8LHistogramSetClear(out);
 out->size = out_size;

 for (i = 0; i < in_size; ++i) {
   int idx;
   if (in_histo[i] == NULL) continue;
   idx = symbols[i];
   HistogramAdd(in_histo[i], out_histo[idx], out_histo[idx]);
 }
}

static int32_t GetCombineCostFactor(int histo_size, int quality) {
 int32_t combine_cost_factor = 16;
 if (quality < 90) {
   if (histo_size > 256) combine_cost_factor /= 2;
   if (histo_size > 512) combine_cost_factor /= 2;
   if (histo_size > 1024) combine_cost_factor /= 2;
   if (quality <= 50) combine_cost_factor /= 2;
 }
 return combine_cost_factor;
}

int VP8LGetHistoImageSymbols(int xsize, int ysize,
                            const VP8LBackwardRefs* const refs, int quality,
                            int low_effort, int histogram_bits, int cache_bits,
                            VP8LHistogramSet* const image_histo,
                            VP8LHistogram* const tmp_histo,
                            uint32_t* const histogram_symbols,
                            const WebPPicture* const pic, int percent_range,
                            int* const percent) {
 const int histo_xsize =
     histogram_bits ? VP8LSubSampleSize(xsize, histogram_bits) : 1;
 const int histo_ysize =
     histogram_bits ? VP8LSubSampleSize(ysize, histogram_bits) : 1;
 const int image_histo_raw_size = histo_xsize * histo_ysize;
 VP8LHistogramSet* const orig_histo =
     VP8LAllocateHistogramSet(image_histo_raw_size, cache_bits);
 // Don't attempt linear bin-partition heuristic for
 // histograms of small sizes (as bin_map will be very sparse) and
 // maximum quality q==100 (to preserve the compression gains at that level).
 const int entropy_combine_num_bins = low_effort ? NUM_PARTITIONS : BIN_SIZE;
 int entropy_combine;
 if (orig_histo == NULL) {
   WebPEncodingSetError(pic, VP8_ENC_ERROR_OUT_OF_MEMORY);
   goto Error;
 }

 // Construct the histograms from backward references.
 HistogramBuild(xsize, histogram_bits, refs, orig_histo);
 HistogramCopyAndAnalyze(orig_histo, image_histo);
 entropy_combine =
     (image_histo->size > entropy_combine_num_bins * 2) && (quality < 100);

 if (entropy_combine) {
   const int32_t combine_cost_factor =
       GetCombineCostFactor(image_histo_raw_size, quality);

   HistogramAnalyzeEntropyBin(image_histo, low_effort);
   // Collapse histograms with similar entropy.
   HistogramCombineEntropyBin(image_histo, tmp_histo, entropy_combine_num_bins,
                              combine_cost_factor, low_effort);
 }

 // Don't combine the histograms using stochastic and greedy heuristics for
 // low-effort compression mode.
 if (!low_effort || !entropy_combine) {
   // cubic ramp between 1 and MAX_HISTO_GREEDY:
   const int threshold_size =
       (int)(1 + DivRound(quality * quality * quality * (MAX_HISTO_GREEDY - 1),
                          100 * 100 * 100));
   int do_greedy;
   if (!HistogramCombineStochastic(image_histo, threshold_size, &do_greedy)) {
     WebPEncodingSetError(pic, VP8_ENC_ERROR_OUT_OF_MEMORY);
     goto Error;
   }
   if (do_greedy) {
     if (!HistogramCombineGreedy(image_histo)) {
       WebPEncodingSetError(pic, VP8_ENC_ERROR_OUT_OF_MEMORY);
       goto Error;
     }
   }
 }

 // Find the optimal map from original histograms to the final ones.
 HistogramRemap(orig_histo, image_histo, histogram_symbols);

 if (!WebPReportProgress(pic, *percent + percent_range, percent)) {
   goto Error;
 }

Error:
 VP8LFreeHistogramSet(orig_histo);
 return (pic->error_code == VP8_ENC_OK);
}