tiling.c

/*
    This file is part of darktable,
    Copyright (C) 2011-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 "develop/tiling.h"
#include "common/opencl.h"
#include "control/control.h"
#include "develop/blend.h"
#include "develop/pixelpipe.h"

#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <strings.h>
#include <unistd.h>


/* parameter RESERVE for extended roi_in sizes due to inaccuracies when doing
   roi_out -> roi_in estimations.
   Needs to be increased if tiling fails due to insufficient buffer sizes. */
#define RESERVE 5

#ifdef HAVE_OPENCL
/* greatest common divisor */
static unsigned _gcd(unsigned a, unsigned b)
{
  unsigned t;
  while(b != 0)
  {
    t = b;
    b = a % b;
    a = t;
  }
  return MAX(a, 1);
}

/* least common multiple */
static unsigned _lcm(const unsigned a, const unsigned b)
{
  return (((unsigned long)a * b) / _gcd(a, b));
}
#endif

static inline int _align_up(const int n, const int a)
{
  return n + a - (n % a);
}
static inline int _align_down(const int n, const int a)
{
  return n - (n % a);
}
static inline int _align_close(const int n, const int a)
{
  const int off = n % a;
  const int shift = (off > a/2) ? a - off : -off;
  return n + shift;
}

/*
  _maximum_number_tiles is the assumed maximum sane number of tiles
  if during tiling this number is exceeded darktable assumes that tiling is not possible and falls back
  to untiled processing - with all system memory limits taking full effect.
  For huge images like stitched panos the user might choose resourcelevel="unrestricted", in that
  case the allowed number of tiles is practically unlimited
*/
static inline int _maximum_number_tiles()
{
  return (darktable.dtresources.level == 3) ? 0x40000000 : 10000;
}

static inline void _print_roi(const dt_iop_roi_t *roi, const char *label)
{
  dt_print(DT_DEBUG_TILING | DT_DEBUG_VERBOSE,"     {%5d %5d ->%5d %5d (%5dx%5d)  %.6f } %s",
           roi->x, roi->y, roi->x + roi->width, roi->y + roi->height,
           roi->width, roi->height, roi->scale, label);
}


static double _nm_fitness(double x[], void *rest[])
{
  struct dt_iop_module_t *self = (struct dt_iop_module_t *)rest[0];
  struct dt_dev_pixelpipe_iop_t *piece = (struct dt_dev_pixelpipe_iop_t *)rest[1];
  struct dt_iop_roi_t *iroi = (struct dt_iop_roi_t *)rest[2];
  struct dt_iop_roi_t *oroi = (struct dt_iop_roi_t *)rest[3];

  dt_iop_roi_t oroi_test = *oroi;
  oroi_test.x = x[0] * piece->iwidth;
  oroi_test.y = x[1] * piece->iheight;
  oroi_test.width = x[2] * piece->iwidth;
  oroi_test.height = x[3] * piece->iheight;

  dt_iop_roi_t iroi_probe = *iroi;
  self->modify_roi_in(self, piece, &oroi_test, &iroi_probe);

  double fitness = 0.0;

  fitness += (double)(iroi_probe.x - iroi->x) * (iroi_probe.x - iroi->x);
  fitness += (double)(iroi_probe.y - iroi->y) * (iroi_probe.y - iroi->y);
  fitness += (double)(iroi_probe.width - iroi->width) * (iroi_probe.width - iroi->width);
  fitness += (double)(iroi_probe.height - iroi->height) * (iroi_probe.height - iroi->height);

  return fitness;
}


/* We use a Nelder-Mead simplex algorithm based on an implementation of Michael F. Hutt.
   It is covered by the following copyright notice: */
/*
 * Program: nmsimplex.c
 * Author : Michael F. Hutt
 * http://www.mikehutt.com
 * 11/3/97
 *
 * An implementation of the Nelder-Mead simplex method.
 *
 * Copyright (c) 1997-2011 <Michael F. Hutt>
 *
 * Permission is hereby granted, free of charge, to any person obtaining
 * a copy of this software and associated documentation files (the
 * "Software"), to deal in the Software without restriction, including
 * without limitation the rights to use, copy, modify, merge, publish,
 * distribute, sublicense, and/or sell copies of the Software, and to
 * permit persons to whom the Software is furnished to do so, subject to
 * the following conditions:
 *
 * The above copyright notice and this permission notice shall be
 * included in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 *
 */

#define ALPHA 1.0   /* reflection coefficient */
#define BETA 0.5    /* contraction coefficient */
#define GAMMA 2.0   /* expansion coefficient */

// #define DT_TILING_DEBUG

static int _simplex(double (*objfunc)(double[], void *[]),
                    const double start[],
                    const int n,
                    const double EPSILON,
                    const double scale,
                    const int maxiter,
                    void (*constrain)(double[], int n),
                    void *rest[])
{

  int vs; /* vertex with smallest value */
  int vh; /* vertex with next smallest value */
  int vg; /* vertex with largest value */

  int i, j = 0, m, row;
  int itr; /* track the number of iterations */

  double **v;    /* holds vertices of simplex */
  double pn, qn; /* values used to create initial simplex */
  double *f;     /* value of function at each vertex */
  double fr;     /* value of function at reflection point */
  double fe;     /* value of function at expansion point */
  double fc;     /* value of function at contraction point */
  double *vr;    /* reflection - coordinates */
  double *ve;    /* expansion - coordinates */
  double *vc;    /* contraction - coordinates */
  double *vm;    /* centroid - coordinates */

  double fsum, favg, s, cent;

  /* dynamically allocate arrays */

  /* allocate the rows of the arrays */
  v = (double **)malloc(sizeof(double *) * (n + 1));
  f = (double *)malloc(sizeof(double) * (n + 1));
  vr = (double *)malloc(sizeof(double) * n);
  ve = (double *)malloc(sizeof(double) * n);
  vc = (double *)malloc(sizeof(double) * n);
  vm = (double *)malloc(sizeof(double) * n);

  /* allocate the columns of the arrays */
  for(i = 0; i <= n; i++)
  {
    v[i] = (double *)malloc(sizeof(double) * n);
  }

  /* create the initial simplex */
  /* assume one of the vertices is 0,0 */

  pn = scale * (sqrt(n + 1) - 1 + n) / (n * sqrt(2));
  qn = scale * (sqrt(n + 1) - 1) / (n * sqrt(2));

  for(i = 0; i < n; i++)
  {
    v[0][i] = start[i];
  }

  for(i = 1; i <= n; i++)
  {
    for(j = 0; j < n; j++)
    {
      if(i - 1 == j)
      {
        v[i][j] = pn + start[j];
      }
      else
      {
        v[i][j] = qn + start[j];
      }
    }
  }

  if(constrain != NULL)
  {
    constrain(v[j], n);
  }
  /* find the initial function values */
  for(j = 0; j <= n; j++)
  {
    f[j] = objfunc(v[j], rest);
  }

#ifdef DT_TILING_DEBUG
  /* print out the initial values */
  printf ("Initial Values\n");
  for(j = 0; j <= n; j++)
  {
    for(i = 0; i < n; i++)
    {
      printf ("%f %f\n", v[j][i], f[j]);
    }
  }
#endif

  /* begin the main loop of the minimization */
  for(itr = 1; itr <= maxiter; itr++)
  {
    /* find the index of the largest value */
    vg = 0;
    for(j = 0; j <= n; j++)
    {
      if(f[j] > f[vg])
      {
        vg = j;
      }
    }

    /* find the index of the smallest value */
    vs = 0;
    for(j = 0; j <= n; j++)
    {
      if(f[j] < f[vs])
      {
        vs = j;
      }
    }

    /* find the index of the second largest value */
    vh = vs;
    for(j = 0; j <= n; j++)
    {
      if(f[j] > f[vh] && f[j] < f[vg])
      {
        vh = j;
      }
    }

    /* calculate the centroid */
    for(j = 0; j <= n - 1; j++)
    {
      cent = 0.0;
      for(m = 0; m <= n; m++)
      {
        if(m != vg)
        {
          cent += v[m][j];
        }
      }
      vm[j] = cent / n;
    }

    /* reflect vg to new vertex vr */
    for(j = 0; j <= n - 1; j++)
    {
      /*vr[j] = (1+ALPHA)*vm[j] - ALPHA*v[vg][j]; */
      vr[j] = vm[j] + ALPHA * (vm[j] - v[vg][j]);
    }
    if(constrain != NULL)
    {
      constrain(vr, n);
    }
    fr = objfunc(vr, rest);

    if(fr < f[vh] && fr >= f[vs])
    {
      for(j = 0; j <= n - 1; j++)
      {
        v[vg][j] = vr[j];
      }
      f[vg] = fr;
    }

    /* investigate a step further in this direction */
    if(fr < f[vs])
    {
      for(j = 0; j <= n - 1; j++)
      {
        /*ve[j] = GAMMA*vr[j] + (1-GAMMA)*vm[j]; */
        ve[j] = vm[j] + GAMMA * (vr[j] - vm[j]);
      }
      if(constrain != NULL)
      {
        constrain(ve, n);
      }
      fe = objfunc(ve, rest);

      /* by making fe < fr as opposed to fe < f[vs],
         Rosenbrocks function takes 63 iterations as opposed
         to 64 when using double variables. */

      if(fe < fr)
      {
        for(j = 0; j <= n - 1; j++)
        {
          v[vg][j] = ve[j];
        }
        f[vg] = fe;
      }
      else
      {
        for(j = 0; j <= n - 1; j++)
        {
          v[vg][j] = vr[j];
        }
        f[vg] = fr;
      }
    }

    /* check to see if a contraction is necessary */
    if(fr >= f[vh])
    {
      if(fr < f[vg] && fr >= f[vh])
      {
        /* perform outside contraction */
        for(j = 0; j <= n - 1; j++)
        {
          /*vc[j] = BETA*v[vg][j] + (1-BETA)*vm[j]; */
          vc[j] = vm[j] + BETA * (vr[j] - vm[j]);
        }
        if(constrain != NULL)
        {
          constrain(vc, n);
        }
        fc = objfunc(vc, rest);
      }
      else
      {
        /* perform inside contraction */
        for(j = 0; j <= n - 1; j++)
        {
          /*vc[j] = BETA*v[vg][j] + (1-BETA)*vm[j]; */
          vc[j] = vm[j] - BETA * (vm[j] - v[vg][j]);
        }
        if(constrain != NULL)
        {
          constrain(vc, n);
        }
        fc = objfunc(vc, rest);
      }


      if(fc < f[vg])
      {
        for(j = 0; j <= n - 1; j++)
        {
          v[vg][j] = vc[j];
        }
        f[vg] = fc;
      }
      /* at this point the contraction is not successful,
         we must halve the distance from vs to all the
         vertices of the simplex and then continue.
         10/31/97 - modified to account for ALL vertices.
       */
      else
      {
        for(row = 0; row <= n; row++)
        {
          if(row != vs)
          {
            for(j = 0; j <= n - 1; j++)
            {
              v[row][j] = v[vs][j] + (v[row][j] - v[vs][j]) / 2.0;
            }
          }
        }
        if(constrain != NULL)
        {
          constrain(v[vg], n);
        }
        f[vg] = objfunc(v[vg], rest);
        if(constrain != NULL)
        {
          constrain(v[vh], n);
        }
        f[vh] = objfunc(v[vh], rest);
      }
    }

#ifdef DT_TILING_DEBUG
    /* print out the value at each iteration */
    printf ("Iteration %d\n", itr);
    for(j = 0; j <= n; j++)
    {
      for(i = 0; i < n; i++)
      {
        printf ("%f %f\n", v[j][i], f[j]);
      }
    }
#endif

    /* test for convergence */
    fsum = 0.0;
    for(j = 0; j <= n; j++)
    {
      fsum += f[j];
    }
    favg = fsum / (n + 1);
    s = 0.0;
    for(j = 0; j <= n; j++)
    {
      s += pow((f[j] - favg), 2.0) / (n);
    }
    s = sqrt(s);
    if(s < EPSILON) break;
  }
  /* end main loop of the minimization */

  /* find the index of the smallest value */
  vs = 0;
  for(j = 0; j <= n; j++)
  {
    if(f[j] < f[vs])
    {
      vs = j;
    }
  }

#ifdef DT_TILING_DEBUG
  printf ("The minimum was found at\n");
  for(j = 0; j < n; j++)
  {
    printf ("%e\n", v[vs][j]);
    start[j] = v[vs][j];
  }
  double min = objfunc (v[vs], rest);
  printf ("Function value at minimum %f\n", min);
  k++;
  printf ("%d Function Evaluations\n", k);
  printf ("%d Iterations through program\n", itr);
#endif

  free(f);
  free(vr);
  free(ve);
  free(vc);
  free(vm);
  for(i = 0; i <= n; i++)
  {
    free(v[i]);
  }
  free(v);
  return itr;
}


