blurs.c

Go to the documentation of this file.

/*
    This file is part of darktable,
    Copyright (C) 2021, 2023, 2025-2026 Aurélien PIERRE.
    Copyright (C) 2021 Hubert Kowalski.
    Copyright (C) 2021-2022 Pascal Obry.
    Copyright (C) 2021 Ralf Brown.
    Copyright (C) 2022 Diederik Ter Rahe.
    Copyright (C) 2022 Hanno Schwalm.
    Copyright (C) 2022 Martin Bařinka.
    Copyright (C) 2022 Philipp Lutz.
    Copyright (C) 2023 Luca Zulberti.
    Copyright (C) 2024 Alynx Zhou.
    Copyright (C) 2026 Guillaume Stutin.
    
    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/>.
*/
#ifdef HAVE_CONFIG_H
#include "system/macros.h"
#include "system/mem_alloc.h"
#include "common/module_versioning.h"
#include "common/logging.h"
#include "system/openmp.h"
#include "system/simd.h"
#include "system/target_clones.h"
#include "config.h"
#endif
// our includes go first:
#include "widgets/bauhaus.h"
#include "develop/imageop.h"
#include "develop/imageop_gui.h"
#include "widgets/drawingarea.h"
#include "iop/iop_api.h"
 
// #include <fftw3.h> // one day, include FFT convolution
#include <gtk/gtk.h>
#include <stdlib.h>
 
DT_MODULE_INTROSPECTION(1, dt_iop_blurs_params_t)
 
typedef enum dt_iop_blur_type_t
{
  DT_BLUR_LENS = 0,     // $DESCRIPTION: "lens"
  DT_BLUR_MOTION = 1,   // $DESCRIPTION: "motion"
  DT_BLUR_GAUSSIAN = 2, // $DESCRIPTION: "gaussian"
} dt_iop_blur_type_t;
 
 
typedef struct dt_iop_blurs_params_t
{
  dt_iop_blur_type_t type; // $DEFAULT: DT_BLUR_LENS $DESCRIPTION: "blur type"
  int radius;              // $MIN: 4 $MAX: 128 $DEFAULT: 8 $DESCRIPTION: "blur radius"
 
  // lens blur params
  int blades;              // $MIN: 3 $MAX: 11 $DEFAULT: 5 $DESCRIPTION: "diaphragm blades"
  float concavity;         // $MIN: 1. $MAX: 9.  $DEFAULT: 1. $DESCRIPTION: "concavity"
  float linearity;         // $MIN: 0. $MAX: 1.  $DEFAULT: 1. $DESCRIPTION: "linearity"
  float rotation;          // $MIN: -1.57 $MAX: 1.57 $DEFAULT: 0. $DESCRIPTION: "rotation"
 
  // motion blur params
  float angle;             // $MIN: -3.14 $MAX: 3.14 $DEFAULT: 0. $DESCRIPTION: "direction"
  float curvature;         // $MIN: -2.   $MAX: 2.   $DEFAULT: 0. $DESCRIPTION: "curvature"
  float offset;            // $MIN: -1.   $MAX: 1.   $DEFAULT: 0  $DESCRIPTION: "offset"
 
} dt_iop_blurs_params_t;
 
 
typedef struct dt_iop_blurs_gui_data_t
{
  GtkWidget *type, *radius, *blades, *concavity, *linearity, *rotation, *angle, *curvature, *offset;
  GtkDrawingArea *area;
  unsigned char *img;
  int img_cached;
  float img_width;
} dt_iop_blurs_gui_data_t;
 
 
typedef struct dt_iop_blurs_global_data_t
{
  int kernel_blurs_convolve;
} dt_iop_blurs_global_data_t;
 
 
const char *name()
{
  return _("blurs");
}
 
const char *aliases()
{
  return _("blur|lens|motion");
}
 
const char **description(struct dt_iop_module_t *self)
{
  return dt_iop_set_description(self,
                                _("simulate physically-accurate lens and motion blurs"),
                                _("creative"), _("linear, RGB, scene-referred"), _("linear, RGB"),
                                _("linear, RGB, scene-referred"));
}
 
