/* 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/pixel-region.h"

#include "paint-funcs.h"

#include "scale-region.h"

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

static void  scale_pixel_region (PixelRegion *srcPR,
                                 guchar      *dst,
                                 gint         bytes,
                                 ScaleMode    mode,
                                 GimpInterpolationType  interpolation);
static void  scale_pixel_buffer (const guchar *src,
                                 guchar       *dst,
                                 gint          bytes,
                                 gint          width,
                                 gint          height,
                                 ScaleMode    mode,
                                 GimpInterpolationType  interpolation);
static void  scale_nearest_pr   (PixelRegion  *dstPR,
                                 PixelRegion  *srcPR);
static void  scale_nearest      (PixelRegion  *dstPR,
                                 const guchar *src,
                                 gint          source_w,
                                 gint          source_h,
                                 gint          bytes);
static void  scale_bilinear_pr  (PixelRegion  *dstPR,
                                 PixelRegion  *srcPR);
static void  scale_bilinear     (PixelRegion  *dstPR,
                                 const guchar *src,
                                 gint          source_w,
                                 gint          source_h,
                                 gint          bytes);
static void  scale_cubic_pr     (PixelRegion  *dstPR,
                                 PixelRegion  *srcPR);
static void  scale_cubic        (PixelRegion  *dstPR,
                                 const guchar *src,
                                 gint          source_w,
                                 gint          source_h,
                                 gint          bytes);
static void  scale_lanczos_pr   (PixelRegion  *dstPR,
                                 PixelRegion  *srcPR);
static void  scale_lanczos      (PixelRegion  *dstPR,
                                 const guchar *src,
                                 gint          source_w,
                                 gint          source_h,
                                 gint          bytes);


static void reduce_x_nearest  (guchar            *src[5],
                               guchar            *dst,
                               gint               width,
                               gint               bytes);
static void reduce_x_linear   (guchar            *src[5],
                               guchar            *dst,
                               gint               width,
                               gint               bytes);
static void reduce_x_lanczos2 (guchar            *src[5],
                               guchar            *dst,
                               gint               width,
                               gint               bytes);
static void reduce_y_nearest  (guchar            *src[5],
                               guchar            *dst,
                               gint               width,
                               gint               bytes);
static void reduce_y_linear   (guchar            *src[5],
                               guchar            *dst,
                               gint               width,
                               gint               bytes);
static void reduce_y_lanczos2 (guchar            *src[5],
                               guchar            *dst,
                               gint               width,
                               gint               bytes);


static void expand_x_nearest  (guchar           *src[5],
                               guchar           *dst,
                               gint              width,
                               gint              bytes);
static void expand_x_linear   (guchar           *src[5],
                               guchar           *dst,
                               gint              width,
                               gint              bytes);
static void expand_x_lanczos2 (guchar           *src[5],
                               guchar           *dst,
                               gint              width,
                               gint              bytes);
                              

static void expand_x_nearest  (guchar           *src[5],
                               guchar           *dst,
                               gint              width,
                               gint              bytes);
static void expand_x_linear   (guchar           *src[5],
                               guchar           *dst,
                               gint              width,
                               gint              bytes);
static void expand_x_lanczos2 (guchar           *src[5],
                               guchar           *dst,
                               gint              width,
                               gint              bytes);

static void expand_y_nearest  (guchar           *src[5],
                               guchar           *dst,
                               gint              width,
                               gint              bytes);
static void expand_y_linear   (guchar           *src[5],
                               guchar           *dst,
                               gint              width,
                               gint              bytes);
static void expand_y_lanczos2 (guchar           *src[5],
                               guchar           *dst,
                               gint              width,
                               gint              bytes);
                      

static gboolean inv_lin_trans (const gdouble *t,
                               gdouble       *it);
                               
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 void
rotate_pointers (guchar  **p,
                 guint32   n);

/* Cubic interpolation */
static void expand_x_cubic             (guchar       *src[5],
                                        guchar       *dst,
                                        gint          width,
                                        gint          bytes);
static void expand_y_cubic             (guchar       *src[5],
                                        guchar       *dst,
                                        gint          width,
                                        gint          bytes);
static inline gdouble cubic_spline_fit (gdouble dx,
                                        gint    pt0,
                                        gint    pt1,
                                        gint    pt2,
                                        gint    pt3);


inline void rotate_pointer_1     (guchar *ptr[5]);
inline void rotate_pointer_2     (guchar *ptr[5]);
static inline gdouble sinc (gdouble x);


void
scale_region  (PixelRegion           *srcPR,
               PixelRegion           *dstPR,
               GimpInterpolationType  interpolation,
               GimpProgressFunc       progress_callback,
               gpointer               progress_data)
{
  const gint bytes = dstPR->bytes;
  gint       levelx = 0;
  gint       levely = 0;
  guchar    *dst;
  guchar    *src;
  gint       srcwidth;
  gint       srcheight;
  gint       dstwidth  = 0;
  gint       dstheight = 0;
 /*
  InterpolateKernel expand_nearest    = {1, { 1,  0,  0,  0,  0},  1};
  InterpolateKernel expand_linear     = {2, { 1,  1,  0,  0,  0},  2};
  InterpolateKernel expand_lanczos2   = {4, {-1,  9,  9, -1,  0}, 16};
  InterpolateKernel reduce_nearest    = {1, { 1,  0,  0,  0,  0},  1};
  InterpolateKernel reduce_linear     = {2, { 1,  1,  0,  0,  0},  2};
  InterpolateKernel reduce_lanczos2   = {5, {-1,  9,  16, 9, -1}, 32};
  InterpolateKernel reduce_kernel     = {0, { 0,  0,  0,  0,  0},  1};
  InterpolateKernel expand_kernel     = {0, { 0,  0,  0,  0,  0},  1};
*/  
  gdouble    scalex = (gdouble) dstPR->w / (gdouble) srcPR->w;
  gdouble    scaley = (gdouble) dstPR->h / (gdouble) srcPR->h;

  /* 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++;
    }

  /* Quick step  */
  if (levelx == 0 && levely == 0)
    {
      switch (interpolation)
        {
          case GIMP_INTERPOLATION_NONE:
            scale_nearest_pr (dstPR, srcPR);
            break;      
          case GIMP_INTERPOLATION_LINEAR:
            scale_bilinear_pr (dstPR, srcPR);
            break;           
          case GIMP_INTERPOLATION_CUBIC:
            scale_cubic_pr (dstPR, srcPR);
            break;
          case GIMP_INTERPOLATION_LANCZOS:
            scale_lanczos_pr (dstPR, srcPR);
            break;
        }
      return;
    }



  srcwidth  = srcPR->w;
  srcheight = srcPR->h;

  dst = NULL;                                                     
  /* First expand in x or y if level(x|y) is negative */
  if (levelx < 0 && dst == NULL)
    {
      dstwidth  = srcwidth * 2;
      dstheight = srcheight;
      dst = g_new (guchar, dstwidth * dstheight * bytes);

      scale_pixel_region (srcPR, dst, bytes, EXPAND_X, interpolation);

      levelx++;
      srcwidth  = dstwidth;
      srcheight = dstheight;
    }

  if (levely < 0 && dst == NULL)
    {
      dstwidth  = srcwidth;
      dstheight = srcheight * 2;
      dst = g_new (guchar, dstwidth * dstheight * bytes);

      scale_pixel_region (srcPR, dst, bytes, EXPAND_Y, interpolation);

      levely++;
      srcwidth  = dstwidth;
      srcheight = dstheight;
    }
  /* 
   * While level(x|y) is negative expand alternating x/y direction 
   * Start with y direction because first step could have 
   * been executed for x direction
   */
  while ( levelx < 0 || levely < 0 ) 
    {
      while (levely < 0)
        {
          src = dst;
          dstwidth  = srcwidth;
          dstheight = srcheight * 2;

          dst = g_new (guchar, dstwidth * dstheight * bytes);
          scale_pixel_buffer (src, dst, bytes, srcwidth, srcheight, EXPAND_Y, interpolation);
          g_free(src);

          levely++;
          srcwidth  = dstwidth;
          srcheight = dstheight;
        }
  
      while (levelx < 0)
        {
          src = dst;
          dstwidth  = srcwidth * 2;
          dstheight = srcheight;

          dst = g_new (guchar, dstwidth * dstheight * bytes);
          scale_pixel_buffer (src, dst, bytes, srcwidth, srcheight, EXPAND_X, interpolation);
          g_free(src);

          levelx++;
          srcwidth  = dstwidth;
          srcheight = dstheight;
       }
    }

  /* Reduce in x or y if level(x|y) is positive */
  if (levelx > 0 && dst == NULL)
    {
      dstwidth  = srcwidth / 2;
      dstheight = srcheight;
      dst = g_new (guchar, dstwidth * dstheight * bytes);

      scale_pixel_region (srcPR, dst, bytes, REDUCE_X, interpolation);

      levelx--;
      srcwidth  = dstwidth;
      srcheight = dstheight;
    }
  if (levely > 0 && dst == NULL)
    {
      dstwidth  = srcwidth;
      dstheight = srcheight / 2;
      dst       = g_new (guchar, dstwidth * dstheight * bytes);

      scale_pixel_region (srcPR, dst, bytes, REDUCE_Y, interpolation);

      levely--;        
      srcwidth  = dstwidth;
      srcheight = dstheight;
    }
  /* 
   * While level(x|y) is positive reduce alternating x/y direction 
   * Start with y direction because first step could have 
   * been executed for x direction
   */
    
  while ( levelx > 0 || levely > 0 ) 
    {
      while (levelx > 0)
        {
          src = dst;

          dstwidth  = srcwidth  / 2;
          dstheight = srcheight;

          dst = g_new (guchar, dstwidth * dstheight * bytes);
          scale_pixel_buffer (src, dst, bytes, srcwidth, srcheight, REDUCE_X, interpolation);
          g_free(src);

          levelx--;
          srcwidth  = dstwidth;
          srcheight = dstheight;
       }
    
     while (levely > 0)
       {
         src = dst;
         dstwidth  = srcwidth;
         dstheight = srcheight / 2;

         dst = g_new (guchar, dstwidth * dstheight * bytes);
         scale_pixel_buffer (src, dst, bytes, srcwidth, srcheight, REDUCE_Y, interpolation);
         g_free(src);

         levely--;
         srcwidth  = dstwidth;
         srcheight = dstheight;
       }
    }

  switch (interpolation)
    {
    case GIMP_INTERPOLATION_NONE:
      scale_nearest (dstPR, dst, dstwidth, dstheight, bytes);
      break;      
    case GIMP_INTERPOLATION_LINEAR:
      scale_bilinear (dstPR, dst, dstwidth, dstheight, bytes);
      break;           
    case GIMP_INTERPOLATION_CUBIC:
      scale_cubic (dstPR, dst, dstwidth, dstheight, bytes);
      break;
    case GIMP_INTERPOLATION_LANCZOS:
      scale_lanczos (dstPR, dst, dstwidth, dstheight, bytes);
      break;
    }
  g_free (dst);

}

inline void 
rotate_pointer_1(guchar *ptr[5])
{
  guchar *tmp;
  tmp = ptr[0];
  ptr[0] = ptr[1];
  ptr[1] = ptr[2];
  ptr[2] = ptr[3];
  ptr[3] = ptr[4];
  ptr[4] = tmp;
}

inline void 
rotate_pointer_2(guchar *ptr[5])
{
  guchar *tmp1;
  guchar *tmp2;
  tmp1 = ptr[0];
  tmp2 = ptr[1];
  ptr[0] = ptr[2];
  ptr[1] = ptr[3];
  ptr[2] = ptr[4];
  ptr[3] = tmp1;
  ptr[4] = tmp2;
}


