/*
Copyright (c) 2014 Timur Gafarov

Boost Software License - Version 1.0 - August 17th, 2003

Permission is hereby granted, free of charge, to any person or organization
obtaining a copy of the software and accompanying documentation covered by
this license (the "Software") to use, reproduce, display, distribute,
execute, and transmit the Software, and to prepare derivative works of the
Software, and to permit third-parties to whom the Software is furnished to
do so, all subject to the following:

The copyright notices in the Software and this entire statement, including
the above license grant, this restriction and the following disclaimer,
must be included in all copies of the Software, in whole or in part, and
all derivative works of the Software, unless such copies or derivative
works are solely in the form of machine-executable object code generated by
a source language processor.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
DEALINGS IN THE SOFTWARE.
*/

// dimage is actually stripped out part of dlib - just to support reading PNG and JPEG
module dimage.idct; //dlib.image.io.idct

import std.math;

/*
 * Inverse discrete cosine transform (DCT) for 64x64 blocks
 */

enum blockSize = 64; // A DCT block is 8x8.

enum w1 = 2841; // 2048*sqrt(2)*cos(1*pi/16)
enum w2 = 2676; // 2048*sqrt(2)*cos(2*pi/16)
enum w3 = 2408; // 2048*sqrt(2)*cos(3*pi/16)
enum w5 = 1609; // 2048*sqrt(2)*cos(5*pi/16)
enum w6 = 1108; // 2048*sqrt(2)*cos(6*pi/16)
enum w7 = 565; // 2048*sqrt(2)*cos(7*pi/16)

enum w1pw7 = w1 + w7;
enum w1mw7 = w1 - w7;
enum w2pw6 = w2 + w6;
enum w2mw6 = w2 - w6;
enum w3pw5 = w3 + w5;
enum w3mw5 = w3 - w5;

enum r2 = 181; // 256/sqrt(2)

// idct performs a 2-D Inverse Discrete Cosine Transformation.
//
// The input coefficients should already have been multiplied by the
// appropriate quantization table. We use fixed-point computation, with the
// number of bits for the fractional component varying over the intermediate
// stages.
//
// For more on the actual algorithm, see Z. Wang, "Fast algorithms for the
// discrete W transform and for the discrete Fourier transform", IEEE Trans. on
// ASSP, Vol. ASSP- 32, pp. 803-816, Aug. 1984.
void idct64(int* src)
{
    // Horizontal 1-D IDCT.
    for (uint y = 0; y < 8; y++)
    {
        int y8 = y * 8;
        // If all the AC components are zero, then the IDCT is trivial.
        if (src[y8+1] == 0 && src[y8+2] == 0 && src[y8+3] == 0 &&
            src[y8+4] == 0 && src[y8+5] == 0 && src[y8+6] == 0 && src[y8+7] == 0)
        {
            int dc = src[y8+0] << 3;
            src[y8+0] = dc;
            src[y8+1] = dc;
            src[y8+2] = dc;
            src[y8+3] = dc;
            src[y8+4] = dc;
            src[y8+5] = dc;
            src[y8+6] = dc;
            src[y8+7] = dc;
            continue;
        }

        // Prescale.
        int x0 = (src[y8+0] << 11) + 128;
        int x1 = src[y8+4] << 11;
        int x2 = src[y8+6];
        int x3 = src[y8+2];
        int x4 = src[y8+1];
        int x5 = src[y8+7];
        int x6 = src[y8+5];
        int x7 = src[y8+3];

        // Stage 1.
        int x8 = w7 * (x4 + x5);
        x4 = x8 + w1mw7*x4;
        x5 = x8 - w1pw7*x5;
        x8 = w3 * (x6 + x7);
        x6 = x8 - w3mw5*x6;
        x7 = x8 - w3pw5*x7;

        // Stage 2.
        x8 = x0 + x1;
        x0 -= x1;
        x1 = w6 * (x3 + x2);
        x2 = x1 - w2pw6*x2;
        x3 = x1 + w2mw6*x3;
        x1 = x4 + x6;
        x4 -= x6;
        x6 = x5 + x7;
        x5 -= x7;

        // Stage 3.
        x7 = x8 + x3;
        x8 -= x3;
        x3 = x0 + x2;
        x0 -= x2;
        x2 = (r2*(x4+x5) + 128) >> 8;
        x4 = (r2*(x4-x5) + 128) >> 8;

        // Stage 4.
        src[y8+0] = (x7 + x1) >> 8;
        src[y8+1] = (x3 + x2) >> 8;
        src[y8+2] = (x0 + x4) >> 8;
        src[y8+3] = (x8 + x6) >> 8;
        src[y8+4] = (x8 - x6) >> 8;
        src[y8+5] = (x0 - x4) >> 8;
        src[y8+6] = (x3 - x2) >> 8;
        src[y8+7] = (x7 - x1) >> 8;
    }

    // Vertical 1-D IDCT.
    for (uint x = 0; x < 8; x++)
    {
        // Similar to the horizontal 1-D IDCT case, if all the AC components are zero, then the IDCT is trivial.
        // However, after performing the horizontal 1-D IDCT, there are typically non-zero AC components, so
        // we do not bother to check for the all-zero case.

        // Prescale.
        int y0 = (src[8*0+x] << 8) + 8192;
        int y1 = src[8*4+x] << 8;
        int y2 = src[8*6+x];
        int y3 = src[8*2+x];
        int y4 = src[8*1+x];
        int y5 = src[8*7+x];
        int y6 = src[8*5+x];
        int y7 = src[8*3+x];

        // Stage 1.
        int y8 = w7*(y4+y5) + 4;
        y4 = (y8 + w1mw7*y4) >> 3;
        y5 = (y8 - w1pw7*y5) >> 3;
        y8 = w3*(y6+y7) + 4;
        y6 = (y8 - w3mw5*y6) >> 3;
        y7 = (y8 - w3pw5*y7) >> 3;

        // Stage 2.
        y8 = y0 + y1;
        y0 -= y1;
        y1 = w6*(y3+y2) + 4;
        y2 = (y1 - w2pw6*y2) >> 3;
        y3 = (y1 + w2mw6*y3) >> 3;
        y1 = y4 + y6;
        y4 -= y6;
        y6 = y5 + y7;
        y5 -= y7;

        // Stage 3.
        y7 = y8 + y3;
        y8 -= y3;
        y3 = y0 + y2;
        y0 -= y2;
        y2 = (r2*(y4+y5) + 128) >> 8;
        y4 = (r2*(y4-y5) + 128) >> 8;

        // Stage 4.
        src[8*0+x] = (y7 + y1) >> 14;
        src[8*1+x] = (y3 + y2) >> 14;
        src[8*2+x] = (y0 + y4) >> 14;
        src[8*3+x] = (y8 + y6) >> 14;
        src[8*4+x] = (y8 - y6) >> 14;
        src[8*5+x] = (y0 - y4) >> 14;
        src[8*6+x] = (y3 - y2) >> 14;
        src[8*7+x] = (y7 - y1) >> 14;
    }
}

/*
void idct(float* inMat, float* outMat)
{
    uint i, j, u, v;
    float s;

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

        for (u = 0; u < 8; u++)
        for (v = 0; v < 8; v++)
        {
            s += inMat[u * 8 + v] 
                * cos((2 * i + 1) * u * PI / 16.0f)
                * cos((2 * j + 1) * v * PI / 16.0f)
                * ((u == 0) ? 1 / sqrt(2.0f) : 1.0f)
                * ((v == 0) ? 1 / sqrt(2.0f) : 1.0f);
        }

        outMat[i * 8 + j] = s / 4.0f;
    }
}
*/