int flags()
{
  return IOP_FLAGS_INCLUDE_IN_STYLES | IOP_FLAGS_SUPPORTS_BLENDING;
}
 
 
int default_group()
{
  return IOP_GROUP_SHARPNESS;
}
 
 
int default_colorspace(dt_iop_module_t *self, dt_dev_pixelpipe_t *pipe, const dt_dev_pixelpipe_iop_t *piece)
{
  return IOP_CS_RGB;
}
 
 
void commit_params(dt_iop_module_t *self, dt_iop_params_t *p1, dt_dev_pixelpipe_t *pipe, dt_dev_pixelpipe_iop_t *piece)
{
  memcpy(piece->data, p1, self->params_size);
}
 
// B spline filter
#define FSIZE 5
 
inline static void blur_2D_Bspline(const float *const restrict in, float *const restrict out,
                                   const size_t width, const size_t height)
{
  __OMP_PARALLEL_FOR__( collapse(2))
  for(size_t i = 0; i < height; i++)
  {
    for(size_t j = 0; j < width; j++)
    {
      const size_t index = (i * width + j);
      float acc = 0.f;
 
      for(size_t ii = 0; ii < FSIZE; ++ii)
        for(size_t jj = 0; jj < FSIZE; ++jj)
        {
          const size_t row = CLAMP((int)i + (int)(ii - (FSIZE - 1) / 2), (int)0, (int)height - 1);
          const size_t col = CLAMP((int)j + (int)(jj - (FSIZE - 1) / 2), (int)0, (int)width - 1);
          const size_t k_index = (row * width + col);
 
          static const float DT_ALIGNED_ARRAY filter[FSIZE]
              = { 1.0f / 16.0f, 4.0f / 16.0f, 6.0f / 16.0f, 4.0f / 16.0f, 1.0f / 16.0f };
 
          acc += filter[ii] * filter[jj] * in[k_index];
        }
 
      out[index] = acc;
    }
  }
}
 
 
__DT_CLONE_TARGETS__
static inline void init_kernel(float *const restrict buffer, const size_t width, const size_t height)
{
  // init an empty kernel with zeros
  __OMP_PARALLEL_FOR_SIMD__(aligned(buffer:64))
  for(size_t k = 0; k < height * width; k++) buffer[k] = 0.f;
}
 
 
__DT_CLONE_TARGETS__
static inline void create_lens_kernel(float *const restrict buffer, const size_t width,
                                      const size_t height, const float n, const float m,
                                      const float k, const float rotation)
{
  // n is number of diaphragm blades
  // m is the concavity, aka the number of vertices on straight lines (?)
  // k is the roundness vs. linearity factor
  //   see https://math.stackexchange.com/a/4160104/498090
  // buffer sizes need to be odd
 
  // Spatial coordinates rounding error
  const float eps = 1.f / (float)width;
  const float radius = (float)(width - 1) / 2.f - 1;
  __OMP_PARALLEL_FOR_SIMD__(aligned(buffer:64) collapse(2))
  for(size_t i = 0; i < height; i++)
    for(size_t j = 0; j < width; j++)
    {
      // get normalized kernel coordinates in [-1 ; 1]
      const float x = (float)(i - 1) / radius - 1;
      const float y = (float)(j - 1) / radius - 1;
 
      // get current radial distance from kernel center
      const float r = dt_fast_hypotf(x, y);
 
      // get the radial distance at current angle of the shape envelope
      const float M = cosf((2.f * asinf(k) + M_PI_F * m) / (2.f * n))
                      / cosf((2.f * asinf(k * cosf(n * (atan2f(y, x) + rotation))) + M_PI_F * m) / (2.f * n));
 
      // write 1 if we are inside the envelope of the shape, else 0
      buffer[i * width + j] = (M >= r + eps);
    }
}
 
 
__DT_CLONE_TARGETS__
static inline void create_motion_kernel(float *const restrict buffer, const size_t width,
                                        const size_t height, const float angle,
                                        const float curvature, const float offset)
{
  // Compute the polynomial params from user params
  const float A = curvature / 2.f;
  const float B = 1.f;
  const float C = -A * offset * offset + B * offset;
  // Note : C ensures the polynomial arc always goes through the central pixel
  // so we don't shift pixels. This is meant to allow seamless connection
  // with unmasked areas when using masked blur.
 
  // Spatial coordinates rounding error
  const float eps = 1.f / (float)width;
 
  const float radius = (float)(width - 1) / 2.f - 1;
  const float corr_angle = -M_PI_F / 4.f - angle;
 
  // Matrix of rotation
  const float M[2][2] = { { cosf(corr_angle), -sinf(corr_angle) },
                          { sinf(corr_angle), cosf(corr_angle) } };
  __OMP_PARALLEL_FOR_SIMD__(aligned(buffer:64))
  for(size_t i = 0; i < 8 * width; i++)
  {
    // Note : for better smoothness of the polynomial discretization,
    // we oversample 8 times, meaning we evaluate the polynomial
    // every eighth of pixel
 
    // get normalized kernel coordinates in [-1 ; 1]
    const float x = (float)(i / 8.f - 1) / radius - 1;
    //const float y = (j - 1) / radius - 1; // not used here
 
    // build the motion path : 2nd order polynomial
    const float X = x - offset;
    const float y = X * X * A + X * B + C;
 
    // rotate the motion path around the kernel center
    const float rot_x = x * M[0][0] + y * M[0][1];
    const float rot_y = x * M[1][0] + y * M[1][1];
 
    // convert back to kernel absolute coordinates ± eps
    const int y_f[2] = { roundf((rot_y + 1) * radius - eps),
                         roundf((rot_y + 1) * radius + eps) };
    const int x_f[2] = { roundf((rot_x + 1) * radius - eps),
                         roundf((rot_x + 1) * radius + eps) };
 
    // write 1 if we are inside the envelope of the shape, else 0
    // leave 1px padding on each border of the kernel for the anti-aliasing
    for(int l = 0; l < 2; l++)
      for(int m = 0; m < 2; m++)
      {
        if(x_f[l] > 0 && x_f[l] < width - 1 && y_f[m] > 0 && y_f[m] < width - 1)
          buffer[y_f[m] * width + x_f[l]] = 1.f;
      }
  }
}
 
 
__DT_CLONE_TARGETS__
static inline void create_gauss_kernel(float *const restrict buffer, const size_t width,
                                       const size_t height)
{
  // This is not optimized. Gauss kernel is separable and can be turned into
  // 2 x 1D convolutions.
  const float radius = (width - 1) / 2.f - 1;
  __OMP_PARALLEL_FOR_SIMD__(aligned(buffer:64) collapse(2))
  for(size_t i = 0; i < height; i++)
    for(size_t j = 0; j < width; j++)
    {
      // get normalized kernel coordinates in [-1 ; 1]
      const float x = (float)(i - 1) / radius - 1;
      const float y = (float)(j - 1) / radius - 1;
 
      // get current square radial distance from kernel center
      const float r_2 = x * x + y * y;
      buffer[i * width + j] = expf(-4.f * r_2);
    }
}
 
 
 
