/*
    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 "common/gdk_event_utils.h"
#include "bauhaus/bauhaus.h"
#include "common/darktable.h"
#include "common/imagebuf.h"
#include "common/dwt.h"
#include "control/control.h"
#include "develop/imageop.h"
#include "develop/imageop_math.h"
#include "develop/imageop_gui.h"
#include "develop/tiling.h"
#include "develop/openmp_maths.h"
#include "dtgtk/drawingarea.h"
#include "gui/accelerators.h"
#include "gui/gtk.h"
#include "iop/iop_api.h"
#include "common/opencl.h"

#include <gtk/gtk.h>
#include <stdlib.h>
#include <strings.h>

DT_MODULE_INTROSPECTION(2, dt_iop_rawdenoise_params_t)

#define DT_IOP_RAWDENOISE_INSET DT_PIXEL_APPLY_DPI(5)
#define DT_IOP_RAWDENOISE_RES 64
#define DT_IOP_RAWDENOISE_BANDS 5

typedef enum dt_iop_rawdenoise_channel_t
{
  DT_RAWDENOISE_ALL = 0,
  DT_RAWDENOISE_R = 1,
  DT_RAWDENOISE_G = 2,
  DT_RAWDENOISE_B = 3,
  DT_RAWDENOISE_NONE = 4
} dt_iop_rawdenoise_channel_t;

typedef struct dt_iop_rawdenoise_params_t
{
  float threshold; // $MIN: 0.0 $MAX: 1.0 $DEFAULT: 0.01 $DESCRIPTION: "noise threshold"
  float x[DT_RAWDENOISE_NONE][DT_IOP_RAWDENOISE_BANDS];
  float y[DT_RAWDENOISE_NONE][DT_IOP_RAWDENOISE_BANDS]; // $DEFAULT: 0.5
} dt_iop_rawdenoise_params_t;

typedef struct dt_iop_rawdenoise_gui_data_t
{
  dt_draw_curve_t *transition_curve; // curve for gui to draw

  GtkWidget *threshold;
  GtkDrawingArea *area;
  GtkNotebook *channel_tabs;
  double mouse_x, mouse_y, mouse_pick;
  float mouse_radius;
  dt_iop_rawdenoise_params_t drag_params;
  int dragging;
  int x_move;
  dt_iop_rawdenoise_channel_t channel;
  float draw_xs[DT_IOP_RAWDENOISE_RES], draw_ys[DT_IOP_RAWDENOISE_RES];
  float draw_min_xs[DT_IOP_RAWDENOISE_RES], draw_min_ys[DT_IOP_RAWDENOISE_RES];
  float draw_max_xs[DT_IOP_RAWDENOISE_RES], draw_max_ys[DT_IOP_RAWDENOISE_RES];
} dt_iop_rawdenoise_gui_data_t;

typedef struct dt_iop_rawdenoise_data_t
{
  float threshold;
  dt_draw_curve_t *curve[DT_RAWDENOISE_NONE];
  dt_iop_rawdenoise_channel_t channel;
  float force[DT_RAWDENOISE_NONE][DT_IOP_RAWDENOISE_BANDS];
} dt_iop_rawdenoise_data_t;

typedef struct dt_iop_rawdenoise_global_data_t
{
  int kernel_rawdenoise_extract;
  int kernel_rawdenoise_scatter;
  int kernel_rawdenoise_clear;
  int kernel_rawdenoise_dwt_vert;
  int kernel_rawdenoise_dwt_horiz;
} dt_iop_rawdenoise_global_data_t;

int legacy_params(dt_iop_module_t *self,
                  const void *const old_params,
                  const int old_version,
                  void **new_params,
                  int32_t *new_params_size,
                  int *new_version)
{
  typedef struct dt_iop_rawdenoise_params_v2_t
  {
    float threshold;
    float x[DT_RAWDENOISE_NONE][DT_IOP_RAWDENOISE_BANDS];
    float y[DT_RAWDENOISE_NONE][DT_IOP_RAWDENOISE_BANDS];
  } dt_iop_rawdenoise_params_v2_t;

  if(old_version == 1)
  {
    typedef struct dt_iop_rawdenoise_params_v1_t
    {
      float threshold;
    } dt_iop_rawdenoise_params_v1_t;

    const dt_iop_rawdenoise_params_v1_t *o = (dt_iop_rawdenoise_params_v1_t *)old_params;
    dt_iop_rawdenoise_params_v2_t *n = malloc(sizeof(dt_iop_rawdenoise_params_v2_t));

    n->threshold = o->threshold;
    for(int k = 0; k < DT_IOP_RAWDENOISE_BANDS; k++)
    {
      for(int ch = 0; ch < DT_RAWDENOISE_NONE; ch++)
      {
        n->x[ch][k] = k / (DT_IOP_RAWDENOISE_BANDS - 1.0);
        n->y[ch][k] = 0.5f;
      }
    }

    *new_params = n;
    *new_params_size = sizeof(dt_iop_rawdenoise_params_v2_t);
    *new_version = 2;
    return 0;
  }
  return 1;
}


const char *name()
{
  return _("raw denoise");
}

const char **description(dt_iop_module_t *self)
{
  return dt_iop_set_description(self, _("denoise the raw image early in the pipeline"),
                                      _("corrective"),
                                      _("linear, raw, scene-referred"),
                                      _("linear, raw"),
                                      _("linear, raw, scene-referred"));
}

int flags()
{
  return IOP_FLAGS_SUPPORTS_BLENDING;
}

int default_group()
{
  return IOP_GROUP_CORRECT | IOP_GROUP_TECHNICAL;
}

dt_iop_colorspace_type_t default_colorspace(dt_iop_module_t *self,
                                            dt_dev_pixelpipe_t *pipe,
                                            dt_dev_pixelpipe_iop_t *piece)
{
  return IOP_CS_RAW;
}

#define BIT16 65536.0

static void compute_channel_noise(float *const noise, int color, const dt_iop_rawdenoise_data_t *const data)
{
  // note that these constants are the same for X-Trans and Bayer, as they are proportional to image detail on
  // each channel, not the sensor pattern
  static const float noise_all[] = { 0.8002, 0.2735, 0.1202, 0.0585, 0.0291, 0.0152, 0.0080, 0.0044 };
  for(int i = 0; i < DT_IOP_RAWDENOISE_BANDS; i++)
  {
    // scale the value from [0,1] to [0,16],
    // and makes the "0.5" neutral value become 1
    float chan_threshold_exp_4;
    switch(color)
    {
    case 0:
      chan_threshold_exp_4 = data->force[DT_RAWDENOISE_R][DT_IOP_RAWDENOISE_BANDS - i - 1];
      break;
    case 2:
      chan_threshold_exp_4 = data->force[DT_RAWDENOISE_B][DT_IOP_RAWDENOISE_BANDS - i - 1];
      break;
    default:
      chan_threshold_exp_4 = data->force[DT_RAWDENOISE_G][DT_IOP_RAWDENOISE_BANDS - i - 1];
      break;
    }
    chan_threshold_exp_4 *= chan_threshold_exp_4;
    chan_threshold_exp_4 *= chan_threshold_exp_4;
    // repeat for the overall all-channels thresholds
    float all_threshold_exp_4 = data->force[DT_RAWDENOISE_ALL][DT_IOP_RAWDENOISE_BANDS - i - 1];
    all_threshold_exp_4 *= all_threshold_exp_4;
    all_threshold_exp_4 *= all_threshold_exp_4;
    noise[i] = noise_all[i] * all_threshold_exp_4 * chan_threshold_exp_4 * 16.0f * 16.0f;
    // the following multiplication needs to stay separate from the above line, because merging the two changes
    // the results on the integration test!
    noise[i] *= data->threshold;
  }
}

static void wavelet_denoise(const float *const restrict in, float *const restrict out, const dt_iop_roi_t *const roi,
                            const dt_iop_rawdenoise_data_t * const data, const uint32_t filters)
{
  const size_t size = (size_t)(roi->width / 2 + 1) * (roi->height / 2 + 1);
  float *const restrict fimg = dt_alloc_align_float(size);
  if(!fimg)
    return;

  const int nc = 4;
  for(int c = 0; c < nc; c++) /* denoise R,G1,B,G3 individually */
  {
    const int color = FC(c % 2, c / 2, filters);
    float noise[DT_IOP_RAWDENOISE_BANDS];
    compute_channel_noise(noise,color,data);

    // adjust for odd width and height
    const int halfwidth = roi->width / 2 + (roi->width & (~(c >> 1)) & 1);
    const int halfheight = roi->height / 2 + (roi->height & (~c) & 1);

    // collect one of the R/G1/G2/B channels into a monochrome image, applying sqrt() to the values as a
    // variance-stabilizing transform
    DT_OMP_FOR()
    for(int row = c & 1; row < roi->height; row += 2)
    {
      float *const restrict fimgp = fimg + (size_t)row / 2 * halfwidth;
      const int offset = (c & 2) >> 1;
      const float *const restrict inp = in + (size_t)row * roi->width + offset;
      const int senselwidth = (roi->width-offset+1)/2;
      for(int col = 0; col < senselwidth; col++)
        fimgp[col] = sqrtf(MAX(0.0f, inp[2*col]));
    }

    // perform the wavelet decomposition and denoising
    dwt_denoise(fimg,halfwidth,halfheight,DT_IOP_RAWDENOISE_BANDS,noise);

    // distribute the denoised data back out to the original R/G1/G2/B channel, squaring the resulting values to
    // undo the original transform
    DT_OMP_FOR()
    for(int row = c & 1; row < roi->height; row += 2)
    {
      const float *const restrict fimgp = fimg + (size_t)row / 2 * halfwidth;
      const int offset = (c & 2) >> 1;
      float *const restrict outp = out + (size_t)row * roi->width + offset;
      const int senselwidth = (roi->width-offset+1)/2;
      for(int col = 0; col < senselwidth; col++)
      {
        float d = fimgp[col];
        outp[2*col] = d * d;
      }
    }
  }
