/**
 * Copyright (c) 2003 Billy Biggs <vektor@dumbterm.net>
 *
 * Permission is hereby granted, free of charge, to any person obtaining
 * a copy of this software and associated documentation files (the
 * "Software"), to deal in the Software without restriction, including
 * without limitation the rights to use, copy, modify, merge, publish,
 * distribute, sublicense, and/or sell copies of the Software, and to
 * permit persons to whom the Software is furnished to do so, subject to
 * the following conditions:
 *
 * The above copyright notice and this permission notice shall be
 * included in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
 * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
 * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
 * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 *
 * Last updated Mon Jun 16 13:12:17 ADT 2003 by vektor
 *  - Initial release.
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include "huerotate.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 );
}

static 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 );
}

static 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 );
}

static 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 );
}

static 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 );
}

static void xyz_to_cielab( double x, double y, double z,
                           double *l, double *a, double *b )
{
    double ypow;

    ypow = pow( y, 1.0 / 3.0 );
    *l = ( 116.0 * ypow ) - 16.0;
    *a = 500.0 * ( pow( x, 1.0 / 3.0 ) - ypow );
    *b = 200.0 * ( ypow - pow( z, 1.0 / 3.0 ) );
}

static void cielab_to_xyz( double l, double a, double b,
                           double *x, double *y, double *z )
{
    double p;

    p = ( l + 16.0 ) / 116.0;
    *y = pow( p, 3.0 );
    *x = pow( p + ( a / 500.0 ), 3.0 );
    *z = pow( p - ( b / 200.0 ), 3.0 );
}

#ifndef M_PI
#define M_PI 3.14159265358979323846264338327
#endif

static void rotate_hue( double rp, double gp, double bp, double angle,
                        double *rrp, double *rgp, double *rbp )
{
    double r, g, b, X, Y, Z, Lstar, astar, bstar, hueab, chroma;

    nonlinear_to_linear_srgb( rp, gp, bp, &r, &g, &b );
    srgb_to_xyz( r, g, b, &X, &Y, &Z );
    xyz_to_cielab( X, Y, Z, &Lstar, &astar, &bstar );

    hueab = atan2( bstar, astar );
    chroma = sqrt( (astar * astar) + (bstar * bstar) );

    hueab += angle;
    if( hueab > M_PI ) hueab -= (2.0 * M_PI);

    astar = chroma * cos( hueab );
    bstar = chroma * sin( hueab );

    cielab_to_xyz( Lstar, astar, bstar, &X, &Y, &Z );
    xyz_to_srgb( X, Y, Z, &r, &g, &b );
    linear_to_nonlinear_srgb( r, g, b, rrp, rgp, rbp );
}

void hue_rotate_scanline_rgb24( uint8_t *dst, uint8_t *src, int width, double angle )
{
    int x;

    for( x = 0; x < width; x++ ) {
        double rp, gp, bp, rrp, rgp, rbp;

        rp = (double) src[ (x * 3) + 0 ] / 255.0;
        gp = (double) src[ (x * 3) + 1 ] / 255.0;
        bp = (double) src[ (x * 3) + 2 ] / 255.0;

        rotate_hue( rp, gp, bp, angle, &rrp, &rgp, &rbp );

        dst[ (x * 3) + 0 ] = (int) ( ( rrp * 255.0 ) + 0.5 );
        dst[ (x * 3) + 1 ] = (int) ( ( rgp * 255.0 ) + 0.5 );
        dst[ (x * 3) + 2 ] = (int) ( ( rbp * 255.0 ) + 0.5 );
    }
}

void hue_rotate_scanline_rgba32( uint8_t *dst, uint8_t *src, int width, double angle )
{
    int x;

    for( x = 0; x < width; x++ ) {
        double rp, gp, bp, rrp, rgp, rbp;

        rp = (double) src[ (x * 4) + 0 ] / 255.0;
        gp = (double) src[ (x * 4) + 1 ] / 255.0;
        bp = (double) src[ (x * 4) + 2 ] / 255.0;

        rotate_hue( rp, gp, bp, angle, &rrp, &rgp, &rbp );

        dst[ (x * 4) + 0 ] = (int) ( ( rrp * 255.0 ) + 0.5 );
        dst[ (x * 4) + 1 ] = (int) ( ( rgp * 255.0 ) + 0.5 );
        dst[ (x * 4) + 2 ] = (int) ( ( rbp * 255.0 ) + 0.5 );
        dst[ (x * 4) + 3 ] = src[ (x * 4) + 3 ];
    }
}