static void
scale_pixel_region (PixelRegion *srcPR,
                    guchar      *dst,
                    gint         bytes,
                    ScaleMode    mode,
                    GimpInterpolationType  interpolation)
{
  guchar *buf;
  guchar *bufptr[5];
  guchar *dstptr;
  void   (*scalefunc)(guchar *[5], guchar *, gint, gint);
  void   (*reducefuncx)(guchar *[5], guchar *, gint, gint);
  void   (*expandfuncx)(guchar *[5], guchar *, gint, gint);
  void   (*reducefuncy)(guchar *[5], guchar *, gint, gint);
  void   (*expandfuncy)(guchar *[5], guchar *, gint, gint);
      
  gint pos, i,y;
  gint step      = 1;
  gint srcwidth  = srcPR->w;
  gint srcstride = srcwidth * bytes;  
  gint dstwidth  = 0;
  gint dststride = 0;

  scalefunc   = (*reduce_x_linear);
  reducefuncx = (*reduce_x_linear);
  reducefuncy = (*reduce_y_linear);
  expandfuncx = (*expand_x_linear);
  expandfuncy = (*expand_y_linear);


  switch (interpolation)
    {
    case GIMP_INTERPOLATION_NONE:
      reducefuncx = (*reduce_x_nearest);
      reducefuncy = (*reduce_y_nearest);
      expandfuncx = (*expand_x_nearest);
      expandfuncy = (*expand_y_nearest);
      break;
      
    case GIMP_INTERPOLATION_LINEAR:
      reducefuncx = (*reduce_x_linear);
      reducefuncy = (*reduce_y_linear);
      expandfuncx = (*expand_x_linear);
      expandfuncy = (*expand_y_linear);
      break;
           
    case GIMP_INTERPOLATION_CUBIC:
      reducefuncx = (*reduce_x_linear);
      reducefuncy = (*reduce_y_linear);
      expandfuncx = (*expand_x_cubic);
      expandfuncy = (*expand_y_cubic);
      break;

    case GIMP_INTERPOLATION_LANCZOS:
      reducefuncx = (*reduce_x_lanczos2);
      reducefuncy = (*reduce_y_lanczos2);
      expandfuncx = (*expand_x_lanczos2);
      expandfuncy = (*expand_y_lanczos2);
      break;
    }
  
  switch (mode)
    {
    case REDUCE_X:
      step = 1;
      dstwidth  = srcwidth/2;
      dststride = dstwidth*bytes;
      scalefunc = (*reducefuncx);
      break;
    case REDUCE_Y:
      step = 2;
      dstwidth  = srcwidth;
      dststride = dstwidth*bytes;
      scalefunc = (*reducefuncy);
      break;
    case EXPAND_X:
      step = 1;
      dstwidth   = srcwidth * 2;
      dststride  = dstwidth*bytes;
      scalefunc = (*expandfuncx);
      break;
    case EXPAND_Y:
      step = 1;
      dstwidth   = srcwidth;
      dststride  = dstwidth*bytes*2;
      scalefunc = (*expandfuncy);
      break;
    }

  buf  = g_new (guchar, srcstride * 5);
  for (i=0 ; i<5 ; i++)
    {
      bufptr[i]=buf+(srcstride*i);
    }
       
  pixel_region_get_row (srcPR, 0, 2, srcwidth, bufptr[0], 1);
  pixel_region_get_row (srcPR, 0, 1, srcwidth, bufptr[1], 1);
  pixel_region_get_row (srcPR, 0, 0, srcwidth, bufptr[2], 1);
  pixel_region_get_row (srcPR, 0, 1, srcwidth, bufptr[3], 1);
  pixel_region_get_row (srcPR, 0, 2, srcwidth, bufptr[4], 1);
  
  dstptr = dst;
  for (y=0 ; y < srcPR->h ; y+=step)
    {  
      (*scalefunc)(bufptr, dstptr, srcwidth, bytes);
      dstptr += dststride;

      if ( step == 2 )
        {
          rotate_pointer_2(bufptr);
          pos = ( y+1 < srcPR->h ) ? y+1 : y-1;  
          pixel_region_get_row (srcPR, 0, pos, srcwidth, bufptr[3], 1);
          pos = ( y+2 < srcPR->h ) ? y+2 : y-2;
          pixel_region_get_row (srcPR, 0, pos, srcwidth, bufptr[4], 1);
        }
      else 
        {
          rotate_pointer_1(bufptr);
          pos = ( y+2 < srcPR->h ) ? y+2 : y-2;  
          pixel_region_get_row (srcPR, 0, pos, srcwidth, bufptr[4], 1);
        }
    }
  g_free (buf);
}


static void
scale_pixel_buffer (const guchar *src,
                    guchar       *dst,
                    gint          bytes,
                    gint          width,
                    gint          height,
                    ScaleMode    mode,
                    GimpInterpolationType  interpolation)
{
  guchar *buf;
  guchar *bufptr[5];
  guchar *dstptr;
  void   (*scalefunc)(guchar *[5], guchar *, gint, gint);
  void   (*reducefuncx)(guchar *[5], guchar *, gint, gint);
  void   (*expandfuncx)(guchar *[5], guchar *, gint, gint);
  void   (*reducefuncy)(guchar *[5], guchar *, gint, gint);
  void   (*expandfuncy)(guchar *[5], guchar *, gint, gint);

  
  gint pos, i, y;
  gint step      = 1;
  gint dstheight = 0;
  gint dstwidth  = 0;
  gint srcstride = width * bytes;
  gint dststride = 0;

  scalefunc   = (*reduce_x_linear);
  reducefuncx = (*reduce_x_linear);
  reducefuncy = (*reduce_y_linear);
  expandfuncx = (*expand_x_linear);
  expandfuncy = (*expand_y_linear);

  switch (interpolation)
    {
    case GIMP_INTERPOLATION_NONE:
      reducefuncx = (*reduce_x_nearest);
      reducefuncy = (*reduce_y_nearest);
      expandfuncx = (*expand_x_nearest);
      expandfuncy = (*expand_y_nearest);
      break;
      
    case GIMP_INTERPOLATION_LINEAR:
      reducefuncx = (*reduce_x_linear);
      reducefuncy = (*reduce_y_linear);
      expandfuncx = (*expand_x_linear);
      expandfuncy = (*expand_y_linear);
      break;
           
    case GIMP_INTERPOLATION_CUBIC:
      reducefuncx = (*reduce_x_linear);
      reducefuncy = (*reduce_y_linear);
      expandfuncx = (*expand_x_cubic);
      expandfuncy = (*expand_y_cubic);
      break;

    case GIMP_INTERPOLATION_LANCZOS:
      reducefuncx = (*reduce_x_lanczos2);
      reducefuncy = (*reduce_y_lanczos2);
      expandfuncx = (*expand_x_lanczos2);
      expandfuncy = (*expand_y_lanczos2);
      break;
    }
  
  switch (mode)
    {
    case REDUCE_X:
      step = 1;
      dstheight = height;
      dstwidth  = width / 2;
      dststride = dstwidth*bytes;
      scalefunc = (*reducefuncx);
      break;
    case REDUCE_Y:
      step = 2;
      dstheight = height / 2;
      dstwidth  = width;
      dststride = dstwidth*bytes;
      scalefunc = (*reducefuncy);
      break;
    case EXPAND_X:
      step = 1;
      dstheight = height;
      dstwidth  = width * 2;
      dststride = dstwidth*bytes;
      scalefunc = (*expandfuncx);
      break;
    case EXPAND_Y:
      step = 1;
      dstheight = height * 2;
      dstwidth  = width;
      dststride = dstwidth*bytes*2;
      scalefunc = (*expandfuncy);
      break;
    }
  
  buf  = g_new (guchar, srcstride * 5);
  for (i=0 ; i<5 ; i++)
    {
      bufptr[i]=buf+(srcstride*i);
    }
  
  memcpy (bufptr[0], src+srcstride*2, srcstride);
  memcpy (bufptr[1], src+srcstride*1, srcstride);
  memcpy (bufptr[2], src,             srcstride);
  memcpy (bufptr[3], src+srcstride*1, srcstride);
  memcpy (bufptr[4], src+srcstride*2, srcstride);
  
  dstptr = dst;
  for (y=0 ; y < height ; y+=step)
    {  
      (*scalefunc)(bufptr, dstptr, width, bytes);
      dstptr += dststride;
      
      if ( step == 2 )
        {
          rotate_pointer_2(bufptr);
          pos = ( y+1 < height ) ? y+1 : y-1;  
          memcpy (bufptr[3], src+srcstride*pos, srcstride);
          pos = ( y+2 < height ) ? y+2 : y-2;
          memcpy (bufptr[4], src+srcstride*pos, srcstride);         
        }
      else 
        {
          rotate_pointer_1(bufptr);
          pos = ( y+2 < height ) ? y+2 : y-2;
          memcpy (bufptr[4], src+srcstride*pos, srcstride);
        }        
    }
  g_free (buf);
}

/* Interpolation dependent expanton functions */
/* Cubic interpolation */
static void
expand_x_cubic (guchar       *src[5],
                guchar       *dst,
                gint          width,
                gint          bytes)
{

  gint          x, b;
  gint          g1,g2,g3,g4;
  gdouble       alphasum;
  gdouble       sum;
  const guchar *src0;
  const guchar *src1;
  const guchar *src2;
  const guchar *src3;
  const guchar *src4;

  src0 = src[0];
  src1 = src[1];
  src2 = src[2];
  src3 = src[3];
  src4 = src[4];

  for (x=0; x < width ; x++)
    {
      g1 = ( x-1 >= 0 ) ? (x-1): (x+1);
      g2 =  x;
      g3 = ( x+1 < width ) ? (x+1): (x-1);
      g4 = ( x+2 < width ) ? (x+2): (x-2);
      
      switch (bytes)
        {                
        case 1:
          dst[(x*bytes*2)] = src2[g2*bytes];
          sum  = cubic_spline_fit (0.5, src2[g1*bytes],
                                        src2[g2*bytes],
                                        src2[g3*bytes],
                                        src2[g4*bytes]);
          dst[(x*bytes*2) + bytes] = (guchar) CLAMP (sum, 0, 255);
          break;
        case 2:
          dst[(x*bytes*2)  ] = src2[g2*bytes];
          dst[(x*bytes*2)+1] = src2[g2*bytes+1];
          alphasum = cubic_spline_fit (0.5, src2[g1*bytes+1],
                                            src2[g2*bytes+1],
                                            src2[g3*bytes+1],
                                            src2[g4*bytes+1]);
          sum      = cubic_spline_fit (0.5, src2[g1*bytes]*src2[g1*bytes+1],
                                            src2[g2*bytes]*src2[g2*bytes+1],
                                            src2[g3*bytes]*src2[g3*bytes+1],
                                            src2[g4*bytes]*src2[g4*bytes+1]);
          if (alphasum > 0) 
            {
              sum /= alphasum;
              dst[(x*bytes*2) + bytes]     = (guchar) CLAMP (sum, 0, 255);
              dst[(x*bytes*2) + bytes + 1] = (guchar) CLAMP (alphasum, 0, 255);
            }
          else
            {            
              dst[(x*bytes*2) + bytes    ] = 0;
              dst[(x*bytes*2) + bytes + 1] = 0;              
            }        
          break;
        case 3:
          for (b=0 ; b < 3 ; b++)
            {
              dst[(x*bytes*2) + b] = src2[g2*bytes+b];
              sum = cubic_spline_fit (0.5, src2[g1*bytes+b],
                                           src2[g2*bytes+b],
                                           src2[g3*bytes+b],
                                           src2[g4*bytes+b]);
              dst[(x*bytes*2) + bytes + b] = (guchar) CLAMP (sum, 0, 255);
            }
          break;
        case 4:
          dst[(x*bytes*2)] = src2[g2*bytes];
          dst[(x*bytes*2) + 3] = src2[g2*bytes+3];
          alphasum = cubic_spline_fit (0.5, src2[g1*bytes+3],
                                            src2[g2*bytes+3],
                                            src2[g3*bytes+3],
                                            src2[g4*bytes+3]);
          if (alphasum > 0) 
            {
              dst[(x*bytes*2) + bytes + 3] = (guchar) CLAMP (alphasum, 0, 255);
              for (b=0 ; b < 3 ; b++)
                {
                  sum = cubic_spline_fit (0.5, src2[g1*bytes+b]*src2[g1*bytes+3],
                                               src2[g2*bytes+b]*src2[g2*bytes+3],
                                               src2[g3*bytes+b]*src2[g3*bytes+3],
                                               src2[g4*bytes+b]*src2[g4*bytes+3]);
                  sum /= alphasum;
                  dst[(x*bytes*2) + bytes + b] = (guchar) CLAMP (sum, 0, 255);
                }
            }
          else
            {
              dst[(x*bytes*2) + bytes + 0] = 0;
              dst[(x*bytes*2) + bytes + 1] = 0;
              dst[(x*bytes*2) + bytes + 2] = 0;
              dst[(x*bytes*2) + bytes + 3] = 0;
            }
          break;
        }
    }
}

static void
expand_y_cubic (guchar       *src[5],
                guchar       *dst,
                gint          width,
                gint          bytes)
{

  gint          x, b;
  gdouble       alphasum;
  gdouble       sum;
  const guchar *src0;
  const guchar *src1;
  const guchar *src2;
  const guchar *src3;
  const guchar *src4;

  src0 = src[0];
  src1 = src[1];
  src2 = src[2];
  src3 = src[3];
  src4 = src[4];
  
  for (x=0; x < width ; x++)
    {
      switch (bytes)
        {                
        case 1:
          dst[x*bytes        ] = src2[x*bytes];
          sum                  = cubic_spline_fit (0.5,
                                                   src1[x*bytes],
                                                   src2[x*bytes],
                                                   src3[x*bytes],
                                                   src4[x*bytes]);
          dst[(x+width)*bytes] = (guchar) CLAMP (sum, 0, 255);
          break;
        case 2:
          dst[x*bytes  ] = src2[x*bytes];
          dst[x*bytes+1] = src2[x*bytes+1];
          alphasum = cubic_spline_fit (0.5, src1[x*bytes+1],
                                            src2[x*bytes+1],
                                            src3[x*bytes+1],
                                            src4[x*bytes+1]);
          sum      = cubic_spline_fit (0.5, src1[x*bytes]*src1[x*bytes+1],
                                            src2[x*bytes]*src2[x*bytes+1],
                                            src3[x*bytes]*src3[x*bytes+1],
                                            src4[x*bytes]*src4[x*bytes+1]);
          if (alphasum > 0) 
            {
              sum /= alphasum;
              dst[(x+width)*bytes]     = (guchar) CLAMP (sum, 0, 255);
              dst[(x+width)*bytes + 1] = (guchar) CLAMP (alphasum, 0, 255);
            }
          else
            {            
              dst[(x+width)*bytes] = 0;
              dst[(x+width)*bytes + 1] = 0;              
            }        
          break;
        case 3:
          for (b=0 ; b < 3 ; b++)
            {
              dst[x*bytes         + b] = src2[x*bytes+b];
              sum                      = cubic_spline_fit (0.5,
                                                           src1[x*bytes+b],
                                                           src2[x*bytes+b],
                                                           src3[x*bytes+b],
                                                           src4[x*bytes+b]);
              dst[(x+width)*bytes + b] = (guchar) CLAMP (sum, 0, 255);
            }
          break;
        case 4:
          dst[x*bytes  ] = src2[x*bytes];
          dst[x*bytes+1] = src2[x*bytes+1];
          alphasum = cubic_spline_fit (0.5, src1[x*bytes+3],
                                            src2[x*bytes+3],
                                            src3[x*bytes+3],
                                            src4[x*bytes+3]);
          sum      = cubic_spline_fit (0.5, src1[x*bytes]*src1[x*bytes+3],
                                            src2[x*bytes]*src2[x*bytes+3],
                                            src3[x*bytes]*src3[x*bytes+3],
                                            src4[x*bytes]*src4[x*bytes+3]);
          if (alphasum > 0) 
            {
              dst[(x+width)*bytes + 3] = (guchar) CLAMP (alphasum, 0, 255);
              for (b=0 ; b < 3 ; b++)
                {
                  sum      = cubic_spline_fit (0.5, src1[x*bytes+b]*src1[x*bytes+3],
                                                    src2[x*bytes+b]*src2[x*bytes+3],
                                                    src3[x*bytes+b]*src3[x*bytes+3],
                                                    src4[x*bytes+b]*src4[x*bytes+3]);
                  sum /= alphasum;
                  dst[(x+width)*bytes + b] = (guchar) CLAMP (sum, 0, 255);
                }
            }
          else
            { 
              for (b=0 ; b<bytes ; b++)
                {
                  dst[(x+width)*bytes + b] = 0;
                }
            }        

          break;
        }
    }
}

