]> git.sesse.net Git - narabu/commitdiff
Encoder with 4x4 blocks (using TF switching).
authorSteinar H. Gunderson <sgunderson@bigfoot.com>
Sat, 16 Sep 2017 13:01:46 +0000 (15:01 +0200)
committerSteinar H. Gunderson <sgunderson@bigfoot.com>
Sat, 16 Sep 2017 13:57:57 +0000 (15:57 +0200)
curve.pl [new file with mode: 0644]
qdc.cpp

diff --git a/curve.pl b/curve.pl
new file mode 100644 (file)
index 0000000..ce322d0
--- /dev/null
+++ b/curve.pl
@@ -0,0 +1,26 @@
+#! /usr/bin/perl
+use strict;
+use warnings;
+
+my @alts = ();
+while (<>) {
+       chomp;
+       my @x = split / /;
+       push @x, 0;
+       push @alts, \@x; 
+}
+
+for my $i (0..$#alts) {
+       next if ($alts[$i]->[4]);
+       for my $j (($i+1)..$#alts) {
+               if ($alts[$i]->[2] <= $alts[$j]->[2] && $alts[$i]->[3] >= $alts[$j]->[3]) {
+                       $alts[$j]->[4] = 1;
+               }
+       }
+}
+
+@alts = sort { $a->[2] <=> $b->[2] } @alts;
+for my $i (0..$#alts) {
+       next if ($alts[$i]->[4]);
+       print join(' ', @{$alts[$i]}), "\n";
+}
diff --git a/qdc.cpp b/qdc.cpp
index 1af42d2d94b79be478b4702871b2a176d7c96500..b97f9a93fb7624523aedf311302144d50504745b 100644 (file)
--- a/qdc.cpp
+++ b/qdc.cpp
@@ -1,6 +1,7 @@
 #include <stdio.h>
 #include <stdint.h>
 #include <stdlib.h>
+#include <string.h>
 #include <assert.h>
 #include <math.h>
 
 
 using namespace std;
 
+static constexpr int dc_scalefac = 8;  // Matches the FDCT's gain.
+static double quant_scalefac = 5.0;  // whatever?
+static double lambda = 0.1;
+
 void fdct_int32(short *const In);
 void idct_int32(short *const In);
 
-unsigned char pix[WIDTH * HEIGHT];
-short coeff[WIDTH * HEIGHT];
+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)];
 
-static const unsigned char std_luminance_quant_tbl[64] = {
+static const unsigned char quant_8x8[64] = {
 #if 0
        16,  11,  10,  16,  24,  40,  51,  61,
        12,  12,  14,  19,  26,  58,  60,  55,
@@ -32,15 +38,31 @@ static const unsigned char std_luminance_quant_tbl[64] = {
        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
+       // ff_mpeg1_default_intra_matrix
+         8, 16, 19, 22, 26, 27, 29, 34,
+        16, 16, 22, 24, 27, 29, 34, 37,                                                 
+        19, 22, 26, 27, 29, 34, 34, 38,                                                 
+        22, 22, 26, 27, 29, 34, 37, 40,
+        22, 26, 27, 29, 32, 35, 40, 48,
+        26, 27, 29, 32, 35, 40, 48, 58,
+        26, 27, 29, 34, 38, 46, 56, 69,
+        27, 29, 35, 38, 46, 56, 69, 83
 #endif
-    16, 16, 19, 22, 26, 27, 29, 34,
-    16, 16, 22, 24, 27, 29, 34, 37,
-    19, 22, 26, 27, 29, 34, 34, 38,
-    22, 22, 26, 27, 29, 34, 37, 40,
-    22, 26, 27, 29, 32, 35, 40, 48,
-    26, 27, 29, 32, 35, 40, 48, 58,
-    26, 27, 29, 34, 38, 46, 56, 69,
-    27, 29, 35, 38, 46, 56, 69, 83
+};
+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
@@ -135,13 +157,16 @@ void SymbolStats::normalize_freqs(uint32_t target_total)
 
 SymbolStats stats[64];
 
-int pick_stats_for(int y, int x)
+int pick_stats_for(int y, int x, bool is_4x4)
 {
+       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;
+       //return y * 8 + x;
 #if 0
        if (y == 0 && x == 0) {
                return 0;
@@ -173,11 +198,11 @@ public:
                clear();
        }
 
-       void init_prob(const SymbolStats &s)
+       void init_prob(const SymbolStats &s1, const SymbolStats &s2)
        {
                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);
+                       RansEncSymbolInit(&esyms[i], s1.cum_freqs[i], s1.freqs[i], prob_bits);
                }
        }
 
@@ -260,8 +285,156 @@ private:
        RansEncSymbol esyms[NUM_SYMS];
 };
 
