/*
    This file is part of darktable,
    Copyright (C) 2009-2026 darktable developers.

    darktable 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 3 of the License, or
    (at your option) any later version.

    darktable 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 darktable.  If not, see <http://www.gnu.org/licenses/>.
*/

#include "colorspace.h"
#include "color_conversion.h"
#include "common.h"
#include "rgb_norms.h"
#include "lut3d.h"

#include "diffuse.cl"

int
BL(const int row, const int col)
{
  return (((row & 1) << 1) + (col & 1));
}

kernel void
rawprepare_1f(read_only image2d_t in, write_only image2d_t out,
              const int width, const int height,
              const int cx, const int cy,
              global const float *sub, global const float *div, global const float *rdiv,
              const int rx, const int ry)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  const float pix = read_imageui(in, sampleri, (int2)(x + cx, y + cy)).x;
  const int id = BL(ry+cy+y, rx+cx+x);
  const float n = pix - sub[id];
  const float q = n * rdiv[id];
  const float rem = fma(-q, div[id], n);
  const float oval = fma(rem, rdiv[id], q);

  write_imagef(out, (int2)(x, y), oval);
}

// the gain at map_pt, interpolated as in rawprepare's process(). The
// linear sampler would do the same, but GPUs may compute its weights with
// as little as 8 bits of precision, which shows in the image
static float _gainmap_gain(read_only image2d_t map, const int2 map_size, const float2 map_pt)
{
  const float2 pt = clamp(map_pt, (float2)(0.0f, 0.0f), convert_float2(map_size));
  const int2 i0 = min(convert_int2(pt), map_size - 1);
  const int2 i1 = min(i0 + 1, map_size - 1);
  const float2 frac = pt - convert_float2(i0);
  const float top = (1.0f - frac.x) * read_imagef(map, sampleri, (int2)(i0.x, i0.y)).x
                    + frac.x * read_imagef(map, sampleri, (int2)(i1.x, i0.y)).x;
  const float bottom = (1.0f - frac.x) * read_imagef(map, sampleri, (int2)(i0.x, i1.y)).x
                       + frac.x * read_imagef(map, sampleri, (int2)(i1.x, i1.y)).x;
  return (1.0f - frac.y) * top + frac.y * bottom;
}

kernel void
rawprepare_1f_gainmap(read_only image2d_t in, write_only image2d_t out,
              const int width, const int height,
              const int cx, const int cy,
              global const float *sub, global const float *div, global const float *rdiv,
              const int rx, const int ry,
              read_only image2d_t map0, read_only image2d_t map1,
              read_only image2d_t map2, read_only image2d_t map3,
              const int2 map_size, const float2 im_to_rel,
              const float2 rel_to_map, const float2 map_origin)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  const float pix = read_imageui(in, sampleri, (int2)(x + cx, y + cy)).x;
  const int id = BL(ry+cy+y, rx+cx+x);
  const float n = pix - sub[id];
  const float q = n * rdiv[id];
  const float rem = fma(-q, div[id], n);
  float pixel_scaled = fma(rem, rdiv[id], q);

  const float2 map_pt = ((float2)(rx+cx+x,ry+cy+y) * im_to_rel - map_origin) * rel_to_map;
  switch(id)
  {
    case 0:
      pixel_scaled *= _gainmap_gain(map0, map_size, map_pt);
      break;
    case 1:
      pixel_scaled *= _gainmap_gain(map1, map_size, map_pt);
      break;
    case 2:
      pixel_scaled *= _gainmap_gain(map2, map_size, map_pt);
      break;
    case 3:
      pixel_scaled *= _gainmap_gain(map3, map_size, map_pt);
      break;
  }

  write_imagef(out, (int2)(x, y), pixel_scaled);
}

kernel void
rawprepare_1f_unnormalized(read_only image2d_t in, write_only image2d_t out,
                           const int width, const int height,
                           const int cx, const int cy,
                           global const float *sub, global const float *div, global const float *rdiv,
                           const int rx, const int ry)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width  || y >= height) return;

  const float pix = read_imagef(in, sampleri, (int2)(x + cx, y + cy)).x;
  const int id = BL(ry+cy+y, rx+cx+x);
  const float n = pix - sub[id];
  const float q = n * rdiv[id];
  const float rem = fma(-q, div[id], n);
  const float oval = fma(rem, rdiv[id], q);

  write_imagef(out, (int2)(x, y), oval);
}

kernel void
rawprepare_1f_unnormalized_gainmap(read_only image2d_t in, write_only image2d_t out,
                           const int width, const int height,
                           const int cx, const int cy,
                           global const float *sub, global const float *div, global const float *rdiv,
                           const int rx, const int ry,
                           read_only image2d_t map0, read_only image2d_t map1,
                           read_only image2d_t map2, read_only image2d_t map3,
                           const int2 map_size, const float2 im_to_rel,
                           const float2 rel_to_map, const float2 map_origin)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width  || y >= height) return;

  const float pix = read_imagef(in, sampleri, (int2)(x + cx, y + cy)).x;
  const int id = BL(ry+cy+y, rx+cx+x);
  const float n = pix - sub[id];
  const float q = n * rdiv[id];
  const float rem = fma(-q, div[id], n);
  float pixel_scaled = fma(rem, rdiv[id], q);

  const float2 map_pt = ((float2)(rx+cx+x,ry+cy+y) * im_to_rel - map_origin) * rel_to_map;
  switch(id)
  {
    case 0:
      pixel_scaled *= _gainmap_gain(map0, map_size, map_pt);
      break;
    case 1:
      pixel_scaled *= _gainmap_gain(map1, map_size, map_pt);
      break;
    case 2:
      pixel_scaled *= _gainmap_gain(map2, map_size, map_pt);
      break;
    case 3:
      pixel_scaled *= _gainmap_gain(map3, map_size, map_pt);
      break;
  }

  write_imagef(out, (int2)(x, y), pixel_scaled);
}

kernel void
rawprepare_4f(read_only image2d_t in, write_only image2d_t out,
              const int width, const int height,
              const int cx, const int cy,
              global const float *sub, global const float *div, global const float *rdiv,
              const int rx, const int ry)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  const float4 sub4 = (const float4)(sub[0], sub[1], sub[2], sub[3]);
  const float4 div4 = (const float4)(div[0], div[1], div[2], div[3]);
  const float4 rdiv4 = (const float4)(rdiv[0], rdiv[1], rdiv[2], rdiv[3]);
  const float4 pix = readpixel(in, x + cx, y + cy);
  const float4 n = pix - sub4;
  const float4 q = n * rdiv4;
  const float4 rem = fma(-q, div4, n);
  float4 oval = fma(rem, rdiv4, q);
  oval.w = 0.0f;
  write_imagef(out, (int2)(x, y), oval);
}

kernel void
invert_1f(read_only image2d_t in, write_only image2d_t out, const int width, const int height, global float *color,
          const unsigned int filters, const int rx, const int ry)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);
  if(x >= width || y >= height) return;
  const float pixel = read_imagef(in, sampleri, (int2)(x, y)).x;
  const float inv_pixel = color[FC(ry+y, rx+x, filters)] - pixel;

  write_imagef (out, (int2)(x, y), (float4)(clipf(inv_pixel), 0.0f, 0.0f, 0.0f));
}

kernel void
invert_4f(read_only image2d_t in, write_only image2d_t out, const int width, const int height, global float *color,
                const unsigned int filters, const int rx, const int ry)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);
  if(x >= width || y >= height) return;
  float4 pixel = read_imagef(in, sampleri, (int2)(x, y));
  pixel.x = color[0] - pixel.x;
  pixel.y = color[1] - pixel.y;
  pixel.z = color[2] - pixel.z;
  pixel.xyz = clamp(pixel.xyz, 0.0f, 1.0f);

  write_imagef (out, (int2)(x, y), pixel);
}

kernel void
whitebalance_1f(read_only image2d_t in, write_only image2d_t out, const int width, const int height, global float *coeffs,
    const unsigned int filters, global const unsigned char (*const xtrans)[6])
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);
  if(x >= width || y >= height) return;
  float pixel = readsingle(in, x, y);
  write_imagef(out, (int2)(x, y), pixel * coeffs[FC(y, x, filters)]);
}

kernel void
whitebalance_1f_xtrans(read_only image2d_t in, write_only image2d_t out, const int width, const int height, global float *coeffs,
    const unsigned int filters, global const unsigned char (*const xtrans)[6])
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);
  if(x >= width || y >= height) return;
  float pixel = readsingle(in, x, y);
  write_imagef(out, (int2)(x, y), pixel * coeffs[FCxtrans(y, x, xtrans)]);
}


kernel void
whitebalance_4f(read_only image2d_t in, write_only image2d_t out, const int width, const int height, global float *coeffs,
    const unsigned int filters, global const unsigned char (*const xtrans)[6])
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);
  if(x >= width || y >= height) return;
  float4 pixel = Areadpixel(in, x, y);
  write_imagef (out, (int2)(x, y), (float4)(pixel.x * coeffs[0], pixel.y * coeffs[1], pixel.z * coeffs[2], pixel.w));
}

/* kernel for the exposure plugin. should work transparently with float4 and float image2d. */
kernel void
exposure (read_only image2d_t in, write_only image2d_t out, const int width, const int height, const float black, const float scale)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;
  float4 pixel = Areadpixel(in, x, y);
  pixel.xyz = ((pixel - black ) * scale).xyz;
  write_imagef (out, (int2)(x, y), pixel);
}

/* kernel for the highlights plugin. */
kernel void
highlights_4f_clip (read_only image2d_t in, write_only image2d_t out, const int width, const int height, const float clip)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  // 4f/pixel means that this has been debayered already.
  // it's thus hopeless to recover highlights here (this code path is just used for preview and non-raw images)
  float4 pixel = Areadpixel(in, x, y);
  // default: // 0, DT_IOP_HIGHLIGHTS_CLIP
  pixel.x = fmin(clip, pixel.x);
  pixel.y = fmin(clip, pixel.y);
  pixel.z = fmin(clip, pixel.z);
  write_imagef (out, (int2)(x, y), pixel);
}

kernel void
highlights_1f_clip (read_only image2d_t in, write_only image2d_t out,
                    const int iwidth, const int iheight,
                    const int owidth, const int oheight,
                    const float4 clip, const int dx, const int dy,
                    const int filters, global const unsigned char (*const xtrans)[6])
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= owidth || y >= oheight) return;

  const int irow = y + dy;
  const int icol = x + dx;
  float pixel = 0.0f;
  float clips[4] = { clip.x, clip.y, clip.z, clip.w};
  if((icol >= 0) && (irow >= 0) && (irow < iheight) && (icol < iwidth))
  {
    const int color = fcol(irow, icol, filters, xtrans);
    pixel = Areadsingle(in, icol, irow);
    pixel = fmin(clips[color], pixel);
  }
  write_imagef(out, (int2)(x, y), pixel);
}

kernel void highlights_false_color(read_only image2d_t in,
                                   write_only image2d_t out,
                                   const int owidth,
                                   const int oheight,
                                   const int iwidth,
                                   const int iheight,
                                   const int dx,
                                   const int dy,
                                   const unsigned int filters,
                                   global const unsigned char (*const xtrans)[6],
                                   const float4 clip)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= owidth || y >= oheight) return;

  const int irow = y + dy;
  const int icol = x + dx;
  float oval = 0.0f;
  float clips[4] = { clip.x, clip.y, clip.z, clip.w};

  if((irow >= 0) && (icol >= 0) && (icol < iwidth) && (irow < iheight))
  {
    const float ival = Areadsingle(in, icol, irow);
    const int c = fcol(irow, icol, filters, xtrans);
    oval = (ival < clips[c]) ? 0.2f * ival : 1.0f;
  }
  write_imagef (out, (int2)(x, y), oval);
}

static float _calc_refavg(read_only image2d_t in,
                          global const unsigned char (*const xtrans)[6],
                          const unsigned int filters,
                          int row,
                          int col,
                          int maxrow,
                          int maxcol,
                          global const float *correction)
{
  float mean[3] = { 0.0f, 0.0f, 0.0f };
  int cnt[3]  = { 0, 0, 0 };

  for(int dy = max(0, row - 1); dy < min(maxrow - 1, row + 2); dy++)
  {
    for(int dx = max(0, col - 1); dx < min(maxcol - 1, col + 2); dx++)
    {
      const float val = fmax(0.0f, Areadsingle(in, dx, dy));
      const int c = fcol(dy, dx, filters, xtrans);
      mean[c] += val;
      cnt[c] += 1;
    }
  }

  for(int c = 0; c < 3; c++)
    mean[c] = cnt[c] > 0 ? cbrt((correction[c] * mean[c]) / (float)cnt[c]) : 0.0f;

  float croot_refavg[3] = { 0.5f * (mean[1] + mean[2]), 0.5f * (mean[0] + mean[2]), 0.5f * (mean[0] + mean[1])};
  int color = fcol(row, col, filters, xtrans);
  return fcube(croot_refavg[color]);
}

kernel void highlights_initmask(read_only image2d_t in,
                                global char *inmask,
                                const int msize,
                                const int mwidth,
                                const int mheight,
                                const int iwidth,
                                const int iheight,
                                const unsigned int filters,
                                global const unsigned char (*const xtrans)[6],
                                global const float *clips)
{
  const int mcol = get_global_id(0);
  const int mrow = get_global_id(1);

  if(mcol >= mwidth || mrow >= mheight)
    return;

  int mdx = mad24(mrow, mwidth, mcol);
  char mbuff[3] = { 0, 0, 0 };
  int irow = 3*mrow;
  int icol = 3*mcol;

  for(int y = irow; y < irow + 3; y++)
  {
    for(int x = icol; x < icol + 3; x++)
    {
      if(y < iheight && x < iwidth)
      {
        int color = fcol(y, x, filters, xtrans);
        float val = Areadsingle(in, x, y);
        mbuff[color] += (val >= clips[color]) ? 1 : 0;
      }
    }
  }

  for(int c = 0; c < 3; c++)
    inmask[c*msize + mdx] = mbuff[c] ? 1 : 0;
}

static char _mask_dilated(global char *in, const int w1)
{
  if(in[0])
    return 1;

  if(in[-w1-1] | in[-w1] | in[-w1+1] | in[-1] | in[1] | in[w1-1] | in[w1] | in[w1+1])
    return 1;

  const int w2 = 2*w1;
  const int w3 = 3*w1;
  return (in[-w3-2] | in[-w3-1] | in[-w3]   | in[-w3+1] | in[-w3+2] |
          in[-w2-3] | in[-w2-2] | in[-w2-1] | in[-w2]   | in[-w2+1] | in[-w2+2] | in[-w2+3] |
          in[-w1-3] | in[-w1-2] | in[-w1+2] | in[-w1+3] |
          in[-3]    | in[-2]    | in[2]     | in[3]     |
          in[w1-3]  | in[w1-2]  | in[w1+2]  | in[w1+3]  |
          in[w2-3]  | in[w2-2]  | in[w2-1]  | in[w2]    | in[w2+1]  | in[w2+2]  | in[w2+3] |
          in[w3-2]  | in[w3-1]  | in[w3]    | in[w3+1]  | in[w3+2]) ? 1 : 0;
}

kernel void highlights_dilatemask(global char *in,
                                  global char *out,
                                  const int msize,
                                  const int mwidth,
                                  const int mheight)
{
  const int col = get_global_id(0);
  const int row = get_global_id(1);

  if(col >= mwidth || row >= mheight)
    return;

  int mx = mad24(row, mwidth, col);
  bool safe = col >= 3 && row >= 3 && col < mwidth - 3 && row < mheight - 3;

  out[mx] =         safe ? _mask_dilated(in + mx, mwidth)           : in[mx];
  out[mx+msize] =   safe ? _mask_dilated(in + mx + msize, mwidth)   : in[mx + msize];
  out[mx+2*msize] = safe ? _mask_dilated(in + mx + 2*msize, mwidth) : in[mx + 2*msize];
}

kernel void highlights_opposed(read_only image2d_t in,
                               write_only image2d_t out,
                               const int owidth,
                               const int oheight,
                               const int iwidth,
                               const int iheight,
                               const int dx,
                               const int dy,
                               const unsigned int filters,
                               global const unsigned char (*const xtrans)[6],
                               global const float *clips,
                               global const float *chroma,
                               global const float *correction)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);
  if(x >= owidth || y >= oheight) return;

  int irow = y + dy;
  int icol = x + dx;
  float val = 0.0f;

  if((icol >= 0) && (icol < iwidth) && (irow >= 0) && (irow < iheight))
  {
    val = Areadsingle(in, icol, irow);
    const int color = fcol(irow, icol, filters, xtrans);
    if(val >= clips[color])
    {
      const float ref = _calc_refavg(in, xtrans, filters, irow, icol, iheight, iwidth, correction);
      val = fmax(val, ref + chroma[color]);
    }
  }
  write_imagef(out, (int2)(x, y), val);
}

// linear raws, see _process_linear_opposed() in src/iop/hlreconstruct/opposed.c
static float _calc_linear_refavg(const float4 pixel, const int color)
{
  const float ins[3] = { cbrt(fmax(0.0f, pixel.x)), cbrt(fmax(0.0f, pixel.y)), cbrt(fmax(0.0f, pixel.z)) };
  const float opp[3] = { 0.5f*(ins[1]+ins[2]), 0.5f*(ins[0]+ins[2]), 0.5f*(ins[0]+ins[1]) };
  return fcube(opp[color]);
}

// a mask cell covers 3x3 pixels; the last row and column of cells also take
// the pixels left over by the division, as _raw_to_cmap() clamps to them
kernel void highlights_initmask_linear(read_only image2d_t in,
                                       global char *inmask,
                                       const int msize,
                                       const int mwidth,
                                       const int mheight,
                                       const int iwidth,
                                       const int iheight,
                                       const float4 clips)
{
  const int mcol = get_global_id(0);
  const int mrow = get_global_id(1);

  if(mcol >= mwidth || mrow >= mheight)
    return;

  const int irow = 3*mrow;
  const int icol = 3*mcol;
  const int rowend = mrow == mheight - 1 ? iheight : irow + 3;
  const int colend = mcol == mwidth - 1 ? iwidth : icol + 3;
  char mbuff[3] = { 0, 0, 0 };

  for(int y = irow; y < rowend; y++)
  {
    for(int x = icol; x < colend; x++)
    {
      const float4 val = readpixel(in, x, y);
      mbuff[0] |= val.x >= clips.x;
      mbuff[1] |= val.y >= clips.y;
      mbuff[2] |= val.z >= clips.z;
    }
  }

  const int mdx = mad24(mrow, mwidth, mcol);
  for(int c = 0; c < 3; c++)
    inmask[c*msize + mdx] = mbuff[c] ? 1 : 0;
}

kernel void highlights_opposed_linear(read_only image2d_t in,
                                      write_only image2d_t out,
                                      const int width,
                                      const int height,
                                      const float4 clips,
                                      const float4 chroma)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);
  if(x >= width || y >= height) return;

  const float4 pixel = readpixel(in, x, y);
  const float clip[3] = { clips.x, clips.y, clips.z };
  const float chr[3] = { chroma.x, chroma.y, chroma.z };
  const float inval[3] = { fmax(0.0f, pixel.x), fmax(0.0f, pixel.y), fmax(0.0f, pixel.z) };
  float oval[3];
  for(int c = 0; c < 3; c++)
    oval[c] = inval[c] >= clip[c] ? fmax(inval[c], _calc_linear_refavg(pixel, c) + chr[c]) : inval[c];

  write_imagef(out, (int2)(x, y), (float4)(oval[0], oval[1], oval[2], pixel.w));
}