/* 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;
}


/* Expand / Reduce */
static void
reduce_x_nearest (guchar            *src[5],
                  guchar            *dst,
                  gint               width,
                  gint               bytes)
{
  gint          x;
  guchar       *dstptr = dst;
  guchar       *srcptr = src[2];

  for (x=0; x < width ; x+=2)
    {
      memcpy(dstptr, srcptr, bytes);
      srcptr += bytes * 2;
      dstptr += bytes;
    }
}

static void
reduce_y_nearest (guchar            *src[5],
                  guchar            *dst,
                  gint               width,
                  gint               bytes)
{
  const guchar *src2;
  src2 = src[2];
  memcpy(dst, src2, width * bytes);
}

static void
expand_x_nearest (guchar           *src[5],
                  guchar           *dst,
                  gint              width,
                  gint              bytes)
{
  gint          x;
  const guchar *src2   = src[2];
  guchar       *dstptr1 = dst;
  guchar       *dstptr2 = dst + bytes;

  for (x=0; x < width; x++, dstptr1+=bytes*2, dstptr2+=bytes*2, src2+=bytes)
    {
      memcpy(dstptr1, src2, bytes);
      memcpy(dstptr2, src2, bytes);
    }
}

static void
expand_y_nearest (guchar       *src[5],
                  guchar       *dst,
                  gint          width,
                  gint          bytes)
{
  gint          x;
  const guchar *src2   = src[2];
  guchar       *dstptr1 = dst;
  guchar       *dstptr2 = dst + width * bytes;
 
  for (x=0; x < width*bytes ; x++)
    {
      *(dstptr1++)=*src2;
      *(dstptr2++)=*(src2++);
    }
}

static void
reduce_x_linear (guchar            *src[5],
                 guchar            *dst,
                 gint               width,
                 gint               bytes)
{

  gint          x, b;
  gint          g2,g3;
  gdouble       alphasum;
  gdouble       sum;
  const guchar *src2;
  guchar       *dstptr;

  dstptr = dst;

  src2 = src[2];

  for (x=0; x < width ; x+=2)
    {
      g2 =  x;
      g3 = ( x+1 < width ) ? (x+1): (x-1);
      switch (bytes)
        {                
          case 1:
            sum =  (src2[g2*bytes] + src2[g3*bytes]) / 2 ;
            *dstptr = (guchar) CLAMP (sum, 0, 255);
          break;
        case 2:
            alphasum =  (src2[g2*bytes+1] + src2[g3*bytes+1]);
            sum      =  (src2[g2*bytes] * src2[g2*bytes+1] +
                         src2[g3*bytes] * src2[g3*bytes+1]);
            if (alphasum > 0) 
              {
                sum /= alphasum;
                alphasum /= 2;
                dstptr[0] = (guchar) CLAMP (sum, 0, 255);
                dstptr[1] = (guchar) CLAMP (alphasum, 0, 255);
              }
            else 
              {
                dstptr[0] = 0;
                dstptr[1] = 0;
              }
          break;
        case 3:
          for (b=0 ; b < 3 ; b++)
            {
              sum =  (src2[g2*bytes+b] + src2[g3*bytes+b]) / 2 ;
              dstptr[b] = (guchar) CLAMP (sum, 0, 255);
            }
          break;
        case 4:
            alphasum =  (src2[g2*bytes+3] + src2[g3*bytes+3]);
            if (alphasum > 0) 
              {
                for (b=0 ; b < 3 ; b++)
                  {
                    sum = (src2[g2*bytes+b] * src2[g2*bytes+3] +
                           src2[g3*bytes+b] * src2[g3*bytes+3]);
                    sum /= alphasum;
                    dstptr[b] = (guchar) CLAMP (sum, 0, 255);
                  }
                alphasum /= 2;
                dstptr[3] = (guchar) CLAMP (alphasum, 0, 255);
              }
            else 
              {
                for (b=0 ; b < bytes ; b++)
                  {
                    dstptr[b] = 0;
                  }
              }
          break;
        }
      dstptr += bytes;
    }
}

static void
reduce_y_linear (guchar            *src[5],
                 guchar            *dst,
                 gint               width,
                 gint               bytes)
{

  gint          x, b;
  gdouble       alphasum;
  gdouble       sum;
  const guchar *src2;
  const guchar *src3;

  src2 = src[2];
  src3 = src[3];
  
  for (x=0; x < width ; x++)
    {
      switch (bytes)
        {                
        case 1:
          sum  = (src2[x*bytes] + src3[x*bytes]) / 2;
          dst[x*bytes] = (guchar) CLAMP (sum, 0, 255);
          break;
        case 2:
          alphasum  = (src2[x*bytes+1] + src3[x*bytes+1]) / 2;
          if (alphasum > 0) 
            {
               sum  = (src2[x*bytes] * src2[x*bytes+1] + 
                       src3[x*bytes] * src3[x*bytes+1]) / alphasum;
               dst[x*bytes] = (guchar) CLAMP (sum, 0, 255);            
               dst[x*bytes+1] = (guchar) CLAMP (alphasum, 0, 255);
            }
          else
            {
               dst[x*bytes]   = 0;            
               dst[x*bytes+1] = 0;
            }
          break;
        case 3:
          for (b=0 ; b < 3 ; b++)
            {
              sum  = (src2[x*bytes+b] + src3[x*bytes+b]) / 2;
              dst[x*bytes+b] = (guchar) CLAMP (sum, 0, 255);
            }
          break;
        case 4:
          alphasum  = (src2[x*bytes+1] + src3[x*bytes+1]);
          if (alphasum > 0) 
            {
              for (b=0 ; b < 3 ; b++)
                {
                  sum  = (src2[x*bytes+b] * src2[x*bytes+3] + 
                          src3[x*bytes+b] * src3[x*bytes+3]) / alphasum;
                  dst[x*bytes] = (guchar) CLAMP (sum, 0, 255);            
                }
              alphasum /= 2;
              dst[x*bytes+3] = (guchar) CLAMP (alphasum, 0, 255);
            }
          else
            {
              for (b=0 ; b < bytes ; b++)
                {
                  dst[x*bytes+b] = 0;
                }
            }
          break;
        }
    }
}

static void
reduce_x_lanczos2 (guchar            *src[5],
                   guchar            *dst,
                   gint               width,
                   gint               bytes)
{
  InterpolateKernel kernel  = {5, {-1,  9,  16, 9, -1}, 32}; /* Turkowski decimator */
  gint          x, b;
  gint          g0,g1,g2,g3,g4;
  gdouble       alphasum;
  gdouble       sum;
  const guchar *src2;
  guchar       *dstptr;
  
  dstptr = dst;

  src2 = src[2];

  for (x=0; x < width ; x+=2)
    {
      g0 = ( x-2 >= 0 ) ? (x-2): (x+2);
      g1 = ( x-1 >= 0 ) ? (x-1): (x+1);
      g2 =  x;
      g3 = ( x+1 < width ) ? (x+1): (x-1);
      g4 = ( x+2 < width ) ? (x+2): (x-2);
      switch (bytes)
        {                
        case 1:
          sum =  src2[g0*bytes] * kernel.q[0];
          sum += src2[g1*bytes] * kernel.q[1];
          sum += src2[g2*bytes] * kernel.q[2];
          sum += src2[g3*bytes] * kernel.q[3];
          sum += src2[g4*bytes] * kernel.q[4];
          sum /= kernel.weight;
          dstptr[0] = (guchar) CLAMP (sum, 0, 255);
          break;
        case 2:
          alphasum =  src2[g0*bytes+1] * kernel.q[0];
          alphasum += src2[g1*bytes+1] * kernel.q[1];
          alphasum += src2[g2*bytes+1] * kernel.q[2];
          alphasum += src2[g3*bytes+1] * kernel.q[3];
          alphasum += src2[g4*bytes+1] * kernel.q[4];
          sum =  src2[g0*bytes] * kernel.q[0];
          sum += src2[g1*bytes] * kernel.q[1];
          sum += src2[g2*bytes] * kernel.q[2];
          sum += src2[g3*bytes] * kernel.q[3];
          sum += src2[g4*bytes] * kernel.q[4];
          if (alphasum > 0) 
            {
              sum /= alphasum;
              alphasum /= kernel.weight;
              dstptr[0]     = (guchar) CLAMP (sum, 0, 255);
              dstptr[1] = (guchar) CLAMP (alphasum, 0, 255);
            }
          else
            {            
              dstptr[0] = 0;
              dstptr[1] = 0;              
            }        
          break;
        case 3:
          for (b=0 ; b < 3 ; b++)
            {
              sum =  src2[g0*bytes+b] * kernel.q[0];
              sum += src2[g1*bytes+b] * kernel.q[1];
              sum += src2[g2*bytes+b] * kernel.q[2];
              sum += src2[g3*bytes+b] * kernel.q[3];
              sum += src2[g4*bytes+b] * kernel.q[4];
              sum /= kernel.weight;
              dstptr[b] = (guchar) CLAMP (sum, 0, 255);
            }
          break;
        case 4:
          alphasum =  src2[g0*bytes+3] * kernel.q[0];
          alphasum += src2[g1*bytes+3] * kernel.q[1];
          alphasum += src2[g2*bytes+3] * kernel.q[2];
          alphasum += src2[g3*bytes+3] * kernel.q[3];
          alphasum += src2[g4*bytes+3] * kernel.q[4];
          if (alphasum > 0) 
            {
              for (b=0 ; b < 3 ; b++)
                {
                  sum =  src2[g0*bytes+b] * kernel.q[0];
                  sum += src2[g1*bytes+b] * kernel.q[1];
                  sum += src2[g2*bytes+b] * kernel.q[2];
                  sum += src2[g3*bytes+b] * kernel.q[3];
                  sum += src2[g4*bytes+b] * kernel.q[4];
                  sum /= alphasum;
                  dstptr[b] = (guchar) CLAMP (sum, 0, 255);
                }
              alphasum /= kernel.weight;
              dstptr[4] = (guchar) CLAMP (alphasum, 0, 255);
            }
          else
            { 
              for ( b=0 ; b < bytes ; b++)
                {
                  dstptr[b] = 0;
                }
            }        
          break;
        }
      dstptr += bytes;
    }
}

static void
reduce_y_lanczos2 (guchar            *src[5],
                   guchar            *dst,
                   gint               width,
                   gint               bytes)
{
  InterpolateKernel kernel  = {5, {-1,  9,  16, 9, -1}, 32}; /* Turkowski decimator */
  gint          x, b;
  gdouble       alphasum;
  gdouble       sum;
  const guchar *src0;
  const guchar *src1;
  const guchar *src2;
  const guchar *src3;
  const guchar *src4;

  src0 = src[0];
  src1 = src[1];
  src2 = src[2];
  src3 = src[3];
  src4 = src[4];
  
  for (x=0; x < width ; x++)
    {
      switch (bytes)
        {                
        case 1:
          sum  = src0[x*bytes] * kernel.q[0];
          sum += src1[x*bytes] * kernel.q[1];
          sum += src2[x*bytes] * kernel.q[2];
          sum += src3[x*bytes] * kernel.q[3];
          sum += src4[x*bytes] * kernel.q[4];
          sum /= kernel.weight;
          dst[x*bytes] = (guchar) CLAMP (sum, 0, 255);
          break;
        case 2:
          alphasum =  src0[x*bytes+1] * kernel.q[0];
          alphasum += src1[x*bytes+1] * kernel.q[1];
          alphasum += src2[x*bytes+1] * kernel.q[2];
          alphasum += src3[x*bytes+1] * kernel.q[3];
          alphasum += src4[x*bytes+1] * kernel.q[4];
          sum  = src0[x*bytes] * src0[x*bytes+1] * kernel.q[0];
          sum += src1[x*bytes] * src1[x*bytes+1] * kernel.q[1];
          sum += src2[x*bytes] * src2[x*bytes+1] * kernel.q[2];
          sum += src3[x*bytes] * src3[x*bytes+1] * kernel.q[3];
          sum += src4[x*bytes] * src4[x*bytes+1] * kernel.q[4];
          if (alphasum > 0) 
            {
              sum /= alphasum;
              alphasum /= kernel.weight;
              dst[x*bytes]    = (guchar) CLAMP (sum, 0, 255);
              dst[x*bytes+ 1] = (guchar) CLAMP (alphasum, 0, 255);
            }
          else
            {            
              dst[x*bytes] = 0;
              dst[x*bytes + 1] = 0;              
            }
          break;
        case 3:
          for (b=0 ; b < 3 ; b++)
            {
              sum  = src0[x*bytes+b] * kernel.q[0];
              sum += src1[x*bytes+b] * kernel.q[1];
              sum += src2[x*bytes+b] * kernel.q[2];
              sum += src3[x*bytes+b] * kernel.q[3];
              sum += src4[x*bytes+b] * kernel.q[4];
              sum /= kernel.weight;
              dst[x*bytes+b] = (guchar) CLAMP (sum, 0, 255);
            }
          break;
        case 4:
          alphasum =  src0[x*bytes+3] * kernel.q[0];
          alphasum += src1[x*bytes+3] * kernel.q[1];
          alphasum += src2[x*bytes+3] * kernel.q[2];
          alphasum += src3[x*bytes+3] * kernel.q[3];
          alphasum += src4[x*bytes+3] * kernel.q[4];
          if (alphasum > 0) 
            {
              for (b=0 ; b<3 ; b++)
                {
                  sum  = src0[x*bytes] * src0[x*bytes+3] * kernel.q[0];
                  sum += src1[x*bytes] * src1[x*bytes+3] * kernel.q[1];
                  sum += src2[x*bytes] * src2[x*bytes+3] * kernel.q[2];
                  sum += src3[x*bytes] * src3[x*bytes+3] * kernel.q[3];
                  sum += src4[x*bytes] * src4[x*bytes+3] * kernel.q[4];
                  sum /= alphasum;
                  dst[x*bytes+b]    = (guchar) CLAMP (sum, 0, 255);
                }
              alphasum /= kernel.weight;
              dst[x*bytes+ 3] = (guchar) CLAMP (alphasum, 0, 255);
            }
          else
            {            
              for ( b=0 ; b < bytes ; b++)
                {
                  dst[x*bytes+b] = 0;
                }
            }
          break;
        }
    }
}

