/* This file is an image processing operation for GEGL

 *

 * GEGL is free software; you can redistribute it and/or

 * modify it under the terms of the GNU Lesser General Public

 * License as published by the Free Software Foundation; either

 * version 3 of the License, or (at your option) any later version.

 *

 * GEGL 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 GEGL; if not, see <http://www.gnu.org/licenses/>.

 *

 * Copyright 2011 Michael Muré <batolettre@gmail.com>

 */

#include "config.h"

#include <glib/gi18n-lib.h>

#include <math.h>

#include <stdio.h> /* for test only */

#ifdef GEGL_CHANT_PROPERTIES

gegl_chant_double (strength, _("Strength"), 0.0, 100.0, 50,

                   _("Effect Strength"))

gegl_chant_double (size, _("Size"), 1.0, 10000.0, 40.0,

                   _("Effect Size"))

gegl_chant_double (hardness, _("Hardness"), 0.0, 1.0, 0.5,

                   _("Effect Hardness"))

gegl_chant_path (stroke, _("Stroke"), _("Effect Strength"))

gegl_chant_enum (behavior, _("Behavior"), GeglWarpBehavior, GEGL_TYPE_WARP_BEHAVIOR,

                 GEGL_WARP_BEHAVIOR_MOVE, _("Behavior of the op"))

#else

#define GEGL_CHANT_TYPE_FILTER

#define GEGL_CHANT_C_FILE       "warp.c"

#include "gegl-plugin.h"

#include "gegl-path.h"

static void path_changed (GeglPath            *path,

                          const GeglRectangle *roi,

                          gpointer             userdata);

#include "gegl-chant.h"

typedef struct {

  gdouble     *lookup;

  GeglBuffer  *buffer;

  gdouble      last_x;

  gdouble      last_y;

  gboolean     last_point_set;

} WarpPrivate;

static void

path_changed (GeglPath            *path,

              const GeglRectangle *roi,

              gpointer             userdata)

{

  GeglRectangle rect = *roi;

  GeglChantO    *o   = GEGL_CHANT_PROPERTIES (userdata);

  /* invalidate the incoming rectangle */

  rect.x -= o->size/2;

  rect.y -= o->size/2;

  rect.width += o->size;

  rect.height += o->size;

  gegl_operation_invalidate (userdata, &rect, FALSE);

}

static void

prepare (GeglOperation *operation)

{

  GeglChantO  *o = GEGL_CHANT_PROPERTIES (operation);

  WarpPrivate *priv;

  const Babl *format = babl_format_n (babl_type ("float"), 2);

  gegl_operation_set_format (operation, "input", format);

  gegl_operation_set_format (operation, "output", format);

  if (!o->chant_data)

    {

      o->chant_data = g_slice_new (WarpPrivate);

    }

  priv = (WarpPrivate*) o->chant_data;

  priv->last_point_set = FALSE;

  priv->lookup = NULL;

  priv->buffer = NULL;

}

static void

finalize (GObject *object)

{

  GeglChantO *o = GEGL_CHANT_PROPERTIES (object);

  if (o->chant_data)

    {

      g_slice_free (WarpPrivate, o->chant_data);

      o->chant_data = NULL;

    }

  G_OBJECT_CLASS (gegl_chant_parent_class)->finalize (object);

}

static gdouble

gauss (gdouble f)

{

  /* This is not a real gauss function. */

  /* Approximation is valid if -1 < f < 1 */

  if (f < -0.5)

    {

      f = -1.0 - f;

      return (2.0 * f*f);

    }

  if (f < 0.5)

    return (1.0 - 2.0 * f*f);

  f = 1.0 - f;

  return (2.0 * f*f);

}

/* set up lookup table */

static void

calc_lut (GeglChantO  *o)

{

  WarpPrivate  *priv = (WarpPrivate*) o->chant_data;

  gint          length;

  gint          x;

  gdouble       exponent;

  length = ceil (0.5 * o->size + 1.0);

  priv->lookup = g_malloc (length * sizeof (gdouble));

  if ((1.0 - o->hardness) < 0.0000004)

    exponent = 1000000.0;

  else

    exponent = 0.4 / (1.0 - o->hardness);

  for (x = 0; x < length; x++)

    {

      priv->lookup[x] = gauss (pow (2.0 * x / o->size, exponent));

    }

}

static gdouble

get_stamp_force (GeglChantO *o,

                 gdouble     x,

                 gdouble     y)

{

  WarpPrivate  *priv = (WarpPrivate*) o->chant_data;

  gfloat        radius;

  if (!priv->lookup)

    {

      calc_lut (o);

    }

  radius = sqrt(x*x+y*y);

  if (radius < 0.5 * o->size + 1)

    {

      /* linear interpolation */

      gdouble a, ratio;

      gdouble before, after;

      a = floor (radius);

      ratio = (radius - a);

      before = priv->lookup[(gint) a];

      after = priv->lookup[(gint) a + 1];

      return ratio * before + (1.0 - ratio) * after;

    }

  return 0.0;

}

static void

stamp (GeglChantO          *o,

       const GeglRectangle *result,

       gdouble              x,

       gdouble              y)

