gimp/app/core/gimp-transform-region.c

/* GIMP - The GNU Image Manipulation Program
 * Copyright (C) 1995-2003 Spencer Kimball, Peter Mattis, and others
 *
 * 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 <stdlib.h>

#include <glib-object.h>

#include "libgimpbase/gimpbase.h"
#include "libgimpmath/gimpmath.h"

#include "core-types.h"

#include "base/pixel-region.h"
#include "base/pixel-surround.h"
#include "base/tile-manager.h"
#include "base/tile.h"

#include "paint-funcs/scale-funcs.h"

#include "gimp-transform-region.h"
#include "gimpchannel.h"
#include "gimpcontext.h"
#include "gimpimage.h"
#include "gimppickable.h"
#include "gimpprogress.h"


/*  forward function prototypes  */

static void  gimp_transform_region_nearest (TileManager       *orig_tiles,
                                            PixelRegion       *destPR,
                                            gint               dest_x1,
                                            gint               dest_y1,
                                            gint               dest_x2,
                                            gint               dest_y2,
                                            gint               u1,
                                            gint               v1,
                                            gint               u2,
                                            gint               v2,
                                            const GimpMatrix3 *m,
                                            gint               alpha,
                                            const guchar      *bg_color,
                                            GimpProgress      *progress);
static void  gimp_transform_region_linear  (TileManager       *orig_tiles,
                                            PixelRegion       *destPR,
                                            gint               dest_x1,
                                            gint               dest_y1,
                                            gint               dest_x2,
                                            gint               dest_y2,
                                            gint               u1,
                                            gint               v1,
                                            gint               u2,
                                            gint               v2,
                                            const GimpMatrix3 *m,
                                            gint               alpha,
                                            gboolean           supersample,
                                            gint               recursion_level,
                                            const guchar      *bg_color,
                                            GimpProgress      *progress);
static void  gimp_transform_region_cubic   (TileManager       *orig_tiles,
                                            PixelRegion       *destPR,
                                            gint               dest_x1,
                                            gint               dest_y1,
                                            gint               dest_x2,
                                            gint               dest_y2,
                                            gint               u1,
                                            gint               v1,
                                            gint               u2,
                                            gint               v2,
                                            const GimpMatrix3 *m,
                                            gint               alpha,
                                            gboolean           supersample,
                                            gint               recursion_level,
                                            const guchar      *bg_color,
                                            GimpProgress      *progress);
static void  gimp_transform_region_lanczos (TileManager       *orig_tiles,
                                            PixelRegion       *destPR,
                                            gint               dest_x1,
                                            gint               dest_y1,
                                            gint               dest_x2,
                                            gint               dest_y2,
                                            gint               u1,
                                            gint               v1,
                                            gint               u2,
                                            gint               v2,
                                            const GimpMatrix3 *m,
                                            gint               alpha,
                                            const guchar      *bg_color,
                                            GimpProgress      *progress);

static inline void  normalize_coordinates (const gint     coords,
                                           const gdouble *tu,
                                           const gdouble *tv,
                                           const gdouble *tw,
                                           gdouble       *u,
                                           gdouble       *v);

static inline gboolean supersample_dtest (gdouble u0, gdouble v0,
                                          gdouble u1, gdouble v1,
                                          gdouble u2, gdouble v2,
                                          gdouble u3, gdouble v3);

static void     sample_adapt      (TileManager   *tm,
                                   gdouble        uc,     gdouble     vc,
                                   gdouble        u0,     gdouble     v0,
                                   gdouble        u1,     gdouble     v1,
                                   gdouble        u2,     gdouble     v2,
                                   gdouble        u3,     gdouble     v3,
                                   gint           level,
                                   guchar        *color,
                                   const guchar  *bg_color,
                                   gint           bpp,
                                   gint           alpha);

static void     sample_cubic      (PixelSurround *surround,
                                   gdouble        u,
                                   gdouble        v,
                                   guchar        *color,
                                   gint           bytes,
                                   gint           alpha);

static void     sample_linear     (PixelSurround *surround,
                                   gdouble        u,
                                   gdouble        v,
                                   guchar        *color,
                                   gint           bytes,
                                   gint           alpha);


/*  public functions  */

