/*
 *  Copyright (c) 2010 The WebM project authors. All Rights Reserved.
 *
 *  Use of this source code is governed by a BSD-style license
 *  that can be found in the LICENSE file in the root of the source
 *  tree. An additional intellectual property rights grant can be found
 *  in the file PATENTS.  All contributing project authors may
 *  be found in the AUTHORS file in the root of the source tree.
 */


#include <float.h>
#include <math.h>
#include <stdio.h>
#include "vpx_mem/vpx_mem.h"
#include "vpxscale_arbitrary.h"

#define FIXED_POINT

#define MAX_IN_WIDTH        800
#define MAX_IN_HEIGHT       600
#define MAX_OUT_WIDTH       800
#define MAX_OUT_HEIGHT      600
#define MAX_OUT_DIMENSION   ((MAX_OUT_WIDTH > MAX_OUT_HEIGHT) ? \
                             MAX_OUT_WIDTH : MAX_OUT_HEIGHT)

BICUBIC_SCALER_STRUCT g_b_scaler;
static int g_first_time = 1;

#pragma DATA_SECTION(g_hbuf, "VP6_HEAP")
#pragma DATA_ALIGN (g_hbuf, 32);
unsigned char g_hbuf[MAX_OUT_DIMENSION];

#pragma DATA_SECTION(g_hbuf_uv, "VP6_HEAP")
#pragma DATA_ALIGN (g_hbuf_uv, 32);
unsigned char g_hbuf_uv[MAX_OUT_DIMENSION];


#ifdef FIXED_POINT
static int a_i = 0.6 * 65536;
#else
static float a = -0.6;
#endif

#ifdef FIXED_POINT
//         3     2
// C0 = a*t - a*t
//
static short c0_fixed(unsigned int t)
{
    // put t in Q16 notation
    unsigned short v1, v2;

    // Q16
    v1 = (a_i * t) >> 16;
    v1 = (v1 * t) >> 16;

    // Q16
    v2 = (a_i * t) >> 16;
    v2 = (v2 * t) >> 16;
    v2 = (v2 * t) >> 16;

    // Q12
    return -((v1 - v2) >> 4);
}

//                     2          3
// C1 = a*t + (3-2*a)*t  - (2-a)*t
//
static short c1_fixed(unsigned int t)
{
    unsigned short v1, v2, v3;
    unsigned short two, three;

    // Q16
    v1 = (a_i * t) >> 16;

    // Q13
    two = 2 << 13;
    v2 = two - (a_i >> 3);
    v2 = (v2 * t) >> 16;
    v2 = (v2 * t) >> 16;
    v2 = (v2 * t) >> 16;

    // Q13
    three = 3 << 13;
    v3 = three - (2 * (a_i >> 3));
    v3 = (v3 * t) >> 16;
    v3 = (v3 * t) >> 16;

    // Q12
    return (((v1 >> 3) - v2 + v3) >> 1);

}

//                 2          3
// C2 = 1 - (3-a)*t  + (2-a)*t
//
static short c2_fixed(unsigned int t)
{
    unsigned short v1, v2, v3;
    unsigned short two, three;

    // Q13
    v1 = 1 << 13;

    // Q13
    three = 3 << 13;
    v2 = three - (a_i >> 3);
    v2 = (v2 * t) >> 16;
    v2 = (v2 * t) >> 16;

    // Q13
    two = 2 << 13;
    v3 = two - (a_i >> 3);
    v3 = (v3 * t) >> 16;
    v3 = (v3 * t) >> 16;
    v3 = (v3 * t) >> 16;

    // Q12
    return (v1 - v2 + v3) >> 1;
}

//                 2      3
// C3 = a*t - 2*a*t  + a*t
//
static short c3_fixed(unsigned int t)
{
    int v1, v2, v3;

    // Q16
    v1 = (a_i * t) >> 16;

    // Q15
    v2 = 2 * (a_i >> 1);
    v2 = (v2 * t) >> 16;
    v2 = (v2 * t) >> 16;

    // Q16
    v3 = (a_i * t) >> 16;
    v3 = (v3 * t) >> 16;
    v3 = (v3 * t) >> 16;

    // Q12
    return ((v2 - (v1 >> 1) - (v3 >> 1)) >> 3);
}
#else
//          3     2
// C0 = -a*t + a*t
//
float C0(float t)
{
    return -a * t * t * t + a * t * t;
}