#if 0
  /* FIXME: Haven't ported this part yet */
  float maximum = 1.0;		/* FIXME */
  float black = 0.0;		/* FIXME */
  maximum *= BIT16;
  black *= BIT16;
  for(c=0; c<4; c++)
    cblack[c] *= BIT16;
  if(filters && colors == 3)	/* pull G1 and G3 closer together */
  {
    float *window[4];
    int wlast, blk[2];
    float mul[2];
    float thold = threshold/512;
    for(row=0; row < 2; row++)
    {
      mul[row] = 0.125 * pre_mul[FC(row+1,0) | 1] / pre_mul[FC(row,0) | 1];
      blk[row] = cblack[FC(row,0) | 1];
    }
    for(i=0; i < 4; i++)
      window[i] = fimg + width*i;
    for(wlast=-1, row=1; row < height-1; row++)
    {
      while(wlast < row+1)
      {
        for(wlast++, i=0; i < 4; i++)
          window[(i+3) & 3] = window[i];
        for(col = FC(wlast,1) & 1; col < width; col+=2)
          window[2][col] = BAYER(wlast,col);
      }
      for(col = (FC(row,0) & 1)+1; col < width-1; col+=2)
      {
        float avg = ( window[0][col-1] + window[0][col+1] +
                      window[2][col-1] + window[2][col+1] - blk[~row & 1]*4 )
                    * mul[row & 1] + (window[1][col] + blk[row & 1]) * 0.5;
        avg = avg > 0 ? sqrtf(avg) : 0;
        float diff = sqrtf(BAYER(row,col)) - avg;
        if      (diff < -thold) diff += thold;
        else if(diff >  thold) diff -= thold;
        else diff = 0;
        BAYER(row,col) = SQR(avg+diff);
      }
    }
  }
#endif
  dt_free_align(fimg);
}

static inline float vstransform(const float value)
{
  return sqrtf(fmaxf(0.0f, value));
}

