/* GIMP - The GNU Image Manipulation Program
 * Copyright (C) 1995 Spencer Kimball and Peter Mattis
 *
 * gimpoperationlchlightnessmode.c
 * Copyright (C) 2015 Thomas Manni <thomas.manni@free.fr>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program 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 General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

 #include "config.h"

#include <cairo.h>
#include <gegl-plugin.h>
#include <gdk-pixbuf/gdk-pixbuf.h>

#include "libgimpcolor/gimpcolor.h"

#include "libgimpmath/gimpmath.h"

#include "operations-types.h"

#include "gimpoperationlchlightnessmode.h"

#define D50_WHITE_REF_X   0.964202880
#define D50_WHITE_REF_Y   1.000000000
#define D50_WHITE_REF_Z   0.824905400
#define LAB_EPSILON (216.0 / 24389.0)
#define LAB_KAPPA   (24389.0 / 27.0)

static inline void gimp_rgb_to_xyz          (const GimpRGB  *rgb,
                                             GimpXYZ        *xyz);
static inline void gimp_xyz_to_rgb          (const GimpXYZ  *xyz,
                                             GimpRGB        *rgb);

static inline void gimp_xyz_to_lab          (const GimpXYZ *xyz,
                                             GimpLab       *lab);
static inline void gimp_lab_to_xyz          (const GimpLab *lab,
                                             GimpXYZ       *xyz);

static inline void gimp_rgb_to_lab          (const GimpRGB  *rgb,
                                             GimpLab        *lab);
static inline void gimp_lab_to_rgb          (const GimpLab  *lab,
                                             GimpRGB        *rgb);

//
static inline void
gimp_rgb_to_xyz (const GimpRGB *rgb,
                  GimpXYZ       *xyz)
 {
   g_return_if_fail (rgb != NULL);
   g_return_if_fail (xyz != NULL);

   xyz->x = rgb->r * 0.43603516 + rgb->g * 0.38511658 + rgb->b * 0.14305115;
   xyz->y = rgb->r * 0.22248840 + rgb->g * 0.71690369 + rgb->b * 0.06060791;
   xyz->z = rgb->r * 0.01391602 + rgb->g * 0.09706116 + rgb->b * 0.71392822;
   xyz->a = rgb->a;
 }

static inline void
gimp_xyz_to_rgb (const GimpXYZ *xyz,
                  GimpRGB       *rgb)
 {
   g_return_if_fail (xyz != NULL);
   g_return_if_fail (rgb != NULL);   

   rgb->r = xyz->x *  3.134274799724 + xyz->y * -1.617275708956 + xyz->z * -0.490724283042;
   rgb->g = xyz->x * -0.978795575994 + xyz->y *  1.916161689117 + xyz->z *  0.033453331711;
   rgb->b = xyz->x *  0.071976988401 + xyz->y * -0.228984974402 + xyz->z *  1.405718224383;
   rgb->a = xyz->a;
 }

static inline void
gimp_rgb_to_lab (const GimpRGB *rgb,
                 GimpLab       *lab)
{
  GimpXYZ xyz;
  gimp_rgb_to_xyz (rgb, &xyz);
  gimp_xyz_to_lab (&xyz, lab);
}

static inline void
gimp_lab_to_rgb (const GimpLab *lab,
                 GimpRGB       *rgb)
{
  GimpXYZ xyz;
  gimp_lab_to_xyz (lab, &xyz);
  gimp_xyz_to_rgb (&xyz, rgb);
}

 /* intermediate conversions */

static inline void
gimp_xyz_to_lab (const GimpXYZ *xyz,
                 GimpLab       *lab)
{
  gdouble xr, yr, zr;
  gdouble fx, fy, fz;

  g_return_if_fail (xyz != NULL);
  g_return_if_fail (lab != NULL);

  xr = xyz->x / D50_WHITE_REF_X;
  yr = xyz->y / D50_WHITE_REF_Y;
  zr = xyz->z / D50_WHITE_REF_Z;

  if (xr > LAB_EPSILON)
    fx = cbrt (xr);
  else
    fx = (LAB_KAPPA * xr + 16) / 116.0;

  if (yr > LAB_EPSILON)
    fy = cbrt (yr);
  else
    fy = (LAB_KAPPA * yr + 16) / 116.0;

  if (zr > LAB_EPSILON)
    fz = cbrt (zr);
  else
    fz = (LAB_KAPPA * zr + 16) / 116.0;

  lab->l = 116.0 * fy - 16;
  lab->A = 500.0 * (fx - fy);
  lab->B = 200.0 * (fy - fz);
  lab->a = xyz->a;
}

