#include "develop/tiling.h"

#include "develop/pixelpipe.h"

#include "develop/blend.h"

#include "common/opencl.h"

#include "control/control.h"

#include <string.h>

#include <strings.h>

#include <stdlib.h>

#include <math.h>

#include <unistd.h>

#include <assert.h>

#define CLAMPI(a, mn, mx) ((a) < (mn) ? (mn) : ((a) > (mx) ? (mx) : (a)))

#define CL_ALIGNMENT 4

#define RESERVE 5

static unsigned

_gcd(unsigned a, unsigned b)

{

unsigned t;

while(b != 0)

{

t = b;

b = a % b;

a = t;

}

return a;

}

static unsigned

_lcm(unsigned a, unsigned b)

{

return (((unsigned long)a * b) / _gcd(a, b));

}

static inline int

_min(int a, int b)

{

return a < b ? a : b;

}

static inline int

_max(int a, int b)

{

return a > b ? a : b;

}

static inline int

_align_up(int n, int a)

{

return n % a !=0 ? (n/a + 1) * a : n;

}

static inline int

_align_down(int n, int a)

{

return n % a !=0 ? (n/a) * a : n;

}

_print_roi(const dt_iop_roi_t *roi, const char *label)

{

printf("{ %5d %5d %5d %5d %.6f } %s\n", roi->x, roi->y, roi->width, roi->height, roi->scale, label);

}

#if 0

static void

_nm_constraints(double x[], int n)

{

x[0] = fabs(x[0]);

x[1] = fabs(x[1]);

x[2] = fabs(x[2]);

x[3] = fabs(x[3]);

if(x[0] > 1.0) x[0] = 1.0 - x[0];

if(x[1] > 1.0) x[1] = 1.0 - x[1];

if(x[2] > 1.0) x[2] = 1.0 - x[2];

if(x[3] > 1.0) x[3] = 1.0 - x[3];

}

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;

}

#define MAX_IT 1000 /* maximum number of iterations */

#define ALPHA 1.0 /* reflection coefficient */

#define BETA 0.5 /* contraction coefficient */

#define GAMMA 2.0 /* expansion coefficient */

static int

_simplex(double (*objfunc) (double[], void*[]), double start[], int n, double EPSILON,

double scale, 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, m, row;

int k; /* track the number of function evaluations */

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 ((n + 1) * sizeof (double *));

f = (double *) malloc ((n + 1) * sizeof (double));

vr = (double *) malloc (n * sizeof (double));

ve = (double *) malloc (n * sizeof (double));

vc = (double *) malloc (n * sizeof (double));

vm = (double *) malloc (n * sizeof (double));

/* allocate the columns of the arrays */

for (i = 0; i <= n; i++)

{

v[i] = (double *) malloc (n * sizeof (double));

}

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

}

k = n + 1;

#if 0

/* 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]);

}

}

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

k++;

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

k++;

/* by making fe < fr as opposed to fe < f[vs],

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

k++;

}

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

k++;

}

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,

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

k++;

if (constrain != NULL)

{

constrain (v[vh], n);

}

f[vh] = objfunc (v[vh], rest);

k++;

}

}

#if 0

/* 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]);

}

}

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

}

}

#if 0

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

free (f);

free (vr);

free (ve);

free (vc);

free (vm);

for (i = 0; i <= n; i++)

{

free (v[i]);

}

free (v);

return itr;

}

static int

_nm_fit_output_to_input_roi(struct dt_iop_module_t *self, struct dt_dev_pixelpipe_iop_t *piece, const dt_iop_roi_t *iroi, dt_iop_roi_t *oroi, 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 };

double epsilon = (double)delta / MIN(piece->iwidth, piece->iheight);

int maxiter = 1000;

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

//printf("_simplex: %d, delta: %d, epsilon: %f\n", 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);

}

static int

_fit_output_to_input_roi(struct dt_iop_module_t *self, struct dt_dev_pixelpipe_iop_t *piece, const dt_iop_roi_t *iroi, dt_iop_roi_t *oroi, 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.

int fit = _nm_fit_output_to_input_roi(self, piece, iroi, oroi, delta);

return fit;

}

static void

_default_process_tiling_ptp (struct dt_iop_module_t *self, struct dt_dev_pixelpipe_iop_t *piece, void *ivoid, void *ovoid, const dt_iop_roi_t *roi_in, const dt_iop_roi_t *roi_out, const int in_bpp)

{

void *input = NULL;

void *output = NULL;

const int out_bpp = self->output_bpp(self, piece->pipe, piece);

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

self->tiling_callback(self, piece, roi_in, roi_out, &tiling);

/* 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_DEV, "[default_process_tiling_ptp] no need to use tiling for module '%s' as no real memory saving to be expected\n", self->op);

goto fallback;

}

/* calculate optimal size of tiles */

float available = (float)dt_conf_get_int("host_memory_limit")*1024.0f*1024.0f;

assert(available >= 500.0f*1024.0f*1024.0f);

/* correct for size of ivoid and ovoid which are needed on top of tiling */