#define SQRT3 1.7320508075688772935274463415058723669f
#define SQRT12 3.4641016151377545870548926830117447339f // 2*SQRT3
kernel void
highlights_1f_lch_bayer (read_only image2d_t in, write_only image2d_t out, const int width, const int height,
                         const float clip, const int filters)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  int clipped = 0;
  float R = 0.0f;
  float Gmin = FLT_MAX;
  float Gmax = -FLT_MAX;
  float B = 0.0f;
  float pixel = 0.0f;

  // sample 1 bayer block. thus we will have 2 green values.
  for(int jj = 0; jj <= 1; jj++)
  {
    for(int ii = 0; ii <= 1; ii++)
    {
      const float val = readsingle(in, x+ii, y+jj);

      pixel = (ii == 0 && jj == 0) ? val : pixel;

      clipped = (clipped || (val > clip));

      const int c = FC(y + jj, x + ii, filters);

      switch(c)
      {
        case 0:
          R = val;
          break;
        case 1:
          Gmin = fmin(Gmin, val);
          Gmax = fmax(Gmax, val);
          break;
        case 2:
          B = val;
          break;
      }
    }
  }

  if(clipped)
  {
    const float Ro = fmin(R, clip);
    const float Go = fmin(Gmin, clip);
    const float Bo = fmin(B, clip);

    const float L = (R + Gmax + B) / 3.0f;

    float C = SQRT3 * (R - Gmax);
    float H = 2.0f * B - Gmax - R;

    const float Co = SQRT3 * (Ro - Go);
    const float Ho = 2.0f * Bo - Go - Ro;

    const float ratio = (R != Gmax && Gmax != B) ? sqrt((Co * Co + Ho * Ho) / (C * C + H * H)) : 1.0f;

    C *= ratio;
    H *= ratio;

    /*
     * backtransform proof, sage:
     *
     * R,G,B,L,C,H = var('R,G,B,L,C,H')
     * solve([L==(R+G+B)/3, C==sqrt(3)*(R-G), H==2*B-G-R], R, G, B)
     *
     * result:
     * [[R == 1/6*sqrt(3)*C - 1/6*H + L, G == -1/6*sqrt(3)*C - 1/6*H + L, B == 1/3*H + L]]
     */
    const int c = FC(y, x, filters);
    C = (c == 1) ? -C : C;

    pixel = L;
    pixel += (c == 2) ? H / 3.0f : -H / 6.0f + C / SQRT12;
  }

  write_imagef (out, (int2)(x, y), pixel);
}


kernel void
highlights_1f_lch_xtrans (read_only image2d_t in, write_only image2d_t out, const int width, const int height,
                         const float clip, global const unsigned char (*const xtrans)[6],
                         local float *buffer)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);
  const int xlsz = get_local_size(0);
  const int ylsz = get_local_size(1);
  const int xlid = get_local_id(0);
  const int ylid = get_local_id(1);
  const int xgid = get_group_id(0);
  const int ygid = get_group_id(1);

  // individual control variable in this work group and the work group size
  const int l = mad24(ylid, xlsz, xlid);
  const int lsz = mul24(xlsz, ylsz);

  // stride and maximum capacity of local buffer
  // cells of 1*float per pixel with a surrounding border of 2 cells
  const int stride = xlsz + 2*2;
  const int maxbuf = mul24(stride, ylsz + 2*2);

  // coordinates of top left pixel of buffer
  // this is 2 pixel left and above of the work group origin
  const int xul = mul24(xgid, xlsz) - 2;
  const int yul = mul24(ygid, ylsz) - 2;

  // populate local memory buffer
  for(int n = 0; n <= maxbuf/lsz; n++)
  {
    const int bufidx = mad24(n, lsz, l);
    if(bufidx >= maxbuf) continue;
    const int xx = xul + bufidx % stride;
    const int yy = yul + bufidx / stride;
    buffer[bufidx] = readsingle(in, xx, yy);
  }

  // center buffer around current x,y-Pixel
  buffer += mad24(ylid + 2, stride, xlid + 2);

  barrier(CLK_LOCAL_MEM_FENCE);

  if(x >= width || y >= height) return;

  float pixel = 0.0f;

  if(x < 2 || x > width - 3 || y < 2 || y > height - 3)
  {
    // fast path for border
    pixel = fmin(clip, buffer[0]);
  }
  else
  {
    // if current pixel is clipped, always reconstruct
    int clipped = (buffer[0] > clip);

    if(!clipped)
    {
      clipped = 1;
      // check if there is any 3x3 block touching the current
      // pixel which has no clipping, as then we don't need to
      // reconstruct the current pixel. This avoids zippering in
      // edge transitions from clipped to unclipped areas. The
      // X-Trans sensor seems prone to this, unlike Bayer, due
      // to its irregular pattern.
      for(int offset_j = -2; offset_j <= 0; offset_j++)
      {
        for(int offset_i = -2; offset_i <= 0; offset_i++)
        {
          if(clipped)
          {
            clipped = 0;
            for(int jj = offset_j; jj <= offset_j + 2; jj++)
            {
              for(int ii = offset_i; ii <= offset_i + 2; ii++)
              {
                const float val = buffer[mad24(jj, stride, ii)];
                clipped = (clipped || (val > clip));
              }
            }
          }
        }
      }
    }

    if(clipped)
    {
      float mean[3] = { 0.0f, 0.0f, 0.0f };
      int cnt[3] = { 0, 0, 0 };
      float RGBmax[3] = { -FLT_MAX, -FLT_MAX, -FLT_MAX };

      for(int jj = -1; jj <= 1; jj++)
      {
        for(int ii = -1; ii <= 1; ii++)
        {
          const float val = buffer[mad24(jj, stride, ii)];
          const int c = FCxtrans(y + jj, x + ii, xtrans);
          mean[c] += val;
          cnt[c]++;
          RGBmax[c] = fmax(RGBmax[c], val);
        }
      }

      const float Ro = fmin(mean[0]/cnt[0], clip);
      const float Go = fmin(mean[1]/cnt[1], clip);
      const float Bo = fmin(mean[2]/cnt[2], clip);

      const float R = RGBmax[0];
      const float G = RGBmax[1];
      const float B = RGBmax[2];

      const float L = (R + G + B) / 3.0f;
      float C = SQRT3 * (R - G);
      float H = 2.0f * B - G - R;

      const float Co = SQRT3 * (Ro - Go);
      const float Ho = 2.0f * Bo - Go - Ro;

      if(R != G && G != B)
      {
        const float ratio = sqrt((Co * Co + Ho * Ho) / (C * C + H * H));
        C *= ratio;
        H *= ratio;
      }

      float RGB[3] = { 0.0f, 0.0f, 0.0f };

      RGB[0] = L - H / 6.0f + C / SQRT12;
      RGB[1] = L - H / 6.0f - C / SQRT12;
      RGB[2] = L + H / 3.0f;

      pixel = RGB[FCxtrans(y, x, xtrans)];
    }
    else
      pixel = buffer[0];
  }

  write_imagef (out, (int2)(x, y), pixel);
}
#undef SQRT3
#undef SQRT12


kernel void
interpolate_and_mask(read_only image2d_t input,
                     write_only image2d_t interpolated,
                     write_only image2d_t clipping_mask,
                     global const float *clips,
                     global const float *wb,
                     const unsigned int filters,
                     const int width,
                     const int height)
{
  // Bilinear interpolation
  const int j = get_global_id(0); // = x
  const int i = get_global_id(1); // = y

  if(j >= width || i >= height) return;

  /** On current mesa we have to use one sampler in one kernel!! */
  const float center = fmax(0.0f, Areadsingle(input, j, i));

  const int c = FC(i, j, filters);

  float R = 0.f;
  float G = 0.f;
  float B = 0.f;

  int R_clipped = 0;
  int G_clipped = 0;
  int B_clipped = 0;

  if(i == 0 || j == 0 || i == height - 1 || j == width - 1)
  {
    // We are on the image edges. We don't need to demosaic,
    // just set R = G = B = center and record clipping.
    // This will introduce a marginal error close to edges, mostly irrelevant
    // because we are dealing with local averages anyway, later on.
    // Also we remosaic the image at the end, so only the relevant channel gets picked.
    // Finally, it's unlikely that the borders of the image get clipped due to vignetting.
    R = G = B = center;
    R_clipped = G_clipped = B_clipped = (center > clips[c]);
  }
  else
  {
    // fetch neighbours and cache them for perf
    const size_t i_prev = (i - 1);
    const size_t i_next = (i + 1);
    const size_t j_prev = (j - 1);
    const size_t j_next = (j + 1);

    const float north = Areadsingle(input, j, i_prev);
    const float south = Areadsingle(input, j, i_next);
    const float west = Areadsingle(input, j_prev, i);
    const float east = Areadsingle(input, j_next, i);

    const float north_east = Areadsingle(input, j_next, i_prev);
    const float north_west = Areadsingle(input, j_prev, i_prev);
    const float south_east = Areadsingle(input, j_next, i_next);
    const float south_west = Areadsingle(input, j_prev, i_next);

    if(c == GREEN) // green pixel
    {
      G = center;
      G_clipped = (center > clips[GREEN]);
    }
    else // non-green pixel
    {
      // interpolate inside an X/Y cross
      G = (north + south + east + west) / 4.f;
      G_clipped = (north > clips[GREEN] || south > clips[GREEN] || east > clips[GREEN] || west > clips[GREEN]);
    }

    if(c == RED ) // red pixel
    {
      R = center;
      R_clipped = (center > clips[RED]);
    }
    else // non-red pixel
    {
      if(FC(i - 1, j, filters) == RED && FC(i + 1, j, filters) == RED)
      {
        // we are on a red column, so interpolate column-wise
        R = (north + south) / 2.f;
        R_clipped = (north > clips[RED] || south > clips[RED]);
      }
      else if(FC(i, j - 1, filters) == RED && FC(i, j + 1, filters) == RED)
      {
        // we are on a red row, so interpolate row-wise
        R = (west + east) / 2.f;
        R_clipped = (west > clips[RED] || east > clips[RED]);
      }
      else
      {
        // we are on a blue row, so interpolate inside a square
        R = (north_west + north_east + south_east + south_west) / 4.f;
        R_clipped = (north_west > clips[RED] || north_east > clips[RED] || south_west > clips[RED]
                      || south_east > clips[RED]);
      }
    }

    if(c == BLUE ) // blue pixel
    {
      B = center;
      B_clipped = (center > clips[BLUE]);
    }
    else // non-blue pixel
    {
      if(FC(i - 1, j, filters) == BLUE && FC(i + 1, j, filters) == BLUE)
      {
        // we are on a blue column, so interpolate column-wise
        B = (north + south) / 2.f;
        B_clipped = (north > clips[BLUE] || south > clips[BLUE]);
      }
      else if(FC(i, j - 1, filters) == BLUE && FC(i, j + 1, filters) == BLUE)
      {
        // we are on a red row, so interpolate row-wise
        B = (west + east) / 2.f;
        B_clipped = (west > clips[BLUE] || east > clips[BLUE]);
      }
      else
      {
        // we are on a red row, so interpolate inside a square
        B = (north_west + north_east + south_east + south_west) / 4.f;

        B_clipped = (north_west > clips[BLUE] || north_east > clips[BLUE] || south_west > clips[BLUE]
                    || south_east > clips[BLUE]);
      }
    }
  }

  const float4 RGB = {R, G, B, dtcl_sqrt(R * R + G * G + B * B) };
  const float4 clipped = { R_clipped ? 1.0f : 0.0f, G_clipped ? 1.0f : 0.0f, B_clipped ? 1.0f : 0.0f, (R_clipped || G_clipped || B_clipped) ? 1.0f : 0.0f };
  const float4 WB4 = { wb[0], wb[1], wb[2], 1.0f };
  write_imagef(interpolated, (int2)(j, i), RGB / WB4);
  write_imagef(clipping_mask, (int2)(j, i), clipped);
}


kernel void
remosaic_and_replace(read_only image2d_t input,
                     read_only image2d_t interpolated,
                     read_only image2d_t clipping_mask,
                     write_only image2d_t output,
                     global const float *wb,
                     const unsigned int filters,
                     const int width,
                     const int height)
{
  // Take RGB ratios and norm, reconstruct RGB and remosaic the image
  const int j = get_global_id(0); // = x
  const int i = get_global_id(1); // = y

  if(j >= width || i >= height) return;

  const int c = FC(i, j, filters);
  const float4 center = Areadpixel(interpolated, j, i);
  const float4 clipmask = Areadpixel(clipping_mask, j, i);
  const float opacity = isfinite(clipmask.w) ? clipf(clipmask.w) : 0.0f;
  const float pix_in = readsingle(input, j, i);

  const float corrections[4] = { center.x * wb[0], center.y * wb[1], center.z * wb[2], 0.0f};
  const float correction = fmax(corrections[c], 0.0f);
  write_imagef(output, (int2)(j, i), (float4)(mix(pix_in, correction, opacity), 0.0f, 0.0f, 0.0f));
}

kernel void
box_blur_5x5(read_only image2d_t in,
             write_only image2d_t out,
             const int width,
             const int height)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  float4 acc = 0.f;

  for(int ii = -2; ii < 3; ++ii)
  {
    for(int jj = -2; jj < 3; ++jj)
    {
      const int row = clamp(y + ii, 0, height - 1);
      const int col = clamp(x + jj, 0, width - 1);
      acc += fmax(0.0f, Areadpixel(in, col, row));
    }
  }
  acc *= 0.04f;
  write_imagef(out, (int2)(x, y), acc);
}


enum wavelets_scale_t
{
  ANY_SCALE   = 1 << 0, // any wavelets scale   : reconstruct += HF
  FIRST_SCALE = 1 << 1, // first wavelets scale : reconstruct = 0
  LAST_SCALE  = 1 << 2, // last wavelets scale  : reconstruct += residual
};


kernel void
guide_laplacians(read_only image2d_t HF,
                 read_only image2d_t LF,
                 read_only image2d_t mask,
                 read_only image2d_t output_r,
                 write_only image2d_t output_w,
                 const int width,
                 const int height,
                 const int mult,
                 const float noise_level,
                 const int salt,
                 const unsigned int scale,
                 const float radius_sq)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  const float alpha = read_imagef(mask, sampleri, (int2)(x, y)).w;
  const float alpha_comp = 1.f - alpha;

  float4 high_frequency = readpixel(HF, x, y);

  float4 out;

  if(alpha > 0.f) // reconstruct
  {
    // non-local neighbours coordinates
    const int j_neighbours[3] = { max(x - mult, 0), x, min(x + mult, width - 1) };
    const int i_neighbours[3] = { max(y - mult, 0), y, min(y + mult, height - 1) };

    // fetch non-local pixels and store them locally and contiguously
    float4 neighbour_pixel_HF[9];
    neighbour_pixel_HF[3 * 0 + 0] = read_imagef(HF, samplerA, (int2)(j_neighbours[0], i_neighbours[0]));
    neighbour_pixel_HF[3 * 0 + 1] = read_imagef(HF, samplerA, (int2)(j_neighbours[1], i_neighbours[0]));
    neighbour_pixel_HF[3 * 0 + 2] = read_imagef(HF, samplerA, (int2)(j_neighbours[2], i_neighbours[0]));

    neighbour_pixel_HF[3 * 1 + 0] = read_imagef(HF, samplerA, (int2)(j_neighbours[0], i_neighbours[1]));
    neighbour_pixel_HF[3 * 1 + 1] = read_imagef(HF, samplerA, (int2)(j_neighbours[1], i_neighbours[1]));
    neighbour_pixel_HF[3 * 1 + 2] = read_imagef(HF, samplerA, (int2)(j_neighbours[2], i_neighbours[1]));

    neighbour_pixel_HF[3 * 2 + 0] = read_imagef(HF, samplerA, (int2)(j_neighbours[0], i_neighbours[2]));
    neighbour_pixel_HF[3 * 2 + 1] = read_imagef(HF, samplerA, (int2)(j_neighbours[1], i_neighbours[2]));
    neighbour_pixel_HF[3 * 2 + 2] = read_imagef(HF, samplerA, (int2)(j_neighbours[2], i_neighbours[2]));

    // Compute the linear fit of the laplacian of chromaticity against the laplacian of the norm
    // that is the chromaticity filter guided by the norm

    // Get the local average per channel
    const float ninth = 1.0f / 9.0f;
    float4 means_HF = 0.f;
    for(int k = 0; k < 9; k++)
    {
      means_HF += neighbour_pixel_HF[k];
    }
    means_HF *= ninth;

    // Get the local variance per channel
    float4 variance_HF = 0.f;
    for(int k = 0; k < 9; k++)
    {
      variance_HF += sqf(neighbour_pixel_HF[k] - means_HF);
    }
    variance_HF *= ninth;

    // Find the channel most likely to contain details = max( variance(HF) )
    // But since OpenCL is not designed to iterate over float4,
    // we need to check each channel in sequence
    int guiding_channel_HF = ALPHA;
    float guiding_value_HF = 0.f;

    if(variance_HF.x > guiding_value_HF)
    {
      guiding_value_HF = variance_HF.x;
      guiding_channel_HF = RED;
    }
    if(variance_HF.y > guiding_value_HF)
    {
      guiding_value_HF = variance_HF.y;
      guiding_channel_HF = GREEN;
    }
    if(variance_HF.z > guiding_value_HF)
    {
      guiding_value_HF = variance_HF.z;
      guiding_channel_HF = BLUE;
    }

    // Extract the guiding values for HF and LF now
    // so we can proceed after with vectorized code
    float means_HF_guide = 0.f;
    float variance_HF_guide = 0.f;
    float channel_guide_HF[9];
    float high_frequency_guide = 0.f;

    if(guiding_channel_HF == RED)
    {
      means_HF_guide = means_HF.x;
      variance_HF_guide = variance_HF.x;
      high_frequency_guide = high_frequency.x;
      for(int k = 0; k < 9; k++) channel_guide_HF[k] = neighbour_pixel_HF[k].x;
    }
    else if(guiding_channel_HF == GREEN)
    {
      means_HF_guide = means_HF.y;
      variance_HF_guide = variance_HF.y;
      high_frequency_guide = high_frequency.y;
      for(int k = 0; k < 9; k++) channel_guide_HF[k] = neighbour_pixel_HF[k].y;
    }
    else // BLUE
    {
      means_HF_guide = means_HF.z;
      variance_HF_guide = variance_HF.z;
      high_frequency_guide = high_frequency.z;
      for(int k = 0; k < 9; k++) channel_guide_HF[k] = neighbour_pixel_HF[k].z;
    }

    // Compute the linear regression channel = f(guide)
    float4 covariance_HF = 0.f;
    for(int k = 0; k < 9; k++)
    {
      covariance_HF += (neighbour_pixel_HF[k] - means_HF)
                       * (channel_guide_HF[k] - means_HF_guide);
    }
    covariance_HF *= ninth;
    const float scale_multiplier = 1.f / radius_sq;
    const float4 alpha_ch = readpixel(mask, x, y);

    const float4 a_HF = fmax(covariance_HF / variance_HF_guide, 0.f);
    const float4 b_HF = means_HF - a_HF * means_HF_guide;

    // Guide all channels by the norms
    high_frequency = alpha_ch * scale_multiplier * (a_HF * high_frequency_guide + b_HF)
                   + (1.f - alpha_ch * scale_multiplier) * high_frequency;

  }

  if(scale & FIRST_SCALE)
  {
    // out is not inited yet
    out = high_frequency;
  }
  else
  {
    // just accumulate HF
    out = read_imagef(output_r, sampleri, (int2)(x, y)) + high_frequency;
  }

  if(scale & LAST_SCALE)
  {
    // add the residual and clamp
    out = fmax(out + readpixel(LF, x, y), (float4)0.f);
  }

  // Last step of RGB reconstruct : add noise
  if((scale & LAST_SCALE) && salt && alpha > 0.f)
  {
    // Init random number generator
    unsigned int state[4] = { splitmix32(x + 1), splitmix32((x + 1) * (y + 3)), splitmix32(1337), splitmix32(666) };
    xoshiro128plus(state);
    xoshiro128plus(state);
    xoshiro128plus(state);
    xoshiro128plus(state);

    // Model noise on the max RGB
    const float4 sigma = out * noise_level;
    float4 noise = dt_noise_generator_simd(DT_NOISE_POISSONIAN, out, sigma, state);

    // Ensure the noise only brightens the image, since it's clipped
    noise = out + fabs(noise - out);
    out = fmax(alpha * noise + alpha_comp * out, 0.f);
  }

  if(scale & LAST_SCALE)
  {
    // Break the RGB channels into ratios/norm for the next step of reconstruction
    const float4 out_2 = out * out;
    const float norm = fmax(sqrt(out_2.x + out_2.y + out_2.z), 1e-6f);
    out /= norm;
    out.w = norm;
  }

  write_imagef(output_w, (int2)(x, y), out);
}

