/*
  Copyright (C) 2016 Mark Oteiza <mvoteiza@udel.edu>

  Author: Mark Oteiza <mvoteiza@udel.edu>
  Created: 01 March 2016

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License version 2 as published by the Free Software Foundation.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public
  License along with this library. If not, see
  <http://www.gnu.org/licenses/>.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <X11/Xlib.h>
#include <lcms2.h>

#include "emacs-module.h"

#define UNUSED __attribute__((unused))

int plugin_is_GPL_compatible;

static emacs_value Flcms2_ciede2000 (emacs_env *env, ptrdiff_t UNUSED argc,
                                     emacs_value argv[], void UNUSED *data) {
    emacs_value ret;
    double L1 = env->extract_float(env, argv[0]);
    double a1 = env->extract_float(env, argv[1]);
    double b1 = env->extract_float(env, argv[2]);
    double L2 = env->extract_float(env, argv[3]);
    double a2 = env->extract_float(env, argv[4]);
    double b2 = env->extract_float(env, argv[5]);
    double kL = env->extract_float(env, argv[6]);
    double kC = env->extract_float(env, argv[7]);
    double kH = env->extract_float(env, argv[8]);

    const cmsCIELab lab1 = { .L = L1, .a = a1, .b = b1 };
    const cmsCIELab lab2 = { .L = L2, .a = a2, .b = b2 };
    ret = env->make_float(env, cmsCIE2000DeltaE (&lab1, &lab2, kL, kC, kH));
    return ret;
}

static emacs_value Flcms2_ciecamde02 (emacs_env *env, ptrdiff_t UNUSED argc,
                                      emacs_value argv[], void UNUSED *data) {
    emacs_value ret;
    cmsViewingConditions vc;
    cmsJCh jch1, jch2;
    cmsHANDLE h1, h2;
    /* scale XYZ because emacs funs expect these correlates to be in
       the unit line segment [0,1] */
    double X1 = 100 * env->extract_float(env, argv[0]);
    double Y1 = 100 * env->extract_float(env, argv[1]);
    double Z1 = 100 * env->extract_float(env, argv[2]);
    double X2 = 100 * env->extract_float(env, argv[3]);
    double Y2 = 100 * env->extract_float(env, argv[4]);
    double Z2 = 100 * env->extract_float(env, argv[5]);
    /* printf("(%f, %f, %f) <-> (%f, %f, %f)\n", X1, Y1, Z1, X2, Y2, Z2); */
    double Mp1, Mp2, FL, k;
    double Jp1, ap1, bp1, Jp2, ap2, bp2;
    /* UCS coefficients */
    /* double KL = 0.77; */
    /* double c1 = 0.007; */
    /* double c2 = 0.0228; */

    vc.whitePoint.X = 95.047;
    vc.whitePoint.Y = 100.00;
    vc.whitePoint.Z = 108.883;
    vc.Yb = 20;
    vc.La = 100;
    vc.surround = AVG_SURROUND;
    vc.D_value = 1.0;

    h1 = cmsCIECAM02Init(0, &vc);
    h2 = cmsCIECAM02Init(0, &vc);
    const cmsCIEXYZ xyz1 = { .X = X1, .Y = Y1, .Z = Z1 };
    const cmsCIEXYZ xyz2 = { .X = X2, .Y = Y2, .Z = Z2 };

    cmsCIECAM02Forward(h1, &xyz1, &jch1);
    cmsCIECAM02Forward(h2, &xyz2, &jch2);

    cmsCIECAM02Done(h1);
    cmsCIECAM02Done(h2);
    /* Now have JCh, need to calculate Jab

       M = C * F_L^0.25
       J' = 1.7 J / (1 + 0.007 J)
       M' = 43.86 ln(1 + 0.0228 M)
       a' = M' cos(h)
       b' = M' sin(h)

       where

       F_L = 0.2 k^4 (5 L_A) + 0.1 (1 - k^4)^2 (5 L_A)^(1/3),
       k = 1/(5 L_A + 1)
     */
    k = 1.0 / (1.0 + (5.0 * vc.La));
    FL = pow(k, 4) * vc.La + 0.1 * pow(1 - pow(k, 4), 2) * pow(5 * vc.La, 1.0 / 3.0);
    Mp1 = 43.86 * log(1.0 + 0.0228 * (jch1.C * pow(FL, 0.25)));
    Mp2 = 43.86 * log(1.0 + 0.0228 * (jch2.C * pow(FL, 0.25)));
    Jp1 = 1.7 * jch1.J / (1.0 + 0.007 * jch1.J);
    Jp2 = 1.7 * jch2.J / (1.0 + 0.007 * jch2.J);
    ap1 = Mp1 * cos(jch1.h);
    ap2 = Mp2 * cos(jch2.h);
    bp1 = Mp1 * sin(jch1.h);
    bp2 = Mp2 * sin(jch2.h);
    ret = env->make_float(env, sqrt(pow(Jp2 - Jp1, 2.0) +
                                    pow(ap2 - ap1, 2.0) +
                                    pow(bp2 - bp1, 2.0)));
    return ret;
}

static void bind_function(emacs_env *env, const char *name, emacs_value Sfun) {
    emacs_value Qfset = env->intern(env, "fset");
    emacs_value Qsym = env->intern(env, name);
    emacs_value args[] = { Qsym, Sfun };

    env->funcall(env, Qfset, 2, args);
}

static void provide(emacs_env *env, const char *feature) {
    emacs_value Qfeat = env->intern(env, feature);
    emacs_value Qprovide = env->intern (env, "provide");
    emacs_value args[] = { Qfeat };

    env->funcall(env, Qprovide, 1, args);
}

int emacs_module_init(struct emacs_runtime *ert) {
    emacs_env *env = ert->get_environment(ert);

    bind_function(env, "lcms2-ciede2000-internal",
                  env->make_function(env, 9, 9, Flcms2_ciede2000,
                                     "Compute CIEDE2000 between two colors.\
\n(fn L1 A1 B1 L2 A2 B2 kL kC kH)", NULL));
    bind_function(env, "lcms2-ciecamde02-internal",
                  env->make_function(env, 6, 6, Flcms2_ciecamde02,
                                     "Compute CIECAM02 between two colors.\
\n(fn X1 Y1 Z1 X2 Y2 Z2)", NULL));
    provide(env, "lcms2");
    return EXIT_SUCCESS;
}