static void wavelet_denoise_xtrans(const float *const restrict in,
                                   float *const restrict out,
                                   const dt_iop_roi_t *const restrict roi,
                                   const dt_iop_rawdenoise_data_t *const data,
                                   const uint8_t (*const xtrans)[6])
{
  const int width = roi->width;
  const int height = roi->height;
  const size_t size = (size_t)width * height;
  // regions smaller than a CFA period may lack matching sensels for the border interpolation
  if(width < 6 || height < 6)
  {
    memcpy(out, in, sizeof(float) * size);
    return;
  }
  // allocate a buffer for the particular color channel to be denoise; we add two rows to simplify the
  // channel-extraction code (no special case for top/bottom row)
  float *const img = dt_alloc_align_float((size_t)width * (height+2));
  if(!img)
  {
    // we ran out of memory, so just pass through the image without denoising
    memcpy(out, in, sizeof(float) * size);
    return;
  }
  float *const fimg = img + width;	// point at the actual color channel contents in the buffer

  for(int c = 0; c < 3; c++)
  {
    float noise[DT_IOP_RAWDENOISE_BANDS];
    compute_channel_noise(noise, c, data);

    // negative values mark border positions that no matching sensel fills; vstransform() is nonnegative
    // so these positions can be filled from the completed plane before wavelet filtering
    for(size_t col = 0; col < width; col++)
    {
      fimg[col] = -1.0f;
      fimg[(size_t)(height-1)*width + col] = -1.0f;
    }
    // build the plane in one pass, top to bottom: sensels are also copied into the neighboring rows, so
    // the order of the writes decides which value a position keeps; splitting the rows among threads
    // would make the result depend on the thread count
    for(size_t row = 0; row < height; row++)
    {
      const float *const restrict inp = in + row * width;
      float *const restrict fimgp = fimg + row * width;
      // handle red/blue pixel in first column
      if(c != 1 && FCNxtrans(row, 0, xtrans) == c)
      {
        // copy to neighbors above and right
        const float d = vstransform(inp[0]);
        fimgp[0] = fimgp[-width] = fimgp[-width+1] = d;
      }
      for(size_t col = (c != 1); col < width-1; col++)
      {
        if(FCNxtrans(row, col, xtrans) == c)
        {
          // the pixel at the current location has the desired color, so apply sqrt() as a variance-stablizing
          // transform, and then do cheap nearest-neighbor interpolation by copying it to appropriate neighbors
          const float d = vstransform(inp[col]);
          fimgp[col] = d;
          if(c == 1) // green pixel
          {
            // Copy to the right and down.  The X-Trans color layout is such that copying to those two neighbors
            // results in all positions being filled except in the left-most and right-most columns and sometimes
            // the topmost and bottom-most rows (depending on how the ROI aligns with the CFA).
            fimgp[col+1] = fimgp[col+width] = d;
          }
          else // red or blue pixel
          {
            // Copy value to all eight neighbors; it's OK to copy to the row above even when we're in row 0 (or
            // the row below when in the last row) because the destination is then one of the two extra rows of
            // the buffer, which dwt_denoise() never reads.  We need to copy to all adjacent positions because
            // there may be two green pixels between nearest red/red or blue/blue, so each will cover one of the
            // greens.
            fimgp[col-width-1] = fimgp[col-width] = fimgp[col-width+1] = d; // row above
            fimgp[col-1] = fimgp[col+1] = d;                                // left and right
            fimgp[col+width-1] = fimgp[col+width] = fimgp[col+width+1] = d; // row below
          }
        }
      }
      // leftmost and rightmost pixel in the row may still need to be filled in from a neighbor
      if(FCNxtrans(row, 0, xtrans) != c)
      {
        int src = 0;	// fallback is current sensel even if it has the wrong color
        if(row > 0 && FCNxtrans(row-1, 0, xtrans) == c)
          src = -width;
        else if(FCNxtrans(row, 1, xtrans) == c)
          src = 1;
        else if(row > 0 && FCNxtrans(row-1, 1, xtrans) == c)
          src = -width + 1;
        fimgp[0] = (src == 0 && (row == 0 || row == height-1)) ? -1.0f : vstransform(inp[src]);
      }
      // check the right-most pixel; if it's the desired color and not green, copy it to the neighbors
      if(c != 1 && FCNxtrans(row, width-1, xtrans) == c)
      {
        // copy to neighbors above and left
        const float d = vstransform(inp[width-1]);
        fimgp[width-2] = fimgp[width-1] = fimgp[-1] = d;
      }
      else if(FCNxtrans(row, width-1, xtrans) != c)
      {
        int src = width-1;	// fallback is current sensel even if it has the wrong color
        if(FCNxtrans(row, width-2, xtrans) == c)
          src = width-2;
        else if(row > 0 && FCNxtrans(row-1, width-1, xtrans) == c)
          src = -1;
        else if(row > 0 && FCNxtrans(row-1, width-2, xtrans) == c)
          src = -2;
        fimgp[width-1] = (src == width-1 && (row == 0 || row == height-1)) ? -1.0f : vstransform(inp[src]);
      }
      else // green right-most pixel in the green plane
      {
        // the loop above stops one column short, so it never writes this sensel into the plane. If the left
        // neighbor is green, the position would keep that neighbor's value. Otherwise it would keep the value
        // the red pass left there, or an unset position in the top and bottom rows
        fimgp[width-1] = vstransform(inp[width-1]);
      }
    }

    // missing neighbors beyond the border leave holes; extend the completed plane one row outward
    float *const last = fimg + (size_t)(height-1)*width;
    for(size_t col = 0; col < width; col++)
    {
      if(fimg[col] < 0.0f) fimg[col] = fimg[width + col];
      if(last[col] < 0.0f) last[col] = fimg[(size_t)(height-2)*width + col];
    }

    // perform the wavelet decomposition and denoising
    dwt_denoise(fimg,width,height,DT_IOP_RAWDENOISE_BANDS,noise);

    // distribute the denoised data back out to the original R/G/B channel, squaring the resulting values to
    // undo the original transform
    DT_OMP_FOR()
    for(int row = 0; row < height; row++)
    {
      const float *const restrict fimgp = fimg + (size_t)row * width;
      float *const restrict outp = out + (size_t)row * width;
      for(int col = 0; col < width; col++)
        if(FCNxtrans(row, col, xtrans) == c)
        {
          float d = fimgp[col];
          outp[col] = d * d;
        }
    }
  }

  dt_free_align(img);
}