available = fmax(available - (roi_out->width*roi_out->height*out_bpp) - (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.

float singlebuffer = (float)dt_conf_get_int("singlebuffer_limit")*1024.0f*1024.0f;

singlebuffer = fmax(singlebuffer, 2.0f*1024.0f*1024.0f);

float factor = fmax(tiling.factor, 1.0f);

float maxbuf = fmax(tiling.maxbuf, 1.0f);

singlebuffer = fmax(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/(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 * sqrt(scale));

height = floorf(height * sqrt(scale));

}

}

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

}

/* Alignment rules: we need to make sure that alignment requirements of module are fulfilled.

const unsigned int xyalign = _lcm(tiling.xalign, tiling.yalign);

assert(xyalign != 0);

/* properly align tile width and height by making them smaller if needed */

if(width < roi_in->width) width = (width / xyalign) * xyalign;

if(height < roi_in->height) height = (height / xyalign) * xyalign;

/* also make sure that overlap follows alignment rules by making it wider when needed */

const int overlap = tiling.overlap % xyalign != 0 ? (tiling.overlap / xyalign + 1) * xyalign : 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 > DT_TILING_MAXTILES)

{

dt_print(DT_DEBUG_DEV, "[default_process_tiling_ptp] gave up tiling for module '%s'. too many tiles: %d x %d\n", self->op, tiles_x, tiles_y);

goto error;

}

dt_print(DT_DEBUG_DEV, "[default_process_tiling_ptp] use tiling on module '%s' for image with full size %d x %d\n", self->op, roi_in->width, roi_in->height);

dt_print(DT_DEBUG_DEV, "[default_process_tiling_ptp] (%d x %d) tiles with max dimensions %d x %d and overlap %d\n", tiles_x, tiles_y, width, height, overlap);

/* reserve input and output buffers for tiles */

input = dt_alloc_align(64, width*height*in_bpp);

if(input == NULL)

{

dt_print(DT_DEBUG_DEV, "[default_process_tiling_ptp] could not alloc input buffer for module '%s'\n", self->op);

goto error;

}

output = dt_alloc_align(64, width*height*out_bpp);

if(output == NULL)

{

dt_print(DT_DEBUG_DEV, "[default_process_tiling_ptp] could not alloc output buffer for module '%s'\n", self->op);

goto error;

}

/* store processed_maximum to be re-used and aggregated */

float processed_maximum_saved[3];

float processed_maximum_new[3] = { 1.0f };

for(int k=0; k<3; k++)

processed_maximum_saved[k] = piece->pipe->processed_maximum[k];

/* iterate over tiles */

for(int tx=0; tx<tiles_x; tx++)

for(int ty=0; ty<tiles_y; ty++)

{

piece->pipe->tiling = 1;

size_t wd = tx * tile_wd + width > roi_in->width ? roi_in->width - tx * tile_wd : width;

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

if((wd <= overlap && tx > 0) || (ht <= overlap && ty > 0)) continue;

/* origin and region of effective part of tile, which we want to store later */

size_t origin[] = { 0, 0, 0 };

size_t region[] = { wd, ht, 1 };

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

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(DT_DEBUG_DEV, "[default_process_tiling_ptp] tile (%d, %d) with %d x %d at origin [%d, %d]\n", tx, ty, wd, ht, tx*tile_wd, ty*tile_ht);

/* prepare input tile buffer */

#ifdef _OPENMP

#pragma omp parallel for default(none) shared(input,width,ivoid,ioffs,wd,ht) schedule(static)

for(int j=0; j<ht; j++)

memcpy((char *)input+j*wd*in_bpp, (char *)ivoid+ioffs+j*ipitch, wd*in_bpp);

/* take original processed_maximum as starting point */

for(int k=0; k<3; k++)

piece->pipe->processed_maximum[k] = processed_maximum_saved[k];

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

for(int k=0; k<3; k++)

{

if(tx+ty > 0 && fabs(processed_maximum_new[k] - piece->pipe->processed_maximum[k]) > 1.0e-6f)

dt_print(DT_DEBUG_DEV, "[default_process_tiling_ptp] processed_maximum[%d] differs between tiles in module '%s'\n", k, self->op);

processed_maximum_new[k] = piece->pipe->processed_maximum[k];

}

/* correct origin and region of tile for overlap.

if(tx > 0)

{

origin[0] += overlap;

region[0] -= overlap;

ooffs += overlap*out_bpp;

}

if(ty > 0)

{

origin[1] += overlap;

region[1] -= overlap;

ooffs += overlap*opitch;

}

/* copy "good" part of tile to output buffer */

#ifdef _OPENMP

#pragma omp parallel for default(none) shared(ovoid,ooffs,output,width,origin,region,wd) schedule(static)

for(int j=0; j<region[1]; j++)

memcpy((char *)ovoid+ooffs+j*opitch, (char *)output+((j+origin[1])*wd+origin[0])*out_bpp, region[0]*out_bpp);

}

/* copy back final processed_maximum */

for(int k=0; k<3; k++)

