git.sesse.net Git - narabu/blob - qdd.cpp

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <stdint.h>
   4 #include <assert.h>
   5 #include <memory>
   6
   7 #define WIDTH 1280
   8 #define HEIGHT 720
   9 #define NUM_SYMS 256
  10 #define ESCAPE_LIMIT (NUM_SYMS - 1)
  11
  12 #include "ryg_rans/rans_byte.h"
  13
  14 using namespace std;
  15
  16 void fdct_int32(short *const In);
  17 void idct_int32(short *const In);
  18
  19 unsigned char pix[WIDTH * HEIGHT];
  20 short coeff[WIDTH * HEIGHT];
  21
  22 struct RansDecoder {
  23 };
  24
  25 static constexpr uint32_t prob_bits = 12;
  26 static constexpr uint32_t prob_scale = 1 << prob_bits;
  27
  28 struct RansDecodeTable {
  29         int cum2sym[prob_scale];
  30         RansDecSymbol dsyms[NUM_SYMS];
  31 };
  32 #define NUM_TABLES 8
  33 RansDecodeTable decode_tables[NUM_TABLES];
  34
  35 static const unsigned char std_luminance_quant_tbl[64] = {
  36 #if 0
  37         16,  11,  10,  16,  24,  40,  51,  61,
  38         12,  12,  14,  19,  26,  58,  60,  55,
  39         14,  13,  16,  24,  40,  57,  69,  56,
  40         14,  17,  22,  29,  51,  87,  80,  62,
  41         18,  22,  37,  56,  68, 109, 103,  77,
  42         24,  35,  55,  64,  81, 104, 113,  92,
  43         49,  64,  78,  87, 103, 121, 120, 101,
  44         72,  92,  95,  98, 112, 100, 103,  99
  45 #else
  46         // ff_mpeg1_default_intra_matrix
  47          8, 16, 19, 22, 26, 27, 29, 34,
  48         16, 16, 22, 24, 27, 29, 34, 37,
  49         19, 22, 26, 27, 29, 34, 34, 38,
  50         22, 22, 26, 27, 29, 34, 37, 40,
  51         22, 26, 27, 29, 32, 35, 40, 48,
  52         26, 27, 29, 32, 35, 40, 48, 58,
  53         26, 27, 29, 34, 38, 46, 56, 69,
  54         27, 29, 35, 38, 46, 56, 69, 83
  55 #endif
  56 };
  57
  58
  59 const int luma_mapping[64] = {
  60         0, 0, 1, 1, 2, 2, 3, 3,
  61         0, 0, 1, 2, 2, 2, 3, 3,
  62         1, 1, 2, 2, 2, 3, 3, 3,
  63         1, 1, 2, 2, 2, 3, 3, 3,
  64         1, 2, 2, 2, 2, 3, 3, 3,
  65         2, 2, 2, 2, 3, 3, 3, 3,
  66         2, 2, 3, 3, 3, 3, 3, 3,
  67         3, 3, 3, 3, 3, 3, 3, 3,
  68 };
  69 const int chroma_mapping[64] = {
  70         0, 1, 1, 2, 2, 2, 3, 3,
  71         1, 1, 2, 2, 2, 3, 3, 3,
  72         2, 2, 2, 2, 3, 3, 3, 3,
  73         2, 2, 2, 3, 3, 3, 3, 3,
  74         2, 3, 3, 3, 3, 3, 3, 3,
  75         3, 3, 3, 3, 3, 3, 3, 3,
  76         3, 3, 3, 3, 3, 3, 3, 3,
  77         3, 3, 3, 3, 3, 3, 3, 3,
  78 };
  79
  80 int pick_stats_for(int x, int y, bool is_chroma)
  81 {
  82         if (is_chroma) {
  83                 return chroma_mapping[y * 8 + x] + 4;
  84         } else {
  85                 return luma_mapping[y * 8 + x];
  86         }
  87 }
  88
  89 uint32_t read_varint(FILE *fp)
  90 {
  91         uint32_t x = 0;
  92         int shift = 0;
  93         for ( ;; ) {
  94                 int ch = getc(fp);
  95                 if (ch == -1) {
  96                         fprintf(stderr, "Premature EOF\n");
  97                         exit(1);
  98                 }
  99
 100                 x |= (ch & 0x7f) << shift;
 101                 if ((ch & 0x80) == 0) return x;
 102                 shift += 7;
 103         }
 104 }
 105
 106 static constexpr int dc_scalefac = 8;  // Matches the FDCT's gain.
 107 static constexpr double quant_scalefac = 4.0;  // whatever?
 108
 109 static inline int unquantize(int qf, int coeff_idx)
 110 {
 111         if (coeff_idx == 0) {
 112                 return qf * dc_scalefac;
 113         }
 114         if (qf == 0) {
 115                 return 0;
 116         }
 117
 118         const int w = std_luminance_quant_tbl[coeff_idx];
 119         const int s = quant_scalefac;
 120         return (2 * qf * w * s) / 32;
 121 }
 122
 123 int main(void)
 124 {
 125         FILE *fp = fopen("coded.dat", "rb");
 126         if (fp == nullptr) {
 127                 perror("coded.dat");
 128                 exit(1);
 129         }
 130
 131         uint32_t sign_bias[NUM_TABLES];
 132         for (unsigned table = 0; table < NUM_TABLES; ++table) {
 133                 uint32_t cum_freq = 0;
 134                 for (unsigned sym = 0; sym < NUM_SYMS; ++sym) {
 135                         uint32_t freq = read_varint(fp);
 136                         fprintf(stderr, "sym=%u/%u: freq=%u\n", sym, NUM_SYMS, freq);
 137                         RansDecSymbolInit(&decode_tables[table].dsyms[(sym + 1) & 255], cum_freq, freq);
 138                         for (uint32_t i = 0; i < freq; ++i) {
 139                                 if (cum_freq < prob_scale)
 140                                         decode_tables[table].cum2sym[cum_freq] = (sym + 1) & 255;
 141                                 ++cum_freq;
 142                         }
 143                 }
 144                 sign_bias[table] = cum_freq;
 145                 printf("sign_bias=%u (of %d)\n", sign_bias[table], prob_scale * 2);
 146         }
 147
 148         // loop over all coefficients
 149         for (unsigned y = 0; y < 8; ++y) {
 150                 for (unsigned x = 0; x < 8; ++x) {
 151                         unsigned tbl = pick_stats_for(x, y, false);
 152
 153                         RansState rans = 0;
 154
 155                         unique_ptr<uint8_t[]> rans_bytes;
 156                         uint8_t *rans_ptr = nullptr;
 157
 158                         // loop over all DCT blocks
 159                         for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
 160                                 if (yb % 16 == 0) {
 161                                         // read a block
 162                                         uint32_t num_rans_bytes = read_varint(fp);
 163                                         rans_bytes.reset(new uint8_t[num_rans_bytes]);
 164                                         fread(rans_bytes.get(), 1, num_rans_bytes, fp);
 165
 166                                         printf("%d,%d: read %d rANS bytes\n", x, y, num_rans_bytes);
 167                                         //printf("first bytes: %02x %02x %02x %02x %02x %02x %02x %02x\n", rans_bytes[0], rans_bytes[1], rans_bytes[2], rans_bytes[3], rans_bytes[4], rans_bytes[5], rans_bytes[6], rans_bytes[7]);
 168
 169
 170                                         // init rANS decoder
 171                                         rans_ptr = rans_bytes.get();
 172                                         RansDecInit(&rans, &rans_ptr);
 173                                 }
 174                                 for (unsigned xb = 0; xb < WIDTH; xb += 8) {
 175                                         uint32_t bottom_bits = RansDecGet(&rans, prob_bits + 1);
 176                                         uint32_t sign = 0;
 177                                         if (bottom_bits >= sign_bias[tbl]) {
 178                                                 bottom_bits -= sign_bias[tbl];
 179                                                 rans -= sign_bias[tbl];
 180                                                 sign = 1;
 181                                         }
 182                                         uint32_t k = decode_tables[tbl].cum2sym[std::min(bottom_bits, prob_scale - 1)];
 183                                         RansDecAdvanceSymbol(&rans, &rans_ptr, &decode_tables[tbl].dsyms[k], prob_bits + 1);
 184                                         assert(k <= ESCAPE_LIMIT);
 185                                         if (k == ESCAPE_LIMIT) {
 186                                                 k = RansDecGet(&rans, prob_bits);
 187                                                 assert(k >= ESCAPE_LIMIT);
 188                                                 RansDecAdvance(&rans, &rans_ptr, k, 1, prob_bits);
 189                                         }
 190                                         if (sign) {
 191                                                 assert(k != 0);
 192                                                 k = -k;
 193                                         }
 194
 195                                         // reverse
 196                                         int reversed_yb = yb ^ 8;
 197                                         int reversed_xb = WIDTH - 8 - xb;
 198                                         coeff[(reversed_yb + y) * WIDTH + (reversed_xb + x)] = k;
 199                 //                      printf("coeff %d xb,yb=%d,%d: decoded %d\n", y * 8 + x, reversed_xb, reversed_yb, k);
 200                                 }
 201                         }
 202
 203 #if 0
 204                         for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
 205                                 for (unsigned xb = 0; xb < WIDTH; xb += 8) {
 206 #endif
 207                 }
 208         }
 209         fclose(fp);
 210
 211         // DC coefficient pred from the right to left
 212         for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
 213                 for (int xb = WIDTH - 16; xb >= 0; xb -= 8) {
 214                         coeff[yb * WIDTH + xb] += coeff[yb * WIDTH + (xb + 8)];
 215                 }
 216         }
 217
 218         // IDCT
 219         for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
 220                 for (unsigned xb = 0; xb < WIDTH; xb += 8) {
 221                         // Read one block
 222                         short in[64];
 223                         for (unsigned y = 0; y < 8; ++y) {
 224                                 for (unsigned x = 0; x < 8; ++x) {
 225                                         int k = coeff[(yb + y) * WIDTH + (xb + x)];
 226                                         in[y * 8 + x] = unquantize(k, y * 8 + x);
 227                                         printf("%3d ", in[y * 8 + x]);
 228                                 }
 229                                 printf("\n");
 230                         }
 231                         printf("\n");
 232
 233                         idct_int32(in);
 234
 235                         // Clamp and move back
 236                         for (unsigned y = 0; y < 8; ++y) {
 237                                 for (unsigned x = 0; x < 8; ++x) {
 238                                         int k = in[y * 8 + x];
 239                                         printf("%3d ", k);
 240                                         if (k < 0) k = 0;
 241                                         if (k > 255) k = 255;
 242                                         pix[(yb + y) * WIDTH + (xb + x)] = k;
 243                                 }
 244                                 printf("\n");
 245                         }
 246                         printf("\n");
 247                 }
 248         }
 249
 250         fp = fopen("output.pgm", "wb");
 251         fprintf(fp, "P5\n%d %d\n255\n", WIDTH, HEIGHT);
 252         fwrite(pix, 1, WIDTH * HEIGHT, fp);
 253         fclose(fp);
 254 }