#include <algorithm>
#include <cmath>
#include <cstdint>
#include <tuple>
#include <iostream>
#include <utility>

std::tuple<float, float, float> rgb_to_xyz(const uint8_t p_red, const uint8_t p_green, const uint8_t p_blue) {
    auto normalize_uint8_t = [](const uint8_t p_value) -> float {
        return std::min(static_cast<float>(p_value) / 255.0F, 1.0F);
    };

    auto gamma_correction = [](const float p_value) -> float {
        float gamma_corrected_value = 0.0F;

        if (p_value <= 0.0405F) {
            gamma_corrected_value = p_value / 12.92F;
        } else {
            gamma_corrected_value = std::pow((p_value + 0.055F) / 1.055F, 2.4F);
        }

        return gamma_corrected_value;
    };

    float red = gamma_correction(normalize_uint8_t(p_red));
    float green = gamma_correction(normalize_uint8_t(p_green));
    float blue = gamma_correction(normalize_uint8_t(p_blue));

    float x = 0.4124564F * red + 0.3575761F * green + 0.1804375F * blue;
    float y = 0.2126729F * red + 0.7151522F * green + 0.0721750F * blue;
    float z = 0.0193339F * red + 0.1191920F * green + 0.9503041F * blue;

    return {x * 100.0F, y * 100.0F, z * 100.0F};
}

std::tuple<float, float, float> xyz_to_rgb_float(float p_x, float p_y, float p_z) {
    auto reverse_gamma_correction = [](const float p_value) -> float {
        float reverse_gamma_corrected_value = 0.0F;

        if (p_value <= 0.0031308F) {
            reverse_gamma_corrected_value = p_value * 12.92F;
        } else {
            reverse_gamma_corrected_value = std::pow(p_value, 1.0F / 2.4F) * 1.055F - 0.055F;
        }

        return reverse_gamma_corrected_value;
    };

    auto clamp_if_within_error_margin = [](const float p_value) -> float {
        float value = p_value;
        constexpr float MARGIN = 0.001F;

        if ((p_value >= -MARGIN) and (p_value < 0.0F)) {
            value = 0.0F;
        } else if ((p_value <= 1.0F + MARGIN) and (p_value > 1.0F)) {
            value = 1.0F;
        }

        return value;
    };

    float red = 3.2404542F * p_x / 100.0F - 1.5371385F * p_y / 100.0F - 0.4985314F * p_z / 100.0F;
    float green = -0.9692660F * p_x / 100.0F + 1.8760108F * p_y / 100.0F + 0.0415560F * p_z / 100.0F;
    float blue = 0.0556434F * p_x / 100.0F - 0.2040259F * p_y / 100.0F + 1.0572252F * p_z / 100.0F;

    red = clamp_if_within_error_margin(reverse_gamma_correction(red));
    green = clamp_if_within_error_margin(reverse_gamma_correction(green));
    blue = clamp_if_within_error_margin(reverse_gamma_correction(blue));

    return {red, green, blue};
}

std::tuple<uint8_t, uint8_t, uint8_t> xyz_to_rgb(float p_x, float p_y, float p_z) {
    auto uint8_t_scale = [](const float p_value) -> uint8_t {
        return static_cast<uint8_t>(std::clamp(p_value, 0.0F, 1.0F) * 255.0F + 0.5F);
    };

    auto [red, green, blue] = xyz_to_rgb_float(p_x, p_y, p_z);

    uint8_t red_uint8_t = uint8_t_scale(red);
    uint8_t green_uint8_t = uint8_t_scale(green);
    uint8_t blue_uint8_t = uint8_t_scale(blue);

    return {red_uint8_t, green_uint8_t, blue_uint8_t};
}

std::tuple<float, float> xyz_to_uv(const float p_x, const float p_y, const float p_z) {
    float denominator = (p_x + 15.0F * p_y + 3.0F * p_z);

    if (denominator == 0.0F) {
        denominator = 1.0F;
    }

    float u = (4.0F * p_x) / denominator;
    float v = (9.0F * p_y) / denominator;

    return {u, v};
}

