#include "math_lib.h"

#include <math.h>


/**
 * @brief The function f(x) to define
 *
 * @param x The x param for the function
 * @return double return the result of f(x)
 */
double my_function_x(double x)
{
    return -((2 * x - 1) / (exp(pow(x, 2) - x + 2)));
}

/**
 * @brief Intergration of my_function_x
 *
 * @param a the a
 * @param b the b
 * @param N the N
 * @param f the function to use
 *
 * @return double the result of integrations
 */
double interg(double a, double b, double N, double (*f)(double))
{
    double delta_x = ((b - a) / N);
    double res = 0.0;

    for (int32_t i = 0; i < N; i += 1)
    {
        res += (*f)(a + (i * delta_x)) * delta_x;
    }
    return res;
}

/**
 * @brief Compute E(n) with function f(x)
 *
 * @param N The N number
 * @return double return the result of E(N)
 */
double E_n(double N)
{
    return fabs((VALUE_I - interg(VALUE_A, VALUE_B, N, my_function_x)) / VALUE_I);
}

/**
 * @brief Compute the singal s(x)
 *
 * @param x The param x for the signal
 * @return double return the result of s(x)
 */
double s_x(double x, double w_1, double w_2)
{
    return sin(2 * M_PI * w_1 * x) + sin(2 * M_PI * w_2 * x);
}

/**
 * @brief
 *
 * @param x
 * @param psi
 * @return double
 */
double h_x(double x, double psi)
{
    return (1 / sqrt(2 * M_PI * psi)) * exp(-pow(x, 2) / (2 * psi));
}

matrix* convolve(matrix img, matrix kernel,matrix* new,const double factor)
{
    // matrix new = NULL;
    // matrix_alloc(&new,img->col,img->row);
    // matrix_init(&new);
    matrix_clone(new,img);
    matrix_init(new,img.m,img.n);
    uint32_t kernel_size = kernel.col * kernel.row;
    uint32_t pixels = 0;
    int32_t k_size = (sqrt(kernel_size) / 2);
    for (int32_t i = 0; i < img.col; i += 1)
    {
        for (int32_t j = 0; i < img.row; j += 1)
        {
            for (int32_t k = 0; k < kernel.col; k += 1)
            {
                for (int32_t l = 0; l < kernel.row; l += 1)
                {
                    if ((i - k) < 0 || (i + k) >= img.col)
                    {
                        continue;
                    }
                    if ((j - l) < 0 || (j - l) >= img.row)
                    {
                        continue;
                    }
                    pixels += img.data[i + k][j + l] * kernel.data[k_size + k][k_size + l];
                }
                new->data[i][j] = (pixels / kernel_size) * factor;
                pixels = 0;
            }
        }
    }
    return new;
}


matrix* convolve_matrix(matrix* conv,matrix* orig, const matrix *const kernel, const double factor)
{
    matrix tmp;
    //matrix_alloc(&tmp, kernel->col, kernel->row);

    tmp = matrix_creat(kernel->col,kernel->row);
    //tmp = matrix_create(kernel->col,kernel->row,0);
    //matrix_alloc(conv,orig->col,orig->row);

    // Browse PGM struct orig and create sub matrix but with default value
    // int32_t because (res < 0)
    int32_t res = 0;
    for (int i = 0; i < orig->col; i += 1)
    {
        for (int j = 0; j < orig->row; j += 1)
        {
            res = 0;
            matrix_reset(&tmp);
            
            for (int k = 0, base_i = i - (int)kernel->col / 2; k < kernel->col; k += 1, base_i += 1)
            {
                for (int l = 0, base_j = j - (int)kernel->row / 2; l < kernel->row; l += 1, base_j += 1)
                {
                    if (base_i < 0 || base_j < 0 || base_i > orig->col - 1 || base_j > orig->row - 1)
                    {
                        tmp.data[k][l] = 0; // Default value for out of bounds
                    }
                    else
                    {
                        tmp.data[k][l] = orig->data[base_i][base_j];
                    }
                }
            }

            // uint32_t tm, ker;
            // Use matrix tmp with sub matrix
            for (int k = 0; k < tmp.col; k += 1)
            {
                for (int l = 0; l < tmp.row; l += 1)
                {
                    // tm = tmp.data[k][l];
                    // ker = kernel->data[k][l];
                    res +=  tmp.data[k][l] * kernel->data[k][l];
                }
            }

            //printf("Res(before) : %d / %lf\n",res,factor);
            res *= factor; // Divide by factor of kernel
            //printf("Res(after) : %d\n",res);

            if(res > 65535)
            {
                res = 65535;
            }
            if (res < 0)
            {
                res = 0;
            }

            conv->data[i][j] = res;
        }
    }
    //free(&tmp.data);
    //matrix_destroy(&tmp);       // PROBLEME DE FREEE
    matrix_destroy(&tmp);

    return NULL;
}

double function_for_trapez(double a, uint32_t N, double (*f)(double), double (*O)(double))
{
    double delta_x = 0.0;
    double res = 0.0;
    for (uint32_t i = 0; i < N; i += 1)
    {
        res += ((*f)(a + i * delta_x) + (*f)(a + (i + 1) * delta_x) / 2) + O(delta_x * delta_x);
    }
    return res;
}

// matrix *convolv_2d(matrix *m1, matrix *kernel)
// {
//     matrix *mat = NULL;
//     matrix_alloc(mat, m1->col, m1->row);
//     matrix_init(mat);
//     matrix_print(*mat);
// }