void
gimp_transform_region (GimpPickable          *pickable,
                       GimpContext           *context,
                       TileManager           *orig_tiles,
                       PixelRegion           *destPR,
                       gint                   dest_x1,
                       gint                   dest_y1,
                       gint                   dest_x2,
                       gint                   dest_y2,
                       const GimpMatrix3     *matrix,
                       GimpInterpolationType  interpolation_type,
                       gboolean               supersample,
                       gint                   recursion_level,
                       GimpProgress          *progress)
{
  GimpImageType  pickable_type;
  GimpMatrix3    m;
  gint           u1, v1, u2, v2;       /* source bounding box */
  gint           alpha;
  guchar         bg_color[MAX_CHANNELS];

  g_return_if_fail (GIMP_IS_PICKABLE (pickable));

  tile_manager_get_offsets (orig_tiles, &u1, &v1);

  u2 = u1 + tile_manager_width (orig_tiles);
  v2 = v1 + tile_manager_height (orig_tiles);

  m = *matrix;
  gimp_matrix3_invert (&m);

  alpha = 0;

  /*  turn interpolation off for simple transformations (e.g. rot90)  */
  if (gimp_matrix3_is_simple (matrix))
    interpolation_type = GIMP_INTERPOLATION_NONE;

  pickable_type = gimp_pickable_get_image_type (pickable);

  /*  Get the background color  */
  gimp_image_get_background (gimp_pickable_get_image (pickable), context,
                             pickable_type, bg_color);

  switch (GIMP_IMAGE_TYPE_BASE_TYPE (pickable_type))
    {
    case GIMP_RGB:
      bg_color[ALPHA_PIX] = TRANSPARENT_OPACITY;
      alpha = ALPHA_PIX;
      break;

    case GIMP_GRAY:
      bg_color[ALPHA_G_PIX] = TRANSPARENT_OPACITY;
      alpha = ALPHA_G_PIX;
      break;

    case GIMP_INDEXED:
      bg_color[ALPHA_I_PIX] = TRANSPARENT_OPACITY;
      alpha = ALPHA_I_PIX;
      /*  If the image is indexed color, ignore interpolation value  */
      interpolation_type = GIMP_INTERPOLATION_NONE;
      break;

    default:
      g_assert_not_reached ();
      break;
    }

  /*  "Outside" a channel is transparency, not the bg color  */
  if (GIMP_IS_CHANNEL (pickable))
    bg_color[0] = TRANSPARENT_OPACITY;

  /*  setting alpha = 0 will cause the channel's value to be treated
   *  as alpha and the color channel loops never to be entered
   */
  if (tile_manager_bpp (orig_tiles) == 1)
    alpha = 0;

  switch (interpolation_type)
    {
    case GIMP_INTERPOLATION_NONE:
      gimp_transform_region_nearest (orig_tiles, destPR,
                                     dest_x1, dest_y1, dest_x2, dest_y2,
                                     u1, v1, u2, v2,
                                     &m, alpha, bg_color, progress);
      break;

    case GIMP_INTERPOLATION_LINEAR:
      gimp_transform_region_linear (orig_tiles, destPR,
                                    dest_x1, dest_y1, dest_x2, dest_y2,
                                    u1, v1, u2, v2,
                                    &m, alpha, supersample, recursion_level,
                                    bg_color, progress);
      break;

    case GIMP_INTERPOLATION_CUBIC:
      gimp_transform_region_cubic (orig_tiles, destPR,
                                   dest_x1, dest_y1, dest_x2, dest_y2,
                                   u1, v1, u2, v2,
                                   &m, alpha, supersample, recursion_level,
                                   bg_color, progress);
      break;

    case GIMP_INTERPOLATION_LANCZOS:
      gimp_transform_region_lanczos (orig_tiles, destPR,
                                     dest_x1, dest_y1, dest_x2, dest_y2,
                                     u1, v1, u2, v2,
                                     &m, alpha, bg_color, progress);
      break;
    }
}

static void
gimp_transform_region_nearest (TileManager        *orig_tiles,
                               PixelRegion        *destPR,
                               gint                dest_x1,
                               gint                dest_y1,
                               gint                dest_x2,
                               gint                dest_y2,
                               gint                u1,
                               gint                v1,
                               gint                u2,
                               gint                v2,
                               const GimpMatrix3  *m,
                               gint                alpha,
                               const guchar       *bg_color,
                               GimpProgress       *progress)
{
  gdouble   uinc, vinc, winc;  /* increments in source coordinates  */
  gint      pixels;
  gint      total;
  gint      n;
  gpointer  pr;

  uinc = m->coeff[0][0];
  vinc = m->coeff[1][0];
  winc = m->coeff[2][0];

  total = destPR->w * destPR->h;

  for (pr = pixel_regions_register (1, destPR), pixels = 0, n = 0;
       pr != NULL;
       pr = pixel_regions_process (pr), n++)
    {
      guchar *dest = destPR->data;
      gint    y;

      for (y = destPR->y; y < destPR->y + destPR->h; y++)
        {
          gint     x     = dest_x1 + destPR->x;
          gint     width = destPR->w;
          guchar  *d     = dest;
          gdouble  tu, tv, tw;   /* undivided source coordinates and divisor */

          /* set up inverse transform steps */
          tu = uinc * x + m->coeff[0][1] * (dest_y1 + y) + m->coeff[0][2];
          tv = vinc * x + m->coeff[1][1] * (dest_y1 + y) + m->coeff[1][2];
          tw = winc * x + m->coeff[2][1] * (dest_y1 + y) + m->coeff[2][2];

          while (width--)
            {
              gdouble u, v; /* source coordinates */
              gint    iu, iv;

              /*  normalize homogeneous coords  */
              normalize_coordinates (1, &tu, &tv, &tw, &u, &v);

              iu = (gint) u;
              iv = (gint) v;

              /*  Set the destination pixels  */
              if (iu >= u1 && iu < u2 &&
                  iv >= v1 && iv < v2)
                {
                  read_pixel_data_1 (orig_tiles, iu - u1, iv - v1, d);

                  d += destPR->bytes;
                }
              else /* not in source range */
                {
                  gint b;

                  for (b = 0; b < destPR->bytes; b++)
                    *d++ = bg_color[b];
                }

              tu += uinc;
              tv += vinc;
              tw += winc;
            }

          dest += destPR->rowstride;
        }

      if (progress)
        {
          pixels += destPR->w * destPR->h;

          if (n % 16 == 0)
            gimp_progress_set_value (progress,
                                     (gdouble) pixels / (gdouble) total);
        }
    }
}

