#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
+#include <string.h>
#include <assert.h>
#include <math.h>
//#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)
+#define BLOCKS_PER_STREAM 320
+
+// 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 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[WIDTH * HEIGHT];
-short coeff[WIDTH * HEIGHT];
+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 std_luminance_quant_tbl[64] = {
#if 0
24, 35, 55, 64, 81, 104, 113, 92,
49, 64, 78, 87, 103, 121, 120, 101,
72, 92, 95, 98, 112, 100, 103, 99
+#else
+ // 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
};
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)
+#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)
{
- //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 &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[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 (num_rans_bytes == last_block.size() &&
+ memcmp(last_block.data(), ptr, last_block.size()) == 0) {
+ write_varint(0, codedfp);
+ clear();
+ return 1;
+ } else {
+ last_block = string((const char *)ptr, num_rans_bytes);
+ }
+
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;
+
+ std::string last_block;
};
-int main(void)
+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;
+ }
+ if (f == 0) {
+ return 0;
+ }
+
+ 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 unquantize(int qf, int coeff_idx)
{
- FILE *fp = fopen("pic.pgm", "rb");
- fread(pix, 1, WIDTH * HEIGHT, fp);
+ if (coeff_idx == 0) {
+ return qf * dc_scalefac;
+ }
+ if (qf == 0) {
+ return 0;
+ }
+
+ const int w = std_luminance_quant_tbl[coeff_idx];
+ const int s = quant_scalefac;
+ return (2 * qf * w * s) / 32;
+}
+
+void readpix(unsigned char *ptr, const char *filename)
+{
+ FILE *fp = fopen(filename, "rb");
+ if (fp == nullptr) {
+ perror(filename);
+ exit(1);
+ }
+
+ fseek(fp, 0, SEEK_END);
+ long len = ftell(fp);
+ assert(len >= WIDTH * HEIGHT * 3);
+ fseek(fp, len - WIDTH * HEIGHT * 3, SEEK_SET);
+
+ fread(ptr, 1, WIDTH * HEIGHT * 3, fp);
fclose(fp);
+}
+
+void convert_ycbcr()
+{
+ 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));
+ }
+ }
+
+ // 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);
+ }
+ }
+}
+
+#if FIND_OPTIMAL_STREAM_ASSIGNMENT
+double find_best_assignment(const int *medoids, int *assignment)
+{
+ 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;
+ }
+ return current_score;
+}
+
+void find_optimal_stream_assignment(int base)
+{
+ 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;
+ }
+ inv_sum[i] = 1.0 / s;
+ }
+
+ 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)
+ 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;
+ int max_val_x[8] = {0}, min_val_x[8] = {0};
+ int max_val_y[8] = {0}, min_val_y[8] = {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];
+ 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[y * 8 + x] = pix_y[(yb + y) * WIDTH + (xb + x)];
}
}
// FDCT it
- fdct_int32(in);
+ fdct_int32(in_y);
- //constexpr int extra_deadzone = 64;
- constexpr int extra_deadzone = 4;
+ for (unsigned y = 0; y < 8; ++y) {
+ for (unsigned x = 0; x < 8; ++x) {
+ int coeff_idx = y * 8 + x;
+ int k = quantize(in_y[coeff_idx], coeff_idx);
+ coeff_y[(yb + y) * WIDTH + (xb + x)] = k;
+
+ max_val_x[x] = std::max(max_val_x[x], k);
+ min_val_x[x] = std::min(min_val_x[x], k);
+ max_val_y[y] = std::max(max_val_y[y], k);
+ min_val_y[y] = std::min(min_val_y[y], k);
+
+ // Store back for reconstruction / PSNR calculation
+ 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) {
- 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 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
- 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
- }
- }
+ // 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) {
+ 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");
+ }
+
+ 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) {
+ short a = in_y[y * 8 + x];
+ printf(" %4d", a);
+ }
+ printf("\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) {
+ in_y[y * 8 + x] = pix_y[(yb + y) * (WIDTH) + (xb + x) * 2];
+ }
+ }
+ fdct_int32(in_y);
#endif
+
+ // 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)];
}
}
- // Store it
+ // 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;
+
for (unsigned y = 0; y < 8; ++y) {
for (unsigned x = 0; x < 8; ++x) {
- coeff[(yb + y) * WIDTH + (xb + x)] = in[y * 8 + x];
+ 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]);
}
+ //in_cb[0] += 1024;
+ //in_cr[0] += 1024;
+ //in_cb[0] -= in_y[0];
+ //in_cr[0] -= in_y[0];
+#endif
- // and back
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 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
+ in_cb[coeff_idx] = unquantize(k_cb, coeff_idx);
+ in_cr[coeff_idx] = unquantize(k_cr, coeff_idx);
}
}
- idct_int32(in);
+ idct_int32(in_y); // DEBUG
+ idct_int32(in_cb);
+ idct_int32(in_cr);
for (unsigned y = 0; y < 8; ++y) {
for (unsigned x = 0; x < 8; ++x) {
- int k = in[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;
+ 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];
}
}
+
+#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);
- // 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)];
+#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 (within each slice)
+ for (unsigned block_idx = 0; block_idx < NUM_BLOCKS; block_idx += BLOCKS_PER_STREAM) {
+ int prev_k = 128;
+
+ for (unsigned subblock_idx = BLOCKS_PER_STREAM; subblock_idx --> 0; ) {
+ unsigned yb = (block_idx + subblock_idx) / WIDTH_BLOCKS;
+ unsigned xb = (block_idx + subblock_idx) % WIDTH_BLOCKS;
+ int k = coeff_y[(yb * 8) * WIDTH + (xb * 8)];
+
+ coeff_y[(yb * 8) * WIDTH + (xb * 8)] = k - prev_k;
+
+ prev_k = k;
+ }
+ }
+ for (unsigned block_idx = 0; block_idx < NUM_BLOCKS_CHROMA; block_idx += BLOCKS_PER_STREAM) {
+ int prev_k_cb = 0;
+ int prev_k_cr = 0;
+
+ for (unsigned subblock_idx = BLOCKS_PER_STREAM; subblock_idx --> 0; ) {
+ unsigned yb = (block_idx + subblock_idx) / WIDTH_BLOCKS_CHROMA;
+ unsigned xb = (block_idx + subblock_idx) % WIDTH_BLOCKS_CHROMA;
+ int k_cb = coeff_cb[(yb * 8) * WIDTH/2 + (xb * 8)];
+ int k_cr = coeff_cr[(yb * 8) * WIDTH/2 + (xb * 8)];
+
+ coeff_cb[(yb * 8) * WIDTH/2 + (xb * 8)] = k_cb - prev_k_cb;
+ coeff_cr[(yb * 8) * WIDTH/2 + (xb * 8)] = k_cr - prev_k_cr;
+
+ prev_k_cb = k_cb;
+ prev_k_cr = k_cr;
}
}
- fp = fopen("reconstructed.pgm", "wb");
+ FILE *fp = fopen("reconstructed.pgm", "wb");
fprintf(fp, "P5\n%d %d\n255\n", WIDTH, HEIGHT);
- fwrite(pix, 1, WIDTH * HEIGHT, fp);
+ 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)];
+ 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(coeff[(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)];
+ }
+ }
+ // 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
+ }
+ ++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);
+ write_varint(stats[i].freqs[j], codedfp);
}
-#endif
}
RansEncoder rans_encoder;
size_t tot_bytes = 0;
+
+ // Luma
for (unsigned y = 0; y < 8; ++y) {
for (unsigned x = 0; x < 8; ++x) {
- SymbolStats &s = stats[pick_stats_for(x, y)];
+ SymbolStats &s_luma = stats[pick_stats_for(x, y, false)];
+ rans_encoder.init_prob(s_luma);
- rans_encoder.init_prob(s);
+ // Luma
+ std::vector<int> lens;
// need to reverse later
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);
+ 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 % BLOCKS_PER_STREAM == (BLOCKS_PER_STREAM - 1) || block_idx == NUM_BLOCKS - 1) {
+ int l = rans_encoder.save_block(codedfp);
+ num_bytes += l;
+ lens.push_back(l);
}
- if (yb % 16 == 8) {
+ }
+ 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);
+ }
+ }
+
+ // 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 % BLOCKS_PER_STREAM == (BLOCKS_PER_STREAM - 1) || block_idx == NUM_BLOCKS - 1) {
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 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 % BLOCKS_PER_STREAM == (BLOCKS_PER_STREAM - 1) || block_idx == NUM_BLOCKS - 1) {
+ num_bytes += rans_encoder.save_block(codedfp);
+ }
}
tot_bytes += num_bytes;
- printf("coeff %d: %ld bytes\n", y * 8 + x, num_bytes);
+ printf("coeff %d Cr: %ld bytes\n", y * 8 + x, num_bytes);
}
}
- 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);
+
+#if 0
+ printf("Max coefficient ranges (as a function of x):\n\n");
+ for (unsigned x = 0; x < 8; ++x) {
+ int range = std::max(max_val_x[x], -min_val_x[x]);
+ printf(" [%4d, %4d] (%.2f bits)\n", min_val_x[x], max_val_x[x], log2(range * 2 + 1));
+ }
+
+ printf("Max coefficient ranges (as a function of y):\n\n");
+ for (unsigned y = 0; y < 8; ++y) {
+ int range = std::max(max_val_y[y], -min_val_y[y]);
+ printf(" [%4d, %4d] (%.2f bits)\n", min_val_y[y], max_val_y[y], log2(range * 2 + 1));
+ }
+#endif
}
projects / narabu / blobdiff