void process(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 dt_iop_rawdenoise_data_t *const restrict d = piece->data;

  if(!(d->threshold > 0.0f))
  {
    dt_iop_image_copy_by_size(ovoid, ivoid, roi_in->width, roi_in->height, piece->colors);
  }
  else
  {
    const uint32_t filters = piece->filters;
    const uint8_t(*const xtrans)[6] = piece->xtrans;
    if(filters != 9u)
      wavelet_denoise(ivoid, ovoid, roi_in, d, filters);
    else
      wavelet_denoise_xtrans(ivoid, ovoid, roi_in, d, xtrans);
  }
}

#ifdef HAVE_OPENCL
void 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,
                     dt_develop_tiling_t *tiling)
{
  default_tiling_callback(self, piece, roi_in, roi_out, tiling);
  // process_cl() adds three quarter-size planes for the Bayer wavelets
  tiling->factor_cl += 0.75f;
}

int process_cl(dt_iop_module_t *self,
               dt_dev_pixelpipe_iop_t *piece,
               cl_mem dev_in,
               cl_mem dev_out,
               const dt_iop_roi_t *const roi_in,
               const dt_iop_roi_t *const roi_out)
{
  const dt_iop_rawdenoise_data_t *const d = piece->data;
  const dt_iop_rawdenoise_global_data_t *const gd = self->global_data;
  const int devid = piece->pipe->devid;
  const int width = roi_in->width;
  const int height = roi_in->height;

  if(!(d->threshold > 0.0f))
  {
    const size_t region[2] = { width, height };
    return dt_opencl_enqueue_copy_image(devid, dev_in, dev_out, CLIMG_ORIGIN, CLIMG_ORIGIN, region);
  }

  // the X-Trans planes are built sequentially, see wavelet_denoise_xtrans()
  const uint32_t filters = piece->filters;
  if(filters == 9u) return DT_OPENCL_PROCESS_CL;

  const int maxhalfwidth = width / 2 + (width & 1);
  const int maxhalfheight = height / 2 + (height & 1);
  const size_t size = (size_t)maxhalfwidth * maxhalfheight * sizeof(float);

  cl_int err = CL_MEM_OBJECT_ALLOCATION_FAILURE;
  cl_mem dev_fimg = dt_opencl_alloc_device_buffer(devid, size);
  cl_mem dev_interm = dt_opencl_alloc_device_buffer(devid, size);
  cl_mem dev_accum = dt_opencl_alloc_device_buffer(devid, size);
  if(!dev_fimg || !dev_interm || !dev_accum) goto error;

  for(int c = 0; c < 4; c++) /* denoise R,G1,B,G3 individually */
  {
    const int color = FC(c % 2, c / 2, filters);
    float noise[DT_IOP_RAWDENOISE_BANDS];
    compute_channel_noise(noise, color, d);

    // as in wavelet_denoise()
    const int halfwidth = width / 2 + (width & (~(c >> 1)) & 1);
    const int halfheight = height / 2 + (height & (~c) & 1);

    err = dt_opencl_enqueue_kernel_2d_args(devid, gd->kernel_rawdenoise_extract, halfwidth, halfheight,
            CLARG(dev_in), CLARG(dev_fimg), CLARG(c), CLARG(halfwidth), CLARG(halfheight));
    if(err != CL_SUCCESS) goto error;
    err = dt_opencl_enqueue_kernel_2d_args(devid, gd->kernel_rawdenoise_clear, halfwidth, halfheight,
            CLARG(dev_accum), CLARG(halfwidth), CLARG(halfheight));
    if(err != CL_SUCCESS) goto error;

    for(int lev = 0; lev < DT_IOP_RAWDENOISE_BANDS; lev++)
    {
      const int last = lev + 1 == DT_IOP_RAWDENOISE_BANDS;
      const int vscale = MIN(1 << lev, halfheight - 1);
      const int hscale = MIN(1 << lev, halfwidth - 1);
      err = dt_opencl_enqueue_kernel_2d_args(devid, gd->kernel_rawdenoise_dwt_vert, halfwidth, halfheight,
              CLARG(dev_fimg), CLARG(dev_interm), CLARG(halfwidth), CLARG(halfheight), CLARG(vscale));
      if(err != CL_SUCCESS) goto error;
      err = dt_opencl_enqueue_kernel_2d_args(devid, gd->kernel_rawdenoise_dwt_horiz, halfwidth, halfheight,
              CLARG(dev_interm), CLARG(dev_fimg), CLARG(dev_accum),
              CLARG(halfwidth), CLARG(halfheight), CLARG(hscale), CLARG(noise[lev]), CLARG(last));
      if(err != CL_SUCCESS) goto error;
    }

    err = dt_opencl_enqueue_kernel_2d_args(devid, gd->kernel_rawdenoise_scatter, halfwidth, halfheight,
            CLARG(dev_fimg), CLARG(dev_out), CLARG(c), CLARG(halfwidth), CLARG(halfheight));
    if(err != CL_SUCCESS) goto error;
  }

error:
  dt_opencl_release_mem_object(dev_fimg);
  dt_opencl_release_mem_object(dev_interm);
  dt_opencl_release_mem_object(dev_accum);
  return err;
}
#endif

void init_global(dt_iop_module_so_t *self)
{
  const int program = 20; // dwt.cl, from programs.conf
  dt_iop_rawdenoise_global_data_t *gd = malloc(sizeof(dt_iop_rawdenoise_global_data_t));
  self->data = gd;
  gd->kernel_rawdenoise_extract = dt_opencl_create_kernel(program, "rawdenoise_extract");
  gd->kernel_rawdenoise_scatter = dt_opencl_create_kernel(program, "rawdenoise_scatter");
  gd->kernel_rawdenoise_clear = dt_opencl_create_kernel(program, "rawdenoise_clear");
  gd->kernel_rawdenoise_dwt_vert = dt_opencl_create_kernel(program, "rawdenoise_dwt_vert");
  gd->kernel_rawdenoise_dwt_horiz = dt_opencl_create_kernel(program, "rawdenoise_dwt_horiz");
}

