/* GIMP - The GNU Image Manipulation Program
 * Copyright (C) 1995 Spencer Kimball and Peter Mattis
 *
 * 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 2 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, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 */

#include "config.h"

#include <string.h>

#include <glib-object.h>

#include "libgimpmath/gimpmath.h"

#include "paint-funcs-types.h"

#include "base/tile.h"

#include "base/tile-manager.h"

#include "base/pixel-region.h"

#include "paint-funcs.h"

#include "scale-region.h"

#define EPSILON             (0.0001)  /* arbitary small number for avoiding zero */

/* Tile based scaling */
static void scale  (TileManager           *dstTM, 
                    TileManager           *srcTM,
                    GimpInterpolationType  interpolation);                 


static void inline interpolate_bilinear   (TileManager *srcTM,
                                           gint         x0,
                                           gint         y0,
                                           gint         x1,
                                           gint         y1,
                                           gdouble      xfrac,
                                           gdouble      yfrac,
                                           guchar      *pixel);

static void inline interpolate_nearest    (TileManager *srcTM,
                                           gint         x0,
                                           gint         y0,
                                           gint         x1,
                                           gint         y1,
                                           gdouble      xfrac,
                                           gdouble      yfrac,
                                           guchar      *pixel);

static inline gdouble cubic_spline_fit    (gdouble      dx,
                                           gint         pt0,
                                           gint         pt1,
                                           gint         pt2,
                                           gint         pt3);
static void inline    interpolate_cubic   (TileManager *srcTM,
                                           gint         x0,
                                           gint         y0,
                                           gint         x1,
                                           gint         y1,
                                           gdouble      xfrac,
                                           gdouble      yfrac,
                                           guchar      *pixel);

static gfloat *  create_lanczos2_lookup   (void);
static void inline interpolate_lanczos2   (TileManager *srcTM,
                                           gint         x1,
                                           gint         y1,
                                           gint         x2,
                                           gint         y2,
                                           gdouble      xfrac,
                                           gdouble      yfrac,
                                           guchar      *pixel,
                                           gfloat      *kernel_lookup);
static gfloat *    create_lanczos3_lookup (void);
static void inline interpolate_lanczos3   (TileManager *srcTM,
                                           gint         x1,
                                           gint         y1,
                                           gint         x2,
                                           gint         y2,
                                           gdouble      xfrac,
                                           gdouble      yfrac,
                                           guchar      *pixel,
                                           gfloat      *kernel_lookup);
                                         


                     

static inline gdouble lanczos_sum_mul (guchar         *ptr,
                                       const gdouble  *kernel,    /* 1-D kernel of transform coeffs */
                                       gint            u,
                                       gint            bytes,
                                       gint            byte,
                                       gint            alpha);
static inline gdouble
lanczos_sum (guchar         *ptr,
             const gdouble  *kernel,  /* 1-D kernel of transform coeffs */
             gint            u,
             gint            bytes,
             gint            byte);
static inline gdouble sinc (gdouble x);