std::tuple<float, float> uv_to_xy(const float p_u, const float p_v) {
    float denominator = 6.0F * p_u - 16.0F * p_v + 12.0F;

    if (denominator == 0.0F) {
        denominator = 1.0F;
    }

    float x = 9.0F * p_u / denominator;
    float y = 4.0F * p_v / denominator;

    return {x, y};
}

std::tuple<float, float, float> xyz_d65_reference_white_point() {
    // Standard D64 reference white point values
    constexpr float X_N = 95.047F;
    constexpr float Y_N = 100.0F;
    constexpr float Z_N = 108.833F;

    return {X_N, Y_N, Z_N};
}

std::tuple<float, float, float> xyz_to_cieluv(const float p_x, const float p_y, const float p_z) {
    auto [X_N, Y_N, Z_N] = xyz_d65_reference_white_point();
    auto [U_N, V_N] = xyz_to_uv(X_N, Y_N, Z_N);

    auto [u, v] = xyz_to_uv(p_x, p_y, p_z);

    float l_star = 0.0F;

    if (p_y == 0.0F) {
        l_star = 0.0F;
    } else if ((p_y / Y_N) > 0.0088565F) {
        l_star = 116.0F * std::cbrt(p_y / Y_N) - 16.0F;
    } else {
        l_star = p_y / Y_N * 903.3F;
    }

    float u_star = 13.0F * l_star * (u - U_N);
    float v_star = 13.0F * l_star * (v - V_N);

    return {l_star, u_star, v_star};
}

std::tuple<float, float, float> cieluv_to_xyz(float p_l_star, float p_u_star, float p_v_star) {
    auto [X_N, Y_N, Z_N] = xyz_d65_reference_white_point();
    auto [U_N, V_N] = xyz_to_uv(X_N, Y_N, Z_N);

    float denominator = 13.0F * p_l_star;

    if (denominator == 0.0F) {
        denominator = 13.0F;
    }

    float u = p_u_star / denominator + U_N;
    float v = p_v_star / denominator + V_N;

    // auto [x, y] = uv_to_xy(u, v);

    float y = 0.0F;
    if (p_l_star <= 8.0F) {
        y = Y_N * p_l_star * 0.0011071F;
    } else {
        y = Y_N * std::pow((p_l_star + 16.F) / 116.0F, 3.0F);
    }

    if (v == 0.0F) {
        denominator = 4.0F;
    } else {
        denominator = 4.0F * v;
    }
    float x = y * 9.0F * u / denominator;
    float z = y * (12.0F - 3.0F * u - 20.0F * v) / denominator;

    return {x, y, z};
}

std::tuple<float, float, float> cieluv_to_cielchuv(float p_l_star, float p_u_star, float p_v_star) {
    float c_star = std::sqrt(p_u_star * p_u_star + p_v_star * p_v_star);
    float h_star = std::atan2(p_v_star, p_u_star);

    if (h_star < 0.0F) {
        h_star += 2.0F * M_PI;
    }

    return {p_l_star, c_star, h_star};
}

std::tuple<float, float, float> cielchuv_to_cieluv(float p_l_star, float p_c_star, float p_h_star) {
    float u_star = p_c_star * std::cos(p_h_star);
    float v_star = p_c_star * std::sin(p_h_star);

    return {p_l_star, u_star, v_star};
}

bool is_xyz_in_rgb_range(float p_x, float p_y, float p_z) {
    auto is_in_range_inclusive = [](const float p_value, const float p_min, const float p_max) -> bool {
        return (p_value >= p_min and p_value <= p_max);
    };

    auto [red, green, blue] = xyz_to_rgb_float(p_x, p_y, p_z);

    bool is_red_in_range = is_in_range_inclusive(red, 0.0F, 1.0F);
    bool is_green_in_range = is_in_range_inclusive(green, 0.0F, 1.0F);
    bool is_blue_in_range = is_in_range_inclusive(blue, 0.0F, 1.0F);

    return is_red_in_range and is_green_in_range and is_blue_in_range;
}