static void
expand_x_linear (guchar           *src[5],
                 guchar           *dst,
                 gint              width,
                 gint              bytes)
{

  gint          x, b;
  gint          g2,g3;
  gdouble       alphasum;
  gdouble       sum;
  const guchar *src2;

  src2 = src[2];

  for (x=0; x < width ; x++)
    {
      g2 =  x;
      g3 = ( x+1 < width ) ? (x+1): (x-1);    
      switch (bytes)
        {                
        case 1:
          dst[(x*bytes*2)] = src2[g2*bytes];
          sum  = ( src2[g2*bytes] + src2[g3*bytes] ) / 2;
          dst[(x*bytes*2) + bytes] = (guchar) CLAMP (sum, 0, 255);
          break;
        case 2:
          dst[(x*bytes*2)]   = src2[g2*bytes];
          dst[(x*bytes*2)+1] = src2[g2*bytes+1];
          alphasum  = ( src2[g2*bytes+1] + src2[g3*bytes+1] );
          sum  = ( src2[g2*bytes] * src2[g2*bytes+1] + 
                   src2[g3*bytes] * src2[g2*bytes+1] );
          if ( alphasum > 0 )
            {
              sum/=alphasum;
              dst[(x*bytes*2) + bytes] = (guchar) CLAMP (sum, 0, 255);
              alphasum /= 2;
              dst[(x*bytes*2) + bytes+1] = (guchar) CLAMP (alphasum, 0, 255);
            }
          else
            {            
              dst[(x*bytes*2) + bytes  ] = 0;
              dst[(x*bytes*2) + bytes+1] = 0;           
            }
          break;
        case 3:
          for (b=0 ; b < 3 ; b++)
            {
              sum  = ( src2[g2*bytes+b] + src2[g3*bytes+b] ) / 2;
              dst[(x*bytes*2)         + b] = src2[g2*bytes+b];
              dst[(x*bytes*2) + bytes + b] = (guchar) CLAMP (sum, 0, 255);
            }
          break;
        case 4:
          dst[(x*bytes*2)  ] = src2[g2*bytes];
          dst[(x*bytes*2)+3] = src2[g2*bytes+3];
          alphasum  = ( src2[g2*bytes+3] + src2[g3*bytes+3] );
          if ( alphasum > 0 )
            {
              for (b=0 ; b < 3 ; b++)
                {
                  sum  = ( src2[g2*bytes+b] + src2[g3*bytes+b] ) / 2;
                  sum/=alphasum;
                  dst[(x*bytes*2) + bytes + b] = (guchar) CLAMP (sum, 0, 255);
                }
              alphasum /= 2;
              dst[(x*bytes*2) + bytes + 3] = (guchar) CLAMP (alphasum, 0, 255);
            }
          else
            {
              for ( b=0 ; b < bytes ; b++)
                {
                  dst[x*bytes+b] = 0;
                }
            }
          break;
        }
    }
}

static void
expand_x_lanczos2 (guchar           *src[5],
                   guchar           *dst,
                   gint              width,
                   gint              bytes)
{
  InterpolateKernel kernel   = {4, {-1,  9,  9, -1,  0}, 16}; /*Lanczos2 phase (1/2) */ 
  gint          x, b;
  gint          g0,g1,g2,g3,g4;
  gdouble       alphasum;
  gdouble       sum;
  const guchar *src2;

  src2 = src[2];

  for (x=0; x < width ; x++)
    {
      g0 = ( x-2 >= 0 ) ? (x-2): (x+2);
      g1 = ( x-1 >= 0 ) ? (x-1): (x+1);
      g2 =  x;
      g3 = ( x+1 < width ) ? (x+1): (x-1);
      g4 = ( x+2 < width ) ? (x+2): (x-2);
      
      switch (bytes)
        {                
        case 1:
          sum  = src2[g1*bytes] * kernel.q[0];
          sum += src2[g2*bytes] * kernel.q[1];
          sum += src2[g3*bytes] * kernel.q[2];
          sum += src2[g4*bytes] * kernel.q[3];
          sum /= kernel.weight;
          dst[(x*bytes*2)        ] = src2[g2*bytes];
          dst[(x*bytes*2) + bytes] = (guchar) CLAMP (sum, 0, 255);
          break;
        case 2:
          dst[(x*bytes*2)]   = src2[g2*bytes];
          dst[(x*bytes*2)+1] = src2[g2*bytes+1];
          alphasum  = src2[g1*bytes+1] * kernel.q[0];
          alphasum += src2[g2*bytes+1] * kernel.q[1];
          alphasum += src2[g3*bytes+1] * kernel.q[2];
          alphasum += src2[g4*bytes+1] * kernel.q[3];          
          sum  = src2[g1*bytes] * src2[g1*bytes+1] * kernel.q[0];
          sum += src2[g2*bytes] * src2[g2*bytes+1] * kernel.q[1];
          sum += src2[g3*bytes] * src2[g3*bytes+1] * kernel.q[2];
          sum += src2[g4*bytes] * src2[g4*bytes+1] * kernel.q[3];          
          if ( alphasum > 0 )
            {              
              sum/=alphasum;
              dst[(x*bytes*2) + bytes] = (guchar) CLAMP (sum, 0, 255);
              alphasum /= 2;
              dst[(x*bytes*2) + bytes+1] = (guchar) CLAMP (alphasum, 0, 255);
            }
          else
            {            
              dst[(x*bytes*2) + bytes  ] = 0;
              dst[(x*bytes*2) + bytes+1] = 0;           
            }
          break;
        case 3:
          for (b=0 ; b < 3 ; b++)
            {
              sum  = src2[g1*bytes+b] * kernel.q[0];
              sum += src2[g2*bytes+b] * kernel.q[1];
              sum += src2[g3*bytes+b] * kernel.q[2];
              sum += src2[g4*bytes+b] * kernel.q[3];
              sum /= kernel.weight;
              dst[(x*bytes*2)         + b] = src2[g2*bytes+b];
              dst[(x*bytes*2) + bytes + b] = (guchar) CLAMP (sum, 0, 255);
            }
          break;
        case 4:
          dst[(x*bytes*2)]   = src2[g2*bytes];
          dst[(x*bytes*2)+1] = src2[g2*bytes+1];
          alphasum  = src2[g1*bytes+1] * kernel.q[0];
          alphasum += src2[g2*bytes+1] * kernel.q[1];
          alphasum += src2[g3*bytes+1] * kernel.q[2];
          alphasum += src2[g4*bytes+1] * kernel.q[3];          
          if ( alphasum > 0 )
            {
              for (b=0 ; b<3 ; b++)
                {
                  sum  = src2[g1*bytes+b] * src2[g1*bytes+1] * kernel.q[0];
                  sum += src2[g2*bytes+b] * src2[g2*bytes+1] * kernel.q[1];
                  sum += src2[g3*bytes+b] * src2[g3*bytes+1] * kernel.q[2];
                  sum += src2[g4*bytes+b] * src2[g4*bytes+1] * kernel.q[3];          
                  sum/=alphasum;
                  dst[(x*bytes*2) + bytes] = (guchar) CLAMP (sum, 0, 255);
                }
              alphasum /= kernel.weight;
              dst[(x*bytes*2) + bytes+1] = (guchar) CLAMP (alphasum, 0, 255);
            }
          else
            {            
              for (b=0 ; b<bytes ; b++)
                {
                  dst[(x*bytes*2) + b] = 0;
                }
            }
          break;
        }
    }
}


static void
expand_y_linear (guchar       *src[5],
                 guchar       *dst,
                 gint          width,
                 gint          bytes)
{
  gint          x, b;
  gdouble       alphasum;
  gdouble       sum;
  const guchar *src2;
  const guchar *src3;

  src2 = src[2];
  src3 = src[3];
  
  for (x=0; x < width ; x++)
    {
      switch (bytes)
        {                
        case 1:
          sum  = (src2[x*bytes] + src3[x*bytes]) / 2;
          dst[x*bytes] = src2[x*bytes];
          dst[(x+width)*bytes] = (guchar) CLAMP (sum, 0, 255);
          break;
        case 2:
          dst[x*bytes]   = src2[x*bytes];
          dst[x*bytes+1] = src2[x*bytes+1];
          alphasum  = (src2[x*bytes+1] + src3[x*bytes+1]);
          if (alphasum > 0 )
            {
              sum  = (src2[x*bytes] * src2[x*bytes+1] + 
                      src3[x*bytes] * src3[x*bytes+1]) / alphasum;
              dst[(x+width)*bytes] = (guchar) CLAMP (sum, 0, 255);
              alphasum/=2;
              dst[(x+width)*bytes+1]= (guchar) CLAMP (alphasum, 0, 255);
            }
          else
            {
              dst[(x+width)*bytes] = 0;
              dst[(x+width)*bytes+1]= 0;
            }            
          break;
        case 3:
          for (b=0 ; b < 3 ; b++)
            {
              sum  = (src2[x*bytes+b] + src3[x*bytes+b]) / 2;
              dst[ x       *bytes   + b] = src2[x*bytes+b];
              dst[(x+width)*bytes   + b] = (guchar) CLAMP (sum, 0, 255);
            }
          break;
        case 4:
          dst[x*bytes]   = src2[x*bytes];
          dst[x*bytes+3] = src2[x*bytes+3];
          alphasum  = (src2[x*bytes+3] + src3[x*bytes+3]);
          if (alphasum > 0 )
            {
              for (b=0 ; b<3 ; b++)
                {
                  sum  = (src2[x*bytes+b] * src2[x*bytes+3] + 
                          src3[x*bytes+b] * src3[x*bytes+3]) / alphasum;
                  dst[(x+width)*bytes+b] = (guchar) CLAMP (sum, 0, 255);
                }
              alphasum /= 2;
              dst[(x+width)*bytes+3]= (guchar) CLAMP (alphasum, 0, 255);
            }
          else
            {
              for (b=0 ; b<bytes ; b++)
                {
                  dst[(x+width)*bytes+b] = 0;
                }
            }            
          break;
        }
    }
}

