#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
+#include <string.h>
#include <assert.h>
#include <math.h>
#include "ryg_rans/renormalize.h"
#include <memory>
+#include <vector>
#define WIDTH 1280
#define HEIGHT 720
uint32_t cum_freqs[NUM_SYMS + 1];
void clear();
- void count_freqs(uint8_t const* in, size_t nbytes);
void calc_cum_freqs();
void normalize_freqs(uint32_t target_total);
};
freqs[i] = 0;
}
-void SymbolStats::count_freqs(uint8_t const* in, size_t nbytes)
-{
- clear();
-
- for (size_t i=0; i < nbytes; i++)
- freqs[in[i]]++;
-}
-
void SymbolStats::calc_cum_freqs()
{
cum_freqs[0] = 0;
cum_freqs[i+1] = cum_freqs[i] + freqs[i];
}
-static double cache[NUM_SYMS + 1][prob_scale + 1];
-static double log2cache[prob_scale + 1];
-static int64_t cachefill = 0;
-
-double find_optimal_cost(const uint32_t *cum_freqs, int sym_to, int available_slots)
-{
- assert(sym_to >= 0);
-
- while (sym_to > 0 && cum_freqs[sym_to] == cum_freqs[sym_to - 1]) { --sym_to; }
- if (cache[sym_to][available_slots] >= 0.0) {
- //printf("CACHE: %d,%d\n", sym_to, available_slots);
- return cache[sym_to][available_slots];
- }
- if (sym_to == 0) {
- return 0.0;
- }
- if (sym_to == 1) {
- return cum_freqs[0] * log2cache[available_slots];
- }
- if (available_slots == 1) {
- return cum_freqs[0] * log2cache[1] + find_optimal_cost(cum_freqs, sym_to - 1, 0);
- }
-
-// printf("UNCACHE: %d,%d\n", sym_to, available_slots);
-#if 0
- // ok, test all possible options for the last symbol (TODO: save the choice)
- double best_so_far = HUGE_VAL;
- //for (int i = num_syms - 1; i < available_slots; ++i) {
- double f = freqs[sym_to - 1];
- for (int i = available_slots; i --> 0; ) {
- double cost1 = f * log2cache[available_slots - i];
- double cost2 = find_optimal_cost(freqs, sym_to - 1, i);
-
- if (sym_to == 3 && available_slots == 838) {
- printf("%d %f\n", i, cost1 + cost2);
- } else
- if (cost1 + cost2 > best_so_far) {
- break;
- }
- best_so_far = cost1 + cost2;
- }
-#elif 1
- // Minimize the number of total bits spent as a function of how many slots
- // we assign to this symbol.
- //
- // The cost function is convex (I don't know how to prove it, but it makes
- // intuitively a lot of sense). Find a reasonable guess and see what way
- // we should search, then iterate until we either hit the end or we start
- // increasing again.
- double f = cum_freqs[sym_to - 1] - cum_freqs[sym_to - 2];
- double start = lrint(available_slots * f / cum_freqs[sym_to - 1]);
-
- int x1 = std::max<int>(floor(start), 1);
- int x2 = x1 + 1;
-
- double f1 = f * log2cache[x1] + find_optimal_cost(cum_freqs, sym_to - 1, available_slots - x1);
- double f2 = f * log2cache[x2] + find_optimal_cost(cum_freqs, sym_to - 1, available_slots - x2);
-
- int x, direction; // -1 or +1
- double best_so_far = std::min(f1, f2);
- if (isinf(f1) && isinf(f2)) {
- // The cost isn't infinite due to the first term, so we need to go downwards
- // to give the second term more room to breathe.
- x = x1;
- direction = -1;
- } else if (f1 < f2) {
- x = x1;
- direction = -1;
- } else {
- x = x2;
- direction = 1;
- }
-
- //printf("[%d,%d] freq=%ld cumfreq=%d From %d and %d, chose %d [%f] and direction=%d\n",
- // sym_to, available_slots, freqs[sym_to - 1], cum_freqs[sym_to - 1], x1, x2, x, best_so_far, direction);
-
- while ((x + direction) > 0 && (x + direction) <= available_slots) {
- x += direction;
- double fn = f * log2cache[x] + find_optimal_cost(cum_freqs, sym_to - 1, available_slots - x);
- // printf("[%d,%d] %d is %f\n", sym_to, available_slots, x, fn);
- if (fn > best_so_far) {
- break;
- }
- best_so_far = fn;
- }
-#endif
- if (++cachefill % 131072 == 0) {
- // printf("%d,%d = %f (cachefill = %.2f%%)\n", sym_to, available_slots, best_so_far,
- // 100.0 * (cachefill / double((NUM_SYMS + 1) * (prob_scale + 1))));
- }
- assert(best_so_far >= 0.0);
- assert(cache[sym_to][available_slots] < 0.0);
- cache[sym_to][available_slots] = best_so_far;
- return best_so_far;
-}
-
-double find_optimal_cost(const uint32_t *cum_freqs, const uint64_t *freqs)
-{
- for (int j = 0; j <= NUM_SYMS; ++j) {
- for (unsigned k = 0; k <= prob_scale; ++k) {
- cache[j][k] = -1.0;
- }
- }
- for (unsigned k = 0; k <= prob_scale; ++k) {
- log2cache[k] = -log2(k * (1.0 / prob_scale));
- //printf("log2cache[%d] = %f\n", k, log2cache[k]);
- }
- cachefill = 0;
- double ret = find_optimal_cost(cum_freqs, NUM_SYMS, prob_scale);
- printf("Used %ld function invocations\n", cachefill);
- return ret;
-}
-
void SymbolStats::normalize_freqs(uint32_t target_total)
{
uint64_t real_freq[NUM_SYMS + 1]; // hack
OptimalRenormalize(cum_freqs, NUM_SYMS, prob_scale);
-#if 0
- double optimal_cost = find_optimal_cost(cum_freqs + 1, real_freq + 1);
-
- // resample distribution based on cumulative freqs
- for (int i = 1; i <= NUM_SYMS; i++)
- //cum_freqs[i] = ((uint64_t)target_total * cum_freqs[i])/cur_total;
- cum_freqs[i] = lrint(cum_freqs[i] * double(target_total) / cur_total);
-
- // if we nuked any non-0 frequency symbol to 0, we need to steal
- // the range to make the frequency nonzero from elsewhere.
- //
- // this is not at all optimal, i'm just doing the first thing that comes to mind.
- for (int i=0; i < NUM_SYMS; i++) {
- if (freqs[i] && cum_freqs[i+1] == cum_freqs[i]) {
- // symbol i was set to zero freq
-
- // find best symbol to steal frequency from (try to steal from low-freq ones)
- uint32_t best_freq = ~0u;
- int best_steal = -1;
- for (int j=0; j < NUM_SYMS; j++) {
- uint32_t freq = cum_freqs[j+1] - cum_freqs[j];
- if (freq > 1 && freq < best_freq) {
- best_freq = freq;
- best_steal = j;
- }
- }
- assert(best_steal != -1);
-
- // and steal from it!
- if (best_steal < i) {
- for (int j = best_steal + 1; j <= i; j++)
- cum_freqs[j]--;
- } else {
- assert(best_steal > i);
- for (int j = i + 1; j <= best_steal; j++)
- cum_freqs[j]++;
- }
- }
- }
-#endif
-
// calculate updated freqs and make sure we didn't screw anything up
assert(cum_freqs[0] == 0 && cum_freqs[NUM_SYMS] == target_total);
for (int i=0; i < NUM_SYMS; i++) {
int pick_stats_for(int y, int x)
{
+ //return 0;
//return std::min<int>(hypot(x, y), 7);
return std::min<int>(x + y, 7);
//if (x + y >= 7) return 7;
RansEncoder()
{
out_buf.reset(new uint8_t[out_max_size]);
- sign_buf.reset(new uint8_t[max_num_sign]);
clear();
}
- void init_prob(const SymbolStats &s)
+ void init_prob(SymbolStats &s)
{
for (int i = 0; i < NUM_SYMS; i++) {
- printf("%d: cumfreqs=%d freqs=%d prob_bits=%d\n", i, s.cum_freqs[i], s.freqs[i], prob_bits);
- RansEncSymbolInit(&esyms[i], s.cum_freqs[i], s.freqs[i], prob_bits);
+ printf("%d: cumfreqs=%d freqs=%d prob_bits=%d\n", i, s.cum_freqs[i], s.freqs[i], prob_bits + 1);
+ RansEncSymbolInit(&esyms[i], s.cum_freqs[i], s.freqs[i], prob_bits + 1);
}
+ sign_bias = s.cum_freqs[256];
}
void clear()
{
out_end = out_buf.get() + out_max_size;
- sign_end = sign_buf.get() + max_num_sign;
ptr = out_end; // *end* of output buffer
- sign_ptr = sign_end; // *end* of output buffer
RansEncInit(&rans);
- free_sign_bits = 0;
}
uint32_t save_block(FILE *codedfp) // Returns number of bytes.
//printf("post-flush = %08x\n", rans);
uint32_t num_rans_bytes = out_end - ptr;
+#if 0
+ if (num_rans_bytes == 4) {
+ uint32_t block;
+ memcpy(&block, ptr, 4);
+
+ if (block == last_block) {
+ write_varint(0, codedfp);
+ clear();
+ return 1;
+ }
+
+ last_block = block;
+ } else {
+ last_block = 0;
+ }
+#endif
+
write_varint(num_rans_bytes, codedfp);
//fwrite(&num_rans_bytes, 1, 4, codedfp);
fwrite(ptr, 1, num_rans_bytes, codedfp);
//printf("first rANS bytes: %02x %02x %02x %02x %02x %02x %02x %02x\n", ptr[0], ptr[1], ptr[2], ptr[3], ptr[4], ptr[5], ptr[6], ptr[7]);
- if (free_sign_bits > 0) {
- *sign_ptr <<= free_sign_bits;
- }
-
-#if 1
- uint32_t num_sign_bytes = sign_end - sign_ptr;
- write_varint((num_sign_bytes << 3) | free_sign_bits, codedfp);
- fwrite(sign_ptr, 1, num_sign_bytes, codedfp);
-#endif
clear();
- printf("Saving block: %d rANS bytes, %d sign bytes\n", num_rans_bytes, num_sign_bytes);
- return num_rans_bytes + num_sign_bytes;
+ printf("Saving block: %d rANS bytes\n", num_rans_bytes);
+ return num_rans_bytes;
//return num_rans_bytes;
}
void encode_coeff(short signed_k)
{
- printf("encoding coeff %d\n", signed_k);
- short k = abs(signed_k);
+ //printf("encoding coeff %d (sym %d), rans before encoding = %08x\n", signed_k, ((abs(signed_k) - 1) & 255), rans);
+ unsigned short k = abs(signed_k);
if (k >= ESCAPE_LIMIT) {
// Put the coefficient as a 1/(2^12) symbol _before_
// the 255 coefficient, since the decoder will read the
RansEncPut(&rans, &ptr, k, 1, prob_bits);
k = ESCAPE_LIMIT;
}
- if (k != 0) {
-#if 1
- if (free_sign_bits == 0) {
- --sign_ptr;
- *sign_ptr = 0;
- free_sign_bits = 8;
- }
- *sign_ptr <<= 1;
- *sign_ptr |= (signed_k < 0);
- --free_sign_bits;
-#else
- RansEncPut(&rans, &ptr, (k < 0) ? prob_scale / 2 : 0, prob_scale / 2, prob_bits);
-#endif
+ RansEncPutSymbol(&rans, &ptr, &esyms[(k - 1) & 255]);
+ if (signed_k < 0) {
+ rans += sign_bias;
}
- RansEncPutSymbol(&rans, &ptr, &esyms[k]);
}
private:
static constexpr size_t out_max_size = 32 << 20; // 32 MB.
static constexpr size_t max_num_sign = 1048576; // Way too big. And actually bytes.
- unique_ptr<uint8_t[]> out_buf, sign_buf;
- uint8_t *out_end, *sign_end;
- uint8_t *ptr, *sign_ptr;
+ unique_ptr<uint8_t[]> out_buf;
+ uint8_t *out_end;
+ uint8_t *ptr;
RansState rans;
- size_t free_sign_bits;
RansEncSymbol esyms[NUM_SYMS];
+ uint32_t sign_bias;
+
+ uint32_t last_block = 0; // Not a valid 4-byte rANS block (?)
};
static constexpr int dc_scalefac = 8; // Matches the FDCT's gain.
fclose(fp);
// For each coefficient, make some tables.
- size_t extra_bits = 0, sign_bits = 0;
+ size_t extra_bits = 0;
for (unsigned i = 0; i < 64; ++i) {
stats[i].clear();
}
// Luma
for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
for (unsigned xb = 0; xb < WIDTH; xb += 8) {
- short k = abs(coeff_y[(yb + y) * WIDTH + (xb + x)]);
+ unsigned short k = abs(coeff_y[(yb + y) * WIDTH + (xb + x)]);
if (k >= ESCAPE_LIMIT) {
k = ESCAPE_LIMIT;
extra_bits += 12; // escape this one
}
- if (k != 0) ++sign_bits;
- ++s_luma.freqs[k];
+ ++s_luma.freqs[(k - 1) & 255];
}
}
// Chroma
for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
for (unsigned xb = 0; xb < WIDTH/2; xb += 8) {
- short k_cb = abs(coeff_cb[(yb + y) * WIDTH/2 + (xb + x)]);
- short k_cr = abs(coeff_cr[(yb + y) * WIDTH/2 + (xb + x)]);
+ unsigned short k_cb = abs(coeff_cb[(yb + y) * WIDTH/2 + (xb + x)]);
+ unsigned short k_cr = abs(coeff_cr[(yb + y) * WIDTH/2 + (xb + x)]);
if (k_cb >= ESCAPE_LIMIT) {
k_cb = ESCAPE_LIMIT;
extra_bits += 12; // escape this one
k_cr = ESCAPE_LIMIT;
extra_bits += 12; // escape this one
}
- if (k_cb != 0) ++sign_bits;
- if (k_cr != 0) ++sign_bits;
- ++s_chroma.freqs[k_cb];
- ++s_chroma.freqs[k_cr];
+ ++s_chroma.freqs[(k_cb - 1) & 255];
+ ++s_chroma.freqs[(k_cr - 1) & 255];
}
}
}
}
for (unsigned i = 0; i < 64; ++i) {
+ stats[i].freqs[255] /= 2; // zero, has no sign bits (yes, this is trickery)
stats[i].normalize_freqs(prob_scale);
+ stats[i].cum_freqs[256] += stats[i].freqs[255];
+ stats[i].freqs[255] *= 2;
}
FILE *codedfp = fopen("coded.dat", "wb");
rans_encoder.init_prob(s_luma);
// Luma
+ std::vector<int> lens;
// need to reverse later
rans_encoder.clear();
rans_encoder.encode_coeff(k);
}
if (yb % 16 == 8) {
- num_bytes += rans_encoder.save_block(codedfp);
+ int l = rans_encoder.save_block(codedfp);
+ num_bytes += l;
+ lens.push_back(l);
}
}
if (HEIGHT % 16 != 0) {
}
tot_bytes += num_bytes;
printf("coeff %d Y': %ld bytes\n", y * 8 + x, num_bytes);
+
+ double sum_l = 0.0;
+ for (int l : lens) {
+ sum_l += l;
+ }
+ double avg_l = sum_l / lens.size();
+
+ double sum_sql = 0.0;
+ for (int l : lens) {
+ sum_sql += (l - avg_l) * (l - avg_l);
+ }
+ double stddev_l = sqrt(sum_sql / (lens.size() - 1));
+ printf("coeff %d: avg=%.2f bytes, stddev=%.2f bytes\n", y*8+x, avg_l, stddev_l);
}
}
}
}
- printf("%ld bytes + %ld sign bits (%ld) + %ld escape bits (%ld) = %ld total bytes\n",
- tot_bytes - sign_bits / 8 - extra_bits / 8,
- sign_bits,
- sign_bits / 8,
+ printf("%ld bytes + %ld escape bits (%ld) = %ld total bytes\n",
+ tot_bytes - extra_bits / 8,
extra_bits,
extra_bits / 8,
tot_bytes);