kernel void
diffuse_color(read_only image2d_t HF,
              read_only image2d_t LF,
              read_only image2d_t mask,
              read_only image2d_t output_r,
              write_only image2d_t output_w,
              const int width,
              const int height,
              const int mult,
              const unsigned int scale,
              const float first_order_factor)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  const float4 alpha = readpixel(mask, x, y);
  float4 high_frequency = readpixel(HF, x, y);

  const float norm_backup = high_frequency.w;

  float4 out;

  if(alpha.w > 0.f) // reconstruct
  {
    // non-local neighbours coordinates
    const int j_neighbours[3] = { max(x - mult, 0), x, min(x + mult, width - 1) };
    const int i_neighbours[3] = { max(y - mult, 0), y, min(y + mult, height - 1) };

    // fetch non-local pixels and store them locally and contiguously
    float4 neighbour_pixel_HF[9];
    neighbour_pixel_HF[3 * 0 + 0] = read_imagef(HF, samplerA, (int2)(j_neighbours[0], i_neighbours[0]));
    neighbour_pixel_HF[3 * 0 + 1] = read_imagef(HF, samplerA, (int2)(j_neighbours[1], i_neighbours[0]));
    neighbour_pixel_HF[3 * 0 + 2] = read_imagef(HF, samplerA, (int2)(j_neighbours[2], i_neighbours[0]));

    neighbour_pixel_HF[3 * 1 + 0] = read_imagef(HF, samplerA, (int2)(j_neighbours[0], i_neighbours[1]));
    neighbour_pixel_HF[3 * 1 + 1] = read_imagef(HF, samplerA, (int2)(j_neighbours[1], i_neighbours[1]));
    neighbour_pixel_HF[3 * 1 + 2] = read_imagef(HF, samplerA, (int2)(j_neighbours[2], i_neighbours[1]));

    neighbour_pixel_HF[3 * 2 + 0] = read_imagef(HF, samplerA, (int2)(j_neighbours[0], i_neighbours[2]));
    neighbour_pixel_HF[3 * 2 + 1] = read_imagef(HF, samplerA, (int2)(j_neighbours[1], i_neighbours[2]));
    neighbour_pixel_HF[3 * 2 + 2] = read_imagef(HF, samplerA, (int2)(j_neighbours[2], i_neighbours[2]));

    float4 update = 0.f;

    // Compute the laplacian in the direction parallel to the steepest gradient on the norm
    float anisotropic_kernel_isophote[9] = { 0.25f, 0.5f, 0.25f, 0.5f, -3.f, 0.5f, 0.25f, 0.5f, 0.25f };

    // Convolve the filter to get the laplacian
    float4 laplacian_HF = 0.f;
    for(int k = 0; k < 9; k++)
    {
      laplacian_HF += neighbour_pixel_HF[k] * anisotropic_kernel_isophote[k];
    }

    // Diffuse
    const float4 multipliers_HF = { 1.f / B_SPLINE_TO_LAPLACIAN, 1.f / B_SPLINE_TO_LAPLACIAN, 1.f / B_SPLINE_TO_LAPLACIAN, 0.f };
    high_frequency += alpha * multipliers_HF * (laplacian_HF - first_order_factor * high_frequency);

    high_frequency.w = norm_backup;
  }

  if(scale & FIRST_SCALE)
  {
    // out is not inited yet
    out = high_frequency;
  }
  else
  {
    // just accumulate HF
    out = readpixel(output_r, x, y) + high_frequency;
  }

  if(scale & LAST_SCALE)
  {
    // add the residual and clamp
    out = fmax(out + readpixel(LF, x, y), (float4)0.f);

    // renormalize ratios
    if(alpha.w > 0.f)
    {
      const float4 out_sq = sqf(out);
      const float norm = sqrt(out_sq.x + out_sq.y + out_sq.z);
      if(norm > 1e-4f) out.xyz /= norm;
    }

    // Last scale : reconstruct RGB from ratios and norm - norm stays in the 4th channel
    // we need it to evaluate the gradient
    out.xyz *= out.w;
  }

  write_imagef(output_w, (int2)(x, y), out);
}

float
lookup_unbounded_twosided(read_only image2d_t lut, const float x, constant float *a)
{
  // in case the tone curve is marked as linear, return the fast
  // path to linear unbounded (does not clip x at 1)
  if(a[0] >= 0.0f)
  {
    const float ar = 1.0f/a[0];
    const float al = 1.0f - 1.0f/a[3];
    if(x < ar && x >= al)
    {
      // lut lookup
      const int xi = clamp((int)(x * 0x10000ul), 0, 0xffff);
      const int2 p = (int2)((xi & 0xff), (xi >> 8));
      return read_imagef(lut, sampleri, p).x;
    }
    else
    {
      // two-sided extrapolation (with inverted x-axis for left side)
      const float xx = (x >= ar) ? x : 1.0f - x;
      constant float *aa = (x >= ar) ? a : a + 3;
      return aa[1] * dtcl_pow(xx*aa[0], aa[2]);
    }
  }
  else return x;
}

float
lerp_lookup_unbounded0(read_only image2d_t lut, const float x, global const float *a)
{
  // in case the tone curve is marked as linear, return the fast
  // path to linear unbounded (does not clip x at 1)
  if(a[0] >= 0.0f)
  {
    if(x < 1.0f)
    {
      const float ft = clamp(x * (float)0xffff, 0.0f, (float)0xffff);
      const int t = ft < 0xfffe ? ft : 0xfffe;
      const float f = ft - t;
      const int2 p1 = (int2)((t & 0xff), (t >> 8));
      const int2 p2 = (int2)(((t + 1) & 0xff), ((t + 1) >> 8));
      const float l1 = read_imagef(lut, sampleri, p1).x;
      const float l2 = read_imagef(lut, sampleri, p2).x;
      return l1 * (1.0f - f) + l2 * f;
    }
    else return a[1] * dtcl_pow(x*a[0], a[2]);
  }
  else return x;
}

static float4 _lookup_hsm(const float4 hsv,
                          const int hue_div,
                          const int sat_div,
                          const int val_div,
                          global const float *hsm)
{
  const float h = (hsv.x - floor(hsv.x)) * (float)hue_div;
  const float s = clipf(hsv.y) * (float)(sat_div-1);
  const float v = clipf(hsv.z) * (float)(val_div-1);
  const int hi[2] = { min((int)h, hue_div-1), (min((int)h,   hue_div-1)+1) % hue_div };
  const int si[2] = { min((int)s, sat_div-1),  min((int)s+1, sat_div-1) };
  const int vi[2] = { min((int)v, val_div-1),  min((int)v+1, val_div-1) };
  const float hf = clipf(h - (float)hi[0]);
  const float sf = clipf(s - (float)si[0]);
  const float vf = clipf(v - (float)vi[0]);

  float4 correction = (float4)0.0f;
  for(int z = 0; z < 2; z++)
    for(int y = 0; y < 2; y++)
      for(int x = 0; x < 2; x++)
      {
        // DNG stores saturation first, then hue, then value, with three floats per cell
        const int index = 3 * ((vi[z] * hue_div + hi[y]) * sat_div + si[x]);
        const float weight = (z ? vf : 1.0f - vf)
                           * (y ? hf : 1.0f - hf)
                           * (x ? sf : 1.0f - sf);
        correction.x += weight * hsm[index + 0];
        correction.y += weight * hsm[index + 1];
        correction.z += weight * hsm[index + 2];
      }
  return correction;
}

static float _srgb_to_linear(const float x)
{
  return x <= 0.04045f ? x / 12.92f : dtcl_pow((x + 0.055f) / 1.055f, 2.4f);
}
static float _linear_to_srgb(const float x)
{
  return x <= 0.0031308f ? 12.92f * x : 1.055f * dtcl_pow(x, 1.0f / 2.4f) - 0.055f;
}

/* kernel for the plugin colorin: unbound processing */
kernel void
colorin_unbound (read_only image2d_t in, write_only image2d_t out, const int width, const int height,
                 global float *cmat, global float *lmat,
                 read_only image2d_t lutr, read_only image2d_t lutg, read_only image2d_t lutb,
                 const int blue_mapping, global const float (*const a)[3], global const float *corr,
                 const float scale, const int hue_div, const int sat_div, const int val_div, global const float *dev_hsm, const int encoding, const int use_hsm, const int hdr_transfer)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  const float4 corval = (const float4)(corr[0], corr[1], corr[2], corr[3]);
  float4 pixel = corval * readpixel(in, x, y);
  const float alpha = pixel.w;

  float cam[3], XYZ[3];
  if(hdr_transfer)
  {
    cam[0] = pixel.x;
    cam[1] = pixel.y;
    cam[2] = pixel.z;
    dt_hdr_decode(cam, hdr_transfer);
  }
  else
  {
    cam[0] = lerp_lookup_unbounded0(lutr, pixel.x, a[0]);
    cam[1] = lerp_lookup_unbounded0(lutg, pixel.y, a[1]);
    cam[2] = lerp_lookup_unbounded0(lutb, pixel.z, a[2]);
  }

  if(blue_mapping)
  {
    const float YY = cam[0] + cam[1] + cam[2];
    if(YY > 0.0f)
    {
      // manual gamut mapping. these values cause trouble when converting back from Lab to sRGB:
      const float zz = cam[2] / YY;
      // lower amount and higher bound_z make the effect smaller.
      // the effect is weakened the darker input values are, saturating at bound_Y
      const float bound_z = 0.5f, bound_Y = 0.8f;
      const float amount = 0.11f;
      if (zz > bound_z)
      {
        const float t = (zz - bound_z) / (1.0f - bound_z) * fmin(1.0f, YY / bound_Y);
        cam[1] += t * amount;
        cam[2] -= t * amount;
      }
    }
  }

  // now convert camera to XYZ using the color matrix
  for(int j=0;j<3;j++)
  {
    XYZ[j] = 0.0f;
    for(int i=0;i<3;i++) XYZ[j] += cmat[3*j+i] * cam[i];
  }
  float4 xyz = (float4)(XYZ[0], XYZ[1], XYZ[2], 0.0f);
  pixel.xyz = XYZ_to_Lab(xyz).xyz;
  if(use_hsm)
  {
    float4 prophoto_rgb = Lab_to_prophotorgb(pixel) / (float4)scale;
    float4 hsv = RGB_2_HSV(prophoto_rgb);
    // Encode valueV for table lookup if needed
    float v_encoded = hsv.z;
    if(encoding == 1) // sRGB
      v_encoded = _linear_to_srgb(clipf(v_encoded));

    // Use encoded value for table lookup (need to temporarily set hsv[2])
    float4 hsv_for_lookup = hsv;
    hsv_for_lookup.z = v_encoded;
    float4 correction = _lookup_hsm(hsv_for_lookup, hue_div, sat_div, val_div, dev_hsm);

    hsv.x += correction.x / 360.0f;
    hsv.x -= floor(hsv.x); // stay within the 360° hue circle
    hsv.y = clipf(hsv.y * correction.y);

    float v_corrected = v_encoded * correction.z;
    if(encoding == 1)  // sRGB
      v_corrected = _srgb_to_linear(v_corrected);
    hsv.z = clipf(v_corrected);
    prophoto_rgb = HSV_2_RGB(hsv) * (float4)scale;
    pixel = prophotorgb_to_Lab(prophoto_rgb);
  }
  pixel.w = alpha;
  write_imagef(out, (int2)(x, y), pixel);
}

/* kernel for the plugin colorin: the LittleCMS transform baked into a 3D LUT
   over input values 0..1, see _bake_clut() */
kernel void
colorin_clut(read_only image2d_t in, write_only image2d_t out, const int width, const int height,
             global const float *clut, const int level, global const float *corr)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  const float4 corval = (const float4)(corr[0], corr[1], corr[2], corr[3]);
  const float4 raw = readpixel(in, x, y);
  float4 pixel = dt_lut3d_tetrahedral(clip4(corval * raw), clut, level);
  pixel.w = raw.w;
  write_imagef(out, (int2)(x, y), pixel);
}

/* kernel for the plugin colorin: with clipping */
kernel void
colorin_clipping (read_only image2d_t in, write_only image2d_t out, const int width, const int height,
                  global float *cmat, global float *lmat,
                  read_only image2d_t lutr, read_only image2d_t lutg, read_only image2d_t lutb,
                  const int blue_mapping, global const float (*const a)[3], global const float *corr,
                  const float scale, const int hue_div, const int sat_div, const int val_div, global const float *dev_hsm, const int encoding, const int use_hsm, const int hdr_transfer)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  const float4 corval = (const float4)(corr[0], corr[1], corr[2], corr[3]);
  float4 pixel = corval * readpixel(in, x, y);
  const float alpha = pixel.w;

  float cam[3], RGB[3], XYZ[3];
  if(hdr_transfer)
  {
    cam[0] = pixel.x;
    cam[1] = pixel.y;
    cam[2] = pixel.z;
    dt_hdr_decode(cam, hdr_transfer);
  }
  else
  {
    cam[0] = lerp_lookup_unbounded0(lutr, pixel.x, a[0]);
    cam[1] = lerp_lookup_unbounded0(lutg, pixel.y, a[1]);
    cam[2] = lerp_lookup_unbounded0(lutb, pixel.z, a[2]);
  }

  if(blue_mapping)
  {
    const float YY = cam[0] + cam[1] + cam[2];
    if(YY > 0.0f)
    {
      // manual gamut mapping. these values cause trouble when converting back from Lab to sRGB:
      const float zz = cam[2] / YY;
      // lower amount and higher bound_z make the effect smaller.
      // the effect is weakened the darker input values are, saturating at bound_Y
      const float bound_z = 0.5f, bound_Y = 0.8f;
      const float amount = 0.11f;
      if (zz > bound_z)
      {
        const float t = (zz - bound_z) / (1.0f - bound_z) * fmin(1.0f, YY / bound_Y);
        cam[1] += t * amount;
        cam[2] -= t * amount;
      }
    }
  }

  // convert camera to RGB using the first color matrix
  for(int j=0;j<3;j++)
  {
    RGB[j] = 0.0f;
    for(int i=0;i<3;i++) RGB[j] += cmat[3*j+i] * cam[i];
  }

  // clamp at this stage
  for(int i=0; i<3; i++) RGB[i] = clipf(RGB[i]);

  // convert clipped RGB to XYZ
  for(int j=0;j<3;j++)
  {
    XYZ[j] = 0.0f;
    for(int i=0;i<3;i++) XYZ[j] += lmat[3*j+i] * RGB[i];
  }

  float4 xyz = (float4)(XYZ[0], XYZ[1], XYZ[2], 0.0f);
  pixel.xyz = XYZ_to_Lab(xyz).xyz;
  if(use_hsm)
  {
    float4 prophoto_rgb = Lab_to_prophotorgb(pixel) / (float4)scale;
    float4 hsv = RGB_2_HSV(prophoto_rgb);
    // Encode valueV for table lookup if needed
    float v_encoded = hsv.z;
    if(encoding == 1)  // sRGB
      v_encoded = _linear_to_srgb(clipf(v_encoded));

    // Use encoded value for table lookup (need to temporarily set hsv[2])
    float4 hsv_for_lookup = hsv;
    hsv_for_lookup.z = v_encoded;
    float4 correction = _lookup_hsm(hsv_for_lookup, hue_div, sat_div, val_div, dev_hsm);

    hsv.x += correction.x / 360.0f;
    hsv.x -= floor(hsv.x); // stay within the 360° hue circle
    hsv.y = clipf(hsv.y * correction.y);

    float v_corrected = v_encoded * correction.z;
    if(encoding == 1)  // sRGB
      v_corrected = _srgb_to_linear(v_corrected);
    hsv.z = clipf(v_corrected);
    prophoto_rgb = HSV_2_RGB(hsv) * (float4)scale;

    pixel = prophotorgb_to_Lab(prophoto_rgb);
  }
  pixel.w = alpha;
  write_imagef(out, (int2)(x, y), pixel);
}

