#include<iostream>
#include <fstream>
#include <vector>
#include <array>
#include<cmath>
#include<cstdlib>

#if defined(_MSC_VER) && !defined(__INTEL_COMPILER)
#pragma warning(disable:4996)
#endif

using namespace std;
constexpr double pi = 3.1415926535897932384626433832;

constexpr int block_size = 8;
constexpr int qualtiy = 20;




#define _DO_QUANTIZATION
//#define _USE_YCBCR



const uint8_t ZigZagInv[8 * 8] = { 0, 1, 8,16, 9, 2, 3,10,

        17,24,32,25,18,11, 4, 5,
                  12,19,26,33,40,48,41,34,
                   27,20,13, 6, 7,14,21,28,
                   35,42,49,56,57,50,43,36,
                 29,22,15,23,30,37,44,51,
                   58,59,52,45,38,31,39,46,
                53,60,61,54,47,55,62,63 };
typedef  vector<vector< long long int>> map;

void readppm(const char* filename, int* height, int* width, map& r1, map& g1, map& b1, bool* vaild = nullptr)
{
    *vaild = false;
    FILE* fd;
    int  k;
    char c;
    long int i, j;
    char b[100];
    int red, green, blue;

    int n;
    int m;
    int bm;
    // map a1, b1, c1;

    fd = fopen(filename, "rb");
    if (fd == NULL)
    {
        printf("Could not open %s\n", filename);
        return;
    }
    else {
        *vaild = true;
    }
    c = getc(fd);
    if (c == 'P' || c == 'p')
        c = getc(fd);

    if (c == '3')
    {
        printf("%s is a PPM file (plain text version)\n", filename);

        // NOTE: This is not very good PPM code! Comments are not allowed
        // except immediately after the magic number.
        c = getc(fd);
        if (c == '\n' || c == '\r') // Skip any line break and comments
        {
            c = getc(fd);
            while (c == '#')
            {
                fscanf(fd, "%[^\n\r] ", b);
                //   printf("%s\n", b);
                c = getc(fd);
            }
            ungetc(c, fd);
        }
        fscanf(fd, "%d %d %d", &m, &n, &k);
        bm = fmaxf(n, m);
        //   r1.resize(bm);
      // /    g1.resize(bm);
       //    b1.resize(bm);

        map a1(bm, vector<long long>(bm));
        map t1(bm, vector<long long>(bm));
        map c1(bm, vector<long long>(bm));

        //   cout << a1.size() << endl;
           //  printf("%d rows  %d columns  max value= %d\n", n, m, k);



         ///  cout << m << " " << n << endl;

        for (i = 0; i < n; i++) {
            for (j = 0; j < m; j++) // Important bug fix here!
            { // i = row, j = column
                fscanf(fd, "%d %d %d", &red, &green, &blue);

                //   a1[j].resize(bm);
                  // g1[j].resize(bm);
                //   b1[j].resize(bm);

                a1[i][j] = (double(red) * 255.0 / double(k));
                t1[i][j] = (double(green) * 255.0 / double(k));
                c1[i][j] = (double(blue) * 255.0 / double(k));
            }
        }
        r1 = a1;
        g1 = t1;
        b1 = c1;
    }
    else
        if (c == '6')
        {
            printf("%s is a PPM file (raw version)!\n", filename);

            c = getc(fd);
            if (c == '\n' || c == '\r') // Skip any line break and comments
            {
                c = getc(fd);
                while (c == '#')
                {
                    fscanf(fd, "%[^\n\r] ", b);
                    printf("%s\n", b);
                    c = getc(fd);
                }
                ungetc(c, fd);
            }
            fscanf(fd, "%d %d %d", &m, &n, &k);
            bm = fmaxf(n, m);

            map a1(bm, vector<long long>(bm));
            map t1(bm, vector<long long>(bm));
            map c1(bm, vector<long long>(bm));

            //     printf("%d rows  %d columns  max value= %d\n", m, n, k);
            c = getc(fd); // Skip the last whitespace


            // Read and re-order as necessary
            for (i = 0; i < n; i++) {
                for (j = 0; j < m; j++) // Important bug fix here!
                {

                    /// r1[j].resize(n);
                  //   g1[j].resize(n);
                   //  b1[j].resize(n);


                    a1[i][j] = double(getc(fd));
                    t1[i][j] = double(getc(fd));
                    c1[i][j] = double(getc(fd));
                }
            }
            r1 = a1;
            g1 = t1;
            b1 = c1;
        }
        else
        {
            printf("%s is not a PPM file!\n", filename);
            return;
        }

    // printf("read image\n");

    *height = m;
    *width = n;
    //   return image;
}