static void
expand_y_lanczos2 (guchar       *src[5],
                   guchar       *dst,
                   gint          width,
                   gint          bytes)
{
  InterpolateKernel kernel   = {4, {-1,  9,  9, -1,  0}, 16}; /*Lanczos2 phase (1/2) */ 
  gint          x, b;
  gdouble       alphasum;
  gdouble       sum;
  const guchar *src0;
  const guchar *src1;
  const guchar *src2;
  const guchar *src3;
  const guchar *src4;

  src0 = src[0];
  src1 = src[1];
  src2 = src[2];
  src3 = src[3];
  src4 = src[4];
  
  for (x=0; x < width ; x++)
    {
      switch (bytes)
        {                
        case 1:
          sum  = src1[x*bytes] * kernel.q[0];
          sum += src2[x*bytes] * kernel.q[1];
          sum += src3[x*bytes] * kernel.q[2];
          sum += src4[x*bytes] * kernel.q[3];
          sum /= kernel.weight;              
          dst[x*bytes        ] = src2[x*bytes];
          dst[(x+width)*bytes] = (guchar) CLAMP (sum, 0, 255);
          break;
        case 2:
          dst[x*bytes]   = src2[x*bytes];
          dst[x*bytes+1] = src2[x*bytes+1];
          alphasum  = src1[x*bytes+1] * kernel.q[0];
          alphasum += src2[x*bytes+1] * kernel.q[1];
          alphasum += src3[x*bytes+1] * kernel.q[2];
          alphasum += src4[x*bytes+1] * kernel.q[3];
          if (alphasum > 0 )
            {
              sum  = src1[x*bytes] * src1[x*bytes+1] * kernel.q[0];
              sum += src2[x*bytes] * src2[x*bytes+1] * kernel.q[1];
              sum += src3[x*bytes] * src3[x*bytes+1] * kernel.q[2];
              sum += src4[x*bytes] * src4[x*bytes+1] * kernel.q[3];
              sum /= alphasum;
              alphasum /= kernel.weight;
              dst[(x+width)*bytes] = (guchar) CLAMP (sum, 0, 255);
              dst[(x+width)*bytes+1] = (guchar) CLAMP (alphasum, 0, 255);
            }
          else
            {
              dst[(x+width)*bytes]   = 0;
              dst[(x+width)*bytes+1] = 0;
            }             
          break;
        case 3:
          for (b=0 ; b < 3 ; b++)
            {
                  sum  = src1[x*bytes+b] * kernel.q[0];
                  sum += src2[x*bytes+b] * kernel.q[1];
                  sum += src3[x*bytes+b] * kernel.q[2];
                  sum += src4[x*bytes+b] * kernel.q[3];
                  sum /= kernel.weight;              
                  dst[x*bytes           + b] = src2[x*bytes+b];
                  dst[(x+width)*bytes   + b] = (guchar) CLAMP (sum, 0, 255);
            }
          break;
        case 4:
          dst[x*bytes+3] = src2[x*bytes+3];
          alphasum  = src1[x*bytes+3] * kernel.q[0];
          alphasum += src2[x*bytes+3] * kernel.q[1];
          alphasum += src3[x*bytes+3] * kernel.q[2];
          alphasum += src4[x*bytes+3] * kernel.q[3];
          if (alphasum > 0 )
            {
              for (b=0 ; b<3 ;b++)
                {
                  sum  = src1[x*bytes+b] * src1[x*bytes+3] * kernel.q[0];
                  sum += src2[x*bytes+b] * src2[x*bytes+3] * kernel.q[1];
                  sum += src3[x*bytes+b] * src3[x*bytes+3] * kernel.q[2];
                  sum += src4[x*bytes+b] * src4[x*bytes+3] * kernel.q[3];
                  sum /= alphasum;
                  dst[(x+width)*bytes] = (guchar) CLAMP (sum, 0, 255);
                }
              alphasum /= kernel.weight;
              dst[(x+width)*bytes+3] = (guchar) CLAMP (alphasum, 0, 255);
            }
          else
            {
              for (b=0 ; b<bytes ; b++)
                {
                  dst[(x+width)*bytes+b] = 0;
                }
            }
          break;
        }
    }
}

/* Scale functions */
static void
scale_nearest_pr     (PixelRegion  *dstPR,
                      PixelRegion  *srcPR)
{
  const gint bytes = dstPR->bytes;
  gint       x, y, b;
  gint       x2, x3;
  gint       y2, y3;


  gdouble    scalex, scaley;
  gdouble    dx, dy;
 
  guchar *dst;
  guchar *src2;
 
  scaley = (gdouble) dstPR->h / (gdouble) srcPR->h;
  scalex = (gdouble) dstPR->w / (gdouble) srcPR->w;

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

  /* Interpolate */  
  src2 = g_new0 (guchar, srcPR->w * bytes);
  dst  = g_new0 (guchar, dstPR->w * bytes);

  for (y = 0; y < dstPR->h; y++)
    {
      dy = y / scaley;    
      y2 = (gint) floor (dy);
      y3 = (gint) ceil (dy);
      y3 = ( y3 >= srcPR->h) ? y2 : y3;
      dy =  dy - y2;
      y2 = ( dy < 0.5) ? y2 : y3;
 
      pixel_region_get_row (srcPR, 0, y2, srcPR->w, src2, 1);
 
      for (x = 0; x < dstPR->w; x++)
        {
          dx = (x / scalex);        
          x2 = (gint)floor(dx);
          x3 = (gint)ceil(dx);
          x3 = ( x3 >= srcPR->w) ? x2 : x3;
          dx =  dx - x2;
          x2 = ( dy < 0.5) ? x2 : x3;
          
          for (b=0 ; b < bytes ; b++)
            {
              dst[x*bytes+b] = src2[x2*bytes+b];
            }
        }

      pixel_region_set_row (dstPR, 0, y, dstPR->w, dst);
    }

  g_free(src2);
  g_free(dst);
}
                                 
static void
scale_nearest      (PixelRegion  *dstPR,
                    const guchar *src,
                    gint          source_w,
                    gint          source_h,
                    gint          bytes)
{
  gint       x, y, b;
  gint       x2, x3;
  gint       y2, y3;

  gdouble    scalex, scaley;
  gdouble    dx, dy;
 
  guchar *dst;
  guchar *src2;
 
  scaley = (gdouble) dstPR->h / (gdouble) source_h;
  scalex = (gdouble) dstPR->w / (gdouble) source_w;

  /* Interpolate */  
  src2 = g_new0 (guchar, source_w * bytes); 
  dst  = g_new0 (guchar, dstPR->w * bytes);

  for (y = 0; y < dstPR->h; y++)
    {
      dy = y / scaley;    
      y2 = (gint) floor (dy);
      y3 = (gint) ceil (dy);
      y3 = ( y3 >= source_h) ? y2 : y3;      
      dy =  dy - y2;
      y2 = ( dy < 0.5) ? y2 : y3;
      
      memcpy(src2, src+y2*source_w*bytes,source_w*bytes);

      for (x = 0; x < dstPR->w; x++)
        {
          dx = (x / scalex);        
          x2 = (gint)floor(dx);
          x3 = (gint)ceil(dx);
          x3 = ( x3 >= source_w) ? x2 : x3;      
          dx =  dx - x2;
          x2 = ( dy < 0.5) ? x2 : x3;         
          for (b=0 ; b < bytes ; b++)
            {
              dst[x*bytes+b] = src2[x2*bytes+b];
            }
        }

      pixel_region_set_row (dstPR, 0, y, dstPR->w, dst);
    }

  g_free(src2);
  g_free(dst);
}                                 

static void
scale_bilinear (PixelRegion  *dstPR,
                const guchar *src,
                gint          source_w,
                gint          source_h,
                gint          bytes)
{
  const gdouble  scaley = (gdouble) dstPR->h / (gdouble) source_h;
  const gdouble  scalex = (gdouble) dstPR->w / (gdouble) source_w;

  gint     x, y, b;
  gint     x0, x1, y0, y1;
  gdouble  xfrac, yfrac;
  gdouble  s00, s01, s10, s11;
  gdouble  alphasum;
  gdouble  sum;
  guchar  *dst;

  dst = g_new (guchar, dstPR->w * dstPR->h * bytes);

  for (y = 0; y < dstPR->h; y++)
    {
      yfrac = ( y / scaley );
      y0 = (gint) floor (yfrac);
      y1 = (gint) ceil (yfrac);
      y1 = ( y1 >= source_h) ? y0 : y1;      
      yfrac =  yfrac - y0;

      for (x = 0; x < dstPR->w; x++)
        {
           xfrac = (x / scalex);
           x0 = (gint) floor (xfrac);
           x1 = (gint) ceil (xfrac);
           x1 = ( x1 >= source_w) ? x0 : x1;           
           xfrac =  xfrac - x0;           
           switch (bytes)
             {                
             case 1:
               s00 = src[(x0 + y0 * source_w) * bytes];
               s10 = src[(x1 + y0 * source_w) * bytes];
               s01 = src[(x0 + y1 * source_w) * bytes];
               s11 = src[(x1 + y1 * source_w) * bytes];
               sum =  (gint)((1 - yfrac) * ((1 - xfrac) * s00 + xfrac * s01)
                                  + yfrac  * ((1 - xfrac) * s10 + xfrac * s11));
               dst[x * bytes] = sum;
               break;
             case 2:
               s00 = src[(x0 + y0 * source_w) * bytes+1];
               s10 = src[(x1 + y0 * source_w) * bytes+1];
               s01 = src[(x0 + y1 * source_w) * bytes+1];
               s11 = src[(x1 + y1 * source_w) * bytes+1];
               alphasum =  (gint)((1 - yfrac) * ((1 - xfrac) * s00 + xfrac * s01)
                                       + yfrac  * ((1 - xfrac) * s10 + xfrac * s11));
               if (alphasum > 0)
                 {
                   s00 *= src[(x0 + y0 * source_w) * bytes];
                   s10 *= src[(x1 + y0 * source_w) * bytes];
                   s01 *= src[(x0 + y1 * source_w) * bytes];
                   s11 *= src[(x1 + y1 * source_w) * bytes];
                   sum =  (gint)((1 - yfrac) * ((1 - xfrac) * s00 + xfrac * s01)
                                      + yfrac  * ((1 - xfrac) * s10 + xfrac * s11));
                   sum /= alphasum;
                   dst[x * bytes]   = (guchar) CLAMP (sum, 0, 255);
                   dst[x * bytes+1] = (guchar) CLAMP (alphasum, 0, 255);
                 }
               else
                {
                   dst[x * bytes]   = 0;
                   dst[x * bytes+1] = 0;                
                }
               break;
             case 3:
               for (b = 0 ; b < bytes ; b++)
                 {
                   s00 = src[(x0 + y0 * source_w) * bytes + b];
                   s10 = src[(x1 + y0 * source_w) * bytes + b];
                   s01 = src[(x0 + y1 * source_w) * bytes + b];
                   s11 = src[(x1 + y1 * source_w) * bytes + b];
                   sum =  (gint)((1 - yfrac) * ((1 - xfrac) * s00 + xfrac * s01)
                                  + yfrac  * ((1 - xfrac) * s10 + xfrac * s11));
                   dst[x * bytes + b] = sum;
                 }
               break;
             case 4:
               s00 = src[(x0 + y0 * source_w) * bytes+3];
               s10 = src[(x1 + y0 * source_w) * bytes+3];
               s01 = src[(x0 + y1 * source_w) * bytes+3];
               s11 = src[(x1 + y1 * source_w) * bytes+3];
               alphasum =  (gint)((1 - yfrac) * ((1 - xfrac) * s00 + xfrac * s01)
                                       + yfrac  * ((1 - xfrac) * s10 + xfrac * s11));
               if (alphasum > 0)
                 {
                   for (b = 0 ; b < bytes ; b++)
                     {
                       s00 = src[(x0 + y0 * source_w) * bytes + b] * src[(x0 + y0 * source_w) * bytes+3];
                       s10 = src[(x1 + y0 * source_w) * bytes + b] * src[(x1 + y0 * source_w) * bytes+3];
                       s01 = src[(x0 + y1 * source_w) * bytes + b] * src[(x0 + y1 * source_w) * bytes+3];
                       s11 = src[(x1 + y1 * source_w) * bytes + b] * src[(x1 + y1 * source_w) * bytes+3];
                       sum =  (gint)((1 - yfrac) * ((1 - xfrac) * s00 + xfrac * s01)
                                        + yfrac  * ((1 - xfrac) * s10 + xfrac * s11));
                       sum /= alphasum;
                       dst[x * bytes+b]   = (guchar) CLAMP (sum, 0, 255);
                     }
                    dst[x * bytes+3] = (guchar) CLAMP (alphasum, 0, 255);
                 }
               else
                {
                   for (b=0 ; b<bytes ; b++) 
                     {
                       dst[x * bytes+b]   = 0;
                     }
                }
               break;
             }
        }
        pixel_region_set_row (dstPR, 0, y, dstPR->w, dst);
    }
}