/* kernel for input profile conversion directly to working RGB */
kernel void
colorin_rgb (read_only image2d_t in, write_only image2d_t out, const int width, const int height,
                 global float *cmat, global float *lmat,
                 read_only image2d_t lutr, read_only image2d_t lutg, read_only image2d_t lutb,
                 const int blue_mapping, global const float (*const a)[3], global const float *corr,
                 const float scale, const int hue_div, const int sat_div, const int val_div, global const float *dev_hsm, const int encoding, const int use_hsm, const int hdr_transfer,
                 const int clipping, constant dt_colorspaces_iccprofile_info_cl_t *work,
                 read_only image2d_t work_lut, const int lab_input)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  const float4 corval = (const float4)(corr[0], corr[1], corr[2], corr[3]);
  const float4 input = readpixel(in, x, y);
  float4 pixel = corval * input;
  const float alpha = input.w;

  if(lab_input)
  {
    pixel = matrix_product(Lab_to_XYZ(input), work->matrix_out);
    if(work->nonlinearlut || work->hdr_transfer)
      pixel = apply_trc_out(pixel, work, work_lut);
    pixel.w = alpha;
    write_imagef(out, (int2)(x, y), pixel);
    return;
  }

  float cam[3];
  if(hdr_transfer)
  {
    cam[0] = pixel.x;
    cam[1] = pixel.y;
    cam[2] = pixel.z;
    dt_hdr_decode(cam, hdr_transfer);
  }
  else
  {
    cam[0] = lerp_lookup_unbounded0(lutr, pixel.x, a[0]);
    cam[1] = lerp_lookup_unbounded0(lutg, pixel.y, a[1]);
    cam[2] = lerp_lookup_unbounded0(lutb, pixel.z, a[2]);
  }

  if(blue_mapping)
  {
    const float YY = cam[0] + cam[1] + cam[2];
    if(YY > 0.0f)
    {
      // manual gamut mapping. these values cause trouble when converting back from Lab to sRGB:
      const float zz = cam[2] / YY;
      // lower amount and higher bound_z make the effect smaller.
      // the effect is weakened the darker input values are, saturating at bound_Y
      const float bound_z = 0.5f, bound_Y = 0.5f;
      const float amount = 0.11f;
      if(zz > bound_z)
      {
        const float t = (zz - bound_z) / (1.0f - bound_z) * fmin(1.0f, YY / bound_Y);
        cam[1] += t * amount;
        cam[2] -= t * amount;
      }
    }
  }

  // the first matrix maps directly to working RGB, or to the clipping gamut
  float4 linear = (float4)(0.0f);
  for(int j = 0; j < 3; j++)
    for(int i = 0; i < 3; i++) linear[j] += cmat[3*j+i] * cam[i];
  if(clipping)
  {
    float4 clipped = clamp(linear, 0.0f, 1.0f);
    linear = (float4)(0.0f);
    for(int j = 0; j < 3; j++)
      for(int i = 0; i < 3; i++) linear[j] += lmat[3*j+i] * clipped[i];
  }
  pixel = linear;
  if(use_hsm)
  {
    float4 prophoto_rgb = XYZ_to_prophotorgb(matrix_product(pixel, work->matrix_in)) / (float4)scale;
    float4 hsv = RGB_2_HSV(prophoto_rgb);
    // encode value for table lookup if needed
    float v_encoded = hsv.z;
    if(encoding == 1) // sRGB
      v_encoded = _linear_to_srgb(clipf(v_encoded));

    // use the encoded value for table lookup
    float4 hsv_for_lookup = hsv;
    hsv_for_lookup.z = v_encoded;
    float4 correction = _lookup_hsm(hsv_for_lookup, hue_div, sat_div, val_div, dev_hsm);

    hsv.x += correction.x / 360.0f;
    hsv.x -= floor(hsv.x); // stay within the 360° hue circle
    hsv.y = clipf(hsv.y * correction.y);

    float v_corrected = v_encoded * correction.z;
    if(encoding == 1)  // sRGB
      v_corrected = _srgb_to_linear(v_corrected);
    hsv.z = clipf(v_corrected);
    prophoto_rgb = HSV_2_RGB(hsv) * (float4)scale;
    pixel = matrix_product(prophotorgb_to_XYZ(prophoto_rgb), work->matrix_out);
  }
  if(work->nonlinearlut || work->hdr_transfer)
    pixel = apply_trc_out(pixel, work, work_lut);
  pixel.w = alpha;
  write_imagef(out, (int2)(x, y), pixel);
}

/* kernel for the tonecurve plugin. */
kernel void
tonecurve (read_only image2d_t in, write_only image2d_t out, const int width, const int height,
           read_only image2d_t table_L, read_only image2d_t table_a, read_only image2d_t table_b,
           const int autoscale_ab, const int unbound_ab, constant float *coeffs_L, constant float *coeffs_ab,
           const float low_approximation, const int preserve_colors,
           constant dt_colorspaces_iccprofile_info_cl_t *profile_info, read_only image2d_t lut)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  float4 pixel = readpixel(in, x, y);
  const float L_in = pixel.x/100.0f;
  // use lut or extrapolation:
  const float L = lookup_unbounded_last_node(table_L, L_in, coeffs_L);

  if (autoscale_ab == 0)
  {
    const float a_in = (pixel.y + 128.0f) / 256.0f;
    const float b_in = (pixel.z + 128.0f) / 256.0f;

    if (unbound_ab == 0)
    {
      pixel.y = lookup(table_a, a_in);
      pixel.z = lookup(table_b, b_in);
    }
    else
    {
      // use lut or two-sided extrapolation
      pixel.y = lookup_unbounded_twosided(table_a, a_in, coeffs_ab);
      pixel.z = lookup_unbounded_twosided(table_b, b_in, coeffs_ab + 6);
    }
    pixel.x = L;
  }
  else if(autoscale_ab == 1)
  {
    if(L_in > 0.01f)
    {
      pixel.y *= L/pixel.x;
      pixel.z *= L/pixel.x;
    }
    else
    {
      pixel.y *= low_approximation;
      pixel.z *= low_approximation;
    }
    pixel.x = L;
  }
  else if(autoscale_ab == 2)
  {
    float4 xyz = Lab_to_XYZ(pixel);
    xyz.x = lookup_unbounded_last_node(table_L, xyz.x, coeffs_L);
    xyz.y = lookup_unbounded_last_node(table_L, xyz.y, coeffs_L);
    xyz.z = lookup_unbounded_last_node(table_L, xyz.z, coeffs_L);
    pixel.xyz = XYZ_to_Lab(xyz).xyz;
  }
  else if(autoscale_ab == 3)
  {
    float4 rgb = Lab_to_prophotorgb(pixel);

    if (preserve_colors == DT_RGB_NORM_NONE)
    {
      rgb.x = lookup_unbounded_last_node(table_L, rgb.x, coeffs_L);
      rgb.y = lookup_unbounded_last_node(table_L, rgb.y, coeffs_L);
      rgb.z = lookup_unbounded_last_node(table_L, rgb.z, coeffs_L);
    }
    else
    {
      float ratio = 1.f;
      float lum = dt_rgb_norm(rgb, preserve_colors, 1, profile_info, lut);
      if(lum > 0.f)
      {
        const float curve_lum = lookup_unbounded_last_node(table_L, lum, coeffs_L);
        ratio = curve_lum / lum;
      }
      rgb.xyz *= ratio;
    }
    pixel.xyz = prophotorgb_to_Lab(rgb).xyz;
  }

  write_imagef (out, (int2)(x, y), pixel);
}


/* kernel for the colorcorrection plugin. */
__kernel void
colorcorrection (read_only image2d_t in, write_only image2d_t out, const int width, const int height,
                 const float saturation, const float a_scale, const float a_base,
                 const float b_scale, const float b_base)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  float4 pixel = readpixel(in, x, y);
  pixel.y = saturation*(pixel.y + pixel.x * a_scale + a_base);
  pixel.z = saturation*(pixel.z + pixel.x * b_scale + b_base);
  write_imagef (out, (int2)(x, y), pixel);
}


void
mul_mat_vec_2(const float4 m, const float2 *p, float2 *o)
{
  (*o).x = (*p).x*m.x + (*p).y*m.y;
  (*o).y = (*p).x*m.z + (*p).y*m.w;
}

void
backtransform(float2 *p, float2 *o, const float4 m, const float2 t)
{
  (*p).y /= (1.0f + (*p).x*t.x);
  (*p).x /= (1.0f + (*p).y*t.y);
  mul_mat_vec_2(m, p, o);
}

void
keystone_backtransform(float2 *i, const float4 k_space, const float2 ka, const float4 ma, const float2 mb)
{
  float xx = (*i).x - k_space.x;
  float yy = (*i).y - k_space.y;

  /*float u = ka.x-kb.x+kc.x-kd.x;
  float v = ka.x-kb.x;
  float w = ka.x-kd.x;
  float z = ka.x;
  //(*i).x = (xx/k_space.z)*(yy/k_space.w)*(ka.x-kb.x+kc.x-kd.x) - (xx/k_space.z)*(ka.x-kb.x) - (yy/k_space.w)*(ka.x-kd.x) + ka.x + k_space.x;
  (*i).x = (xx/k_space.z)*(yy/k_space.w)*u - (xx/k_space.z)*v - (yy/k_space.w)*w + z + k_space.x;
  u = ka.y-kb.y+kc.y-kd.y;
  v = ka.y-kb.y;
  w = ka.y-kd.y;
  z = ka.y;
  //(*i).y = (xx/k_space.z)*(yy/k_space.w)*(ka.y-kb.y+kc.y-kd.y) - (xx/k_space.z)*(ka.y-kb.y) - (yy/k_space.w)*(ka.y-kd.y) + ka.y + k_space.y;
  (*i).y = (xx/k_space.z)*(yy/k_space.w)*u - (xx/k_space.z)*v - (yy/k_space.w)*w + z + k_space.y;*/
  float div = ((ma.z*xx-ma.x*yy)*mb.y+(ma.y*yy-ma.w*xx)*mb.x+ma.x*ma.w-ma.y*ma.z);

  (*i).x = (ma.w*xx-ma.y*yy)/div + ka.x;
  (*i).y =-(ma.z*xx-ma.x*yy)/div + ka.y;
}


float
interpolation_func_bicubic(float t)
{
  float r;
  t = fabs(t);

  r = (t >= 2.0f) ? 0.0f : ((t > 1.0f) ? (0.5f*(t*(-t*t + 5.0f*t - 8.0f) + 4.0f)) : (0.5f*(t*(3.0f*t*t - 5.0f*t) + 2.0f)));

  return r;
}

#define DT_LANCZOS_EPSILON (1e-9f)

#if 0
float
interpolation_func_lanczos(float width, float t)
{
  float ta = fabs(t);

  float r = (ta > width) ? 0.0f : ((ta < DT_LANCZOS_EPSILON) ? 1.0f : width*native_sin(M_PI_F*t)*native_sin(M_PI_F*t/width)/(M_PI_F*M_PI_F*t*t));

  return r;
}
#else
float
sinf_fast(float t)
{
  /***** if you change this function, you must also change the copy in src/common/math.h *****/
  const float a = 4.0f/(M_PI_F*M_PI_F);
  const float p = 0.225f;

  t = a*t*(M_PI_F - fabs(t));

  return p*(t*fabs(t) - t) + t;
}

float
interpolation_func_lanczos(float width, float t)
{
  /* Compute a value for sinf(pi.t) in [-pi pi] for which the value will be
   * correct */
  int a = (int)t;
  float r = t - (float)a;

  // Compute the correct sign for sinf(pi.r)
  union { float f; unsigned int i; } sign;
  sign.i = ((a&1)<<31) | 0x3f800000;

  return (DT_LANCZOS_EPSILON + width*sign.f*sinf_fast(M_PI_F*r)*sinf_fast(M_PI_F*t/width))/(DT_LANCZOS_EPSILON + M_PI_F*M_PI_F*t*t);
}
#endif


// the taps of _maketaps_bilinear(), _maketaps_bicubic() and _maketaps_lanczos()
// in src/common/interpolation.c, by the dt_interpolation_type of the interpolator
static inline float _interpolation_tap(const int type, const float width, const float t)
{
  if(type == 0) return 1.0f - fabs(t);
  if(type == 1)
  {
    const float ta = fabs(t);
    return ta <= 1.0f ? 0.5f * (ta * (3.0f * ta * ta - 5.0f * ta) + 2.0f)
                      : 0.5f * (ta * (-ta * ta + 5.0f * ta - 8.0f) + 4.0f);
  }
  return interpolation_func_lanczos(width, t);
}

// one pixel by dt_interpolation_compute_pixel4c() at (px, py): the interpolator's
// taps and norm, mirrored borders and zero outside the image
static inline float4 _interpolation_pixel4(read_only image2d_t in,
                                           const float px,
                                           const float py,
                                           const int in_width,
                                           const int in_height,
                                           const int itype,
                                           const int iwidth)
{
  const int ix = (int)floor(px);
  const int iy = (int)floor(py);

  if(ix < -1 || iy < -1 || ix >= in_width || iy >= in_height)
    return (float4)0.0f;

  const int fx = ix - iwidth + 1;
  const int fy = iy - iwidth + 1;
  const float tx = px - (float)fx;
  const float ty = py - (float)fy;

  float kh[6], kv[6];
  float normh = 0.0f, normv = 0.0f;
  for(int k = 0; k < 2 * iwidth; k++)
  {
    kh[k] = _interpolation_tap(itype, iwidth, tx - k);
    kv[k] = _interpolation_tap(itype, iwidth, ty - k);
    normh += kh[k];
    normv += kv[k];
  }
  // the CPU kernels report a norm of 1 for bilinear and bicubic
  if(itype < 2) normh = normv = 1.0f;

  float4 pixel = (float4)0.0f;
  for(int i = 0; i < 2 * iwidth; i++)
  {
    const int yy = clip_mirror(fy + i, in_height - 1);
    float4 h = (float4)0.0f;
    for(int j = 0; j < 2 * iwidth; j++)
      h += kh[j] * readpixel(in, clip_mirror(fx + j, in_width - 1), yy);
    pixel += kv[i] * h;
  }

  return pixel * (1.0f / (normh * normv));
}

// rotate pixels: one output pixel at the back transformed position
kernel void
rotatepixels(read_only image2d_t in,
             write_only image2d_t out,
             const int width,
             const int height,
             const int in_width,
             const int in_height,
             const int2 roi_out,
             const int2 roi_in,
             const float4 m,
             const float2 center,
             const int itype,
             const int iwidth)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);
  if(x >= width || y >= height) return;

  const float pix = roi_out.x + x;
  const float piy = roi_out.y + y;
  const float px = pix * m.x + piy * m.y + center.x - roi_in.x;
  const float py = pix * m.z + piy * m.w + center.y - roi_in.y;

  write_imagef(out, (int2)(x, y), _interpolation_pixel4(in, px, py, in_width, in_height, itype, iwidth));
}

// scale pixels: one output pixel at ((x, y) + shift) * scale - offset of the
// input, see _sampling_grid() in src/iop/scalepixels.c
kernel void
scalepixels(read_only image2d_t in,
            write_only image2d_t out,
            const int width,
            const int height,
            const int in_width,
            const int in_height,
            const float2 scale,
            const int2 shift,
            const float2 offset,
            const int itype,
            const int iwidth)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);
  if(x >= width || y >= height) return;

  const float px = (x + shift.x) * scale.x - offset.x;
  const float py = (y + shift.y) * scale.y - offset.y;

  write_imagef(out, (int2)(x, y), _interpolation_pixel4(in, px, py, in_width, in_height, itype, iwidth));
}

/* kernel for clip&rotate: bilinear interpolation */
__kernel void
clip_rotate_bilinear(read_only image2d_t in, write_only image2d_t out, const int width, const int height,
            const int in_width, const int in_height,
            const int2 roi_in, const float2 roi_out, const float scale_in, const float scale_out,
            const int flip, const float2 t, const float2 k, const float4 mat,
            const float4 k_space, const float2 ka, const float4 ma, const float2 mb)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  float2 pi, po;

  pi.x = roi_out.x + x + 0.5f;
  pi.y = roi_out.y + y + 0.5f;

  pi.x -= flip ? t.y * scale_out : t.x * scale_out;
  pi.y -= flip ? t.x * scale_out : t.y * scale_out;

  pi /= scale_out;
  backtransform(&pi, &po, mat, k);
  po *= scale_in;

  po.x += t.x * scale_in;
  po.y += t.y * scale_in;

  if (k_space.z > 0.0f) keystone_backtransform(&po,k_space,ka,ma,mb);

  po.x -= roi_in.x + 0.5f;
  po.y -= roi_in.y + 0.5f;

  const int ii = (int)po.x;
  const int jj = (int)po.y;

  float4 o = (ii >=0 && jj >= 0 && ii < in_width && jj < in_height) ? read_imagef(in, samplerf, po) : (float4)0.0f;

  write_ipixel(out, (int2)(x, y), o);
}



/* kernel for clip&rotate: bicubic interpolation */
__kernel void
clip_rotate_bicubic(read_only image2d_t in,
                    write_only image2d_t out,
                    const int width,
                    const int height,
                    const int in_width,
                    const int in_height,
                    const int2 roi_in,
                    const float2 roi_out,
                    const float scale_in,
                    const float scale_out,
                    const int flip,
                    const float2 t,
                    const float2 k,
                    const float4 mat,
                    const float4 k_space,
                    const float2 ka,
                    const float4 ma,
                    const float2 mb)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  const int kwidth = 2;

  if(x >= width || y >= height) return;

  float2 pi, po;

  pi.x = roi_out.x + x + 0.5f;
  pi.y = roi_out.y + y + 0.5f;

  pi.x -= flip ? t.y * scale_out : t.x * scale_out;
  pi.y -= flip ? t.x * scale_out : t.y * scale_out;

  pi /= scale_out;
  backtransform(&pi, &po, mat, k);
  po *= scale_in;

  po.x += t.x * scale_in;
  po.y += t.y * scale_in;

  if (k_space.z > 0.0f) keystone_backtransform(&po,k_space,ka,ma,mb);

  po.x -= roi_in.x + 0.5f;
  po.y -= roi_in.y + 0.5f;

  int tx = po.x;
  int ty = po.y;

  float4 pixel = (float4)0.0f;
  float weight = 0.0f;

  for(int jj = 1 - kwidth; jj <= kwidth; jj++)
    for(int ii= 1 - kwidth; ii <= kwidth; ii++)
  {
    const int i = tx + ii;
    const int j = ty + jj;

    float wx = interpolation_func_bicubic((float)i - po.x);
    float wy = interpolation_func_bicubic((float)j - po.y);
    float w = wx * wy;

    pixel += Areadpixel(in, clip_mirror(i, in_width-1), clip_mirror(j, in_height-1)) * (float4)w;
    weight += w;
  }

  pixel = (tx >= 0 && ty >= 0 && tx < in_width && ty < in_height) ? pixel / weight : (float4)0.0f;

  write_ipixel(out, (int2)(x, y), pixel);
}


/* kernel for clip&rotate: lanczos2 interpolation */
__kernel void
clip_rotate_lanczos2(read_only image2d_t in, write_only image2d_t out, const int width, const int height,
            const int in_width, const int in_height,
            const int2 roi_in, const float2 roi_out, const float scale_in, const float scale_out,
            const int flip, const float2 t, const float2 k, const float4 mat,
            const float4 k_space, const float2 ka, const float4 ma, const float2 mb)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  const int kwidth = 2;

  if(x >= width || y >= height) return;

  float2 pi, po;

  pi.x = roi_out.x + x + 0.5f;
  pi.y = roi_out.y + y + 0.5f;

  pi.x -= flip ? t.y * scale_out : t.x * scale_out;
  pi.y -= flip ? t.x * scale_out : t.y * scale_out;

  pi /= scale_out;
  backtransform(&pi, &po, mat, k);
  po *= scale_in;

  po.x += t.x * scale_in;
  po.y += t.y * scale_in;

  if (k_space.z > 0.0f) keystone_backtransform(&po,k_space,ka,ma,mb);

  po.x -= roi_in.x + 0.5f;
  po.y -= roi_in.y + 0.5f;

  int tx = po.x;
  int ty = po.y;

  float4 pixel = (float4)0.0f;
  float weight = 0.0f;

  for(int jj = 1 - kwidth; jj <= kwidth; jj++)
    for(int ii= 1 - kwidth; ii <= kwidth; ii++)
  {
    const int i = tx + ii;
    const int j = ty + jj;

    float wx = interpolation_func_lanczos(2, (float)i - po.x);
    float wy = interpolation_func_lanczos(2, (float)j - po.y);
    float w = wx * wy;

    pixel += Areadpixel(in, clip_mirror(i, in_width-1), clip_mirror(j, in_height-1)) * (float4)w;
    weight += w;
  }

  pixel = (tx >= 0 && ty >= 0 && tx < in_width && ty < in_height) ? pixel / weight : (float4)0.0f;

  write_ipixel(out, (int2)(x, y), pixel);
}