void writeppm(string fileOut, int height, int width, map red, map green, map blue)
{

    ofstream out(fileOut);

    out << "P3" << endl;

    out << width << " " << height << endl;

    out << 255 << endl;

    for (int i = 0; i < height; ++i)
    {
        for (int j = 0; j < width; ++j)
        {
            out << abs(red[i][j]) << " " << abs(green[i][j]) << " " << abs(blue[i][j]) << endl;
        }
    }

    out.close();
}

array<array<double, block_size>, block_size> Compute8x8Dct(array<array<double, block_size>, block_size> in)
{
    int i, j, u, v;
    double s;

    array<array<double, block_size>, block_size> out;

    double cof5 = sqrt(2.0 * double(block_size));

    for (i = 0; i < block_size; i++)
        for (j = 0; j < block_size; j++)
        {
            s = 0;

            for (u = 0; u < block_size; u++)
                for (v = 0; v < block_size; v++)
                    s += in[u][v] * cos((2.0 * (double)u + 1.0) * (double)i * pi / ((double)block_size * 2.0)) *
                    cos((2.0 * (double)v + 1.0) * (double)j * pi / ((double)block_size * 2.0)) *
                    ((i == 0) ? 1.0 / sqrt(2) : 1.0) *
                    ((j == 0) ? 1.0 / sqrt(2) : 1.0);
            s = s / cof5;

            out[i][j] = s;
            //  cout << s / 4 << " " << out[i][j] << endl;
        }
    return out;
}

array<array<double, block_size>, block_size> Compute8x8Idct(array<array<double, block_size>, block_size> in)
{
    int i, j, u, v;
    double s;
    array<array<double, block_size>, block_size> out;
    double cof5 = sqrt(2.0 * double(block_size));

    for (i = 0; i < block_size; i++)
        for (j = 0; j < block_size; j++)
        {
            s = 0;

            for (u = 0; u < block_size; u++)
                for (v = 0; v < block_size; v++)
                    s += in[u][v] * cos((2.0 * (double)i + 1.0) * (double)u * pi / (double(block_size) * 2.0)) *
                    cos((2.0 * (double)j + 1.0) * (double)v * pi / (double(block_size) * 2.0)) *
                    ((u == 0) ? 1.0 / sqrtf(2) : 1.0) *
                    ((v == 0) ? 1.0 / sqrtf(2) : 1.0);

            s = s / cof5;

            if (s > 255)
            {
                s = 255;
            }

            out[i][j] = s;


        }

    return out;
}


array<array<double, block_size>, block_size> qu(array<array<double, block_size>, block_size> a, vector<double> quantization)

{
    ///std::vector<double> quantization = { 16 , 11 , 10 , 16 , 24 , 40 , 51 , 61 , 12 , 12 , 14 , 19 , 26 , 58 , 60 , 55 , 14 , 13 , 16 , 24 , 40 , 57 , 69 , 56 , 14 , 17 , 22 , 29 , 51 , 87 , 80 , 62 , 18 , 22 , 37 , 56 , 68 , 109 , 103 , 77 , 24 , 35 , 55 , 64 , 81 , 104 , 113 , 92 , 49 , 64 , 78 , 87 , 103 , 121 , 120 , 101 , 72 , 92 , 95 , 98 , 112 , 100 , 103 , 99 };

    for (int i = 0; i < block_size * block_size; ++i)
    {
        int col = i % block_size;
        int row = int(i / block_size);
        double olda = a[row][col];

        a[row][col] = round((a[row][col]) / (quantization[i]));

    }


    return a;
}



array<array<double, block_size>, block_size> qu2(array<array<double, block_size>, block_size> a, vector<double> quantization)

