/* 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 (C) 1997 Hirotsuna Mizuno <s1041150@u-aizu.ac.jp>
 * Copyright (C) 2011 Robert Sasu <sasu.robert@gmail.com>
 */

#include "config.h"
#include <glib/gi18n-lib.h>

#ifdef GEGL_CHANT_PROPERTIES

gegl_chant_string (fractal, _("Fractal"), "fractal",
                   _("Type of fractal to use. "
                     "Choices are julia, mandelbrot. Default is mandelbrot."))
gegl_chant_double (X1, _("X1"), -50.0, 50.0, -1.00,
                   _("X1 value, position"))
gegl_chant_double (X2, _("X2"), -50.0, 50.0, 0.50,
                   _("X2 value, position"))
gegl_chant_double (Y1, _("Y1"), -50.0, 50.0, -1.00,
                   _("X2 value, position"))
gegl_chant_double (Y2, _("Y2"), -50.0, 50.0, 1.00,
                   _("Y2 value, position"))
gegl_chant_double (JX, _("JX"), -50.0, 50.0, 0.5,
                   _("Julia seed X value, position"))
gegl_chant_double (JY, _("JY"), -50.0, 50.0, 0.5,
                   _("Julia seed Y value, position"))
gegl_chant_int    (depth, _("Depth"), 1, 65536, 3,
                   _("Depth value"))
gegl_chant_double (bailout, _("Bailout"), 0.0, G_MAXDOUBLE, G_MAXDOUBLE,
                   _("Bailout length"))
gegl_chant_string (background, _("Background"), "wrap",
                   _("Optional parameter to override automatic selection of wrap background. "
                     "Choices are wrap, black, white and transparent."))

#else

#define GEGL_CHANT_TYPE_FILTER
#define GEGL_CHANT_C_FILE       "fractal-trace.c"

#include "gegl-chant.h"
#include <math.h>

typedef enum
{
  BACKGROUND_TYPE_WRAP,
  BACKGROUND_TYPE_TRANSPARENT,
  BACKGROUND_TYPE_BLACK,
  BACKGROUND_TYPE_WHITE
} BackgroundType;

typedef enum
{
  FRACTAL_TYPE_MANDELBROT,
  FRACTAL_TYPE_JULIA
} FractalType;

static void prepare (GeglOperation *operation)
{
  gegl_operation_set_format (operation, "input", babl_format ("RGBA float"));
  gegl_operation_set_format (operation, "output", babl_format ("RGBA float"));
}

static void
julia (gdouble  x,
       gdouble  y,
       gdouble  jx,
       gdouble  jy,
       gdouble *u,
       gdouble *v,
       gint     depth,
       gdouble  bailout2)
{
  gint    i;
  gdouble xx = x;
  gdouble yy = y;

  for (i = 0; i < depth; i++)
    {
      gdouble x2, y2, tmp;

      x2 = xx * xx;
      y2 = yy * yy;
      tmp = x2 - y2 + jx;
      yy  = 2 * xx * yy + jy;
      xx  = tmp;

      if ((x2 + y2) > bailout2)
	break;
    }

  *u = xx;
  *v = yy;
}

