//#include "ryg_rans/rans64.h"
#include "ryg_rans/rans_byte.h"
+#include "ryg_rans/renormalize.h"
+#include <algorithm>
#include <memory>
+#include <numeric>
+#include <random>
+#include <vector>
+#include <unordered_map>
#define WIDTH 1280
#define HEIGHT 720
+#define WIDTH_BLOCKS (WIDTH/8)
+#define WIDTH_BLOCKS_CHROMA (WIDTH/16)
+#define HEIGHT_BLOCKS (HEIGHT/8)
+#define NUM_BLOCKS (WIDTH_BLOCKS * HEIGHT_BLOCKS)
+#define NUM_BLOCKS_CHROMA (WIDTH_BLOCKS_CHROMA * HEIGHT_BLOCKS)
#define NUM_SYMS 256
#define ESCAPE_LIMIT (NUM_SYMS - 1)
-using namespace std;
+// If you set this to 1, the program will try to optimize the placement
+// of coefficients to rANS probability distributions. This is randomized,
+// so you might want to run it a few times.
+#define FIND_OPTIMAL_STREAM_ASSIGNMENT 0
+#define NUM_CLUSTERS 4
-static constexpr int dc_scalefac = 8; // Matches the FDCT's gain.
-static double quant_scalefac = 5.0; // whatever?
-static double lambda = 0.1;
+static constexpr uint32_t prob_bits = 12;
+static constexpr uint32_t prob_scale = 1 << prob_bits;
+
+using namespace std;
void fdct_int32(short *const In);
void idct_int32(short *const In);
-unsigned char pix_4x4[WIDTH * HEIGHT], pix_8x8[WIDTH * HEIGHT], pix[WIDTH * HEIGHT];
-short global_coeff8x8[WIDTH * HEIGHT], global_coeff4x4[WIDTH * HEIGHT];
-double err_8x8[(WIDTH/8) * (HEIGHT/8)], err_4x4[(WIDTH/8) * (HEIGHT/8)];
+unsigned char rgb[WIDTH * HEIGHT * 3];
+unsigned char pix_y[WIDTH * HEIGHT];
+unsigned char pix_cb[(WIDTH/2) * HEIGHT];
+unsigned char pix_cr[(WIDTH/2) * HEIGHT];
+unsigned char full_pix_cb[WIDTH * HEIGHT];
+unsigned char full_pix_cr[WIDTH * HEIGHT];
+short coeff_y[WIDTH * HEIGHT], coeff_cb[(WIDTH/2) * HEIGHT], coeff_cr[(WIDTH/2) * HEIGHT];
+
+int clamp(int x)
+{
+ if (x < 0) return 0;
+ if (x > 255) return 255;
+ return x;
+}
-static const unsigned char quant_8x8[64] = {
+static const unsigned char std_luminance_quant_tbl[64] = {
#if 0
16, 11, 10, 16, 24, 40, 51, 61,
12, 12, 14, 19, 26, 58, 60, 55,
24, 35, 55, 64, 81, 104, 113, 92,
49, 64, 78, 87, 103, 121, 120, 101,
72, 92, 95, 98, 112, 100, 103, 99
-#elif 1
+#else
// ff_mpeg1_default_intra_matrix
8, 16, 19, 22, 26, 27, 29, 34,
16, 16, 22, 24, 27, 29, 34, 37,
27, 29, 35, 38, 46, 56, 69, 83
#endif
};
-static const unsigned char quant_4x4[16] = {
- 8, 17, 27, 37,
- 17, 27, 37, 43,
- 27, 37, 43, 49,
- 37, 43, 49, 56
- //8, 8, 8, 8,
- //8, 8, 8, 8,
- //8, 8, 8, 8,
- //8, 8, 8, 8
- // 8, 19, 26, 29,
- //19, 26, 29, 34,
- //22, 27, 32, 40,
- //26, 29, 38, 56,
-};
struct SymbolStats
{
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;
void SymbolStats::normalize_freqs(uint32_t target_total)
{
+ uint64_t real_freq[NUM_SYMS + 1]; // hack
+
assert(target_total >= NUM_SYMS);
calc_cum_freqs();
if (cur_total == 0) return;
- // resample distribution based on cumulative freqs
+ double ideal_cost = 0.0;
for (int i = 1; i <= NUM_SYMS; i++)
- cum_freqs[i] = ((uint64_t)target_total * cum_freqs[i])/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]++;
- }
- }
+ {
+ real_freq[i] = cum_freqs[i] - cum_freqs[i - 1];
+ if (real_freq[i] > 0)
+ ideal_cost -= real_freq[i] * log2(real_freq[i] / double(cur_total));
}
+ OptimalRenormalize(cum_freqs, NUM_SYMS, prob_scale);
+
// 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++) {
// calc updated freq
freqs[i] = cum_freqs[i+1] - cum_freqs[i];
}
+
+ double calc_cost = 0.0;
+ for (int i = 1; i <= NUM_SYMS; i++)
+ {
+ uint64_t freq = cum_freqs[i] - cum_freqs[i - 1];
+ if (real_freq[i] > 0)
+ calc_cost -= real_freq[i] * log2(freq / double(target_total));
+ }
+
+ static double total_loss = 0.0;
+ total_loss += calc_cost - ideal_cost;
+ static double total_loss_with_dp = 0.0;
+ double optimal_cost = 0.0;
+ //total_loss_with_dp += optimal_cost - ideal_cost;
+ printf("ideal cost = %.0f bits, DP cost = %.0f bits, calc cost = %.0f bits (loss = %.2f bytes, total loss = %.2f bytes, total loss with DP = %.2f bytes)\n",
+ ideal_cost, optimal_cost,
+ calc_cost, (calc_cost - ideal_cost) / 8.0, total_loss / 8.0, total_loss_with_dp / 8.0);
}
-SymbolStats stats[64];
+SymbolStats stats[128];
-int pick_stats_for(int y, int x, bool is_4x4)
+#if FIND_OPTIMAL_STREAM_ASSIGNMENT
+// Distance from one stream to the other, based on a hacked-up K-L divergence.
+float kl_dist[64][64];
+#endif
+
+const int luma_mapping[64] = {
+ 0, 0, 1, 1, 2, 2, 3, 3,
+ 0, 0, 1, 2, 2, 2, 3, 3,
+ 1, 1, 2, 2, 2, 3, 3, 3,
+ 1, 1, 2, 2, 2, 3, 3, 3,
+ 1, 2, 2, 2, 2, 3, 3, 3,
+ 2, 2, 2, 2, 3, 3, 3, 3,
+ 2, 2, 3, 3, 3, 3, 3, 3,
+ 3, 3, 3, 3, 3, 3, 3, 3,
+};
+const int chroma_mapping[64] = {
+ 0, 1, 1, 2, 2, 2, 3, 3,
+ 1, 1, 2, 2, 2, 3, 3, 3,
+ 2, 2, 2, 2, 3, 3, 3, 3,
+ 2, 2, 2, 3, 3, 3, 3, 3,
+ 2, 3, 3, 3, 3, 3, 3, 3,
+ 3, 3, 3, 3, 3, 3, 3, 3,
+ 3, 3, 3, 3, 3, 3, 3, 3,
+ 3, 3, 3, 3, 3, 3, 3, 3,
+};
+
+int pick_stats_for(int x, int y, bool is_chroma)
{
- if (is_4x4) {
- return 8 + std::min<int>(x + y, 3);
- }
- //return std::min<int>(hypot(x, y), 7);
- return std::min<int>(x + y, 7);
- //if (x + y >= 7) return 7;
- //return x + y;
- //return y * 8 + x;
-#if 0
- if (y == 0 && x == 0) {
- return 0;
+#if FIND_OPTIMAL_STREAM_ASSIGNMENT
+ return y * 8 + x + is_chroma * 64;
+#else
+ if (is_chroma) {
+ return chroma_mapping[y * 8 + x] + 4;
} else {
- return 1;
+ return luma_mapping[y * 8 + x];
}
#endif
}
class RansEncoder {
public:
- static constexpr uint32_t prob_bits = 12;
- static constexpr uint32_t prob_scale = 1 << prob_bits;
-
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 &s1, const SymbolStats &s2)
+ 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], s1.cum_freqs[i], s1.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[NUM_SYMS];
}
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) & (NUM_SYMS - 1)]);
+ 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 inline int quantize8x8(int f, int coeff_idx)
+static constexpr int dc_scalefac = 8; // Matches the FDCT's gain.
+static constexpr double quant_scalefac = 4.0; // whatever?
+
+static inline int quantize(int f, int coeff_idx)
{
if (coeff_idx == 0) {
return f / dc_scalefac;
return 0;
}
- const int w = quant_8x8[coeff_idx];
+ const int w = std_luminance_quant_tbl[coeff_idx];
const int s = quant_scalefac;
int sign_f = (f > 0) ? 1 : -1;
return (32 * f + sign_f * w * s) / (2 * w * s);
}
-static inline int unquantize8x8(int qf, int coeff_idx)
+static inline int unquantize(int qf, int coeff_idx)
{
if (coeff_idx == 0) {
return qf * dc_scalefac;
return 0;
}
- const int w = quant_8x8[coeff_idx];
+ const int w = std_luminance_quant_tbl[coeff_idx];
const int s = quant_scalefac;
return (2 * qf * w * s) / 32;
}
-static inline int quantize4x4(int f, int coeff_idx)
+void readpix(unsigned char *ptr, const char *filename)
{
- if (coeff_idx == 0) {
- return f / (dc_scalefac/2);
- }
- if (f == 0) {
- return 0;
- }
-
- const int w = quant_4x4[coeff_idx];
- const int s = quant_scalefac;
- int sign_f = (f > 0) ? 1 : -1;
- return (64 * f + sign_f * w * s) / (2 * w * s);
-}
-
-static inline int unquantize4x4(int qf, int coeff_idx)
-{
- if (coeff_idx == 0) {
- return qf * (dc_scalefac/2);
- }
- if (qf == 0) {
- return 0;
+ FILE *fp = fopen(filename, "rb");
+ if (fp == nullptr) {
+ perror(filename);
+ exit(1);
}
- const int w = quant_4x4[coeff_idx];
- const int s = quant_scalefac;
- return (2 * qf * w * s) / 64;
-}
-
-// https://people.xiph.org/~xiphmont/demo/daala/demo3.shtml
+ fseek(fp, 0, SEEK_END);
+ long len = ftell(fp);
+ assert(len >= WIDTH * HEIGHT * 3);
+ fseek(fp, len - WIDTH * HEIGHT * 3, SEEK_SET);
-static inline void tf_switch(short *a, short *b, short *c, short *d)
-{
- *b = *a - *b;
- *c = *c + *d;
- short e = (*c - *b)/2;
- *a = *a + e;
- *d = *d - e;
- *c = *a - *c;
- *b = *b - *d;
+ fread(ptr, 1, WIDTH * HEIGHT * 3, fp);
+ fclose(fp);
}
-static inline void tf_switch_second_stage(short *b, short *d, short *f, short *h)
+void convert_ycbcr()
{
- *b += *d / 2;
- *d -= *b / 2;
- *d += *f / 2;
- *f -= *d / 2;
- *f += *h / 2;
- *h -= *f / 2;
-}
+ double coeff[3] = { 0.2126, 0.7152, 0.0722 }; // sum = 1.0
+ double cb_fac = 1.0 / (coeff[0] + coeff[1] + 1.0f - coeff[2]); // 0.539
+ double cr_fac = 1.0 / (1.0f - coeff[0] + coeff[1] + coeff[2]); // 0.635
+
+ unique_ptr<float[]> temp_cb(new float[WIDTH * HEIGHT]);
+ unique_ptr<float[]> temp_cr(new float[WIDTH * HEIGHT]);
+ for (unsigned yb = 0; yb < HEIGHT; ++yb) {
+ for (unsigned xb = 0; xb < WIDTH; ++xb) {
+ int r = rgb[((yb * WIDTH) + xb) * 3 + 0];
+ int g = rgb[((yb * WIDTH) + xb) * 3 + 1];
+ int b = rgb[((yb * WIDTH) + xb) * 3 + 2];
+ double y = std::min(std::max(coeff[0] * r + coeff[1] * g + coeff[2] * b, 0.0), 255.0);
+ double cb = (b - y) * cb_fac + 128.0;
+ double cr = (r - y) * cr_fac + 128.0;
+ pix_y[(yb * WIDTH) + xb] = lrint(y);
+ temp_cb[(yb * WIDTH) + xb] = cb;
+ temp_cr[(yb * WIDTH) + xb] = cr;
+ full_pix_cb[(yb * WIDTH) + xb] = lrint(std::min(std::max(cb, 0.0), 255.0));
+ full_pix_cr[(yb * WIDTH) + xb] = lrint(std::min(std::max(cr, 0.0), 255.0));
+ }
+ }
-static inline void tf_switch_second_stage_inv(short *b, short *d, short *f, short *h)
-{
- *h += *f / 2;
- *f -= *h / 2;
- *f += *d / 2;
- *d -= *f / 2;
- *d += *b / 2;
- *b -= *d / 2;
+ // Simple 4:2:2 subsampling with left convention.
+ for (unsigned yb = 0; yb < HEIGHT; ++yb) {
+ for (unsigned xb = 0; xb < WIDTH / 2; ++xb) {
+ int c0 = yb * WIDTH + std::max(int(xb) * 2 - 1, 0);
+ int c1 = yb * WIDTH + xb * 2;
+ int c2 = yb * WIDTH + xb * 2 + 1;
+
+ double cb = 0.25 * temp_cb[c0] + 0.5 * temp_cb[c1] + 0.25 * temp_cb[c2];
+ double cr = 0.25 * temp_cr[c0] + 0.5 * temp_cr[c1] + 0.25 * temp_cr[c2];
+ cb = std::min(std::max(cb, 0.0), 255.0);
+ cr = std::min(std::max(cr, 0.0), 255.0);
+ pix_cb[(yb * WIDTH/2) + xb] = lrint(cb);
+ pix_cr[(yb * WIDTH/2) + xb] = lrint(cr);
+ }
+ }
}
-static void convert_8x8to4x4(short *c)
+#if FIND_OPTIMAL_STREAM_ASSIGNMENT
+double find_best_assignment(const int *medoids, int *assignment)
{
- for (unsigned x = 0; x < 8; ++x) {
- tf_switch_second_stage_inv(&c[1 * 8 + x], &c[3 * 8 + x], &c[5 * 8 + x], &c[7 * 8 + x]);
- }
- for (unsigned y = 0; y < 8; ++y) {
- tf_switch_second_stage_inv(&c[y * 8 + 1], &c[y * 8 + 3], &c[y * 8 + 5], &c[y * 8 + 7]);
- }
- for (unsigned y = 0; y < 4; ++y) {
- for (unsigned x = 0; x < 4; ++x) {
- tf_switch(&c[(y*2) * 8 + x*2], &c[(y*2) * 8 + (x*2+1)], &c[(y*2+1)*8 + x*2], &c[(y*2+1)*8 + (x*2+1)]);
+ double current_score = 0.0;
+ for (int i = 0; i < 64; ++i) {
+ int best_medoid = medoids[0];
+ float best_medoid_score = kl_dist[i][medoids[0]];
+ for (int j = 1; j < NUM_CLUSTERS; ++j) {
+ if (kl_dist[i][medoids[j]] < best_medoid_score) {
+ best_medoid = medoids[j];
+ best_medoid_score = kl_dist[i][medoids[j]];
+ }
}
+ assignment[i] = best_medoid;
+ current_score += best_medoid_score;
}
- short d[64] = {
- c[0*8 + 0], c[0*8 + 2], c[0*8 + 4], c[0*8 + 6], c[0*8 + 1], c[0*8 + 3], c[0*8 + 5], c[0*8 + 7],
- c[2*8 + 0], c[2*8 + 2], c[2*8 + 4], c[2*8 + 6], c[2*8 + 1], c[2*8 + 3], c[2*8 + 5], c[2*8 + 7],
- c[4*8 + 0], c[4*8 + 2], c[4*8 + 4], c[4*8 + 6], c[4*8 + 1], c[4*8 + 3], c[4*8 + 5], c[4*8 + 7],
- c[6*8 + 0], c[6*8 + 2], c[6*8 + 4], c[6*8 + 6], c[6*8 + 1], c[6*8 + 3], c[6*8 + 5], c[6*8 + 7],
- c[1*8 + 0], c[1*8 + 2], c[1*8 + 4], c[1*8 + 6], c[1*8 + 1], c[1*8 + 3], c[1*8 + 5], c[1*8 + 7],
- c[3*8 + 0], c[3*8 + 2], c[3*8 + 4], c[3*8 + 6], c[3*8 + 1], c[3*8 + 3], c[3*8 + 5], c[3*8 + 7],
- c[5*8 + 0], c[5*8 + 2], c[5*8 + 4], c[5*8 + 6], c[5*8 + 1], c[5*8 + 3], c[5*8 + 5], c[5*8 + 7],
- c[7*8 + 0], c[7*8 + 2], c[7*8 + 4], c[7*8 + 6], c[7*8 + 1], c[7*8 + 3], c[7*8 + 5], c[7*8 + 7]
- };
- memcpy(c, d, sizeof(d));
+ return current_score;
}
-static void convert_4x4to8x8(short *c)
+void find_optimal_stream_assignment(int base)
{
- short d[64] = {
- c[0*8 + 0], c[0*8 + 4], c[0*8 + 1], c[0*8 + 5], c[0*8 + 2], c[0*8 + 6], c[0*8 + 3], c[0*8 + 7],
- c[4*8 + 0], c[4*8 + 4], c[4*8 + 1], c[4*8 + 5], c[4*8 + 2], c[4*8 + 6], c[4*8 + 3], c[4*8 + 7],
- c[1*8 + 0], c[1*8 + 4], c[1*8 + 1], c[1*8 + 5], c[1*8 + 2], c[1*8 + 6], c[1*8 + 3], c[1*8 + 7],
- c[5*8 + 0], c[5*8 + 4], c[5*8 + 1], c[5*8 + 5], c[5*8 + 2], c[5*8 + 6], c[5*8 + 3], c[5*8 + 7],
- c[2*8 + 0], c[2*8 + 4], c[2*8 + 1], c[2*8 + 5], c[2*8 + 2], c[2*8 + 6], c[2*8 + 3], c[2*8 + 7],
- c[6*8 + 0], c[6*8 + 4], c[6*8 + 1], c[6*8 + 5], c[6*8 + 2], c[6*8 + 6], c[6*8 + 3], c[6*8 + 7],
- c[3*8 + 0], c[3*8 + 4], c[3*8 + 1], c[3*8 + 5], c[3*8 + 2], c[3*8 + 6], c[3*8 + 3], c[3*8 + 7],
- c[7*8 + 0], c[7*8 + 4], c[7*8 + 1], c[7*8 + 5], c[7*8 + 2], c[7*8 + 6], c[7*8 + 3], c[7*8 + 7]
- };
-
- for (unsigned y = 0; y < 4; ++y) {
- for (unsigned x = 0; x < 4; ++x) {
- tf_switch(&d[(y*2) * 8 + x*2], &d[(y*2) * 8 + (x*2+1)], &d[(y*2+1)*8 + x*2], &d[(y*2+1)*8 + (x*2+1)]);
+ double inv_sum[64];
+ for (unsigned i = 0; i < 64; ++i) {
+ double s = 0.0;
+ for (unsigned k = 0; k < NUM_SYMS; ++k) {
+ s += stats[i + base].freqs[k] + 0.5;
}
- }
- for (unsigned y = 0; y < 8; ++y) {
- tf_switch_second_stage(&d[y * 8 + 1], &d[y * 8 + 3], &d[y * 8 + 5], &d[y * 8 + 7]);
- }
- for (unsigned x = 0; x < 8; ++x) {
- tf_switch_second_stage(&d[1 * 8 + x], &d[3 * 8 + x], &d[5 * 8 + x], &d[7 * 8 + x]);
+ inv_sum[i] = 1.0 / s;
}
- memcpy(c, d, sizeof(d));
+ for (unsigned i = 0; i < 64; ++i) {
+ for (unsigned j = 0; j < 64; ++j) {
+ double d = 0.0;
+ for (unsigned k = 0; k < NUM_SYMS; ++k) {
+ double p1 = (stats[i + base].freqs[k] + 0.5) * inv_sum[i];
+ double p2 = (stats[j + base].freqs[k] + 0.5) * inv_sum[j];
+
+ // K-L divergence is asymmetric; this is a hack.
+ d += p1 * log(p1 / p2);
+ d += p2 * log(p2 / p1);
+ }
+ kl_dist[i][j] = d;
+ //printf("%.3f ", d);
+ }
+ //printf("\n");
+ }
+
+ // k-medoids init
+ int medoids[64]; // only the first NUM_CLUSTERS matter
+ bool is_medoid[64] = { false };
+ std::iota(medoids, medoids + 64, 0);
+ std::random_device rd;
+ std::mt19937 g(rd());
+ std::shuffle(medoids, medoids + 64, g);
+ for (int i = 0; i < NUM_CLUSTERS; ++i) {
+ printf("%d ", medoids[i]);
+ is_medoid[medoids[i]] = true;
+ }
+ printf("\n");
+
+ // assign each data point to the closest medoid
+ int assignment[64];
+ double current_score = find_best_assignment(medoids, assignment);
+
+ for (int i = 0; i < 1000; ++i) {
+ printf("iter %d\n", i);
+ bool any_changed = false;
+ for (int m = 0; m < NUM_CLUSTERS; ++m) {
+ for (int o = 0; o < 64; ++o) {
+ if (is_medoid[o]) continue;
+ int old_medoid = medoids[m];
+ medoids[m] = o;
+
+ int new_assignment[64];
+ double candidate_score = find_best_assignment(medoids, new_assignment);
+
+ if (candidate_score < current_score) {
+ current_score = candidate_score;
+ memcpy(assignment, new_assignment, sizeof(assignment));
+
+ is_medoid[old_medoid] = false;
+ is_medoid[medoids[m]] = true;
+ printf("%f: ", current_score);
+ for (int i = 0; i < 64; ++i) {
+ printf("%d ", assignment[i]);
+ }
+ printf("\n");
+ any_changed = true;
+ } else {
+ medoids[m] = old_medoid;
+ }
+ }
+ }
+ if (!any_changed) break;
+ }
+ printf("\n");
+ std::unordered_map<int, int> rmap;
+ for (int i = 0; i < 64; ++i) {
+ if (i % 8 == 0) printf("\n");
+ if (!rmap.count(assignment[i])) {
+ rmap.emplace(assignment[i], rmap.size());
+ }
+ printf("%d, ", rmap[assignment[i]]);
+ }
+ printf("\n");
}
+#endif
int main(int argc, char **argv)
{
- if (argc >= 2) quant_scalefac = atof(argv[1]);
- if (argc >= 3) lambda = atof(argv[2]);
-
- FILE *fp = fopen("pic.pgm", "rb");
- fread(pix, 1, WIDTH * HEIGHT, fp);
- fclose(fp);
+ if (argc >= 2)
+ readpix(rgb, argv[1]);
+ else
+ readpix(rgb, "color.pnm");
+ convert_ycbcr();
double sum_sq_err = 0.0;
+ //double last_cb_cfl_fac = 0.0;
+ //double last_cr_cfl_fac = 0.0;
+ // DCT and quantize luma
for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
for (unsigned xb = 0; xb < WIDTH; xb += 8) {
// Read one block
- short in[64], reconstructed8x8[64];
- short reconstructed4x4[64];
+ short in_y[64];
for (unsigned y = 0; y < 8; ++y) {
for (unsigned x = 0; x < 8; ++x) {
- in[y * 8 + x] = pix[(yb + y) * WIDTH + (xb + x)];
- // in[y * 8 + x] = 128;
+ in_y[y * 8 + x] = pix_y[(yb + y) * WIDTH + (xb + x)];
}
}
// FDCT it
- fdct_int32(in);
+ fdct_int32(in_y);
- // quant 8x8
for (unsigned y = 0; y < 8; ++y) {
for (unsigned x = 0; x < 8; ++x) {
int coeff_idx = y * 8 + x;
- int k = quantize8x8(in[coeff_idx], coeff_idx);
- global_coeff8x8[(yb + y) * WIDTH + (xb + x)] = k;
+ int k = quantize(in_y[coeff_idx], coeff_idx);
+ coeff_y[(yb + y) * WIDTH + (xb + x)] = k;
// Store back for reconstruction / PSNR calculation
- reconstructed8x8[coeff_idx] = unquantize8x8(k, coeff_idx);
+ in_y[coeff_idx] = unquantize(k, coeff_idx);
}
}
+ idct_int32(in_y);
+
+ for (unsigned y = 0; y < 8; ++y) {
+ for (unsigned x = 0; x < 8; ++x) {
+ int k = clamp(in_y[y * 8 + x]);
+ uint8_t *ptr = &pix_y[(yb + y) * WIDTH + (xb + x)];
+ sum_sq_err += (*ptr - k) * (*ptr - k);
+ *ptr = k;
+ }
+ }
+ }
+ }
+ double mse = sum_sq_err / double(WIDTH * HEIGHT);
+ double psnr_db = 20 * log10(255.0 / sqrt(mse));
+ printf("psnr = %.2f dB\n", psnr_db);
+
+ //double chroma_energy = 0.0, chroma_energy_pred = 0.0;
+
+ // DCT and quantize chroma
+ //double last_cb_cfl_fac = 0.0, last_cr_cfl_fac = 0.0;
+ for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
+ for (unsigned xb = 0; xb < WIDTH/2; xb += 8) {
#if 0
- printf("before TF switch:\n");
+ // TF switch: Two 8x8 luma blocks -> one 16x8 block, then drop high frequencies
+ printf("in blocks:\n");
for (unsigned y = 0; y < 8; ++y) {
for (unsigned x = 0; x < 8; ++x) {
- printf("%4d ", in[y * 8 + x]);
+ short a = coeff_y[(yb + y) * WIDTH + (xb*2 + x)];
+ printf(" %4d", a);
+ }
+ printf(" | ");
+ for (unsigned x = 0; x < 8; ++x) {
+ short b = coeff_y[(yb + y) * WIDTH + (xb*2 + x + 8)];
+ printf(" %4d", b);
}
printf("\n");
}
- convert_8x8to4x4(in);
- printf("after TF switch:\n");
+
+ short in_y[64];
+ for (unsigned y = 0; y < 8; ++y) {
+ for (unsigned x = 0; x < 4; ++x) {
+ short a = coeff_y[(yb + y) * WIDTH + (xb*2 + x)];
+ short b = coeff_y[(yb + y) * WIDTH + (xb*2 + x + 8)];
+ b = a - b;
+ a = 2 * a - b;
+ in_y[y * 8 + x * 2 + 0] = a;
+ in_y[y * 8 + x * 2 + 1] = b;
+ }
+ }
+
+ printf("tf-ed block:\n");
for (unsigned y = 0; y < 8; ++y) {
for (unsigned x = 0; x < 8; ++x) {
- printf("%4d ", in[y * 8 + x]);
+ short a = in_y[y * 8 + x];
+ printf(" %4d", a);
}
printf("\n");
}
- convert_4x4to8x8(in);
- printf("after TF switch and back:\n");
+#else
+ // Read Y block with no tf switch (from reconstructed luma)
+ short in_y[64];
for (unsigned y = 0; y < 8; ++y) {
for (unsigned x = 0; x < 8; ++x) {
- printf("%4d ", in[y * 8 + x]);
+ in_y[y * 8 + x] = pix_y[(yb + y) * (WIDTH) + (xb + x) * 2];
}
- printf("\n");
}
+ fdct_int32(in_y);
#endif
- // reconstruct 8x8
- idct_int32(reconstructed8x8);
+ // Read one block
+ short in_cb[64], in_cr[64];
+ for (unsigned y = 0; y < 8; ++y) {
+ for (unsigned x = 0; x < 8; ++x) {
+ in_cb[y * 8 + x] = pix_cb[(yb + y) * (WIDTH/2) + (xb + x)];
+ in_cr[y * 8 + x] = pix_cr[(yb + y) * (WIDTH/2) + (xb + x)];
+ }
+ }
+
+ // FDCT it
+ fdct_int32(in_cb);
+ fdct_int32(in_cr);
+
+#if 0
+ // Chroma from luma
+ double x0 = in_y[1];
+ double x1 = in_y[8];
+ double x2 = in_y[9];
+ double denom = (x0 * x0 + x1 * x1 + x2 * x2);
+ //double denom = (x1 * x1);
+
+ double y0 = in_cb[1];
+ double y1 = in_cb[8];
+ double y2 = in_cb[9];
+ double cb_cfl_fac = (x0 * y0 + x1 * y1 + x2 * y2) / denom;
+ //double cb_cfl_fac = (x1 * y1) / denom;
- double sum_sq_err8x8 = 0.0;
for (unsigned y = 0; y < 8; ++y) {
for (unsigned x = 0; x < 8; ++x) {
- int k = reconstructed8x8[y * 8 + x];
- if (k < 0) k = 0;
- if (k > 255) k = 255;
- uint8_t *ptr = &pix[(yb + y) * WIDTH + (xb + x)];
- sum_sq_err8x8 += (*ptr - k) * (*ptr - k);
- pix_8x8[(yb + y) * WIDTH + (xb + x)] = k;
-// *ptr = k;
+ short a = in_y[y * 8 + x];
+ printf(" %4d", a);
}
+ printf(" | ");
+ for (unsigned x = 0; x < 8; ++x) {
+ short a = in_cb[y * 8 + x];
+ printf(" %4d", a);
+ }
+ printf("\n");
+ }
+ printf("(%d,%d,%d) -> (%d,%d,%d) gives %f\n",
+ in_y[1], in_y[8], in_y[9],
+ in_cb[1], in_cb[8], in_cb[9],
+ cb_cfl_fac);
+
+ y0 = in_cr[1];
+ y1 = in_cr[8];
+ y2 = in_cr[9];
+ double cr_cfl_fac = (x0 * y0 + x1 * y1 + x2 * y2) / denom;
+ //double cr_cfl_fac = (x1 * y1) / denom;
+ printf("cb CfL = %7.3f dc = %5d cr CfL = %7.3f dc = %d\n",
+ cb_cfl_fac, in_cb[0] - in_y[0],
+ cr_cfl_fac, in_cr[0] - in_y[0]);
+
+ if (denom == 0.0) { cb_cfl_fac = cr_cfl_fac = 0.0; }
+
+ // CHEAT
+ //last_cb_cfl_fac = cb_cfl_fac;
+ //last_cr_cfl_fac = cr_cfl_fac;
+
+ for (unsigned coeff_idx = 1; coeff_idx < 64; ++coeff_idx) {
+ //printf("%2d: cb = %3d prediction = %f * %3d = %7.3f\n", coeff_idx, in_cb[coeff_idx], last_cb_cfl_fac, in_y[coeff_idx], last_cb_cfl_fac * in_y[coeff_idx]);
+ //printf("%2d: cr = %3d prediction = %f * %3d = %7.3f\n", coeff_idx, in_cr[coeff_idx], last_cr_cfl_fac, in_y[coeff_idx], last_cr_cfl_fac * in_y[coeff_idx]);
+ double cb_pred = last_cb_cfl_fac * in_y[coeff_idx];
+ chroma_energy += in_cb[coeff_idx] * in_cb[coeff_idx];
+ chroma_energy_pred += (in_cb[coeff_idx] - cb_pred) * (in_cb[coeff_idx] - cb_pred);
+
+ //in_cb[coeff_idx] -= lrint(last_cb_cfl_fac * in_y[coeff_idx]);
+ //in_cr[coeff_idx] -= lrint(last_cr_cfl_fac * in_y[coeff_idx]);
+ //in_cr[coeff_idx] -= lrint(last_cr_cfl_fac * in_y[coeff_idx]);
+ //in_cb[coeff_idx] = lrint(in_y[coeff_idx] * (1.0 / sqrt(2)));
+ //in_cr[coeff_idx] = lrint(in_y[coeff_idx] * (1.0 / sqrt(2)));
+ //in_cb[coeff_idx] = lrint(in_y[coeff_idx]);
+ //in_cr[coeff_idx] = lrint(in_y[coeff_idx]);
}
- sum_sq_err += sum_sq_err8x8;
+ //in_cb[0] += 1024;
+ //in_cr[0] += 1024;
+ //in_cb[0] -= in_y[0];
+ //in_cr[0] -= in_y[0];
+#endif
- // now let's try 4x4
- convert_8x8to4x4(in);
for (unsigned y = 0; y < 8; ++y) {
for (unsigned x = 0; x < 8; ++x) {
int coeff_idx = y * 8 + x;
- int subcoeff_idx = (y%4) * 4 + (x%4);
- int k = quantize4x4(in[coeff_idx], subcoeff_idx);
- global_coeff4x4[(yb + y) * WIDTH + (xb + x)] = k;
+ int k_cb = quantize(in_cb[coeff_idx], coeff_idx);
+ coeff_cb[(yb + y) * (WIDTH/2) + (xb + x)] = k_cb;
+ int k_cr = quantize(in_cr[coeff_idx], coeff_idx);
+ coeff_cr[(yb + y) * (WIDTH/2) + (xb + x)] = k_cr;
// Store back for reconstruction / PSNR calculation
- reconstructed4x4[coeff_idx] = unquantize4x4(k, subcoeff_idx);
+ in_cb[coeff_idx] = unquantize(k_cb, coeff_idx);
+ in_cr[coeff_idx] = unquantize(k_cr, coeff_idx);
}
}
- // reconstruct 4x4
- convert_4x4to8x8(reconstructed4x4);
- idct_int32(reconstructed4x4);
+ idct_int32(in_y); // DEBUG
+ idct_int32(in_cb);
+ idct_int32(in_cr);
- double sum_sq_err4x4 = 0.0;
for (unsigned y = 0; y < 8; ++y) {
for (unsigned x = 0; x < 8; ++x) {
- int k = reconstructed4x4[y * 8 + x];
- if (k < 0) k = 0;
- if (k > 255) k = 255;
- uint8_t *ptr = &pix[(yb + y) * WIDTH + (xb + x)];
- sum_sq_err4x4 += (*ptr - k) * (*ptr - k);
- //*ptr = k;
- pix_4x4[(yb + y) * WIDTH + (xb + x)] = k;
+ pix_cb[(yb + y) * (WIDTH/2) + (xb + x)] = clamp(in_cb[y * 8 + x]);
+ pix_cr[(yb + y) * (WIDTH/2) + (xb + x)] = clamp(in_cr[y * 8 + x]);
+
+ // pix_cb[(yb + y) * (WIDTH/2) + (xb + x)] = in_y[y * 8 + x];
+ // pix_cr[(yb + y) * (WIDTH/2) + (xb + x)] = in_y[y * 8 + x];
}
}
- err_8x8[(yb/8) * (WIDTH/8) + (xb/8)] = sum_sq_err8x8;
- err_4x4[(yb/8) * (WIDTH/8) + (xb/8)] = sum_sq_err4x4;
- //printf("err 8x8 = %6.2f err 4x4 = %6.2f win = %d\n", sum_sq_err8x8, sum_sq_err4x4, sum_sq_err4x4 < sum_sq_err8x8);
- //sum_sq_err += sum_sq_err4x4;
+#if 0
+ last_cb_cfl_fac = cb_cfl_fac;
+ last_cr_cfl_fac = cr_cfl_fac;
+#endif
}
}
- double mse = sum_sq_err / double(WIDTH * HEIGHT);
- double psnr_db = 20 * log10(255.0 / sqrt(mse));
- printf("psnr = %.2f dB\n", psnr_db);
+
+#if 0
+ printf("chroma_energy = %f, with_pred = %f\n",
+ chroma_energy / (WIDTH * HEIGHT), chroma_energy_pred / (WIDTH * HEIGHT));
+#endif
// DC coefficient pred from the right to left
for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
for (unsigned xb = 0; xb < WIDTH - 8; xb += 8) {
- global_coeff8x8[yb * WIDTH + xb] -= global_coeff8x8[yb * WIDTH + (xb + 8)];
+ coeff_y[yb * WIDTH + xb] -= coeff_y[yb * WIDTH + (xb + 8)];
}
}
- for (unsigned yb = 0; yb < HEIGHT; yb += 4) {
- for (unsigned xb = 0; xb < WIDTH - 4; xb += 4) {
- global_coeff4x4[yb * WIDTH + xb] -= global_coeff4x4[yb * WIDTH + (xb + 4)];
+ for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
+ for (unsigned xb = 0; xb < WIDTH/2 - 8; xb += 8) {
+ coeff_cb[yb * WIDTH/2 + xb] -= coeff_cb[yb * WIDTH/2 + (xb + 8)];
+ coeff_cr[yb * WIDTH/2 + xb] -= coeff_cr[yb * WIDTH/2 + (xb + 8)];
+ }
+ }
+
+ FILE *fp = fopen("reconstructed.pgm", "wb");
+ fprintf(fp, "P5\n%d %d\n255\n", WIDTH, HEIGHT);
+ fwrite(pix_y, 1, WIDTH * HEIGHT, fp);
+ fclose(fp);
+
+ fp = fopen("reconstructed.pnm", "wb");
+ fprintf(fp, "P6\n%d %d\n255\n", WIDTH, HEIGHT);
+ for (unsigned yb = 0; yb < HEIGHT; ++yb) {
+ for (unsigned xb = 0; xb < WIDTH; ++xb) {
+ int y = pix_y[(yb * WIDTH) + xb];
+ int cb, cr;
+ int c0 = yb * (WIDTH/2) + xb/2;
+ if (xb % 2 == 0) {
+ cb = pix_cb[c0] - 128.0;
+ cr = pix_cr[c0] - 128.0;
+ } else {
+ int c1 = yb * (WIDTH/2) + std::min<int>(xb/2 + 1, WIDTH/2 - 1);
+ cb = 0.5 * (pix_cb[c0] + pix_cb[c1]) - 128.0;
+ cr = 0.5 * (pix_cr[c0] + pix_cr[c1]) - 128.0;
+ }
+
+ double r = y + 1.5748 * cr;
+ double g = y - 0.1873 * cb - 0.4681 * cr;
+ double b = y + 1.8556 * cb;
+
+ putc(clamp(lrint(r)), fp);
+ putc(clamp(lrint(g)), fp);
+ putc(clamp(lrint(b)), fp);
}
}
+ 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();
}
for (unsigned y = 0; y < 8; ++y) {
for (unsigned x = 0; x < 8; ++x) {
- SymbolStats &s = stats[pick_stats_for(x, y, false)];
+ SymbolStats &s_luma = stats[pick_stats_for(x, y, false)];
+ SymbolStats &s_chroma = stats[pick_stats_for(x, y, true)];
+ // Luma
for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
for (unsigned xb = 0; xb < WIDTH; xb += 8) {
- short k = abs(global_coeff8x8[(yb + y) * WIDTH + (xb + x)]);
+ unsigned short k = abs(coeff_y[(yb + y) * WIDTH + (xb + x)]);
if (k >= ESCAPE_LIMIT) {
- //printf("coeff (%d,%d) had value %d\n", y, x, k);
k = ESCAPE_LIMIT;
extra_bits += 12; // escape this one
}
- //if (y != 0 || x != 0) ++sign_bits;
- if (k != 0) ++sign_bits;
- ++s.freqs[k];
+ ++s_luma.freqs[(k - 1) & (NUM_SYMS - 1)];
}
}
- }
- }
- for (unsigned y = 0; y < 4; ++y) {
- for (unsigned x = 0; x < 4; ++x) {
- SymbolStats &s = stats[pick_stats_for(x, y, true)];
-
- for (unsigned yb = 0; yb < HEIGHT; yb += 4) {
- for (unsigned xb = 0; xb < WIDTH; xb += 4) {
- short k = abs(global_coeff4x4[(yb + y) * WIDTH + (xb + x)]);
- if (k >= ESCAPE_LIMIT) {
- k = ESCAPE_LIMIT;
+ // Chroma
+ for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
+ for (unsigned xb = 0; xb < WIDTH/2; xb += 8) {
+ 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
+ }
+ if (k_cr >= ESCAPE_LIMIT) {
+ k_cr = ESCAPE_LIMIT;
extra_bits += 12; // escape this one
}
- if (k != 0) ++sign_bits;
- ++s.freqs[k];
+ ++s_chroma.freqs[(k_cb - 1) & (NUM_SYMS - 1)];
+ ++s_chroma.freqs[(k_cr - 1) & (NUM_SYMS - 1)];
}
}
}
}
- for (unsigned i = 0; i < 64; ++i) {
-#if 0
- printf("coeff %i:", i);
- for (unsigned j = 0; j <= ESCAPE_LIMIT; ++j) {
- printf(" %d", stats[i].freqs[j]);
- }
- printf("\n");
+
+#if FIND_OPTIMAL_STREAM_ASSIGNMENT
+ printf("Luma:\n");
+ find_optimal_stream_assignment(0);
+ printf("Chroma:\n");
+ find_optimal_stream_assignment(64);
+ exit(0);
#endif
- stats[i].normalize_freqs(RansEncoder::prob_scale);
+
+ for (unsigned i = 0; i < 64; ++i) {
+ stats[i].freqs[NUM_SYMS - 1] /= 2; // zero, has no sign bits (yes, this is trickery)
+ stats[i].normalize_freqs(prob_scale);
+ stats[i].cum_freqs[NUM_SYMS] += stats[i].freqs[NUM_SYMS - 1];
+ stats[i].freqs[NUM_SYMS - 1] *= 2;
}
FILE *codedfp = fopen("coded.dat", "wb");
exit(1);
}
- // TODO: varint or something on the freqs
+ // TODO: rather gamma-k or something
for (unsigned i = 0; i < 64; ++i) {
if (stats[i].cum_freqs[NUM_SYMS] == 0) {
continue;
}
printf("writing table %d\n", i);
-#if 0
- for (unsigned j = 0; j <= NUM_SYMS; ++j) {
- uint16_t freq = stats[i].cum_freqs[j];
- fwrite(&freq, 1, sizeof(freq), codedfp);
- printf("%d: %d\n", j, stats[i].freqs[j]);
- }
-#else
- // TODO: rather gamma-k or something
for (unsigned j = 0; j < NUM_SYMS; ++j) {
write_varint(stats[i].freqs[j], codedfp);
}
-#endif
}
- RansEncoder rans_encoder_8x8, rans_encoder_4x4;
-
- double total_bits_8x8 = 0.0, total_bits_4x4 = 0.0, total_bits_chosen = 0.0;
- int num_chosen = 0, tot_chosen = 0;
+ RansEncoder rans_encoder;
size_t tot_bytes = 0;
- for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
- for (unsigned xb = 0; xb < WIDTH; xb += 8) {
- double bits_8x8 = 0.0, bits_4x4 = 0.0;
- //rans_encoder.init_prob(s1, s2);
- //rans_encoder.clear();
+ // Luma
+ for (unsigned y = 0; y < 8; ++y) {
+ for (unsigned x = 0; x < 8; ++x) {
+ SymbolStats &s_luma = stats[pick_stats_for(x, y, false)];
+ rans_encoder.init_prob(s_luma);
+
+ // Luma
+ std::vector<int> lens;
+
+ // need to reverse later
+ rans_encoder.clear();
size_t num_bytes = 0;
- for (unsigned y = 0; y < 8; ++y) {
- for (unsigned x = 0; x < 8; ++x) {
- SymbolStats &s1 = stats[pick_stats_for(x, y, false)];
- SymbolStats &s2 = stats[pick_stats_for(x%4, y%4, true)];
+ for (unsigned block_idx = 0; block_idx < NUM_BLOCKS; ++block_idx) {
+ unsigned yb = block_idx / WIDTH_BLOCKS;
+ unsigned xb = block_idx % WIDTH_BLOCKS;
+
+ int k = coeff_y[(yb * 8 + y) * WIDTH + (xb * 8 + x)];
+ //printf("encoding coeff %d xb,yb=%d,%d: %d\n", y*8+x, xb, yb, k);
+ rans_encoder.encode_coeff(k);
+
+ if (block_idx % 320 == 319 || block_idx == NUM_BLOCKS - 1) {
+ int l = rans_encoder.save_block(codedfp);
+ num_bytes += l;
+ lens.push_back(l);
+ }
+ }
+ tot_bytes += num_bytes;
+ printf("coeff %d Y': %ld bytes\n", y * 8 + x, num_bytes);
- int k8 = global_coeff8x8[(yb + y) * WIDTH + (xb + x)];
- int k4 = global_coeff4x4[(yb + y) * WIDTH + (xb + x)];
+ double sum_l = 0.0;
+ for (int l : lens) {
+ sum_l += l;
+ }
+ double avg_l = sum_l / lens.size();
- if (k8 != 0) ++bits_8x8; // sign bits
- if (k4 != 0) ++bits_4x4;
- k8 = abs(k8); k4 = abs(k4);
-
- if (k8 >= ESCAPE_LIMIT) { k8 = ESCAPE_LIMIT; bits_8x8 += 12.0; }
- if (k4 >= ESCAPE_LIMIT) { k4 = ESCAPE_LIMIT; bits_4x4 += 12.0; }
-
- bits_8x8 -= log2(s1.freqs[k8] / 4096.0);
- bits_4x4 -= log2(s2.freqs[k4] / 4096.0);
+ 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);
+ }
+ }
+
+ // Cb
+ for (unsigned y = 0; y < 8; ++y) {
+ for (unsigned x = 0; x < 8; ++x) {
+ SymbolStats &s_chroma = stats[pick_stats_for(x, y, true)];
+ rans_encoder.init_prob(s_chroma);
+
+ rans_encoder.clear();
+ size_t num_bytes = 0;
+ for (unsigned block_idx = 0; block_idx < NUM_BLOCKS_CHROMA; ++block_idx) {
+ unsigned yb = block_idx / WIDTH_BLOCKS_CHROMA;
+ unsigned xb = block_idx % WIDTH_BLOCKS_CHROMA;
+
+ int k = coeff_cb[(yb * 8 + y) * WIDTH/2 + (xb * 8 + x)];
+ //printf("encoding coeff %d xb,yb=%d,%d: %d\n", y*8+x, xb, yb, k);
+ rans_encoder.encode_coeff(k);
+
+ if (block_idx % 320 == 319 || block_idx == NUM_BLOCKS - 1) {
+ num_bytes += rans_encoder.save_block(codedfp);
}
-// if (yb % 16 == 8) {
-// num_bytes += rans_encoder.save_block(codedfp);
-// }
}
-// if (HEIGHT % 16 != 0) {
-// num_bytes += rans_encoder.save_block(codedfp);
-// }
tot_bytes += num_bytes;
- total_bits_8x8 += bits_8x8;
- total_bits_4x4 += bits_4x4;
- auto e8 = err_8x8[(yb/8)*(WIDTH/8) + (xb/8)];
- auto e4 = err_4x4[(yb/8)*(WIDTH/8) + (xb/8)];
- double rd8 = sqrt(e8) + lambda * bits_8x8;
- double rd4 = sqrt(e4) + lambda * bits_4x4;
- const unsigned char *spix = (rd4 < rd8) ? pix_4x4 : pix_8x8;
- unsigned char col = (rd4 < rd8) ? 255 : 0;
- total_bits_chosen += (rd4 < rd8) ? bits_4x4 : bits_8x8;
- num_chosen += (rd4 < rd8);
- ++tot_chosen;
- for (unsigned y = 0; y < 8; ++y) {
- for (unsigned x = 0; x < 8; ++x) {
- pix[(yb + y) * WIDTH + (xb + x)] = spix[(yb + y) * WIDTH + (xb + x)];
- //pix[(yb + y) * WIDTH + (xb + x)] = col;
+ printf("coeff %d Cb: %ld bytes\n", y * 8 + x, num_bytes);
+ }
+ }
+
+ // Cr
+ for (unsigned y = 0; y < 8; ++y) {
+ for (unsigned x = 0; x < 8; ++x) {
+ SymbolStats &s_chroma = stats[pick_stats_for(x, y, true)];
+ rans_encoder.init_prob(s_chroma);
+
+ rans_encoder.clear();
+ size_t num_bytes = 0;
+ for (unsigned block_idx = 0; block_idx < NUM_BLOCKS_CHROMA; ++block_idx) {
+ unsigned yb = block_idx / WIDTH_BLOCKS_CHROMA;
+ unsigned xb = block_idx % WIDTH_BLOCKS_CHROMA;
+
+ int k = coeff_cr[(yb * 8 + y) * WIDTH/2 + (xb * 8 + x)];
+ //printf("encoding coeff %d xb,yb=%d,%d: %d\n", y*8+x, xb, yb, k);
+ rans_encoder.encode_coeff(k);
+
+ if (block_idx % 320 == 319 || block_idx == NUM_BLOCKS - 1) {
+ num_bytes += rans_encoder.save_block(codedfp);
}
- }
- printf("block (%d,%d): 8x8 %.2f bits [err=%.2f], 4x4 %.2f bits [err=%.2f], win_bits = %d, win_err = %d, win = %d\n",
- yb, xb,
- bits_8x8, sqrt(e8),
- bits_4x4, sqrt(e4),
- bits_4x4 < bits_8x8,
- e4 < e8,
- rd4 < rd8);
+ }
+ tot_bytes += num_bytes;
+ printf("coeff %d Cr: %ld bytes\n", y * 8 + x, num_bytes);
}
}
- printf("4x4: %.2f bits, 8x8: %.2f bits, chosen: %d/%d times, %.2f bits (%.0f bytes)\n", total_bits_4x4, total_bits_8x8,
- num_chosen, tot_chosen, total_bits_chosen, total_bits_chosen / 8.0);
- 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);
-
- fp = fopen("reconstructed.pgm", "wb");
- fprintf(fp, "P5\n%d %d\n255\n", WIDTH, HEIGHT);
- fwrite(pix, 1, WIDTH * HEIGHT, fp);
- fclose(fp);
-
}