piece->pipe->processed_maximum[k] = processed_maximum_new[k];

if(input != NULL) free(input);

if(output != NULL) free(output);

piece->pipe->tiling = 0;

return;

error:

dt_control_log(_("tiling failed for module '%s'. output might be garbled."), self->op);

// fall through

fallback:

if(input != NULL) free(input);

if(output != NULL) free(output);

piece->pipe->tiling = 0;

dt_print(DT_DEBUG_DEV, "[default_process_tiling_ptp] fall back to standard processing for module '%s'\n", self->op);

self->process(self, piece, ivoid, ovoid, roi_in, roi_out);

return;

}

static void

_default_process_tiling_roi (struct dt_iop_module_t *self, struct dt_dev_pixelpipe_iop_t *piece, void *ivoid, void *ovoid, const dt_iop_roi_t *roi_in, const dt_iop_roi_t *roi_out, const int in_bpp)

{

void *input = NULL;

void *output = NULL;

//_print_roi(roi_in, "module roi_in");

//_print_roi(roi_out, "module roi_out");

const int out_bpp = self->output_bpp(self, piece->pipe, piece);

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 = fmax(roi_in->scale / roi_out->scale, sqrt((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 };

self->tiling_callback(self, piece, roi_in, roi_out, &tiling);

/* 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_DEV, "[default_process_tiling_roi] no need to use tiling for module '%s' as no real memory saving to be expected\n", self->op);

goto fallback;

}

/* calculate optimal size of tiles */

float available = (float)dt_conf_get_int("host_memory_limit")*1024.0f*1024.0f;

assert(available >= 500.0f*1024.0f*1024.0f);

/* correct for size of ivoid and ovoid which are needed on top of tiling */

available = fmax(available - (roi_out->width*roi_out->height*out_bpp) - (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.

float singlebuffer = (float)dt_conf_get_int("singlebuffer_limit")*1024.0f*1024.0f;

singlebuffer = fmax(singlebuffer, 2.0f*1024.0f*1024.0f);

float factor = fmax(tiling.factor, 1.0f);

float maxbuf = fmax(tiling.maxbuf, 1.0f);

singlebuffer = fmax(available / factor, singlebuffer);

int width = _max(roi_in->width, roi_out->width);

int height = _max(roi_in->height, roi_out->height);

/* shrink tile size in case it would exceed singlebuffer size */

if((float)width*height*max_bpp*maxbuf > singlebuffer)

{

const float scale = singlebuffer/(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 * sqrt(scale));

height = floorf(height * sqrt(scale));

}

}

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

}

/* Alignment rules: we need to make sure that alignment requirements of module are fulfilled.

/* for simplicity reasons we use only one alignment that fits to x and y requirements at the same time */

unsigned int xyalign = _lcm(tiling.xalign, tiling.yalign);

assert(xyalign != 0);

/* make sure that overlap follows alignment rules by making it wider when needed.

const int overlap_in = _align_up(tiling.overlap, xyalign);

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.

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 > DT_TILING_MAXTILES)

{

dt_print(DT_DEBUG_DEV, "[default_process_tiling_roi] gave up tiling for module '%s'. too many tiles: %d x %d\n", self->op, tiles_x, tiles_y);

goto error;

}

/* calculate tile width and height excl. overlap (i.e. the good part) for output.

const int tile_wd = _align_up(roi_out->width % tiles_x == 0 ? roi_out->width / tiles_x : roi_out->width / tiles_x + 1, xyalign);

const int tile_ht = _align_up(roi_out->height % tiles_y == 0 ? roi_out->height / tiles_y : roi_out->height / tiles_y + 1, xyalign);

dt_print(DT_DEBUG_DEV, "[default_process_tiling_roi] use tiling on module '%s' for image with full input size %d x %d\n", self->op, roi_in->width, roi_in->height);

dt_print(DT_DEBUG_DEV, "[default_process_tiling_roi] (%d x %d) tiles with max dimensions %d x %d\n", tiles_x, tiles_y, width, height);

/* store processed_maximum to be re-used and aggregated */

float processed_maximum_saved[3];

float processed_maximum_new[3] = { 1.0f };

for(int k=0; k<3; k++)

processed_maximum_saved[k] = piece->pipe->processed_maximum[k];

/* iterate over tiles */

for(int tx=0; tx<tiles_x; tx++)

for(int ty=0; ty<tiles_y; ty++)

{

piece->pipe->tiling = 1;

/* the output dimensions of the good part of this specific tile */

size_t wd = (tx + 1) * tile_wd > roi_out->width ? roi_out->width - tx * tile_wd : tile_wd;

size_t ht = (ty + 1) * tile_ht > roi_out->height ? 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.

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

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_down(x_in - overlap_in - delta, xyalign), roi_in->x);

const int new_y_in = _max(_align_down(y_in - overlap_in - delta, xyalign), roi_in->y);

const int new_width_in = _min(_align_up(width_in + overlap_in + delta + (x_in - new_x_in), xyalign), 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), xyalign), 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))