void
scale_region  (PixelRegion           *srcPR,
               PixelRegion           *dstPR,
               GimpInterpolationType  interpolation,
               GimpProgressFunc       progress_callback,
               gpointer               progress_data)
{
  TileManager       *tmpTM=NULL;
  TileManager       *srcTM = srcPR->tiles;
  TileManager       *dstTM = dstPR->tiles;
  gint       levelx = 0;
  gint       levely = 0;
  gdouble    scalex = (gdouble) dstPR->w / (gdouble) srcPR->w;
  gdouble    scaley = (gdouble) dstPR->h / (gdouble) srcPR->h;
  gint       width  = srcPR->w;
  gint       height = srcPR->h;
  gint       bytes  = srcPR->bytes;

  /* Copy and return if scale = 1.0 */
  if (srcPR->h == dstPR->h && srcPR->w == dstPR->w)
    {
      copy_region (srcPR, dstPR);
      return;
    }

  /* determine scaling levels */
  while (scalex >= 2)
    {
      scalex  /= 2;
      levelx--;
    }

  while (scaley >= 2)
    {
      scaley  /= 2;
      levely--;
    }

  while (scalex <= 0.5)
    {
      scalex  *= 2;
      levelx++;
    }

  while (scaley <= 0.5)
    {
      scaley  *= 2;
      levely++;
    }

  srcTM = srcPR->tiles;
  while ( levelx < 0 && levely < 0 )   
   {
      width <<=1;
      height <<=1;
      tmpTM = tile_manager_new (width, height, bytes);

      scale (tmpTM, srcTM, interpolation);
      if (srcTM != srcPR->tiles)
        tile_manager_unref (srcTM);

      srcTM = tmpTM;

      levelx++;
      levely++;

   }
  while ( levelx < 0 )
   {
      width <<=1;
      tmpTM = tile_manager_new (width, height, bytes);
      scale (tmpTM, srcTM, interpolation);
      if (srcTM != srcPR->tiles)
        tile_manager_unref (srcTM);
      levelx++;
   }  
  
  while ( levely < 0 )
   {
      height <<=1;
      tmpTM = tile_manager_new (width, height, bytes);
      scale (tmpTM, srcTM, interpolation);
      if (srcTM != srcPR->tiles)
        tile_manager_unref (srcTM);
      levely++;
   }  
  
  while ( levelx > 0 && levely > 0 )   
   {
      width >>=1;
      height >>=1;
      tmpTM = tile_manager_new (width, height, bytes);

      scale (tmpTM, srcTM, interpolation);
      if (srcTM != srcPR->tiles)
        tile_manager_unref (srcTM);
      srcTM = tmpTM;

      levelx--;
      levely--;

   }
  while ( levelx > 0 )
   {
      width <<=1;
      tmpTM = tile_manager_new (width, height, bytes);
      scale (tmpTM, srcTM, interpolation);
      if (srcTM != srcPR->tiles)
        tile_manager_unref (srcTM);
      levelx--;
   }  
  
  while ( levely > 0 )
   {
      height <<=1;
      tmpTM = tile_manager_new (width, height, bytes);
      scale (tmpTM, srcTM, interpolation);
      if (srcTM != srcPR->tiles)
        tile_manager_unref (srcTM);
      levely--;
   }  

  if ( tmpTM == NULL )
    {
       tmpTM = srcTM;
    }

  scale(dstTM, tmpTM, interpolation);
  tile_manager_unref (tmpTM);
  return;
  
}


static inline gdouble
sinc (gdouble x)
{
  gdouble y = x * G_PI;

  if (ABS (x) < LANCZOS_MIN)
    return 1.0;

  return sin (y) / y;
}

/*
 * allocate and fill lookup table of Lanczos windowed sinc function
 * use gfloat since errors due to granularity of array far exceed data precision
 */
gfloat *
create_lanczos_lookup (void)
{
  const gdouble dx = LANCZOS_WIDTH / (gdouble) (LANCZOS_SAMPLES - 1);

  gfloat  *lookup = g_new (gfloat, LANCZOS_SAMPLES);
  gdouble  x      = 0.0;
  gint     i;

  for (i = 0; i < LANCZOS_SAMPLES; i++)
    {
      lookup[i] = ((ABS (x) < LANCZOS_WIDTH) ?
                   (sinc (x) * sinc (x / LANCZOS_WIDTH)) : 0.0);
      x += dx;
    }

  return lookup;
}


gfloat *
create_lanczos2_lookup (void)
{
  const gdouble dx = 2.0 / (gdouble) (LANCZOS_SAMPLES - 1);

  gfloat  *lookup = g_new (gfloat, LANCZOS_SAMPLES);
  gdouble  x      = 0.0;
  gint     i;

  for (i = 0; i < LANCZOS_SAMPLES; i++)
    {
      lookup[i] = ((ABS (x) < 2.0) ?
                   (sinc (x) * sinc (x / 2.0)) : 0.0);
      x += dx;
    }

  return lookup;
}