static gboolean _nm_fit_output_to_input_roi(const dt_iop_module_t *self,
                                            const dt_dev_pixelpipe_iop_t *piece,
                                            const dt_iop_roi_t *iroi,
                                            dt_iop_roi_t *oroi,
                                            const int delta)
{
  void *rest[4] = { (void *)self, (void *)piece, (void *)iroi, (void *)oroi };
  double start[4] = { (float)oroi->x / piece->iwidth, (float)oroi->y / piece->iheight,
                      (float)oroi->width / piece->iwidth, (float)oroi->height / piece->iheight };
  const double epsilon = (double)delta / MIN(piece->iwidth, piece->iheight);
  const int maxiter = 1000;

  const int iter = _simplex(_nm_fitness, start, 4, epsilon, 1.0, maxiter, NULL, rest);

  dt_print(DT_DEBUG_TILING | DT_DEBUG_VERBOSE,
           "[_nm_fit_output_to_input_roi] _simplex: %d, delta: %d, epsilon: %f",
           iter, delta, epsilon);

  oroi->x = start[0] * piece->iwidth;
  oroi->y = start[1] * piece->iheight;
  oroi->width = start[2] * piece->iwidth;
  oroi->height = start[3] * piece->iheight;

  return (iter <= maxiter);
}



/* find a matching oroi_full by probing start value of oroi and get corresponding input roi into iroi_probe.
   We search in two steps. first by a simplicistic iterative search which will succeed in most cases.
   If this does not converge, we do a downhill simplex (nelder-mead) fitting */
static gboolean _fit_output_to_input_roi(dt_iop_module_t *self,
                                         dt_dev_pixelpipe_iop_t *piece,
                                         dt_iop_roi_t *iroi,
                                         dt_iop_roi_t *oroi,
                                         const int delta,
                                         int iter)
{
  dt_iop_roi_t iroi_probe = *iroi;
  dt_iop_roi_t save_oroi = *oroi;

  // try to go the easy way. this works in many cases where output is
  // just like input, only scaled down
  self->modify_roi_in(self, piece, oroi, &iroi_probe);
  while((abs((int)iroi_probe.x - (int)iroi->x) > delta || abs((int)iroi_probe.y - (int)iroi->y) > delta
         || abs((int)iroi_probe.width - (int)iroi->width) > delta
         || abs((int)iroi_probe.height - (int)iroi->height) > delta) && iter > 0)
  {
    _print_roi(&iroi_probe, "tile iroi_probe");
    _print_roi(oroi, "tile oroi old");

    oroi->x += (iroi->x - iroi_probe.x) * oroi->scale / iroi->scale;
    oroi->y += (iroi->y - iroi_probe.y) * oroi->scale / iroi->scale;
    oroi->width += (iroi->width - iroi_probe.width) * oroi->scale / iroi->scale;
    oroi->height += (iroi->height - iroi_probe.height) * oroi->scale / iroi->scale;

    _print_roi(oroi, "tile oroi new");

    self->modify_roi_in(self, piece, oroi, &iroi_probe);
    iter--;
  }

  if(iter > 0) return TRUE;

  *oroi = save_oroi;

  // simplicistic approach did not converge.
  // try simplex downhill fitting now.
  // it's crucial that we have a good starting point in oroi, else this
  // will not converge as well.
  return _nm_fit_output_to_input_roi(self, piece, iroi, oroi, delta);
}