static void
gimp_transform_region_linear (TileManager        *orig_tiles,
                              PixelRegion        *destPR,
                              gint                dest_x1,
                              gint                dest_y1,
                              gint                dest_x2,
                              gint                dest_y2,
                              gint                u1,
                              gint                v1,
                              gint                u2,
                              gint                v2,
                              const GimpMatrix3  *m,
                              gint                alpha,
                              gboolean            supersample,
                              gint                recursion_level,
                              const guchar       *bg_color,
                              GimpProgress       *progress)
{
  PixelSurround *surround;
  gdouble        uinc, vinc, winc;  /* increments in source coordinates  */
  gint           pixels;
  gint           total;
  gint           n;
  gpointer       pr;

  surround = pixel_surround_new (orig_tiles, 2, 2, bg_color);

  uinc = m->coeff[0][0];
  vinc = m->coeff[1][0];
  winc = m->coeff[2][0];

  total = destPR->w * destPR->h;

  for (pr = pixel_regions_register (1, destPR), pixels = 0, n = 0;
       pr != NULL;
       pr = pixel_regions_process (pr), n++)
    {
      guchar *dest = destPR->data;
      gint    y;

      for (y = destPR->y; y < destPR->y + destPR->h; y++)
        {
          gint     x     = dest_x1 + destPR->x;
          gint     yy    = dest_y1 + y;
          gint     width = destPR->w;
          guchar  *d     = dest;

          gdouble  tu[5], tv[5];   /* undivided source coordinates */
          gdouble  tw[5];          /* divisor                      */

          /* set up inverse transform steps */
          tu[0] = uinc * (x    ) + m->coeff[0][1] * (yy    ) + m->coeff[0][2];
          tv[0] = vinc * (x    ) + m->coeff[1][1] * (yy    ) + m->coeff[1][2];
          tw[0] = winc * (x    ) + m->coeff[2][1] * (yy    ) + m->coeff[2][2];

          tu[1] = uinc * (x - 1) + m->coeff[0][1] * (yy    ) + m->coeff[0][2];
          tv[1] = vinc * (x - 1) + m->coeff[1][1] * (yy    ) + m->coeff[1][2];
          tw[1] = winc * (x - 1) + m->coeff[2][1] * (yy    ) + m->coeff[2][2];

          tu[2] = uinc * (x    ) + m->coeff[0][1] * (yy - 1) + m->coeff[0][2];
          tv[2] = vinc * (x    ) + m->coeff[1][1] * (yy - 1) + m->coeff[1][2];
          tw[2] = winc * (x    ) + m->coeff[2][1] * (yy - 1) + m->coeff[2][2];

          tu[3] = uinc * (x + 1) + m->coeff[0][1] * (yy    ) + m->coeff[0][2];
          tv[3] = vinc * (x + 1) + m->coeff[1][1] * (yy    ) + m->coeff[1][2];
          tw[3] = winc * (x + 1) + m->coeff[2][1] * (yy    ) + m->coeff[2][2];

          tu[4] = uinc * (x    ) + m->coeff[0][1] * (yy + 1) + m->coeff[0][2];
          tv[4] = vinc * (x    ) + m->coeff[1][1] * (yy + 1) + m->coeff[1][2];
          tw[4] = winc * (x    ) + m->coeff[2][1] * (yy + 1) + m->coeff[2][2];

          while (width--)
            {
              gdouble u[5], v[5]; /* source coordinates */
              gint    i;

              /*  normalize homogeneous coords  */
              normalize_coordinates (5, tu, tv, tw, u, v);

              /*  Set the destination pixels  */
              if (u[0] <  u1 || v[0] <  v1 ||
                  u[0] >= u2 || v[0] >= v2 )
                {
                  /* not in source range */
                  /* increment the destination pointers  */

                  for (i = 0; i < destPR->bytes; i++)
                    *d++ = bg_color[i];
                }
              else
                {
                  /* clamp texture coordinates */
                  for (i = 0; i < 5; i++)
                    {
                      u[i] = CLAMP (u[i], u1, u2 - 1);
                      v[i] = CLAMP (v[i], v1, v2 - 1);
                    }

                  if (supersample &&
                      supersample_dtest (u[1], v[1], u[2], v[2],
                                         u[3], v[3], u[4], v[4]))
                    {
                      sample_adapt (orig_tiles,
                                    u[0] - u1, v[0] - v1,
                                    u[1] - u1, v[1] - v1,
                                    u[2] - u1, v[2] - v1,
                                    u[3] - u1, v[3] - v1,
                                    u[4] - u1, v[4] - v1,
                                    recursion_level,
                                    d, bg_color, destPR->bytes, alpha);
                    }
                  else
                    {
                      sample_linear (surround, u[0] - u1, v[0] - v1,
                                     d, destPR->bytes, alpha);
                    }

                  d += destPR->bytes;
                }

              for (i = 0; i < 5; i++)
                {
                  tu[i] += uinc;
                  tv[i] += vinc;
                  tw[i] += winc;
                }
            }

          dest += destPR->rowstride;
        }

      if (progress)
        {
          pixels += destPR->w * destPR->h;

          if (n % 16 == 0)
            gimp_progress_set_value (progress,
                                     (gdouble) pixels / (gdouble) total);
        }
    }

  pixel_surround_destroy (surround);
}