static void
scale_bilinear_pr (PixelRegion *dstPR,
                   PixelRegion *srcPR)
{
  const gint bytes = dstPR->bytes;
  gint       x, y, b;
  gint       x0, x1, y0, y1;

  gdouble    scale;
  gdouble    xfrac, yfrac;
  gdouble    s00, s01, s10, s11;

  gdouble       alphasum;
  gdouble       sum;

  guchar    *dst;
  guchar    *src0;
  guchar    *src1;

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

  scale = (gdouble) dstPR->h / (gdouble) srcPR->h;

  /* Interpolate */
  src0 = g_new0 (guchar, srcPR->w * bytes);
  src1 = g_new0 (guchar, srcPR->w * bytes);
  dst  = g_new0 (guchar, dstPR->w * bytes);

  for (y = 0; y < dstPR->h; y++)
    {
      yfrac = y / scale;
      y0 = (gint) floor (yfrac);
      y1 = (gint) ceil (yfrac);
      y1 = ( y1 >= srcPR->h) ? y0 : y1;      
      yfrac =  yfrac - y0;

      pixel_region_get_row (srcPR, 0, y0, srcPR->w, src0, 1);
      pixel_region_get_row (srcPR, 0, y1, srcPR->w, src1, 1);

      for (x = 0; x < dstPR->w; x++)
        {
          xfrac = (x / scale);
          x0 = (gint)floor(xfrac);
          x1 = (gint)ceil(xfrac);
          x1 = ( x1 >= srcPR->w) ? x0 : x1; 
          xfrac =  xfrac - x0;
          switch (bytes)
            {                
             case 1:
               s00 = src0[x0 * bytes];
               s01 = src0[x1 * bytes];
               s10 = src1[x0 * bytes];
               s11 = src1[x1 * bytes];
               sum =  (1 - yfrac) * ( (1 - xfrac) * s00 + xfrac * s01 )
                        +  yfrac  * ( (1 - xfrac) * s10 + xfrac * s11 );
               dst[x * bytes] = (guchar) sum;
               break;
             case 2:
               s00 = src0[x0 * bytes + 1];
               s01 = src0[x1 * bytes + 1];
               s10 = src1[x0 * bytes + 1];
               s11 = src1[x1 * bytes + 1];
               alphasum =  (gint)((1 - yfrac) * ((1 - xfrac) * s00 + xfrac * s01)
                                       + yfrac  * ((1 - xfrac) * s10 + xfrac * s11));
               if (alphasum > 0)
                 {
                   s00 *= src0[x0 * bytes];
                   s10 *= src0[x1 * bytes];
                   s01 *= src1[x0 * bytes];
                   s11 *= src1[x1 * bytes];
                   sum =  (gint)((1 - yfrac) * ((1 - xfrac) * s00 + xfrac * s01)
                                      + yfrac  * ((1 - xfrac) * s10 + xfrac * s11));
                   sum /= alphasum;
                   dst[x * bytes]   = (guchar) CLAMP (sum, 0, 255);
                   dst[x * bytes+1] = (guchar) CLAMP (alphasum, 0, 255);
                 }
               else
                {
                   dst[x * bytes]   = 0;
                   dst[x * bytes+1] = 0;                
                }
               break;
             case 3:
               for (b = 0 ; b < bytes ; b++)
                 {
                   s00 = src0[x0 * bytes + b];
                   s01 = src0[x1 * bytes + b];
                   s10 = src1[x0 * bytes + b];
                   s11 = src1[x1 * bytes + b];
                   sum =  (gint)((1 - yfrac) * ((1 - xfrac) * s00 + xfrac * s01)
                                  + yfrac  * ((1 - xfrac) * s10 + xfrac * s11));
                   dst[x * bytes + b] = sum;
                 }
               break;
             case 4:
               s00 = src0[x0 * bytes + 3];
               s01 = src0[x1 * bytes + 3];
               s10 = src1[x0 * bytes + 3];
               s11 = src1[x1 * bytes + 3];
               alphasum =  (gint)((1 - yfrac) * ((1 - xfrac) * s00 + xfrac * s01)
                                       + yfrac  * ((1 - xfrac) * s10 + xfrac * s11));
               if (alphasum > 0)
                 {
                   for (b = 0 ; b < bytes ; b++)
                     {
                       s00 = src0[x0 * bytes + b] * src0[x0 * bytes+3];
                       s10 = src0[x1 * bytes + b] * src0[x1 * bytes+3];
                       s01 = src1[x0 * bytes + b] * src1[x0 * bytes+3];
                       s11 = src1[x1 * bytes + b] * src1[x1 * bytes+3];
                       sum =  (gint)((1 - yfrac) * ((1 - xfrac) * s00 + xfrac * s01)
                                        + yfrac  * ((1 - xfrac) * s10 + xfrac * s11));
                       sum /= alphasum;
                       dst[x * bytes+b]   = (guchar) CLAMP (sum, 0, 255);
                     }
                    dst[x * bytes+3] = (guchar) CLAMP (alphasum, 0, 255);
                 }
               else
                {
                   for (b=0 ; b<bytes ; b++) 
                     {
                       dst[x * bytes+b]   = 0;
                     }
                }
               break;
           }
        }

      pixel_region_set_row (dstPR, 0, y, dstPR->w, dst);
    }

  g_free(src0);
  g_free(src1);
  g_free(dst);
}

static void
scale_cubic_pr     (PixelRegion  *dstPR,
                    PixelRegion  *srcPR)
{
  const gint bytes = dstPR->bytes;
  gint       x, y, b;
  gint       x1, x2, x3, x4;
  gint       y1, y2, y3, y4;


  gdouble    scalex, scaley;
  gdouble    dx, dy;
  gdouble    p0, p1, p2, p3;

  gdouble       alphasum;
  gdouble       sum;

  guchar *dst;
  guchar *src0;
  guchar *src1;
  guchar *src2;
  guchar *src3;
  guchar *src4;

  scaley = (gdouble) dstPR->h / (gdouble) srcPR->h;
  scalex = (gdouble) dstPR->w / (gdouble) srcPR->w;

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

  /* Interpolate */  
  src0 = g_new0 (guchar, srcPR->w * bytes * 5);
  src1 = src0 + srcPR->w * bytes;
  src2 = src1 + srcPR->w * bytes;
  src3 = src2 + srcPR->w * bytes;
  src4 = src3 + srcPR->w * bytes; 
  
  dst  = g_new0 (guchar, dstPR->w * bytes);

  for (y = 0; y < dstPR->h; y++)
    {
      dy = (y / scaley);
      y2 = (gint) floor (dy);     
      y1 = ( y2-1 >= 0 )       ? y2-1: y2+1;
      y3 = (gint) ceil (dy);
      y3 = ( y3 < srcPR->h)    ? y3  : y2;
      y4 = ( y3+1 < srcPR->h ) ? y3+1: y3-1;
      dy =  dy - y2;

      pixel_region_get_row (srcPR, 0, y1, srcPR->w, src1, 1);
      pixel_region_get_row (srcPR, 0, y2, srcPR->w, src2, 1);
      pixel_region_get_row (srcPR, 0, y3, srcPR->w, src3, 1);
      pixel_region_get_row (srcPR, 0, y4, srcPR->w, src4, 1);

      for (x = 0; x < dstPR->w; x++)
        {
          dx = (x / scalex);        
          x2 = (gint)floor(dx);
          x1 = ( x2-1 >= 0 )       ? x2-1: x2+1;
          x3 = (gint)ceil(dx);
          x3 = ( x3 < srcPR->w)    ? x3  : x2;
          x4 = ( x3+1 < srcPR->w ) ? x3+1: x3-1;
          dx =  dx - x2;
          
          switch (bytes)
            {                
             case 1:
               p0  = cubic_spline_fit (dx, src1[x1*bytes],
                                           src1[x2*bytes],
                                           src1[x3*bytes],
                                           src1[x4*bytes]);
               p1  = cubic_spline_fit (dx, src2[x1*bytes],
                                           src2[x2*bytes],
                                           src2[x3*bytes],
                                           src2[x4*bytes]);
               p2  = cubic_spline_fit (dx, src3[x1*bytes],
                                           src3[x2*bytes],
                                           src3[x3*bytes],
                                           src3[x4*bytes]);
               p3  = cubic_spline_fit (dx, src4[x1*bytes],
                                           src4[x2*bytes],
                                           src4[x3*bytes],
                                           src4[x4*bytes]);
               sum = cubic_spline_fit (dy, p0, p1, p2, p3);                                          
               dst[x*bytes] = (guchar) CLAMP (sum, 0, 255);
               break;
             case 2:
               p0  = cubic_spline_fit (dx, src1[x1*bytes+1],
                                           src1[x2*bytes+1],
                                           src1[x3*bytes+1],
                                           src1[x4*bytes+1]);
               p1  = cubic_spline_fit (dx, src2[x1*bytes+1],
                                           src2[x2*bytes+1],
                                           src2[x3*bytes+1],
                                           src2[x4*bytes+1]);
               p2  = cubic_spline_fit (dx, src3[x1*bytes+1],
                                           src3[x2*bytes+1],
                                           src3[x3*bytes+1],
                                           src3[x4*bytes+1]);
               p3  = cubic_spline_fit (dx, src4[x1*bytes+1],
                                           src4[x2*bytes+1],
                                           src4[x3*bytes+1],
                                           src4[x4*bytes+1]);
               alphasum = cubic_spline_fit (dy, p0, p1, p2, p3);                                          
               if (alphasum > 0)
                 {
                   p0  = cubic_spline_fit (dx, src1[x1*bytes] * src1[x1*bytes+1],
                                               src1[x2*bytes] * src1[x2*bytes+1],
                                               src1[x3*bytes] * src1[x3*bytes+1],
                                               src1[x4*bytes] * src1[x4*bytes+1]);
                   p1  = cubic_spline_fit (dx, src2[x1*bytes] * src2[x1*bytes+1],
                                               src2[x2*bytes] * src2[x2*bytes+1],
                                               src2[x3*bytes] * src2[x3*bytes+1],
                                               src2[x4*bytes] * src2[x4*bytes+1]);
                   p2  = cubic_spline_fit (dx, src3[x1*bytes] * src3[x1*bytes+1],
                                               src3[x2*bytes] * src3[x2*bytes+1],
                                               src3[x3*bytes] * src3[x3*bytes+1],
                                               src3[x4*bytes] * src3[x4*bytes+1]);
                   p3  = cubic_spline_fit (dx, src4[x1*bytes] * src4[x1*bytes+1],
                                               src4[x2*bytes] * src4[x2*bytes+1],
                                               src4[x3*bytes] * src4[x3*bytes+1],
                                               src4[x4*bytes] * src4[x4*bytes+1]);
                   sum = cubic_spline_fit (dy, p0, p1, p2, p3);                                          
                   sum /= alphasum;
                   dst[x * bytes]   = (guchar) CLAMP (sum, 0, 255);
                   dst[x * bytes+1] = (guchar) CLAMP (alphasum, 0, 255);
                 }
               else
                {
                   dst[x * bytes]   = 0;
                   dst[x * bytes+1] = 0;                
                }
               break;
             case 3:
               for (b = 0 ; b < bytes ; b++)
                 {
                   p0  = cubic_spline_fit (dx, src1[x1*bytes+b],
                                               src1[x2*bytes+b],
                                               src1[x3*bytes+b],
                                               src1[x4*bytes+b]);
                   p1  = cubic_spline_fit (dx, src2[x1*bytes+b],
                                               src2[x2*bytes+b],
                                               src2[x3*bytes+b],
                                               src2[x4*bytes+b]);
                   p2  = cubic_spline_fit (dx, src3[x1*bytes+b],
                                               src3[x2*bytes+b],
                                               src3[x3*bytes+b],
                                               src3[x4*bytes+b]);
                   p3  = cubic_spline_fit (dx, src4[x1*bytes+b],
                                               src4[x2*bytes+b],
                                               src4[x3*bytes+b],
                                               src4[x4*bytes+b]);
                   sum = cubic_spline_fit (dy, p0, p1, p2, p3);                                          
                   dst[x*bytes+b] = (guchar) CLAMP (sum, 0, 255);
                 }
               break;
             case 4:
               p0  = cubic_spline_fit (dx, src1[x1*bytes+3],
                                           src1[x2*bytes+3],
                                           src1[x3*bytes+3],
                                           src1[x4*bytes+3]);
               p1  = cubic_spline_fit (dx, src2[x1*bytes+3],
                                           src2[x2*bytes+3],
                                           src2[x3*bytes+3],
                                           src2[x4*bytes+3]);
               p2  = cubic_spline_fit (dx, src3[x1*bytes+3],
                                           src3[x2*bytes+3],
                                           src3[x3*bytes+3],
                                           src3[x4*bytes+3]);
               p3  = cubic_spline_fit (dx, src4[x1*bytes+3],
                                           src4[x2*bytes+3],
                                           src4[x3*bytes+3],
                                           src4[x4*bytes+3]);
               alphasum = cubic_spline_fit (dy, p0, p1, p2, p3);                                          
               if (alphasum > 0)
                 {
                   for (b = 0 ; b < 3 ; b++)
                     {
                       p0  = cubic_spline_fit (dx, src1[x1*bytes+b] * src1[x1*bytes+3],
                                                   src1[x2*bytes+b] * src1[x2*bytes+3],
                                                   src1[x3*bytes+b] * src1[x3*bytes+3],
                                                   src1[x4*bytes+b] * src1[x4*bytes+3]);
                       p1  = cubic_spline_fit (dx, src2[x1*bytes+b] * src2[x1*bytes+3],
                                                   src2[x2*bytes+b] * src2[x2*bytes+3],
                                                   src2[x3*bytes+b] * src2[x3*bytes+3],
                                                   src2[x4*bytes+b] * src2[x4*bytes+3]);
                       p2  = cubic_spline_fit (dx, src3[x1*bytes+b] * src3[x1*bytes+3],
                                                   src3[x2*bytes+b] * src3[x2*bytes+3],
                                                   src3[x3*bytes+b] * src3[x3*bytes+3],
                                                   src3[x4*bytes+b] * src3[x4*bytes+3]);
                       p3  = cubic_spline_fit (dx, src4[x1*bytes+b] * src4[x1*bytes+3],
                                                   src4[x2*bytes+b] * src4[x2*bytes+3],
                                                   src4[x3*bytes+b] * src4[x3*bytes+3],
                                                   src4[x4*bytes+b] * src4[x4*bytes+3]);
                       sum = cubic_spline_fit (dy, p0, p1, p2, p3);                                          
                       sum /= alphasum;
                       dst[x * bytes+b] = (guchar) CLAMP (sum, 0, 255);
                     }
                     dst[x * bytes+3] = (guchar) CLAMP (alphasum, 0, 255);
                 }
               else
                {
                   for (b = 0 ; b < 3 ; b++)
                     {
                       dst[x * bytes+b] = 0;
                     }
                }
               break;
           }
        }

      pixel_region_set_row (dstPR, 0, y, dstPR->w, dst);
    }

  g_free(src0);
  g_free(dst);
}
                                 