gfloat *
create_lanczos3_lookup (void)
{
  const gdouble dx = 3.0 / (gdouble) (LANCZOS_SAMPLES - 1);

  gfloat  *lookup = g_new (gfloat, LANCZOS_SAMPLES);
  gdouble  x      = 0.0;
  gint     i;

  for (i = 0; i < LANCZOS_SAMPLES; i++)
    {
      lookup[i] = ((ABS (x) < 3.0) ?
                   (sinc (x) * sinc (x / 3.0)) : 0.0);
      x += dx;
    }

  return lookup;
}


static inline gdouble
lanczos_sum (guchar         *ptr,
             const gdouble  *kernel,  /* 1-D kernel of transform coeffs */
             gint            u,
             gint            bytes,
             gint            byte)
{
  gdouble sum = 0;
  gint    i;

  for (i = 0; i < LANCZOS_WIDTH2 ; i++)
    sum += kernel[i] * ptr[ (u + i - LANCZOS_WIDTH) * bytes + byte ];

  return sum;
}

static inline gdouble
lanczos_sum_mul (guchar         *ptr,
                 const gdouble  *kernel,    /* 1-D kernel of transform coeffs */
                 gint            u,
                 gint            bytes,
                 gint            byte,
                 gint            alpha)
{
  gdouble sum = 0;
  gint    i;

  for (i = 0; i < LANCZOS_WIDTH2 ; i++ )
    sum += kernel[i] * ptr[ (u + i - LANCZOS_WIDTH) * bytes + byte ]
                 * ptr[ (u + i - LANCZOS_WIDTH) * bytes + alpha];

  return sum;
}

/* Catmull-Rom spline - not bad
  * basic intro http://www.mvps.org/directx/articles/catmull/
  * This formula will calculate an interpolated point between pt1 and pt2
  * dx=0 returns pt1; dx=1 returns pt2
  */

static inline gdouble
cubic_spline_fit (gdouble dx,
                  gint    pt0,
                  gint    pt1,
                  gint    pt2,
                  gint    pt3)
{

  return (gdouble) ((( ( - pt0 + 3 * pt1 - 3 * pt2 + pt3 ) *   dx +
      ( 2 * pt0 - 5 * pt1 + 4 * pt2 - pt3 ) ) * dx +
      ( - pt0 + pt2 ) ) * dx + (pt1 + pt1) ) / 2.0;
}


static void inline
interpolate_bilinear (TileManager *srcTM,
                      gint         x0,
                      gint         y0,
                      gint         x1,
                      gint         y1,
                      gdouble      xfrac,
                      gdouble      yfrac,
                      guchar      *pixel)
{
  gint   b;  
  gint   src_bpp         = tile_manager_bpp  (srcTM);
  guchar pixel1[src_bpp];
  guchar pixel2[src_bpp];
  guchar pixel3[src_bpp];
  guchar pixel4[src_bpp];  
  gdouble sum;
  read_pixel_data_1 (srcTM, x0, y0, pixel1);
  read_pixel_data_1 (srcTM, x1, y0, pixel2);
  read_pixel_data_1 (srcTM, x0, y1, pixel3);
  read_pixel_data_1 (srcTM, x1, y1, pixel4);

  for (b=0 ; b < src_bpp ; b++)
    {
      sum =  (gint)((1 - yfrac) * ((1 - xfrac) * pixel1[b] + xfrac * pixel2[b])
                                + yfrac  * ((1 - xfrac) * pixel3[b] + xfrac * pixel4[b]));
      pixel[b] = (guchar) CLAMP (sum, 0, 255);
    }
}

static void inline
interpolate_nearest (TileManager *srcTM,
                     gint         x0,
                     gint         y0,
                     gint         x1,
                     gint         y1,
                     gdouble      xfrac,
                     gdouble      yfrac,
                     guchar      *pixel)
{
  gint   x = ( xfrac <= 0.5) ? x0 : x1;
  gint   y = ( yfrac <= 0.5) ? y0 : y1;           
  read_pixel_data_1 (srcTM, x, y, pixel);
}