static void
gimp_transform_region_cubic (TileManager        *orig_tiles,
                             PixelRegion        *destPR,
                             gint                dest_x1,
                             gint                dest_y1,
                             gint                dest_x2,
                             gint                dest_y2,
                             gint                u1,
                             gint                v1,
                             gint                u2,
                             gint                v2,
                             const GimpMatrix3  *m,
                             gint                alpha,
                             gboolean            supersample,
                             gint                recursion_level,
                             const guchar       *bg_color,
                             GimpProgress       *progress)
{
  PixelSurround *surround;
  gdouble        uinc, vinc, winc;  /* increments in source coordinates  */
  gint           pixels;
  gint           total;
  gint           n;
  gpointer       pr;

  surround = pixel_surround_new (orig_tiles, 4, 4, bg_color);

  uinc = m->coeff[0][0];
  vinc = m->coeff[1][0];
  winc = m->coeff[2][0];

  total = destPR->w * destPR->h;

  for (pr = pixel_regions_register (1, destPR), pixels = 0, n = 0;
       pr != NULL;
       pr = pixel_regions_process (pr), n++)
    {
      guchar *dest = destPR->data;
      gint    y;

      for (y = destPR->y; y < destPR->y + destPR->h; y++)
        {
          gint     x     = dest_x1 + destPR->x;
          gint     yy    = dest_y1 + y;
          gint     width = destPR->w;
          guchar  *d     = dest;

          gdouble  tu[5], tv[5];   /* undivided source coordinates */
          gdouble  tw[5];          /* divisor                      */

          /* set up inverse transform steps */
          tu[0] = uinc * (x    ) + m->coeff[0][1] * (yy    ) + m->coeff[0][2];
          tv[0] = vinc * (x    ) + m->coeff[1][1] * (yy    ) + m->coeff[1][2];
          tw[0] = winc * (x    ) + m->coeff[2][1] * (yy    ) + m->coeff[2][2];

          tu[1] = uinc * (x - 1) + m->coeff[0][1] * (yy    ) + m->coeff[0][2];
          tv[1] = vinc * (x - 1) + m->coeff[1][1] * (yy    ) + m->coeff[1][2];
          tw[1] = winc * (x - 1) + m->coeff[2][1] * (yy    ) + m->coeff[2][2];

          tu[2] = uinc * (x    ) + m->coeff[0][1] * (yy - 1) + m->coeff[0][2];
          tv[2] = vinc * (x    ) + m->coeff[1][1] * (yy - 1) + m->coeff[1][2];
          tw[2] = winc * (x    ) + m->coeff[2][1] * (yy - 1) + m->coeff[2][2];

          tu[3] = uinc * (x + 1) + m->coeff[0][1] * (yy    ) + m->coeff[0][2];
          tv[3] = vinc * (x + 1) + m->coeff[1][1] * (yy    ) + m->coeff[1][2];
          tw[3] = winc * (x + 1) + m->coeff[2][1] * (yy    ) + m->coeff[2][2];

          tu[4] = uinc * (x    ) + m->coeff[0][1] * (yy + 1) + m->coeff[0][2];
          tv[4] = vinc * (x    ) + m->coeff[1][1] * (yy + 1) + m->coeff[1][2];
          tw[4] = winc * (x    ) + m->coeff[2][1] * (yy + 1) + m->coeff[2][2];

          while (width--)
            {
              gdouble u[5], v[5]; /* source coordinates */
              gint    i;

              /*  normalize homogeneous coords  */
              normalize_coordinates (5, tu, tv, tw, u, v);

              /*  Set the destination pixels  */
              if (u[0] <  u1 || v[0] <  v1 ||
                  u[0] >= u2 || v[0] >= v2 )
                {
                  /* not in source range */
                  /* increment the destination pointers  */

                  for (i = 0; i < destPR->bytes; i++)
                    *d++ = bg_color[i];
                }
              else
                {
                  /* clamp texture coordinates */
                  for (i = 0; i < 5; i++)
                    {
                      u[i] = CLAMP (u[i], u1, u2 - 1);
                      v[i] = CLAMP (v[i], v1, v2 - 1);
                    }

                  if (supersample &&
                      supersample_dtest (u[1], v[1], u[2], v[2],
                                         u[3], v[3], u[4], v[4]))
                    {
                      sample_adapt (orig_tiles,
                                    u[0] - u1, v[0] - v1,
                                    u[1] - u1, v[1] - v1,
                                    u[2] - u1, v[2] - v1,
                                    u[3] - u1, v[3] - v1,
                                    u[4] - u1, v[4] - v1,
                                    recursion_level,
                                    d, bg_color, destPR->bytes, alpha);
                    }
                  else
                    {
                      sample_cubic (surround, u[0] - u1, v[0] - v1,
                                    d, destPR->bytes, alpha);
                    }

                  d += destPR->bytes;
                }

              for (i = 0; i < 5; i++)
                {
                  tu[i] += uinc;
                  tv[i] += vinc;
                  tw[i] += winc;
                }
            }

          dest += destPR->rowstride;
        }

      if (progress)
        {
          pixels += destPR->w * destPR->h;

          if (n % 16 == 0)
            gimp_progress_set_value (progress,
                                     (gdouble) pixels / (gdouble) total);
        }
    }

  pixel_surround_destroy (surround);
}