static void  scale_cubic        (PixelRegion  *dstPR,
                                 const guchar *src,
                                 gint          source_w,
                                 gint          source_h,
                                 gint          bytes)
{
  gint       x, y, b;
  gint       x1, x2, x3, x4;
  gint       y1, y2, y3, y4;

  gdouble    scalex, scaley;
  gdouble    dx, dy;
  gdouble    p0, p1, p2, p3;

  gdouble       alphasum;
  gdouble       sum;

  guchar *dst;
  guchar *src0;
  guchar *src1;
  guchar *src2;
  guchar *src3;
  guchar *src4;

  scaley = (gdouble) dstPR->h / (gdouble) source_h;
  scalex = (gdouble) dstPR->w / (gdouble) source_w;

  /* Interpolate */  
  src0 = g_new0 (guchar, source_w * bytes * 5);
  src1 = src0 + source_w * bytes;
  src2 = src1 + source_w * bytes;
  src3 = src2 + source_w * bytes;
  src4 = src3 + source_w * bytes; 
  
  dst  = g_new0 (guchar, dstPR->w * bytes);

  for (y = 0; y < dstPR->h; y++)
    {
      dy = (y / scaley);
      y2 = (gint) floor (dy);     
      y1 = ( y2-1 >= 0 )       ? y2-1: y2+1;
      y3 = (gint) ceil (dy);
      y3 = ( y3 < source_h)    ? y3  : y2;
      y4 = ( y3+1 < source_h ) ? y3+1: y3-1;
      dy =  dy - y2;
      
      memcpy(src1, src+y1*source_w*bytes,source_w*bytes);
      memcpy(src2, src+y2*source_w*bytes,source_w*bytes);
      memcpy(src3, src+y3*source_w*bytes,source_w*bytes);
      memcpy(src4, src+y4*source_w*bytes,source_w*bytes);

      for (x = 0; x < dstPR->w; x++)
        {
          dx = (x / scalex);        
          x2 = (gint)floor(dx);
          x1 = ( x2-1 >= 0 )       ? x2-1: x2+1;
          x3 = (gint)ceil(dx);
          x3 = ( x3 < source_w)    ? x3  : x2;
          x4 = ( x3+1 < source_w ) ? x3+1: x3-1;
          dx =  dx - x2;
          
          switch (bytes)
            {                
             case 1:
               p0  = cubic_spline_fit (dx, src1[x1*bytes],
                                           src1[x2*bytes],
                                           src1[x3*bytes],
                                           src1[x4*bytes]);
               p1  = cubic_spline_fit (dx, src2[x1*bytes],
                                           src2[x2*bytes],
                                           src2[x3*bytes],
                                           src2[x4*bytes]);
               p2  = cubic_spline_fit (dx, src3[x1*bytes],
                                           src3[x2*bytes],
                                           src3[x3*bytes],
                                           src3[x4*bytes]);
               p3  = cubic_spline_fit (dx, src4[x1*bytes],
                                           src4[x2*bytes],
                                           src4[x3*bytes],
                                           src4[x4*bytes]);
               sum = cubic_spline_fit (dy, p0, p1, p2, p3);                                          
               dst[x*bytes] = (guchar) CLAMP (sum, 0, 255);
               break;
             case 2:
               p0  = cubic_spline_fit (dx, src1[x1*bytes+1],
                                           src1[x2*bytes+1],
                                           src1[x3*bytes+1],
                                           src1[x4*bytes+1]);
               p1  = cubic_spline_fit (dx, src2[x1*bytes+1],
                                           src2[x2*bytes+1],
                                           src2[x3*bytes+1],
                                           src2[x4*bytes+1]);
               p2  = cubic_spline_fit (dx, src3[x1*bytes+1],
                                           src3[x2*bytes+1],
                                           src3[x3*bytes+1],
                                           src3[x4*bytes+1]);
               p3  = cubic_spline_fit (dx, src4[x1*bytes+1],
                                           src4[x2*bytes+1],
                                           src4[x3*bytes+1],
                                           src4[x4*bytes+1]);
               alphasum = cubic_spline_fit (dy, p0, p1, p2, p3);                                          
               if (alphasum > 0)
                 {
                   p0  = cubic_spline_fit (dx, src1[x1*bytes] * src1[x1*bytes+1],
                                               src1[x2*bytes] * src1[x2*bytes+1],
                                               src1[x3*bytes] * src1[x3*bytes+1],
                                               src1[x4*bytes] * src1[x4*bytes+1]);
                   p1  = cubic_spline_fit (dx, src2[x1*bytes] * src2[x1*bytes+1],
                                               src2[x2*bytes] * src2[x2*bytes+1],
                                               src2[x3*bytes] * src2[x3*bytes+1],
                                               src2[x4*bytes] * src2[x4*bytes+1]);
                   p2  = cubic_spline_fit (dx, src3[x1*bytes] * src3[x1*bytes+1],
                                               src3[x2*bytes] * src3[x2*bytes+1],
                                               src3[x3*bytes] * src3[x3*bytes+1],
                                               src3[x4*bytes] * src3[x4*bytes+1]);
                   p3  = cubic_spline_fit (dx, src4[x1*bytes] * src4[x1*bytes+1],
                                               src4[x2*bytes] * src4[x2*bytes+1],
                                               src4[x3*bytes] * src4[x3*bytes+1],
                                               src4[x4*bytes] * src4[x4*bytes+1]);
                   sum = cubic_spline_fit (dy, p0, p1, p2, p3);                                          
                   sum /= alphasum;
                   dst[x * bytes]   = (guchar) CLAMP (sum, 0, 255);
                   dst[x * bytes+1] = (guchar) CLAMP (alphasum, 0, 255);
                 }
               else
                {
                   dst[x * bytes]   = 0;
                   dst[x * bytes+1] = 0;                
                }
               break;
             case 3:
               for (b = 0 ; b < bytes ; b++)
                 {
                   p0  = cubic_spline_fit (dx, src1[x1*bytes+b],
                                               src1[x2*bytes+b],
                                               src1[x3*bytes+b],
                                               src1[x4*bytes+b]);
                   p1  = cubic_spline_fit (dx, src2[x1*bytes+b],
                                               src2[x2*bytes+b],
                                               src2[x3*bytes+b],
                                               src2[x4*bytes+b]);
                   p2  = cubic_spline_fit (dx, src3[x1*bytes+b],
                                               src3[x2*bytes+b],
                                               src3[x3*bytes+b],
                                               src3[x4*bytes+b]);
                   p3  = cubic_spline_fit (dx, src4[x1*bytes+b],
                                               src4[x2*bytes+b],
                                               src4[x3*bytes+b],
                                               src4[x4*bytes+b]);
                   sum = cubic_spline_fit (dy, p0, p1, p2, p3);                                          
                   dst[x*bytes+b] = (guchar) CLAMP (sum, 0, 255);
                 }
               break;
             case 4:
               p0  = cubic_spline_fit (dx, src1[x1*bytes+3],
                                           src1[x2*bytes+3],
                                           src1[x3*bytes+3],
                                           src1[x4*bytes+3]);
               p1  = cubic_spline_fit (dx, src2[x1*bytes+3],
                                           src2[x2*bytes+3],
                                           src2[x3*bytes+3],
                                           src2[x4*bytes+3]);
               p2  = cubic_spline_fit (dx, src3[x1*bytes+3],
                                           src3[x2*bytes+3],
                                           src3[x3*bytes+3],
                                           src3[x4*bytes+3]);
               p3  = cubic_spline_fit (dx, src4[x1*bytes+3],
                                           src4[x2*bytes+3],
                                           src4[x3*bytes+3],
                                           src4[x4*bytes+3]);
               alphasum = cubic_spline_fit (dy, p0, p1, p2, p3);                                          
               if (alphasum > 0)
                 {
                   for (b = 0 ; b < 3 ; b++)
                     {
                       p0  = cubic_spline_fit (dx, src1[x1*bytes+b] * src1[x1*bytes+3],
                                                   src1[x2*bytes+b] * src1[x2*bytes+3],
                                                   src1[x3*bytes+b] * src1[x3*bytes+3],
                                                   src1[x4*bytes+b] * src1[x4*bytes+3]);
                       p1  = cubic_spline_fit (dx, src2[x1*bytes+b] * src2[x1*bytes+3],
                                                   src2[x2*bytes+b] * src2[x2*bytes+3],
                                                   src2[x3*bytes+b] * src2[x3*bytes+3],
                                                   src2[x4*bytes+b] * src2[x4*bytes+3]);
                       p2  = cubic_spline_fit (dx, src3[x1*bytes+b] * src3[x1*bytes+3],
                                                   src3[x2*bytes+b] * src3[x2*bytes+3],
                                                   src3[x3*bytes+b] * src3[x3*bytes+3],
                                                   src3[x4*bytes+b] * src3[x4*bytes+3]);
                       p3  = cubic_spline_fit (dx, src4[x1*bytes+b] * src4[x1*bytes+3],
                                                   src4[x2*bytes+b] * src4[x2*bytes+3],
                                                   src4[x3*bytes+b] * src4[x3*bytes+3],
                                                   src4[x4*bytes+b] * src4[x4*bytes+3]);
                       sum = cubic_spline_fit (dy, p0, p1, p2, p3);                                          
                       sum /= alphasum;
                       dst[x * bytes+b] = (guchar) CLAMP (sum, 0, 255);
                     }
                     dst[x * bytes+3] = (guchar) CLAMP (alphasum, 0, 255);
                 }
               else
                {
                   for (b = 0 ; b < 3 ; b++)
                     {
                       dst[x * bytes+b] = 0;
                     }
                }
               break;
           }
        }

      pixel_region_set_row (dstPR, 0, y, dstPR->w, dst);
    }

  g_free(src0);
  g_free(dst);
}                                 


static void
scale_lanczos_pr (PixelRegion      *srcPR,
                  PixelRegion      *dstPR)
{
  gfloat        *kernel_lookup  = NULL;       /* Lanczos lookup table                */
  gdouble        x_kernel[LANCZOS_WIDTH2],    /* 1-D kernels of Lanczos window coeffs */
                 y_kernel[LANCZOS_WIDTH2];

  gdouble        newval;                      /* new interpolated RGB value          */

  guchar        *win_buf = NULL;              /* Sliding window buffer               */
    guchar        *win_ptr[LANCZOS_WIDTH2];     /* Ponters to sliding window rows      */

  guchar        *dst_buf = NULL;              /* Pointer to destination image data   */

  gint           x, y;                        /* Position in destination image       */
  gint           i, byte;                     /* loop vars  */
  gint           row;

  gdouble        trans[6], itrans[6];         /* Scale transformations               */

  const gint     dst_width     = dstPR->w;
  const gint     dst_height    = dstPR->h;
  const gint     bytes         = dstPR->bytes;
  const gint     src_width     = srcPR->w;
  const gint     src_height    = srcPR->h;

  const gint     src_row_span  = src_width * bytes;
  const gint     dst_row_span  = dst_width * bytes;
  const gint     win_row_span  = (src_width + LANCZOS_WIDTH2) * bytes;

  const gdouble  scale_x       = dst_width / (gdouble) src_width;
  const gdouble  scale_y       = dst_height / (gdouble) src_height;

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

  for (i = 0; i < 6; i++)
    trans[i] = 0.0;

  trans[0] = scale_x;
  trans[4] = scale_y;

  if (! inv_lin_trans (trans, itrans))
    {
      g_warning ("transformation matrix is not invertible");
      return;
    }

  /* allocate buffer for destination row */
  dst_buf = g_new0 (guchar, dst_row_span);

  /* allocate and fill kernel_lookup lookup table */
  kernel_lookup = create_lanczos_lookup ();

  /* allocate buffer for source rows */
  win_buf = g_new0 (guchar, win_row_span * LANCZOS_WIDTH2);

  /* Set the window pointers */
  for ( i = 0 ; i < LANCZOS_WIDTH2 ; i++ )
    win_ptr[i] = win_buf + (win_row_span * i) + LANCZOS_WIDTH * bytes;

  /* fill the data for the first loop */
  for ( i = 0 ; i <= LANCZOS_WIDTH && i < src_height ; i++)
    pixel_region_get_row (srcPR, 0, i, src_width, win_ptr[i + LANCZOS_WIDTH], 1);

  for (row = y = 0; y < dst_height; y++)
    {
      pixel_region_get_row (dstPR, 0, y, dst_width, dst_buf, 1);
      for  (x = 0; x < dst_width; x++)
        {
          gdouble  dsrc_x ,dsrc_y;        /* corresponding scaled position in source image */
          gint     int_src_x, int_src_y;  /* integer part of coordinates in source image   */
          gint     x_shift, y_shift;      /* index into Lanczos lookup                     */
          gdouble  kx_sum, ky_sum;        /* sums of Lanczos kernel coeffs                 */

          /* -0.5 corrects image drift.due to average offset used in lookup */
          dsrc_x = x / scale_x - 0.5;
          dsrc_y = y / scale_y - 0.5;

          /* avoid (int) on negative*/
          if (dsrc_x > 0)
            int_src_x = (gint) (dsrc_x);
          else
            int_src_x = (gint) (dsrc_x + 1) - 1;

          if (dsrc_y > 0)
            int_src_y = (gint) (dsrc_y);
          else
            int_src_y = (gint) (dsrc_y + 1) - 1;

          /* calc lookup offsets for non-interger remainders */
          x_shift = (gint) ((dsrc_x - int_src_x) * LANCZOS_SPP + 0.5);
          y_shift = (gint) ((dsrc_y - int_src_y) * LANCZOS_SPP + 0.5);

          /*  Fill x_kernel[] and y_kernel[] with lanczos coeffs
           *
           *  kernel_lookup = Is a lookup table that contains half of the symetrical windowed-sinc func.
           *
           *  x_shift, y_shift = shift from kernel center due to fractional part
           *           of interpollation
           *
           *  The for-loop creates two 1-D kernels for convolution.
           *    - If the center position +/- LANCZOS_WIDTH is out of
           *      the source image coordinates set the value to 0.0
           *      FIXME => partial kernel. Define a more rigourous border mode.
           *    - If the kernel index is out of range set value to 0.0
           *      ( caused by offset coeff. obselete??)
           */
          kx_sum = ky_sum = 0.0;

          for (i = LANCZOS_WIDTH; i >= -LANCZOS_WIDTH; i--)
            {
              gint pos = i * LANCZOS_SPP;

              if ( int_src_x + i >= 0 && int_src_x + i < src_width)
                kx_sum += x_kernel[LANCZOS_WIDTH + i] = kernel_lookup[ABS (x_shift - pos)];
              else
                x_kernel[LANCZOS_WIDTH + i] = 0.0;

              if ( int_src_y + i >= 0 && int_src_y + i < src_height)
                ky_sum += y_kernel[LANCZOS_WIDTH + i] = kernel_lookup[ABS (y_shift - pos)];
              else
                y_kernel[LANCZOS_WIDTH + i] = 0.0;
            }

          /* normalise the kernel arrays */
          for (i = -LANCZOS_WIDTH; i <= LANCZOS_WIDTH; i++)
          {
             x_kernel[LANCZOS_WIDTH + i] /= kx_sum;
             y_kernel[LANCZOS_WIDTH + i] /= ky_sum;
          }

          /*
            Scaling up
            New determined source row is > than last read row
            rotate the pointers and get next source row from region.
            If no more source rows are available fill buffer with 0
            ( Probably not necessary because multipliers should be 0).
          */
          for ( ; row < int_src_y ; )
            {
              row++;
              rotate_pointers (win_ptr, LANCZOS_WIDTH2);
              if ( row + LANCZOS_WIDTH < src_height)
                pixel_region_get_row (srcPR, 0,
                                      row + LANCZOS_WIDTH, src_width,
                                      win_ptr[LANCZOS_WIDTH2 - 1], 1);
              else
                memset (win_ptr[LANCZOS_WIDTH2 - 1], 0,
                        sizeof (guchar) * src_row_span);
            }
           /*
              Scaling down
            */
            for ( ; row > int_src_y ; )
            {
              row--;
              for ( i = 0 ; i < LANCZOS_WIDTH2 - 1 ; i++ )
                 rotate_pointers (win_ptr, LANCZOS_WIDTH2);
              if ( row >=  0)
                pixel_region_get_row (srcPR, 0,
                                      row, src_width,
                                      win_ptr[0], 1);
              else
                memset (win_ptr[0], 0,
                        sizeof (guchar) * src_row_span);

            }


           if (pixel_region_has_alpha (srcPR))
             {
               const gint alpha = bytes - 1;
               gint       byte;
               gdouble    arecip;
               gdouble    aval;

               aval = 0.0;
               for (i = 0; i < LANCZOS_WIDTH2 ; i++ )
                 aval += y_kernel[i] * lanczos_sum (win_ptr[i], x_kernel,
                         int_src_x, bytes, alpha);

               if (aval <= 0.0)
                 {
                   arecip = 0.0;
                   dst_buf[x * bytes + alpha] = 0;
                 }
               else if (aval > 255.0)
                 {
                   arecip = 1.0 / aval;
                   dst_buf[x * bytes + alpha] = 255;
                 }
               else
                 {
                   arecip = 1.0 / aval;
                   dst_buf[x * bytes + alpha] = RINT (aval);
                 }

               for (byte = 0; byte < alpha; byte++)
                 {
                   newval = 0.0;
                   for (i = 0; i < LANCZOS_WIDTH2; i++ )
                     newval += y_kernel[i] * lanczos_sum_mul (win_ptr[i], x_kernel,
                               int_src_x, bytes, byte, alpha);
                   newval *= arecip;
                   dst_buf[x * bytes + byte] = CLAMP (newval, 0, 255);
                 }
             }
           else
             {
               for (byte = 0; byte < bytes; byte++)
                 {
                   /* Calculate new value */
                   newval = 0.0;
                   for (i = 0; i < LANCZOS_WIDTH2; i++ )
                     newval += y_kernel[i] * lanczos_sum (win_ptr[i], x_kernel,
                               int_src_x, bytes, byte);
                   dst_buf[x * bytes + byte] = CLAMP ((gint) newval, 0, 255);
                 }
             }
        }

      pixel_region_set_row (dstPR, 0, y , dst_width, dst_buf);
    }

  g_free (dst_buf);
  g_free (win_buf);
  g_free (kernel_lookup);
}

