// Copyright 2023 Google LLC

// SPDX-License-Identifier: BSD-2-Clause

#include "avif/internal.h"

#include <math.h>

static const double epsilon = 1e-12;

static avifBool avifXyToXYZ(const float xy[2], double XYZ[3])

{

    if (fabsf(xy[1]) < epsilon) {

        return AVIF_FALSE;

    }

    const double factor = 1.0 / xy[1];

    XYZ[0] = xy[0] * factor;

    XYZ[1] = 1;

    XYZ[2] = (1 - xy[0] - xy[1]) * factor;

    return AVIF_TRUE;

}

// Computes I = M^-1. Returns false if M seems to be singular.

static avifBool avifMatInv(double M[3][3], double I[3][3])

{

    double det = M[0][0] * (M[1][1] * M[2][2] - M[2][1] * M[1][2]) - M[0][1] * (M[1][0] * M[2][2] - M[1][2] * M[2][0]) +

                 M[0][2] * (M[1][0] * M[2][1] - M[1][1] * M[2][0]);

    if (fabs(det) < epsilon) {

        return AVIF_FALSE;

    }

    det = 1.0 / det;

    I[0][0] = (M[1][1] * M[2][2] - M[2][1] * M[1][2]) * det;

    I[0][1] = (M[0][2] * M[2][1] - M[0][1] * M[2][2]) * det;

    I[0][2] = (M[0][1] * M[1][2] - M[0][2] * M[1][1]) * det;

    I[1][0] = (M[1][2] * M[2][0] - M[1][0] * M[2][2]) * det;

    I[1][1] = (M[0][0] * M[2][2] - M[0][2] * M[2][0]) * det;

    I[1][2] = (M[1][0] * M[0][2] - M[0][0] * M[1][2]) * det;

    I[2][0] = (M[1][0] * M[2][1] - M[2][0] * M[1][1]) * det;

    I[2][1] = (M[2][0] * M[0][1] - M[0][0] * M[2][1]) * det;

    I[2][2] = (M[0][0] * M[1][1] - M[1][0] * M[0][1]) * det;

    return AVIF_TRUE;

}

// Computes C = A*B

static void avifMatMul(double A[3][3], double B[3][3], double C[3][3])

{

    C[0][0] = A[0][0] * B[0][0] + A[0][1] * B[1][0] + A[0][2] * B[2][0];

    C[0][1] = A[0][0] * B[0][1] + A[0][1] * B[1][1] + A[0][2] * B[2][1];

    C[0][2] = A[0][0] * B[0][2] + A[0][1] * B[1][2] + A[0][2] * B[2][2];

    C[1][0] = A[1][0] * B[0][0] + A[1][1] * B[1][0] + A[1][2] * B[2][0];

    C[1][1] = A[1][0] * B[0][1] + A[1][1] * B[1][1] + A[1][2] * B[2][1];

    C[1][2] = A[1][0] * B[0][2] + A[1][1] * B[1][2] + A[1][2] * B[2][2];

    C[2][0] = A[2][0] * B[0][0] + A[2][1] * B[1][0] + A[2][2] * B[2][0];

    C[2][1] = A[2][0] * B[0][1] + A[2][1] * B[1][1] + A[2][2] * B[2][1];

    C[2][2] = A[2][0] * B[0][2] + A[2][1] * B[1][2] + A[2][2] * B[2][2];

}

// Set M to have values of d on the leading diagonal, and zero elsewhere.

static void avifMatDiag(const double d[3], double M[3][3])

{

    M[0][0] = d[0];

    M[0][1] = 0;

    M[0][2] = 0;

    M[1][0] = 0;

    M[1][1] = d[1];

    M[1][2] = 0;

    M[2][0] = 0;

    M[2][1] = 0;

    M[2][2] = d[2];

}

// Computes y = M.x

static void avifVecMul(double M[3][3], const double x[3], double y[3])

{

    y[0] = M[0][0] * x[0] + M[0][1] * x[1] + M[0][2] * x[2];

    y[1] = M[1][0] * x[0] + M[1][1] * x[1] + M[1][2] * x[2];

    y[2] = M[2][0] * x[0] + M[2][1] * x[1] + M[2][2] * x[2];

}

// Bradford chromatic adaptation matrix

// from https://www.researchgate.net/publication/253799640_A_uniform_colour_space_based_upon_CIECAM97s

static double avifBradford[3][3] = {

    { 0.8951, 0.2664, -0.1614 },

    { -0.7502, 1.7135, 0.0367 },

    { 0.0389, -0.0685, 1.0296 },

};