__DT_CLONE_TARGETS__
static inline int build_gui_kernel(unsigned char *const buffer, const size_t width,
                                   const size_t height, dt_iop_blurs_params_t *p)
{
  float *const restrict kernel_1 = dt_alloc_align_float(width * height);
  float *const restrict kernel_2 = dt_alloc_align_float(width * height);
  if(IS_NULL_PTR(kernel_1) || IS_NULL_PTR(kernel_2)) goto error;
 
  if(p->type == DT_BLUR_LENS)
  {
    create_lens_kernel(kernel_1, width, height, p->blades, p->concavity, p->linearity, p->rotation);
 
    // anti-aliasing step
    blur_2D_Bspline(kernel_1, kernel_2, width, height);
  }
  else if(p->type == DT_BLUR_MOTION)
  {
    init_kernel(kernel_1, width, height);
    create_motion_kernel(kernel_1, width, height, p->angle, p->curvature, p->offset);
 
    // anti-aliasing step
    blur_2D_Bspline(kernel_1, kernel_2, width, height);
  }
  else if(p->type == DT_BLUR_GAUSSIAN)
  {
    create_gauss_kernel(kernel_2, width, height);
  }
 
  // Convert to Gtk/Cairo RGBA 8x4 bits
  __OMP_PARALLEL_FOR_SIMD__(aligned(buffer, kernel_2:64))
  for(size_t k = 0; k < height * width; k++)
  {
    buffer[k * 4] = buffer[k * 4 + 1] = buffer[k * 4 + 2] = buffer[k * 4 + 3] = roundf(255.f * kernel_2[k]);
  }
 
error:;
  int err = (IS_NULL_PTR(kernel_1) || IS_NULL_PTR(kernel_2));
  dt_free_align(kernel_1);
  dt_free_align(kernel_2);
  return err;
}
 
 
__DT_CLONE_TARGETS__
static inline float compute_norm(float *const buffer, const size_t width, const size_t height)
{
  float norm = 0.f;
  __OMP_PARALLEL_FOR_SIMD__(aligned(buffer:64) reduction(+:norm))
  for(size_t i = 0; i < width * height; i++)
  {
    norm += buffer[i];
  }
 
  return norm;
}
 
 
__DT_CLONE_TARGETS__
static inline void normalize(float *const buffer, const size_t width, const size_t height,
                             const float norm)
{
  __OMP_PARALLEL_FOR_SIMD__(aligned(buffer:64))
  for(size_t i = 0; i < width * height; i++)
  {
    buffer[i] /= norm;
  }
}
 
 
static inline int build_pixel_kernel(float *const buffer, const size_t width, const size_t height,
                                     dt_iop_blurs_params_t *p)
{
  float *const restrict kernel_1 = dt_alloc_align_float(width * height);
  if(IS_NULL_PTR(kernel_1)) return 1;
 
  if(p->type == DT_BLUR_LENS)
  {
    create_lens_kernel(kernel_1, width, height, p->blades, p->concavity, p->linearity, p->rotation + M_PI_F);
 
    // anti-aliasing step
    blur_2D_Bspline(kernel_1, buffer, width, height);
  }
  else if(p->type == DT_BLUR_MOTION)
  {
    init_kernel(kernel_1, width, height);
    create_motion_kernel(kernel_1, width, height, p->angle + M_PI_F, p->curvature, p->offset);
 
    // anti-aliasing step
    blur_2D_Bspline(kernel_1, buffer, width, height);
  }
  else if(p->type == DT_BLUR_GAUSSIAN)
  {
    create_gauss_kernel(buffer, width, height);
  }
 
  // normalize to respect the conservation of energy law
  const float norm = compute_norm(buffer, width, height);
  normalize(buffer, width, height, norm);
 
  dt_free_align(kernel_1);
  return 0;
}
 