void cleanup_global(dt_iop_module_so_t *self)
{
  dt_iop_rawdenoise_global_data_t *gd = self->data;
  if(gd)
  {
    dt_opencl_free_kernel(gd->kernel_rawdenoise_extract);
    dt_opencl_free_kernel(gd->kernel_rawdenoise_scatter);
    dt_opencl_free_kernel(gd->kernel_rawdenoise_clear);
    dt_opencl_free_kernel(gd->kernel_rawdenoise_dwt_vert);
    dt_opencl_free_kernel(gd->kernel_rawdenoise_dwt_horiz);
  }
  free(self->data);
  self->data = NULL;
}

void init(dt_iop_module_t *self)
{
  dt_iop_default_init(self);

  dt_iop_rawdenoise_params_t *d = self->default_params;

  for(int k = 0; k < DT_IOP_RAWDENOISE_BANDS; k++)
  {
    for(int ch = 0; ch < DT_RAWDENOISE_NONE; ch++)
    {
      d->x[ch][k] = k / (DT_IOP_RAWDENOISE_BANDS - 1.f);
    }
  }
}

void reload_defaults(dt_iop_module_t *self)
{
  // can't be switched on for non-raw images:
  self->hide_enable_button = !dt_image_is_raw(&self->dev->image_storage);

  if(self->widget)
  {
    gtk_stack_set_visible_child_name(GTK_STACK(self->widget), self->hide_enable_button ? "non_raw" : "raw");
  }

  self->default_enabled = FALSE;
}

void commit_params(dt_iop_module_t *self, dt_iop_params_t *params, dt_dev_pixelpipe_t *pipe,
                   dt_dev_pixelpipe_iop_t *piece)
{
  dt_iop_rawdenoise_params_t *p = (dt_iop_rawdenoise_params_t *)params;
  dt_iop_rawdenoise_data_t *d = piece->data;

  d->threshold = p->threshold;

  for(int ch = 0; ch < DT_RAWDENOISE_NONE; ch++)
  {
    dt_draw_curve_set_point(d->curve[ch], 0, p->x[ch][DT_IOP_RAWDENOISE_BANDS - 2] - 1.0, p->y[ch][0]);
    for(int k = 0; k < DT_IOP_RAWDENOISE_BANDS; k++)
      dt_draw_curve_set_point(d->curve[ch], k, p->x[ch][k], p->y[ch][k]);
    dt_draw_curve_set_point(d->curve[ch], DT_IOP_RAWDENOISE_BANDS + 1, p->x[ch][1] + 1.0,
                            p->y[ch][DT_IOP_RAWDENOISE_BANDS - 1]);
    dt_draw_curve_calc_values(d->curve[ch], 0.0, 1.0, DT_IOP_RAWDENOISE_BANDS, NULL, d->force[ch]);
  }

  if(!(dt_image_is_raw(&pipe->image)))
    piece->enabled = FALSE;
}

void init_pipe(dt_iop_module_t *self, dt_dev_pixelpipe_t *pipe, dt_dev_pixelpipe_iop_t *piece)
{
  dt_iop_rawdenoise_data_t *d = malloc(sizeof(dt_iop_rawdenoise_data_t));
  const dt_iop_rawdenoise_params_t *const default_params = self->default_params;

  piece->data = (void *)d;
  for(int ch = 0; ch < DT_RAWDENOISE_NONE; ch++)
  {
    d->curve[ch] = dt_draw_curve_new(0.0, 1.0, CATMULL_ROM);
    for(int k = 0; k < DT_IOP_RAWDENOISE_BANDS; k++)
      (void)dt_draw_curve_add_point(d->curve[ch], default_params->x[ch][k], default_params->y[ch][k]);
  }
}

void cleanup_pipe(dt_iop_module_t *self, dt_dev_pixelpipe_t *pipe, dt_dev_pixelpipe_iop_t *piece)
{
  dt_iop_rawdenoise_data_t *d = piece->data;
  for(int ch = 0; ch < DT_RAWDENOISE_NONE; ch++) dt_draw_curve_destroy(d->curve[ch]);
  free(piece->data);
  piece->data = NULL;
}

void gui_update(dt_iop_module_t *self)
{
  dt_iop_rawdenoise_gui_data_t *g = self->gui_data;
  gtk_widget_queue_draw(GTK_WIDGET(g->area));
}

static void dt_iop_rawdenoise_get_params(dt_iop_rawdenoise_params_t *p, const int ch, const double mouse_x,
                                         const double mouse_y, const float rad)
{
  for(int k = 0; k < DT_IOP_RAWDENOISE_BANDS; k++)
  {
    const float f = expf(-(mouse_x - p->x[ch][k]) * (mouse_x - p->x[ch][k]) / (rad * rad));
    p->y[ch][k] = (1 - f) * p->y[ch][k] + f * mouse_y;
  }
}