-int main(void)
+static inline int quantize8x8(int f, int coeff_idx)
+{
+       if (coeff_idx == 0) {
+               return f / dc_scalefac;
+       }
+       if (f == 0) {
+               return 0;
+       }
+
+       const int w = quant_8x8[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)
+{
+       if (coeff_idx == 0) {
+               return qf * dc_scalefac;
+       }
+       if (qf == 0) {
+               return 0;
+       }
+
+       const int w = quant_8x8[coeff_idx];
+       const int s = quant_scalefac;
+       return (2 * qf * w * s) / 32;
+}
+
+static inline int quantize4x4(int f, int coeff_idx)
+{
+       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;
+       }
+
+       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
+
+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;
+}
+
+static inline void tf_switch_second_stage(short *b, short *d, short *f, short *h)
+{
+       *b += *d / 2;
+       *d -= *b / 2;
+       *d += *f / 2;
+       *f -= *d / 2;
+       *f += *h / 2;
+       *h -= *f / 2;
+}
+
+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;
+}
+
+static void convert_8x8to4x4(short *c)
 {
+       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)]);
+               }
+       }
+       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));
+}
+
+static void convert_4x4to8x8(short *c)
+{
+       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)]);
+               }
+       }
+       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]);
+       }
+
+       memcpy(c, d, sizeof(d));
+}
+
+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);
@@ -271,73 +444,108 @@ int main(void)
        for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
                for (unsigned xb = 0; xb < WIDTH; xb += 8) {
                        // Read one block
-                       short in[64];
+                       short in[64], reconstructed8x8[64];
+                       short reconstructed4x4[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;
                                }
                        }
 
                        // FDCT it
                        fdct_int32(in);
 
-                       //constexpr int extra_deadzone = 64;
-                       constexpr int extra_deadzone = 4;
-
+                       // quant 8x8
                        for (unsigned y = 0; y < 8; ++y) {
                                for (unsigned x = 0; x < 8; ++x) {
-                                       short *c = &in[y * 8 + x];
-                                       *c <<= 3;
-                                       *c = copysign(std::max(abs(*c) - extra_deadzone, 0), *c);
-                                       //*c /= std_luminance_quant_tbl[y * 8 + x];
-                                       *c = (int)(double(*c) / std_luminance_quant_tbl[y * 8 + x]);
+                                       int coeff_idx = y * 8 + x;
+                                       int k = quantize8x8(in[coeff_idx], coeff_idx);
+                                       global_coeff8x8[(yb + y) * WIDTH + (xb + x)] = k;
+
+                                       // Store back for reconstruction / PSNR calculation
+                                       reconstructed8x8[coeff_idx] = unquantize8x8(k, coeff_idx);
+                               }
+                       }
+
 #if 0
-                                       if (x != 0 || y != 0) {
-                                               int ss = 1;
-                                               if (::abs(int(*c)) <= ss) {
-                                                       *c = 0; // eeh
-                                               } else if (*c > 0) {
-                                                       *c -= ss;  // eeh
-                                               } else {
-                                                       *c += ss;  // eeh
-                                               }
-                                       }
-#endif
+                       printf("before TF switch:\n");
+                       for (unsigned y = 0; y < 8; ++y) {
+                               for (unsigned x = 0; x < 8; ++x) {
+                                       printf("%4d ", in[y * 8 + x]);
+                               }
+                               printf("\n");
+                       }
+                       convert_8x8to4x4(in);
+                       printf("after TF switch:\n");
+                       for (unsigned y = 0; y < 8; ++y) {
+                               for (unsigned x = 0; x < 8; ++x) {
+                                       printf("%4d ", in[y * 8 + x]);
+                               }
+                               printf("\n");
+                       }
+                       convert_4x4to8x8(in);
+                       printf("after TF switch and back:\n");
+                       for (unsigned y = 0; y < 8; ++y) {
+                               for (unsigned x = 0; x < 8; ++x) {
+                                       printf("%4d ", in[y * 8 + x]);
                                }
+                               printf("\n");
                        }
+#endif
+
+                       // reconstruct 8x8
+                       idct_int32(reconstructed8x8);
 
-                       // Store it
+                       double sum_sq_err8x8 = 0.0;
                        for (unsigned y = 0; y < 8; ++y) {
                                for (unsigned x = 0; x < 8; ++x) {
-                                       coeff[(yb + y) * WIDTH + (xb + x)] = in[y * 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;
                                }
                        }
+                       sum_sq_err += sum_sq_err8x8;
 
-                       // and back
+                       // now let's try 4x4
+                       convert_8x8to4x4(in);
                        for (unsigned y = 0; y < 8; ++y) {
                                for (unsigned x = 0; x < 8; ++x) {
-                                       in[y * 8 + x] *= std_luminance_quant_tbl[y * 8 + x];
-                                       if (in[y * 8 + x] > 0) {
-                                               in[y * 8 + x] += extra_deadzone;
-                                       } else if (in[y * 8 + x] < 0) {
-                                               in[y * 8 + x] -= extra_deadzone;
-                                       }
-                                       in[y * 8 + x] >>= 3;
+                                       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;
+
+                                       // Store back for reconstruction / PSNR calculation
+                                       reconstructed4x4[coeff_idx] = unquantize4x4(k, subcoeff_idx);
                                }
                        }
 
-                       idct_int32(in);
+                       // reconstruct 4x4
+                       convert_4x4to8x8(reconstructed4x4);
+                       idct_int32(reconstructed4x4);
 
+                       double sum_sq_err4x4 = 0.0;
                        for (unsigned y = 0; y < 8; ++y) {
                                for (unsigned x = 0; x < 8; ++x) {
-                                       int k = in[y * 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_err += (*ptr - k) * (*ptr - k);
-                                       *ptr = k;
+                                       sum_sq_err4x4 += (*ptr - k) * (*ptr - k);
+                                       //*ptr = k;
+                                       pix_4x4[(yb + y) * WIDTH + (xb + x)] = k;
                                }
                        }
+
+                       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;
                }
        }
        double mse = sum_sq_err / double(WIDTH * HEIGHT);