{

    for (int i = 0; i < block_size * block_size; ++i)
    {
        int col = i % block_size;
        int row = int(i / block_size);

        a[row][col] = ((a[row][col]) * (quantization[i]));


    }
    return a;
}

map jpg(map r, int height, int width, vector<double> quantization)
{
    vector<array<array<double, block_size>, block_size>> block = vector<array<array<double, block_size>, block_size>>();


    int wth = ((width * height) / (block_size * block_size));


    for (int i = 0; i < ((width * height) / (block_size * block_size)); ++i)
    {
        block.push_back(array<array<double, block_size>, block_size>());
    }

    // cout << " h\n";

    for (int i = 0; i < height / block_size; i++)
    {
        for (int j = 0; j < width / block_size; j++)
        {
            for (int k = i * block_size, l = 0; k < (i * block_size) + block_size; k++, l++)
            {
                for (int m = j * block_size, n = 0; m < (j * block_size) + block_size; m++, n++)
                {
                    block[(i * (width / block_size)) + j][n][l] = r[m][k];
                }

            }
        }
    }

    int bidx = block.size();


    for (int i = 0; i < bidx; ++i)
    {
        block[i] = Compute8x8Dct(block[i]);

    }
#ifdef _DO_QUANTIZATION
    for (int i = 0; i < bidx; ++i)
    {
        block[i] = qu(block[i], quantization);
    }
#endif

    for (int i = 0; i < bidx; ++i)
    {
        for (int j = 0; j < block_size * block_size; ++j)
        {
            //  int de = ZigZagInv[j];
              //block[bidx][(int)de / 8][de % 8] = block[bidx][(int)j / 8][j % 8];
        }

    }








    for (int i = 0; i < height / block_size; i++)
    {
        for (int j = 0; j < width / block_size; j++)
        {
            for (int k = i * block_size, l = 0; k < (i * block_size) + block_size; k++, l++)
            {
                for (int m = j * block_size, n = 0; m < (j * block_size) + block_size; m++, n++)
                {
                    r[m][k] = ((block[(i * (width / block_size)) + j][n][l]));
                }
            }
        }
    }


    cout << "done" << endl;

    return r;
}




map jpgde(map r, int height, int width, vector<double> quantization)
{

    //read

    vector<array<array<double, block_size>, block_size>> block = vector<array<array<double, block_size>, block_size>>();
    /// block.push_back(array<array<double, 8>, 8>());
    // block.resize(width * height / 64);


    int wth = ((width * height) / (block_size * block_size));


    for (int i = 0; i < wth; ++i)
    {
        block.push_back(array<array<double, block_size>, block_size>());
    }

    // cout << " h\n";

    for (int i = 0; i < height / block_size; i++)
    {
        for (int j = 0; j < width / block_size; j++)
        {
            for (int k = i * block_size, l = 0; k < (i * block_size) + block_size; k++, l++)
            {
                for (int m = j * block_size, n = 0; m < (j * block_size) + block_size; m++, n++)
                {
                    block[(i * (width / block_size)) + j][n][l] = r[m][k];
                }
            }
        }
    }

    int bidx = block.size();



#ifdef _DO_QUANTIZATION
    for (int i = 0; i < bidx; ++i)
    {
        block[i] = qu2(block[i], quantization);
    }
#endif // _DO_QUANTIZATION

    for (int i = 0; i < bidx; ++i)
    {
        block[i] = Compute8x8Idct(block[i]);
    }

    for (int i = 0; i < bidx; ++i)
    {
        for (int j = 0; j < 64; ++j)
        {
            //int de = ZigZagInv[j];
             //block[bidx][(int)de / 8][de % 8] = block[bidx][(int)j / 8][j % 8];
        }

    }

    for (int i = 0; i < bidx; ++i)
    {
        for (int j = 0; j < 64; ++j)
        {
            //  int de = ZigZagInv[j];
       ///s       block[bidx][(int)j / 8][j % 8] = block[bidx][(int)de / 8][de % 8];
        }

    }

    double val_scale = 1.0 / (double(block_size) / 8.0);

    for (int i = 0; i < height / block_size; i++)
    {
        for (int j = 0; j < width / block_size; j++)
        {
            for (int k = i * block_size, l = 0; k < (i * block_size) + block_size; k++, l++)
            {
                for (int m = j * block_size, n = 0; m < (j * block_size) + block_size; m++, n++)
                {
                    r[m][k] = ((block[(i * (width / block_size)) + j][n][l])) * val_scale;
                }
            }
        }
    }


    bool deb = 0;

    int bidx2 = 0;


    return r;
}