static gboolean rawdenoise_draw(GtkWidget *widget, cairo_t *crf, dt_iop_module_t *self)
{
  dt_iop_rawdenoise_gui_data_t *g = self->gui_data;
  dt_iop_rawdenoise_params_t p = *(dt_iop_rawdenoise_params_t *)self->params;

  int ch = (int)g->channel;
  dt_draw_curve_set_point(g->transition_curve, 0, p.x[ch][DT_IOP_RAWDENOISE_BANDS - 2] - 1.0, p.y[ch][0]);
  for(int k = 0; k < DT_IOP_RAWDENOISE_BANDS; k++)
    dt_draw_curve_set_point(g->transition_curve, k + 1, p.x[ch][k], p.y[ch][k]);
  dt_draw_curve_set_point(g->transition_curve, DT_IOP_RAWDENOISE_BANDS + 1, p.x[ch][1] + 1.0,
                          p.y[ch][DT_IOP_RAWDENOISE_BANDS - 1]);

  const int inset = DT_IOP_RAWDENOISE_INSET;
  GtkAllocation allocation;
  gtk_widget_get_allocation(widget, &allocation);
  int width = allocation.width, height = allocation.height;
  cairo_surface_t *cst = dt_cairo_image_surface_create(CAIRO_FORMAT_ARGB32, width, height);
  cairo_t *cr = cairo_create(cst);
  cairo_set_source_rgb(cr, .2, .2, .2);

  cairo_paint(cr);

  cairo_translate(cr, inset, inset);
  width -= 2 * inset;
  height -= 2 * inset;

  cairo_set_line_width(cr, DT_PIXEL_APPLY_DPI(1.0));
  cairo_set_source_rgb(cr, .1, .1, .1);
  cairo_rectangle(cr, 0, 0, width, height);
  cairo_stroke(cr);

  cairo_set_source_rgb(cr, .3, .3, .3);
  cairo_rectangle(cr, 0, 0, width, height);
  cairo_fill(cr);

  // draw grid
  cairo_set_line_width(cr, DT_PIXEL_APPLY_DPI(.4));
  cairo_set_source_rgb(cr, .1, .1, .1);
  dt_draw_grid(cr, 8, 0, 0, width, height);

  if(g->mouse_y > 0 || g->dragging)
  {
    // draw min/max curves:
    dt_iop_rawdenoise_get_params(&p, g->channel, g->mouse_x, 1., g->mouse_radius);
    dt_draw_curve_set_point(g->transition_curve, 0, p.x[ch][DT_IOP_RAWDENOISE_BANDS - 2] - 1.0, p.y[ch][0]);
    for(int k = 0; k < DT_IOP_RAWDENOISE_BANDS; k++)
      dt_draw_curve_set_point(g->transition_curve, k + 1, p.x[ch][k], p.y[ch][k]);
    dt_draw_curve_set_point(g->transition_curve, DT_IOP_RAWDENOISE_BANDS + 1, p.x[ch][1] + 1.0,
                            p.y[ch][DT_IOP_RAWDENOISE_BANDS - 1]);
    dt_draw_curve_calc_values(g->transition_curve, 0.0, 1.0, DT_IOP_RAWDENOISE_RES, g->draw_min_xs, g->draw_min_ys);

    p = *(dt_iop_rawdenoise_params_t *)self->params;
    dt_iop_rawdenoise_get_params(&p, g->channel, g->mouse_x, .0, g->mouse_radius);
    dt_draw_curve_set_point(g->transition_curve, 0, p.x[ch][DT_IOP_RAWDENOISE_BANDS - 2] - 1.0, p.y[ch][0]);
    for(int k = 0; k < DT_IOP_RAWDENOISE_BANDS; k++)
      dt_draw_curve_set_point(g->transition_curve, k + 1, p.x[ch][k], p.y[ch][k]);
    dt_draw_curve_set_point(g->transition_curve, DT_IOP_RAWDENOISE_BANDS + 1, p.x[ch][1] + 1.0,
                            p.y[ch][DT_IOP_RAWDENOISE_BANDS - 1]);
    dt_draw_curve_calc_values(g->transition_curve, 0.0, 1.0, DT_IOP_RAWDENOISE_RES, g->draw_max_xs, g->draw_max_ys);
  }

  cairo_save(cr);

  // draw selected cursor
  cairo_translate(cr, 0, height);

  cairo_set_operator(cr, CAIRO_OPERATOR_OVER);
  cairo_set_line_width(cr, DT_PIXEL_APPLY_DPI(2.));

  for(int i = 0; i < DT_RAWDENOISE_NONE; i++)
  {
    // draw curves, selected last
    ch = ((int)g->channel + i + 1) % DT_RAWDENOISE_NONE;
    float alpha = 0.3;
    if(i == DT_RAWDENOISE_NONE - 1) alpha = 1.0;
    switch(ch)
    {
      case 0:
        cairo_set_source_rgba(cr, .7, .7, .7, alpha);
        break;
      case 1:
        cairo_set_source_rgba(cr, .7, .1, .1, alpha);
        break;
      case 2:
        cairo_set_source_rgba(cr, .1, .7, .1, alpha);
        break;
      case 3:
        cairo_set_source_rgba(cr, .1, .1, .7, alpha);
        break;
    }

    p = *(dt_iop_rawdenoise_params_t *)self->params;
    dt_draw_curve_set_point(g->transition_curve, 0, p.x[ch][DT_IOP_RAWDENOISE_BANDS - 2] - 1.0, p.y[ch][0]);
    for(int k = 0; k < DT_IOP_RAWDENOISE_BANDS; k++)
      dt_draw_curve_set_point(g->transition_curve, k + 1, p.x[ch][k], p.y[ch][k]);
    dt_draw_curve_set_point(g->transition_curve, DT_IOP_RAWDENOISE_BANDS + 1, p.x[ch][1] + 1.0,
                            p.y[ch][DT_IOP_RAWDENOISE_BANDS - 1]);
    dt_draw_curve_calc_values(g->transition_curve, 0.0, 1.0, DT_IOP_RAWDENOISE_RES, g->draw_xs, g->draw_ys);
    cairo_move_to(cr, 0 * width / (float)(DT_IOP_RAWDENOISE_RES - 1), -height * g->draw_ys[0]);
    for(int k = 1; k < DT_IOP_RAWDENOISE_RES; k++)
      cairo_line_to(cr, k * width / (float)(DT_IOP_RAWDENOISE_RES - 1), -height * g->draw_ys[k]);
    cairo_stroke(cr);
  }

  ch = g->channel;
  // draw dots on knots
  cairo_set_source_rgb(cr, 0.7, 0.7, 0.7);
  cairo_set_line_width(cr, DT_PIXEL_APPLY_DPI(1.));
  for(int k = 0; k < DT_IOP_RAWDENOISE_BANDS; k++)
  {
    cairo_arc(cr, width * p.x[ch][k], -height * p.y[ch][k], DT_PIXEL_APPLY_DPI(3.0), 0.0, 2.0 * M_PI);
    if(g->x_move == k)
      cairo_fill(cr);
    else
      cairo_stroke(cr);
  }

  if(g->mouse_y > 0 || g->dragging)
  {
    // draw min/max, if selected
    cairo_set_source_rgba(cr, .7, .7, .7, .6);
    cairo_move_to(cr, 0, -height * g->draw_min_ys[0]);
    for(int k = 1; k < DT_IOP_RAWDENOISE_RES; k++)
      cairo_line_to(cr, k * width / (float)(DT_IOP_RAWDENOISE_RES - 1), -height * g->draw_min_ys[k]);
    for(int k = DT_IOP_RAWDENOISE_RES - 1; k >= 0; k--)
      cairo_line_to(cr, k * width / (float)(DT_IOP_RAWDENOISE_RES - 1), -height * g->draw_max_ys[k]);
    cairo_close_path(cr);
    cairo_fill(cr);
    // draw mouse focus circle
    cairo_set_source_rgba(cr, .9, .9, .9, .5);
    const float pos = DT_IOP_RAWDENOISE_RES * g->mouse_x;
    int k = (int)pos;
    const float f = k - pos;
    if(k >= DT_IOP_RAWDENOISE_RES - 1) k = DT_IOP_RAWDENOISE_RES - 2;
    float ht = -height * (f * g->draw_ys[k] + (1 - f) * g->draw_ys[k + 1]);
    cairo_arc(cr, g->mouse_x * width, ht, g->mouse_radius * width, 0, 2. * M_PI);
    cairo_stroke(cr);
  }

  cairo_restore(cr);

  cairo_set_operator(cr, CAIRO_OPERATOR_SOURCE);

  // draw labels:
  PangoLayout *layout;
  PangoRectangle ink;
  PangoFontDescription *desc = dt_gui_get_font();
  pango_font_description_set_weight(desc, PANGO_WEIGHT_BOLD);
  pango_font_description_set_absolute_size(desc, (.08 * height) * PANGO_SCALE);
  layout = pango_cairo_create_layout(cr);
  pango_layout_set_font_description(layout, desc);
  cairo_set_source_rgb(cr, .1, .1, .1);

  pango_layout_set_text(layout, C_("graph", "coarse"), -1);
  pango_layout_get_pixel_extents(layout, &ink, NULL);
  cairo_move_to(cr, .02 * width - ink.y, .5 * (height + ink.width));
  cairo_save(cr);
  cairo_rotate(cr, -M_PI * .5f);
  pango_cairo_show_layout(cr, layout);
  cairo_restore(cr);

  pango_layout_set_text(layout, _("fine"), -1);
  pango_layout_get_pixel_extents(layout, &ink, NULL);
  cairo_move_to(cr, .98 * width - ink.height, .5 * (height + ink.width));
  cairo_save(cr);
  cairo_rotate(cr, -M_PI * .5f);
  pango_cairo_show_layout(cr, layout);
  cairo_restore(cr);


  pango_layout_set_text(layout, _("smooth"), -1);
  pango_layout_get_pixel_extents(layout, &ink, NULL);
  cairo_move_to(cr, .5 * (width - ink.width), .08 * height - ink.height);
  pango_cairo_show_layout(cr, layout);

  pango_layout_set_text(layout, _("noisy"), -1);
  pango_layout_get_pixel_extents(layout, &ink, NULL);
  cairo_move_to(cr, .5 * (width - ink.width), .97 * height - ink.height);
  pango_cairo_show_layout(cr, layout);

  pango_font_description_free(desc);
  g_object_unref(layout);
  cairo_destroy(cr);
  cairo_set_source_surface(crf, cst, 0, 0);
  cairo_paint(crf);
  cairo_surface_destroy(cst);
  return FALSE;
}

