histogram.c revision 2a99a7e74a7f215066514fe81d2bfa6639d9eddd
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 22static 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 31void 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 39void 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 49void VP8LHistogramInit(VP8LHistogram* const p, int palette_code_bits) { 50 p->palette_code_bits_ = palette_code_bits; 51 HistogramClear(p); 52} 53 54VP8LHistogramSet* 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 80void 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 101static 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. 152static 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 189static double PopulationCost(const int* const population, int length) { 190 return BitsEntropy(population, length) + HuffmanCost(population, length); 191} 192 193static 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. 201double 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 211double 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' 225static 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. 247static 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. 294static 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 333static 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 353static uint32_t MyRand(uint32_t *seed) { 354 *seed *= 16807U; 355 if (*seed == 0) { 356 *seed = 1; 357 } 358 return *seed; 359} 360 361static 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. 442static 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. 450static 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 479int 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 509Error: 510 free(image_out); 511 return ok; 512} 513