//                      2          3
// C1 = -a*t + (2*a+3)*t  - (a+2)*t
//
float C1(float t)
{
    return -(a + 2.0f) * t * t * t + (2.0f * a + 3.0f) * t * t - a * t;
}

//                 2          3
// C2 = 1 - (a+3)*t  + (a+2)*t
//
float C2(float t)
{
    return (a + 2.0f) * t * t * t - (a + 3.0f) * t * t + 1.0f;
}

//                 2      3
// C3 = a*t - 2*a*t  + a*t
//
float C3(float t)
{
    return a * t * t * t - 2.0f * a * t * t + a * t;
}
#endif

#if 0
int compare_real_fixed()
{
    int i, errors = 0;
    float mult = 1.0 / 10000.0;
    unsigned int fixed_mult = mult * 4294967296;//65536;
    unsigned int phase_offset_int;
    float phase_offset_real;

    for (i = 0; i < 10000; i++)
    {
        int fixed0, fixed1, fixed2, fixed3, fixed_total;
        int real0, real1, real2, real3, real_total;

        phase_offset_real = (float)i * mult;
        phase_offset_int = (fixed_mult * i) >> 16;
//      phase_offset_int = phase_offset_real * 65536;

        fixed0 = c0_fixed(phase_offset_int);
        real0 = C0(phase_offset_real) * 4096.0;

        if ((abs(fixed0) > (abs(real0) + 1)) || (abs(fixed0) < (abs(real0) - 1)))
            errors++;

        fixed1 = c1_fixed(phase_offset_int);
        real1 = C1(phase_offset_real) * 4096.0;

        if ((abs(fixed1) > (abs(real1) + 1)) || (abs(fixed1) < (abs(real1) - 1)))
            errors++;

        fixed2 = c2_fixed(phase_offset_int);
        real2 = C2(phase_offset_real) * 4096.0;

        if ((abs(fixed2) > (abs(real2) + 1)) || (abs(fixed2) < (abs(real2) - 1)))
            errors++;

        fixed3 = c3_fixed(phase_offset_int);
        real3 = C3(phase_offset_real) * 4096.0;

        if ((abs(fixed3) > (abs(real3) + 1)) || (abs(fixed3) < (abs(real3) - 1)))
            errors++;

        fixed_total = fixed0 + fixed1 + fixed2 + fixed3;
        real_total = real0 + real1 + real2 + real3;

        if ((fixed_total > 4097) || (fixed_total < 4094))
            errors ++;

        if ((real_total > 4097) || (real_total < 4095))
            errors ++;
    }

    return errors;
}
#endif

// Find greatest common denominator between two integers.  Method used here is
//  slow compared to Euclid's algorithm, but does not require any division.
int gcd(int a, int b)
{
    // Problem with this algorithm is that if a or b = 0 this function
    //  will never exit.  Don't want to return 0 because any computation
    //  that was based on a common denoninator and tried to reduce by
    //  dividing by 0 would fail.  Best solution that could be thought of
    //  would to be fail by returing a 1;
    if (a <= 0 || b <= 0)
        return 1;

    while (a != b)
    {
        if (b > a)
            b = b - a;
        else
        {
            int tmp = a;//swap large and
            a = b; //small
            b = tmp;
        }
    }

    return b;
}

void bicubic_coefficient_init()
{
    vpx_memset(&g_b_scaler, 0, sizeof(BICUBIC_SCALER_STRUCT));
    g_first_time = 0;
}

void bicubic_coefficient_destroy()
{
    if (!g_first_time)
    {
        vpx_free(g_b_scaler.l_w);

        vpx_free(g_b_scaler.l_h);

        vpx_free(g_b_scaler.l_h_uv);

        vpx_free(g_b_scaler.c_w);

        vpx_free(g_b_scaler.c_h);

        vpx_free(g_b_scaler.c_h_uv);

        vpx_memset(&g_b_scaler, 0, sizeof(BICUBIC_SCALER_STRUCT));
    }
}