/* kernel for clip&rotate: lanczos3 interpolation */
__kernel void
clip_rotate_lanczos3(read_only image2d_t in, write_only image2d_t out, const int width, const int height,
            const int in_width, const int in_height,
            const int2 roi_in, const float2 roi_out, const float scale_in, const float scale_out,
            const int flip, const float2 t, const float2 k, const float4 mat,
            const float4 k_space, const float2 ka, const float4 ma, const float2 mb)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  const int kwidth = 3;

  if(x >= width || y >= height) return;

  float2 pi, po;

  pi.x = roi_out.x + x + 0.5f;
  pi.y = roi_out.y + y + 0.5f;

  pi.x -= flip ? t.y * scale_out : t.x * scale_out;
  pi.y -= flip ? t.x * scale_out : t.y * scale_out;

  pi /= scale_out;
  backtransform(&pi, &po, mat, k);
  po *= scale_in;

  po.x += t.x * scale_in;
  po.y += t.y * scale_in;

  if (k_space.z > 0.0f) keystone_backtransform(&po,k_space,ka,ma,mb);

  po.x -= roi_in.x + 0.5f;
  po.y -= roi_in.y + 0.5f;

  int tx = (int)po.x;
  int ty = (int)po.y;

  float4 pixel = (float4)0.0f;
  float weight = 0.0f;

  for(int jj = 1 - kwidth; jj <= kwidth; jj++)
    for(int ii= 1 - kwidth; ii <= kwidth; ii++)
  {
    const int i = tx + ii;
    const int j = ty + jj;

    float wx = interpolation_func_lanczos(3, (float)i - po.x);
    float wy = interpolation_func_lanczos(3, (float)j - po.y);
    float w = wx * wy;

    pixel += Areadpixel(in, clip_mirror(i, in_width-1), clip_mirror(j, in_height-1)) * (float4)w;
    weight += w;
  }

  pixel = (tx >= 0 && ty >= 0 && tx < in_width && ty < in_height) ? pixel / weight : (float4)0.0f;

  write_ipixel(out, (int2)(x, y), pixel);
}

float maketaps_bilinear(float *taps,
                        const int num_taps,
                        const float width,
                        const float first_tap,
                        const float interval)
{
  float iter[4];
  float vt[4];
  for(int c = 0; c < 4; c++)
    iter[c] = 4.0f * interval;
  for(int c = 0; c < 4; c++)
    vt[c] = first_tap + (float)c * interval;

  const int runs = (num_taps + 3) / 4;

  for(int i = 0; i < runs; i++)
  {
    for(int c = 0; c < 4; c++)
      taps[4*i + c] = 1.0f - (vt[c] < 0.0f ? -vt[c] : vt[c]);
    // prepare next iteration
    for(int c = 0; c < 4; c++)
      vt[c] += iter[c];
  }
  return 1.0f; //kernel norm is 1.0f by construction
}

float maketaps_bicubic(float *taps,
                       const int num_taps,
                       const float width,
                       const float first_tap,
                       const float interval)
{
  float iter[4];
  float vt[4];
  for(int c = 0; c < 4; c++)
    iter[c] = 4.0f * interval;
  for(int c = 0; c < 4; c++)
    vt[c] = first_tap + (float)c * interval;

  const int runs = (num_taps + 3) / 4;

  for(int i = 0; i < runs; i++)
  {
    // compute and store the values for the current four taps
    float vt_abs[4];
    float t2[4];   // tap-squared
    for(int c = 0; c < 4; c++)
    {
      vt_abs[c] = vt[c] < 0.0f ? -vt[c] : vt[c];
      t2[c] = vt[c] * vt[c];
    }
    float t5[4];
    float mt2_add_t5_sub_8[4];
    for(int c = 0; c < 4; c++)
    {
      t5[c] = 5.0f * vt_abs[c];
      mt2_add_t5_sub_8[c] = t5[c] - 8.0f - t2[c];
    }
    float b[4];
    float r12[4];
    for(int c = 0; c < 4; c++)
    {
      b[c] = vt_abs[c] * mt2_add_t5_sub_8[c] + 4.0f;
      r12[c] = b[c] * 0.5f; // the value for 1 < t < 2
    }
    float t23[4];
    float e[4];
    float r01[4];
    for(int c = 0; c < 4; c++)
    {
      t23[c] = 3.0f * t2[c] - t5[c];
      e[c] = t23[c] * vt_abs[c] + 2.0f;
      r01[c] = e[c] * 0.5f;
    }
    // combine the values depending on whether abs(tap) is less than one or not
    for(int c = 0; c < 4; c++)
    {
      taps[4*i + c] = vt_abs[c] <= 1.0f ? r01[c] : r12[c];
    }
    // prepare next iteration
    for(int c = 0; c < 4; c++)
      vt[c] += iter[c];
  }
  return 1.0f; //kernel norm is 1.0f by construction
}

void vector_sin(const float *arg, float *sine)
{
  const float a = 4.0f / (M_PI_F * M_PI_F);
  float abs_arg[4];
  for(int c = 0; c < 4; c++)
    abs_arg[c] = (arg[c] < 0.0f) ? -arg[c] : arg[c];

  float scaled[4];
  for(int c = 0; c < 4; c++)
    scaled[c] = a * arg[c] * (M_PI_F - abs_arg[c]);

  float abs_scaled[4];
  for(int c = 0; c < 4; c++)
    abs_scaled[c] = (scaled[c] < 0.0f) ? -scaled[c] : scaled[c];
  for(int c = 0; c < 4; c++)
    sine[c] = scaled[c] * (0.225f * (abs_scaled[c] - 1.0f) + 1.0f);
}

float maketaps_lanczos(float *taps,
                       const int num_taps,
                       const float width,
                       const float first_tap,
                       const float interval)
{
  float iter[4];
  float vt[4];
  for(int c = 0; c < 4; c++)
    iter[c] = 4.0f * interval;
  for(int c = 0; c < 4; c++)
    vt[c] = first_tap + (float)c * interval;
  float vw[4];
  for(int c = 0; c < 4; c++)
    vw[c] = width;

  const int runs = (num_taps + 3) / 4;

  for(int i = 0; i < runs; i++)
  {
    float r[4];
    float sign[4];
    for(int c = 0; c < 4; c++)
    {
      const int a = (int)vt[c];
      r[c] = vt[c] - (float)a;
      sign[c] = (a & 1) ? -1.0f : 1.0f;
    }
    float sine_arg1[4];
    float sine_arg2[4];
    for(int c = 0; c < 4; c++)
    {
      sine_arg1[c] = M_PI_F * r[c];
      sine_arg2[c] = M_PI_F * vt[c] / vw[c];
    }
    float sine1[4];
    float sine2[4];
    vector_sin(sine_arg1, sine1);
    vector_sin(sine_arg2, sine2);
    float num[4];
    float denom[4];
    for(int c = 0; c < 4; c++)
    {
      num[c] = (vw[c] * sign[c] * sine1[c] * sine2[c]) + 1e-9f;
      denom[c] = (M_PI_F*M_PI_F * vt[c] * vt[c]) + 1e-9f;
    }
    for(int c = 0; c < 4; c++)
    {
      taps[4*i + c] = num[c] / denom[c];
    }
    // prepare next iteration
    for(int c = 0; c < 4; c++)
      vt[c] += iter[c];
  }
  float norm = 0.0f;
  for(int i = 0; i < num_taps; i++)
    norm += taps[i];
  return norm;
}

float compute_upsampling_taps(const int itor_mode,
                              const int itor_width,
                              float *taps,
                              float tt)
{
  const int f = (int)floor(tt) - itor_width + 1;
  const float t = tt - (float)f;

  // itor_mode is the dt_interpolation_type: lanczos2 and lanczos3 share the
  // lanczos taps and differ only in width, as in src/common/interpolation.c
  if(itor_mode == 1)
    return maketaps_bicubic(taps, 2*itor_width, (float)itor_width, t, -1.0f);
  else if(itor_mode == 2 || itor_mode == 3)
    return maketaps_lanczos(taps, 2*itor_width, (float)itor_width, t, -1.0f);
  else
    return maketaps_bilinear(taps, 2*itor_width, (float)itor_width, t, -1.0f);
}

static inline float _get_image_channel(read_only image2d_t in,
                                      const int x,
                                      const int y,
                                      const int c)
{
  float4 pixel = Areadpixel(in, x, y);  // we ensure secure locations in callers
  if(c == 0)
    return pixel.x;
  else if(c == 1)
    return pixel.y;
  else if(c == 2)
    return pixel.z;

  return pixel.w;
}

// Keep in sync with defines in interpolation
#define MAX_HALF_FILTER_WIDTH 3
#define MAX_KERNEL_REQ ((2 * (MAX_HALF_FILTER_WIDTH) + 3) & (~3))

float interpolation_compute_single(read_only image2d_t in,
                                   const int itor_mode,
                                   const int itor_width,
                                   const float x,
                                   const float y,
                                   const int width,
                                   const int height,
                                   const int channel)
{
  float kernelh[MAX_KERNEL_REQ];
  float kernelv[MAX_KERNEL_REQ];

  int ix = (int)x;
  int iy = (int)y;

  float s = 0.0f;
  if(ix >= 0 && iy >= 0 && ix < width && iy < height)
  {
    // Compute both horizontal and vertical kernels
    float normh = compute_upsampling_taps(itor_mode, itor_width, kernelh, x);
    float normv = compute_upsampling_taps(itor_mode, itor_width, kernelv, y);

    iy -= itor_width - 1;
    ix -= itor_width - 1;

    // Apply the kernel
    for(int i = 0; i < 2 * itor_width; i++)
    {
      float h = 0.0f;
      for(int j = 0; j < 2 * itor_width; j++)
      {
        h += kernelh[j] * _get_image_channel(in, clip_mirror(ix+j, width-1), clip_mirror(iy+i, height-1), channel);
      }
      s += kernelv[i] * h;
    }
    s = s / (normh * normv);
  }
  return s;
}

#undef MAX_KERNEL_REQ
#undef MAX_HALF_FILTER_WIDTH

/* lens plugin: the lensfun values of pixel (x, y), interpolated bilinearly
   from a grid of every step-th pixel, gwidth points to a row and stride floats
   to a point. This is the OpenCL version of _lens_grid_row() in
   src/iop/lens.cc */
static inline float _lens_grid_lerp(const float a, const float b, const float t)
{
  // mixf() of the CPU path, except that a weight of zero skips b: a grid of
  // step 1 can hold non-finite values, which must not leak into the
  // neighboring pixels
  return t > 0.0f ? fma(t, b - a, a) : a;
}

static inline void _lens_grid_sample(global const float *grid, const int x, const int y,
                                     const int step, const int gwidth, const int stride,
                                     const int count, float *out)
{
  const int i = x / step;
  const int j = y / step;
  const float tx = (float)(x - i * step) / (float)step;
  const float ty = (float)(y - j * step) / (float)step;
  global const float *p00 = grid + (j * gwidth + i) * stride;
  global const float *p01 = p00 + gwidth * stride;
  for(int n = 0; n < count; n++)
  {
    const float left = _lens_grid_lerp(p00[n], p01[n], ty);
    const float right = _lens_grid_lerp(p00[stride + n], p01[stride + n], ty);
    out[n] = _lens_grid_lerp(left, right, tx);
  }
}

/* lens plugin: a lensfun coordinate relative to the input, held inside it as
   in the CPU path of src/iop/lens.cc. Coordinates past the edge sample the
   edge pixel there; mirroring the taps alone would reflect them back into the
   image. fmin() and fmax() also map NaN to an edge as the CPU path does. */
static inline float _lens_coord(const float v, const int offset, const int size)
{
  return fmax(fmin(v - (float)offset, (float)size - 1.0f), 0.0f);
}

/* kernels for the lens plugin: bilinear interpolation */
kernel void
lens_distort_bilinear (read_only image2d_t in, write_only image2d_t out, const int width, const int height,
               const int iwidth, const int iheight, const int roi_in_x, const int roi_in_y, global float *pi,
               const int do_nan_checks, const int step, const int gwidth)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  float4 pixel;

  float rx, ry;
  float ppi[6];
  _lens_grid_sample(pi, x, y, step, gwidth, 6, 6, ppi);

  if(do_nan_checks)
  {
    bool valid = true;

    for(int i = 0; i < 6; i++) valid = valid && isfinite(ppi[i]);

    if(!valid)
    {
      pixel = (float4)0.0f;
      write_imagef (out, (int2)(x, y), pixel);
      return;
    }
  }

  rx = _lens_coord(ppi[0], roi_in_x, iwidth);
  ry = _lens_coord(ppi[1], roi_in_y, iheight);
  {
    const float fx = floor(rx);
    const float fy = floor(ry);
    const float dx = rx - fx;
    const float dy = ry - fy;

    const int x0 = clip_mirror((int)fx,     iwidth - 1);
    const int x1 = clip_mirror((int)fx + 1, iwidth - 1);
    const int y0 = clip_mirror((int)fy,     iheight - 1);
    const int y1 = clip_mirror((int)fy + 1, iheight - 1);

    const float4 px = (1.0f - dy) * ((1.0f - dx) * Areadpixel(in, x0, y0) + dx * Areadpixel(in, x1, y0))
                    + dy  * ((1.0f - dx) * Areadpixel(in, x0, y1) + dx * Areadpixel(in, x1, y1));
    pixel.x = px.x;
  }

  rx = _lens_coord(ppi[2], roi_in_x, iwidth);
  ry = _lens_coord(ppi[3], roi_in_y, iheight);
  {
    const float fx = floor(rx);
    const float fy = floor(ry);
    const float dx = rx - fx;
    const float dy = ry - fy;

    const int x0 = clip_mirror((int)fx,     iwidth - 1);
    const int x1 = clip_mirror((int)fx + 1, iwidth - 1);
    const int y0 = clip_mirror((int)fy,     iheight - 1);
    const int y1 = clip_mirror((int)fy + 1, iheight - 1);

    const float4 px = (1.0f - dy) * ((1.0f - dx) * Areadpixel(in, x0, y0) + dx * Areadpixel(in, x1, y0))
                    + dy  * ((1.0f - dx) * Areadpixel(in, x0, y1) + dx * Areadpixel(in, x1, y1));
    pixel.yw = px.yw;
  }

  rx = _lens_coord(ppi[4], roi_in_x, iwidth);
  ry = _lens_coord(ppi[5], roi_in_y, iheight);
  {
    const float fx = floor(rx);
    const float fy = floor(ry);
    const float dx = rx - fx;
    const float dy = ry - fy;

    const int x0 = clip_mirror((int)fx,     iwidth - 1);
    const int x1 = clip_mirror((int)fx + 1, iwidth - 1);
    const int y0 = clip_mirror((int)fy,     iheight - 1);
    const int y1 = clip_mirror((int)fy + 1, iheight - 1);

    const float4 px = (1.0f - dy) * ((1.0f - dx) * Areadpixel(in, x0, y0) + dx * Areadpixel(in, x1, y0))
                    + dy  * ((1.0f - dx) * Areadpixel(in, x0, y1) + dx * Areadpixel(in, x1, y1));
    pixel.z = px.z;
  }

  pixel = all(isfinite(pixel.xyz)) ? pixel : (float4)0.0f;

  write_ipixel(out, (int2)(x, y), pixel);
}

/* kernels for the lens plugin: bicubic interpolation */
kernel void
lens_distort_bicubic (read_only image2d_t in,
                      write_only image2d_t out,
                      const int width,
                      const int height,
                      const int iwidth,
                      const int iheight,
                      const int roi_in_x,
                      const int roi_in_y,
                      global float *pi,
                      const int do_nan_checks,
                      const int step,
                      const int gwidth)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  const int kwidth = 2;

  if(x >= width || y >= height) return;

  float4 pixel = (float4)0.0f;

  float rx, ry;
  int tx, ty;
  float sum, weight;
  float2 sum2;
  float ppi[6];
  _lens_grid_sample(pi, x, y, step, gwidth, 6, 6, ppi);

  if(do_nan_checks)
  {
    bool valid = true;

    for(int i = 0; i < 6; i++) valid = valid && isfinite(ppi[i]);

    if(!valid)
    {
      pixel = (float4)0.0f;
      write_imagef (out, (int2)(x, y), pixel);
      return;
    }
  }

  rx = _lens_coord(ppi[0], roi_in_x, iwidth);
  ry = _lens_coord(ppi[1], roi_in_y, iheight);
  tx = floor(rx);
  ty = floor(ry);

  sum = 0.0f;
  weight = 0.0f;
  for(int jj = 1 - kwidth; jj <= kwidth; jj++)
    for(int ii= 1 - kwidth; ii <= kwidth; ii++)
  {
    int i = tx + ii;
    int j = ty + jj;

    float wx = interpolation_func_bicubic((float)i - rx);
    float wy = interpolation_func_bicubic((float)j - ry);
    float w = wx * wy;

    sum += Areadpixel(in, clip_mirror(i, iwidth-1), clip_mirror(j, iheight-1)).x * w;
    weight += w;
  }
  pixel.x = sum/weight;

  rx = _lens_coord(ppi[2], roi_in_x, iwidth);
  ry = _lens_coord(ppi[3], roi_in_y, iheight);
  tx = floor(rx);
  ty = floor(ry);

  sum2 = (float2)0.0f;
  weight = 0.0f;
  for(int jj = 1 - kwidth; jj <= kwidth; jj++)
    for(int ii= 1 - kwidth; ii <= kwidth; ii++)
  {
    int i = tx + ii;
    int j = ty + jj;

    float wx = interpolation_func_bicubic((float)i - rx);
    float wy = interpolation_func_bicubic((float)j - ry);
    float w = wx * wy;

    sum2 += Areadpixel(in, clip_mirror(i, iwidth-1), clip_mirror(j, iheight-1)).yw * w;
    weight += w;
  }
  pixel.yw = sum2/weight;

  rx = _lens_coord(ppi[4], roi_in_x, iwidth);
  ry = _lens_coord(ppi[5], roi_in_y, iheight);
  tx = floor(rx);
  ty = floor(ry);

  sum = 0.0f;
  weight = 0.0f;
  for(int jj = 1 - kwidth; jj <= kwidth; jj++)
    for(int ii= 1 - kwidth; ii <= kwidth; ii++)
  {
    int i = tx + ii;
    int j = ty + jj;

    float wx = interpolation_func_bicubic((float)i - rx);
    float wy = interpolation_func_bicubic((float)j - ry);
    float w = wx * wy;

    sum += Areadpixel(in, clip_mirror(i, iwidth-1), clip_mirror(j, iheight-1)).z * w;
    weight += w;
  }
  pixel.z = sum/weight;

  pixel = all(isfinite(pixel.xyz)) ? pixel : (float4)0.0f;
  write_ipixel(out, (int2)(x, y), pixel);
}


