/* Yliluoma's ordered dithering algorithm 1 d. Tri-tone dithering Usage: ./yliluoma_ordered_dither_1d.exe scene.png scenedither6.png Arbitrary-palette positional dithering algorithm https://bisqwit.iki.fi/story/howto/dither/jy/ */ #include #include #include /* 8x8 threshold map */ #define d(x) x / 64.0 static const double map[8 * 8] = { d(0), d(48), d(12), d(60), d(3), d(51), d(15), d(63), d(32), d(16), d(44), d(28), d(35), d(19), d(47), d(31), d(8), d(56), d(4), d(52), d(11), d(59), d(7), d(55), d(40), d(24), d(36), d(20), d(43), d(27), d(39), d(23), d(2), d(50), d(14), d(62), d(1), d(49), d(13), d(61), d(34), d(18), d(46), d(30), d(33), d(17), d(45), d(29), d(10), d(58), d(6), d(54), d(9), d(57), d(5), d(53), d(42), d(26), d(38), d(22), d(41), d(25), d(37), d(21)}; #undef d /* Palette */ static const unsigned pal[16] = {0x080000, 0x201A0B, 0x432817, 0x492910, 0x234309, 0x5D4F1E, 0x9C6B20, 0xA9220F, 0x2B347C, 0x2B7409, 0xD0CA40, 0xE8A077, 0x6A94AB, 0xD5C4B3, 0xFCE76E, 0xFCFAE2}; // Compare the difference of two RGB values, weigh by CCIR 601 luminosity: double ColorCompare(int r1, int g1, int b1, int r2, int g2, int b2) { double luma1 = (r1 * 299 + g1 * 587 + b1 * 114) / (255.0 * 1000); double luma2 = (r2 * 299 + g2 * 587 + b2 * 114) / (255.0 * 1000); double lumadiff = luma1 - luma2; double diffR = (r1 - r2) / 255.0, diffG = (g1 - g2) / 255.0, diffB = (b1 - b2) / 255.0; return (diffR * diffR * 0.299 + diffG * diffG * 0.587 + diffB * diffB * 0.114) * 0.75 + lumadiff * lumadiff; } double EvaluateMixingError(int r, int g, int b, // Desired color int r0, int g0, int b0, // Mathematical mix product int r1, int g1, int b1, // Mix component 1 int r2, int g2, int b2, // Mix component 2 double ratio) // Mixing ratio { return ColorCompare(r, g, b, r0, g0, b0) + ColorCompare(r1, g1, b1, r2, g2, b2) * 0.1 * (fabs(ratio - 0.5) + 0.5); } struct MixingPlan { unsigned colors[4]; double ratio; /* 0 = always index1, 1 = always index2, 0.5 = 50% of both */ /* 4 = special three or four-color dither */ }; MixingPlan DeviseBestMixingPlan(unsigned color) { const unsigned r = color >> 16, g = (color >> 8) & 0xFF, b = color & 0xFF; MixingPlan result = {{0, 0}, 0.5}; double least_penalty = 1e99; for (unsigned index1 = 0; index1 < 16; ++index1) for (unsigned index2 = index1; index2 < 16; ++index2) // for(int ratio=0; ratio<64; ++ratio) { // Determine the two component colors unsigned color1 = pal[index1], color2 = pal[index2]; unsigned r1 = color1 >> 16, g1 = (color1 >> 8) & 0xFF, b1 = color1 & 0xFF; unsigned r2 = color2 >> 16, g2 = (color2 >> 8) & 0xFF, b2 = color2 & 0xFF; int ratio = 32; if (color1 != color2) { // Determine the ratio of mixing for each channel. // solve(r1 + ratio*(r2-r1)/64 = r, ratio) // Take a weighed average of these three ratios according to the // perceived luminosity of each channel (according to CCIR 601). ratio = ((r2 != r1 ? 299 * 64 * int(r - r1) / int(r2 - r1) : 0) + (g2 != g1 ? 587 * 64 * int(g - g1) / int(g2 - g1) : 0) + (b1 != b2 ? 114 * 64 * int(b - b1) / int(b2 - b1) : 0)) / ((r2 != r1 ? 299 : 0) + (g2 != g1 ? 587 : 0) + (b2 != b1 ? 114 : 0)); if (ratio < 0) ratio = 0; else if (ratio > 63) ratio = 63; } // Determine what mixing them in this proportion will produce unsigned r0 = r1 + ratio * int(r2 - r1) / 64; unsigned g0 = g1 + ratio * int(g2 - g1) / 64; unsigned b0 = b1 + ratio * int(b2 - b1) / 64; double penalty = EvaluateMixingError( r, g, b, r0, g0, b0, r1, g1, b1, r2, g2, b2, ratio / double(64)); if (penalty < least_penalty) { least_penalty = penalty; result.colors[0] = index1; result.colors[1] = index2; result.ratio = ratio / double(64); } if (index1 != index2) for (unsigned index3 = 0; index3 < 16; ++index3) { if (index3 == index2 || index3 == index1) continue; // 50% index3, 25% index2, 25% index1 unsigned color3 = pal[index3]; unsigned r3 = color3 >> 16, g3 = (color3 >> 8) & 0xFF, b3 = color3 & 0xFF; r0 = (r1 + r2 + r3 * 2) / 4; g0 = (g1 + g2 + g3 * 2) / 4; b0 = (b1 + b2 + b3 * 2) / 4; penalty = ColorCompare(r, g, b, r0, g0, b0) + ColorCompare(r1, g1, b1, r2, g2, b2) * 0.025 + ColorCompare((r1 + g1) / 2, (g1 + g2) / 2, (b1 + b2) / 2, r3, g3, b3) * 0.025; if (penalty < least_penalty) { least_penalty = penalty; result.colors[0] = index3; // (0,0) index3 occurs twice result.colors[1] = index1; // (0,1) result.colors[2] = index2; // (1,0) result.colors[3] = index3; // (1,1) result.ratio = 4.0; } } } return result; } int main(int argc, char **argv) { if (argc != 3) { printf("Usage: %s INPUT.png OUTPUT.png\n", argv[0]); return 0; } FILE *fp = fopen(argv[1], "rb"); gdImagePtr srcim = gdImageCreateFromPng(fp); fclose(fp); unsigned w = gdImageSX(srcim), h = gdImageSY(srcim); gdImagePtr im = gdImageCreate(w, h); for (unsigned c = 0; c < 16; ++c) gdImageColorAllocate(im, pal[c] >> 16, (pal[c] >> 8) & 0xFF, pal[c] & 0xFF); #pragma omp parallel for for (unsigned y = 0; y < h; ++y) for (unsigned x = 0; x < w; ++x) { unsigned color = gdImageGetTrueColorPixel(srcim, x, y); MixingPlan plan = DeviseBestMixingPlan(color); if (plan.ratio == 4.0) // Tri-tone or quad-tone dithering { gdImageSetPixel(im, x, y, plan.colors[((y & 1) * 2 + (x & 1))]); } else { double map_value = map[(x & 7) + ((y & 7) << 3)]; gdImageSetPixel(im, x, y, plan.colors[map_value < plan.ratio ? 1 : 0]); } } fp = fopen(argv[2], "wb"); gdImagePng(im, fp); fclose(fp); gdImageDestroy(im); gdImageDestroy(srcim); }