#if 0
 
// This crashes on the FFT step - not sure why and no time to investigate opaque libs now
// FFT convolution should be faster for large blurs because it is o(N log2(N))
// where N is the width of the kernel
// TODO
 
static void process_fft(struct dt_iop_module_t *self, dt_dev_pixelpipe_iop_t *piece,
                        const void *const restrict ivoid, void *const restrict ovoid,
                        const dt_iop_roi_t *const roi_in, const dt_iop_roi_t *const roi_out)
{
  dt_iop_blurs_params_t *p = (dt_iop_blurs_params_t *)piece->data;
  const float scale = dt_dev_get_module_scale(pipe, roi_in);
 
  const float *const restrict in = __builtin_assume_aligned(ivoid, 64);
  //float *const restrict out = __builtin_assume_aligned(ovoid, 64);
 
  // FFT needs odd buffer sizes, so fix that here
  const int is_width_even = (roi_in->width % 2 == 0);
  const int is_height_even = (roi_in->height % 2 == 0);
  const size_t padded_width = roi_in->width + 1 * is_width_even;
  const size_t padded_height = roi_in->height + 1 * is_height_even;
 
  float *const restrict padded_in = dt_alloc_align_float(padded_width * padded_height * 4);
  float *const restrict padded_out = dt_alloc_align_float(padded_width * padded_height * 4);
 
  // Write the image in the padded buffer
  __OMP_PARALLEL_FOR_SIMD__(aligned(in, padded_in:64))
  for(size_t i = 0; i < roi_in->height; i++)
    for(size_t j = 0; j < roi_in->width; j++)
    {
      const size_t index_in = (i * roi_in->width + j) * 4;
      const size_t index_out = (i * padded_width + j) * 4;
      for_four_channels(c, aligned(in, padded_in : 64)) padded_in[index_out + c] = in[index_in + c];
    }
 
  // Write the padding if needed
  if(padded_width > roi_in->width)
  {
  __OMP_PARALLEL_FOR_SIMD__(aligned(in, padded_in:64))
  for(size_t i = 0; i < roi_in->height; i++)
    {
      const size_t index_in = (i * (roi_in->width - 1)) * 4;
      const size_t index_out = (i * (padded_width - 1)) * 4;
      for_four_channels(c, aligned(in, padded_in : 64)) padded_in[index_out + c] = in[index_in + c];
    }
  }
 
  if(padded_height > roi_in->height)
  {
  __OMP_PARALLEL_FOR_SIMD__(aligned(in, padded_in:64))
  for(size_t j = 0; j < roi_in->width; j++)
    {
      const size_t index_in = ((roi_in->height - 1) * roi_in->width + j) * 4;
      const size_t index_out = ((padded_height - 1) * padded_width + j) * 4;
      for_four_channels(c, aligned(in, padded_in : 64)) padded_in[index_out + c] = in[index_in + c];
    }
  }
 
  // Init the blur kernel
  const size_t radius = MAX(roundf(p->radius / scale), 1);
  const size_t kernel_width = 2 * radius + 1;
 
  float *const restrict kernel = dt_alloc_align_float(kernel_width * kernel_width);
  build_pixel_kernel(kernel, kernel_width, kernel_width, p);
 
  // Convert to padded kernel - copy kernel in the center
  float *const restrict padded_kernel = dt_alloc_align_float(padded_width * padded_height);
  const size_t offset_i = (padded_height - 1) / 2 - (kernel_width - 1) / 2;
  const size_t offset_j = (padded_width - 1) / 2 - (kernel_width - 1) / 2;
  const size_t i_reach = offset_i + kernel_width;
  const size_t j_reach = offset_j + kernel_width;
  __OMP_PARALLEL_FOR_SIMD__(aligned(kernel, padded_kernel:64))
  for(size_t i = 0; i < padded_width; i++)
    for(size_t j = 0; j < padded_width; j++)
    {
      const size_t padded_idx = (i * padded_width) + j;
 
      if(i >= offset_i && i < i_reach && j >= offset_j && j < j_reach)
      {
        const size_t i_kern = i - offset_i;
        const size_t j_kern = j - offset_j;
        const size_t kern_idx = i_kern * kernel_width + j_kern;
        padded_kernel[padded_idx] = kernel[kern_idx];
      }
      else
        padded_kernel[padded_idx] = 0.f;
    }
 
  // Init the FFT transforms
  int threads = fftw_init_threads();
  fftw_plan_with_nthreads(threads);
 
  // TODO: things go well until this point
 
  // Plan the dimensions of the FFT
  // notice we use fftwf prefix to use the float 32 variant of fftw
  int rank = 2;     /* 2D FFT */
  int n[2] = { padded_width, padded_height };
  int howmany = 4; /* 4 channels : RGBa */
  int idist = 1;   /* channels are distanced by one float in memory */
  int odist = 1;   /* channels are distanced by one float in memory */
  int istride = 4; /* array is not contiguous in memory */
  int ostride = 4; /* array is not contiguous in memory */
  int *inembed = n;
  int *onembed = n;
 
  fftwf_complex *const restrict kernel_fft = fftwf_alloc_complex(padded_height * padded_height * 4);
  fftwf_complex *const restrict image_fft = fftwf_alloc_complex(padded_height * padded_height * 4);
 
  // FFT convert the kernel
  fftwf_plan kernel_plan = fftwf_plan_many_dft_r2c(rank, n, howmany,
                                 padded_kernel, inembed,
                                 istride, idist,
                                 kernel_fft, onembed,
                                 ostride, odist,
                                 FFTW_ESTIMATE);
  fftwf_execute(kernel_plan);
 
  // Clean FFT
  fftwf_destroy_plan(kernel_plan);
  fftwf_free(kernel_fft);
  fftwf_free(image_fft);
  fftw_cleanup_threads();
 
  dt_free_align(kernel);
  dt_free_align(padded_kernel);
  dt_free_align(padded_in);
  dt_free_align(padded_out);
}
#endif
 
