lossless: entropy clustering for high qualities.

Tested on 1000 pngs corpus with quality 90-100 it gives ~0.15% improvement
in compression density and ~7% speed up.

Change-Id: I460f56c96707edb3c1f0b51a024e5122e10458df

src/enc/histogram.c

@@ -27,6 +27,8 @@
 #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
 
 static void HistogramClear(VP8LHistogram* const p) {
   uint32_t* const literal = p->literal_;
@@ -560,8 +562,235 @@ static uint32_t MyRand(uint32_t *seed) {
   return *seed;
 }
 
-static void HistogramCombine(VP8LHistogramSet* const image_histo,
-                             VP8LHistogramSet* const histos, int quality) {
+// -----------------------------------------------------------------------------
+// Histogram pairs priority queue
+
+// Pair of histograms. Negative idx1 value means that pair is out-of-date.
+typedef struct {
+  int idx1;
+  int idx2;
+  double cost_diff;
+  double cost_combo;
+} HistogramPair;
+
+typedef struct {
+  HistogramPair* heap;
+  int* positions;
+  int size;
+  int max_index;
+} HistoHeap;
+
+static int HistoHeapInit(HistoHeap* const histo_heap, const int max_index) {
+  histo_heap->size = 0;
+  histo_heap->max_index = max_index;
+  histo_heap->heap = WebPSafeMalloc(max_index * max_index,
+                                    sizeof(*histo_heap->heap));
+  histo_heap->positions = WebPSafeMalloc(max_index * max_index,
+                                         sizeof(*histo_heap->heap));
+  return histo_heap->heap != NULL && histo_heap->positions != NULL;
+}
+
+static void HistoHeapClear(HistoHeap* const histo_heap) {
+  assert(histo_heap != NULL);
+  WebPSafeFree(histo_heap->heap);
+  WebPSafeFree(histo_heap->positions);
+}
+
+static void SwapHistogramPairs(HistogramPair *p1,
+                               HistogramPair *p2) {
+  const HistogramPair tmp = *p1;
+  *p1 = *p2;
+  *p2 = tmp;
+}
+
+// Given a valid min-heap in range [0, heap_size-1) this function places value
+// heap[heap_size-1] into right location within heap and sets its position in
+// positions array.
+static void HeapPush(HistoHeap* const histo_heap) {
+  HistogramPair* const heap = histo_heap->heap - 1;
+  int* const positions = histo_heap->positions;
+  const int max_index = histo_heap->max_index;
+  int v;
+  ++histo_heap->size;
+  v = histo_heap->size;
+  while (v > 1 && heap[v].cost_diff < heap[v >> 1].cost_diff) {
+    SwapHistogramPairs(&heap[v], &heap[v >> 1]);
+    // Change position of moved pair in heap.
+    if (heap[v].idx1 >= 0) {
+      const int pos = heap[v].idx1 * max_index + heap[v].idx2;
+      assert(pos >= 0 && pos < max_index * max_index);
+      positions[pos] = v;
+    }
+    v >>= 1;
+  }
+  positions[heap[v].idx1 * max_index + heap[v].idx2] = v;
+}
+
+// Given a valid min-heap in range [0, heap_size) this function shortens heap
+// range by one and places element with the lowest value to (heap_size-1).
+static void HeapPop(HistoHeap* const histo_heap) {
+  HistogramPair* const heap = histo_heap->heap - 1;
+  int* const positions = histo_heap->positions;
+  const int heap_size = histo_heap->size;
+  const int max_index = histo_heap->max_index;
+  int v = 1;
+  if (heap[v].idx1 >= 0) {
+    positions[heap[v].idx1 * max_index + heap[v].idx2] = -1;
+  }
+  SwapHistogramPairs(&heap[v], &heap[heap_size]);
+  while ((v << 1) < heap_size) {
+    int son = (heap[v << 1].cost_diff < heap[v].cost_diff) ? (v << 1) : v;
+    if (((v << 1) + 1) < heap_size &&
+        heap[(v << 1) + 1].cost_diff < heap[son].cost_diff) {
+      son = (v << 1) + 1;
+    }
+    if (son == v) break;
+    SwapHistogramPairs(&heap[v], &heap[son]);
+    // Change position of moved pair in heap.
+    if (heap[v].idx1 >= 0) {
+      positions[heap[v].idx1 * max_index + heap[v].idx2] = v;
+    }
+    v = son;
+  }
+  if (heap[v].idx1 >= 0) {
+    positions[heap[v].idx1 * max_index + heap[v].idx2] = v;
+  }
+  --histo_heap->size;
+}
+
+// -----------------------------------------------------------------------------
+
+static void PreparePair(VP8LHistogram** histograms, int idx1, int idx2,
+                        HistogramPair* const pair,
+                        VP8LHistogram* const histos) {
+  if (idx1 > idx2) {
+    const int tmp = idx2;
+    idx2 = idx1;
+    idx1 = tmp;
+  }
+  pair->idx1 = idx1;
+  pair->idx2 = idx2;
+  pair->cost_diff =
+      HistogramAddEval(histograms[idx1], histograms[idx2], histos, 0);
+  pair->cost_combo = histos->bit_cost_;
+}
+
+#define POSITION_INVALID (-1)
+
+// Invalidates pairs intersecting (idx1, idx2) in heap.
+static void InvalidatePairs(int idx1, int idx2,
+                            const HistoHeap* const histo_heap) {
+  HistogramPair* const heap = histo_heap->heap - 1;
+  int* const positions = histo_heap->positions;
+  const int max_index = histo_heap->max_index;
+  int i;
+  for (i = 0; i < idx1; ++i) {
+    const int pos = positions[i * max_index + idx1];
+    if (pos >= 0) {
+      heap[pos].idx1 = POSITION_INVALID;
+    }
+  }
+  for (i = idx1 + 1; i < max_index; ++i) {
+    const int pos = positions[idx1 * max_index + i];
+    if (pos >= 0) {
+      heap[pos].idx1 = POSITION_INVALID;
+    }
+  }
+  for (i = 0; i < idx2; ++i) {
+    const int pos = positions[i * max_index + idx2];
+    if (pos >= 0) {
+      heap[pos].idx1 = POSITION_INVALID;
+    }
+  }
+  for (i = idx2 + 1; i < max_index; ++i) {
+    const int pos = positions[idx2 * max_index + i];
+    if (pos >= 0) {
+      heap[pos].idx1 = POSITION_INVALID;
+    }
+  }
+}
+
+// Combines histograms by continuously choosing the one with the highest cost
+// reduction.
+static int HistogramCombineGreedy(VP8LHistogramSet* const image_histo,
+                                  VP8LHistogram* const histos) {
+  int ok = 0;
+  int image_histo_size = image_histo->size;
+  int i, j;
+  VP8LHistogram** const histograms = image_histo->histograms;
+  // Indexes of remaining histograms.
+  int* const clusters = WebPSafeMalloc(image_histo_size, sizeof(*clusters));
+  // Heap of histogram pairs.
+  HistoHeap histo_heap;
+
+  if (!HistoHeapInit(&histo_heap, image_histo_size) || clusters == NULL) {
+    goto End;
+  }
+
+  for (i = 0; i < image_histo_size; ++i) {
+    // Initialize clusters indexes.
+    clusters[i] = i;
+    for (j = i + 1; j < image_histo_size; ++j) {
+      // Initialize positions array.
+      histo_heap.positions[i * histo_heap.max_index + j] = POSITION_INVALID;
+      PreparePair(histograms, i, j, &histo_heap.heap[histo_heap.size], histos);
+      if (histo_heap.heap[histo_heap.size].cost_diff < 0) {
+        HeapPush(&histo_heap);
+      }
+    }
+  }
+
+  while (image_histo_size > 1 && histo_heap.size > 0) {
+    const int idx1 = histo_heap.heap[0].idx1;
+    const int idx2 = histo_heap.heap[0].idx2;
+    VP8LHistogramAdd(histograms[idx2], histograms[idx1], histograms[idx1]);
+    histograms[idx1]->bit_cost_ = histo_heap.heap[0].cost_combo;
+    // Remove merged histogram.
+    for (i = 0; i + 1 < image_histo_size; ++i) {
+      if (clusters[i] >= idx2) {
+        clusters[i] = clusters[i + 1];
+      }
+    }
+    --image_histo_size;
+
+    // Invalidate pairs intersecting the just combined best pair.
+    InvalidatePairs(idx1, idx2, &histo_heap);
+
+    // Pop invalid pairs from the top of the heap.
+    while (histo_heap.size > 0 && histo_heap.heap[0].idx1 < 0) {
+      HeapPop(&histo_heap);
+    }
+
+    // Push new pairs formed with combined histogram to the heap.
+    for (i = 0; i < image_histo_size; ++i) {
+      if (clusters[i] != idx1) {
+        PreparePair(histograms, idx1, clusters[i],
+                    &histo_heap.heap[histo_heap.size], histos);
+        if (histo_heap.heap[histo_heap.size].cost_diff < 0) {
+          HeapPush(&histo_heap);
+        }
+      }
+    }
+  }
+  // Move remaining histograms to the beginning of the array.
+  for (i = 0; i < image_histo_size; ++i) {
+    if (i != clusters[i]) {
+      HistogramCopy(histograms[clusters[i]], histograms[i]);
+    }
+  }
+
+  image_histo->size = image_histo_size;
+  ok = 1;
+
+ End:
+  WebPSafeFree(clusters);
+  HistoHeapClear(&histo_heap);
+  return ok;
+}
+
+static void HistogramCombineStochastic(VP8LHistogramSet* const image_histo,
+                                       VP8LHistogramSet* const histos,
+                                       int quality, int min_cluster_size) {
   int iter;
   uint32_t seed = 0;
   int tries_with_no_success = 0;
@@ -570,12 +799,12 @@ static void HistogramCombine(VP8LHistogramSet* const image_histo,
   const int outer_iters = image_histo_size * iter_mult;
   const int num_pairs = image_histo_size / 2;
   const int num_tries_no_success = outer_iters / 2;
-  const int min_cluster_size = 2;
   VP8LHistogram** const histograms = image_histo->histograms;
   VP8LHistogram* cur_combo = histos->histograms[0];   // trial histogram
   VP8LHistogram* best_combo = histos->histograms[1];  // best histogram so far
 
   // Collapse similar histograms in 'image_histo'.
+  ++min_cluster_size;
   for (iter = 0;
        iter < outer_iters && image_histo_size >= min_cluster_size;
        ++iter) {
@@ -734,8 +963,15 @@ int VP8LGetHistoImageSymbols(int xsize, int ysize,
                                bin_map, bin_depth, combine_cost_factor);
   }
 
-  // Collapse similar histograms by random histogram-pair compares.
-  HistogramCombine(image_histo, histos, quality);
+  {
+    const float x = quality / 100.f;
+    const int threshold_size = 1 + (x * x * x) * (MAX_HISTO_GREEDY - 1);
+    HistogramCombineStochastic(image_histo, histos, quality, threshold_size);
+    if ((image_histo->size <= threshold_size) &&
+        !HistogramCombineGreedy(image_histo, histos->histograms[0])) {
+      goto Error;
+    }
+  }
 
   // Find the optimal map from original histograms to the final ones.
   HistogramRemap(orig_histo, image_histo, histogram_symbols);