X-Git-Url: https://git.sesse.net/?a=blobdiff_plain;f=qdc.cpp;h=556a8a5d541315a54d6bddf332cdd746a5c1769a;hb=34138d5b1a1302a7a1050fd46d2fb95c0186140a;hp=c4e1ba4e73f2615addb13f6ca09551e8bb9970e4;hpb=fb83fc30cf33cec1d155b3a63c338bbb64adb4e3;p=narabu diff --git a/qdc.cpp b/qdc.cpp index c4e1ba4..556a8a5 100644 --- a/qdc.cpp +++ b/qdc.cpp @@ -18,6 +18,11 @@ #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) @@ -25,7 +30,7 @@ // 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 8 +#define NUM_CLUSTERS 4 static constexpr uint32_t prob_bits = 12; static constexpr uint32_t prob_scale = 1 << prob_bits; @@ -154,12 +159,37 @@ SymbolStats stats[128]; 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 FIND_OPTIMAL_STREAM_ASSIGNMENT return y * 8 + x + is_chroma * 64; #else - return std::min(x + y, 7) + is_chroma * 8; + if (is_chroma) { + return chroma_mapping[y * 8 + x] + 4; + } else { + return luma_mapping[y * 8 + x]; + } #endif } @@ -187,7 +217,7 @@ public: printf("%d: cumfreqs=%d freqs=%d prob_bits=%d\n", i, s.cum_freqs[i], s.freqs[i], prob_bits + 1); RansEncSymbolInit(&esyms[i], s.cum_freqs[i], s.freqs[i], prob_bits + 1); } - sign_bias = s.cum_freqs[256]; + sign_bias = s.cum_freqs[NUM_SYMS]; } void clear() @@ -245,7 +275,7 @@ public: RansEncPut(&rans, &ptr, k, 1, prob_bits); k = ESCAPE_LIMIT; } - RansEncPutSymbol(&rans, &ptr, &esyms[(k - 1) & 255]); + RansEncPutSymbol(&rans, &ptr, &esyms[(k - 1) & (NUM_SYMS - 1)]); if (signed_k < 0) { rans += sign_bias; } @@ -374,75 +404,80 @@ double find_best_assignment(const int *medoids, int *assignment) return current_score; } -double find_inv_sum(const SymbolStats &stats) -{ - double s = 0.0; - for (unsigned j = 0; j < NUM_SYMS; ++j) { - s += stats.freqs[j] + 0.5; - } - return 1.0 / s; -} - void find_optimal_stream_assignment(int base) { - // k-means init; make random assignments - std::random_device rd; - std::mt19937 g(rd()); - std::uniform_int_distribution<> u(0, NUM_CLUSTERS - 1); - int assignment[64]; - for (unsigned i = 0; i < 64; ++i) { - assignment[i] = u(g); - } - double inv_sum_coeffs[64]; + double inv_sum[64]; for (unsigned i = 0; i < 64; ++i) { - inv_sum_coeffs[i] = find_inv_sum(stats[i + base]); + 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 iter = 0; iter < 1000; ++iter) { - // make new clusters based on the current assignments - SymbolStats clusters[NUM_CLUSTERS]; - for (unsigned i = 0; i < NUM_CLUSTERS; ++i) { - clusters[i].clear(); - } - for (unsigned i = 0; i < 64; ++i) { - for (unsigned j = 0; j < NUM_SYMS; ++j) { - clusters[assignment[i]].freqs[j] += stats[i + base].freqs[j]; + 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"); + } - double inv_sum_clusters[NUM_CLUSTERS]; - for (unsigned i = 0; i < NUM_CLUSTERS; ++i) { - inv_sum_clusters[i] = find_inv_sum(clusters[i]); - } - - // find new assignments based on distance to the clusters + // 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; - double total_d = 0.0; - for (unsigned i = 0; i < 64; ++i) { - int best_assignment = -1; - double best_distance = HUGE_VAL; - for (unsigned j = 0; j < NUM_CLUSTERS; ++j) { - double d = 0.0; - for (unsigned k = 0; k < NUM_SYMS; ++k) { - double p1 = (clusters[j].freqs[k] + 0.5) * inv_sum_clusters[j]; - double p2 = (stats[i + base].freqs[k] + 0.5) * inv_sum_coeffs[i]; - - // K-L divergence is asymmetric; this is a hack. - d += p1 * log(p1 / p2); - d += p2 * log(p2 / p1); - } - if (d < best_distance) { - best_assignment = j; - best_distance = d; + 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 (assignment[i] != best_assignment) { - any_changed = true; - } - assignment[i] = best_assignment; - total_d += best_distance; } - printf("iter %u: %.3f\n", iter, total_d); if (!any_changed) break; } printf("\n"); @@ -690,6 +725,12 @@ int main(int argc, char **argv) coeff_y[yb * WIDTH + xb] -= coeff_y[yb * WIDTH + (xb + 8)]; } } + 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); @@ -741,7 +782,7 @@ int main(int argc, char **argv) k = ESCAPE_LIMIT; extra_bits += 12; // escape this one } - ++s_luma.freqs[(k - 1) & 255]; + ++s_luma.freqs[(k - 1) & (NUM_SYMS - 1)]; } } // Chroma @@ -757,8 +798,8 @@ int main(int argc, char **argv) k_cr = ESCAPE_LIMIT; extra_bits += 12; // escape this one } - ++s_chroma.freqs[(k_cb - 1) & 255]; - ++s_chroma.freqs[(k_cr - 1) & 255]; + ++s_chroma.freqs[(k_cb - 1) & (NUM_SYMS - 1)]; + ++s_chroma.freqs[(k_cr - 1) & (NUM_SYMS - 1)]; } } } @@ -773,10 +814,10 @@ int main(int argc, char **argv) #endif for (unsigned i = 0; i < 64; ++i) { - stats[i].freqs[255] /= 2; // zero, has no sign bits (yes, this is trickery) + 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[256] += stats[i].freqs[255]; - stats[i].freqs[255] *= 2; + 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"); @@ -812,21 +853,20 @@ int main(int argc, char **argv) // 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_y[(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) { + 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); } } - if (HEIGHT % 16 != 0) { - num_bytes += rans_encoder.save_block(codedfp); - } tot_bytes += num_bytes; printf("coeff %d Y': %ld bytes\n", y * 8 + x, num_bytes); @@ -853,18 +893,18 @@ int main(int argc, char **argv) rans_encoder.clear(); size_t num_bytes = 0; - for (unsigned yb = 0; yb < HEIGHT; yb += 8) { - for (unsigned xb = 0; xb < WIDTH/2; xb += 8) { - int k = coeff_cb[(yb + y) * WIDTH/2 + (xb + x)]; - rans_encoder.encode_coeff(k); - } - if (yb % 16 == 8) { + 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 (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); } @@ -878,18 +918,18 @@ int main(int argc, char **argv) rans_encoder.clear(); size_t num_bytes = 0; - for (unsigned yb = 0; yb < HEIGHT; yb += 8) { - for (unsigned xb = 0; xb < WIDTH/2; xb += 8) { - int k = coeff_cr[(yb + y) * WIDTH/2 + (xb + x)]; - rans_encoder.encode_coeff(k); - } - if (yb % 16 == 8) { + 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); } } - if (HEIGHT % 16 != 0) { - num_bytes += rans_encoder.save_block(codedfp); - } tot_bytes += num_bytes; printf("coeff %d Cr: %ld bytes\n", y * 8 + x, num_bytes); }