// Spatial convolution should be slower for large blurs because it is o(N²) where N is the width of the kernel
// but code is much simpler and easier to debug
 
__DT_CLONE_TARGETS__
int process(struct dt_iop_module_t *self, const dt_dev_pixelpipe_t *pipe, const dt_dev_pixelpipe_iop_t *piece,
                    const void *const restrict ivoid, void *const restrict ovoid)
{
  const dt_iop_roi_t *const roi_in = &piece->roi_in;
  const dt_iop_roi_t *const roi_out = &piece->roi_out;
  dt_iop_blurs_params_t *p = (dt_iop_blurs_params_t *)piece->data;
  const float scale = dt_dev_get_module_scale(pipe, roi_in);
 
  const float *const restrict in = __builtin_assume_aligned(ivoid, 64);
  float *const restrict out = __builtin_assume_aligned(ovoid, 64);
 
  // Init the blur kernel
  const int radius = MAX(roundf(p->radius / scale), 2);
  const size_t kernel_width = 2 * radius + 1;
 
  float *restrict kernel = dt_alloc_align_float(kernel_width * kernel_width);
  if(IS_NULL_PTR(kernel)) return 1;
  if(build_pixel_kernel(kernel, kernel_width, kernel_width, p))
  {
    dt_free_align(kernel);
    return 1;
  }
  // The pipe abandoned this run: the rows not started yet are skipped and the module fails. One
  // convolution of up to (2 * 128 + 1)² taps per pixel at full resolution has no other step to stop
  // at than its rows. A new request may lower the killswitch again before the loop ends, hence the
  // latch, shared since the macro makes every variable firstprivate.
  dt_atomic_int abandoned = FALSE;
  __OMP_PARALLEL_FOR__(shared(abandoned))
  for(int i = 0; i < roi_out->height; i++)
  {
    if(dt_atomic_get_int(&abandoned) || dt_dev_pixelpipe_has_shutdown(pipe))
    {
      dt_atomic_set_int(&abandoned, TRUE);
      continue;
    }
    for(int j = 0; j < roi_out->width; j++)
    {
      const size_t index = ((i * roi_out->width) + j) * 4;
      float DT_ALIGNED_PIXEL acc[4] = { 0.f };
 
      if(i >= radius && j >= radius && i < roi_out->height - radius && j < roi_out->width - radius)
      {
        // We are in the safe area, no need to check for out-of-bounds
        for(int l = -radius; l <= radius; l++)
          for(int m = -radius; m <= radius; m++)
          {
            const int ii = i + l;
            const int jj = j + m;
            const size_t idx_shift = ((ii * roi_out->width) + jj) * 4;
 
            const int ik = l + radius;
            const int jk = m + radius;
            const size_t idx_kernel = (ik * kernel_width) + jk;
            const float k = kernel[idx_kernel];
 
            for_four_channels(c, aligned(in : 64)) acc[c] += k * in[idx_shift + c];
          }
      }
      else
      {
        // We are close to borders, we need to clamp indices to bounds
        // assume constant boundary conditions
        for(int l = -radius; l <= radius; l++)
          for(int m = -radius; m <= radius; m++)
          {
            const int ii = CLAMP((int)i + l, (int)0, (int)roi_out->height - 1);
            const int jj = CLAMP((int)j + m, (int)0, (int)roi_out->width - 1);
            const size_t idx_shift = ((ii * roi_out->width) + jj) * 4;
 
            const int ik = l + radius;
            const int jk = m + radius;
            const size_t idx_kernel = (ik * kernel_width) + jk;
            const float k = kernel[idx_kernel];
 
            for_four_channels(c, aligned(in : 64)) acc[c] += k * in[idx_shift + c];
          }
      }
 
      for_each_channel(c, aligned(out : 64) aligned(acc : 16)) out[index + c] = acc[c];
 
      // copy alpha
      out[index + 3] = in[index + 3];
    }
  }
  dt_free_align(kernel);
  return dt_atomic_get_int(&abandoned);
}
 
 
#if HAVE_OPENCL
int process_cl(struct dt_iop_module_t *self, const dt_dev_pixelpipe_t *pipe, const dt_dev_pixelpipe_iop_t *piece, cl_mem dev_in, cl_mem dev_out)
{
  const dt_iop_roi_t *const roi_in = &piece->roi_in;
  const dt_iop_roi_t *const roi_out = &piece->roi_out;
  dt_iop_blurs_params_t *p = (dt_iop_blurs_params_t *)piece->data;
  dt_iop_blurs_global_data_t *const gd = (dt_iop_blurs_global_data_t *)self->global_data;
 
  cl_int err = -999;
 
  const int devid = pipe->devid;
  const int width = roi_in->width;
  const int height = roi_in->height;
 
  size_t sizes[] = { ROUNDUPDWD(width, devid), ROUNDUPDHT(height, devid), 1 };
 
  // Init the blur kernel
  const float scale = dt_dev_get_module_scale(pipe, roi_in);
  const int radius = MAX(roundf(p->radius / scale), 2);
  const size_t kernel_width = 2 * radius + 1;
 
  float *const restrict kernel = dt_alloc_align_float(kernel_width * kernel_width);
  if(IS_NULL_PTR(kernel)) return FALSE;
  if(build_pixel_kernel(kernel, kernel_width, kernel_width, p))
  {
    dt_free_align(kernel);
    return FALSE;
  }
 
  cl_mem kernel_cl = dt_opencl_copy_host_to_device(devid, kernel, kernel_width, kernel_width, sizeof(float));
 
  dt_opencl_set_kernel_arg(devid, gd->kernel_blurs_convolve, 0, sizeof(cl_mem), (void *)&dev_in);
  dt_opencl_set_kernel_arg(devid, gd->kernel_blurs_convolve, 1, sizeof(cl_mem), (void *)&kernel_cl);
  dt_opencl_set_kernel_arg(devid, gd->kernel_blurs_convolve, 2, sizeof(cl_mem), (void *)&dev_out);
  dt_opencl_set_kernel_arg(devid, gd->kernel_blurs_convolve, 3, sizeof(int), (void *)&roi_out->width);
  dt_opencl_set_kernel_arg(devid, gd->kernel_blurs_convolve, 4, sizeof(int), (void *)&roi_out->height);
  dt_opencl_set_kernel_arg(devid, gd->kernel_blurs_convolve, 5, sizeof(int), (void *)&radius);
 
  err = dt_opencl_enqueue_kernel_2d(devid, gd->kernel_blurs_convolve, sizes);
  if(err != CL_SUCCESS) goto error;
 
  // cleanup and exit on success
  dt_free_align(kernel);
  dt_opencl_release_mem_object(kernel_cl);
  return TRUE;
 
error:
  dt_free_align(kernel);
  dt_opencl_release_mem_object(kernel_cl);
  dt_print(DT_DEBUG_OPENCL, "[opencl_blurs] couldn't enqueue kernel! %d\n", err);
  return FALSE;
}
 