bool is_cieluv_in_rgb_range(float p_l_star, float p_u_star, float p_v_star) {
    auto [x, y, z] = cieluv_to_xyz(p_l_star, p_u_star, p_v_star);

    return is_xyz_in_rgb_range(x, y, z);
}

bool is_cielchuv_in_rgb_range(float p_l_star, float p_c_star, float p_h_star) {
    auto [l, u, v] = cielchuv_to_cieluv(p_l_star, p_c_star, p_h_star);

    return is_cieluv_in_rgb_range(l, u, v);
}

std::tuple<float, float, float> map_cielchuv_to_visible_color(float p_l_star, float p_c_star, float p_h_star) {
    float c_star = p_c_star;
    uint32_t count = 0U;

    while (not is_cielchuv_in_rgb_range(p_l_star, c_star, p_h_star)) {
        c_star = c_star * 0.99F;
        count = count + 1U;

        if (count > 1000U) {
            c_star = 0.0F;
            std::cout << "breaky breaky!" << std::endl;
            break;
        }
    }

    return {p_l_star, c_star, p_h_star};
}

std::tuple<float, float, float> rgb_to_cielchuv(uint8_t p_red, uint8_t p_green, uint8_t p_blue) {
    auto [x, y, z] = rgb_to_xyz(p_red, p_green, p_blue);
    auto [l, u, v] = xyz_to_cieluv(x, y, z);
    auto [l_star, c_star, h_star] = cieluv_to_cielchuv(l, u, v);

    return {l_star, c_star, h_star};
}

std::tuple<uint8_t, uint8_t, uint8_t> cielchuv_to_rgb(float p_l_star, float p_c_star, float p_h_star) {
    auto [l, u, v] = cielchuv_to_cieluv(p_l_star, p_c_star, p_h_star);
    auto [x, y, z] = cieluv_to_xyz(l, u, v);
    auto [red, green, blue] = xyz_to_rgb(x, y, z);

    return {red, green, blue};
}