// LMS values for D50 whitepoint

static const double avifLmsD50[3] = { 0.996284, 1.02043, 0.818644 };

avifBool avifColorPrimariesComputeRGBToXYZD50Matrix(avifColorPrimaries colorPrimaries, double coeffs[3][3])

{

    float primaries[8];

    avifColorPrimariesGetValues(colorPrimaries, primaries);

    double whitePointXYZ[3];

    AVIF_CHECK(avifXyToXYZ(&primaries[6], whitePointXYZ));

    double rgbPrimaries[3][3] = {

        { primaries[0], primaries[2], primaries[4] },

        { primaries[1], primaries[3], primaries[5] },

        { 1.0 - primaries[0] - primaries[1], 1.0 - primaries[2] - primaries[3], 1.0 - primaries[4] - primaries[5] }

    };

    double rgbPrimariesInv[3][3];

    AVIF_CHECK(avifMatInv(rgbPrimaries, rgbPrimariesInv));

    double rgbCoefficients[3];

    avifVecMul(rgbPrimariesInv, whitePointXYZ, rgbCoefficients);

    double rgbCoefficientsMat[3][3];

    avifMatDiag(rgbCoefficients, rgbCoefficientsMat);

    double rgbXYZ[3][3];

    avifMatMul(rgbPrimaries, rgbCoefficientsMat, rgbXYZ);

    // ICC stores primaries XYZ under PCS.

    // Adapt using linear bradford transform

    // from https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781119021780.app3

    double lms[3];

    avifVecMul(avifBradford, whitePointXYZ, lms);

    for (int i = 0; i < 3; ++i) {

        if (fabs(lms[i]) < epsilon) {

            return AVIF_FALSE;

        }

        lms[i] = avifLmsD50[i] / lms[i];

    }

    double adaptation[3][3];

    avifMatDiag(lms, adaptation);

    double tmp[3][3];

    avifMatMul(adaptation, avifBradford, tmp);

    double bradfordInv[3][3];

    if (!avifMatInv(avifBradford, bradfordInv)) {

        return AVIF_FALSE;

    }

    avifMatMul(bradfordInv, tmp, adaptation);

    avifMatMul(adaptation, rgbXYZ, coeffs);

    return AVIF_TRUE;

}

avifBool avifColorPrimariesComputeXYZD50ToRGBMatrix(avifColorPrimaries colorPrimaries, double coeffs[3][3])

{

    double rgbToXyz[3][3];

    AVIF_CHECK(avifColorPrimariesComputeRGBToXYZD50Matrix(colorPrimaries, rgbToXyz));

    AVIF_CHECK(avifMatInv(rgbToXyz, coeffs));

    return AVIF_TRUE;

}

avifBool avifColorPrimariesComputeRGBToRGBMatrix(avifColorPrimaries srcColorPrimaries,

                                                 avifColorPrimaries dstColorPrimaries,

                                                 double coeffs[3][3])

{

    // Note: no special casing for srcColorPrimaries == dstColorPrimaries to allow

    // testing that the computation actually produces the identity matrix.

    double srcRGBToXYZ[3][3];

    AVIF_CHECK(avifColorPrimariesComputeRGBToXYZD50Matrix(srcColorPrimaries, srcRGBToXYZ));

    double xyzToDstRGB[3][3];

    AVIF_CHECK(avifColorPrimariesComputeXYZD50ToRGBMatrix(dstColorPrimaries, xyzToDstRGB));

    // coeffs = xyzToDstRGB * srcRGBToXYZ

    // i.e. srcRGB -> XYZ -> dstRGB

    avifMatMul(xyzToDstRGB, srcRGBToXYZ, coeffs);

    return AVIF_TRUE;

}

// Converts a linear RGBA pixel to a different color space. This function actually works for gamma encoded

// RGB as well but linear gives better results. Also, for gamma encoded values, it would be

// better to clamp the output to [0, 1]. Linear values don't need clamping because values

// > 1.0 are valid for HDR transfer curves, and the gamma compression function will do the

// clamping as necessary.

void avifLinearRGBConvertColorSpace(float rgb[4], double coeffs[3][3])

{

    const double rgbDouble[3] = { rgb[0], rgb[1], rgb[2] };

    double converted[3];

    avifVecMul(coeffs, rgbDouble, converted);

    rgb[0] = (float)converted[0];

    rgb[1] = (float)converted[1];

    rgb[2] = (float)converted[2];

}