/* kernels for the lens plugin: lanczos2 interpolation */
kernel void
lens_distort_lanczos2 (read_only image2d_t in,
                       write_only image2d_t out,
                       const int width,
                       const int height,
                       const int iwidth,
                       const int iheight,
                       const int roi_in_x,
                       const int roi_in_y,
                       global float *pi,
                       const int do_nan_checks,
                       const int step,
                       const int gwidth)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  const int kwidth = 2;

  if(x >= width || y >= height) return;

  float4 pixel = (float4)0.0f;

  float rx, ry;
  int tx, ty;
  float sum, weight;
  float2 sum2;
  float ppi[6];
  _lens_grid_sample(pi, x, y, step, gwidth, 6, 6, ppi);

  if(do_nan_checks)
  {
    bool valid = true;

    for(int i = 0; i < 6; i++) valid = valid && isfinite(ppi[i]);

    if(!valid)
    {
      pixel = (float4)0.0f;
      write_imagef (out, (int2)(x, y), pixel);
      return;
    }
  }


  rx = _lens_coord(ppi[0], roi_in_x, iwidth);
  ry = _lens_coord(ppi[1], roi_in_y, iheight);
  tx = floor(rx);
  ty = floor(ry);

  sum = 0.0f;
  weight = 0.0f;
  for(int jj = 1 - kwidth; jj <= kwidth; jj++)
    for(int ii= 1 - kwidth; ii <= kwidth; ii++)
  {
    int i = tx + ii;
    int j = ty + jj;

    float wx = interpolation_func_lanczos(2, (float)i - rx);
    float wy = interpolation_func_lanczos(2, (float)j - ry);
    float w = wx * wy;

    sum += Areadpixel(in, clip_mirror(i, iwidth-1), clip_mirror(j, iheight-1)).x * w;
    weight += w;
  }
  pixel.x = sum/weight;


  rx = _lens_coord(ppi[2], roi_in_x, iwidth);
  ry = _lens_coord(ppi[3], roi_in_y, iheight);
  tx = floor(rx);
  ty = floor(ry);

  sum2 = (float2)0.0f;
  weight = 0.0f;
  for(int jj = 1 - kwidth; jj <= kwidth; jj++)
    for(int ii= 1 - kwidth; ii <= kwidth; ii++)
  {
    int i = tx + ii;
    int j = ty + jj;

    float wx = interpolation_func_lanczos(2, (float)i - rx);
    float wy = interpolation_func_lanczos(2, (float)j - ry);
    float w = wx * wy;

    sum2 += Areadpixel(in, clip_mirror(i, iwidth-1), clip_mirror(j, iheight-1)).yw * w;
    weight += w;
  }
  pixel.yw = sum2/weight;

  rx = _lens_coord(ppi[4], roi_in_x, iwidth);
  ry = _lens_coord(ppi[5], roi_in_y, iheight);
  tx = floor(rx);
  ty = floor(ry);

  sum = 0.0f;
  weight = 0.0f;
  for(int jj = 1 - kwidth; jj <= kwidth; jj++)
    for(int ii= 1 - kwidth; ii <= kwidth; ii++)
  {
    int i = tx + ii;
    int j = ty + jj;

    float wx = interpolation_func_lanczos(2, (float)i - rx);
    float wy = interpolation_func_lanczos(2, (float)j - ry);
    float w = wx * wy;

    sum += Areadpixel(in, clip_mirror(i, iwidth-1), clip_mirror(j, iheight-1)).z * w;
    weight += w;
  }
  pixel.z = sum/weight;

  pixel = all(isfinite(pixel.xyz)) ? pixel : (float4)0.0f;
  write_ipixel(out, (int2)(x, y), pixel);
}


/* kernels for the lens plugin: lanczos3 interpolation */
kernel void
lens_distort_lanczos3 (read_only image2d_t in, write_only image2d_t out, const int width, const int height,
                      const int iwidth, const int iheight, const int roi_in_x, const int roi_in_y, global float *pi,
                      const int do_nan_checks, const int step, const int gwidth)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  const int kwidth = 3;

  if(x >= width || y >= height) return;

  float4 pixel = (float4)0.0f;

  float rx, ry;
  int tx, ty;
  float sum, weight;
  float2 sum2;
  float ppi[6];
  _lens_grid_sample(pi, x, y, step, gwidth, 6, 6, ppi);

  if(do_nan_checks)
  {
    bool valid = true;

    for(int i = 0; i < 6; i++) valid = valid && isfinite(ppi[i]);

    if(!valid)
    {
      pixel = (float4)0.0f;
      write_imagef (out, (int2)(x, y), pixel);
      return;
    }
  }

  rx = _lens_coord(ppi[0], roi_in_x, iwidth);
  ry = _lens_coord(ppi[1], roi_in_y, iheight);
  tx = floor(rx);
  ty = floor(ry);

  sum = 0.0f;
  weight = 0.0f;
  for(int jj = 1 - kwidth; jj <= kwidth; jj++)
    for(int ii= 1 - kwidth; ii <= kwidth; ii++)
  {
    int i = tx + ii;
    int j = ty + jj;

    float wx = interpolation_func_lanczos(3, (float)i - rx);
    float wy = interpolation_func_lanczos(3, (float)j - ry);
    float w = wx * wy;

    sum += Areadpixel(in, clip_mirror(i, iwidth-1), clip_mirror(j, iheight-1)).x * w;
    weight += w;
  }
  pixel.x = sum/weight;

  rx = _lens_coord(ppi[2], roi_in_x, iwidth);
  ry = _lens_coord(ppi[3], roi_in_y, iheight);
  tx = floor(rx);
  ty = floor(ry);

  sum2 = (float2)0.0f;
  weight = 0.0f;
  for(int jj = 1 - kwidth; jj <= kwidth; jj++)
    for(int ii= 1 - kwidth; ii <= kwidth; ii++)
  {
    int i = tx + ii;
    int j = ty + jj;

    float wx = interpolation_func_lanczos(3, (float)i - rx);
    float wy = interpolation_func_lanczos(3, (float)j - ry);
    float w = wx * wy;

    sum2 += Areadpixel(in, clip_mirror(i, iwidth-1), clip_mirror(j, iheight-1)).yw * w;
    weight += w;
  }
  pixel.yw = sum2/weight;

  rx = _lens_coord(ppi[4], roi_in_x, iwidth);
  ry = _lens_coord(ppi[5], roi_in_y, iheight);
  tx = floor(rx);
  ty = floor(ry);

  sum = 0.0f;
  weight = 0.0f;
  for(int jj = 1 - kwidth; jj <= kwidth; jj++)
    for(int ii= 1 - kwidth; ii <= kwidth; ii++)
  {
    int i = tx + ii;
    int j = ty + jj;

    float wx = interpolation_func_lanczos(3, (float)i - rx);
    float wy = interpolation_func_lanczos(3, (float)j - ry);
    float w = wx * wy;

    sum += Areadpixel(in, clip_mirror(i, iwidth-1), clip_mirror(j, iheight-1)).z * w;
    weight += w;
  }
  pixel.z = sum/weight;

  pixel = all(isfinite(pixel.xyz)) ? pixel : (float4)0.0f;
  write_ipixel(out, (int2)(x, y), pixel);
}



/* kernel for the ashift module: bilinear interpolation */
kernel void
ashift_bilinear(read_only image2d_t in, write_only image2d_t out, const int width, const int height,
                const int iwidth, const int iheight, const int2 roi_in, const int2 roi_out,
                const float in_scale, const float out_scale, const float2 clip, global float *homograph)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  float pin[3], pout[3];

  // convert output pixel coordinates to original image coordinates
  pout[0] = roi_out.x + x + clip.x;
  pout[1] = roi_out.y + y + clip.y;
  pout[0] = fp32_div(pout[0], out_scale);
  pout[1] = fp32_div(pout[1], out_scale);
  pout[2] = 1.0f;

  // apply homograph
  for(int i = 0; i < 3; i++)
  {
    pin[i] = 0.0f;
    for(int j = 0; j < 3; j++) pin[i] = fma(homograph[3 * i + j], pout[j], pin[i]);
  }

  // convert to input pixel coordinates
  pin[0] = fp32_div(pin[0], pin[2]);
  pin[1] = fp32_div(pin[1], pin[2]);
  pin[0] *= in_scale;
  pin[1] *= in_scale;
  pin[0] -= roi_in.x;
  pin[1] -= roi_in.y;
  const float rx = pin[0];
  const float ry = pin[1];

  /** don't get output values by interpolation from input image using fast hardware bilinear interpolation
      as that can't mirror and thus doesn't match CPU output.
      Instead do that manually.
  */
  float4 pixel = (float4)0.0f;
  if(rx >= -1.0f && ry >= -1.0f && rx < (float)iwidth && ry < (float)iheight)
  {
    const float fx = floor(rx);
    const float fy = floor(ry);
    const float dx = rx - fx;
    const float dy = ry - fy;

    const int x0 = clip_mirror((int)fx,     iwidth - 1);
    const int x1 = clip_mirror((int)fx + 1, iwidth - 1);
    const int y0 = clip_mirror((int)fy,     iheight - 1);
    const int y1 = clip_mirror((int)fy + 1, iheight - 1);

    pixel = (1.0f - dy) * ((1.0f - dx) * Areadpixel(in, x0, y0) + dx * Areadpixel(in, x1, y0))
                  + dy  * ((1.0f - dx) * Areadpixel(in, x0, y1) + dx * Areadpixel(in, x1, y1));
  }
  write_ipixel(out, (int2)(x, y), pixel);
}

/* kernel for the ashift module: bicubic interpolation */
kernel void
ashift_bicubic (read_only image2d_t in,
                write_only image2d_t out,
                const int width,
                const int height,
                const int iwidth,
                const int iheight,
                const int2 roi_in,
                const int2 roi_out,
                const float in_scale,
                const float out_scale,
                const float2 clip,
                global float *homograph)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  const int kwidth = 2;

  if(x >= width || y >= height) return;

  float pin[3], pout[3];

  // convert output pixel coordinates to original image coordinates
  pout[0] = roi_out.x + x + clip.x;
  pout[1] = roi_out.y + y + clip.y;
  pout[0] = fp32_div(pout[0], out_scale);
  pout[1] = fp32_div(pout[1], out_scale);
  pout[2] = 1.0f;

  // apply homograph
  for(int i = 0; i < 3; i++)
  {
    pin[i] = 0.0f;
    for(int j = 0; j < 3; j++) pin[i] = fma(homograph[3 * i + j], pout[j], pin[i]);
  }

  // convert to input pixel coordinates
  pin[0] = fp32_div(pin[0], pin[2]);
  pin[1] = fp32_div(pin[1], pin[2]);
  pin[0] *= in_scale;
  pin[1] *= in_scale;
  pin[0] -= roi_in.x;
  pin[1] -= roi_in.y;

  // get output values by interpolation from input image
  float rx = pin[0];
  float ry = pin[1];
  int tx = floor(rx);
  int ty = floor(ry);

  float4 pixel = (float4)0.0f;
  for(int jj = 1 - kwidth; jj <= kwidth; jj++)
    for(int ii= 1 - kwidth; ii <= kwidth; ii++)
  {
    const int i = tx + ii;
    const int j = ty + jj;

    float wx = interpolation_func_bicubic((float)i - rx);
    float wy = interpolation_func_bicubic((float)j - ry);
    float w = wx * wy;

    pixel += Areadpixel(in, clip_mirror(i, iwidth-1), clip_mirror(j, iheight-1)) * (float4)w;
  }

  pixel = (tx >= -1 && ty >= -1 && tx < iwidth && ty < iheight)
          ? pixel
          : (float4)0.0f;

  write_ipixel(out, (int2)(x, y), pixel);
}


/* kernel for the ashift module: lanczos2 interpolation */
kernel void
ashift_lanczos2(read_only image2d_t in,
                write_only image2d_t out,
                const int width,
                const int height,
                const int iwidth,
                const int iheight,
                const int2 roi_in,
                const int2 roi_out,
                const float in_scale,
                const float out_scale,
                const float2 clip,
                global float *homograph)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  const int kwidth = 2;

  if(x >= width || y >= height) return;

  float pin[3], pout[3];

  // convert output pixel coordinates to original image coordinates
  pout[0] = roi_out.x + x + clip.x;
  pout[1] = roi_out.y + y + clip.y;
  pout[0] = fp32_div(pout[0], out_scale);
  pout[1] = fp32_div(pout[1], out_scale);
  pout[2] = 1.0f;

  // apply homograph
  for(int i = 0; i < 3; i++)
  {
    pin[i] = 0.0f;
    for(int j = 0; j < 3; j++) pin[i] = fma(homograph[3 * i + j], pout[j], pin[i]);
  }

  // convert to input pixel coordinates
  pin[0] = fp32_div(pin[0], pin[2]);
  pin[1] = fp32_div(pin[1], pin[2]);
  pin[0] *= in_scale;
  pin[1] *= in_scale;
  pin[0] -= roi_in.x;
  pin[1] -= roi_in.y;

  // get output values by interpolation from input image
  float rx = pin[0];
  float ry = pin[1];
  int tx = floor(rx);
  int ty = floor(ry);

  float4 pixel = (float4)0.0f;
  float weight = 0.0f;
  for(int jj = 1 - kwidth; jj <= kwidth; jj++)
    for(int ii= 1 - kwidth; ii <= kwidth; ii++)
  {
    const int i = tx + ii;
    const int j = ty + jj;

    float wx = interpolation_func_lanczos(2, (float)i - rx);
    float wy = interpolation_func_lanczos(2, (float)j - ry);
    float w = wx * wy;

    pixel += Areadpixel(in, clip_mirror(i, iwidth-1), clip_mirror(j, iheight-1)) * (float4)w;
    weight += w;
  }

  pixel = (tx >= -1 && ty >= -1 && tx < iwidth && ty < iheight)
        ? pixel/weight
        : (float4)0.0f;

  write_ipixel(out, (int2)(x, y), pixel);
}


/* kernels for the ashift module: lanczos3 interpolation */
kernel void
ashift_lanczos3(read_only image2d_t in,
                write_only image2d_t out,
                const int width,
                const int height,
                const int iwidth,
                const int iheight,
                const int2 roi_in,
                const int2 roi_out,
                const float in_scale,
                const float out_scale,
                const float2 clip,
                global float *homograph)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  const int kwidth = 3;

  if(x >= width || y >= height) return;

  float pin[3], pout[3];

  // convert output pixel coordinates to original image coordinates
  pout[0] = roi_out.x + x + clip.x;
  pout[1] = roi_out.y + y + clip.y;
  pout[0] = fp32_div(pout[0], out_scale);
  pout[1] = fp32_div(pout[1], out_scale);
  pout[2] = 1.0f;

  // apply homograph
  for(int i = 0; i < 3; i++)
  {
    pin[i] = 0.0f;
    for(int j = 0; j < 3; j++) pin[i] = fma(homograph[3 * i + j], pout[j], pin[i]);
  }

  // convert to input pixel coordinates
  pin[0] = fp32_div(pin[0], pin[2]);
  pin[1] = fp32_div(pin[1], pin[2]);
  pin[0] *= in_scale;
  pin[1] *= in_scale;
  pin[0] -= roi_in.x;
  pin[1] -= roi_in.y;

  // get output values by interpolation from input image
  float rx = pin[0];
  float ry = pin[1];
  int tx = floor(rx);
  int ty = floor(ry);

  float4 pixel = (float4)0.0f;
  float weight = 0.0f;
  for(int jj = 1 - kwidth; jj <= kwidth; jj++)
  {
    for(int ii= 1 - kwidth; ii <= kwidth; ii++)
    {
      const int i = tx + ii;
      const int j = ty + jj;

      float wx = interpolation_func_lanczos(3, (float)i - rx);
      float wy = interpolation_func_lanczos(3, (float)j - ry);
      float w = wx * wy;

      pixel += Areadpixel(in, clip_mirror(i, iwidth-1), clip_mirror(j, iheight-1)) * (float4)w;
      weight += w;
    }
  }
  pixel = (tx >= -1 && ty >= -1 && tx < iwidth && ty < iheight)
          ? pixel/weight
          : (float4)0.0f;

  write_ipixel(out, (int2)(x, y), pixel);
}

float _calc_vignette_spline(const float radius,
                            global float *spline,
                            const int splinesize)
{
  if(radius >= 1.0f) return spline[splinesize-1];

  const float r = radius * (float)(splinesize-1 - 1);
  const int i = (int)r;
  const float frac = r - (float)i;

  const float p0 = spline[i];
  return p0 + (spline[i+1] - p0) * frac;
}

float _interpolate_linear_spline(global float *xi,
                                 global float *yi,
                                 const int ni,
                                 const float x)
{
  if(x < xi[0]) return yi[0];
  for(int i = 1; i < ni; i++)
  {
    if(x >= xi[i - 1] && x <= xi[i])
    {
      const float dydx = (yi[i] - yi[i - 1]) / (xi[i] - xi[i - 1]);
      return yi[i - 1] + (x - xi[i - 1]) * dydx;
    }
  }
  return yi[ni - 1];
}

kernel void md_vignette(read_only image2d_t in,
                        write_only image2d_t out,
                        global float *knots_vig,
                        global float *vig,
                        const int width,
                        const int height,
                        const float w2,
                        const float h2,
                        const float r,
                        const int roix,
                        const int roiy,
                        const int knots)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);
  if(x >= width || y >= height) return;

  const float cx = ((float)(roix + x) - w2);
  const float cy = ((float)(roiy + y) - h2);
  const float4 spline =
    _interpolate_linear_spline(knots_vig, vig, knots, r * hypot(cx, cy));

  float4 pixel  = Areadpixel(in, x, y);
  float mask = pixel.w;
  pixel.x = fp32_div(pixel.x, fmax(1e-4f, spline.x));
  pixel.y = fp32_div(pixel.y, fmax(1e-4f, spline.y));
  pixel.z = fp32_div(pixel.z, fmax(1e-4f, spline.z));
  pixel.w = mask;
  write_imagef (out, (int2)(x, y), pixel);
}

kernel void lens_man_vignette(read_only image2d_t in,
                              write_only image2d_t out,
                              global float *spline,
                              const int width,
                              const int height,
                              const float w2,
                              const float h2,
                              const int roix,
                              const int roiy,
                              const float inv_maxr,
                              const float intensity,
                              const int splinesize,
                              const int vigmask)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);
  if(x >= width || y >= height) return;

  const float dx = ((float)(roix + x) - w2);
  const float dy = ((float)(roiy + y) - h2);
  const float radius = dt_fast_hypot(dx, dy) * inv_maxr;
  const float4 val = fmax(0.0f, intensity * _calc_vignette_spline(radius, spline, splinesize));

  float4 pixel  = Areadpixel(in, x, y);
  const float mask = pixel.w;
  pixel *= (1.0f + val);
  pixel.w = (vigmask) ? val.w : mask;

  write_imagef (out, (int2)(x, y), pixel);
}

kernel void
lens_vignette (read_only image2d_t in,
               write_only image2d_t out,
               const int width,
               const int height,
               global float *pi,
               const int step,
               const int gwidth)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  float4 pixel = Areadpixel(in, x, y);
  float scale[3];
  _lens_grid_sample(pi, x, y, step, gwidth, 4, 3, scale);

  pixel.xyz *= (float3)(scale[0], scale[1], scale[2]);

  write_imagef (out, (int2)(x, y), pixel);
}