//void rgb2ycbr

void rgb2ycbcr(int width, int height, map& r, map& g, map& b)
{
    map r1 = r;
    map g1 = g;
    map b1 = b;

    for (int i = 0; i < height; ++i)
    {
        for (int j = 0; j < width; ++j)
        {
            r1[i][j] = 0.0 + (0.299 * r[i][j]) + (0.587 * g[i][j]) + (0.114 * b[i][j]);
            g1[i][j] = 128.0 - (0.168736 * r[i][j]) - (0.331264 * g[i][j]) + (0.5 * b[i][j]);
            b1[i][j] = 128.0 + (0.5 * r[i][j]) - (0.418688 * g[i][j]) - (0.081312 * b[i][j]);
        }
    }
    r = r1;
    g = g1;
    b = b1;
    //   Y = 0.257R´ + 0.504G´ + 0.098B´ + 16

   //Cb = -0.148R´ - 0.291G´ + 0.439B´ + 128

   //Cr = 0.439R´ - 0.368G´ - 0.071B´ + 128
}

void ycbcr2rgb(int width, int height, map& y, map& cb, map& cr)
{
    map r = y;
    map g = cb;
    map b = cr;

    for (int i = 0; i < height; ++i)
    {
        for (int j = 0; j < width; ++j)
        {
            r[i][j] = (y[i][j]) + 1.402 * (cr[i][j] - 128);
            g[i][j] = (y[i][j]) - 0.344136 * (cb[i][j] - 128) - 0.714136 * (cr[i][j] - 128);
            b[i][j] = (y[i][j]) + 1.772 * (cb[i][j] - 128);
        }
    }
    y = r;
    cb = g;
    cr = b;
    // R´ = 1.164(Y - 16) + 1.596(Cr - 128)

      //   G´ = 1.164(Y - 16) - 0.813(Cr - 128) - 0.392(Cb - 128)

        // B´ = 1.164(Y - 16) + 2.017(Cb - 128)
}



int main()
{
    map r;
    map g;
    map b;

    int width = 0;

    int height = 0;
    bool vaild = false;
    double S = 0;
    vector<double> tb = { 16, 11, 10, 16,24, 40, 51, 61, 12, 12, 14, 19, 26, 58, 60, 55,
        14 ,13, 16, 24, 40, 57, 69, 56, 14, 17, 22, 29, 51, 87, 80, 62,
        18,22 ,37, 56, 68, 109, 103, 77, 24, 35, 55, 64, 81, 104, 113, 92,
        49 ,64, 78 ,87 ,103 ,121, 120, 101, 72, 92 ,95, 98 ,112, 100, 103, 99 };

    ///  %// Determine S
    if (qualtiy < 50)
        S = 5000 / qualtiy;
    else
        S = 200 - 2 * qualtiy;
    //end

    for (int i = 0; i < 64; ++i)
    {
        tb[i] = floor((S * tb[i] + 50.0) / 100.0);
        if (tb[i] == 0)
        {
            tb[i] = 1;
        }
    }
    // Ts = floor((S * Tb + 50) / 100);


// Ts(Ts == 0) = 1;

 //read image
    readppm("test3.ppm", &width, &height, r, g, b, &vaild);
    map c = r;
    if (!vaild)
    {
        return EXIT_FAILURE;

    }

    //do stuff here








#ifdef _USE_YCBCR
    rgb2ycbcr(width, height, r, g, b);
#endif // _USE_YCBCR

    //write to output image

    r = jpg(r, width, height, tb);
    r = jpgde(r, width, height, tb);

    g = jpg(g, width, height, tb);
    g = jpgde(g, width, height, tb);

    b = jpg(b, width, height, tb);
    b = jpgde(b, width, height, tb);

#ifdef _USE_YCBCR
    ycbcr2rgb(width, height, r, g, b);
#endif // _USE_YCBCR


    writeppm("output.ppm", height, width, r, g, b);




    cout << width << "  " << height << endl;





}