static void
gimp_transform_region_lanczos (TileManager       *orig_tiles,
                               PixelRegion       *destPR,
                               gint               dest_x1,
                               gint               dest_y1,
                               gint               dest_x2,
                               gint               dest_y2,
                               gint               u1,
                               gint               v1,
                               gint               u2,
                               gint               v2,
                               const GimpMatrix3 *m,
                               gint               alpha,
                               const guchar      *bg_color,
                               GimpProgress      *progress)
{
  gfloat  *lanczos;                  /* Lanczos lookup table         */
  gdouble  x_kernel[LANCZOS_WIDTH2]; /* 1-D kernels of window coeffs */
  gdouble  y_kernel[LANCZOS_WIDTH2];
  gdouble  x_sum, y_sum;             /* sum of Lanczos weights       */

  gdouble  uw;
  gdouble  ww;
  gdouble  vw;
  gdouble  du;
  gdouble  dv;

  guchar  *dest, *d;
  gint     x, y;
  gint     width;
  gint     bytes;
  gint     pos;
  gdouble  newval;
  gdouble  arecip;
  gdouble  aval;

  guchar   lwin[LANCZOS_WIDTH2 * LANCZOS_WIDTH2][MAX_CHANNELS + 1];
  gint     b, u, v, i, j;
  gint     pu, pv, su, sv;

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

  width = dest_x2 - dest_x1;
  bytes = destPR->bytes;

  dest = g_new (guchar, width * bytes);

  for (y = dest_y1; y < dest_y2; y++)
    {
      if (progress && !(y & 0xf))
        gimp_progress_set_value (progress,
                                 (gdouble) (y - dest_y1) /
                                 (gdouble) (dest_y2 - dest_y1));

      pixel_region_get_row (destPR, 0, (y - dest_y1), width, dest, 1);

      d = dest;

      for (x = dest_x1; x < dest_x2; x++)
        {
          uw = m->coeff[0][0] * x + m->coeff[0][1] * y + m->coeff[0][2];
          vw = m->coeff[1][0] * x + m->coeff[1][1] * y + m->coeff[1][2];
          ww = m->coeff[2][0] * x + m->coeff[2][1] * y + m->coeff[2][2];

          normalize_coordinates (1, &uw, &vw, &ww, &du, &dv);

          u = (gint) du;
          v = (gint) dv;

          /* get weight for fractional error */
          su = (gint) ((du - u) * LANCZOS_SPP);
          sv = (gint) ((dv - v) * LANCZOS_SPP);

          if (u <  u1 || v <  v1 ||
              u >= u2 || v >= v2)
            {
              /* not in source range */
              /* increment the destination pointers  */
              for (b = 0; b < bytes; b++)
                *d++ = bg_color[b];
            }
          else
            {
              pos = 0;

              for (j = 0; j < LANCZOS_WIDTH2 ; j++)
                for (i = 0; i < LANCZOS_WIDTH2 ; i++, pos++)
                  {
                    pu = CLAMP (u + i - LANCZOS_WIDTH, u1, u2 - 1);
                    pv = CLAMP (v + j - LANCZOS_WIDTH, v1, v2 - 1);

                    read_pixel_data_1 (orig_tiles,
                                       pu - u1, pv - v1, &lwin[pos][0]);
                  }

              /* fill 1D kernels */
              for (x_sum = y_sum = 0.0, i = LANCZOS_WIDTH;
                   i >= -LANCZOS_WIDTH;
                   i--)
                {
                  pos = i * LANCZOS_SPP;
                  x_sum += x_kernel[LANCZOS_WIDTH + i] = lanczos[ABS (su - pos)];
                  y_sum += y_kernel[LANCZOS_WIDTH + i] = lanczos[ABS (sv - pos)];
                }

              /* normalise the weighted arrays */
              for (i = 0; i < LANCZOS_WIDTH2 ; i++)
                {
                  x_kernel[i] /= x_sum;
                  y_kernel[i] /= y_sum;
                }

              pos = 0;
              aval = 0.0;

              for (j = 0; j < LANCZOS_WIDTH2 ; j++)
                for (i = 0; i < LANCZOS_WIDTH2 ; i++, pos++)
                  {
                    aval += y_kernel[j] * x_kernel[i] * lwin[pos][alpha];
                  }

              if (aval <= 0.0)
                {
                  arecip = 0.0;
                  aval = 0;
                }
              else if (aval > 255.0)
                {
                  arecip = 1.0 / aval;
                  aval = 255;
                }
              else
                {
                  arecip = 1.0 / aval;
                }

              for (b = 0; b < alpha; b++)
                {
                  pos = 0;
                  newval = 0.0;

                  for (j = 0; j < LANCZOS_WIDTH2 ; j++)
                    for (i = 0; i < LANCZOS_WIDTH2 ; i++, pos++)
                      {
                        newval += y_kernel[j] * x_kernel[i] * lwin[pos][b] * lwin[pos][alpha];
                      }
                  newval *= arecip;
                  *d++ = CLAMP (newval, 0, 255);
                }

              *d++ = RINT (aval);
            }
        }

      /*  set the pixel region row  */
      pixel_region_set_row (destPR, 0, (y - dest_y1), width, dest);
    }

  g_free (dest);
  g_free (lanczos);

  if (progress)
    gimp_progress_set_value (progress, 1.0);
}


/*  private functions  */