void init_global(dt_iop_module_so_t *module)
{
  const int program = 34;
  dt_iop_blurs_global_data_t *gd = (dt_iop_blurs_global_data_t *)malloc(sizeof(dt_iop_blurs_global_data_t));
  module->data = gd;
  gd->kernel_blurs_convolve = dt_opencl_create_kernel(program, "convolve");
}
 
 
void cleanup_global(dt_iop_module_so_t *module)
{
  dt_iop_blurs_global_data_t *gd = (dt_iop_blurs_global_data_t *)module->data;
  dt_opencl_free_kernel(gd->kernel_blurs_convolve);
  dt_free(module->data);
}
#endif
 
 
void gui_changed(dt_iop_module_t *self, GtkWidget *w, void *previous)
{
  dt_iop_blurs_params_t *p = (dt_iop_blurs_params_t *)self->params;
  dt_iop_blurs_gui_data_t *g = (dt_iop_blurs_gui_data_t *)dt_iop_gui_data(self);
 
  if(IS_NULL_PTR(w) || w == g->type)
  {
    if(p->type == DT_BLUR_LENS)
    {
      gtk_widget_hide(g->angle);
      gtk_widget_hide(g->curvature);
      gtk_widget_hide(g->offset);
 
      gtk_widget_show(g->blades);
      gtk_widget_show(g->concavity);
      gtk_widget_show(g->rotation);
      gtk_widget_show(g->linearity);
    }
    else if(p->type == DT_BLUR_MOTION)
    {
      gtk_widget_show(g->angle);
      gtk_widget_show(g->curvature);
      gtk_widget_show(g->offset);
 
      gtk_widget_hide(g->blades);
      gtk_widget_hide(g->concavity);
      gtk_widget_hide(g->rotation);
      gtk_widget_hide(g->linearity);
    }
    else if(p->type == DT_BLUR_GAUSSIAN)
    {
      gtk_widget_hide(g->angle);
      gtk_widget_hide(g->curvature);
      gtk_widget_hide(g->offset);
 
      gtk_widget_hide(g->blades);
      gtk_widget_hide(g->concavity);
      gtk_widget_hide(g->rotation);
      gtk_widget_hide(g->linearity);
    }
  }
 
  // update kernel view
  if(g->img_cached)
  {
    if(build_gui_kernel(g->img, g->img_width, g->img_width, p)) return;
    gtk_widget_queue_draw(GTK_WIDGET(g->area));
  }
}
 