{

  WarpPrivate         *priv = (WarpPrivate*) o->chant_data;

  GeglBufferIterator  *it;

  const Babl          *format;

  gdouble              influence;

  gdouble              x_mean = 0.0;

  gdouble              y_mean = 0.0;

  gint                 x_iter, y_iter;

  GeglRectangle        area = {x - o->size / 2.0,

                               y - o->size / 2.0,

                               o->size,

                               o->size};

  /* first point of the stroke */

  if (!priv->last_point_set)

    {

      priv->last_x = x;

      priv->last_y = y;

      priv->last_point_set = TRUE;

      return;

    }

  /* don't stamp if outside the roi treated */

  if (!gegl_rectangle_intersect (NULL, result, &area))

    return;

  format = babl_format_n (babl_type ("float"), 2);

  /* If needed, compute the mean deformation */

  if (o->behavior == GEGL_WARP_BEHAVIOR_SMOOTH)

    {

      gint pixel_count = 0;

      it = gegl_buffer_iterator_new (priv->buffer, &area, 0, format, GEGL_BUFFER_READ, GEGL_ABYSS_NONE);

      while (gegl_buffer_iterator_next (it))

        {

          gint    n_pixels    = it->length;

          gfloat *coords      = it->data[0];

          while (n_pixels--)

            {

              x_mean += coords[0];

              y_mean += coords[1];

              coords += 2;

            }

          pixel_count += it->roi->width * it->roi->height;

        }

      x_mean /= pixel_count;

      y_mean /= pixel_count;

    }

  it = gegl_buffer_iterator_new (priv->buffer, &area, 0, format, GEGL_BUFFER_READWRITE, GEGL_ABYSS_NONE);

  while (gegl_buffer_iterator_next (it))

    {

      /* iterate inside the stamp roi */

      gint    n_pixels = it->length;

      gfloat *coords   = it->data[0];

      x_iter = it->roi->x; /* initial x         */

      y_iter = it->roi->y; /* and y coordinates */

      while (n_pixels--)

        {

          influence = 0.01 * o->strength * get_stamp_force (o,

                                                            x_iter - x,

                                                            y_iter - y);

          switch (o->behavior)

            {

              case GEGL_WARP_BEHAVIOR_MOVE:

                coords[0] += influence * (priv->last_x - x);

                coords[1] += influence * (priv->last_y - y);

                break;

              case GEGL_WARP_BEHAVIOR_GROW:

                coords[0] -= 2.0 * influence * (x_iter - x) / o->size;

                coords[1] -= 2.0 * influence * (y_iter - y) / o->size;

                break;

              case GEGL_WARP_BEHAVIOR_SHRINK:

                coords[0] += 2.0 * influence * (x_iter - x) / o->size;

                coords[1] += 2.0 * influence * (y_iter - y) / o->size;

                break;

              case GEGL_WARP_BEHAVIOR_SWIRL_CW:

                coords[0] += 3.0 * influence * (y_iter - y) / o->size;

                coords[1] -= 5.0 * influence * (x_iter - x) / o->size;

                break;

              case GEGL_WARP_BEHAVIOR_SWIRL_CCW:

                coords[0] -= 3.0 * influence * (y_iter - y) / o->size;

                coords[1] += 5.0 * influence * (x_iter - x) / o->size;

                break;

              case GEGL_WARP_BEHAVIOR_ERASE:

                coords[0] *= 1.0 - MIN (influence, 1.0);

                coords[1] *= 1.0 - MIN (influence, 1.0);

                break;

              case GEGL_WARP_BEHAVIOR_SMOOTH:

                coords[0] -= influence * (coords[0] - x_mean);

                coords[1] -= influence * (coords[1] - y_mean);

                break;

            }

          coords += 2;

          /* update x and y coordinates */

          x_iter++;

          if (x_iter >= (it->roi->x + it->roi->width))

            {

              x_iter = it->roi->x;

              y_iter++;

            }

        }

    }

  /* Memorize the stamp location for movement dependant behavior like move */

  priv->last_x = x;

  priv->last_y = y;

}

static gboolean

process (GeglOperation       *operation,

         GeglBuffer          *input,

         GeglBuffer          *output,

         const GeglRectangle *result,

         gint                 level)

{

  GeglChantO          *o = GEGL_CHANT_PROPERTIES (operation);

  WarpPrivate         *priv = (WarpPrivate*) o->chant_data;

  gdouble              dist;

  gdouble              stamps;

  gdouble              spacing = MAX (o->size * 0.01, 0.5); /*1% spacing for starters*/

  GeglPathPoint        prev, next, lerp;

  gulong               i;

  GeglPathList        *event;

  printf("Process %p\n", operation);

  priv->buffer = gegl_buffer_dup (input);

  event = gegl_path_get_path (o->stroke);

  prev = *(event->d.point);

  while (event->next)

    {

      event = event->next;

      next = *(event->d.point);

      dist = gegl_path_point_dist (&next, &prev;;

      stamps = dist / spacing;

      if (stamps < 1)

        {

          stamp (o, result, next.x, next.y);

          prev = next;

        }

      else

        {

          for (i = 0; i < stamps; i++)

            {

              gegl_path_point_lerp (&lerp, &prev, &next, (i * spacing) / dist);

              stamp (o, result, lerp.x, lerp.y);

            }

          prev = lerp;

        }

    }

  /* Affect the output buffer */

  gegl_buffer_copy (priv->buffer, result, output, result);

  gegl_buffer_set_extent (output, gegl_buffer_get_extent (input));

  g_object_unref (priv->buffer);

  /* prepare for the recomputing of the op */

  priv->last_point_set = FALSE;

  return TRUE;

}

static void

gegl_chant_class_init (GeglChantClass *klass)

{

  GObjectClass               *object_class    = G_OBJECT_CLASS (klass);

  GeglOperationClass         *operation_class = GEGL_OPERATION_CLASS (klass);

  GeglOperationFilterClass   *filter_class    = GEGL_OPERATION_FILTER_CLASS (klass);

  object_class->finalize = finalize;

  operation_class->prepare = prepare;

  filter_class->process = process;

  gegl_operation_class_set_keys (operation_class,

  "name"       , "gegl:warp",

  "categories"  , "transform",

  "description" , _("Compute a relative displacement mapping from a stroke"),

  NULL);

}

#endif