#define MAXKNOTS 16
kernel void md_lens_correction(read_only image2d_t in,
                               write_only image2d_t out,
                               global float *knots_dist,
                               global float *cor_rgb,
                               const int owidth,
                               const int oheight,
                               const int iwidth,
                               const int iheight,
                               const float w2,
                               const float h2,
                               const float r,
                               const float scale,
                               const int roix,
                               const int roiy,
                               const int roox,
                               const int rooy,
                               const int knots,
                               const int itor_mode,
                               const int itor_width,
                               const int pass_mode)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);
  if(x >= owidth || y >= oheight) return;

  const float iscale = fp32_div(1.0f, scale);
  const float cx = ((float)(roox + x) - w2) * iscale;
  const float cy = ((float)(rooy + y) - h2) * iscale;
  const float radius = r * dt_fast_hypot(cx, cy);
  const float limw = (float)iwidth - 1.0f;
  const float limh = (float)iheight - 1.0f;
  float output[4];
  for(int c = 0; c < 4; c++)
  {
    const int plane = (c == 3 || pass_mode) ? 1 : c;
    const float dr =
      _interpolate_linear_spline(knots_dist, &cor_rgb[plane * MAXKNOTS], knots, radius);
    const float xs = clamp(dr*cx + w2 - roix, 0.0f, limw);
    const float ys = clamp(dr*cy + h2 - roiy, 0.0f, limh);
    output[c] = interpolation_compute_single(in, itor_mode, itor_width,
                                             xs, ys, iwidth, iheight, c);
  }

  float4 pixel = {output[0], output[1], output[2], output[3]};
  write_ipixel(out, (int2)(x, y), pixel);
}
#undef MAXKNOTS


/* kernel for flip */
__kernel void
flip(read_only image2d_t in,
     write_only image2d_t out,
     const int width,
     const int height,
     const int owidth,
     const int oheight,
     const int orientation)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  // ORIENTATION_FLIP_X = 2
  int ox = (orientation & 2) ? width - x - 1 : x;

  // ORIENTATION_FLIP_Y = 1
  int oy = (orientation & 1) ? height - y - 1 : y;

  // ORIENTATION_SWAP_XY = 4
  if(orientation & 4)
  {
     const int tmp = ox;
     ox = oy;
     oy = tmp;
  }

  if(ox >= 0 && oy >= 0 && ox < owidth && oy < oheight)
  {
    const float4 pixel = readpixel(in, x, y);
    write_imagef(out, (int2)(ox, oy), pixel);
  }
}

float
envelope(const float L)
{
  const float x = clipf(L/100.0f);
  // const float alpha = 2.0f;
  const float beta = 0.6f;
  if(x < beta)
  {
    // return 1.0f-fabsf(x/beta-1.0f)^2
    const float tmp = fabs(x/beta-1.0f);
    return 1.0f-tmp*tmp;
  }
  else
  {
    const float tmp1 = (1.0f-x)/(1.0f-beta);
    const float tmp2 = tmp1*tmp1;
    const float tmp3 = tmp2*tmp1;
    return 3.0f*tmp2 - 2.0f*tmp3;
  }
}

/* kernel for monochrome */
kernel void
monochrome_filter(read_only image2d_t in,
                  write_only image2d_t out,
                  const int width,
                  const int height,
                  const float a,
                  const float b,
                  const float size)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  float4 pixel = Areadpixel(in, x, y);
  pixel.x = 100.0f*dt_fast_expf(-clipf((fsquare(pixel.y - a) + fsquare(pixel.z - b)) / (2.0f * size)));
  write_imagef (out, (int2)(x, y), pixel);
}

kernel void
monochrome(read_only image2d_t in,
           read_only image2d_t base,
           write_only image2d_t out,
           const int width,
           const int height,
           const float a,
           const float b,
           const float size,
           float highlights)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  float4 pixel = Areadpixel(in, x, y);
  float4 basep = Areadpixel(base, x, y);
  float filter  = dt_fast_expf(-clipf((fsquare(pixel.y - a) + fsquare(pixel.z - b)) / (2.0f * size)));
  float tt = envelope(pixel.x);
  float t  = tt + (1.0f-tt)*(1.0f-highlights);
  pixel.x = mix(pixel.x, pixel.x*basep.x/100.0f, t);
  pixel.y = pixel.z = 0.0f;
  write_imagef (out, (int2)(x, y), pixel);
}

/* kernel for the plugin colorout, fast matrix + shaper path only */
kernel void
colorout (read_only image2d_t in,
          write_only image2d_t out,
          const int width,
          const int height,
          global float *mat,
          read_only image2d_t lutr,
          read_only image2d_t lutg,
          read_only image2d_t lutb,
          global const float (*const a)[3], const int hdr_transfer)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  float4 pixel = Areadpixel(in, x, y);
  float XYZ[3], rgb[3];
  float4 xyz = Lab_to_XYZ(pixel);
  XYZ[0] = xyz.x;
  XYZ[1] = xyz.y;
  XYZ[2] = xyz.z;
  for(int i=0;i<3;i++)
  {
    rgb[i] = 0.0f;
    for(int j=0;j<3;j++) rgb[i] += mat[3*i+j]*XYZ[j];
  }
  dt_hdr_encode(rgb, hdr_transfer);
  pixel.x = lerp_lookup_unbounded0(lutr, rgb[0], a[0]);
  pixel.y = lerp_lookup_unbounded0(lutg, rgb[1], a[1]);
  pixel.z = lerp_lookup_unbounded0(lutb, rgb[2], a[2]);
  write_imagef (out, (int2)(x, y), pixel);
}


/* kernel for the plugin colorout: the LittleCMS transform baked into a 3D LUT
   over the Lab range LittleCMS encodes floats in, see _bake_clut() */
kernel void
colorout_clut(read_only image2d_t in,
              write_only image2d_t out,
              const int width,
              const int height,
              global const float *clut,
              const int level,
              const int gamutcheck,
              global const float *gamut,
              const int hdr_transfer)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  const float4 Lab = readpixel(in, x, y);
  const float4 coord = clip4((float4)(Lab.x / 100.0f, (Lab.y + 128.0f) / 255.0f, (Lab.z + 128.0f) / 255.0f, 0.0f));
  float4 pixel = dt_lut3d_tetrahedral(coord, clut, level);
  pixel.w = Lab.w;

  // as _transform_lcms(): out of gamut, or negative, shows as cyan
  if(gamutcheck
     && (dt_lut3d_tetrahedral(coord, gamut, level).x > 0.5f
         || pixel.x < 0.0f || pixel.y < 0.0f || pixel.z < 0.0f))
    pixel = (float4)(0.0f, 1.0f, 1.0f, 0.0f);

  float rgb[3] = { pixel.x, pixel.y, pixel.z };
  dt_hdr_encode(rgb, hdr_transfer);
  pixel.x = rgb[0];
  pixel.y = rgb[1];
  pixel.z = rgb[2];
  write_imagef(out, (int2)(x, y), pixel);
}

/* kernel for the levels plugin */
kernel void
levels (read_only image2d_t in,
        write_only image2d_t out,
        const int width,
        const int height,
        read_only image2d_t lut,
        const float in_low,
        const float in_high,
        const float in_inv_gamma)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  float4 pixel = Areadpixel(in, x, y);
  const float L = pixel.x;
  const float L_in = pixel.x/100.0f;

  if(L_in <= in_low)
  {
    pixel.x = 0.0f;
  }
  else if(L_in >= in_high)
  {
    float percentage = (L_in - in_low) / (in_high - in_low);
    pixel.x = 100.0f * dtcl_pow(percentage, in_inv_gamma);
  }
  else
  {
    float percentage = (L_in - in_low) / (in_high - in_low);
    pixel.x = lookup(lut, percentage);
  }

  if(L_in > 0.01f)
  {
    pixel.y *= pixel.x/L;
    pixel.z *= pixel.x/L;
  }
  else
  {
    pixel.y *= pixel.x;
    pixel.z *= pixel.x;
  }

  write_imagef (out, (int2)(x, y), pixel);
}

/* kernel for the colorzones plugin */
enum
{
  DT_IOP_COLORZONES_L = 0,
  DT_IOP_COLORZONES_C = 1,
  DT_IOP_COLORZONES_h = 2
};


kernel void
colorzones_v3 (read_only image2d_t in,
               write_only image2d_t out,
               const int width,
               const int height,
               const int channel,
               read_only image2d_t table_L,
               read_only image2d_t table_a,
               read_only image2d_t table_b)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  float4 pixel = readpixel(in, x, y);

  const float a = pixel.y;
  const float b = pixel.z;
  const float h = fmod(atan2(b, a) + DT_2PI_F, DT_2PI_F) / DT_2PI_F;
  const float C = dt_fast_hypot(b, a);

  float select = 0.0f;
  float blend = 0.0f;

  switch(channel)
  {
    case DT_IOP_COLORZONES_L:
      select = fmin(1.0f, pixel.x/100.0f);
      break;
    case DT_IOP_COLORZONES_C:
      select = fmin(1.0f, C/128.0f);
      break;
    default:
    case DT_IOP_COLORZONES_h:
      select = h;
      blend = fsquare(1.0f - C/128.0f);
      break;
  }

  const float Lm = (blend * 0.5f + (1.0f-blend)*lookup(table_L, select)) - 0.5f;
  const float hm = (blend * 0.5f + (1.0f-blend)*lookup(table_b, select)) - 0.5f;
  blend *= blend; // saturation isn't as prone to artifacts:
  // const float Cm = 2.0f* (blend*0.5f + (1.0f-blend)*lookup(d->lut[1], select));
  const float Cm = 2.0f * lookup(table_a, select);
  const float L = pixel.x * dtcl_pow(2.0f, 4.0f*Lm);

  pixel.x = L;
  pixel.y = dtcl_cos(DT_2PI_F*(h + hm)) * Cm * C;
  pixel.z = dtcl_sin(DT_2PI_F*(h + hm)) * Cm * C;

  write_imagef (out, (int2)(x, y), pixel);
}

kernel void
colorzones (read_only image2d_t in,
            write_only image2d_t out,
            const int width,
            const int height,
            const int channel,
            read_only image2d_t table_L,
            read_only image2d_t table_C,
            read_only image2d_t table_h)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  float4 pixel = readpixel(in, x, y);

  float4 LCh;
  const float normalize_C = 1.f / (128.0f * M_SQRT2_F);

  LCh = Lab_2_LCH(pixel);

  float select = 0.0f;
  switch(channel)
  {
    case DT_IOP_COLORZONES_L:
      select = LCh.x * 0.01f;
      break;
    case DT_IOP_COLORZONES_C:
      select = LCh.y * normalize_C;
      break;
    case DT_IOP_COLORZONES_h:
    default:
      select = LCh.z;
      break;
  }
  select = clipf(select);

  LCh.x *= dtcl_pow(2.0f, 4.0f * (lookup(table_L, select) - .5f));
  LCh.y *= 2.f * lookup(table_C, select);
  LCh.z += lookup(table_h, select) - .5f;

  pixel.xyz = LCH_2_Lab(LCh).xyz;

  write_imagef (out, (int2)(x, y), pixel);
}


/* kernel for the zonesystem plugin */
kernel void
zonesystem (read_only image2d_t in,
            write_only image2d_t out,
            const int width,
            const int height,
            const int size,
            global float *zonemap_offset,
            global float *zonemap_scale)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  float4 pixel = readpixel(in, x, y);

  const float rzscale = (float)(size-1)/100.0f;
  const int rz = clamp((int)(pixel.x*rzscale), 0, size-2);
  const float zs = ((rz > 0) ? (zonemap_offset[rz]/pixel.x) : 0) + zonemap_scale[rz];

  pixel.xyz *= zs;

  write_imagef (out, (int2)(x, y), pixel);
}




/* kernel to fill an image with a color (for the borders plugin). */
kernel void
borders_fill (write_only image2d_t out,
              const int left,
              const int top,
              const int width,
              const int height,
              const float4 color)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x < left || y < top) return;
  if(x >= width + left || y >= height + top) return;

  write_imagef (out, (int2)(x, y), color);
}


/* kernel for the overexposed plugin. */
typedef enum dt_clipping_preview_mode_t
{
  DT_CLIPPING_PREVIEW_GAMUT = 0,
  DT_CLIPPING_PREVIEW_ANYRGB = 1,
  DT_CLIPPING_PREVIEW_LUMINANCE = 2,
  DT_CLIPPING_PREVIEW_SATURATION = 3
} dt_clipping_preview_mode_t;

kernel void
overexposed (read_only image2d_t in,
             write_only image2d_t out,
             read_only image2d_t tmp,
             const int width,
             const int height,
             const float lower,
             const float upper,
             const float4 lower_color,
             const float4 upper_color,
             constant dt_colorspaces_iccprofile_info_cl_t *profile_info,
             read_only image2d_t lut,
             const int use_work_profile,
             dt_clipping_preview_mode_t mode)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  float4 pixel = readpixel(in, x, y);
  float4 pixel_tmp = readpixel(tmp, x, y);

  if(mode == DT_CLIPPING_PREVIEW_ANYRGB)
  {
    if(pixel_tmp.x >= upper || pixel_tmp.y >= upper || pixel_tmp.z >= upper)
      pixel.xyz = upper_color.xyz;

    else if(pixel_tmp.x <= lower && pixel_tmp.y <= lower && pixel_tmp.z <= lower)
      pixel.xyz = lower_color.xyz;

  }
  else if(mode == DT_CLIPPING_PREVIEW_GAMUT && use_work_profile)
  {
    const float luminance = get_rgb_matrix_luminance(pixel, profile_info, profile_info->matrix_in, lut);

    if(luminance >= upper)
    {
      pixel.xyz = upper_color.xyz;
    }
    else if(luminance <= lower)
    {
      pixel.xyz = lower_color.xyz;
    }
    else
    {
      float4 saturation = { 0.f, 0.f, 0.f, 0.f};
      saturation = pixel_tmp - (float4)luminance;
      saturation = dtcl_sqrt(saturation * saturation / ((float4)(luminance * luminance) + pixel_tmp * pixel_tmp));

      if(saturation.x > upper || saturation.y > upper || saturation.z > upper ||
         pixel_tmp.x >= upper || pixel_tmp.y >= upper || pixel_tmp.z >= upper)
        pixel.xyz = upper_color.xyz;

      else if(pixel_tmp.x <= lower && pixel_tmp.y <= lower && pixel_tmp.z <= lower)
        pixel.xyz = lower_color.xyz;
    }
  }
  else if(mode == DT_CLIPPING_PREVIEW_LUMINANCE && use_work_profile)
  {
    const float luminance = get_rgb_matrix_luminance(pixel, profile_info, profile_info->matrix_in, lut);

    if(luminance >= upper)
      pixel.xyz = upper_color.xyz;

    else if(luminance <= lower)
      pixel.xyz = lower_color.xyz;
  }
  else if(mode == DT_CLIPPING_PREVIEW_SATURATION && use_work_profile)
  {
    const float luminance = get_rgb_matrix_luminance(pixel, profile_info, profile_info->matrix_in, lut);

    if(luminance < upper && luminance > lower)
    {
      float4 saturation = { 0.f, 0.f, 0.f, 0.f};
      saturation = pixel_tmp - (float4)luminance;
      saturation = dtcl_sqrt(saturation * saturation / ((float4)(luminance * luminance) + pixel_tmp * pixel_tmp));

      if(saturation.x > upper || saturation.y > upper || saturation.z > upper ||
         pixel_tmp.x >= upper || pixel_tmp.y >= upper || pixel_tmp.z >= upper)
        pixel.xyz = upper_color.xyz;

      else if(pixel_tmp.x <= lower && pixel_tmp.y <= lower && pixel_tmp.z <= lower)
        pixel.xyz = lower_color.xyz;
    }
  }

  write_imagef (out, (int2)(x, y), pixel);
}


/* kernel for the rawoverexposed plugin. */
kernel void
rawoverexposed_mark_cfa (read_only image2d_t in,
                         write_only image2d_t out,
                         global float *pi,
                         const int width,
                         const int height,
                         read_only image2d_t raw,
                         const int raw_width,
                         const int raw_height,
                         const unsigned int filters,
                         global const unsigned char (*const xtrans)[6],
                         global unsigned int *threshold,
                         global float4 *colors)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  const int piwidth = 2*width;
  global float *ppi = pi + mad24(y, piwidth, 2*x);

  const int raw_x = ppi[0];
  const int raw_y = ppi[1];

  if(raw_x < 0 || raw_y < 0 || raw_x >= raw_width || raw_y >= raw_height) return;

  const uint raw_pixel = read_imageui(raw, sampleri, (int2)(raw_x, raw_y)).x;

  const int c = fcol(raw_y, raw_x, filters, xtrans);
  if(raw_pixel < threshold[c]) return;

  float4 pixel = fmax(0.0f, Areadpixel(in, x, y));
  const float4 color = colors[c & 3];

  // cfa color
  pixel.xyz = color.xyz;
  write_imagef (out, (int2)(x, y), pixel);
}

kernel void
rawoverexposed_mark_solid (read_only image2d_t in,
                           write_only image2d_t out,
                           global float *pi,
                           const int width,
                           const int height,
                           read_only image2d_t raw,
                           const int raw_width,
                           const int raw_height,
                           const unsigned int filters,
                           global const unsigned char (*const xtrans)[6],
                           global unsigned int *threshold,
                           const float4 solid_color)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  const int piwidth = 2*width;
  global float *ppi = pi + mad24(y, piwidth, 2*x);

  const int raw_x = ppi[0];
  const int raw_y = ppi[1];

  if(raw_x < 0 || raw_y < 0 || raw_x >= raw_width || raw_y >= raw_height) return;

  const uint raw_pixel = read_imageui(raw, sampleri, (int2)(raw_x, raw_y)).x;

  const int c = fcol(raw_y, raw_x, filters, xtrans);
  if(raw_pixel < threshold[c]) return;

  float4 pixel = fmax(0.0f, Areadpixel(in, x, y));

  // solid color
  pixel.xyz = solid_color.xyz;

  write_imagef (out, (int2)(x, y), pixel);
}

kernel void
rawoverexposed_falsecolor (read_only image2d_t in,
                           write_only image2d_t out,
                           global float *pi,
                           const int width,
                           const int height,
                           read_only image2d_t raw,
                           const int raw_width,
                           const int raw_height,
                           const unsigned int filters,
                           global const unsigned char (*const xtrans)[6],
                           global unsigned int *threshold)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  const int piwidth = 2*width;
  global float *ppi = pi + mad24(y, piwidth, 2*x);

  const int raw_x = ppi[0];
  const int raw_y = ppi[1];

  if(raw_x < 0 || raw_y < 0 || raw_x >= raw_width || raw_y >= raw_height) return;

  const uint raw_pixel = read_imageui(raw, sampleri, (int2)(raw_x, raw_y)).x;

  const int c = fcol(raw_y, raw_x, filters, xtrans);
  if(raw_pixel < threshold[c]) return;

  float4 pixel = fmax(0.0f, Areadpixel(in, x, y));
  if(c == 2)      pixel.z = 0.0f;
  else if(c == 1) pixel.y = 0.0f;
  else if(c == 0) pixel.x = 0.0f;

  write_imagef (out, (int2)(x, y), pixel);
}