static inline void
normalize_coordinates (const gint     coords,
                       const gdouble *tu,
                       const gdouble *tv,
                       const gdouble *tw,
                       gdouble       *u,
                       gdouble       *v)
{
  gint i;

  /*  normalize homogeneous coords  */
  for (i = 0; i < coords; i++)
    {
      if (tw[i] == 1.0)
        {
          u[i] = tu[i];
          v[i] = tv[i];
        }
      else if (tw[i] != 0.0)
        {
          u[i] = tu[i] / tw[i];
          v[i] = tv[i] / tw[i];
        }
      else
        {
          g_warning ("homogeneous coordinate = 0...\n");

          u[i] = tu[i];
          v[i] = tv[i];
        }
    }
}


#define BILINEAR(jk, j1k, jk1, j1k1, dx, dy) \
                ((1 - dy) * (jk  + dx * (j1k  - jk)) + \
                      dy  * (jk1 + dx * (j1k1 - jk1)))

  /*  u & v are the subpixel coordinates of the point in
   *  the original selection's floating buffer.
   *  We need the two pixel coords around them:
   *  iu to iu + 1, iv to iv + 1
   */
static void
sample_linear (PixelSurround *surround,
               gdouble        u,
               gdouble        v,
               guchar        *color,
               gint           bytes,
               gint           alpha)
{
  gdouble       a_val, a_recip;
  gint          i;
  gint          iu = floor (u);
  gint          iv = floor (v);
  gint          rowstride;
  gdouble       du, dv;
  const guchar *alphachan;
  const guchar *data;

  /* lock the pixel surround */
  data = pixel_surround_lock (surround, iu, iv, &rowstride);

  /* the fractional error */
  du = u - iu;
  dv = v - iv;

  /* calculate alpha value of result pixel */
  alphachan = &data[alpha];
  a_val = BILINEAR (alphachan[0],         alphachan[bytes],
                    alphachan[rowstride], alphachan[rowstride + bytes], du, dv);

  if (a_val <= 0.0)
    {
      a_recip = 0.0;
      color[alpha] = 0.0;
    }
  else if (a_val >= 255.0)
    {
      a_recip = 1.0 / a_val;
      color[alpha] = 255;
    }
  else
    {
      a_recip = 1.0 / a_val;
      color[alpha] = RINT (a_val);
    }

  /*  for colour channels c,
   *  result = bilinear (c * alpha) / bilinear (alpha)
   *
   *  never entered for alpha == 0
   */
  for (i = 0; i < alpha; i++)
    {
      gint newval;

      newval = (a_recip *
                BILINEAR (alphachan[0]                 * data[i],
                          alphachan[bytes]             * data[bytes + i],
                          alphachan[rowstride]         * data[rowstride + i],
                          alphachan[rowstride + bytes] * data[rowstride + bytes + i],
                          du, dv));

      color[i] = CLAMP (newval, 0, 255);
    }
}

/* macros to handle conversion to/from fixed point, this fixed point code
 * uses signed integers, by using 8 bits for the fractional part we have
 *
 *  1 bit  sign
 * 21 bits integer part
 *  8 bit  fractional part
 *
 * 1023 discrete subpixel sample positions should be enough for the needs
 * of the supersampling algorithm, drawables where the dimensions have a need
 * exceeding 2^21 ( 2097152px, will typically use terabytes of memory, when
 * that is the common need, we can probably assume 64 bit integers and adjust
 * FIXED_SHIFT accordingly.
 */
#define FIXED_SHIFT          10
#define FIXED_UNIT           (1 << FIXED_SHIFT)
#define DOUBLE2FIXED(val)    ((val) * FIXED_UNIT)
#define FIXED2DOUBLE(val)    ((val) / FIXED_UNIT)

/*
    bilinear interpolation of a fixed point pixel
*/
static void
sample_bi (TileManager  *tm,
           gint          x,
           gint          y,
           guchar       *color,
           const guchar *bg_color,
           gint          bpp,
           gint          alpha)
{
  guchar C[4][4];
  gint   i;
  gint   xscale = (x & (FIXED_UNIT-1));
  gint   yscale = (y & (FIXED_UNIT-1));

  gint   x0 = x >> FIXED_SHIFT;
  gint   y0 = y >> FIXED_SHIFT;
  gint   x1 = x0 + 1;
  gint   y1 = y0 + 1;


  /*  fill the color with default values, since read_pixel_data_1
   *  does nothing, when accesses are out of bounds.
   */
  for (i = 0; i < 4; i++)
    *(guint*) (&C[i]) = *(guint*) (bg_color);

  read_pixel_data_1 (tm, x0, y0, C[0]);
  read_pixel_data_1 (tm, x1, y0, C[2]);
  read_pixel_data_1 (tm, x0, y1, C[1]);
  read_pixel_data_1 (tm, x1, y1, C[3]);

#define lerp(v1,v2,r) \
        (((guint)(v1) * (FIXED_UNIT - (guint)(r)) + \
          (guint)(v2) * (guint)(r)) >> FIXED_SHIFT)

  color[alpha]= lerp (lerp (C[0][alpha], C[1][alpha], yscale),
                      lerp (C[2][alpha], C[3][alpha], yscale), xscale);

  if (color[alpha])
    { /* to avoid problems, calculate with premultiplied alpha */
      for (i=0; i<alpha; i++)
        {
          C[0][i] = (C[0][i] * C[0][alpha] / 255);
          C[1][i] = (C[1][i] * C[1][alpha] / 255);
          C[2][i] = (C[2][i] * C[2][alpha] / 255);
          C[3][i] = (C[3][i] * C[3][alpha] / 255);
        }

      for (i = 0; i < alpha; i++)
        color[i] = lerp (lerp (C[0][i], C[1][i], yscale),
                         lerp (C[2][i], C[3][i], yscale), xscale);
    }
  else
    {
      for (i = 0; i < alpha; i++)
        color[i] = 0;
    }
#undef lerp
}