static inline void
gimp_lab_to_xyz (const GimpLab *lab,
                 GimpXYZ       *xyz)
{
  gdouble xr, yr, zr;
  gdouble fx, fy, fz;
  gdouble fxcube, fzcube;

  g_return_if_fail (lab != NULL);
  g_return_if_fail (xyz != NULL);

  fy = (lab->l + 16) / 116.0;
  fx = lab->A / 500.0 + fy;
  fz = fy - lab->B / 200;

  fxcube = pow (fx, 3);
  fzcube = pow (fz, 3);

  if (fxcube > LAB_EPSILON)
    xr = fxcube;
  else
    xr = (116.0 * fx - 16) / LAB_KAPPA;

  if (lab->l > LAB_KAPPA * LAB_EPSILON)
    yr = pow (((lab->l + 16) / 116.0), 3);
  else
    yr = lab->l / LAB_KAPPA;

  if (fzcube > LAB_EPSILON)
    zr = fzcube;
  else
    zr = (116.0 * fz - 16) / LAB_KAPPA;  
  
  xyz->x = xr * D50_WHITE_REF_X;
  xyz->y = yr * D50_WHITE_REF_Y;
  xyz->z = zr * D50_WHITE_REF_Z;
  xyz->a = lab->a;
}
//



static gboolean gimp_operation_lch_lightness_mode_process (GeglOperation       *operation,
                                                           void                *in_buf,
        	                                                 void                *aux_buf,
        	                                                 void                *aux2_buf,
        	                                                 void                *out_buf,
        	                                                 glong                samples,
        	                                                 const GeglRectangle *roi,
        	                                                 gint                 level);
static void     gimp_operation_lch_lightness_mode_prepare (GeglOperation *operation);

G_DEFINE_TYPE (GimpOperationLchLightnessMode, gimp_operation_lch_lightness_mode,
               GIMP_TYPE_OPERATION_POINT_LAYER_MODE)


static void
gimp_operation_lch_lightness_mode_class_init (GimpOperationLchLightnessModeClass *klass)
{
  GeglOperationClass               *operation_class;
  GeglOperationPointComposer3Class *point_class;

  operation_class = GEGL_OPERATION_CLASS (klass);
  point_class     = GEGL_OPERATION_POINT_COMPOSER3_CLASS (klass);

  operation_class->prepare = gimp_operation_lch_lightness_mode_prepare;

  gegl_operation_class_set_keys (operation_class,
                                 "name",        "gimp:lch-lightness-mode",
                                 "description", "GIMP LCH lightness mode operation",
                                 NULL);

  point_class->process = gimp_operation_lch_lightness_mode_process;
}

static void
gimp_operation_lch_lightness_mode_init (GimpOperationLchLightnessMode *self)
{
}

static void
gimp_operation_lch_lightness_mode_prepare (GeglOperation *operation)
{
  /* force linear RGB values (required for XYZ color space conversion) */
  gegl_operation_set_format (operation, "input",  babl_format ("RGBA float"));
  gegl_operation_set_format (operation, "output", babl_format ("RGBA float"));
  gegl_operation_set_format (operation, "aux",    babl_format ("RGBA float"));
  gegl_operation_set_format (operation, "aux2",   babl_format ("Y float"));
}

static gboolean
gimp_operation_lch_lightness_mode_process (GeglOperation       *operation,
	                                 void                *in_buf,
	                                 void                *aux_buf,
	                                 void                *aux2_buf,
	                                 void                *out_buf,
	                                 glong                samples,
	                                 const GeglRectangle *roi,
	                                 gint                 level)
{
  gfloat opacity = GIMP_OPERATION_POINT_LAYER_MODE (operation)->opacity;

  return gimp_operation_lch_lightness_mode_process_pixels (in_buf, aux_buf, aux2_buf, out_buf, opacity, samples, roi, level);
}

gboolean
gimp_operation_lch_lightness_mode_process_pixels (gfloat              *in,
  	                                        gfloat              *layer,
  	                                        gfloat              *mask,
  	                                        gfloat              *out,
  	                                        gfloat               opacity,
  	                                        glong                samples,
  	                                        const GeglRectangle *roi,
  	                                        gint                 level)
{
  const gboolean has_mask = mask != NULL;

  while (samples--)
    {
      GimpLab layer_lab, out_lab;
      GimpRGB layer_rgb = {layer[0], layer[1], layer[2]};
      GimpRGB out_rgb   = {in[0], in[1], in[2]};
      gfloat  comp_alpha, new_alpha;

      comp_alpha = MIN (in[ALPHA], layer[ALPHA]) * opacity;
      if (has_mask)
        comp_alpha *= *mask;

      new_alpha = in[ALPHA] + (1.0 - in[ALPHA]) * comp_alpha;

      if (comp_alpha && new_alpha)
        {
          gint   b;
          gfloat out_tmp[3];
          gfloat ratio = comp_alpha / new_alpha;

          gimp_rgb_to_lab (&layer_rgb, &layer_lab);
          gimp_rgb_to_lab (&out_rgb, &out_lab);

          out_lab.l = layer_lab.l;
          
          gimp_lab_to_rgb (&out_lab, &out_rgb);

          out_tmp[0] = out_rgb.r;
          out_tmp[1] = out_rgb.g;
          out_tmp[2] = out_rgb.b;

          for (b = RED; b < ALPHA; b++)
            {
              out[b] = out_tmp[b] * ratio + in[b] * (1.0 - ratio);
            }
        }
      else
        {
          gint b;

          for (b = RED; b < ALPHA; b++)
            {
              out[b] = in[b];
            }
        }

      out[ALPHA] = in[ALPHA];

      in    += 4;
      layer += 4;
      out   += 4;

      if (has_mask)
        mask++;
    }

  return TRUE;
}