/* kernel for the lowlight plugin. */
kernel void
lowlight (read_only image2d_t in,
          write_only image2d_t out,
          const int width,
          const int height,
          const float4 XYZ_sw,
          read_only image2d_t lut)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);

  if(x >= width || y >= height) return;

  const float c = 0.5f;
  const float threshold = 0.01f;

  float V;
  float w;

  float4 pixel = readpixel(in, x, y);

  float4 XYZ = Lab_to_XYZ(pixel);

  // calculate scotopic luminance
  if (XYZ.x > threshold)
  {
    // normal flow
    V = XYZ.y * ( 1.33f * ( 1.0f + (XYZ.y+XYZ.z)/XYZ.x) - 1.68f );
  }
  else
  {
    // low red flow, avoids "snow" on dark noisy areas
    V = XYZ.y * ( 1.33f * ( 1.0f + (XYZ.y+XYZ.z)/threshold) - 1.68f );
  }

  // scale using empiric coefficient and fit inside limits
  V = clipf(c*V);

  // blending coefficient from curve
  w = lookup(lut, pixel.x/100.0f);

  XYZ = w * XYZ + (1.0f - w) * V * XYZ_sw;

  pixel = XYZ_to_Lab(XYZ);

  write_imagef (out, (int2)(x, y), pixel);
}


/* kernel for the contrast lightness saturation module */
kernel void
colisa (read_only image2d_t in,
        write_only image2d_t out,
        unsigned int width,
        unsigned int height,
        const float saturation,
        read_only image2d_t ctable,
        constant float *ca,
        read_only image2d_t ltable,
        constant float *la)
{
  const unsigned int x = get_global_id(0);
  const unsigned int y = get_global_id(1);

  if(x >= width || y >= height) return;

  float4 i = readpixel(in, x, y);
  float4 o;

  o.x = lookup_unbounded(ctable, i.x/100.0f, ca);
  o.x = lookup_unbounded(ltable, o.x/100.0f, la);
  o.y = i.y*saturation;
  o.z = i.z*saturation;
  o.w = i.w;

  write_imagef(out, (int2)(x, y), o);
}

/* kernel for the unbreak input profile module - gamma version */

kernel void
profilegamma (read_only image2d_t in,
              write_only image2d_t out,
              int width,
              int height,
              read_only image2d_t table,
              constant float *ta)
{
  const unsigned int x = get_global_id(0);
  const unsigned int y = get_global_id(1);

  if(x >= width || y >= height) return;

  float4 i = readpixel(in, x, y);
  float4 o;

  o.x = lookup_unbounded(table, i.x, ta);
  o.y = lookup_unbounded(table, i.y, ta);
  o.z = lookup_unbounded(table, i.z, ta);
  o.w = i.w;

  write_imagef(out, (int2)(x, y), o);
}

/* kernel for the unbreak input profile module - log version */
kernel void
profilegamma_log (read_only image2d_t in,
                  write_only image2d_t out,
                  int width,
                  int height,
                  const float dynamic_range,
                  const float shadows_range,
                  const float grey)
{
  const unsigned int x = get_global_id(0);
  const unsigned int y = get_global_id(1);

  if(x >= width || y >= height) return;

  float4 i = readpixel(in, x, y);
  const float4 noise = pow((float4)2.0f, (float4)-16.0f);
  const float4 dynamic4 = dynamic_range;
  const float4 shadows4 = shadows_range;
  const float4 grey4 = grey;

  float4 o;

  o = (i < noise) ? noise : i / grey4;
  o = (dtcl_log2(o) - shadows4) / dynamic4;
  o = (o < noise) ? noise : o;
  i.xyz = o.xyz;

  write_imagef(out, (int2)(x, y), i);
}

/* kernel for the interpolation resample helper */
kernel void
interpolation_resample (read_only image2d_t in,
                        write_only image2d_t out,
                        const int width,
                        const int height,
                        const global int *hmeta,
                        const global int *vmeta,
                        const global int *hlength,
                        const global int *vlength,
                        const global int *hindex,
                        const global int *vindex,
                        const global float *hkernel,
                        const global float *vkernel,
                        const int htaps,
                        const int vtaps,
                        local float *lkernel,
                        local int *lindex,
                        local float4 *buffer)
{
  const int x = get_global_id(0);
  const int yi = get_global_id(1);
  const int ylsz = get_local_size(1);
  const int xlid = get_local_id(0);
  const int ylid = get_local_id(1);
  const int y = yi / vtaps;
  const int iy = yi % vtaps;

  // Initialize resampling indices
  const int xm = min(x, width - 1);
  const int ym = min(y, height - 1);
  const int hlidx = hmeta[xm*3];   // H(orizontal) L(ength) I(n)d(e)x
  const int hkidx = hmeta[xm*3+1]; // H(orizontal) K(ernel) I(n)d(e)x
  const int hiidx = hmeta[xm*3+2]; // H(orizontal) I(ndex) I(n)d(e)x
  const int vlidx = vmeta[ym*3];   // V(ertical) L(ength) I(n)d(e)x
  const int vkidx = vmeta[ym*3+1]; // V(ertical) K(ernel) I(n)d(e)x
  const int viidx = vmeta[ym*3+2]; // V(ertical) I(ndex) I(n)d(e)x

  const int hl = hlength[hlidx];   // H(orizontal) L(ength)
  const int vl = vlength[vlidx];   // V(ertical) L(ength)

  // generate local copy of horizontal index field and kernel
  for(int n = 0; n <= htaps/ylsz; n++)
  {
    int k = mad24(n, ylsz, ylid);
    if(k >= hl) continue;
    lindex[k] = hindex[hiidx+k];
    lkernel[k] = hkernel[hkidx+k];
  }

  barrier(CLK_LOCAL_MEM_FENCE);

  // horizontal convolution kernel; store intermediate result in local buffer
  if(x < width && y < height)
  {
    const int yvalid = iy < vl;

    const int yy = yvalid ? vindex[viidx+iy] : -1;

    float4 vpixel = (float4)0.0f;

    for (int ix = 0; ix < hl && yvalid; ix++)
    {
      const int xx = lindex[ix];
      float4 hpixel = readpixel(in, xx, yy);
      vpixel += hpixel * lkernel[ix];
    }

    buffer[ylid] = yvalid ? vpixel * vkernel[vkidx+iy] : (float4)0.0f;
  }
  else
    buffer[ylid] = (float4)0.0f;

  barrier(CLK_LOCAL_MEM_FENCE);

  // recursively reduce local buffer (vertical convolution kernel)
  for(int offset = vtaps / 2; offset > 0; offset >>= 1)
  {
    if (iy < offset)
    {
      buffer[ylid] += buffer[ylid + offset];
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  }

  // store final result
  if (iy == 0 && x < width && y < height)
  {
    write_ipixel(out, (int2)(x, y), buffer[ylid]);
  }
}

/* the statistics of dt_color_picker_helper() in src/common/color_picker.c,
   reduced on the device. Each work item reduces one tile of the box and
   writes 16 floats: the sums, minima, maxima and pixel counts of the four
   channels. The host adds up the tiles, so only those come back. */
#define DT_PICK_PLAIN 0
#define DT_PICK_LCH 1
#define DT_PICK_HSL 2
#define DT_PICK_JZCZHZ 3

typedef struct dt_picker_stats_t
{
  float4 sum, low, high, count;
} dt_picker_stats_t;

static inline void _picker_stats_store(const dt_picker_stats_t *s, global float *stats, const int item)
{
  vstore4(s->sum, 4 * item, stats);
  vstore4(s->low, 4 * item + 1, stats);
  vstore4(s->high, 4 * item + 2, stats);
  vstore4(s->count, 4 * item + 3, stats);
}

// blur_2D_Bspline() at one scale and without clipping, as the picker's
// denoise does over the whole image, with the edges clamped to the image
static inline float4 _picker_denoised(read_only image2d_t in, const int x, const int y,
                                      const int width, const int height)
{
  const float filter[5] = { 1.0f / 16.0f, 4.0f / 16.0f, 6.0f / 16.0f, 4.0f / 16.0f, 1.0f / 16.0f };
  float4 acc = (float4)0.0f;
  for(int ii = 0; ii < 5; ii++)
  {
    const int xx = clamp(x + ii - 2, 0, width - 1);
    float4 v = (float4)0.0f;
    for(int jj = 0; jj < 5; jj++)
      v += filter[jj] * read_imagef(in, sampleri, (int2)(xx, clamp(y + jj - 2, 0, height - 1)));
    acc += filter[ii] * v;
  }
  return acc;
}

// rgb_to_JzCzhz() of color_picker.c: the profile's matrix to XYZ D50, then
// the Bradford adaptation to D65
static inline float4 _picker_jzczhz(const float4 rgb,
                                    constant const dt_colorspaces_iccprofile_info_cl_t *const profile_info,
                                    read_only image2d_t lut)
{
  const float4 D50 = rgb_matrix_to_xyz(rgb, profile_info, profile_info->matrix_in, lut);
  const float4 D65 = (float4)(0.9555766f * D50.x - 0.0230393f * D50.y + 0.0631636f * D50.z,
                              -0.0282895f * D50.x + 1.0099416f * D50.y + 0.0210077f * D50.z,
                              0.0122982f * D50.x - 0.0204830f * D50.y + 1.3299098f * D50.z,
                              0.0f);
  const float4 JzAzBz = XYZ_to_JzAzBz(D65);
  const float h = atan2(JzAzBz.z, JzAzBz.y) / DT_2PI_F;
  return (float4)(JzAzBz.x, hypot(JzAzBz.y, JzAzBz.z), h >= 0.0f ? h : 1.0f + h, 0.0f);
}

kernel void
color_picker_stats_4ch(read_only image2d_t in,
                       global float *stats,
                       const int width,
                       const int height,
                       const int bx0,
                       const int by0,
                       const int bx1,
                       const int by1,
                       const int tile_w,
                       const int tile_h,
                       const int tiles_x,
                       const int tiles_y,
                       const int mode,
                       const int denoise,
                       constant const dt_colorspaces_iccprofile_info_cl_t *const profile_info,
                       read_only image2d_t lut)
{
  const int tx = get_global_id(0);
  const int ty = get_global_id(1);
  if(tx >= tiles_x || ty >= tiles_y) return;

  const int x0 = bx0 + tx * tile_w;
  const int y0 = by0 + ty * tile_h;
  const int x1 = min(x0 + tile_w, bx1);
  const int y1 = min(y0 + tile_h, by1);

  dt_picker_stats_t s = { (float4)0.0f, (float4)FLT_MAX, (float4)(-FLT_MAX), (float4)0.0f };
  for(int y = y0; y < y1; y++)
    for(int x = x0; x < x1; x++)
    {
      const float4 p = denoise ? _picker_denoised(in, x, y, width, height)
                               : read_imagef(in, sampleri, (int2)(x, y));
      float4 pick = p;
      if(mode != DT_PICK_PLAIN)
      {
        pick = mode == DT_PICK_LCH ? Lab_2_LCH(p)
             : mode == DT_PICK_HSL ? RGB_2_HSL(p)
             : _picker_jzczhz(p, profile_info, lut);
        // the hue shifted by half a turn, for sensible minima and maxima,
        // as _update_stats_4ch() does
        pick.w = pick.z < 0.5f ? pick.z + 0.5f : pick.z - 0.5f;
      }
      s.sum += pick;
      s.low = fmin(s.low, pick);
      s.high = fmax(s.high, pick);
      s.count += 1.0f;
    }
  _picker_stats_store(&s, stats, ty * tiles_x + tx);
}

// a Bayer or X-Trans mosaic, per CFA color
kernel void
color_picker_stats_mosaic(read_only image2d_t in,
                          global float *stats,
                          const int bx0,
                          const int by0,
                          const int bx1,
                          const int by1,
                          const int tile_w,
                          const int tile_h,
                          const int tiles_x,
                          const int tiles_y,
                          const int rx,
                          const int ry,
                          const unsigned int filters,
                          global const unsigned char (*const xtrans)[6])
{
  const int tx = get_global_id(0);
  const int ty = get_global_id(1);
  if(tx >= tiles_x || ty >= tiles_y) return;

  const int x0 = bx0 + tx * tile_w;
  const int y0 = by0 + ty * tile_h;
  const int x1 = min(x0 + tile_w, bx1);
  const int y1 = min(y0 + tile_h, by1);

  float sum[4] = { 0.0f, 0.0f, 0.0f, 0.0f };
  float low[4] = { FLT_MAX, FLT_MAX, FLT_MAX, FLT_MAX };
  float high[4] = { -FLT_MAX, -FLT_MAX, -FLT_MAX, -FLT_MAX };
  float count[4] = { 0.0f, 0.0f, 0.0f, 0.0f };
  for(int y = y0; y < y1; y++)
    for(int x = x0; x < x1; x++)
    {
      const int c = fcol(y + ry, x + rx, filters, xtrans);
      const float v = read_imagef(in, sampleri, (int2)(x, y)).x;
      sum[c] += v;
      low[c] = fmin(low[c], v);
      high[c] = fmax(high[c], v);
      count[c] += 1.0f;
    }
  const dt_picker_stats_t s = { vload4(0, sum), vload4(0, low), vload4(0, high), vload4(0, count) };
  _picker_stats_store(&s, stats, ty * tiles_x + tx);
}

/* kernel for the interpolation copy helper */
kernel void
interpolation_copy(read_only image2d_t dev_in,
                   write_only image2d_t dev_out,
                   const int owidth,
                   const int oheight,
                   const int iwidth,
                   const int iheight,
                   const int dx,
                   const int dy)
{
  const int ocol = get_global_id(0);
  const int orow = get_global_id(1);

  if(ocol >= owidth || orow >= oheight) return;

  float4 pix = (float4)0.0f;

  const int irow = orow + dy;
  const int icol = ocol + dx;

  if(irow < iheight && irow >= 0 && icol < iwidth && icol >= 0)
  {
    pix = Areadpixel(dev_in, icol, irow);
  }
  write_imagef(dev_out, (int2)(ocol, orow), pix);
}


// process() of src/iop/hotpixels.c: each work item writes its own sensel,
// either copied or replaced by the largest of its four like-colored
// neighbors. mode 0 is Bayer, with the neighbors 2 sensels away; 1 is a
// monochrome sensor, with the direct neighbors; 2 is X-Trans, with the
// neighbors from offsets[(row % 6) * 6 + col % 6], as (x, y) pairs
kernel void
hotpixels(read_only image2d_t in,
          write_only image2d_t out,
          const int width,
          const int height,
          const float threshold,
          const float multiplier,
          const int min_neighbours,
          const int mode,
          global const int4 *offsets,
          global unsigned int *fixed)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);
  if(x >= width || y >= height) return;

  const float pixel = Areadsingle(in, x, y);
  float o = pixel;
  const int margin = mode == 1 ? 1 : 2;
  if(pixel > threshold && x >= margin && y >= margin && x < width - margin && y < height - margin)
  {
    int2 n[4];
    if(mode == 2)
    {
      const int k = (y % 6) * 6 + x % 6;
      const int4 a = offsets[2 * k];
      const int4 b = offsets[2 * k + 1];
      n[0] = a.xy;
      n[1] = a.zw;
      n[2] = b.xy;
      n[3] = b.zw;
    }
    else
    {
      const int s = mode == 1 ? 1 : 2;
      n[0] = (int2)(-s, 0);
      n[1] = (int2)(0, -s);
      n[2] = (int2)(s, 0);
      n[3] = (int2)(0, s);
    }

    const float mid = pixel * multiplier;
    int count = 0;
    float maxin = 0.0f;
    for(int i = 0; i < 4; i++)
    {
      const float other = Areadsingle(in, x + n[i].x, y + n[i].y);
      if(mid > other)
      {
        count++;
        maxin = fmax(maxin, other);
      }
    }
    if(count >= min_neighbours)
    {
      o = maxin;
      atomic_inc(fixed);
    }
  }
  write_imagef(out, (int2)(x, y), o);
}


// the module histograms of dt_histogram_helper() in src/common/histogram.c,
// counted on the device so that only the bins are read back. bins holds 4
// counts per bin as on the CPU and must be zeroed by the caller. With
// use_local each work group counts into lhist (3 * bins_count) first and
// adds that to bins; without, for bin counts too large for local memory,
// every pixel goes to bins directly.
#define DT_HISTOGRAM_RGB 0
#define DT_HISTOGRAM_RGB_COMPENSATED 1
#define DT_HISTOGRAM_LAB 2
#define DT_HISTOGRAM_LCH 3

// dt_ioppr_compensate_middle_grey(): the lightness of a gray of value x
static inline float _histogram_compensate(const float x,
                                          constant const dt_colorspaces_iccprofile_info_cl_t *const profile_info,
                                          read_only image2d_t lut)
{
  const float4 xyz = rgb_matrix_to_xyz((float4)(x, x, x, 0.0f), profile_info, profile_info->matrix_in, lut);
  return XYZ_to_Lab(xyz).x * 0.01f;
}

kernel void
histogram_collect(read_only image2d_t in,
                  const int x0,
                  const int y0,
                  const int x1,
                  const int y1,
                  const int bins_count,
                  const int mode,
                  const int use_local,
                  constant const dt_colorspaces_iccprofile_info_cl_t *const profile_info,
                  read_only image2d_t lut,
                  global unsigned int *bins,
                  local unsigned int *lhist)
{
  const int x = x0 + get_global_id(0);
  const int y = y0 + get_global_id(1);
  const int l = mad24((int)get_local_id(1), (int)get_local_size(0), (int)get_local_id(0));
  const int lsz = mul24((int)get_local_size(0), (int)get_local_size(1));

  if(use_local)
  {
    for(int i = l; i < 3 * bins_count; i += lsz) lhist[i] = 0;
    barrier(CLK_LOCAL_MEM_FENCE);
  }

  if(x < x1 && y < y1)
  {
    const float4 p = read_imagef(in, sampleri, (int2)(x, y));
    const float max_bin = bins_count - 1;
    float4 b;
    if(mode == DT_HISTOGRAM_RGB)
      b = max_bin * p;
    else if(mode == DT_HISTOGRAM_RGB_COMPENSATED)
      b = max_bin * (float4)(_histogram_compensate(p.x, profile_info, lut),
                             _histogram_compensate(p.y, profile_info, lut),
                             _histogram_compensate(p.z, profile_info, lut), 0.0f);
    else if(mode == DT_HISTOGRAM_LAB)
      b = (float4)(max_bin / 100.0f, max_bin / 256.0f, max_bin / 256.0f, 0.0f)
          * (p + (float4)(0.0f, 128.0f, 128.0f, 0.0f));
    else
      b = (float4)(max_bin / 100.0f, max_bin / (128.0f * M_SQRT2_F), max_bin, 0.0f) * Lab_2_LCH(p);

    const int4 bin = convert_int4(clamp(b, 0.0f, max_bin));
    if(use_local)
    {
      atomic_inc(&lhist[bin.x]);
      atomic_inc(&lhist[bins_count + bin.y]);
      atomic_inc(&lhist[2 * bins_count + bin.z]);
    }
    else
    {
      atomic_inc(&bins[4 * bin.x]);
      atomic_inc(&bins[4 * bin.y + 1]);
      atomic_inc(&bins[4 * bin.z + 2]);
    }
  }
  if(!use_local) return;
  barrier(CLK_LOCAL_MEM_FENCE);

  for(int i = l; i < 3 * bins_count; i += lsz)
  {
    const unsigned int count = lhist[i];
    if(count) atomic_add(&bins[4 * (i % bins_count) + i / bins_count], count);
  }
}
// clang-format off
// modelines: These editor modelines have been set for all relevant files by tools/update_modelines.py
// vim: shiftwidth=2 expandtab tabstop=2 cindent
// kate: tab-indents: off; indent-width 2; replace-tabs on; indent-mode cstyle; remove-trailing-spaces modified;
// clang-format on