/*
 * Returns TRUE if one of the deltas of the
 * quad edge is > 1.0 (16.16 fixed values).
 */
static inline gboolean
supersample_test (gint x0, gint y0,
                  gint x1, gint y1,
                  gint x2, gint y2,
                  gint x3, gint y3)
{
  return (abs (x0 - x1) > FIXED_UNIT ||
          abs (x1 - x2) > FIXED_UNIT ||
          abs (x2 - x3) > FIXED_UNIT ||
          abs (x3 - x0) > FIXED_UNIT ||

          abs (y0 - y1) > FIXED_UNIT ||
          abs (y1 - y2) > FIXED_UNIT ||
          abs (y2 - y3) > FIXED_UNIT ||
          abs (y3 - y0) > FIXED_UNIT);
}

/*
 *  Returns TRUE if one of the deltas of the
 *  quad edge is > 1.0 (double values).
 */
static inline gboolean
supersample_dtest (gdouble x0, gdouble y0,
                   gdouble x1, gdouble y1,
                   gdouble x2, gdouble y2,
                   gdouble x3, gdouble y3)
{
  return (fabs (x0 - x1) > 1.0 ||
          fabs (x1 - x2) > 1.0 ||
          fabs (x2 - x3) > 1.0 ||
          fabs (x3 - x0) > 1.0 ||

          fabs (y0 - y1) > 1.0 ||
          fabs (y1 - y2) > 1.0 ||
          fabs (y2 - y3) > 1.0 ||
          fabs (y3 - y0) > 1.0);
}

/*
    sample a grid that is spaced according to the quadraliteral's edges,
    it subdivides a maximum of level times before sampling.
    0..3 is a cycle around the quad
*/
static void
get_sample (TileManager  *tm,
            gint          xc,
            gint          yc,
            gint          x0,
            gint          y0,
            gint          x1,
            gint          y1,
            gint          x2,
            gint          y2,
            gint          x3,
            gint          y3,
            gint         *cc,
            gint          level,
            guint        *color,
            const guchar *bg_color,
            gint          bpp,
            gint          alpha)
{
  if (!level || !supersample_test (x0, y0, x1, y1, x2, y2, x3, y3))
    {
      gint   i;
      guchar C[4];

      sample_bi (tm, xc, yc, C, bg_color, bpp, alpha);

      for (i = 0; i < bpp; i++)
        color[i]+= C[i];

      (*cc)++;  /* increase number of samples taken */
    }
  else
    {
      gint tx, lx, rx, bx, tlx, trx, blx, brx;
      gint ty, ly, ry, by, tly, try, bly, bry;

      /* calculate subdivided corner coordinates (including centercoords
         thus using a bilinear interpolation,. almost as good as
         doing the perspective transform for each subpixel coordinate*/

      tx  = (x0 + x1) / 2;
      tlx = (x0 + xc) / 2;
      trx = (x1 + xc) / 2;
      lx  = (x0 + x3) / 2;
      rx  = (x1 + x2) / 2;
      blx = (x3 + xc) / 2;
      brx = (x2 + xc) / 2;
      bx  = (x3 + x2) / 2;

      ty  = (y0 + y1) / 2;
      tly = (y0 + yc) / 2;
      try = (y1 + yc) / 2;
      ly  = (y0 + y3) / 2;
      ry  = (y1 + y2) / 2;
      bly = (y3 + yc) / 2;
      bry = (y2 + yc) / 2;
      by  = (y3 + y2) / 2;

      get_sample (tm,
                  tlx,tly,
                  x0,y0, tx,ty, xc,yc, lx,ly,
                  cc, level-1, color, bg_color, bpp, alpha);

      get_sample (tm,
                  trx,try,
                  tx,ty, x1,y1, rx,ry, xc,yc,
                  cc, level-1, color, bg_color, bpp, alpha);

      get_sample (tm,
                  brx,bry,
                  xc,yc, rx,ry, x2,y2, bx,by,
                  cc, level-1, color, bg_color, bpp, alpha);

      get_sample (tm,
                  blx,bly,
                  lx,ly, xc,yc, bx,by, x3,y3,
                  cc, level-1, color, bg_color, bpp, alpha);
    }
}

