#include <math.h>
#include "colourspace.h"

static double intensity_to_voltage_srgb( double val )
{
    /* Handle invalid values before doing a gamma transform. */
    if( val < 0.0 ) return 0.0;
    if( val > 1.0 ) return 1.0;

    /* sRGB uses an exponent of 2.4 in its nonstandard curve, but it
     * still advertises a display gamma of 2.2.
     */
    if( val <= 0.0031308 ) {
        return 12.92 * val;
    }
    return ( ( 1.055 * pow( val, 1.0 / 2.4 ) ) - 0.055 );
}

static double voltage_to_intensity_srgb( double val )
{
    /* Handle invalid values before doing a gamma transform. */
    if( val < 0.0 ) return 0.0;
    if( val > 1.0 ) return 1.0;

    if( val <= 0.04045 ) {
        return val / 12.92;
    }
    return pow( ( val + 0.055 ) / 1.055, 2.4 );
}


/**
 * Rec-709 specifies a gamma of 2.2 and the following curve with linear
 * segment.
 */

static double intensity_to_voltage_rec709( double val )
{
    /* Handle invalid values before doing a gamma transform. */
    if( val < 0.0 ) return 0.0;
    if( val > 1.0 ) return 1.0;

    if( val <= 0.018 ) {
        return 4.5 * val;
    }
    return ( ( 1.099 * pow( val, 1.0 / 2.2 ) ) - 0.099 );
}

static double voltage_to_intensity_rec709( double val )
{
    /* Handle invalid values before doing a gamma transform. */
    if( val < 0.0 ) return 0.0;
    if( val > 1.0 ) return 1.0;

    if( val <= 0.081 ) {
        return val / 4.5;
    }
    return pow( ( val + 0.099 ) / 1.099, 2.2 );
}

void linear_to_nonlinear_srgb( double   r, double   g, double   b,
                               double *rp, double *gp, double *bp )
{
    *rp = intensity_to_voltage_srgb( r );
    *gp = intensity_to_voltage_srgb( g );
    *bp = intensity_to_voltage_srgb( b );
}

void nonlinear_to_linear_srgb( double  rp, double  gp, double  bp,
                               double  *r, double  *g, double  *b )
{
    *r = voltage_to_intensity_srgb( rp );
    *g = voltage_to_intensity_srgb( gp );
    *b = voltage_to_intensity_srgb( bp );
}

void linear_to_nonlinear_rec709( double   r, double   g, double   b,
                                 double *rp, double *gp, double *bp )
{
    *rp = intensity_to_voltage_rec709( r );
    *gp = intensity_to_voltage_rec709( g );
    *bp = intensity_to_voltage_rec709( b );
}

void nonlinear_to_linear_rec709( double  rp, double  gp, double  bp,
                                 double  *r, double  *g, double  *b )
{
    *r = voltage_to_intensity_rec709( rp );
    *g = voltage_to_intensity_rec709( gp );
    *b = voltage_to_intensity_rec709( bp );
}

void xyz_to_rgb_rec709( double  x, double  y, double  z,
                        double *r, double *g, double *b )
{
    *r = (  3.240479 * x ) + ( -1.537150 * y ) + ( -0.498535 * z );
    *g = ( -0.969256 * x ) + (  1.875992 * y ) + (  0.041556 * z );
    *b = (  0.055648 * x ) + ( -0.204043 * y ) + (  1.057311 * z );
}

void rgb_to_xyz_rec709( double  r, double  g, double  b,
                        double *x, double *y, double *z )
{
    *x = ( 0.412453 * r ) + ( 0.357580 * g ) + ( 0.180423 * b );
    *y = ( 0.212671 * r ) + ( 0.715160 * g ) + ( 0.072169 * b );
    *z = ( 0.019334 * r ) + ( 0.119193 * g ) + ( 0.950227 * b );
}

void xyz_to_srgb( double  x, double  y, double  z,
                  double *r, double *g, double *b )
{
    *r = (  3.2406 * x ) + ( -1.5372 * y ) + ( -0.4986 * z );
    *g = ( -0.9689 * x ) + (  1.8758 * y ) + (  0.0415 * z );
    *b = (  0.0557 * x ) + ( -0.2040 * y ) + (  1.0570 * z );
}

void srgb_to_xyz( double  r, double  g, double  b,
                  double *x, double *y, double *z )
{
    *x = ( 0.4124 * r ) + ( 0.3576 * g ) + ( 0.1805 * b );
    *y = ( 0.2126 * r ) + ( 0.7152 * g ) + ( 0.0722 * b );
    *z = ( 0.0193 * r ) + ( 0.1192 * g ) + ( 0.9505 * b );
}

void rgb_to_fakelab_rec709( double r, double g, double b,
                            double *lstar, double *astar, double *bstar )
{
    double x, y, z;
    rgb_to_xyz_rec709( r, g, b, &x, &y, &z );

    *lstar = pow( y, 1.0 / 3.0 );
    *astar = pow( x, 1.0 / 3.0 );
    *bstar = pow( z, 1.0 / 3.0 );
}

void fakelab_to_rgb_rec709( double lstar, double astar, double bstar,
                            double *r, double *g, double *b )
{
    double x, y, z;

    y = lstar * lstar * lstar;
    x = astar * astar * astar;
    z = bstar * bstar * bstar;

    xyz_to_rgb_rec709( x, y, z, r, g, b );
}

void xyy_to_xyz( double x, double y, double Y,
                 double *rx, double *ry, double *rz )
{
/*
    *ry = Y;
    *rx = ( x / y ) * Y;
    *rz = ( ( 1.0 - x - y ) / y ) * Y;
*/
    *ry = Y;
    *rx = ( x / y ) * Y;
    *rz = ( *rx / x ) - *rx - Y;
}

void uvy_to_xyz( double up, double vp, double Y,
                 double *rx, double *ry, double *rz )
{
    *ry = Y;
    *rx = ( 9.0 * Y * up ) / ( 4.0 * vp );
    *rz = ( ( ( 4.0 * *rx ) / up ) - *rx - ( 15.0 * Y ) ) / 3.0;
}

void closest_in_gamut_srgb( double r, double g, double b,
                            double wr, double wg, double wb,
                            double *ri, double *gi, double *bi )
{
    *ri = r;
    *gi = g;
    *bi = b;

    if( *ri >= 0.0 && *gi >= 0.0 && *bi >= 0.0 ) return;

    for(;;) {
        double vr = wr - *ri;
        double vg = wg - *gi;
        double vb = wb - *bi;
        double scale;

        if( *ri >= 0.0 && *gi >= 0.0 && *bi >= 0.0 ) return;

        if( *ri < *gi && *ri < *bi ) {
            scale = ( -(*ri) ) / vr;
        } else if( *gi < *ri && *gi < *bi ) {
            scale = ( -(*gi) ) / vg;
        } else {
            scale = ( -(*bi) ) / vb;
        }
        *ri = *ri + ( scale * vr );
        *gi = *gi + ( scale * vg );
        *bi = *bi + ( scale * vb );
    }
}