/* simple tiling algorithm for roi_in == roi_out, i.e. for pixel to pixel modules/operations */
static void _default_process_tiling_ptp(dt_iop_module_t *self,
                                        dt_dev_pixelpipe_iop_t *piece,
                                        const void *const ivoid,
                                        void *const ovoid,
                                        const dt_iop_roi_t *const roi_in,
                                        const dt_iop_roi_t *const roi_out,
                                        const int in_bpp)
{
  void *input = NULL;
  void *output = NULL;

  dt_iop_buffer_dsc_t dsc;
  self->output_format(self, piece->pipe, piece, &dsc);
  const int out_bpp = dt_iop_buffer_dsc_to_bpp(&dsc);

  const int ipitch = roi_in->width * in_bpp;
  const int opitch = roi_out->width * out_bpp;
  const int max_bpp = MAX(in_bpp, out_bpp);

  /* get tiling requirements of module */
  dt_develop_tiling_t tiling = { 0 };
  tiling.factor_cl = tiling.maxbuf_cl = -1;
  self->tiling_callback(self, piece, roi_in, roi_out, &tiling);
  if(tiling.factor_cl < 0) tiling.factor_cl = tiling.factor;
  if(tiling.maxbuf_cl < 0) tiling.maxbuf_cl = tiling.maxbuf;

  /* tiling really does not make sense in these cases. standard process() is not better or worse than we are
   */
  if((tiling.factor < 2.2f)
     && (tiling.overhead < 0.2f * roi_in->width * roi_in->height * max_bpp))
  {
    dt_print(DT_DEBUG_PIPE | DT_DEBUG_TILING,
             "[default_process_tiling_ptp] [%s]  no need to use tiling for module '%s%s' "
             "as no real memory saving to be expected",
             dt_dev_pixelpipe_type_to_str(piece->pipe->type), self->op, dt_iop_get_instance_id(self));
    goto fallback;
  }

  /* calculate optimal size of tiles */
  float available = dt_get_available_pipe_mem(piece->pipe);
  assert(available >= 500.0f * 1024.0f * 1024.0f);
  /* correct for size of ivoid and ovoid which are needed on top of tiling */
  available = fmaxf(available - ((float)roi_out->width * roi_out->height * out_bpp)
                   - ((float)roi_in->width * roi_in->height * in_bpp) - tiling.overhead,
                   0);

  /* we ignore the above value if singlebuffer_limit (is defined and) is higher than available/tiling.factor.
     this will mainly allow tiling for modules with high and "unpredictable" memory demand which is
     reflected in high values of tiling.factor (take bilateral noise reduction as an example). */
  float singlebuffer = dt_get_singlebuffer_mem();
  const float factor = fmaxf(tiling.factor, 1.0f);
  const float maxbuf = fmaxf(tiling.maxbuf, 1.0f);
  singlebuffer = fmaxf(available / factor, singlebuffer);

  int width = roi_in->width;
  int height = roi_in->height;

  /* shrink tile size in case it would exceed singlebuffer size */
  if((float)width * height * max_bpp * maxbuf > singlebuffer)
  {
    const float scale = singlebuffer / ((float)width * height * max_bpp * maxbuf);

    /* TODO: can we make this more efficient to minimize total overlap between tiles? */
    if(width < height && scale >= 0.333f)
    {
      height = floorf(height * scale);
    }
    else if(height <= width && scale >= 0.333f)
    {
      width = floorf(width * scale);
    }
    else
    {
      width = floorf(width * sqrtf(scale));
      height = floorf(height * sqrtf(scale));
    }
    dt_print(DT_DEBUG_TILING | DT_DEBUG_VERBOSE,
             "[default_process_tiling_ptp] buffer exceeds singlebuffer, corrected to %dx%d",
             width, height);
  }

  /* make sure we have a reasonably effective tile dimension. if not try square tiles */
  if(3 * tiling.overlap > width || 3 * tiling.overlap > height)
  {
    width = height = floorf(sqrtf((float)width * height));
    dt_print(DT_DEBUG_TILING | DT_DEBUG_VERBOSE,
             "[default_process_tiling_roi] use squares because of overlap, corrected to %dx%d",
             width, height);
  }

  /* Alignment rules: we need to make sure that alignment requirements of module are fulfilled.
     Modules will report alignment requirements via align within tiling_callback().
     We guarantee alignment by selecting image width/height and overlap accordingly. For a tile width/height
     that is identical to image width/height no special alignment is needed. */

  const unsigned int align = tiling.align;
  assert(align != 0);

  /* properly align tile width and height by making them smaller if needed */
  if(width < roi_in->width) width = (width / align) * align;
  if(height < roi_in->height) height = (height / align) * align;

  /* also make sure that overlap follows alignment rules by making it wider when needed */
  const int overlap = tiling.overlap % align != 0 ? (tiling.overlap / align + 1) * align
                                                    : tiling.overlap;

  /* calculate effective tile size */
  const int tile_wd = width - 2 * overlap > 0 ? width - 2 * overlap : 1;
  const int tile_ht = height - 2 * overlap > 0 ? height - 2 * overlap : 1;

  /* calculate number of tiles */
  const int tiles_x = width < roi_in->width ? ceilf(roi_in->width / (float)tile_wd) : 1;
  const int tiles_y = height < roi_in->height ? ceilf(roi_in->height / (float)tile_ht) : 1;

  /* sanity check: don't run wild on too many tiles */
  if(tiles_x * tiles_y > _maximum_number_tiles())
  {
    dt_print(DT_DEBUG_PIPE | DT_DEBUG_TILING,
             "[default_process_tiling_ptp] [%s] gave up tiling for module '%s%s'. too many tiles: %d x %d",
             dt_dev_pixelpipe_type_to_str(piece->pipe->type),
             self->op, dt_iop_get_instance_id(self), tiles_x, tiles_y);
    goto error;
  }

  /* reserve input and output buffers for tiles */
  input = dt_alloc_aligned((size_t)width * height * in_bpp);
  if(input == NULL)
  {
    dt_print(DT_DEBUG_TILING,
             "[default_process_tiling_ptp] [%s] could not alloc input buffer for module '%s%s'",
             dt_dev_pixelpipe_type_to_str(piece->pipe->type), self->op, dt_iop_get_instance_id(self));
    goto error;
  }
  output = dt_alloc_aligned((size_t)width * height * out_bpp);
  if(output == NULL)
  {
    dt_print(DT_DEBUG_TILING,
             "[default_process_tiling_ptp] [%s]  could not alloc output buffer for module '%s%s'",
             dt_dev_pixelpipe_type_to_str(piece->pipe->type), self->op, dt_iop_get_instance_id(self));
    goto error;
  }

  /* store processed_maximum to be re-used and aggregated */
  dt_aligned_pixel_t processed_maximum_saved;
  dt_aligned_pixel_t processed_maximum_new = { 1.0f };
  for_four_channels(k) processed_maximum_saved[k] = piece->pipe->dsc.processed_maximum[k];

  piece->pipe->tiling = TRUE;
  dt_print_pipe(DT_DEBUG_PIPE | DT_DEBUG_TILING,
                        "default *tiled* ptp", piece->pipe, piece->module, DT_DEVICE_CPU, roi_in, roi_out,
                        "%dx%d tiles, size=%dx%d, overlap=%d",
                        tiles_x, tiles_y, tile_wd, tile_ht, overlap);

  /* iterate over tiles */
  for(size_t tx = 0; tx < tiles_x; tx++)
  {
    const size_t wd = tx * tile_wd + width > roi_in->width ? roi_in->width - tx * tile_wd : width;
    for(size_t ty = 0; ty < tiles_y; ty++)
    {
      const size_t ht = ty * tile_ht + height > roi_in->height ? roi_in->height - ty * tile_ht : height;

      /* no need to process end-tiles that are smaller than the total overlap area */
      const gboolean skipped = (wd <= 2 * overlap && tx > 0) || (ht <= 2 * overlap && ty > 0);

      /* origin and region of effective part of tile, which we want to store later */
      size_t origin[2] = { 0, 0 };
      size_t region[2] = { wd, ht };

      /* roi_in and roi_out for process_cl on subbuffer */
      dt_iop_roi_t iroi = { roi_in->x + tx * tile_wd, roi_in->y + ty * tile_ht, wd, ht, roi_in->scale };
      dt_iop_roi_t oroi = { roi_out->x + tx * tile_wd, roi_out->y + ty * tile_ht, wd, ht, roi_out->scale };

      /* offsets of tile into ivoid and ovoid */
      const size_t ioffs = (ty * tile_ht) * ipitch + (tx * tile_wd) * in_bpp;
      size_t ooffs = (ty * tile_ht) * opitch + (tx * tile_wd) * out_bpp;

      dt_print_pipe(DT_DEBUG_TILING,
               skipped ? "  tile ptp skipped" : "  tile ptp", piece->pipe, piece->module, DT_DEVICE_CPU, &iroi, &oroi,
               "tile (%zu,%zu)", tx, ty);
      if(skipped) continue;

      /* prepare input tile buffer */
      DT_OMP_FOR()
      for(size_t j = 0; j < ht; j++)
        memcpy((char *)input + j * wd * in_bpp, (char *)ivoid + ioffs + j * ipitch, (size_t)wd * in_bpp);

      /* take original processed_maximum as starting point */
      for(int k = 0; k < 4; k++) piece->pipe->dsc.processed_maximum[k] = processed_maximum_saved[k];
      dt_dev_prepare_piece_cfa(piece, &iroi);

      /* call process() of module */
      self->process(self, piece, input, output, &iroi, &oroi);

      /* aggregate resulting processed_maximum */
      /* TODO: check if there really can be differences between tiles and take
               appropriate action (calculate minimum, maximum, average, ...?) */
      for(int k = 0; k < 4; k++)
      {
        if(tx + ty > 0 && fabs(processed_maximum_new[k] - piece->pipe->dsc.processed_maximum[k]) > 1.0e-6f)
          dt_print(DT_DEBUG_TILING,
                   "[default_process_tiling_ptp] [%s] processed_maximum[%d] differs between tiles in module '%s%s'",
                   dt_dev_pixelpipe_type_to_str(piece->pipe->type), k,
                   self->op, dt_iop_get_instance_id(self));
        processed_maximum_new[k] = piece->pipe->dsc.processed_maximum[k];
      }

      /* correct origin and region of tile for overlap.
         make sure that we only copy back the "good" part. */
      if(tx > 0)
      {
        origin[0] += overlap;
        region[0] -= overlap;
        ooffs += (size_t)overlap * out_bpp;
      }
      if(ty > 0)
      {
        origin[1] += overlap;
        region[1] -= overlap;
        ooffs += (size_t)overlap * opitch;
      }

      /* copy "good" part of tile to output buffer */
      DT_OMP_FOR(shared(origin, region))
      for(size_t j = 0; j < region[1]; j++)
        memcpy((char *)ovoid + ooffs + j * opitch,
               (char *)output + ((j + origin[1]) * wd + origin[0]) * out_bpp, (size_t)region[0] * out_bpp);
    }
  }

  /* copy back final processed_maximum */
  for(int k = 0; k < 4; k++) piece->pipe->dsc.processed_maximum[k] = processed_maximum_new[k];

  dt_free_align(input);
  dt_free_align(output);
  piece->pipe->tiling = FALSE;
  return;

error:
  dt_control_log(_("tiling failed for module '%s%s'. the output most likely will be OK, but you might want to check."),
    self->op, dt_iop_get_instance_id(self));
// fall through

fallback:
  dt_free_align(input);
  dt_free_align(output);
  piece->pipe->tiling = FALSE;
  dt_print(DT_DEBUG_TILING,
           "[default_process_tiling_ptp] [%s] fall back to standard processing for module '%s%s'",
           dt_dev_pixelpipe_type_to_str(piece->pipe->type), self->op, dt_iop_get_instance_id(self));
  self->process(self, piece, ivoid, ovoid, roi_in, roi_out);
  return;
}



/* more elaborate tiling algorithm for roi_in != roi_out: slower than the ptp variant,
   more tiles and larger overlap */