static void rawdenoise_motion_notify(GtkEventControllerMotion *controller, gdouble x, gdouble y, dt_iop_module_t *self)
{
  dt_iop_rawdenoise_gui_data_t *g = self->gui_data;
  dt_iop_rawdenoise_params_t *p = self->params;
  GtkWidget *widget = dt_gui_get_widget(controller);
  const int inset = DT_IOP_RAWDENOISE_INSET;
  GtkAllocation allocation;
  gtk_widget_get_allocation(widget, &allocation);
  int height = allocation.height - 2 * inset, width = allocation.width - 2 * inset;
  if(!g->dragging) g->mouse_x = CLAMP(x - inset, 0, width) / (float)width;
  g->mouse_y = 1.0 - CLAMP(y - inset, 0, height) / (float)height;
  if(g->dragging)
  {
    *p = g->drag_params;
    if(g->x_move < 0)
    {
      dt_iop_rawdenoise_get_params(p, g->channel, g->mouse_x, g->mouse_y + g->mouse_pick, g->mouse_radius);
    }
    gtk_widget_queue_draw(widget);
    dt_dev_add_history_item_target(darktable.develop, self, TRUE, widget + g->channel);
  }
  else
  {
    g->x_move = -1;
    gtk_widget_queue_draw(widget);
  }
}

static void rawdenoise_button_press(GtkGestureSingle *gesture, gint n_press, gdouble x, gdouble y, dt_iop_module_t *self)
{
  dt_iop_rawdenoise_gui_data_t *g = self->gui_data;
  const int ch = g->channel;
  if(gtk_gesture_single_get_current_button(gesture) != GDK_BUTTON_PRIMARY)
    return;

  GtkWidget *widget = dt_gui_get_widget(gesture);

  if(n_press >= 2)
  {
    // reset current curve
    dt_iop_rawdenoise_params_t *p = self->params;
    const dt_iop_rawdenoise_params_t *const d = self->default_params;
    for(int k = 0; k < DT_IOP_RAWDENOISE_BANDS; k++)
    {
      p->x[ch][k] = d->x[ch][k];
      p->y[ch][k] = d->y[ch][k];
    }
    dt_dev_add_history_item_target(darktable.develop, self, TRUE, widget + ch);
    gtk_widget_queue_draw(GTK_WIDGET(g->area));
  }
  else
  {
    g->drag_params = *(dt_iop_rawdenoise_params_t *)self->params;
    const int inset = DT_IOP_RAWDENOISE_INSET;
    GtkAllocation allocation;
    gtk_widget_get_allocation(widget, &allocation);
    int height = allocation.height - 2 * inset, width = allocation.width - 2 * inset;
    g->mouse_pick
        = dt_draw_curve_calc_value(g->transition_curve, CLAMP(x - inset, 0, width) / (float)width);
    g->mouse_pick -= 1.0 - CLAMP(y - inset, 0, height) / (float)height;
    g->dragging = 1;
  }
}

