/*
This file contains source code for

Yliluoma's ordered dithering algorithm 2

by Joel Yliluoma 2011.
Originally published at https://bisqwit.iki.fi/story/howto/dither/jy/

Modifications by Pekka Väänänen in 2026.

Original description
--------------------

An altogetherly different dithering algorithm can be devised by discarding the initial assumption that the dithering mixes two colortones together, and instead, assuming that each matrix value corresponds to a particular color tone. A 8x8 matrix has 64 color tones, a 2x2 matrix has 4 color tones, and so on.

An algorithm for populating such a color array will need to find the N-term expression of color values that, when combined, will produce the closest approximation of the input color.

One such algorithm is to start with a guess (the closest color), and then find out how much it went wrong, and then find out by experimentation which terms are needed to improve the result.

To solve the issue about pixel orientations changing, the colors will be sorted by luma. They will still change relative orientation, but such action is relatively minor.

In pseudo code, the process of converting the input bitmap into a target bitmap goes like this: 

  For each pixel, Input, in the original picture:
    Achieved      = 0      // Total color sum achieved so far
    CandidateList.clear()
    LoopWhile CandidateList.Size < (X * Y)
      Count     = 1
      Candidate = 0 // Candidate color from palette
      Comparison.reset()
      Max = CandidateList.Size, Or 1 if empty
      For each Color in palette:
        AddingCount = 1
        LoopWhile AddingCount <= Max
          Sum         = Achieved + Color * AddingCount
          Divide      = CandidateList.Size + AddingCount
          Test        = Sum / Divide
          Compare Test to Input using CIEDE2000 or RGB;
          If it was the best match since Comparison was reset:
            Candidate := Color
            Count     := AddingCount
          EndCompare
          AddingCount = AddingCount * 2 // Faster version
          // AddingCount = AddingCount + 1 // Slower version
        EndWhile
      CandidateList.Add(Candidate, Count times)
      Achieved = Achieved + Candidate * Count
    LoopEnd
    CandidateList.Sort( by: luminance )
    Index = ThresholdMatrix[xcoordinate % X][ycoordinate % Y]
    Draw pixel using CandidateList[Index * CandidateList.Size() / (X*Y)]
  EndFor
*/
#include <gd.h>
#include <stdio.h>
#include <math.h>
#include <algorithm>
#include <map>
#include <iostream>

#define USE_BAYER8x8 0
#define USE_ACCELERATED_SEARCH 0
#define USE_LUMA 0

#if USE_BAYER8x8
static const unsigned char map[8*8] = {
     0,48,12,60, 3,51,15,63,
    32,16,44,28,35,19,47,31,
     8,56, 4,52,11,59, 7,55,
    40,24,36,20,43,27,39,23,
     2,50,14,62, 1,49,13,61,
    34,18,46,30,33,17,45,29,
    10,58, 6,54, 9,57, 5,53,
    42,26,38,22,41,25,37,21 };
#else
static const unsigned char map[4*4] = {
     0,  8,  2, 10,
    12,  4, 14,  6,
     3, 11,  1, 9,
    15,  7, 13,  5 };
#endif

#if 1
// chronocross
static const unsigned pal[16] = { 
    0x080000,0x201A0B,0x432817,0x492910, 0x234309,0x5D4F1E,0x9C6B20,0xA9220F,
    0x2B347C,0x2B7409,0xD0CA40,0xE8A077, 0x6A94AB,0xD5C4B3,0xFCE76E,0xFCFAE2 };
#elif 0
static const unsigned pal[16] = // bigbuckbunny_bird
{
    0x8fa9ee,
    0x2c2e4b,
    0xffb801,
    0xd960c,
    0x994909,
    0xebd2cf,
    0x58e54,
    0x514b68,
    0xb88681,
    0x25f3e,
    0x686488,
    0xbdbae3,
    0x5e3314,
    0xe77906,
    0x8c5f62,
    0xe9ae83
};
#elif 0
static const unsigned pal[16] = // kodim21
{
    0x868373,
    0xfffed2,
    0x2d251f,
    0xc1b999,
    0x844838,
    0x859799,
    0x5a6060,
    0x9f8f6e,
    0xfaf1bc,
    0x87785c,
    0x443b2d,
    0x6a8699,
    0x5d503a,
    0xded5ad,
    0x73654b,
    0xaca489
};
#elif 0
static const unsigned pal[16] = { // kodim23
    0x69734b,
    0xefe6e9,
    0xf1bf06,
    0xc8ab9e,
    0xa2a4a5,
    0x3b392a,
    0xc74339,
    0x999f70,
    0xc09d13,
    0x56878d,
    0x83322a,
    0xd6c7c7,
    0x4f573e,
    0xebc224,
    0x83876e,
    0x557a2a
};

#endif


// { 0x080000,0x201A0B,0x432817,0x492910, 0x234309,0x5D4F1E,0x9C6B20,0xA9220F,
//   0x2B347C,0x2B7409,0xD0CA40,0xE8A077, 0x6A94AB,0xD5C4B3,0xFCE76E,0xFCFAE2 };