static void _default_process_tiling_roi(dt_iop_module_t *self,
                                        dt_dev_pixelpipe_iop_t *piece,
                                        const void *const ivoid,
                                        void *const ovoid,
                                        const dt_iop_roi_t *const roi_in,
                                        const dt_iop_roi_t *const roi_out,
                                        const int in_bpp)
{
  void *input = NULL;
  void *output = NULL;

  dt_iop_buffer_dsc_t dsc;
  self->output_format(self, piece->pipe, piece, &dsc);
  const int out_bpp = dt_iop_buffer_dsc_to_bpp(&dsc);

  const int ipitch = roi_in->width * in_bpp;
  const int opitch = roi_out->width * out_bpp;
  const int max_bpp = MAX(in_bpp, out_bpp);

  float fullscale = fmaxf(roi_in->scale / roi_out->scale, sqrtf(((float)roi_in->width * roi_in->height)
                                                              / ((float)roi_out->width * roi_out->height)));

  /* inaccuracy for roi_in elements in roi_out -> roi_in calculations */
  const int delta = ceilf(fullscale);

  /* estimate for additional (space) requirement in buffer dimensions due to inaccuracies */
  const int inacc = RESERVE * delta;

  /* get tiling requirements of module */
  dt_develop_tiling_t tiling = { 0 };
  tiling.factor_cl = tiling.maxbuf_cl = -1;
  self->tiling_callback(self, piece, roi_in, roi_out, &tiling);
  if(tiling.factor_cl < 0) tiling.factor_cl = tiling.factor;
  if(tiling.maxbuf_cl < 0) tiling.maxbuf_cl = tiling.maxbuf;

  /* tiling really does not make sense in these cases. standard process() is not better or worse than we are
   */
  if((tiling.factor < 2.2f && tiling.overhead < 0.2f * roi_in->width * roi_in->height * max_bpp))
  {
    dt_print(DT_DEBUG_TILING,
             "[default_process_tiling_roi] [%s] no need to use tiling for module "
             "'%s%s' as no memory saving is expected",
             dt_dev_pixelpipe_type_to_str(piece->pipe->type), self->op, dt_iop_get_instance_id(self));
    goto fallback;
  }

  /* calculate optimal size of tiles */
  float available = dt_get_available_pipe_mem(piece->pipe);
  assert(available >= 500.0f * 1024.0f * 1024.0f);
  /* correct for size of ivoid and ovoid which are needed on top of tiling */
  available = fmaxf(available - ((float)roi_out->width * roi_out->height * out_bpp)
                   - ((float)roi_in->width * roi_in->height * in_bpp) - tiling.overhead,
                   0);

  /* we ignore the above value if singlebuffer_limit (is defined and) is higher than available/tiling.factor.
     this will mainly allow tiling for modules with high and "unpredictable" memory demand which is
     reflected in high values of tiling.factor (take bilateral noise reduction as an example). */
  float singlebuffer = dt_get_singlebuffer_mem();
  const float factor = fmaxf(tiling.factor, 1.0f);
  const float maxbuf = fmaxf(tiling.maxbuf, 1.0f);
  singlebuffer = fmaxf(available / factor, singlebuffer);

  int width = MAX(roi_in->width, roi_out->width);
  int height = MAX(roi_in->height, roi_out->height);

  /* Alignment rules: we need to make sure that alignment requirements of module are fulfilled.
     Modules will report alignment requirements via align within tiling_callback(). */

  const unsigned int align = tiling.align;
  assert(align != 0);

  /* shrink tile size in case it would exceed singlebuffer size */
  if((float)width * height * max_bpp * maxbuf > singlebuffer)
  {
    const float scale = singlebuffer / ((float)width * height * max_bpp * maxbuf);

    /* TODO: can we make this more efficient to minimize total overlap between tiles? */
    if(width < height && scale >= 0.333f)
    {
      height = _align_down((int)floorf(height * scale), align);
    }
    else if(height <= width && scale >= 0.333f)
    {
      width = _align_down((int)floorf(width * scale), align);
    }
    else
    {
      width = _align_down((int)floorf(width * sqrtf(scale)), align);
      height = _align_down((int)floorf(height * sqrtf(scale)), align);
    }
    dt_print(DT_DEBUG_TILING | DT_DEBUG_VERBOSE,
             "[default_process_tiling_roi] [%s] buffer exceeds singlebuffer, corrected to %dx%d",
             dt_dev_pixelpipe_type_to_str(piece->pipe->type), width, height);
  }

  /* make sure we have a reasonably effective tile dimension. if not try square tiles */
  if(3 * tiling.overlap > width || 3 * tiling.overlap > height)
  {
    width = height = _align_down((int)floorf(sqrtf((float)width * height)), align);
    dt_print(DT_DEBUG_TILING | DT_DEBUG_VERBOSE,
             "[default_process_tiling_roi] [%s] use squares because of overlap, corrected to %dx%d",
             dt_dev_pixelpipe_type_to_str(piece->pipe->type), width, height);
  }

  /* make sure that overlap follows alignment rules by making it wider when needed.
     overlap_in needs to be aligned, overlap_out is only here to calculate output buffer size */
  const int overlap_in = _align_up(tiling.overlap, align);
  const int overlap_out = ceilf((float)overlap_in / fullscale);

  int tiles_x = 1, tiles_y = 1;

  /* calculate number of tiles taking the larger buffer (input or output) as a guiding one.
     normally it is roi_in > roi_out; but let's be prepared */
  if(roi_in->width > roi_out->width)
    tiles_x = width < roi_in->width
                  ? ceilf((float)roi_in->width / (float)MAX(width - 2 * overlap_in - inacc, 1))
                  : 1;
  else
    tiles_x = width < roi_out->width ? ceilf((float)roi_out->width / (float)MAX(width - 2 * overlap_out, 1))
                                     : 1;

  if(roi_in->height > roi_out->height)
    tiles_y = height < roi_in->height
                  ? ceilf((float)roi_in->height / (float)MAX(height - 2 * overlap_in - inacc, 1))
                  : 1;
  else
    tiles_y = height < roi_out->height
                  ? ceilf((float)roi_out->height / (float)MAX(height - 2 * overlap_out, 1))
                  : 1;

  /* sanity check: don't run wild on too many tiles */
  if(tiles_x * tiles_y > _maximum_number_tiles())
  {
    dt_print(DT_DEBUG_TILING,
             "[default_process_tiling_roi] [%s] gave up tiling for module '%s%s'. too many tiles: %d x %d",
             dt_dev_pixelpipe_type_to_str(piece->pipe->type),
             self->op, dt_iop_get_instance_id(self), tiles_x, tiles_y);
    goto error;
  }


  /* calculate tile width and height excl. overlap (i.e. the good part) for output.
     values are important for all following processing steps. */
  const int tile_wd = _align_up(
      roi_out->width % tiles_x == 0 ? roi_out->width / tiles_x : roi_out->width / tiles_x + 1, align);
  const int tile_ht = _align_up(
      roi_out->height % tiles_y == 0 ? roi_out->height / tiles_y : roi_out->height / tiles_y + 1, align);

  dt_print_pipe(DT_DEBUG_PIPE | DT_DEBUG_TILING,
                        "process *tiled* roi", piece->pipe, piece->module, DT_DEVICE_CPU, roi_in, roi_out,
                        "%dx%d tiles, size=%dx%d",
                        tiles_x, tiles_y, tile_wd, tile_ht);

  /* store processed_maximum to be re-used and aggregated */
  dt_aligned_pixel_t processed_maximum_saved;
  dt_aligned_pixel_t processed_maximum_new = { 1.0f };
  for_four_channels(k) processed_maximum_saved[k] = piece->pipe->dsc.processed_maximum[k];

  /* iterate over tiles */
  for(size_t tx = 0; tx < tiles_x; tx++)
    for(size_t ty = 0; ty < tiles_y; ty++)
    {
      piece->pipe->tiling = TRUE;

      /* the output dimensions of the good part of this specific tile */
      const size_t wd = (tx + 1) * tile_wd > roi_out->width ? (size_t)roi_out->width - tx * tile_wd : tile_wd;
      const size_t ht = (ty + 1) * tile_ht > roi_out->height ? (size_t)roi_out->height - ty * tile_ht : tile_ht;

      /* roi_in and roi_out of good part: oroi_good easy to calculate based on number and dimension of tile.
         iroi_good is calculated by modify_roi_in() of respective module */
      dt_iop_roi_t iroi_good = { roi_in->x  + tx * tile_wd, roi_in->y  + ty * tile_ht, wd, ht, roi_in->scale };
      dt_iop_roi_t oroi_good = { roi_out->x + tx * tile_wd, roi_out->y + ty * tile_ht, wd, ht, roi_out->scale };

      self->modify_roi_in(self, piece, &oroi_good, &iroi_good);

      /* clamp iroi_good to not exceed roi_in */
      iroi_good.x = MAX(iroi_good.x, roi_in->x);
      iroi_good.y = MAX(iroi_good.y, roi_in->y);
      iroi_good.width = MIN(iroi_good.width, roi_in->width + roi_in->x - iroi_good.x);
      iroi_good.height = MIN(iroi_good.height, roi_in->height + roi_in->y - iroi_good.y);

      _print_roi(&iroi_good, "tile iroi_good");
      _print_roi(&oroi_good, "tile oroi_good");

      /* now we need to calculate full region of this tile: increase input roi to take care of overlap
         requirements
         and alignment and add additional delta to correct for possible rounding errors in modify_roi_in()
         -> generates first estimate of iroi_full */
      const int x_in = iroi_good.x;
      const int y_in = iroi_good.y;
      const int width_in = iroi_good.width;
      const int height_in = iroi_good.height;
      const int new_x_in = MAX(_align_close(x_in - overlap_in - delta, align), roi_in->x);
      const int new_y_in = MAX(_align_close(y_in - overlap_in - delta, align), roi_in->y);
      const int new_width_in = MIN(_align_up(width_in + overlap_in + delta + (x_in - new_x_in), align),
                                    roi_in->width + roi_in->x - new_x_in);
      const int new_height_in = MIN(_align_up(height_in + overlap_in + delta + (y_in - new_y_in), align),
                                     roi_in->height + roi_in->y - new_y_in);

      /* iroi_full based on calculated numbers and dimensions. oroi_full just set as a starting point for the
       * following iterative search */
      dt_iop_roi_t iroi_full = { new_x_in, new_y_in, new_width_in, new_height_in, iroi_good.scale };
      dt_iop_roi_t oroi_full = oroi_good; // a good starting point for optimization

      _print_roi(&iroi_full, "tile iroi_full before optimization");
      _print_roi(&oroi_full, "tile oroi_full before optimization");

      /* try to find a matching oroi_full */
      if(!_fit_output_to_input_roi(self, piece, &iroi_full, &oroi_full, delta, 10))
      {
        dt_print(DT_DEBUG_TILING,
                 "[default_process_tiling_roi] [%s] can not handle requested roi's. "
                 "tiling for module '%s%s' not possible",
                 dt_dev_pixelpipe_type_to_str(piece->pipe->type), self->op, dt_iop_get_instance_id(self));
        goto error;
      }

      _print_roi(&iroi_full, "tile iroi_full after optimization");
      _print_roi(&oroi_full, "tile oroi_full after optimization");

      /* make sure that oroi_full at least covers the range of oroi_good.
         this step is needed due to the possibility of rounding errors */
      oroi_full.x = MIN(oroi_full.x, oroi_good.x);
      oroi_full.y = MIN(oroi_full.y, oroi_good.y);
      oroi_full.width = MAX(oroi_full.width, oroi_good.x + oroi_good.width - oroi_full.x);
      oroi_full.height = MAX(oroi_full.height, oroi_good.y + oroi_good.height - oroi_full.y);

      /* clamp oroi_full to not exceed roi_out */
      oroi_full.x = MAX(oroi_full.x, roi_out->x);
      oroi_full.y = MAX(oroi_full.y, roi_out->y);
      oroi_full.width = MIN(oroi_full.width, roi_out->width + roi_out->x - oroi_full.x);
      oroi_full.height = MIN(oroi_full.height, roi_out->height + roi_out->y - oroi_full.y);

      /* calculate final iroi_full */
      self->modify_roi_in(self, piece, &oroi_full, &iroi_full);

      /* clamp iroi_full to not exceed roi_in */
      iroi_full.x = MAX(iroi_full.x, roi_in->x);
      iroi_full.y = MAX(iroi_full.y, roi_in->y);
      iroi_full.width = MIN(iroi_full.width, roi_in->width + roi_in->x - iroi_full.x);
      iroi_full.height = MIN(iroi_full.height, roi_in->height + roi_in->y - iroi_full.y);

      _print_roi(&iroi_full, "tile iroi_full final");
      _print_roi(&oroi_full, "tile oroi_full final");

      /* offsets of tile into ivoid and ovoid */
      const size_t ioffs = ((size_t)iroi_full.y - roi_in->y)  * ipitch + ((size_t)iroi_full.x - roi_in->x) * in_bpp;
            size_t ooffs = ((size_t)oroi_good.y - roi_out->y) * opitch + ((size_t)oroi_good.x - roi_out->x) * out_bpp;

      dt_print_pipe(DT_DEBUG_TILING,
               "  tile roi", piece->pipe, piece->module, DT_DEVICE_CPU, &iroi_full, &oroi_full,
               "tile (%zu,%zu)", tx, ty);

      /* prepare input tile buffer */
      input = dt_alloc_aligned((size_t)iroi_full.width * iroi_full.height * in_bpp);
      if(input == NULL)
      {
        dt_print(DT_DEBUG_TILING,
                 "[default_process_tiling_roi] [%s] could not alloc input buffer for module '%s%s'",
                 dt_dev_pixelpipe_type_to_str(piece->pipe->type), self->op, dt_iop_get_instance_id(self));
        goto error;
      }
      output = dt_alloc_aligned((size_t)oroi_full.width * oroi_full.height * out_bpp);
      if(output == NULL)
      {
        dt_print(DT_DEBUG_TILING,
                 "[default_process_tiling_roi] [%s] could not alloc output buffer for module '%s%s'",
                 dt_dev_pixelpipe_type_to_str(piece->pipe->type), self->op, dt_iop_get_instance_id(self));
        goto error;
      }

      DT_OMP_FOR(shared(iroi_full))
      for(size_t j = 0; j < iroi_full.height; j++)
        memcpy((char *)input + j * iroi_full.width * in_bpp, (char *)ivoid + ioffs + j * ipitch,
               (size_t)iroi_full.width * in_bpp);

      /* take original processed_maximum as starting point */
      for(int k = 0; k < 4; k++) piece->pipe->dsc.processed_maximum[k] = processed_maximum_saved[k];
      dt_dev_prepare_piece_cfa(piece, &iroi_full);

      /* call process() of module */
      self->process(self, piece, input, output, &iroi_full, &oroi_full);

      /* aggregate resulting processed_maximum */
      /* TODO: check if there really can be differences between tiles and take
               appropriate action (calculate minimum, maximum, average, ...?) */
      for(int k = 0; k < 4; k++)
      {
        if(tx + ty > 0 && fabs(processed_maximum_new[k] - piece->pipe->dsc.processed_maximum[k]) > 1.0e-6f)
          dt_print(DT_DEBUG_TILING,
                   "[default_process_tiling_roi] processed_maximum[%d] differs between tiles in module '%s%s'",
                   k, self->op, dt_iop_get_instance_id(self));
        processed_maximum_new[k] = piece->pipe->dsc.processed_maximum[k];
      }

      /* copy "good" part of tile to output buffer */
      const int origin_x = oroi_good.x - oroi_full.x;
      const int origin_y = oroi_good.y - oroi_full.y;
      DT_OMP_FOR(shared(oroi_good, oroi_full))
      for(size_t j = 0; j < oroi_good.height; j++)
        memcpy((char *)ovoid + ooffs + j * opitch,
               (char *)output + ((j + origin_y) * oroi_full.width + origin_x) * out_bpp,
               (size_t)oroi_good.width * out_bpp);

      dt_free_align(input);
      dt_free_align(output);
      input = output = NULL;
    }

  /* copy back final processed_maximum */
  for(int k = 0; k < 4; k++) piece->pipe->dsc.processed_maximum[k] = processed_maximum_new[k];

  dt_free_align(input);
  dt_free_align(output);
  piece->pipe->tiling = FALSE;
  return;

error:
  dt_control_log(_("tiling failed for module '%s%s'. the output most likely will be OK, but you might want to check."),
    self->op,dt_iop_get_instance_id(self));
// fall through

fallback:
  dt_free_align(input);
  dt_free_align(output);
  piece->pipe->tiling = FALSE;
  dt_print(DT_DEBUG_TILING,
           "[default_process_tiling_roi] [%s] fall back to standard processing for module '%s%s'",
           dt_dev_pixelpipe_type_to_str(piece->pipe->type), self->op, dt_iop_get_instance_id(self));
  self->process(self, piece, ivoid, ovoid, roi_in, roi_out);
  return;
}



