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-region.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,
                                            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,
                                            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,
                                            gint               recursion_level,
                                            const guchar      *bg_color,
                                            GimpProgress      *progress);

static inline void  untransform_coords     (const GimpMatrix3 *m,
                                            gint               x,
                                            gint               y,
                                            gdouble           *tu,
                                            gdouble           *tv,
                                            gdouble           *tw);
static inline void  normalize_coords       (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      (PixelSurround *surround,
                                   gdouble        uc,
                                   gdouble        vc,
                                   gdouble        u0,
                                   gdouble        v0,
                                   gdouble        u1,
                                   gdouble        v1,
                                   gdouble        u2,
                                   gdouble        v2,
                                   gdouble        u3,
                                   gdouble        v3,
                                   gint           level,
                                   guchar        *color,
                                   gint           bpp,
                                   gint           alpha);

static void     sample_linear     (PixelSurround *surround,
                                   gdouble        u,
                                   gdouble        v,
                                   guchar        *color,
                                   gint           bytes,
                                   gint           alpha);
static void     sample_cubic      (PixelSurround *surround,
                                   gdouble        u,
                                   gdouble        v,
                                   guchar        *color,
                                   gint           bytes,
                                   gint           alpha);
static void     sample_lanczos    (PixelSurround *surround,
                                   const gfloat  *lanczos,
                                   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,
                       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, 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, 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, recursion_level,
                                     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_coords (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,
                              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, PIXEL_SURROUND_BACKGROUND);
  pixel_surround_set_bg (surround, 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++)
        {
          guchar  *d     = dest;
          gint     width = destPR->w;
          gdouble  tu[5], tv[5];   /* undivided source coordinates */
          gdouble  tw[5];          /* divisor                      */

          /* set up inverse transform steps */
          untransform_coords (m, dest_x1 + destPR->x, dest_y1 + y, tu, tv, tw);

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

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

              /*  Set the destination pixels  */
              if (supersample_dtest (u[1], v[1], u[2], v[2],
                                     u[3], v[3], u[4], v[4]))
                {
                  sample_adapt (surround,
                                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, 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,
                             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, PIXEL_SURROUND_BACKGROUND);
  pixel_surround_set_bg (surround, 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++)
        {
          guchar  *d     = dest;
          gint     width = destPR->w;
          gdouble  tu[5], tv[5];   /* undivided source coordinates */
          gdouble  tw[5];          /* divisor                      */

          /* set up inverse transform steps */
          untransform_coords (m, dest_x1 + destPR->x, dest_y1 + y, tu, tv, tw);

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

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

              if (supersample_dtest (u[1], v[1], u[2], v[2],
                                     u[3], v[3], u[4], v[4]))
                {
                  sample_adapt (surround,
                                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, 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,
                               gint               recursion_level,
                               const guchar      *bg_color,
                               GimpProgress      *progress)
{
  PixelSurround *surround;
  gfloat        *lanczos;           /* Lanczos lookup table              */
  gdouble        uinc, vinc, winc;  /* increments in source coordinates  */
  gint           pixels;
  gint           total;
  gint           n;
  gpointer       pr;

  surround = pixel_surround_new (orig_tiles,
                                 LANCZOS_WIDTH2, LANCZOS_WIDTH2,
                                 PIXEL_SURROUND_BACKGROUND);
  pixel_surround_set_bg (surround, bg_color);

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

  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++)
        {
          guchar  *d     = dest;
          gint     width = destPR->w;
          gdouble  tu[5], tv[5];   /* undivided source coordinates */
          gdouble  tw[5];          /* divisor                      */

          /* set up inverse transform steps */
          untransform_coords (m, dest_x1 + destPR->x, dest_y1 + y, tu, tv, tw);

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

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

              if (supersample_dtest (u[1], v[1], u[2], v[2],
                                     u[3], v[3], u[4], v[4]))
                {
                  sample_adapt (surround,
                                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, destPR->bytes, alpha);
                }
              else
                {
                  sample_lanczos (surround, lanczos, 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);
        }
   }

  g_free (lanczos);
  pixel_surround_destroy (surround);
}


/*  private functions  */

static inline void
untransform_coords (const GimpMatrix3 *m,
                    const gint         x,
                    const gint         y,
                    gdouble           *tu,
                    gdouble           *tv,
                    gdouble           *tw)
{
  tu[0] = m->coeff[0][0] * (x    ) + m->coeff[0][1] * (y    ) + m->coeff[0][2];
  tv[0] = m->coeff[1][0] * (x    ) + m->coeff[1][1] * (y    ) + m->coeff[1][2];
  tw[0] = m->coeff[2][0] * (x    ) + m->coeff[2][1] * (y    ) + m->coeff[2][2];

  tu[1] = m->coeff[0][0] * (x - 1) + m->coeff[0][1] * (y    ) + m->coeff[0][2];
  tv[1] = m->coeff[1][0] * (x - 1) + m->coeff[1][1] * (y    ) + m->coeff[1][2];
  tw[1] = m->coeff[2][0] * (x - 1) + m->coeff[2][1] * (y    ) + m->coeff[2][2];

  tu[2] = m->coeff[0][0] * (x    ) + m->coeff[0][1] * (y - 1) + m->coeff[0][2];
  tv[2] = m->coeff[1][0] * (x    ) + m->coeff[1][1] * (y - 1) + m->coeff[1][2];
  tw[2] = m->coeff[2][0] * (x    ) + m->coeff[2][1] * (y - 1) + m->coeff[2][2];

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

  tu[4] = m->coeff[0][0] * (x    ) + m->coeff[0][1] * (y + 1) + m->coeff[0][2];
  tv[4] = m->coeff[1][0] * (x    ) + m->coeff[1][1] * (y + 1) + m->coeff[1][2];
  tw[4] = m->coeff[2][0] * (x    ) + m->coeff[2][1] * (y + 1) + m->coeff[2][2];
}

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

  for (i = 0; i < coords; i++)
    {
      if (G_LIKELY (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,
               const gdouble  u,
               const gdouble  v,
               guchar        *color,
               const gint     bytes,
               const gint     alpha)
{
  gdouble       a_val, a_recip;
  gint          i;
  const gint    iu = floor (u);
  const 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 =
        ROUND ((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 (PixelSurround *surround,
           const gint     x,
           const gint     y,
           guchar        *color,
           const gint     bpp,
           const gint     alpha)
{
  const gint    xscale = x & (FIXED_UNIT - 1);
  const gint    yscale = y & (FIXED_UNIT - 1);
  const gint    x0     = x >> FIXED_SHIFT;
  const gint    y0     = y >> FIXED_SHIFT;
  gint          rowstride;
  const guchar *src    = pixel_surround_lock (surround, x0, y0, &rowstride);
  const guchar *s0     = src;
  const guchar *s1     = src + bpp;
  const guchar *s2     = src + rowstride;
  const guchar *s3     = src + rowstride + bpp;
  gint          i;

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

  color[alpha]= lerp (lerp (s0[alpha], s1[alpha], yscale),
                      lerp (s2[alpha], s3[alpha], yscale), xscale);

  if (color[alpha])
    {
      /* to avoid problems, calculate with premultiplied alpha */
      for (i = 0; i < alpha; i++)
        {
          color[i] = lerp (lerp (s0[i] * s0[alpha] / 255,
                                 s1[i] * s1[alpha] / 255, yscale),
                           lerp (s2[i] * s2[alpha] / 255,
                                 s3[i] * s3[alpha] / 255, 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). This is the condition used on whether additional recursive
 * subdivision should be used.
 */
static inline gboolean
supersample_test (const gint x0, const gint y0,
                  const gint x1, const gint y1,
                  const gint x2, const gint y2,
                  const gint x3, const 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 > sqrt(2) (double
 * values). This is the condition used on whether supersampling should be used
 * or not. By making this sqrt(2) supersampling will not be triggered by
 * rotations.
 */
static inline gboolean
supersample_dtest (const gdouble x0, const gdouble y0,
                   const gdouble x1, const gdouble y1,
                   const gdouble x2, const gdouble y2,
                   const gdouble x3, const gdouble y3)
{
  return (fabs (x0 - x1) > G_SQRT2 ||
          fabs (x1 - x2) > G_SQRT2 ||
          fabs (x2 - x3) > G_SQRT2 ||
          fabs (x3 - x0) > G_SQRT2 ||

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

/*
    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 (PixelSurround *surround,
            const gint     xc,
            const gint     yc,
            const gint     x0,
            const gint     y0,
            const gint     x1,
            const gint     y1,
            const gint     x2,
            const gint     y2,
            const gint     x3,
            const gint     y3,
            gint          *cc,
            const gint     level,
            guint         *color,
            const gint     bpp,
            const gint     alpha)
{
  if (!level || !supersample_test (x0, y0, x1, y1, x2, y2, x3, y3))
    {
      gint   i;
      guchar C[4];

      sample_bi (surround, xc, yc, C, 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 (surround,
                  tlx,tly,
                  x0,y0, tx,ty, xc,yc, lx,ly,
                  cc, level-1, color, bpp, alpha);

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

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

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

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

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

    get_sample (surround,
                DOUBLE2FIXED (xc), DOUBLE2FIXED (yc),
                DOUBLE2FIXED (x0), DOUBLE2FIXED (y0),
                DOUBLE2FIXED (x1), DOUBLE2FIXED (y1),
                DOUBLE2FIXED (x2), DOUBLE2FIXED (y2),
                DOUBLE2FIXED (x3), DOUBLE2FIXED (y3),
                &cc, level, C, 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 (const gdouble dx,
                               const gint    jm1,
                               const gint    j,
                               const gint    jp1,
                               const 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,
              const gdouble  u,
              const gdouble  v,
              guchar        *color,
              const gint     bytes,
              const gint     alpha)
{
  gdouble       a_val, a_recip;
  gint          i;
  const gint    iu = floor(u);
  const gint    iv = floor(v);
  gdouble       du, dv;
  gint          rowstride;
  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 =
        ROUND ((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);
    }
}

static void
sample_lanczos (PixelSurround *surround,
                const gfloat  *lanczos,
                const gdouble  u,
                const gdouble  v,
                guchar        *color,
                const gint     bytes,
                const gint     alpha)
{
  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       arecip;
  gdouble       aval;
  gint          su, sv;
  gint          b;
  gint          i, j;
  gint          iu, iv;
  gint          rowstride;
  const guchar *data;
  const guchar *src;

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

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

  /* fill 1D kernels */
  for (x_sum = y_sum = 0.0, i = LANCZOS_WIDTH; i >= -LANCZOS_WIDTH; i--)
    {
      gint 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;
    }

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

  src = data + alpha;
  aval = 0.0;

  for (j = 0; j < LANCZOS_WIDTH2 ; j++)
    {
      for (i = 0; i < LANCZOS_WIDTH2 ; i++)
        aval += y_kernel[j] * x_kernel[i] * (gdouble) src[i * bytes];

      src += rowstride;
    }

  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++)
    {
      const guchar *asrc;
      gdouble       newval = 0.0;

      src  = data + b;
      asrc = data + alpha;

      for (j = 0; j < LANCZOS_WIDTH2; j++)
        {
          for (i = 0; i < LANCZOS_WIDTH2; i++)
            newval += (y_kernel[j] * x_kernel[i] *
                       (gdouble) src[i * bytes] * (gdouble) asrc[i * bytes]);

          src += rowstride;
          asrc += rowstride;
        }

      newval *= arecip;
      color[b] = CLAMP (ROUND (newval), 0, 255);
    }

  color[alpha] = RINT (aval);
}