// Create the coeffients that will be used for the cubic interpolation.
//  Because scaling does not have to be equal in the vertical and horizontal
//  regimes the phase offsets will be different.  There are 4 coefficents
//  for each point, two on each side.  The layout is that there are the
//  4 coefficents for each phase in the array and then the next phase.
int bicubic_coefficient_setup(int in_width, int in_height, int out_width, int out_height)
{
    int i;
#ifdef FIXED_POINT
    int phase_offset_int;
    unsigned int fixed_mult;
    int product_val = 0;
#else
    float phase_offset;
#endif
    int gcd_w, gcd_h, gcd_h_uv, d_w, d_h, d_h_uv;

    if (g_first_time)
        bicubic_coefficient_init();


    // check to see if the coefficents have already been set up correctly
    if ((in_width == g_b_scaler.in_width) && (in_height == g_b_scaler.in_height)
        && (out_width == g_b_scaler.out_width) && (out_height == g_b_scaler.out_height))
        return 0;

    g_b_scaler.in_width = in_width;
    g_b_scaler.in_height = in_height;
    g_b_scaler.out_width = out_width;
    g_b_scaler.out_height = out_height;

    // Don't want to allow crazy scaling, just try and prevent a catastrophic
    //  failure here.  Want to fail after setting the member functions so if
    //  if the scaler is called the member functions will not scale.
    if (out_width <= 0 || out_height <= 0)
        return -1;

    // reduce in/out width and height ratios using the gcd
    gcd_w = gcd(out_width, in_width);
    gcd_h = gcd(out_height, in_height);
    gcd_h_uv = gcd(out_height, in_height / 2);

    // the numerator width and height are to be saved in
    //  globals so they can be used during the scaling process
    //  without having to be recalculated.
    g_b_scaler.nw = out_width / gcd_w;
    d_w = in_width / gcd_w;

    g_b_scaler.nh = out_height / gcd_h;
    d_h = in_height / gcd_h;

    g_b_scaler.nh_uv = out_height / gcd_h_uv;
    d_h_uv = (in_height / 2) / gcd_h_uv;

    // allocate memory for the coefficents
    vpx_free(g_b_scaler.l_w);

    vpx_free(g_b_scaler.l_h);

    vpx_free(g_b_scaler.l_h_uv);

    g_b_scaler.l_w = (short *)vpx_memalign(32, out_width * 2);
    g_b_scaler.l_h = (short *)vpx_memalign(32, out_height * 2);
    g_b_scaler.l_h_uv = (short *)vpx_memalign(32, out_height * 2);

    vpx_free(g_b_scaler.c_w);

    vpx_free(g_b_scaler.c_h);

    vpx_free(g_b_scaler.c_h_uv);

    g_b_scaler.c_w = (short *)vpx_memalign(32, g_b_scaler.nw * 4 * 2);
    g_b_scaler.c_h = (short *)vpx_memalign(32, g_b_scaler.nh * 4 * 2);
    g_b_scaler.c_h_uv = (short *)vpx_memalign(32, g_b_scaler.nh_uv * 4 * 2);

    g_b_scaler.hbuf = g_hbuf;
    g_b_scaler.hbuf_uv = g_hbuf_uv;

    // Set up polyphase filter taps.  This needs to be done before
    //  the scaling because of the floating point math required.  The
    //  coefficients are multiplied by 2^12 so that fixed point math
    //  can be used in the main scaling loop.
#ifdef FIXED_POINT
    fixed_mult = (1.0 / (float)g_b_scaler.nw) * 4294967296;

    product_val = 0;

    for (i = 0; i < g_b_scaler.nw; i++)
    {
        if (product_val > g_b_scaler.nw)
            product_val -= g_b_scaler.nw;

        phase_offset_int = (fixed_mult * product_val) >> 16;

        g_b_scaler.c_w[i*4]   = c3_fixed(phase_offset_int);
        g_b_scaler.c_w[i*4+1] = c2_fixed(phase_offset_int);
        g_b_scaler.c_w[i*4+2] = c1_fixed(phase_offset_int);
        g_b_scaler.c_w[i*4+3] = c0_fixed(phase_offset_int);

        product_val += d_w;
    }


    fixed_mult = (1.0 / (float)g_b_scaler.nh) * 4294967296;

    product_val = 0;

    for (i = 0; i < g_b_scaler.nh; i++)
    {
        if (product_val > g_b_scaler.nh)
            product_val -= g_b_scaler.nh;

        phase_offset_int = (fixed_mult * product_val) >> 16;

        g_b_scaler.c_h[i*4]   = c0_fixed(phase_offset_int);
        g_b_scaler.c_h[i*4+1] = c1_fixed(phase_offset_int);
        g_b_scaler.c_h[i*4+2] = c2_fixed(phase_offset_int);
        g_b_scaler.c_h[i*4+3] = c3_fixed(phase_offset_int);

        product_val += d_h;
    }

    fixed_mult = (1.0 / (float)g_b_scaler.nh_uv) * 4294967296;

    product_val = 0;

    for (i = 0; i < g_b_scaler.nh_uv; i++)
    {
        if (product_val > g_b_scaler.nh_uv)
            product_val -= g_b_scaler.nh_uv;

        phase_offset_int = (fixed_mult * product_val) >> 16;

        g_b_scaler.c_h_uv[i*4]   = c0_fixed(phase_offset_int);
        g_b_scaler.c_h_uv[i*4+1] = c1_fixed(phase_offset_int);
        g_b_scaler.c_h_uv[i*4+2] = c2_fixed(phase_offset_int);
        g_b_scaler.c_h_uv[i*4+3] = c3_fixed(phase_offset_int);

        product_val += d_h_uv;
    }

#else

    for (i = 0; i < g_nw; i++)
    {
        phase_offset = (float)((i * d_w) % g_nw) / (float)g_nw;
        g_c_w[i*4]   = (C3(phase_offset) * 4096.0);
        g_c_w[i*4+1] = (C2(phase_offset) * 4096.0);
        g_c_w[i*4+2] = (C1(phase_offset) * 4096.0);
        g_c_w[i*4+3] = (C0(phase_offset) * 4096.0);
    }

    for (i = 0; i < g_nh; i++)
    {
        phase_offset = (float)((i * d_h) % g_nh) / (float)g_nh;
        g_c_h[i*4]   = (C0(phase_offset) * 4096.0);
        g_c_h[i*4+1] = (C1(phase_offset) * 4096.0);
        g_c_h[i*4+2] = (C2(phase_offset) * 4096.0);
        g_c_h[i*4+3] = (C3(phase_offset) * 4096.0);
    }

    for (i = 0; i < g_nh_uv; i++)
    {
        phase_offset = (float)((i * d_h_uv) % g_nh_uv) / (float)g_nh_uv;
        g_c_h_uv[i*4]   = (C0(phase_offset) * 4096.0);
        g_c_h_uv[i*4+1] = (C1(phase_offset) * 4096.0);
        g_c_h_uv[i*4+2] = (C2(phase_offset) * 4096.0);
        g_c_h_uv[i*4+3] = (C3(phase_offset) * 4096.0);
    }

#endif

    // Create an array that corresponds input lines to output lines.
    //  This doesn't require floating point math, but it does require
    //  a division and because hardware division is not present that
    //  is a call.
    for (i = 0; i < out_width; i++)
    {
        g_b_scaler.l_w[i] = (i * d_w) / g_b_scaler.nw;

        if ((g_b_scaler.l_w[i] + 2) <= in_width)
            g_b_scaler.max_usable_out_width = i;

    }

    for (i = 0; i < out_height + 1; i++)
    {
        g_b_scaler.l_h[i] = (i * d_h) / g_b_scaler.nh;
        g_b_scaler.l_h_uv[i] = (i * d_h_uv) / g_b_scaler.nh_uv;
    }

    return 0;
}