static unsigned luma[16];

bool PaletteCompareLuma(unsigned index1, unsigned index2)
{
    return luma[index1] < luma[index2];
}
#if USE_LUMA
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;
}
#else
double ColorCompare(int r1,int g1,int b1, int r2,int g2,int b2)
{
    double diffR = (r1-r2)/255.0, diffG = (g1-g2)/255.0, diffB = (b1-b2)/255.0;
    return diffR*diffR + diffG*diffG + diffB*diffB;
}
#endif
struct MixingPlan
{
    static constexpr unsigned n_colors = 16;
    unsigned colors[n_colors];
};

MixingPlan DeviseBestMixingPlan(unsigned color)
{
    MixingPlan result = {};
    const unsigned src[3] = { color>>16, (color>>8)&0xFF, color&0xFF };
    unsigned proportion_total = 0;
    
    unsigned so_far[3] = {0,0,0};

    while(proportion_total < MixingPlan::n_colors)
    {
        unsigned chosen_amount = 1;
        unsigned chosen        = 0;
        
        const unsigned max_test_count = std::max(1u, proportion_total);

        // Loop over palette colors.
        double least_penalty = -1;
        for(unsigned index=0; index<16; ++index)
        {
            const unsigned color_val = pal[index];
            unsigned sum[3] = { so_far[0], so_far[1], so_far[2] };
            unsigned add[3] = { color_val>>16, (color_val>>8)&0xFF, color_val&0xFF };

            #if USE_ACCELERATED_SEARCH
            for(unsigned p=1; p<=max_test_count; p*=2)
            {
                for(unsigned c=0; c<3; ++c) sum[c] += add[c];
                for(unsigned c=0; c<3; ++c) add[c] += add[c];
            #else
            for(unsigned p=1; p<=max_test_count; p++)
            {
                for(unsigned c=0; c<3; ++c) sum[c] += add[c];
            #endif
                unsigned t_count = proportion_total + p;
                unsigned test[3] = { sum[0] / t_count, sum[1] / t_count, sum[2] / t_count };
                double penalty = ColorCompare(src[0],src[1],src[2],
                                            test[0],test[1],test[2]);

                double t_val = (double)(2*p - 1) / (double)t_count;
                // std::cout << "N=" << proportion_total << " p=" << p << " t=" << t_val 
                //         << " test=" << test[0] << "," << test[1] << "," << test[2]
                //         << " penalty=" << penalty << std::endl;

                if(penalty < least_penalty || least_penalty < 0)
                {
                    least_penalty = penalty;
                    chosen        = index;
                    chosen_amount = p;
                }
            }
        }
        for(unsigned p=0; p<chosen_amount; ++p)
        {
            if(proportion_total >= MixingPlan::n_colors) break;
            result.colors[proportion_total++] = chosen;
        }

        const unsigned color_val = pal[chosen];
        unsigned palcolor[3] = { color_val>>16, (color_val>>8)&0xFF, color_val&0xFF };

        for(unsigned c=0; c<3; ++c)
            so_far[c] += palcolor[c] * chosen_amount;
    }
    std::sort(result.colors, result.colors+MixingPlan::n_colors, PaletteCompareLuma);
    return result;
}

int main(int argc, char**argv)
{
    FILE* fp = fopen(argv[1], "rb");
    printf("fp: %p\n", fp);
    gdImagePtr srcim = gdImageCreateFromPng(fp);
    printf("srcim: %p\n", srcim);
    fclose(fp);

    unsigned w = gdImageSX(srcim), h = gdImageSY(srcim);
    gdImagePtr im = gdImageCreate(w, h);
    for(unsigned c=0; c<16; ++c)
    {
        unsigned r = pal[c]>>16, g = (pal[c]>>8) & 0xFF, b = pal[c] & 0xFF;
        gdImageColorAllocate(im, r,g,b);
        luma[c] = r*299 + g*587 + b*114;
    }
//   #pragma omp parallel for
    for(unsigned y=0; y<h; ++y)
        for(unsigned x=0; x<w; ++x)
        {
            unsigned color = gdImageGetTrueColorPixel(srcim, x, y);
            #if USE_BAYER8x8
            unsigned map_value = map[(x & 7) + ((y & 7) << 3)];
            MixingPlan plan = DeviseBestMixingPlan(color);
            map_value = map_value * MixingPlan::n_colors / 64;
            #else
            unsigned map_value = map[(x & 3) + ((y & 3) << 2)];
            MixingPlan plan = DeviseBestMixingPlan(color);
            map_value = map_value * MixingPlan::n_colors / 16;
            #endif
            gdImageSetPixel(im, x,y, plan.colors[ map_value ]);
        }
    fp = fopen(argv[2], "wb");
    gdImagePng(im, fp);
    fclose(fp); gdImageDestroy(im); gdImageDestroy(srcim);
}