#include <stdio.h>
#include <pslib.h>
#include "tap.h"
#include "pstap.h"


#define NUM_POINTS 32

// Function to define mask values
psMaskType maskFunc(float x, float y)
{
    return 0;
}

// Function to define variance values
double varianceFunc(float x, float y)
{
    return sqrtf(10.0 + 0.1 * x + 0.2 * y);
}

// Function to define image values
double imageFunc(float x, float y)
{
    return 0.01 * x + 0.02 * y;
}

// Generate a dummy mask image
psImage *generateMask(int numCols, int numRows)
{
    psImage *image = psImageAlloc(numCols, numRows, PS_TYPE_MASK);
    for (int j = 0; j < numRows; j++) {
        for (int i = 0; i < numCols; i++) {
            psImageSet(image, i, j, maskFunc(i + 0.5, j + 0.5));
        }
    }
    return image;
}

// Generate a dummy variance image
psImage *generateVariance(int numCols, int numRows, psElemType type)
{
    psImage *image = psImageAlloc(numCols, numRows, type);
    for (int j = 0; j < numRows; j++) {
        for (int i = 0; i < numCols; i++) {
            psImageSet(image, i, j, varianceFunc(i + 0.5, j + 0.5));
        }
    }
    return image;
}

// Generate a dummy image
psImage *generateImage(int numCols, int numRows, psElemType type)
{
    psImage *image = psImageAlloc(numCols, numRows, type);
    for (int j = 0; j < numRows; j++) {
        for (int i = 0; i < numCols; i++) {
            psImageSet(image, i, j, imageFunc(i + 0.5, j + 0.5));
        }
    }
    return image;
}

// Check a particular combination of image and interpolation parameters
// Each has 2 + 2*num tests
void check(int xSize, int ySize,        // Size of image
           psElemType type,             // Type for image
           int num,                     // Number of samples over image
           psImageInterpolateMode mode, // Interpolation mode
           int xKernel, int yKernel,    // Size of interpolation kernel
           float tol                    // Tolerance
    )
{

    psMemId id = psMemGetLastId();

    psImage *image = generateImage(xSize, ySize, type);
    psImage *variance = NULL; // generateVariance(xSize, ySize, type);
    psImage *mask = NULL; // generateMask(xSize, ySize);
    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 12345);

    psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(mode, image,
                                                                       variance, mask, 0, NAN, NAN,
                                                                       0, 0, 0.0);
    ok(interp, "Interpolation options set");
    skip_start(!interp, 2 * num, "Interpolation options failed");

    for (int i = 0; i < num; i++) {
        float x = psRandomUniform(rng) * xSize;
        float y = psRandomUniform(rng) * ySize;

        // Values from interpolation
        double imageVal, varianceVal;
        psMaskType maskVal;

        bool success = psImageInterpolate(&imageVal, &varianceVal, &maskVal, x, y, interp);
        ok(success, "Interpolation at %.2f,%.2f", x, y);
        skip_start(!success, 1, "Interpolation failed");

        int xCentral, yCentral;         // Central pixel of interpolation
        switch (mode) {
          case PS_INTERPOLATE_BICUBE:
          case PS_INTERPOLATE_GAUSS:
            xCentral = x;
            yCentral = y;
            break;
          default:
            xCentral = floor(x - 0.5);
            yCentral = floor(y - 0.5);
            break;
        }

        if (xCentral - (xKernel - 1) / 2 < 0 || xCentral + xKernel / 2 > xSize - 1 ||
            yCentral - (yKernel - 1) / 2 < 0 || yCentral + yKernel / 2 > ySize - 1) {
            ok(isnan(imageVal), "Interpolation = %f vs NAN (border)", imageVal);
        } else {
            is_double_tol(imageVal, imageFunc(x, y), tol, "Interpolation = %f vs %f",
                          imageVal, imageFunc(x, y));
        }

        // Variance is very very very difficult to test without having the answer already
        // Mask might be a bit easier.
        // Defer these for now.

        skip_end();
    }

    psFree(interp);
    skip_end();

    psFree(rng);
    psFree(image);
    psFree(variance);
    psFree(mask);

    ok(!psMemCheckLeaks(id, NULL, NULL, false), "No memory leaks");

    return;
}

int main(int argc, char *argv[])
{
    plan_tests(792);

    check(16, 16, PS_TYPE_F32, 32, PS_INTERPOLATE_BILINEAR, 2, 2, 1.0e-6); // 66 tests
    check(16, 16, PS_TYPE_F64, 32, PS_INTERPOLATE_BILINEAR, 2, 2, 1.0e-6); // 66 tests

    check(16, 16, PS_TYPE_F32, 32, PS_INTERPOLATE_BICUBE, 3, 3, 1.0e-6); // 66 tests
    check(16, 16, PS_TYPE_F64, 32, PS_INTERPOLATE_BICUBE, 3, 3, 1.0e-6); // 66 tests

    check(16, 16, PS_TYPE_F32, 32, PS_INTERPOLATE_GAUSS, 3, 3, 1.0e-3); // 66 tests
    check(16, 16, PS_TYPE_F64, 32, PS_INTERPOLATE_GAUSS, 3, 3, 1.0e-3); // 66 tests

    check(16, 16, PS_TYPE_F32, 32, PS_INTERPOLATE_LANCZOS2, 4, 4, 1.0e-3); // 66 tests
    check(16, 16, PS_TYPE_F64, 32, PS_INTERPOLATE_LANCZOS2, 4, 4, 1.0e-3); // 66 tests

    check(32, 32, PS_TYPE_F32, 32, PS_INTERPOLATE_LANCZOS3, 6, 6, 1.0e-3); // 66 tests
    check(32, 32, PS_TYPE_F64, 32, PS_INTERPOLATE_LANCZOS3, 6, 6, 1.0e-3); // 66 tests

    check(32, 32, PS_TYPE_F32, 32, PS_INTERPOLATE_LANCZOS4, 8, 8, 1.0e-3); // 66 tests
    check(32, 32, PS_TYPE_F64, 32, PS_INTERPOLATE_LANCZOS4, 8, 8, 1.0e-3); // 66 tests

    return exit_status();
}