static void
fractaltrace (GeglBuffer          *input,
              const GeglRectangle *picture,
              gfloat              *dst_buf,
              const GeglRectangle *roi,
              GeglChantO          *o,
              gint                 y,
              FractalType          fractal_type,
              BackgroundType       background_type,
              const Babl          *format)
{
  GeglMatrix2  scale;        /* a matrix indicating scaling factors around the
                                current center pixel.
                              */
  gint         x, i, offset;
  gdouble      scale_x, scale_y;
  gdouble      bailout2;
  gfloat       dest[4];

  scale_x = (o->X2 - o->X1) / picture->width;
  scale_y = (o->Y2 - o->Y1) / picture->height;

  bailout2 = o->bailout * o->bailout;

  offset = (y - roi->y) * roi->width * 4;

  for (x = roi->x; x < roi->x + roi->width; x++)
    {
      gdouble cx, cy;
      gdouble px, py;
      dest[1] = dest[2] = dest[3] = dest[0] = 0.0;

      switch (fractal_type)
        {
        case FRACTAL_TYPE_JULIA:
#define gegl_unmap(u,v,ud,vd) {\
       gdouble rx, ry;\
       cx = o->X1 + ((u) - picture->x) * scale_x; \
       cy = o->Y1 + ((v) - picture->y) * scale_y; \
       julia (cx, cy, o->JX, o->JY, &rx, &ry, o->depth, bailout2);\
       ud = (rx - o->X1) / scale_x + picture->x;\
       vd = (ry - o->Y1) / scale_y + picture->y;\
      }
      gegl_sampler_compute_scale (scale, x, y);
      gegl_unmap(x,y,px,py);
#undef gegl_unmap
          break;

        case FRACTAL_TYPE_MANDELBROT:
#define gegl_unmap(u,v,ud,vd) {\
           gdouble rx, ry;\
           cx = o->X1 + ((u) - picture->x) * scale_x; \
           cy = o->Y1 + ((v) - picture->y) * scale_y; \
           julia (cx, cy, cx, cy, &rx, &ry, o->depth, bailout2);\
           ud = (rx - o->X1) / scale_x + picture->x;\
           vd = (ry - o->Y1) / scale_y + picture->y;\
         }
      gegl_sampler_compute_scale (scale, x, y);
      gegl_unmap(x,y,px,py);
#undef gegl_unmap
          break;

        default:
          g_error (_("Unsupported fractal type"));
        }

      if (0 <= px && px < picture->width && 0 <= py && py < picture->height)
        {
          gegl_buffer_sample (input, px, py, &scale, dest, format,
                              GEGL_SAMPLER_LOHALO, GEGL_ABYSS_NONE);
        }
      else
        {
          switch (background_type)
            {
            case BACKGROUND_TYPE_WRAP:
              px = fmod (px, picture->width);
              py = fmod (py, picture->height);

              if (px < 0.0)
                px += picture->width;
              if (py < 0.0)
                py += picture->height;

              /* Check for NaN */
              if (isnan (px))
                {
                  if (signbit (px))
                    px = 0.0;
                  else
                    px = picture->width - 1.0;
                }

              if (isnan (py))
                {
                  if (signbit (py))
                    py = 0.0;
                  else
                    py = picture->height - 1.0;
                }

              gegl_buffer_sample (input, px, py, &scale, dest, format,
                                  GEGL_SAMPLER_LOHALO, GEGL_ABYSS_NONE);
              break;

            case BACKGROUND_TYPE_TRANSPARENT:
              dest[0] = dest[1] = dest[2] = dest[3] = 0.0;
              break;

            case BACKGROUND_TYPE_BLACK:
              dest[0] = dest[1] = dest[2] = 0.0;
              dest[3] = 1.0;
              break;

            case BACKGROUND_TYPE_WHITE:
              dest[0] = dest[1] = dest[2] = dest[3] = 1.0;
              break;
            }
        }

      for (i = 0; i < 4; i++)
        dst_buf[offset++] = dest[i];
    }
}

static gboolean
process (GeglOperation       *operation,
         GeglBuffer          *input,
         GeglBuffer          *output,
         const GeglRectangle *result,
         gint                 level)
{
  GeglChantO	*o;
  GeglRectangle  boundary;
  const Babl    *format;
  FractalType    fractal_type;
  BackgroundType background_type;
  gfloat        *dst_buf;
  gint           y;

  o = GEGL_CHANT_PROPERTIES (operation);
  boundary = gegl_operation_get_bounding_box (operation);

  fractal_type = FRACTAL_TYPE_MANDELBROT;
  if (!strcmp (o->fractal, "mandelbrot"))
    fractal_type = FRACTAL_TYPE_MANDELBROT;
  else if (!strcmp(o->fractal, "julia"))
    fractal_type = FRACTAL_TYPE_JULIA;

  background_type = BACKGROUND_TYPE_WRAP;
  if (!strcmp (o->background, "wrap"))
    background_type = BACKGROUND_TYPE_WRAP;
  else if (!strcmp (o->background, "transparent"))
    background_type = BACKGROUND_TYPE_TRANSPARENT;
  else if (!strcmp (o->background, "black"))
    background_type = BACKGROUND_TYPE_BLACK;
  else if (!strcmp (o->background, "white"))
    background_type = BACKGROUND_TYPE_WHITE;

  format = babl_format ("RGBA float");
  dst_buf = g_new0 (gfloat, result->width * result->height * 4);

  for (y = result->y; y < result->y + result->height; y++)
    fractaltrace (input, &boundary, dst_buf, result, o, y,
                  fractal_type, background_type, format);

  gegl_buffer_set (output, result, 0, format, dst_buf, GEGL_AUTO_ROWSTRIDE);

  g_free (dst_buf);

  gegl_buffer_sample_cleanup (input);

  return TRUE;
}

static GeglRectangle
get_bounding_box (GeglOperation *operation)
{
  GeglRectangle  result = {0,0,0,0};
  GeglRectangle *in_rect;

  in_rect = gegl_operation_source_get_bounding_box (operation, "input");
  if (!in_rect)
    return result;

  return *in_rect;
}

static GeglRectangle
get_required_for_output (GeglOperation       *operation,
                         const gchar         *input_pad,
                         const GeglRectangle *roi)
{
  return get_bounding_box (operation);
}

static void
gegl_chant_class_init (GeglChantClass *klass)
{
  GeglOperationClass       *operation_class;
  GeglOperationFilterClass *filter_class;

  operation_class = GEGL_OPERATION_CLASS (klass);
  filter_class    = GEGL_OPERATION_FILTER_CLASS (klass);

  filter_class->process                    = process;
  operation_class->prepare                 = prepare;
  operation_class->get_bounding_box        = get_bounding_box;
  operation_class->get_required_for_output = get_required_for_output;

  gegl_operation_class_set_keys (operation_class,
    "categories"  , "map",
    "name"        , "gegl:fractal-trace",
    "description" , _("Performs fractal trace on the image"),
    NULL);
}

#endif