/* if a module does not implement process_tiling() by itself, this function is called instead.
   _default_process_tiling_ptp() is able to handle standard cases where pixels do not change their places.
   _default_process_tiling_roi() takes care of all other cases where image gets distorted and for module
   "clipping",
   "flip" as this may flip or mirror the image. */
void default_process_tiling(dt_iop_module_t *self,
                            dt_dev_pixelpipe_iop_t *piece,
                            const void *const ivoid,
                            void *const ovoid,
                            const dt_iop_roi_t *const roi_in,
                            const dt_iop_roi_t *const roi_out,
                            const int in_bpp)
{
  const gboolean use_roi = memcmp(roi_in, roi_out, sizeof(struct dt_iop_roi_t)) || (self->flags() & IOP_FLAGS_TILING_FULL_ROI);
  if(use_roi)
    _default_process_tiling_roi(self, piece, ivoid, ovoid, roi_in, roi_out, in_bpp);
  else
    _default_process_tiling_ptp(self, piece, ivoid, ovoid, roi_in, roi_out, in_bpp);
  return;
}

#ifdef HAVE_OPENCL
/* simple tiling algorithm for roi_in == roi_out, i.e. for pixel to pixel modules/operations */
static int _default_process_tiling_cl_ptp(dt_iop_module_t *self,
                                          dt_dev_pixelpipe_iop_t *piece,
                                          const void *const ivoid,
                                          void *const ovoid,
                                          const dt_iop_roi_t *const roi_in,
                                          const dt_iop_roi_t *const roi_out,
                                          const int in_bpp)
{
  cl_int err = CL_SUCCESS; // we take care for errors below
  cl_mem input = NULL;
  cl_mem output = NULL;
  size_t cltile_w = 0;
  size_t cltile_h = 0;

  cl_mem pinned_input = NULL;
  cl_mem pinned_output = NULL;
  void *input_buffer = NULL;
  void *output_buffer = NULL;

  dt_iop_buffer_dsc_t dsc;
  self->output_format(self, piece->pipe, piece, &dsc);
  const int out_bpp = dt_iop_buffer_dsc_to_bpp(&dsc);

  const int devid = piece->pipe->devid;
  const int ipitch = roi_in->width * in_bpp;
  const int opitch = roi_out->width * out_bpp;
  const int max_bpp = MAX(in_bpp, out_bpp);

  /* get tiling requirements of module */
  dt_develop_tiling_t tiling = { 0 };
  tiling.factor_cl = tiling.maxbuf_cl = -1.0f;
  self->tiling_callback(self, piece, roi_in, roi_out, &tiling);
  if(tiling.factor_cl < 0.0f) tiling.factor_cl = tiling.factor;
  if(tiling.maxbuf_cl < 0.0f) tiling.maxbuf_cl = tiling.maxbuf;

  /* shall we use pinned memory transfers? */
  gboolean use_pinned_memory = dt_opencl_use_pinned_memory(devid);
  /* If using pinned transfer on devices with dedicated GPU mem there is an additional
     mem pressure as they will allocate also on device as cache for performance
  */
  const float pinned_buffer_overhead = use_pinned_memory && !dt_opencl_unified_memory(devid) ? 2.0f : 0.0f;

  // avoid problems when pinned buffer size gets too close to max_mem_alloc size
  const float pinned_buffer_slack = use_pinned_memory ? 0.85f : 1.0f;
  const float available = (float)dt_opencl_get_device_available(devid);
  const float factor = fmaxf(tiling.factor_cl + pinned_buffer_overhead, 1.0f);
  const float singlebuffer = fminf(fmaxf((available - tiling.overhead) / factor, 0.0f),
                                  pinned_buffer_slack * (float)(dt_opencl_get_device_memalloc(devid)));
  const float maxbuf = fmaxf(tiling.maxbuf_cl, 1.0f);
  int width = MIN(roi_in->width, darktable.opencl->dev[devid].max_image_width);
  int height = MIN(roi_in->height, darktable.opencl->dev[devid].max_image_height);

  /* shrink tile size in case it would exceed singlebuffer size */
  if((float)width * height * max_bpp * maxbuf > singlebuffer)
  {
    const float scale = singlebuffer / ((float)width * height * max_bpp * maxbuf);

    if(width < height && scale >= 0.333f)
    {
      height = floorf(height * scale);
    }
    else if(height <= width && scale >= 0.333f)
    {
      width = floorf(width * scale);
    }
    else
    {
      width = floorf(width * sqrtf(scale));
      height = floorf(height * sqrtf(scale));
    }
    dt_print(DT_DEBUG_TILING | DT_DEBUG_VERBOSE,
             "[default_process_tiling_cl_ptp] [%s] buffer exceeds singlebuffer, corrected to %dx%d",
             dt_dev_pixelpipe_type_to_str(piece->pipe->type), width, height);
  }

  /* make sure we have a reasonably effective tile dimension. if not try square tiles */
  if(3 * tiling.overlap > width || 3 * tiling.overlap > height)
  {
    width = height = floorf(sqrtf((float)width * height));
    dt_print(DT_DEBUG_TILING | DT_DEBUG_VERBOSE,
             "[default_process_tiling_cl_ptp] [%s] use squares because of overlap, corrected to %dx%d",
             dt_dev_pixelpipe_type_to_str(piece->pipe->type), width, height);
  }

  /* Alignment rules: we need to make sure that alignment requirements of module are fulfilled.
     Modules will report alignment requirements via align within tiling_callback().
     Additional alignment requirements are set via dt_opencl_tiling_align().
     We guarantee alignment by selecting image width/height and overlap accordingly. For a tile width/height
     that is identical to image width/height no special alignment is done.
  */
  const unsigned int align = tiling.align;

  const unsigned int walign = _lcm(align, dt_opencl_tiling_align(devid));
  const unsigned int halign = align;
  assert(align != 0 && walign != 0 && halign != 0);

  /* properly align tile width and height by making them smaller if needed */
  if(width < roi_in->width) width = (width / walign) * walign;
  if(height < roi_in->height) height = (height / halign) * halign;

  /* also make sure that overlap follows alignment rules by making it wider when needed */
  const int overlap = tiling.overlap % align != 0 ? (tiling.overlap / align + 1) * align
                                                    : tiling.overlap;

  /* calculate effective tile size */
  const int tile_wd = width - 2 * overlap > 0 ? width - 2 * overlap : 1;
  const int tile_ht = height - 2 * overlap > 0 ? height - 2 * overlap : 1;

  /* calculate number of tiles */
  const int tiles_x = width < roi_in->width ? ceilf(roi_in->width / (float)tile_wd) : 1;
  const int tiles_y = height < roi_in->height ? ceilf(roi_in->height / (float)tile_ht) : 1;

  /* sanity check: don't run wild on too many tiles */
  if(tiles_x * tiles_y > _maximum_number_tiles())
  {
    dt_print(DT_DEBUG_TILING,
             "[default_process_tiling_cl_ptp] [%s] aborted tiling for module '%s%s'. "
             "too many tiles: %d x %d",
             dt_dev_pixelpipe_type_to_str(piece->pipe->type),
             self->op, dt_iop_get_instance_id(self), tiles_x, tiles_y);
    return DT_OPENCL_PROCESS_CL;
  }

  /* store processed_maximum to be re-used and aggregated */
  dt_aligned_pixel_t processed_maximum_saved;
  dt_aligned_pixel_t processed_maximum_new = { 1.0f };
  for_four_channels(k) processed_maximum_saved[k] = piece->pipe->dsc.processed_maximum[k];

  /* reserve pinned input and output memory for host<->device data transfer */
  if(use_pinned_memory)
  {
    const size_t bsize = (size_t)width * height * in_bpp;
    pinned_input = dt_opencl_alloc_device_buffer_with_flags(devid, bsize, CL_MEM_READ_ONLY | CL_MEM_ALLOC_HOST_PTR);

    if(pinned_input)
      input_buffer = dt_opencl_map_buffer(devid, pinned_input, TRUE, CL_MAP_WRITE, 0, bsize);
    if(input_buffer == NULL)
    {
      dt_print(DT_DEBUG_OPENCL | DT_DEBUG_TILING,
               "[default_process_tiling_cl_ptp] [%s] could not map pinned input buffer to host "
               "memory for module '%s%s'",
               dt_dev_pixelpipe_type_to_str(piece->pipe->type), self->op, dt_iop_get_instance_id(self));
      use_pinned_memory = FALSE;
    }
  }

  if(use_pinned_memory)
  {
    const size_t bsize = (size_t)width * height * out_bpp;
    pinned_output = dt_opencl_alloc_device_buffer_with_flags(devid, bsize, CL_MEM_WRITE_ONLY | CL_MEM_ALLOC_HOST_PTR);
    if(pinned_output)
      output_buffer = dt_opencl_map_buffer(devid, pinned_output, TRUE, CL_MAP_READ, 0, bsize);

    if(output_buffer == NULL)
    {
      dt_print(DT_DEBUG_OPENCL | DT_DEBUG_TILING,
               "[default_process_tiling_cl_ptp] [%s] could not map pinned output buffer to host "
               "memory for module '%s%s'",
               dt_dev_pixelpipe_type_to_str(piece->pipe->type), self->op, dt_iop_get_instance_id(self));
      use_pinned_memory = FALSE;
    }
  }

  dt_print_pipe(DT_DEBUG_PIPE | DT_DEBUG_TILING,
                        "default *tiled* cl_ptp", piece->pipe, piece->module, devid, roi_in, roi_out,
                        "%dx%d tiles%s, size=%dx%d, overlap=%d",
                        tiles_x, tiles_y, (use_pinned_memory) ? ", pinned" : "", tile_wd, tile_ht, overlap);
  /* iterate over tiles */
  piece->pipe->tiling = TRUE;
  for(size_t tx = 0; tx < tiles_x; tx++)
  {
    for(size_t ty = 0; ty < tiles_y; ty++)
    {
      const size_t wd = tx * tile_wd + width > roi_in->width ? roi_in->width - tx * tile_wd : width;
      const size_t ht = ty * tile_ht + height > roi_in->height ? roi_in->height - ty * tile_ht : height;

      /* no need to process (end)tiles that are smaller than the total overlap area */
      const gboolean skipped = (wd <= 2 * overlap && tx > 0) || (ht <= 2 * overlap && ty > 0);

      /* origin and region of effective part of tile, which we want to store later */
      size_t origin[2] = { 0, 0 };
      size_t region[2] = { wd, ht };

      /* roi_in and roi_out for process_cl on subbuffer */
      const dt_iop_roi_t iroi = { roi_in->x + tx * tile_wd, roi_in->y + ty * tile_ht, wd, ht, roi_in->scale };
      const dt_iop_roi_t oroi = { roi_out->x + tx * tile_wd, roi_out->y + ty * tile_ht, wd, ht, roi_out->scale };


      /* offsets of tile into ivoid and ovoid */
      const size_t ioffs = (ty * tile_ht) * ipitch + (tx * tile_wd) * in_bpp;
      size_t ooffs = (ty * tile_ht) * opitch + (tx * tile_wd) * out_bpp;

      dt_print_pipe(DT_DEBUG_TILING,
               skipped ? "  tile cl_ptp skipped" : "  tile cl_ptp", piece->pipe, piece->module, devid, &iroi, &oroi,
               "tile (%zu,%zu)", tx, ty);
      if(skipped) continue;

      /* get input and output buffers */
      if(cltile_w != wd || cltile_h != ht)
      {
        dt_opencl_release_mem_object(input);
        dt_opencl_release_mem_object(output);
        input = dt_opencl_alloc_device(devid, wd, ht, in_bpp);
        output = dt_opencl_alloc_device(devid, wd, ht, out_bpp);
        if(output == NULL || input == NULL)
        {
          err = CL_MEM_OBJECT_ALLOCATION_FAILURE;
          goto error;
        }
        cltile_w = wd;
        cltile_h = ht;
      }

      if(use_pinned_memory)
      {
        /* prepare pinned input tile buffer: copy part of input image */
        DT_OMP_FOR()
        for(size_t j = 0; j < ht; j++)
          memcpy((char *)input_buffer + j * wd * in_bpp, (char *)ivoid + ioffs + j * ipitch,
                 (size_t)wd * in_bpp);

        /* blocking memory transfer: pinned host input buffer -> opencl/device tile */
        err = dt_opencl_write_host_to_image_raw(devid, (char *)input_buffer, input, origin, region,
                                                 wd * in_bpp, TRUE);
        if(err != CL_SUCCESS) use_pinned_memory = FALSE;
      }
      else
      {
        /* blocking direct memory transfer: host input image -> opencl/device tile */
        err = dt_opencl_write_host_to_image_raw(devid, (char *)ivoid + ioffs, input, origin, region, ipitch, TRUE);
      }
      if(err != CL_SUCCESS) goto error;

      /* take original processed_maximum as starting point */
      for(int k = 0; k < 4; k++) piece->pipe->dsc.processed_maximum[k] = processed_maximum_saved[k];
      dt_dev_prepare_piece_cfa(piece, &iroi);
      /* call process_cl of module */
      err = self->process_cl(self, piece, input, output, &iroi, &oroi);
      if(err != CL_SUCCESS)
        goto error;

      /* aggregate resulting processed_maximum */
      /* TODO: check if there really can be differences between tiles and take
               appropriate action (calculate minimum, maximum, average, ...?) */
      for(int k = 0; k < 4; k++)
      {
        if(tx + ty > 0 && fabs(processed_maximum_new[k] - piece->pipe->dsc.processed_maximum[k]) > 1.0e-6f)
          dt_print(DT_DEBUG_TILING,
                   "[default_process_tiling_cl_ptp] [%s] processed_maximum[%d] differs between tiles in module '%s%s'",
                   dt_dev_pixelpipe_type_to_str(piece->pipe->type), k, self->op, dt_iop_get_instance_id(self));
        processed_maximum_new[k] = piece->pipe->dsc.processed_maximum[k];
      }

      if(use_pinned_memory)
      {
        /* blocking memory transfer: complete opencl/device tile -> pinned host output buffer */
        err = dt_opencl_read_host_from_image_raw(devid, (char *)output_buffer, output, origin, region,
                                                  wd * out_bpp, TRUE);
        if(err != CL_SUCCESS)
        {
          use_pinned_memory = FALSE;
          goto error;
        }
      }

      /* correct origin and region of tile for overlap.
         makes sure that we only copy back the "good" part. */
      if(tx > 0)
      {
        origin[0] += overlap;
        region[0] -= overlap;
        ooffs += (size_t)overlap * out_bpp;
      }
      if(ty > 0)
      {
        origin[1] += overlap;
        region[1] -= overlap;
        ooffs += (size_t)overlap * opitch;
      }

      if(use_pinned_memory)
      {
        /* copy "good" part of tile from pinned output buffer to output image */
        //        DT_OMP_FOR(shared(origin, region))
        for(size_t j = 0; j < region[1]; j++)
          memcpy((char *)ovoid + ooffs + j * opitch,
                 (char *)output_buffer + ((j + origin[1]) * wd + origin[0]) * out_bpp,
                 (size_t)region[0] * out_bpp);
      }
      else
      {
        /* blocking direct memory transfer: good part of opencl/device tile -> host output image */
        err = dt_opencl_read_host_from_image_raw(devid, (char *)ovoid + ooffs, output, origin, region,
                                                  opitch, TRUE);
        if(err != CL_SUCCESS) goto error;
      }

      /* block until opencl queue has finished to free all used event handlers */
      dt_opencl_finish_sync_pipe(devid, piece->pipe->type);
    }
  }

  /* copy back final processed_maximum */
  for(int k = 0; k < 4; k++) piece->pipe->dsc.processed_maximum[k] = processed_maximum_new[k];

  if(input_buffer != NULL) dt_opencl_unmap_mem_object(devid, pinned_input, input_buffer);
  dt_opencl_release_mem_object(pinned_input);
  if(output_buffer != NULL) dt_opencl_unmap_mem_object(devid, pinned_output, output_buffer);
  dt_opencl_release_mem_object(pinned_output);
  dt_opencl_release_mem_object(input);
  dt_opencl_release_mem_object(output);
  piece->pipe->tiling = FALSE;
  return CL_SUCCESS;

error:
  /* copy back stored processed_maximum */
  for(int k = 0; k < 4; k++) piece->pipe->dsc.processed_maximum[k] = processed_maximum_saved[k];
  if(input_buffer != NULL) dt_opencl_unmap_mem_object(devid, pinned_input, input_buffer);
  dt_opencl_release_mem_object(pinned_input);
  if(output_buffer != NULL) dt_opencl_unmap_mem_object(devid, pinned_output, output_buffer);
  dt_opencl_release_mem_object(pinned_output);
  dt_opencl_release_mem_object(input);
  dt_opencl_release_mem_object(output);
  piece->pipe->tiling = FALSE;
  const gboolean pinning_error = !use_pinned_memory && dt_opencl_use_pinned_memory(devid);
  dt_print(DT_DEBUG_TILING | DT_DEBUG_OPENCL,
           "[default_process_tiling_opencl_ptp] [%s] couldn't run process_cl() for "
           "module '%s%s' in tiling mode:%s %s",
           dt_dev_pixelpipe_type_to_str(piece->pipe->type), self->op, dt_iop_get_instance_id(self),
           (pinning_error) ? " pinning problem" : "", cl_errstr(err));

  if(pinning_error) darktable.opencl->dev[devid].pinned_error = TRUE;
  return err;
}