static void inline
interpolate_cubic    (TileManager *srcTM,
                      gint         x1,
                      gint         y1,
                      gint         x2,
                      gint         y2,
                      gdouble      xfrac,
                      gdouble      yfrac,
                      guchar      *pixel)
{
  gint   src_bpp    = tile_manager_bpp  (srcTM);
  guint  src_width  = tile_manager_width  (srcTM);
  guint  src_height = tile_manager_height (srcTM);
  gint   b, i;  
  gint   x, y;
  gint   x0;
  gint   y0;

  gdouble    p0, p1, p2, p3;

  guchar pixels[16*src_bpp];
  gdouble sum;

  for (y=y1-1, i=0  ; y <= y1+2 ; y++)
   { 
     for (x=x1-1  ; x <= x1+2 ; x++, i++)
      {
        x0 = (x < src_width) ? ABS(x) : 2*src_width-x;
        y0 = (y < src_height) ? ABS(y): 2*src_height-y;        
        read_pixel_data_1 (srcTM, x0, y0, pixels+(i*src_bpp));
      }
    }
   for (b = 0 ; b < src_bpp  ; b++)
     {
        p0  = cubic_spline_fit (xfrac, pixels[0*src_bpp+b],
                                       pixels[1*src_bpp+b],
                                       pixels[2*src_bpp+b],
                                       pixels[3*src_bpp+b]);
        p1  = cubic_spline_fit (xfrac, pixels[4*src_bpp+b],
                                       pixels[5*src_bpp+b],
                                       pixels[6*src_bpp+b],
                                       pixels[7*src_bpp+b]);
        p2  = cubic_spline_fit (xfrac, pixels[8*src_bpp+b],
                                       pixels[9*src_bpp+b],
                                       pixels[10*src_bpp+b],
                                       pixels[11*src_bpp+b]);
        p3  = cubic_spline_fit (xfrac, pixels[12*src_bpp+b],
                                       pixels[13*src_bpp+b],
                                       pixels[14*src_bpp+b],
                                       pixels[15*src_bpp+b]);
        sum = cubic_spline_fit (yfrac, p0, p1, p2, p3);                                          
        pixel[b]  = (guchar) CLAMP (sum, 0, 255);
     }
}

static void inline
interpolate_lanczos2 (TileManager *srcTM,
                      gint         x1,
                      gint         y1,
                      gint         x2,
                      gint         y2,
                      gdouble      xfrac,
                      gdouble      yfrac,
                      guchar      *pixel,
                      gfloat      *kernel_lookup)
{
  gint   src_bpp    = tile_manager_bpp  (srcTM);
  guint  src_width  = tile_manager_width  (srcTM);
  guint  src_height = tile_manager_height (srcTM);
  gint   b, i, pos;  
  gint   x, y;
  gint   x0;
  gint   y0;
  gint     x_shift, y_shift;
  gdouble  kx_sum, ky_sum;
  gdouble        x_kernel[4],
                 y_kernel[4];

  guchar pixels[6*src_bpp];
  gdouble sum, newval;

  for (y=y1-1, i=0  ; y <= y1+2 ; y++)
   { 
     for (x=x1-1  ; x <= x1+2 ; x++, i++)
      {
        x0 = (x < src_width)  ? ABS(x) : 2*src_width-x;
        y0 = (y < src_height) ? ABS(y) : 2*src_height-y;        
        read_pixel_data_1 (srcTM, x0, y0, pixels+(i*src_bpp));
      }
    }
    
  x_shift = (gint)(xfrac * LANCZOS_SPP + 0.5);
  y_shift = (gint)(yfrac * LANCZOS_SPP + 0.5);
  kx_sum = ky_sum = 0.0;
  for (i = 2; i >= -1; i--)
   {
      gint pos = i * LANCZOS_SPP;
      kx_sum += x_kernel[1 + i] = kernel_lookup[ABS (x_shift - pos)];
      ky_sum += y_kernel[1 + i] = kernel_lookup[ABS (y_shift - pos)];
   }

   /* normalise the kernel arrays */
   for (i = -1; i <= 2; i++)
    {
       x_kernel[1 + i] /= kx_sum;
       y_kernel[1 + i] /= ky_sum;
    }        
   for (b = 0 ; b < src_bpp  ; b++)
    {
       newval = 0.0;
       for (y = 0; y < 4; y++ )
        {
          sum = 0;
          for (x = 0; x < 4; x++ )
           {
              pos = (y * 4 + x) * src_bpp;
              sum+=x_kernel[x] * pixels[pos+b];
           }
          sum *= y_kernel[y];
          newval+=sum;
        }                        
       pixel[b]  = (guchar) CLAMP ((gint) newval, 0, 255);
    }
}