static gboolean dt_iop_tonecurve_draw(GtkWidget *widget, cairo_t *crf, gpointer user_data)
{
  dt_iop_module_t *self = (dt_iop_module_t *)user_data;
  dt_iop_blurs_gui_data_t *g = (dt_iop_blurs_gui_data_t *)dt_iop_gui_data(self);
  dt_iop_blurs_params_t *p = (dt_iop_blurs_params_t *)self->params;
 
  GtkAllocation allocation;
  GtkStyleContext *context = gtk_widget_get_style_context(widget);
  gtk_widget_get_allocation(widget, &allocation);
  gtk_render_background(context, crf, 0, 0, allocation.width, allocation.height);
 
  if(allocation.width != g->img_width)
  {
    // Widget size changed, flush the cache buffer and restart
    g->img_cached = FALSE;
    dt_free_align(g->img);
    g->img = NULL;
  }
 
  if(!g->img_cached)
  {
    g->img = dt_alloc_align(sizeof(unsigned char) * 4 * allocation.width * allocation.width);
    if(IS_NULL_PTR(g->img)) return FALSE;
    g->img_width = allocation.width;
    g->img_cached = TRUE;
    if(build_gui_kernel(g->img, g->img_width, g->img_width, p)) 
      return FALSE;
 
    // Note: if params change, we silently recompute the img in the buffer
    // no need to flush the cache. Flush only if buffer size changes,
    // aka GUI widget gets resized.
  }
 
  // Paint the kernel
  const int stride = cairo_format_stride_for_width(CAIRO_FORMAT_ARGB32, g->img_width);
  cairo_surface_t *cst = cairo_image_surface_create_for_data(g->img, CAIRO_FORMAT_ARGB32,
                            g->img_width, g->img_width, stride);
 
  cairo_set_source_surface(crf, cst, 0, 0);
  cairo_paint(crf);
  cairo_surface_destroy(cst);
  return TRUE;
}
 