static void rawdenoise_button_release(GtkGestureSingle *gesture, gint n_press, gdouble x, gdouble y, dt_iop_module_t *self)
{
  if(gtk_gesture_single_get_current_button(gesture) != GDK_BUTTON_PRIMARY)
    return;

  dt_iop_rawdenoise_gui_data_t *g = self->gui_data;
  g->dragging = 0;
}

static void rawdenoise_leave_notify(GtkEventControllerMotion *controller, dt_iop_module_t *self)
{
  dt_iop_rawdenoise_gui_data_t *g = self->gui_data;
  if(!g->dragging) g->mouse_y = -1.0;
  gtk_widget_queue_draw(dt_gui_get_widget(controller));
}

static void rawdenoise_scrolled(GtkEventControllerScroll *controller, gdouble dx, gdouble dy, dt_iop_module_t *self)
{
  dt_iop_rawdenoise_gui_data_t *g = self->gui_data;

  if(dt_modifier_eq(controller, GDK_MOD1_MASK))
  {
    // alt+scroll switches the channel tab (as before the controller conversion)
    dt_gui_forward_scroll(controller, GTK_WIDGET(g->channel_tabs));
    return;
  }

  if(dy != 0.0)
  {
    g->mouse_radius = CLAMP(g->mouse_radius * (1.0 - 0.1 * dy), 0.2 / DT_IOP_RAWDENOISE_BANDS, 1.0);
    gtk_widget_queue_draw(dt_gui_get_widget(controller));
  }
}

static void rawdenoise_tab_switch(GtkNotebook *notebook, GtkWidget *page, guint page_num, dt_iop_module_t *self)
{
  DT_GUARD_GUI_UPDATE();
  dt_iop_rawdenoise_gui_data_t *g = self->gui_data;
  g->channel = (dt_iop_rawdenoise_channel_t)page_num;
  gtk_widget_queue_draw(GTK_WIDGET(g->area));
}

void gui_init(dt_iop_module_t *self)
{
  dt_iop_rawdenoise_gui_data_t *g = IOP_GUI_ALLOC(rawdenoise);
  const dt_iop_rawdenoise_params_t *const p = self->default_params;

  g->channel = dt_conf_get_int("plugins/darkroom/rawdenoise/gui_channel");
  g->channel_tabs = GTK_NOTEBOOK(gtk_notebook_new());
  dt_action_define_iop(self, NULL, N_("channel"), GTK_WIDGET(g->channel_tabs), &dt_action_def_tabs_all_rgb);

  dt_ui_notebook_page(g->channel_tabs, N_("all"), NULL);
  dt_ui_notebook_page(g->channel_tabs, N_("R"), NULL);
  dt_ui_notebook_page(g->channel_tabs, N_("G"), NULL);
  dt_ui_notebook_page(g->channel_tabs, N_("B"), NULL);

  gtk_widget_show(gtk_notebook_get_nth_page(g->channel_tabs, g->channel));
  gtk_notebook_set_current_page(g->channel_tabs, g->channel);
  g_signal_connect(G_OBJECT(g->channel_tabs), "switch_page", G_CALLBACK(rawdenoise_tab_switch), self);

  const int ch = (int)g->channel;
  g->transition_curve = dt_draw_curve_new(0.0, 1.0, CATMULL_ROM);
  (void)dt_draw_curve_add_point(g->transition_curve, p->x[ch][DT_IOP_RAWDENOISE_BANDS - 2] - 1.0,
                                p->y[ch][DT_IOP_RAWDENOISE_BANDS - 2]);
  for(int k = 0; k < DT_IOP_RAWDENOISE_BANDS; k++)
    (void)dt_draw_curve_add_point(g->transition_curve, p->x[ch][k], p->y[ch][k]);
  (void)dt_draw_curve_add_point(g->transition_curve, p->x[ch][1] + 1.0, p->y[ch][1]);

  g->mouse_x = g->mouse_y = g->mouse_pick = -1.0;
  g->dragging = 0;
  g->x_move = -1;
  g->mouse_radius = 1.0 / (DT_IOP_RAWDENOISE_BANDS * 2);

  g->area = GTK_DRAWING_AREA(dt_ui_resize_wrap(NULL,
                                               0,
                                               "plugins/darkroom/rawdenoise/graphheight"));
  g_object_set_data(G_OBJECT(g->area), "iop-instance", self);
  dt_action_define_iop(self, NULL, N_("graph"), GTK_WIDGET(g->area), NULL);

  GtkWidget *box_raw = self->widget = dt_gui_vbox(g->channel_tabs, g->area);

  g_signal_connect(G_OBJECT(g->area), "draw", G_CALLBACK(rawdenoise_draw), self);
  dt_gui_connect_click(g->area, rawdenoise_button_press, rawdenoise_button_release, self);
  dt_gui_connect_motion(g->area, rawdenoise_motion_notify, NULL, rawdenoise_leave_notify, self);
  dt_gui_connect_scroll(g->area, GTK_EVENT_CONTROLLER_SCROLL_BOTH_AXES
                                   | GTK_EVENT_CONTROLLER_SCROLL_DISCRETE,
                        rawdenoise_scrolled, self);

  g->threshold = dt_bauhaus_slider_from_params(self, "threshold");
  dt_bauhaus_slider_set_soft_max(g->threshold, 0.1);
  dt_bauhaus_slider_set_digits(g->threshold, 3);

  // start building top level widget
  self->widget = gtk_stack_new();
  gtk_stack_set_homogeneous(GTK_STACK(self->widget), FALSE);

  GtkWidget *label_non_raw = dt_ui_label_new(_("raw denoising\nonly works for raw images."));

  gtk_stack_add_named(GTK_STACK(self->widget), label_non_raw, "non_raw");
  gtk_stack_add_named(GTK_STACK(self->widget), box_raw, "raw");
}

void gui_cleanup(dt_iop_module_t *self)
{
  dt_iop_rawdenoise_gui_data_t *g = self->gui_data;
  dt_conf_set_int("plugins/darkroom/rawdenoise/gui_channel", g->channel);
  dt_draw_curve_destroy(g->transition_curve);
}
// 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