static void  scale_lanczos      (PixelRegion  *dstPR,
                                 const guchar *src,
                                 gint          source_w,
                                 gint          source_h,
                                 gint          bytes)
{
  gfloat        *kernel_lookup  = NULL;       /* Lanczos lookup table                */
  gdouble        x_kernel[LANCZOS_WIDTH2],    /* 1-D kernels of Lanczos window coeffs */
                 y_kernel[LANCZOS_WIDTH2];

  gdouble        newval;                      /* new interpolated RGB value          */

  guchar        *win_buf = NULL;              /* Sliding window buffer               */
  guchar        *win_ptr[LANCZOS_WIDTH2];     /* Ponters to sliding window rows      */

  guchar        *dst_buf = NULL;              /* Pointer to destination image data   */

  gint           x, y;                        /* Position in destination image       */
  gint           i, byte;                     /* loop vars  */
  gint           row;

  gdouble        trans[6], itrans[6];         /* Scale transformations               */

  const gint     dst_width     = dstPR->w;
  const gint     dst_height    = dstPR->h;
  const gint     src_width     = source_w;
  const gint     src_height    = source_h;

  const gint     src_row_span  = src_width * bytes;
  const gint     dst_row_span  = dst_width * bytes;
  const gint     win_row_span  = (src_width + LANCZOS_WIDTH2) * bytes;

  const gdouble  scale_x       = dst_width / (gdouble) src_width;
  const gdouble  scale_y       = dst_height / (gdouble) src_height;

  for (i = 0; i < 6; i++)
    trans[i] = 0.0;

  trans[0] = scale_x;
  trans[4] = scale_y;

  if (! inv_lin_trans (trans, itrans))
    {
      g_warning ("transformation matrix is not invertible");
      return;
    }

  /* allocate buffer for destination row */
  dst_buf = g_new0 (guchar, dst_row_span);

  /* allocate and fill kernel_lookup lookup table */
  kernel_lookup = create_lanczos_lookup ();

  /* allocate buffer for source rows */
  win_buf = g_new0 (guchar, win_row_span * LANCZOS_WIDTH2);

  /* Set the window pointers */
  for ( i = 0 ; i < LANCZOS_WIDTH2 ; i++ )
    win_ptr[i] = win_buf + (win_row_span * i) + LANCZOS_WIDTH * bytes;

  /* fill the data for the first loop */
  for ( i = 0 ; i <= LANCZOS_WIDTH && i < src_height ; i++)
    memcpy(win_ptr[i + LANCZOS_WIDTH], src+(src_width*bytes*i), src_width*bytes);

  for (row = y = 0; y < dst_height; y++)
    {
      pixel_region_get_row (dstPR, 0, y, dst_width, dst_buf, 1);
      for  (x = 0; x < dst_width; x++)
        {
          gdouble  dsrc_x ,dsrc_y;        /* corresponding scaled position in source image */
          gint     int_src_x, int_src_y;  /* integer part of coordinates in source image   */
          gint     x_shift, y_shift;      /* index into Lanczos lookup                     */
          gdouble  kx_sum, ky_sum;        /* sums of Lanczos kernel coeffs                 */

          /* -0.5 corrects image drift.due to average offset used in lookup */
          dsrc_x = x / scale_x - 0.5;
          dsrc_y = y / scale_y - 0.5;

          /* avoid (int) on negative*/
          if (dsrc_x > 0)
            int_src_x = (gint) (dsrc_x);
          else
            int_src_x = (gint) (dsrc_x + 1) - 1;

          if (dsrc_y > 0)
            int_src_y = (gint) (dsrc_y);
          else
            int_src_y = (gint) (dsrc_y + 1) - 1;

          /* calc lookup offsets for non-interger remainders */
          x_shift = (gint) ((dsrc_x - int_src_x) * LANCZOS_SPP + 0.5);
          y_shift = (gint) ((dsrc_y - int_src_y) * LANCZOS_SPP + 0.5);

          /*  Fill x_kernel[] and y_kernel[] with lanczos coeffs
           *
           *  kernel_lookup = Is a lookup table that contains half of the symetrical windowed-sinc func.
           *
           *  x_shift, y_shift = shift from kernel center due to fractional part
           *           of interpollation
           *
           *  The for-loop creates two 1-D kernels for convolution.
           *    - If the center position +/- LANCZOS_WIDTH is out of
           *      the source image coordinates set the value to 0.0
           *      FIXME => partial kernel. Define a more rigourous border mode.
           *    - If the kernel index is out of range set value to 0.0
           *      ( caused by offset coeff. obselete??)
           */
          kx_sum = ky_sum = 0.0;

          for (i = LANCZOS_WIDTH; i >= -LANCZOS_WIDTH; i--)
            {
              gint pos = i * LANCZOS_SPP;

              if ( int_src_x + i >= 0 && int_src_x + i < src_width)
                kx_sum += x_kernel[LANCZOS_WIDTH + i] = kernel_lookup[ABS (x_shift - pos)];
              else
                x_kernel[LANCZOS_WIDTH + i] = 0.0;

              if ( int_src_y + i >= 0 && int_src_y + i < src_height)
                ky_sum += y_kernel[LANCZOS_WIDTH + i] = kernel_lookup[ABS (y_shift - pos)];
              else
                y_kernel[LANCZOS_WIDTH + i] = 0.0;
            }

          /* normalise the kernel arrays */
          for (i = -LANCZOS_WIDTH; i <= LANCZOS_WIDTH; i++)
          {
             x_kernel[LANCZOS_WIDTH + i] /= kx_sum;
             y_kernel[LANCZOS_WIDTH + i] /= ky_sum;
          }

          /*
            Scaling up
            New determined source row is > than last read row
            rotate the pointers and get next source row from region.
            If no more source rows are available fill buffer with 0
            ( Probably not necessary because multipliers should be 0).
          */
          for ( ; row < int_src_y ; )
            {
              row++;
              rotate_pointers (win_ptr, LANCZOS_WIDTH2);
              if ( row + LANCZOS_WIDTH < src_height)
                memcpy(win_ptr[LANCZOS_WIDTH2 - 1], src+((row + LANCZOS_WIDTH)*src_width*bytes), src_width*bytes);
              else
                memset (win_ptr[LANCZOS_WIDTH2 - 1], 0,
                        sizeof (guchar) * src_row_span);
            }
           /*
              Scaling down
            */
            for ( ; row > int_src_y ; )
            {
              row--;
              for ( i = 0 ; i < LANCZOS_WIDTH2 - 1 ; i++ )
                 rotate_pointers (win_ptr, LANCZOS_WIDTH2);
              if ( row >=  0)
                memcpy(win_ptr[0], src+(row * src_width * bytes), src_width*bytes);
              else
                memset (win_ptr[0], 0,
                        sizeof (guchar) * src_row_span);

            }


           if (pixel_region_has_alpha (dstPR))
             {
               const gint alpha = bytes - 1;
               gint       byte;
               gdouble    arecip;
               gdouble    aval;

               aval = 0.0;
               for (i = 0; i < LANCZOS_WIDTH2 ; i++ )
                 aval += y_kernel[i] * lanczos_sum (win_ptr[i], x_kernel,
                         int_src_x, bytes, alpha);

               if (aval <= 0.0)
                 {
                   arecip = 0.0;
                   dst_buf[x * bytes + alpha] = 0;
                 }
               else if (aval > 255.0)
                 {
                   arecip = 1.0 / aval;
                   dst_buf[x * bytes + alpha] = 255;
                 }
               else
                 {
                   arecip = 1.0 / aval;
                   dst_buf[x * bytes + alpha] = RINT (aval);
                 }

               for (byte = 0; byte < alpha; byte++)
                 {
                   newval = 0.0;
                   for (i = 0; i < LANCZOS_WIDTH2; i++ )
                     newval += y_kernel[i] * lanczos_sum_mul (win_ptr[i], x_kernel,
                               int_src_x, bytes, byte, alpha);
                   newval *= arecip;
                   dst_buf[x * bytes + byte] = CLAMP (newval, 0, 255);
                 }
             }
           else
             {
               for (byte = 0; byte < bytes; byte++)
                 {
                   /* Calculate new value */
                   newval = 0.0;
                   for (i = 0; i < LANCZOS_WIDTH2; i++ )
                     newval += y_kernel[i] * lanczos_sum (win_ptr[i], x_kernel,
                               int_src_x, bytes, byte);
                   dst_buf[x * bytes + byte] = CLAMP ((gint) newval, 0, 255);
                 }
             }
        }

      pixel_region_set_row (dstPR, 0, y , dst_width, dst_buf);
    }

  g_free (dst_buf);
  g_free (win_buf);
  g_free (kernel_lookup);
}




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;
}

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;
}

static gboolean
inv_lin_trans (const gdouble *t,
               gdouble       *it)
{
  gdouble d = (t[0] * t[4]) - (t[1] * t[3]);  /* determinant */

  if (fabs(d) < EPSILON )
    return FALSE;

  it[0] =    t[4] / d;
  it[1] =   -t[1] / d;
  it[2] = (( t[1] * t[5]) - (t[2] * t[4])) / d;
  it[3] =   -t[3] / d;
  it[4] =    t[0] / d;
  it[5] = (( t[2] * t[3]) - (t[0] * t[5])) / d;

  return TRUE;
}
static inline void
rotate_pointers (guchar  **p,
                 guint32   n)
{
  guchar  *tmp = p[0];
  guint32  i;

  for (i = 0; i < n-1; i++)
    p[i] = p[i+1];

  p[i] = tmp;
}