@@ -347,14 +555,14 @@ int main(void)
        // 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) {
-                       coeff[yb * WIDTH + xb] -= coeff[yb * WIDTH + (xb + 8)];
+                       global_coeff8x8[yb * WIDTH + xb] -= global_coeff8x8[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)];
                }
        }
-
-       fp = fopen("reconstructed.pgm", "wb");
-       fprintf(fp, "P5\n%d %d\n255\n", WIDTH, HEIGHT);
-       fwrite(pix, 1, WIDTH * HEIGHT, fp);
-       fclose(fp);
 
        // For each coefficient, make some tables.
        size_t extra_bits = 0, sign_bits = 0;
@@ -363,11 +571,11 @@ int main(void)
        }
        for (unsigned y = 0; y < 8; ++y) {
                for (unsigned x = 0; x < 8; ++x) {
-                       SymbolStats &s = stats[pick_stats_for(x, y)];
+                       SymbolStats &s = stats[pick_stats_for(x, y, false)];
 
                        for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
                                for (unsigned xb = 0; xb < WIDTH; xb += 8) {
-                                       short k = abs(coeff[(yb + y) * WIDTH + (xb + x)]);
+                                       short k = abs(global_coeff8x8[(yb + y) * WIDTH + (xb + x)]);
                                        if (k >= ESCAPE_LIMIT) {
                                                //printf("coeff (%d,%d) had value %d\n", y, x, k);
                                                k = ESCAPE_LIMIT;
@@ -380,6 +588,23 @@ int main(void)
                        }
                }
        }
+       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;
+                                               extra_bits += 12;  // escape this one
+                                       }
+                                       if (k != 0) ++sign_bits;
+                                       ++s.freqs[k];
+                               }
+                       }
+               }
+       }
        for (unsigned i = 0; i < 64; ++i) {
 #if 0
                printf("coeff %i:", i);
@@ -412,40 +637,78 @@ int main(void)
 #else
                // TODO: rather gamma-k or something
                for (unsigned j = 0; j < NUM_SYMS; ++j) {
-               //      write_varint(stats[i].freqs[j], codedfp);
+                       write_varint(stats[i].freqs[j], codedfp);
                }
 #endif
        }
 
-       RansEncoder rans_encoder;
+       RansEncoder rans_encoder_8x8, rans_encoder_4x4;
 
-       size_t tot_bytes = 0;
-       for (unsigned y = 0; y < 8; ++y) {
-               for (unsigned x = 0; x < 8; ++x) {
-                       SymbolStats &s = stats[pick_stats_for(x, y)];
+       double total_bits_8x8 = 0.0, total_bits_4x4 = 0.0, total_bits_chosen = 0.0;
+       int num_chosen = 0, tot_chosen = 0;
 
-                       rans_encoder.init_prob(s);
+       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;
 
-                       // need to reverse later
-                       rans_encoder.clear();
+                       //rans_encoder.init_prob(s1, s2);
+                       //rans_encoder.clear();
                        size_t num_bytes = 0;
-                       for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
-                               for (unsigned xb = 0; xb < WIDTH; xb += 8) {
-                                       int k = coeff[(yb + y) * WIDTH + (xb + x)];
-                                       //printf("encoding coeff %d xb,yb=%d,%d: %d\n", y*8+x, xb, yb, k);
-                                       rans_encoder.encode_coeff(k);
-                               }
-                               if (yb % 16 == 8) {
-                                       num_bytes += rans_encoder.save_block(codedfp);
+                       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)];
+
+                                       int k8 = global_coeff8x8[(yb + y) * WIDTH + (xb + x)];
+                                       int k4 = global_coeff4x4[(yb + y) * WIDTH + (xb + x)];
+
+                                       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);
                                }
+//                             if (yb % 16 == 8) {
+//                                     num_bytes += rans_encoder.save_block(codedfp);
+//                             }
                        }
-                       if (HEIGHT % 16 != 0) {
-                               num_bytes += rans_encoder.save_block(codedfp);
-                       }
+//                     if (HEIGHT % 16 != 0) {
+//                             num_bytes += rans_encoder.save_block(codedfp);
+//                     }
                        tot_bytes += num_bytes;
-                       printf("coeff %d: %ld bytes\n", y * 8 + x, 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("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);
                }
        }
+       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,
@@ -453,4 +716,10 @@ int main(void)
                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);
+
 }