static void inline
interpolate_lanczos3 (TileManager *srcTM,
                      gint         x1,
                      gint         y1,
                      gint         x2,
                      gint         y2,
                      gdouble      xfrac,
                      gdouble      yfrac,
                      guchar      *pixel,
                      gfloat      *kernel_lookup)
{
  gint   src_bpp    = tile_manager_bpp  (srcTM);
  guint  src_width  = tile_manager_width  (srcTM);
  guint  src_height = tile_manager_height (srcTM);
  gint   b, i, pos;  
  gint   x, y;
  gint   x0;
  gint   y0;
  gint     x_shift, y_shift;
  gdouble  kx_sum, ky_sum;
  gdouble        x_kernel[6],
                 y_kernel[6];

  guchar pixels[36*src_bpp];
  gdouble sum, newval;

  for (y=y1-2, i=0  ; y <= y1+3 ; y++)
   { 
     for (x=x1-2  ; x <= x1+3 ; x++, i++)
      {
        x0 = (x < src_width)  ? ABS(x) : 2*src_width-x;
        y0 = (y < src_height) ? ABS(y) : 2*src_height-y;        
        read_pixel_data_1 (srcTM, x0, y0, pixels+(i*src_bpp));
      }
    }
    
  x_shift = (gint)(xfrac * LANCZOS_SPP + 0.5);
  y_shift = (gint)(yfrac * LANCZOS_SPP + 0.5);
  kx_sum = ky_sum = 0.0;
  for (i = 3; i >= -2; i--)
   {
      gint pos = i * LANCZOS_SPP;
      kx_sum += x_kernel[2 + i] = kernel_lookup[ABS (x_shift - pos)];
      ky_sum += y_kernel[2 + i] = kernel_lookup[ABS (y_shift - pos)];
   }

   /* normalise the kernel arrays */
   for (i = -2; i <= 3; i++)
    {
       x_kernel[2 + i] /= kx_sum;
       y_kernel[2 + i] /= ky_sum;
    }        
   for (b = 0 ; b < src_bpp  ; b++)
    {
       newval = 0.0;
       for (y = 0; y < 6; y++ )
        {
          sum = 0;
          for (x = 0; x < 6; x++ )
           {
              pos = (y * 6 + x) * src_bpp;
              sum+=x_kernel[x] * pixels[pos+b];
           }
          sum *= y_kernel[y];
          newval+=sum;
        }                        
       pixel[b]  = (guchar) CLAMP ((gint) newval, 0, 255);
    }
}