/* more elaborate tiling algorithm for roi_in != roi_out: slower than the ptp variant,
   more tiles and larger overlap */
static int _default_process_tiling_cl_roi(dt_iop_module_t *self,
                                          dt_dev_pixelpipe_iop_t *piece,
                                          const void *const ivoid,
                                          void *const ovoid,
                                          const dt_iop_roi_t *const roi_in,
                                          const dt_iop_roi_t *const roi_out,
                                          const int in_bpp)
{
  cl_int err = CL_SUCCESS; // we take care for errors below
  cl_mem input = NULL;
  cl_mem output = NULL;
  size_t cltile_ow = 0;
  size_t cltile_oh = 0;
  size_t cltile_iw = 0;
  size_t cltile_ih = 0;
  cl_mem pinned_input = NULL;
  cl_mem pinned_output = NULL;
  void *input_buffer = NULL;
  void *output_buffer = NULL;

  dt_iop_buffer_dsc_t dsc;
  self->output_format(self, piece->pipe, piece, &dsc);
  const int out_bpp = dt_iop_buffer_dsc_to_bpp(&dsc);

  const int devid = piece->pipe->devid;
  const int ipitch = roi_in->width * in_bpp;
  const int opitch = roi_out->width * out_bpp;
  const int max_bpp = MAX(in_bpp, out_bpp);

  const float fullscale = fmaxf(roi_in->scale / roi_out->scale, sqrtf(((float)roi_in->width * roi_in->height)
                                                              / ((float)roi_out->width * roi_out->height)));

  /* inaccuracy for roi_in elements in roi_out -> roi_in calculations */
  const int delta = ceilf(fullscale);

  /* estimate for additional (space) requirement in buffer dimensions due to inaccuracies */
  const int inacc = RESERVE * delta;

  /* get tiling requirements of module */
  dt_develop_tiling_t tiling = { 0 };
  tiling.factor_cl = tiling.maxbuf_cl = -1;
  self->tiling_callback(self, piece, roi_in, roi_out, &tiling);
  if(tiling.factor_cl < 0) tiling.factor_cl = tiling.factor;
  if(tiling.maxbuf_cl < 0) tiling.maxbuf_cl = tiling.maxbuf;

  /* shall we use pinned memory transfers? */
  gboolean use_pinned_memory = dt_opencl_use_pinned_memory(devid);

  /* If using pinned transfer on devices with dedicated GPU mem there is an additional
     mem pressure as they will allocate also on device as cache for performance
  */
  const float pinned_buffer_overhead = use_pinned_memory && !dt_opencl_unified_memory(devid) ? 2.0f : 0.0f;
  // avoid problems when pinned buffer size gets too close to max_mem_alloc size
  const float pinned_buffer_slack = use_pinned_memory ? 0.85f : 1.0f;
  const float available = (float)dt_opencl_get_device_available(devid);
  const float factor = fmaxf(tiling.factor_cl + pinned_buffer_overhead, 1.0f);
  const float singlebuffer = fminf(fmaxf((available - tiling.overhead) / factor, 0.0f),
                                  pinned_buffer_slack * (float)(dt_opencl_get_device_memalloc(devid)));
  const float maxbuf = fmaxf(tiling.maxbuf_cl, 1.0f);

  int width = MIN(MAX(roi_in->width, roi_out->width), darktable.opencl->dev[devid].max_image_width);
  int height = MIN(MAX(roi_in->height, roi_out->height), darktable.opencl->dev[devid].max_image_height);

  const unsigned int align = _lcm(tiling.align, dt_opencl_tiling_align(devid));
  assert(align != 0);

  /* shrink tile size in case it would exceed singlebuffer size */
  if((float)width * height * max_bpp * maxbuf > singlebuffer)
  {
    const float scale = singlebuffer / ((float)width * height * max_bpp * maxbuf);

    if(width < height && scale >= 0.333f)
    {
      height = _align_down((int)floorf(height * scale), align);
    }
    else if(height <= width && scale >= 0.333f)
    {
      width = _align_down((int)floorf(width * scale), align);
    }
    else
    {
      width = _align_down((int)floorf(width * sqrtf(scale)), align);
      height = _align_down((int)floorf(height * sqrtf(scale)), align);
    }
    dt_print(DT_DEBUG_TILING | DT_DEBUG_OPENCL | DT_DEBUG_VERBOSE,
             "[default_process_tiling_cl_roi] [%s] buffer exceeds singlebuffer, corrected to %dx%d",
             dt_dev_pixelpipe_type_to_str(piece->pipe->type), width, height);
  }

  /* make sure we have a reasonably effective tile dimension. if not try square tiles */
  if(3 * tiling.overlap > width || 3 * tiling.overlap > height)
  {
    width = height = _align_down((int)floorf(sqrtf((float)width * height)), align);
    dt_print(DT_DEBUG_TILING | DT_DEBUG_VERBOSE,
             "[default_process_tiling_cl_roi] [%s] use squares because of overlap, corrected to %dx%d",
             dt_dev_pixelpipe_type_to_str(piece->pipe->type), width, height);
  }

  /* make sure that overlap follows alignment rules by making it wider when needed.
     overlap_in needs to be aligned, overlap_out is only here to calculate output buffer size */
  const int overlap_in = _align_up(tiling.overlap, align);
  const int overlap_out = ceilf((float)overlap_in / fullscale);

  int tiles_x = 1, tiles_y = 1;

  /* calculate number of tiles taking the larger buffer (input or output) as a guiding one.
     normally it is roi_in > roi_out; but let's be prepared */
  if(roi_in->width > roi_out->width)
    tiles_x = width < roi_in->width
                  ? ceilf((float)roi_in->width / (float)MAX(width - 2 * overlap_in - inacc, 1))
                  : 1;
  else
    tiles_x = width < roi_out->width ? ceilf((float)roi_out->width / (float)MAX(width - 2 * overlap_out, 1))
                                     : 1;

  if(roi_in->height > roi_out->height)
    tiles_y = height < roi_in->height
                  ? ceilf((float)roi_in->height / (float)MAX(height - 2 * overlap_in - inacc, 1))
                  : 1;
  else
    tiles_y = height < roi_out->height
                  ? ceilf((float)roi_out->height / (float)MAX(height - 2 * overlap_out, 1))
                  : 1;

  /* sanity check: don't run wild on too many tiles */
  if(tiles_x * tiles_y > _maximum_number_tiles())
  {
    dt_print(DT_DEBUG_TILING,
             "[default_process_tiling_cl_roi] [%s] aborted tiling for module '%s%s'. too many tiles: %dx%d",
             dt_dev_pixelpipe_type_to_str(piece->pipe->type),
             self->op, dt_iop_get_instance_id(self), tiles_x, tiles_y);
    return DT_OPENCL_PROCESS_CL;
  }

  /* calculate tile width and height excl. overlap (i.e. the good part) for output.
     important for all following processing steps. */
  const int tile_wd = _align_up(
      roi_out->width % tiles_x == 0 ? roi_out->width / tiles_x : roi_out->width / tiles_x + 1, align);
  const int tile_ht = _align_up(
      roi_out->height % tiles_y == 0 ? roi_out->height / tiles_y : roi_out->height / tiles_y + 1, align);

  dt_print_pipe(DT_DEBUG_PIPE | DT_DEBUG_TILING,
                        "process *tiled* roi", piece->pipe, piece->module, devid, roi_in, roi_out,
                        "%dx%d tiles%s, size=%dx%d",
                        tiles_x, tiles_y, (use_pinned_memory) ? ", pinned" : "", tile_wd, tile_ht);
  /* store processed_maximum to be re-used and aggregated */
  dt_aligned_pixel_t processed_maximum_saved;
  dt_aligned_pixel_t processed_maximum_new = { 1.0f, 1.0f, 1.0f, 1.0f };
  for_four_channels(k) processed_maximum_saved[k] = piece->pipe->dsc.processed_maximum[k];

  /* reserve pinned input and output memory for host<->device data transfer */
  if(use_pinned_memory)
  {
    const size_t bsize = (size_t)width * height * in_bpp;
    pinned_input = dt_opencl_alloc_device_buffer_with_flags(devid, bsize, CL_MEM_READ_ONLY | CL_MEM_ALLOC_HOST_PTR);
    if(pinned_input)
      input_buffer = dt_opencl_map_buffer(devid, pinned_input, TRUE, CL_MAP_WRITE, 0, bsize);
    if(input_buffer == NULL)
    {
      dt_print(DT_DEBUG_OPENCL | DT_DEBUG_TILING,
               "[default_process_tiling_cl_roi] [%s] could not map pinned input buffer to host "
               "memory for module '%s%s'",
               dt_dev_pixelpipe_type_to_str(piece->pipe->type), self->op, dt_iop_get_instance_id(self));
      use_pinned_memory = FALSE;
    }
  }

  if(use_pinned_memory)
  {
    const size_t bsize = (size_t)width * height * out_bpp;
    pinned_output = dt_opencl_alloc_device_buffer_with_flags(devid, bsize, CL_MEM_WRITE_ONLY | CL_MEM_ALLOC_HOST_PTR);
    if(pinned_output)
      output_buffer = dt_opencl_map_buffer(devid, pinned_output, TRUE, CL_MAP_READ, 0, bsize);
    if(output_buffer == NULL)
    {
      dt_print(DT_DEBUG_OPENCL | DT_DEBUG_TILING,
               "[default_process_tiling_cl_roi] [%s] could not map pinned output buffer to host "
               "memory for module '%s%s'",
               dt_dev_pixelpipe_type_to_str(piece->pipe->type), self->op, dt_iop_get_instance_id(self));
      use_pinned_memory = FALSE;
    }
  }


  /* iterate over tiles */
  piece->pipe->tiling = TRUE;
  for(size_t tx = 0; tx < tiles_x; tx++)
  {
    for(size_t ty = 0; ty < tiles_y; ty++)
    {
      /* the output dimensions of the good part of this specific tile */
      const size_t wd = (tx + 1) * tile_wd > roi_out->width ? (size_t)roi_out->width - tx * tile_wd : tile_wd;
      const size_t ht = (ty + 1) * tile_ht > roi_out->height ? (size_t)roi_out->height - ty * tile_ht : tile_ht;

      /* roi_in and roi_out of good part: oroi_good easy to calculate based on number and dimension of tile.
         iroi_good is calculated by modify_roi_in() of respective module */
      dt_iop_roi_t iroi_good = { roi_in->x  + tx * tile_wd, roi_in->y  + ty * tile_ht, wd, ht, roi_in->scale };
      dt_iop_roi_t oroi_good = { roi_out->x + tx * tile_wd, roi_out->y + ty * tile_ht, wd, ht, roi_out->scale };

      self->modify_roi_in(self, piece, &oroi_good, &iroi_good);

      /* clamp iroi_good to not exceed roi_in */
      iroi_good.x = MAX(iroi_good.x, roi_in->x);
      iroi_good.y = MAX(iroi_good.y, roi_in->y);
      iroi_good.width = MIN(iroi_good.width, roi_in->width + roi_in->x - iroi_good.x);
      iroi_good.height = MIN(iroi_good.height, roi_in->height + roi_in->y - iroi_good.y);

      _print_roi(&iroi_good, "tile iroi_good");
      _print_roi(&oroi_good, "tile oroi_good");

      /* now we need to calculate full region of this tile: increase input roi to take care of overlap
         requirements and alignment and add additional delta to correct for possible rounding errors in
         modify_roi_in()
         -> generates first estimate of iroi_full
      */
      const int x_in = iroi_good.x;
      const int y_in = iroi_good.y;
      const int width_in = iroi_good.width;
      const int height_in = iroi_good.height;
      const int new_x_in = MAX(_align_close(x_in - overlap_in - delta, align), roi_in->x);
      const int new_y_in = MAX(_align_close(y_in - overlap_in - delta, align), roi_in->y);
      const int new_width_in = MIN(_align_up(width_in + overlap_in + delta + (x_in - new_x_in), align),
                                    roi_in->width + roi_in->x - new_x_in);
      const int new_height_in = MIN(_align_up(height_in + overlap_in + delta + (y_in - new_y_in), align),
                                     roi_in->height + roi_in->y - new_y_in);

      /* iroi_full based on calculated numbers and dimensions. oroi_full just set as a starting point for the
       * following iterative search */
      dt_iop_roi_t iroi_full = { new_x_in, new_y_in, new_width_in, new_height_in, iroi_good.scale };
      dt_iop_roi_t oroi_full = oroi_good; // a good starting point for optimization

      _print_roi(&iroi_full, "tile iroi_full before optimization");
      _print_roi(&oroi_full, "tile oroi_full before optimization");

      /* try to find a matching oroi_full */
      if(!_fit_output_to_input_roi(self, piece, &iroi_full, &oroi_full, delta, 10))
      {
        dt_print(DT_DEBUG_OPENCL | DT_DEBUG_TILING,
                 "[default_process_tiling_cl_roi] [%s] can not handle requested roi's tiling "
                 "for module '%s%s' not possible",
                 dt_dev_pixelpipe_type_to_str(piece->pipe->type), self->op, dt_iop_get_instance_id(self));
        err = DT_OPENCL_PROCESS_CL;
        goto error;
      }


      /* make sure that oroi_full at least covers the range of oroi_good.
         this step is needed due to the possibility of rounding errors */
      oroi_full.x = MIN(oroi_full.x, oroi_good.x);
      oroi_full.y = MIN(oroi_full.y, oroi_good.y);
      oroi_full.width = MAX(oroi_full.width, oroi_good.x + oroi_good.width - oroi_full.x);
      oroi_full.height = MAX(oroi_full.height, oroi_good.y + oroi_good.height - oroi_full.y);

      /* clamp oroi_full to not exceed roi_out */
      oroi_full.x = MAX(oroi_full.x, roi_out->x);
      oroi_full.y = MAX(oroi_full.y, roi_out->y);
      oroi_full.width = MIN(oroi_full.width, roi_out->width + roi_out->x - oroi_full.x);
      oroi_full.height = MIN(oroi_full.height, roi_out->height + roi_out->y - oroi_full.y);


      /* calculate final iroi_full */
      self->modify_roi_in(self, piece, &oroi_full, &iroi_full);

      /* clamp iroi_full to not exceed roi_in */
      iroi_full.x = MAX(iroi_full.x, roi_in->x);
      iroi_full.y = MAX(iroi_full.y, roi_in->y);
      iroi_full.width = MIN(iroi_full.width, roi_in->width + roi_in->x - iroi_full.x);
      iroi_full.height = MIN(iroi_full.height, roi_in->height + roi_in->y - iroi_full.y);

      _print_roi(&iroi_full, "tile iroi_full");
      _print_roi(&oroi_full, "tile oroi_full");

      /* offsets of tile into ivoid and ovoid */
      const int in_dx = iroi_full.x - roi_in->x;
      const int in_dy = iroi_full.y - roi_in->y;
      const int out_dx = oroi_good.x - roi_out->x;
      const int out_dy = oroi_good.y - roi_out->y;
      const size_t ioffs = (size_t)(in_dy  * ipitch) + (size_t)(in_dx * in_bpp);
      const size_t ooffs = (size_t)(out_dy * opitch) + (size_t)(out_dx * out_bpp);

      /* region of full input tile */
      const size_t iregion[2] = { iroi_full.width, iroi_full.height };

      /* region of full output tile */
      const size_t ofregion[2] = { oroi_full.width, oroi_full.height };

      /* origin and region of good part of output tile */
      const size_t oorigin[2] = { oroi_good.x - oroi_full.x, oroi_good.y - oroi_full.y };
      const size_t oregion[2] = { oroi_good.width, oroi_good.height };

      dt_print_pipe(DT_DEBUG_TILING,
              "  tile cl_roi", piece->pipe, piece->module, devid, &iroi_full, &oroi_full,
              "tile (%zu,%zu)", tx, ty);
      dt_print_pipe(DT_DEBUG_TILING | DT_DEBUG_VERBOSE,
             "  tile cl_roi", piece->pipe, piece->module, devid, &iroi_full, &oroi_full,
              "tile (%zu,%zu)  offsets=[%i,%i] delta=%i", tx, ty, out_dx, out_dy, delta);

      /* get opencl input and output buffers */
      if(cltile_iw != iroi_full.width || cltile_ih != iroi_full.height)
      {
        dt_opencl_release_mem_object(input);
        input = dt_opencl_alloc_device(devid, iroi_full.width, iroi_full.height, in_bpp);
        cltile_iw = iroi_full.width;
        cltile_ih = iroi_full.height;
      }

      if(cltile_ow != oroi_full.width || cltile_oh != oroi_full.height)
      {
        dt_opencl_release_mem_object(output);
        output = dt_opencl_alloc_device(devid, oroi_full.width, oroi_full.height, out_bpp);
        cltile_ow = oroi_full.width;
        cltile_oh = oroi_full.height;
      }

      if(output == NULL || input == NULL)
      {
        err = CL_MEM_OBJECT_ALLOCATION_FAILURE;
        goto error;
      }

      if(use_pinned_memory)
      {
        /* prepare pinned input tile buffer: copy part of input image */
        DT_OMP_FOR(shared(iroi_full))
        for(size_t j = 0; j < iroi_full.height; j++)
          memcpy((char *)input_buffer + j * iroi_full.width * in_bpp, (char *)ivoid + ioffs + j * ipitch,
                 (size_t)iroi_full.width * in_bpp);

        /* blocking memory transfer: pinned host input buffer -> opencl/device tile */
        err = dt_opencl_write_host_to_image_raw(devid, (char *)input_buffer, input, CLIMG_ORIGIN, iregion,
                                                 (size_t)iroi_full.width * in_bpp, TRUE);
        if(err != CL_SUCCESS) use_pinned_memory = FALSE;
      }
      else
      {
        /* blocking direct memory transfer: host input image -> opencl/device tile */
        err = dt_opencl_write_host_to_image_raw(devid, (char *)ivoid + ioffs, input, CLIMG_ORIGIN, iregion,
                                                 ipitch, TRUE);
      }
      if(err != CL_SUCCESS) goto error;

      /* take original processed_maximum as starting point */
      for(int k = 0; k < 4; k++) piece->pipe->dsc.processed_maximum[k] = processed_maximum_saved[k];

      dt_dev_prepare_piece_cfa(piece, &iroi_full);
      /* call process_cl of module */
      err = self->process_cl(self, piece, input, output, &iroi_full, &oroi_full);
      if(err != CL_SUCCESS)
        goto error;

      /* aggregate resulting processed_maximum */
      /* TODO: check if there really can be differences between tiles and take
               appropriate action (calculate minimum, maximum, average, ...?) */
      for(int k = 0; k < 4; k++)
      {
        if(tx + ty > 0 && fabs(processed_maximum_new[k] - piece->pipe->dsc.processed_maximum[k]) > 1.0e-6f)
          dt_print(DT_DEBUG_TILING,
                   "[default_process_tiling_cl_roi] [%s] processed_maximum[%d] "
                   "differs between tiles in module '%s%s'",
                   dt_dev_pixelpipe_type_to_str(piece->pipe->type), k, self->op, dt_iop_get_instance_id(self));
        processed_maximum_new[k] = piece->pipe->dsc.processed_maximum[k];
      }

      if(use_pinned_memory)
      {
        /* blocking memory transfer: complete opencl/device tile -> pinned host output buffer */
        err = dt_opencl_read_host_from_image_raw(devid, (char *)output_buffer, output, CLIMG_ORIGIN, ofregion,
                                                  (size_t)oroi_full.width * out_bpp, TRUE);
        if(err != CL_SUCCESS)
        {
          use_pinned_memory = FALSE;
          goto error;
        }
        /* copy "good" part of tile from pinned output buffer to output image */
        DT_OMP_FOR(shared(oroi_full, oorigin, oregion))
        for(size_t j = 0; j < oregion[1]; j++)
          memcpy((char *)ovoid + ooffs + j * opitch,
                 (char *)output_buffer + ((j + oorigin[1]) * oroi_full.width + oorigin[0]) * out_bpp,
                 (size_t)oregion[0] * out_bpp);
      }
      else
      {
        /* blocking direct memory transfer: good part of opencl/device tile -> host output image */
        err = dt_opencl_read_host_from_image_raw(devid, (char *)ovoid + ooffs, output, oorigin, oregion,
                                                  opitch, TRUE);
        if(err != CL_SUCCESS) goto error;
      }

      /* block until opencl queue has finished to free all used event handlers */
      dt_opencl_finish_sync_pipe(devid, piece->pipe->type);
    }
  }
  /* copy back final processed_maximum */
  for(int k = 0; k < 4; k++) piece->pipe->dsc.processed_maximum[k] = processed_maximum_new[k];
  if(input_buffer != NULL) dt_opencl_unmap_mem_object(devid, pinned_input, input_buffer);
  dt_opencl_release_mem_object(pinned_input);
  if(output_buffer != NULL) dt_opencl_unmap_mem_object(devid, pinned_output, output_buffer);
  dt_opencl_release_mem_object(pinned_output);
  dt_opencl_release_mem_object(input);
  dt_opencl_release_mem_object(output);
  piece->pipe->tiling = FALSE;
  return CL_SUCCESS;

error:
  /* copy back stored processed_maximum */
  for(int k = 0; k < 4; k++) piece->pipe->dsc.processed_maximum[k] = processed_maximum_saved[k];
  if(input_buffer != NULL) dt_opencl_unmap_mem_object(devid, pinned_input, input_buffer);
  dt_opencl_release_mem_object(pinned_input);
  if(output_buffer != NULL) dt_opencl_unmap_mem_object(devid, pinned_output, output_buffer);
  dt_opencl_release_mem_object(pinned_output);
  dt_opencl_release_mem_object(input);
  dt_opencl_release_mem_object(output);
  piece->pipe->tiling = FALSE;
  const gboolean pinning_error = (use_pinned_memory == FALSE) && dt_opencl_use_pinned_memory(devid);
  dt_print_pipe(DT_DEBUG_OPENCL | DT_DEBUG_TILING,
           "default tiling_cl_roi error", piece->pipe, piece->module, devid, roi_in, roi_out,
           "%serror=%s",
           (pinning_error) ? "pinning problem " : "", cl_errstr(err));

  if(pinning_error) darktable.opencl->dev[devid].pinned_error = TRUE;
  return err;
}