void gui_update(dt_iop_module_t *self)
{
// FIXME check why needed
  gui_changed(self, NULL, NULL);
}
 
#define DEG_TO_RAD 180.f / M_PI_F
 
void gui_init(dt_iop_module_t *self)
{
  dt_iop_blurs_gui_data_t *g = IOP_GUI_ALLOC(blurs);
 
  self->gui->widget = gtk_box_new(GTK_ORIENTATION_VERTICAL, DT_GUI_BOX_SPACING);
 
  // Image buffer to store the kernel look
  // Don't recompute it in the drawing function, only when a param is changed
  // then serve it from cache to the drawing function.
  g->img_cached = FALSE;
  g->img = NULL;
  g->img_width = 0.f;
 
  g->area = GTK_DRAWING_AREA(dtgtk_drawing_area_new_with_aspect_ratio(1.f));
  g_signal_connect(G_OBJECT(g->area), "draw", G_CALLBACK(dt_iop_tonecurve_draw), self);
  gtk_box_pack_start(GTK_BOX(self->gui->widget), GTK_WIDGET(g->area), TRUE, TRUE, 0);
 
  g->radius = dt_bauhaus_slider_from_params(self, "radius");
  dt_bauhaus_slider_set_format(g->radius, " px");
 
  g->type = dt_bauhaus_combobox_from_params(self, "type");
 
  g->blades = dt_bauhaus_slider_from_params(self, "blades");
  g->concavity = dt_bauhaus_slider_from_params(self, "concavity");
  g->linearity = dt_bauhaus_slider_from_params(self, "linearity");
  g->rotation = dt_bauhaus_slider_from_params(self, "rotation");
  dt_bauhaus_slider_set_factor(g->rotation, DEG_TO_RAD);
  dt_bauhaus_slider_set_format(g->rotation, "\302\260");
 
  g->angle = dt_bauhaus_slider_from_params(self, "angle");
  dt_bauhaus_slider_set_factor(g->angle, DEG_TO_RAD);
  dt_bauhaus_slider_set_format(g->angle, "\302\260");
 
 
  g->curvature = dt_bauhaus_slider_from_params(self, "curvature");
  g->offset = dt_bauhaus_slider_from_params(self, "offset");
 
}
 
void gui_cleanup(dt_iop_module_t *self)
{
  dt_iop_blurs_gui_data_t *g = (dt_iop_blurs_gui_data_t *)dt_iop_gui_data(self);
  dt_free_align(g->img);
  g->img = NULL;
  IOP_GUI_FREE;
}
 
 
// 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