static void
scale (TileManager           *dstTM, 
       TileManager           *srcTM,
       GimpInterpolationType  interpolation)
{

  guint              src_width      = tile_manager_width  (srcTM);
  guint              src_height     = tile_manager_height (srcTM);
  Tile              *dst_tile;
  guchar            *dst_data;
  guint              dst_width      = tile_manager_width  (dstTM);
  guint              dst_height     = tile_manager_height (dstTM);
  guint              dst_bpp        = tile_manager_bpp (dstTM);
  guint              dst_tilerows   = tile_manager_tiles_per_row(dstTM);    /*  the number of tiles in each row      */
  guint              dst_tilecols   = tile_manager_tiles_per_col(dstTM);    /*  the number of tiles in each columns  */
  guint              dst_ewidth;
  guint              dst_eheight;
  guint              dst_stride;    
  gdouble            scalex         = (gdouble)dst_width  / (gdouble)src_width;
  gdouble            scaley         = (gdouble)dst_height / (gdouble)src_height;
  gdouble            xfrac;
  gdouble            yfrac;
  gint               x, y, x0, y0, x1, y1;
  gint               sx0, sy0, sx1, sy1;
  gint               col, row;
  guchar             pixel[dst_bpp];
  gfloat            *kernel_lookup=NULL;

  /* if scale is 2^n */  
  if (src_width==dst_width && src_height==dst_height)
    {
      for (row=0 ; row < dst_tilerows ; row++)
       {
         for (col=0 ; col < dst_tilecols ; col++)
          {            
            dst_tile    = tile_manager_get_at (dstTM, col, row, TRUE, TRUE);
            dst_data    = tile_data_pointer (dst_tile, 0, 0);
            dst_bpp     = tile_bpp (dst_tile);
            dst_ewidth  = tile_ewidth (dst_tile);
            dst_eheight = tile_eheight (dst_tile);
            dst_stride  = dst_ewidth * dst_bpp;          
            x0 = col * TILE_WIDTH;
            y0 = row * TILE_HEIGHT;                      
            x1 = x0 + dst_ewidth-1;
            y1 = y0 + dst_eheight-1;
            read_pixel_data (srcTM, x0, y0, x1, y1, dst_data, dst_stride);
          }
        }
        return;
    }

  if (interpolation == GIMP_INTERPOLATION_LANCZOS)
    kernel_lookup = create_lanczos3_lookup();
  
  for (row=0 ; row < dst_tilerows ; row++)
   {
     for (col=0 ; col < dst_tilecols ; col++)
      {
        dst_tile    = tile_manager_get_at (dstTM, col, row, TRUE, FALSE);
        dst_data    = tile_data_pointer (dst_tile, 0, 0);
        dst_bpp     = tile_bpp (dst_tile);
        dst_ewidth  = tile_ewidth (dst_tile);
        dst_eheight = tile_eheight (dst_tile);
        dst_stride  = dst_ewidth * dst_bpp;          

        x0 = col * TILE_WIDTH;
        y0 = row * TILE_HEIGHT;                      
        x1 = x0 + dst_ewidth-1;
        y1 = y0 + dst_eheight-1;

        for (y=y0 ; y <= y1 ; y++)
         {
            yfrac = ( y / scaley );
            sy0 = (gint) floor (yfrac);
            sy1 = (gint) ceil (yfrac);
            sy1 = ( sy1 >= src_height) ? src_height-1 : sy1;      
            yfrac =  yfrac - sy0;

            for (x=x0 ; x <= x1 ; x++)
             {
                xfrac = (x / scalex);
                sx0 = (gint) floor (xfrac);
                sx1 = (gint) ceil (xfrac);
                sx1 = ( sx1 >= src_width) ? src_width-1 : sx1;
                xfrac =  xfrac - sx0;
                
                switch (interpolation)
                 {
                    case GIMP_INTERPOLATION_NONE:
                      interpolate_nearest (srcTM, sx0, sy0,
                                                  sx1, sy1,
                                                  xfrac, yfrac,
                                                  pixel);
                      break;   
                    case GIMP_INTERPOLATION_LINEAR:
                      interpolate_bilinear (srcTM, sx0, sy0,
                                                   sx1, sy1,
                                                   xfrac, yfrac,
                                                   pixel);
                      break;   
                    case GIMP_INTERPOLATION_CUBIC:
                      interpolate_cubic (srcTM, sx0, sy0,
                                                sx1, sy1,
                                                xfrac, yfrac,
                                                pixel);
                      break;
                    case GIMP_INTERPOLATION_LANCZOS:
                      interpolate_lanczos3 (srcTM, sx0, sy0,
                                                   sx1, sy1,
                                                   xfrac, yfrac,
                                                   pixel,
                                                   kernel_lookup);
                      break;   
                 }           
               write_pixel_data_1 (dstTM, x, y, pixel);
            }
         }
      }
   }
  if (interpolation == GIMP_INTERPOLATION_LANCZOS)
    g_free(kernel_lookup);

}