/* if a module does not implement process_tiling_cl() by itself, this function is called instead.
   _default_process_tiling_cl_ptp() is able to handle standard cases where pixels do not change their places.
   _default_process_tiling_cl_roi() takes care of all other cases where image gets distorted. */
int default_process_tiling_cl(dt_iop_module_t *self,
                              dt_dev_pixelpipe_iop_t *piece,
                              const void *const ivoid,
                              void *const ovoid,
                              const dt_iop_roi_t *const roi_in,
                              const dt_iop_roi_t *const roi_out,
                              const int in_bpp)
{
  const gboolean use_roi = memcmp(roi_in, roi_out, sizeof(struct dt_iop_roi_t)) || (self->flags() & IOP_FLAGS_TILING_FULL_ROI);
  if(use_roi)
    return _default_process_tiling_cl_roi(self, piece, ivoid, ovoid, roi_in, roi_out, in_bpp);
  else
    return _default_process_tiling_cl_ptp(self, piece, ivoid, ovoid, roi_in, roi_out, in_bpp);
}

#else
int default_process_tiling_cl(dt_iop_module_t *self,
                              dt_dev_pixelpipe_iop_t *piece,
                              const void *const ivoid,
                              void *const ovoid,
                              const dt_iop_roi_t *const roi_in,
                              const dt_iop_roi_t *const roi_out,
                              const int in_bpp)
{
  return DT_OPENCL_DEFAULT_ERROR;
}
#endif