int bicubic_scale(int in_width, int in_height, int in_stride,
                  int out_width, int out_height, int out_stride,
                  unsigned char *input_image, unsigned char *output_image)
{
    short *RESTRICT l_w, * RESTRICT l_h;
    short *RESTRICT c_w, * RESTRICT c_h;
    unsigned char *RESTRICT ip, * RESTRICT op;
    unsigned char *RESTRICT hbuf;
    int h, w, lw, lh;
    int temp_sum;
    int phase_offset_w, phase_offset_h;

    c_w = g_b_scaler.c_w;
    c_h = g_b_scaler.c_h;

    op = output_image;

    l_w = g_b_scaler.l_w;
    l_h = g_b_scaler.l_h;

    phase_offset_h = 0;

    for (h = 0; h < out_height; h++)
    {
        // select the row to work on
        lh = l_h[h];
        ip = input_image + (in_stride * lh);

        // vp8_filter the row vertically into an temporary buffer.
        //  If the phase offset == 0 then all the multiplication
        //  is going to result in the output equalling the input.
        //  So instead point the temporary buffer to the input.
        //  Also handle the boundry condition of not being able to
        //  filter that last lines.
        if (phase_offset_h && (lh < in_height - 2))
        {
            hbuf = g_b_scaler.hbuf;

            for (w = 0; w < in_width; w++)
            {
                temp_sum =  c_h[phase_offset_h*4+3] * ip[w - in_stride];
                temp_sum += c_h[phase_offset_h*4+2] * ip[w];
                temp_sum += c_h[phase_offset_h*4+1] * ip[w + in_stride];
                temp_sum += c_h[phase_offset_h*4]   * ip[w + 2*in_stride];

                hbuf[w] = temp_sum >> 12;
            }
        }
        else
            hbuf = ip;

        // increase the phase offset for the next time around.
        if (++phase_offset_h >= g_b_scaler.nh)
            phase_offset_h = 0;

        // now filter and expand it horizontally into the final
        //  output buffer
        phase_offset_w = 0;

        for (w = 0; w < out_width; w++)
        {
            // get the index to use to expand the image
            lw = l_w[w];

            temp_sum =  c_w[phase_offset_w*4]   * hbuf[lw - 1];
            temp_sum += c_w[phase_offset_w*4+1] * hbuf[lw];
            temp_sum += c_w[phase_offset_w*4+2] * hbuf[lw + 1];
            temp_sum += c_w[phase_offset_w*4+3] * hbuf[lw + 2];
            temp_sum = temp_sum >> 12;

            if (++phase_offset_w >= g_b_scaler.nw)
                phase_offset_w = 0;

            // boundry conditions
            if ((lw + 2) >= in_width)
                temp_sum = hbuf[lw];

            if (lw == 0)
                temp_sum = hbuf[0];

            op[w] = temp_sum;
        }

        op += out_stride;
    }

    return 0;
}

void bicubic_scale_frame_reset()
{
    g_b_scaler.out_width = 0;
    g_b_scaler.out_height = 0;
}

void bicubic_scale_frame(YV12_BUFFER_CONFIG *src, YV12_BUFFER_CONFIG *dst,
                         int new_width, int new_height)
{

    dst->y_width = new_width;
    dst->y_height = new_height;
    dst->uv_width = new_width / 2;
    dst->uv_height = new_height / 2;

    dst->y_stride = dst->y_width;
    dst->uv_stride = dst->uv_width;

    bicubic_scale(src->y_width, src->y_height, src->y_stride,
                  new_width, new_height, dst->y_stride,
                  src->y_buffer, dst->y_buffer);

    bicubic_scale(src->uv_width, src->uv_height, src->uv_stride,
                  new_width / 2, new_height / 2, dst->uv_stride,
                  src->u_buffer, dst->u_buffer);

    bicubic_scale(src->uv_width, src->uv_height, src->uv_stride,
                  new_width / 2, new_height / 2, dst->uv_stride,
                  src->v_buffer, dst->v_buffer);
}