src/SDL2/external/libwebp-0.3.0/src/enc/histogram.c

   1 // Copyright 2012 Google Inc. All Rights Reserved.
   2 //
   3 // This code is licensed under the same terms as WebM:
   4 //  Software License Agreement:  http://www.webmproject.org/license/software/
   5 //  Additional IP Rights Grant:  http://www.webmproject.org/license/additional/
   6 // -----------------------------------------------------------------------------
   7 //
   8 // Author: Jyrki Alakuijala (jyrki@google.com)
   9 //
  10 #ifdef HAVE_CONFIG_H
  11 #include "config.h"
  12 #endif
  13
  14 #include <math.h>
  15 #include <stdio.h>
  16
  17 #include "./backward_references.h"
  18 #include "./histogram.h"
  19 #include "../dsp/lossless.h"
  20 #include "../utils/utils.h"
  21
  22 static void HistogramClear(VP8LHistogram* const p) {
  23   memset(p->literal_, 0, sizeof(p->literal_));
  24   memset(p->red_, 0, sizeof(p->red_));
  25   memset(p->blue_, 0, sizeof(p->blue_));
  26   memset(p->alpha_, 0, sizeof(p->alpha_));
  27   memset(p->distance_, 0, sizeof(p->distance_));
  28   p->bit_cost_ = 0;
  29 }
  30
  31 void VP8LHistogramStoreRefs(const VP8LBackwardRefs* const refs,
  32                             VP8LHistogram* const histo) {
  33   int i;
  34   for (i = 0; i < refs->size; ++i) {
  35     VP8LHistogramAddSinglePixOrCopy(histo, &refs->refs[i]);
  36   }
  37 }
  38
  39 void VP8LHistogramCreate(VP8LHistogram* const p,
  40                          const VP8LBackwardRefs* const refs,
  41                          int palette_code_bits) {
  42   if (palette_code_bits >= 0) {
  43     p->palette_code_bits_ = palette_code_bits;
  44   }
  45   HistogramClear(p);
  46   VP8LHistogramStoreRefs(refs, p);
  47 }
  48
  49 void VP8LHistogramInit(VP8LHistogram* const p, int palette_code_bits) {
  50   p->palette_code_bits_ = palette_code_bits;
  51   HistogramClear(p);
  52 }
  53
  54 VP8LHistogramSet* VP8LAllocateHistogramSet(int size, int cache_bits) {
  55   int i;
  56   VP8LHistogramSet* set;
  57   VP8LHistogram* bulk;
  58   const uint64_t total_size = sizeof(*set)
  59                             + (uint64_t)size * sizeof(*set->histograms)
  60                             + (uint64_t)size * sizeof(**set->histograms);
  61   uint8_t* memory = (uint8_t*)WebPSafeMalloc(total_size, sizeof(*memory));
  62   if (memory == NULL) return NULL;
  63
  64   set = (VP8LHistogramSet*)memory;
  65   memory += sizeof(*set);
  66   set->histograms = (VP8LHistogram**)memory;
  67   memory += size * sizeof(*set->histograms);
  68   bulk = (VP8LHistogram*)memory;
  69   set->max_size = size;
  70   set->size = size;
  71   for (i = 0; i < size; ++i) {
  72     set->histograms[i] = bulk + i;
  73     VP8LHistogramInit(set->histograms[i], cache_bits);
  74   }
  75   return set;
  76 }
  77
  78 // -----------------------------------------------------------------------------
  79
  80 void VP8LHistogramAddSinglePixOrCopy(VP8LHistogram* const histo,
  81                                      const PixOrCopy* const v) {
  82   if (PixOrCopyIsLiteral(v)) {
  83     ++histo->alpha_[PixOrCopyLiteral(v, 3)];
  84     ++histo->red_[PixOrCopyLiteral(v, 2)];
  85     ++histo->literal_[PixOrCopyLiteral(v, 1)];
  86     ++histo->blue_[PixOrCopyLiteral(v, 0)];
  87   } else if (PixOrCopyIsCacheIdx(v)) {
  88     int literal_ix = 256 + NUM_LENGTH_CODES + PixOrCopyCacheIdx(v);
  89     ++histo->literal_[literal_ix];
  90   } else {
  91     int code, extra_bits_count, extra_bits_value;
  92     PrefixEncode(PixOrCopyLength(v),
  93                  &code, &extra_bits_count, &extra_bits_value);
  94     ++histo->literal_[256 + code];
  95     PrefixEncode(PixOrCopyDistance(v),
  96                  &code, &extra_bits_count, &extra_bits_value);
  97     ++histo->distance_[code];
  98   }
  99 }
 100
 101 static double BitsEntropy(const int* const array, int n) {
 102   double retval = 0.;
 103   int sum = 0;
 104   int nonzeros = 0;
 105   int max_val = 0;
 106   int i;
 107   double mix;
 108   for (i = 0; i < n; ++i) {
 109     if (array[i] != 0) {
 110       sum += array[i];
 111       ++nonzeros;
 112       retval -= VP8LFastSLog2(array[i]);
 113       if (max_val < array[i]) {
 114         max_val = array[i];
 115       }
 116     }
 117   }
 118   retval += VP8LFastSLog2(sum);
 119
 120   if (nonzeros < 5) {
 121     if (nonzeros <= 1) {
 122       return 0;
 123     }
 124     // Two symbols, they will be 0 and 1 in a Huffman code.
 125     // Let's mix in a bit of entropy to favor good clustering when
 126     // distributions of these are combined.
 127     if (nonzeros == 2) {
 128       return 0.99 * sum + 0.01 * retval;
 129     }
 130     // No matter what the entropy says, we cannot be better than min_limit
 131     // with Huffman coding. I am mixing a bit of entropy into the
 132     // min_limit since it produces much better (~0.5 %) compression results
 133     // perhaps because of better entropy clustering.
 134     if (nonzeros == 3) {
 135       mix = 0.95;
 136     } else {
 137       mix = 0.7;  // nonzeros == 4.
 138     }
 139   } else {
 140     mix = 0.627;
 141   }
 142
 143   {
 144     double min_limit = 2 * sum - max_val;
 145     min_limit = mix * min_limit + (1.0 - mix) * retval;
 146     return (retval < min_limit) ? min_limit : retval;
 147   }
 148 }
 149
 150 // Returns the cost encode the rle-encoded entropy code.
 151 // The constants in this function are experimental.
 152 static double HuffmanCost(const int* const population, int length) {
 153   // Small bias because Huffman code length is typically not stored in
 154   // full length.
 155   static const int kHuffmanCodeOfHuffmanCodeSize = CODE_LENGTH_CODES * 3;
 156   static const double kSmallBias = 9.1;
 157   double retval = kHuffmanCodeOfHuffmanCodeSize - kSmallBias;
 158   int streak = 0;
 159   int i = 0;
 160   for (; i < length - 1; ++i) {
 161     ++streak;
 162     if (population[i] == population[i + 1]) {
 163       continue;
 164     }
 165  last_streak_hack:
 166     // population[i] points now to the symbol in the streak of same values.
 167     if (streak > 3) {
 168       if (population[i] == 0) {
 169         retval += 1.5625 + 0.234375 * streak;
 170       } else {
 171         retval += 2.578125 + 0.703125 * streak;
 172       }
 173     } else {
 174       if (population[i] == 0) {
 175         retval += 1.796875 * streak;
 176       } else {
 177         retval += 3.28125 * streak;
 178       }
 179     }
 180     streak = 0;
 181   }
 182   if (i == length - 1) {
 183     ++streak;
 184     goto last_streak_hack;
 185   }
 186   return retval;
 187 }
 188
 189 static double PopulationCost(const int* const population, int length) {
 190   return BitsEntropy(population, length) + HuffmanCost(population, length);
 191 }
 192
 193 static double ExtraCost(const int* const population, int length) {
 194   int i;
 195   double cost = 0.;
 196   for (i = 2; i < length - 2; ++i) cost += (i >> 1) * population[i + 2];
 197   return cost;
 198 }
 199
 200 // Estimates the Entropy + Huffman + other block overhead size cost.
 201 double VP8LHistogramEstimateBits(const VP8LHistogram* const p) {
 202   return PopulationCost(p->literal_, VP8LHistogramNumCodes(p))
 203        + PopulationCost(p->red_, 256)
 204        + PopulationCost(p->blue_, 256)
 205        + PopulationCost(p->alpha_, 256)
 206        + PopulationCost(p->distance_, NUM_DISTANCE_CODES)
 207        + ExtraCost(p->literal_ + 256, NUM_LENGTH_CODES)
 208        + ExtraCost(p->distance_, NUM_DISTANCE_CODES);
 209 }
 210
 211 double VP8LHistogramEstimateBitsBulk(const VP8LHistogram* const p) {
 212   return BitsEntropy(p->literal_, VP8LHistogramNumCodes(p))
 213        + BitsEntropy(p->red_, 256)
 214        + BitsEntropy(p->blue_, 256)
 215        + BitsEntropy(p->alpha_, 256)
 216        + BitsEntropy(p->distance_, NUM_DISTANCE_CODES)
 217        + ExtraCost(p->literal_ + 256, NUM_LENGTH_CODES)
 218        + ExtraCost(p->distance_, NUM_DISTANCE_CODES);
 219 }
 220
 221 // -----------------------------------------------------------------------------
 222 // Various histogram combine/cost-eval functions
 223
 224 // Adds 'in' histogram to 'out'
 225 static void HistogramAdd(const VP8LHistogram* const in,
 226                          VP8LHistogram* const out) {
 227   int i;
 228   for (i = 0; i < PIX_OR_COPY_CODES_MAX; ++i) {
 229     out->literal_[i] += in->literal_[i];
 230   }
 231   for (i = 0; i < NUM_DISTANCE_CODES; ++i) {
 232     out->distance_[i] += in->distance_[i];
 233   }
 234   for (i = 0; i < 256; ++i) {
 235     out->red_[i] += in->red_[i];
 236     out->blue_[i] += in->blue_[i];
 237     out->alpha_[i] += in->alpha_[i];
 238   }
 239 }
 240
 241 // Performs out = a + b, computing the cost C(a+b) - C(a) - C(b) while comparing
 242 // to the threshold value 'cost_threshold'. The score returned is
 243 //  Score = C(a+b) - C(a) - C(b), where C(a) + C(b) is known and fixed.
 244 // Since the previous score passed is 'cost_threshold', we only need to compare
 245 // the partial cost against 'cost_threshold + C(a) + C(b)' to possibly bail-out
 246 // early.
 247 static double HistogramAddEval(const VP8LHistogram* const a,
 248                                const VP8LHistogram* const b,
 249                                VP8LHistogram* const out,
 250                                double cost_threshold) {
 251   double cost = 0;
 252   const double sum_cost = a->bit_cost_ + b->bit_cost_;
 253   int i;
 254
 255   cost_threshold += sum_cost;
 256
 257   // palette_code_bits_ is part of the cost evaluation for literal_.
 258   // TODO(skal): remove/simplify this palette_code_bits_?
 259   out->palette_code_bits_ =
 260       (a->palette_code_bits_ > b->palette_code_bits_) ? a->palette_code_bits_ :
 261                                                         b->palette_code_bits_;
 262   for (i = 0; i < PIX_OR_COPY_CODES_MAX; ++i) {
 263     out->literal_[i] = a->literal_[i] + b->literal_[i];
 264   }
 265   cost += PopulationCost(out->literal_, VP8LHistogramNumCodes(out));
 266   cost += ExtraCost(out->literal_ + 256, NUM_LENGTH_CODES);
 267   if (cost > cost_threshold) return cost;
 268
 269   for (i = 0; i < 256; ++i) out->red_[i] = a->red_[i] + b->red_[i];
 270   cost += PopulationCost(out->red_, 256);
 271   if (cost > cost_threshold) return cost;
 272
 273   for (i = 0; i < 256; ++i) out->blue_[i] = a->blue_[i] + b->blue_[i];
 274   cost += PopulationCost(out->blue_, 256);
 275   if (cost > cost_threshold) return cost;
 276
 277   for (i = 0; i < NUM_DISTANCE_CODES; ++i) {
 278     out->distance_[i] = a->distance_[i] + b->distance_[i];
 279   }
 280   cost += PopulationCost(out->distance_, NUM_DISTANCE_CODES);
 281   cost += ExtraCost(out->distance_, NUM_DISTANCE_CODES);
 282   if (cost > cost_threshold) return cost;
 283
 284   for (i = 0; i < 256; ++i) out->alpha_[i] = a->alpha_[i] + b->alpha_[i];
 285   cost += PopulationCost(out->alpha_, 256);
 286
 287   out->bit_cost_ = cost;
 288   return cost - sum_cost;
 289 }
 290
 291 // Same as HistogramAddEval(), except that the resulting histogram
 292 // is not stored. Only the cost C(a+b) - C(a) is evaluated. We omit
 293 // the term C(b) which is constant over all the evaluations.
 294 static double HistogramAddThresh(const VP8LHistogram* const a,
 295                                  const VP8LHistogram* const b,
 296                                  double cost_threshold) {
 297   int tmp[PIX_OR_COPY_CODES_MAX];  // <= max storage we'll need
 298   int i;
 299   double cost = -a->bit_cost_;
 300
 301   for (i = 0; i < PIX_OR_COPY_CODES_MAX; ++i) {
 302     tmp[i] = a->literal_[i] + b->literal_[i];
 303   }
 304   // note that the tests are ordered so that the usually largest
 305   // cost shares come first.
 306   cost += PopulationCost(tmp, VP8LHistogramNumCodes(a));
 307   cost += ExtraCost(tmp + 256, NUM_LENGTH_CODES);
 308   if (cost > cost_threshold) return cost;
 309
 310   for (i = 0; i < 256; ++i) tmp[i] = a->red_[i] + b->red_[i];
 311   cost += PopulationCost(tmp, 256);
 312   if (cost > cost_threshold) return cost;
 313
 314   for (i = 0; i < 256; ++i) tmp[i] = a->blue_[i] + b->blue_[i];
 315   cost += PopulationCost(tmp, 256);
 316   if (cost > cost_threshold) return cost;
 317
 318   for (i = 0; i < NUM_DISTANCE_CODES; ++i) {
 319     tmp[i] = a->distance_[i] + b->distance_[i];
 320   }
 321   cost += PopulationCost(tmp, NUM_DISTANCE_CODES);
 322   cost += ExtraCost(tmp, NUM_DISTANCE_CODES);
 323   if (cost > cost_threshold) return cost;
 324
 325   for (i = 0; i < 256; ++i) tmp[i] = a->alpha_[i] + b->alpha_[i];
 326   cost += PopulationCost(tmp, 256);
 327
 328   return cost;
 329 }
 330
 331 // -----------------------------------------------------------------------------
 332
 333 static void HistogramBuildImage(int xsize, int histo_bits,
 334                                 const VP8LBackwardRefs* const backward_refs,
 335                                 VP8LHistogramSet* const image) {
 336   int i;
 337   int x = 0, y = 0;
 338   const int histo_xsize = VP8LSubSampleSize(xsize, histo_bits);
 339   VP8LHistogram** const histograms = image->histograms;
 340   assert(histo_bits > 0);
 341   for (i = 0; i < backward_refs->size; ++i) {
 342     const PixOrCopy* const v = &backward_refs->refs[i];
 343     const int ix = (y >> histo_bits) * histo_xsize + (x >> histo_bits);
 344     VP8LHistogramAddSinglePixOrCopy(histograms[ix], v);
 345     x += PixOrCopyLength(v);
 346     while (x >= xsize) {
 347       x -= xsize;
 348       ++y;
 349     }
 350   }
 351 }
 352
 353 static uint32_t MyRand(uint32_t *seed) {
 354   *seed *= 16807U;
 355   if (*seed == 0) {
 356     *seed = 1;
 357   }
 358   return *seed;
 359 }
 360
 361 static int HistogramCombine(const VP8LHistogramSet* const in,
 362                             VP8LHistogramSet* const out, int iter_mult,
 363                             int num_pairs, int num_tries_no_success) {
 364   int ok = 0;
 365   int i, iter;
 366   uint32_t seed = 0;
 367   int tries_with_no_success = 0;
 368   int out_size = in->size;
 369   const int outer_iters = in->size * iter_mult;
 370   const int min_cluster_size = 2;
 371   VP8LHistogram* const histos = (VP8LHistogram*)malloc(2 * sizeof(*histos));
 372   VP8LHistogram* cur_combo = histos + 0;    // trial merged histogram
 373   VP8LHistogram* best_combo = histos + 1;   // best merged histogram so far
 374   if (histos == NULL) goto End;
 375
 376   // Copy histograms from in[] to out[].
 377   assert(in->size <= out->size);
 378   for (i = 0; i < in->size; ++i) {
 379     in->histograms[i]->bit_cost_ = VP8LHistogramEstimateBits(in->histograms[i]);
 380     *out->histograms[i] = *in->histograms[i];
 381   }
 382
 383   // Collapse similar histograms in 'out'.
 384   for (iter = 0; iter < outer_iters && out_size >= min_cluster_size; ++iter) {
 385     double best_cost_diff = 0.;
 386     int best_idx1 = -1, best_idx2 = 1;
 387     int j;
 388     const int num_tries = (num_pairs < out_size) ? num_pairs : out_size;
 389     seed += iter;
 390     for (j = 0; j < num_tries; ++j) {
 391       double curr_cost_diff;
 392       // Choose two histograms at random and try to combine them.
 393       const uint32_t idx1 = MyRand(&seed) % out_size;
 394       const uint32_t tmp = (j & 7) + 1;
 395       const uint32_t diff = (tmp < 3) ? tmp : MyRand(&seed) % (out_size - 1);
 396       const uint32_t idx2 = (idx1 + diff + 1) % out_size;
 397       if (idx1 == idx2) {
 398         continue;
 399       }
 400       // Calculate cost reduction on combining.
 401       curr_cost_diff = HistogramAddEval(out->histograms[idx1],
 402                                         out->histograms[idx2],
 403                                         cur_combo, best_cost_diff);
 404       if (curr_cost_diff < best_cost_diff) {    // found a better pair?
 405         {     // swap cur/best combo histograms
 406           VP8LHistogram* const tmp_histo = cur_combo;
 407           cur_combo = best_combo;
 408           best_combo = tmp_histo;
 409         }
 410         best_cost_diff = curr_cost_diff;
 411         best_idx1 = idx1;
 412         best_idx2 = idx2;
 413       }
 414     }
 415
 416     if (best_idx1 >= 0) {
 417       *out->histograms[best_idx1] = *best_combo;
 418       // swap best_idx2 slot with last one (which is now unused)
 419       --out_size;
 420       if (best_idx2 != out_size) {
 421         out->histograms[best_idx2] = out->histograms[out_size];
 422         out->histograms[out_size] = NULL;   // just for sanity check.
 423       }
 424       tries_with_no_success = 0;
 425     }
 426     if (++tries_with_no_success >= num_tries_no_success) {
 427       break;
 428     }
 429   }
 430   out->size = out_size;
 431   ok = 1;
 432
 433  End:
 434   free(histos);
 435   return ok;
 436 }
 437
 438 // -----------------------------------------------------------------------------
 439 // Histogram refinement
 440
 441 // What is the bit cost of moving square_histogram from cur_symbol to candidate.
 442 static double HistogramDistance(const VP8LHistogram* const square_histogram,
 443                                 const VP8LHistogram* const candidate,
 444                                 double cost_threshold) {
 445   return HistogramAddThresh(candidate, square_histogram, cost_threshold);
 446 }
 447
 448 // Find the best 'out' histogram for each of the 'in' histograms.
 449 // Note: we assume that out[]->bit_cost_ is already up-to-date.
 450 static void HistogramRemap(const VP8LHistogramSet* const in,
 451                            const VP8LHistogramSet* const out,
 452                            uint16_t* const symbols) {
 453   int i;
 454   for (i = 0; i < in->size; ++i) {
 455     int best_out = 0;
 456     double best_bits =
 457         HistogramDistance(in->histograms[i], out->histograms[0], 1.e38);
 458     int k;
 459     for (k = 1; k < out->size; ++k) {
 460       const double cur_bits =
 461           HistogramDistance(in->histograms[i], out->histograms[k], best_bits);
 462       if (cur_bits < best_bits) {
 463         best_bits = cur_bits;
 464         best_out = k;
 465       }
 466     }
 467     symbols[i] = best_out;
 468   }
 469
 470   // Recompute each out based on raw and symbols.
 471   for (i = 0; i < out->size; ++i) {
 472     HistogramClear(out->histograms[i]);
 473   }
 474   for (i = 0; i < in->size; ++i) {
 475     HistogramAdd(in->histograms[i], out->histograms[symbols[i]]);
 476   }
 477 }
 478
 479 int VP8LGetHistoImageSymbols(int xsize, int ysize,
 480                              const VP8LBackwardRefs* const refs,
 481                              int quality, int histo_bits, int cache_bits,
 482                              VP8LHistogramSet* const image_in,
 483                              uint16_t* const histogram_symbols) {
 484   int ok = 0;
 485   const int histo_xsize = histo_bits ? VP8LSubSampleSize(xsize, histo_bits) : 1;
 486   const int histo_ysize = histo_bits ? VP8LSubSampleSize(ysize, histo_bits) : 1;
 487   const int histo_image_raw_size = histo_xsize * histo_ysize;
 488
 489   // Heuristic params for HistogramCombine().
 490   const int num_tries_no_success = 8 + (quality >> 1);
 491   const int iter_mult = (quality < 27) ? 1 : 1 + ((quality - 27) >> 4);
 492   const int num_pairs = (quality < 25) ? 10 : (5 * quality) >> 3;
 493
 494   VP8LHistogramSet* const image_out =
 495       VP8LAllocateHistogramSet(histo_image_raw_size, cache_bits);
 496   if (image_out == NULL) return 0;
 497
 498   // Build histogram image.
 499   HistogramBuildImage(xsize, histo_bits, refs, image_out);
 500   // Collapse similar histograms.
 501   if (!HistogramCombine(image_out, image_in, iter_mult, num_pairs,
 502                         num_tries_no_success)) {
 503     goto Error;
 504   }
 505   // Find the optimal map from original histograms to the final ones.
 506   HistogramRemap(image_out, image_in, histogram_symbols);
 507   ok = 1;
 508
 509 Error:
 510   free(image_out);
 511   return ok;
 512 }