int main(void) {
    uint8_t r = 0;
    uint8_t g = 0;
    uint8_t b = 0;

    for (r = 0; r < 192; r = r + 64) {
        for (g = 0; g < 192; g = g + 64) {
            for (b = 0; b < 192; b = b + 64) {
                auto [x, y, z] = rgb_to_xyz(r, g, b);
                auto [l, u, v] = xyz_to_cieluv(x, y, z);
                auto [x2, y2, z2] = cieluv_to_xyz(l, u, v);
                auto [r2, g2, b2] = xyz_to_rgb(x2, y2, z2);
     
                std::cout << "RGB: [" + std::to_string(r) + ", " + std::to_string(g) + ", " + std::to_string(b) + "]" + 
                        " --> XYZ: [" + std::to_string(x) + ", " + std::to_string(y) + ", " + std::to_string(z) + "]" + 
                        " --> LUV: [" + std::to_string(l) + ", " + std::to_string(u) + ", " + std::to_string(v) + "]" + 
                        " --> XYZ: [" + std::to_string(x2) + ", " + std::to_string(y2) + ", " + std::to_string(z2) + "]" +
                        " --> RGB: [" + std::to_string(r2) + ", " + std::to_string(g2) + ", " + std::to_string(b2) + "]" +
                        " --> E: [" + std::to_string(r - r2) + ", " + std::to_string(g - g2) + ", " + std::to_string(b - b2) + "]" << std::endl;
            }
        }
    }

    std::cout << std::endl;

    for (r = 0; r < 192; r = r + 64) {
        for (g = 0; g < 192; g = g + 64) {
            for (b = 0; b < 192; b = b + 64) {
                auto [x, y, z] = rgb_to_xyz(r, g, b);
                auto [l, u, v] = xyz_to_cieluv(x, y, z);
                auto [l_, c, h] = cieluv_to_cielchuv(l, u, v);

                auto [l2, c2, h2] = map_cielchuv_to_visible_color(l / 2.0F, c, h);

                auto [l2_, u2, v2] = cielchuv_to_cieluv(l2, c2, h2);
                auto [x2, y2, z2] = cieluv_to_xyz(l2, u2, v2);
                auto [r2_f, g2_f, b2_f] = xyz_to_rgb_float(x2, y2, z2);
                auto [r2, g2, b2] = xyz_to_rgb(x2, y2, z2);

                std::cout << "RGB: [" + std::to_string(r) + ", " + std::to_string(g) + ", " + std::to_string(b) + "]" + 
                        " --> XYZ: [" + std::to_string(x) + ", " + std::to_string(y) + ", " + std::to_string(z) + "]" +
                        " --> LUV: [" + std::to_string(l) + ", " + std::to_string(u) + ", " + std::to_string(v) + "]" + 
                        " --> LCH: [" + std::to_string(l) + ", " + std::to_string(c) + ", " + std::to_string(h) + "]" + 
                        " --> / 2: [" + std::to_string(l2) + ", " + std::to_string(c2) + ", " + std::to_string(h2) + "]" + 
                        "; visible: " + std::to_string(is_cielchuv_in_rgb_range(l / 2.0F, c, h)) + 
                        " / " + std::to_string(is_cielchuv_in_rgb_range(l2, c2, h2)) + 
                        " --> LUV: [" + std::to_string(l2) + ", " + std::to_string(u2) + ", " + std::to_string(v2) + "]" + 
                        " --> XYZ: [" + std::to_string(x2) + ", " + std::to_string(y2) + ", " + std::to_string(z2) + "]" + 
                        " --> RGB: [" + std::to_string(r2_f) + ", " + std::to_string(g2_f) + ", " + std::to_string(b2_f) + "]" +
                        " --> RGB: [" + std::to_string(r2) + ", " + std::to_string(g2) + ", " + std::to_string(b2) + "]" << std::endl;
            }
        }
    }

    // std::cout << std::endl;

    // float x = 0.925094F;
    // float y = 0.370037F;
    // float z = 4.872159F;
    // auto [red, green, blue] = xyz_to_rgb_float(x, y, z);
    // std::cout << std::to_string(red) + " " + std::to_string(green) + " " + std::to_string(blue) + " " + std::to_string(is_xyz_in_rgb_range(x, y, z)) << std::endl;

    // float x2 = 0.740333F;
    // float y2 = 0.185027F;
    // float z2 = 0.031254F;
    // auto [red2, green2, blue2] = xyz_to_rgb_float(x2, y2, z2);
    // std::cout << std::to_string(red2) + " " + std::to_string(green2) + " " + std::to_string(blue2) + " " + std::to_string(is_xyz_in_rgb_range(x2, y2, z2)) << std::endl;

    // // float L = 3.342548F;
    // float L = 1.671274F;
    // float C = 13.526361F;
    // float H = 4.640314F;
    // 
    // auto [L_, U, V] = cielchuv_to_cieluv(L, C, H);
    // auto [X, Y, Z] = cieluv_to_xyz(L_, U, V);
    // auto [R, G, B] = xyz_to_rgb_float(X, Y, Z);
    // std::cout << "RGB: " << std::to_string(R) + " " + std::to_string(G) + " " + std::to_string(B) + 
    //            ", XYZ: " + std::to_string(X) + " " + std::to_string(Y) + " " + std::to_string(Z) + 
    //            ", LUV: " + std::to_string(L_) + " " + std::to_string(U) + " " + std::to_string(V) + 
    //            " " + std::to_string(is_xyz_in_rgb_range(X, Y, Z)) << std::endl;
}