+}
+
+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();