/* If a module does not implement tiling_callback() by itself, this function is called instead.
   Default is an image size factor of 2 (i.e. input + output buffer needed), no overhead (1),
   no overlap between tiles, and an pixel alignment of 1 in x and y direction, i.e. no special
   alignment required. Simple pixel to pixel modules (take tonecurve as an example) can happily
   live with that.
   (1) Small overhead like look-up-tables can be ignored safely. */
void default_tiling_callback(dt_iop_module_t *self,
                             dt_dev_pixelpipe_iop_t *piece,
                             const dt_iop_roi_t *roi_in,
                             const dt_iop_roi_t *roi_out,
                             struct dt_develop_tiling_t *tiling)
{
  const float ioratio
      = ((float)roi_out->width * (float)roi_out->height) / ((float)roi_in->width * (float)roi_in->height);

  tiling->factor = 1.0f + ioratio;
  tiling->factor_cl = tiling->factor;
  tiling->maxbuf = 1.0f;
  tiling->maxbuf_cl = tiling->maxbuf;
  tiling->overhead = 0;
  tiling->overlap = 0;
  tiling->align = 1;

  if((self->flags() & IOP_FLAGS_TILING_FULL_ROI) == IOP_FLAGS_TILING_FULL_ROI) tiling->overlap = 4;
}

gboolean dt_tiling_piece_fits_host_memory(const dt_dev_pixelpipe_iop_t *piece,
                                          const size_t width,
                                          const size_t height,
                                          const unsigned bpp,
                                          const float factor,
                                          const size_t overhead)
{
  const size_t available = dt_get_available_pipe_mem(piece->pipe);
  const size_t total = factor * width * height * bpp + overhead;

  return total <= available;
}

// 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