static void
sample_adapt (TileManager  *tm,
              gdouble       xc,
              gdouble       yc,
              gdouble       x0,
              gdouble       y0,
              gdouble       x1,
              gdouble       y1,
              gdouble       x2,
              gdouble       y2,
              gdouble       x3,
              gdouble       y3,
              gint          level,
              guchar       *color,
              const guchar *bg_color,
              gint          bpp,
              gint          alpha)
{
    gint  cc = 0;
    gint  i;
    guint C[MAX_CHANNELS];

    C[0] = C[1] = C[2] = C[3] = 0;

    get_sample (tm,
                DOUBLE2FIXED (xc), DOUBLE2FIXED (yc),
                DOUBLE2FIXED (x0), DOUBLE2FIXED (y0),
                DOUBLE2FIXED (x1), DOUBLE2FIXED (y1),
                DOUBLE2FIXED (x2), DOUBLE2FIXED (y2),
                DOUBLE2FIXED (x3), DOUBLE2FIXED (y3),
                &cc, level, C, bg_color, bpp, alpha);

    if (!cc)
      cc=1;

    color[alpha] = C[alpha] / cc;

    if (color[alpha])
      {
         /* go from premultiplied to postmultiplied alpha */
        for (i = 0; i < alpha; i++)
          color[i] = ((C[i] / cc) * 255) / color[alpha];
      }
    else
      {
        for (i = 0; i < alpha; i++)
          color[i] = 0;
      }
}

/* access interleaved pixels */
#define CUBIC_ROW(dx, row, step) \
  gimp_drawable_transform_cubic(dx,\
            (row)[0], (row)[step], (row)[step+step], (row)[step+step+step])

#define CUBIC_SCALED_ROW(dx, row, arow, step) \
  gimp_drawable_transform_cubic(dx, \
            (arow)[0]              * (row)[0], \
            (arow)[step]           * (row)[step], \
            (arow)[step+step]      * (row)[step+step], \
            (arow)[step+step+step] * (row)[step+step+step])


/*  Note: cubic function no longer clips result. */
/*  Inlining this function makes sample_cubic() run about 10% faster. (Sven) */
static inline gdouble
gimp_drawable_transform_cubic (gdouble dx,
                               gint    jm1,
                               gint    j,
                               gint    jp1,
                               gint    jp2)
{
  gdouble result;

#if 0
  /* Equivalent to Gimp 1.1.1 and earlier - some ringing */
  result = ((( ( - jm1 + j - jp1 + jp2 ) * dx +
               ( jm1 + jm1 - j - j + jp1 - jp2 ) ) * dx +
               ( - jm1 + jp1 ) ) * dx + j );
  /* Recommended by Mitchell and Netravali - too blurred? */
  result = ((( ( - 7 * jm1 + 21 * j - 21 * jp1 + 7 * jp2 ) * dx +
               ( 15 * jm1 - 36 * j + 27 * jp1 - 6 * jp2 ) ) * dx +
               ( - 9 * jm1 + 9 * jp1 ) ) * dx + (jm1 + 16 * j + jp1) ) / 18.0;
#endif

  /* Catmull-Rom - not bad */
  result = ((( ( - jm1 + 3 * j - 3 * jp1 + jp2 ) * dx +
               ( 2 * jm1 - 5 * j + 4 * jp1 - jp2 ) ) * dx +
               ( - jm1 + jp1 ) ) * dx + (j + j) ) / 2.0;

  return result;
}


  /*  u & v are the subpixel coordinates of the point in
   *  the original selection's floating buffer.
   *  We need the four integer pixel coords around them:
   *  iu to iu + 3, iv to iv + 3
   */
static void
sample_cubic (PixelSurround *surround,
              gdouble        u,
              gdouble        v,
              guchar        *color,
              gint           bytes,
              gint           alpha)
{
  gdouble       a_val, a_recip;
  gint          i;
  gint          iu = floor(u);
  gint          iv = floor(v);
  gint          rowstride;
  gdouble       du, dv;
  const guchar *data;

  /* lock the pixel surround */
  data = pixel_surround_lock (surround, iu - 1 , iv - 1, &rowstride);

  /* the fractional error */
  du = u - iu;
  dv = v - iv;

  /* calculate alpha of result */
  a_val = gimp_drawable_transform_cubic
    (dv,
     CUBIC_ROW (du, data + alpha + rowstride * 0, bytes),
     CUBIC_ROW (du, data + alpha + rowstride * 1, bytes),
     CUBIC_ROW (du, data + alpha + rowstride * 2, bytes),
     CUBIC_ROW (du, data + alpha + rowstride * 3, bytes));

  if (a_val <= 0.0)
    {
      a_recip = 0.0;
      color[alpha] = 0;
    }
  else if (a_val > 255.0)
    {
      a_recip = 1.0 / a_val;
      color[alpha] = 255;
    }
  else
    {
      a_recip = 1.0 / a_val;
      color[alpha] = RINT (a_val);
    }

  /*  for colour channels c,
   *  result = bicubic (c * alpha) / bicubic (alpha)
   *
   *  never entered for alpha == 0
   */
  for (i = 0; i < alpha; i++)
    {
      gint newval = (a_recip *
                     gimp_drawable_transform_cubic
                     (dv,
                      CUBIC_SCALED_ROW (du,
                                        i + data + rowstride * 0,
                                        data + alpha + rowstride * 0,
                                        bytes),
                      CUBIC_SCALED_ROW (du,
                                        i + data + rowstride * 1,
                                        data + alpha + rowstride * 1,
                                        bytes),
                      CUBIC_SCALED_ROW (du,
                                        i + data + rowstride * 2,
                                        data + alpha + rowstride * 2,
                                        bytes),
                      CUBIC_SCALED_ROW (du,
                                        i + data + rowstride * 3,
                                        data + alpha + rowstride * 3,
                                        bytes)));

      color[i] = CLAMP (newval, 0, 255);
    }
}