liquify.c

Go to the documentation of this file.

/*
    This file is part of darktable,
    Copyright (C) 2014 Marcello.
    Copyright (C) 2016 Chris Hodapp.
    Copyright (C) 2016, 2018 johannes hanika.
    Copyright (C) 2016-2018, 2020-2022 Pascal Obry.
    Copyright (C) 2016 Roman Lebedev.
    Copyright (C) 2016, 2018-2019 Tobias Ellinghaus.
    Copyright (C) 2016 Ulrich Pegelow.
    Copyright (C) 2017, 2019-2020 Heiko Bauke.
    Copyright (C) 2017, 2019 luzpaz.
    Copyright (C) 2017-2018 Matthieu Moy.
    Copyright (C) 2018-2021, 2023-2026 Aurélien PIERRE.
    Copyright (C) 2018-2019 Edgardo Hoszowski.
    Copyright (C) 2018 Maurizio Paglia.
    Copyright (C) 2018 rawfiner.
    Copyright (C) 2019 Andreas Schneider.
    Copyright (C) 2019 mepi0011.
    Copyright (C) 2020-2021 Aldric Renaudin.
    Copyright (C) 2020-2021 Diederik Ter Rahe.
    Copyright (C) 2020 GrahamByrnes.
    Copyright (C) 2020-2022 Hanno Schwalm.
    Copyright (C) 2020-2021 Hubert Kowalski.
    Copyright (C) 2020 Marco.
    Copyright (C) 2020 Martin Straeten.
    Copyright (C) 2020-2021 Ralf Brown.
    Copyright (C) 2022 Martin Bařinka.
    Copyright (C) 2022 Philipp Lutz.
    Copyright (C) 2024 Alynx Zhou.
    Copyright (C) 2025 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 "caches/pixelpipe_cache_alloc.h"
#include "widgets/gdkkeys.h"
#include "widgets/accelerators.h"
#include "config.h"
#endif
#include "system/macros.h"
#include "system/openmp.h"
#include "system/target_clones.h"
#include "system/mem_alloc.h"
#include "common/logging.h"
#include "common/module_versioning.h"
#include "pixel/interpolation.h"
#include "common/opencl.h"
#include "math/math.h"
#include "common/collection.h"
#include "common/conf.h"
#include "control/control.h"
#include "develop/geometry/geometry.h"
#include "develop/imageop.h"
#include "develop/imageop_gui.h"
#include "develop/develop.h"
 
#include "gui/application.h"
#include "iop/iop_api.h"
#include <assert.h>
#include <cairo.h>
#include <complex.h>
#include <math.h>
#include <stdlib.h>
#include "widgets/label.h"
#include "widgets/widget_settings.h"
 
#ifdef _OPENMP
#include <omp.h>
#endif
 
// this is the version of the modules parameters, and includes version information about compile-time dt
DT_MODULE_INTROSPECTION(1, dt_iop_liquify_params_t)
 
#pragma GCC diagnostic ignored "-Wshadow"
 
#define MAX_NODES 100 // max of nodes in one instance
 
const int   LOOKUP_OVERSAMPLE = 10;
const int   INTERPOLATION_POINTS = 100; // when interpolating bezier
const float STAMP_RELOCATION = 0.1;     // how many radii to move stamp forward when following a path
 
#define CONF_RADIUS "plugins/darkroom/liquify/radius"
#define CONF_ANGLE "plugins/darkroom/liquify/angle"
#define CONF_STRENGTH "plugins/darkroom/liquify/strength"
 
// enum of layers. sorted back to front.
 
typedef enum
{
  DT_LIQUIFY_LAYER_BACKGROUND,
  DT_LIQUIFY_LAYER_RADIUS,
  DT_LIQUIFY_LAYER_HARDNESS1,
  DT_LIQUIFY_LAYER_HARDNESS2,
  DT_LIQUIFY_LAYER_WARPS,
  DT_LIQUIFY_LAYER_PATH,
  DT_LIQUIFY_LAYER_CTRLPOINT1_HANDLE,
  DT_LIQUIFY_LAYER_CTRLPOINT2_HANDLE,
  DT_LIQUIFY_LAYER_RADIUSPOINT_HANDLE,
  DT_LIQUIFY_LAYER_HARDNESSPOINT1_HANDLE,
  DT_LIQUIFY_LAYER_HARDNESSPOINT2_HANDLE,
  DT_LIQUIFY_LAYER_STRENGTHPOINT_HANDLE,
  DT_LIQUIFY_LAYER_CENTERPOINT,
  DT_LIQUIFY_LAYER_CTRLPOINT1,
  DT_LIQUIFY_LAYER_CTRLPOINT2,
  DT_LIQUIFY_LAYER_RADIUSPOINT,
  DT_LIQUIFY_LAYER_HARDNESSPOINT1,
  DT_LIQUIFY_LAYER_HARDNESSPOINT2,
  DT_LIQUIFY_LAYER_STRENGTHPOINT,
  DT_LIQUIFY_LAYER_LAST
} dt_liquify_layer_enum_t;
 
typedef enum
{
  DT_LIQUIFY_LAYER_FLAG_HIT_TEST      =  1,   
  DT_LIQUIFY_LAYER_FLAG_PREV_SELECTED =  2,   
  DT_LIQUIFY_LAYER_FLAG_NODE_SELECTED =  4,   
  DT_LIQUIFY_LAYER_FLAG_POINT_TOOL    =  8,   
  DT_LIQUIFY_LAYER_FLAG_LINE_TOOL     = 16,   
  DT_LIQUIFY_LAYER_FLAG_CURVE_TOOL    = 32,   
  DT_LIQUIFY_LAYER_FLAG_NODE_TOOL     = 64,   
  DT_LIQUIFY_LAYER_FLAG_ANY_TOOL      = 8 + 16 + 32 + 64,
} dt_liquify_layer_flag_enum_t;
 
typedef struct
{
  float red, green, blue, alpha;
} dt_liquify_rgba_t;
 
#define COLOR_NULL                 { 0.0, 0.0, 0.0, 0.8 }
#define GREY                       { 0.3, 0.3, 0.3, 0.8 }
#define LGREY                      { 0.8, 0.8, 0.8, 1.0 }
#define COLOR_DEBUG                { 0.9, 0.9, 0.0, 1.0 }
static const dt_liquify_rgba_t DT_LIQUIFY_COLOR_SELECTED = { 1.0, 1.0, 1.0, 1.0 };
static const dt_liquify_rgba_t DT_LIQUIFY_COLOR_HOVER    = { 1.0, 1.0, 1.0, 0.8 };
 
typedef struct
{
  dt_liquify_layer_enum_t hover_master;    
  dt_liquify_rgba_t fg;                    
  dt_liquify_rgba_t bg;                    
  float opacity;                           
  dt_liquify_layer_flag_enum_t flags;      
  const char *hint;                        
} dt_liquify_layer_t;
 
dt_liquify_layer_t dt_liquify_layers[] =
{
  { DT_LIQUIFY_LAYER_BACKGROUND,     COLOR_NULL,  COLOR_NULL, 0.0,  0,                                                                                                      },
  { DT_LIQUIFY_LAYER_RADIUS,         COLOR_DEBUG, COLOR_NULL, 0.25, DT_LIQUIFY_LAYER_FLAG_ANY_TOOL,                                                                         },
  { DT_LIQUIFY_LAYER_HARDNESS1,      COLOR_DEBUG, COLOR_NULL, 1.0,  0,                                                                                                      },
  { DT_LIQUIFY_LAYER_HARDNESS2,      COLOR_DEBUG, COLOR_NULL, 1.0,  0,                                                                                                      },
  { DT_LIQUIFY_LAYER_WARPS,          COLOR_DEBUG, LGREY,      0.5,  DT_LIQUIFY_LAYER_FLAG_ANY_TOOL,                                                                         },
  { DT_LIQUIFY_LAYER_PATH,           GREY,        LGREY,      1.0,  DT_LIQUIFY_LAYER_FLAG_ANY_TOOL  | DT_LIQUIFY_LAYER_FLAG_HIT_TEST,                                       },
  { DT_LIQUIFY_LAYER_CTRLPOINT1,     GREY,        LGREY,      1.0,  DT_LIQUIFY_LAYER_FLAG_NODE_TOOL,                                                                        },
  { DT_LIQUIFY_LAYER_CTRLPOINT2,     GREY,        LGREY,      1.0,  DT_LIQUIFY_LAYER_FLAG_NODE_TOOL,                                                                        },
  { DT_LIQUIFY_LAYER_RADIUSPOINT,    GREY,        LGREY,      1.0,  DT_LIQUIFY_LAYER_FLAG_NODE_TOOL,                                                                        },
  { DT_LIQUIFY_LAYER_HARDNESSPOINT1, GREY,        LGREY,      1.0,  DT_LIQUIFY_LAYER_FLAG_NODE_TOOL | DT_LIQUIFY_LAYER_FLAG_NODE_SELECTED,                                  },
  { DT_LIQUIFY_LAYER_HARDNESSPOINT2, GREY,        LGREY,      1.0,  DT_LIQUIFY_LAYER_FLAG_NODE_TOOL | DT_LIQUIFY_LAYER_FLAG_NODE_SELECTED,                                  },
  { DT_LIQUIFY_LAYER_STRENGTHPOINT,  GREY,        LGREY,      1.0,  DT_LIQUIFY_LAYER_FLAG_ANY_TOOL,                                                                         },
  { DT_LIQUIFY_LAYER_CENTERPOINT,    GREY,        LGREY,      1.0,  DT_LIQUIFY_LAYER_FLAG_ANY_TOOL  | DT_LIQUIFY_LAYER_FLAG_HIT_TEST,                                       },
  { DT_LIQUIFY_LAYER_CTRLPOINT1,     GREY,        LGREY,      1.0,  DT_LIQUIFY_LAYER_FLAG_NODE_TOOL | DT_LIQUIFY_LAYER_FLAG_HIT_TEST,                                       },
  { DT_LIQUIFY_LAYER_CTRLPOINT2,     GREY,        LGREY,      1.0,  DT_LIQUIFY_LAYER_FLAG_NODE_TOOL | DT_LIQUIFY_LAYER_FLAG_HIT_TEST,                                       },
  { DT_LIQUIFY_LAYER_RADIUSPOINT,    GREY,        LGREY,      1.0,  DT_LIQUIFY_LAYER_FLAG_NODE_TOOL | DT_LIQUIFY_LAYER_FLAG_HIT_TEST,                                       },
  { DT_LIQUIFY_LAYER_HARDNESSPOINT1, GREY,        LGREY,      1.0,  DT_LIQUIFY_LAYER_FLAG_NODE_TOOL | DT_LIQUIFY_LAYER_FLAG_NODE_SELECTED | DT_LIQUIFY_LAYER_FLAG_HIT_TEST, },
  { DT_LIQUIFY_LAYER_HARDNESSPOINT2, GREY,        LGREY,      1.0,  DT_LIQUIFY_LAYER_FLAG_NODE_TOOL | DT_LIQUIFY_LAYER_FLAG_NODE_SELECTED | DT_LIQUIFY_LAYER_FLAG_HIT_TEST, },
  { DT_LIQUIFY_LAYER_STRENGTHPOINT,  GREY,        LGREY,      1.0,  DT_LIQUIFY_LAYER_FLAG_ANY_TOOL  | DT_LIQUIFY_LAYER_FLAG_HIT_TEST,                                       }
};
 
typedef enum
{
  DT_LIQUIFY_UI_WIDTH_THINLINE,
  DT_LIQUIFY_UI_WIDTH_THICKLINE,
  DT_LIQUIFY_UI_WIDTH_DOUBLELINE,
  DT_LIQUIFY_UI_WIDTH_GIZMO,
  DT_LIQUIFY_UI_WIDTH_GIZMO_SMALL,
  DT_LIQUIFY_UI_WIDTH_DEFAULT_RADIUS,
  DT_LIQUIFY_UI_WIDTH_DEFAULT_STRENGTH,
  DT_LIQUIFY_UI_WIDTH_MIN_DRAG,
  DT_LIQUIFY_UI_WIDTH_LAST
} dt_liquify_ui_width_enum_t;
 
float dt_liquify_ui_widths [] =
{
  // value in 1/96 inch (that is: in pixels on a standard 96 dpi screen)
    2.0, // DT_LIQUIFY_UI_WIDTH_THINLINE
    3.0, // DT_LIQUIFY_UI_WIDTH_THICKLINE
    3.0, // DT_LIQUIFY_UI_WIDTH_DOUBLELINE
    9.0, // DT_LIQUIFY_UI_WIDTH_GIZMO
    7.0, // DT_LIQUIFY_UI_WIDTH_GIZMO_SMALL
  100.0, // DT_LIQUIFY_UI_WIDTH_DEFAULT_RADIUS,
   50.0, // DT_LIQUIFY_UI_WIDTH_DEFAULT_STRENGTH,
    3.0  // DT_LIQUIFY_UI_WIDTH_MIN_DRAG
};
 
typedef enum
{
  DT_LIQUIFY_WARP_TYPE_LINEAR,        // $DESCRIPTION: "linear" A linear warp originating from one point.
  DT_LIQUIFY_WARP_TYPE_RADIAL_GROW,   // $DESCRIPTION: "radial grow" A radial warp originating from one point.
  DT_LIQUIFY_WARP_TYPE_RADIAL_SHRINK, // $DESCRIPTION: "radial shrink"
  DT_LIQUIFY_WARP_TYPE_LAST
} dt_liquify_warp_type_enum_t;
 
typedef enum
{
  DT_LIQUIFY_NODE_TYPE_CUSP,        // $DESCRIPTION: "cusp"
  DT_LIQUIFY_NODE_TYPE_SMOOTH,      // $DESCRIPTION: "smooth"
  DT_LIQUIFY_NODE_TYPE_SYMMETRICAL, // $DESCRIPTION: "symmetrical"
  DT_LIQUIFY_NODE_TYPE_AUTOSMOOTH,  // $DESCRIPTION: "autosmooth"
  DT_LIQUIFY_NODE_TYPE_LAST
} dt_liquify_node_type_enum_t;
 
typedef enum
{
  DT_LIQUIFY_STATUS_NONE = 0,         // $DESCRIPTION: "none"
  DT_LIQUIFY_STATUS_NEW = 1,          // $DESCRIPTION: "new"
  DT_LIQUIFY_STATUS_INTERPOLATED = 2, // $DESCRIPTION: "interpolated"
  DT_LIQUIFY_STATUS_PREVIEW = 4,      // $DESCRIPTION: "preview"
  DT_LIQUIFY_STATUS_LAST
} dt_liquify_status_enum_t;
 
// enumerates the shapes types we use.
 
typedef enum
{
  DT_LIQUIFY_PATH_INVALIDATED = 0, // $DESCRIPTION: "invalidated"
  DT_LIQUIFY_PATH_MOVE_TO_V1,      // $DESCRIPTION: "move"
  DT_LIQUIFY_PATH_LINE_TO_V1,      // $DESCRIPTION: "line"
  DT_LIQUIFY_PATH_CURVE_TO_V1,     // $DESCRIPTION: "curve"
} dt_liquify_path_data_enum_t;
 
typedef struct
{
  dt_liquify_path_data_enum_t type;
  dt_liquify_node_type_enum_t node_type;
  dt_liquify_layer_enum_t selected;
  dt_liquify_layer_enum_t hovered;
  int8_t prev;
  int8_t idx;
  int8_t next;
} dt_liquify_path_header_t;
 
// Scalars and vectors are represented here as points because the only
// thing we can reasonably distort_transform are points.
 
typedef struct
{
  float complex point;
  float complex strength;   
  float complex radius;     
  float control1;           
  float control2;           
  dt_liquify_warp_type_enum_t type;
  dt_liquify_status_enum_t status;
} dt_liquify_warp_t;
 
typedef struct
{
  float complex ctrl1;
  float complex ctrl2;
} dt_liquify_node_t;
 
// set up lots of alternative ways to get at the popular members.
 
typedef struct
{
  dt_liquify_path_header_t header;
  dt_liquify_warp_t        warp;
  dt_liquify_node_t        node; // extended node data
} dt_liquify_path_data_t;
 
typedef struct
{
  dt_liquify_layer_enum_t layer;
  dt_liquify_path_data_t *elem;
} dt_liquify_hit_t;
 
static const dt_liquify_hit_t NOWHERE = { DT_LIQUIFY_LAYER_BACKGROUND, NULL };
 
typedef struct
{
  dt_liquify_path_data_t nodes[MAX_NODES];
} dt_iop_liquify_params_t;
 
typedef struct
{
  int warp_kernel;
} dt_iop_liquify_global_data_t;
 
typedef struct
{
  dt_iop_liquify_params_t params;
  int node_index; // last node index inserted
 
  float complex last_mouse_pos;
  float complex last_button1_pressed_pos;
  GdkModifierType last_mouse_mods;  
 
  dt_liquify_hit_t last_hit;      
  dt_liquify_hit_t dragging;      
 
  dt_liquify_path_data_t *temp;    
  dt_liquify_status_enum_t status; 
 
  GtkLabel *label_warp, *label_node;
  GtkToggleButton *btn_point_tool, *btn_line_tool, *btn_curve_tool, *btn_node_tool;
 
  gboolean just_started;
} dt_iop_liquify_gui_data_t;
 
 
// this returns a translatable name
const char *name()
{
  return _("l_iquify");
}
 
const char **description(struct dt_iop_module_t *self)
{
  return dt_iop_set_description(self, _("distort parts of the image"),
                                      _("creative"),
                                      _("linear, RGB, scene-referred"),
                                      _("geometric, RGB"),
                                      _("linear, RGB, scene-referred"));
}
 
 
int default_group()
{
  return IOP_GROUP_EFFECTS;
}
 
int flags()
{
  return IOP_FLAGS_SUPPORTS_BLENDING;
}
 
int operation_tags()
{
   return IOP_TAG_DISTORT;
}
 
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 input_format(dt_iop_module_t *self, dt_dev_pixelpipe_t *pipe, dt_dev_pixelpipe_iop_t *piece,
                  dt_iop_buffer_dsc_t *dsc)
{
  default_input_format(self, pipe, piece, dsc);
  dsc->channels = 4;
  dsc->datatype = TYPE_FLOAT;
}
 
/******************************************************************************/
/* Code common to op-engine and gui.                                          */
/******************************************************************************/
 
static inline float get_rot(const dt_liquify_warp_type_enum_t warp_type)
{
  if(warp_type == DT_LIQUIFY_WARP_TYPE_RADIAL_SHRINK)
    return DT_M_PI_F;
  else
    return 0.0f;
}
 
static dt_liquify_path_data_t *node_alloc(dt_iop_liquify_params_t *p, int *node_index)
{
  for(int k=0; k<MAX_NODES; k++)
    if(p->nodes[k].header.type == DT_LIQUIFY_PATH_INVALIDATED)
    {
      *node_index = k;
      p->nodes[k].header.idx = k;
      p->nodes[k].header.next = p->nodes[k].header.prev = -1;
      p->nodes[k].header.selected = p->nodes[k].header.hovered = 0;
      return &p->nodes[k];
    }
  return NULL;
}
 
static dt_liquify_path_data_t *node_prev(dt_iop_liquify_params_t *p, const dt_liquify_path_data_t *n)
{
  if(n->header.prev == -1)
    return NULL;
  else
    return &p->nodes[n->header.prev];
}
 
static dt_liquify_path_data_t *node_get(dt_iop_liquify_params_t *p, const int index)
{
  if(index > -1 && index < MAX_NODES)
    return &p->nodes[index];
  else
    return NULL;
}
 
static dt_liquify_path_data_t *node_next(dt_iop_liquify_params_t *p, const dt_liquify_path_data_t *n)
{
  if(n->header.next == -1)
    return NULL;
  else
    return &p->nodes[n->header.next];
}
 
static void node_insert_before(dt_iop_liquify_params_t *p, dt_liquify_path_data_t *this, dt_liquify_path_data_t *new)
{
  new->header.next  = this->header.idx;
  new->header.prev  = this->header.prev;
  if(this->header.prev != -1)
    p->nodes[this->header.prev].header.next = new->header.idx;
  this->header.prev = new->header.idx;
}
 
static void node_gc(dt_iop_liquify_params_t *p)
{
  int last=0;
  for(last=MAX_NODES-1; last>0; last--)
    if(p->nodes[last].header.type != DT_LIQUIFY_PATH_INVALIDATED)
      break;
  int k = 0;
 
  while(k<=last)
  {
    if(p->nodes[k].header.type == DT_LIQUIFY_PATH_INVALIDATED)
    {
      for(int e=0; e<last; e++)
      {
        //  then move slot if above position k
        if(e >= k)                       p->nodes[e] = p->nodes[e+1];
        //  update all pointers above position k
        if(e >= k)                       p->nodes[e].header.idx--;
        if(p->nodes[e].header.prev >= k) p->nodes[e].header.prev--;
        if(p->nodes[e].header.next >= k) p->nodes[e].header.next--;
      }
      last--;
    }
    else
      k++;
  }
  //  invalidate all nodes beyond the last moved one
  for(int k=last+1; k<MAX_NODES; k++)
    p->nodes[k].header.type = DT_LIQUIFY_PATH_INVALIDATED;
}
 
static void node_delete(dt_iop_liquify_params_t *p, dt_liquify_path_data_t *this)
{
  dt_liquify_path_data_t *prev = node_prev(p, this);
  dt_liquify_path_data_t *next = node_next(p, this);
 
  if(IS_NULL_PTR(prev) && !IS_NULL_PTR(next))
  {
    next->header.prev = -1;
    next->header.type = DT_LIQUIFY_PATH_MOVE_TO_V1;
  }
  else if(!IS_NULL_PTR(prev))
  {
    prev->header.next = this->header.next;
 
    if(!IS_NULL_PTR(next))
      next->header.prev = prev->header.idx;
  }
 
  this->header.prev = this->header.next = - 1;
  this->header.type = DT_LIQUIFY_PATH_INVALIDATED;
  node_gc(p);
}
 
static void path_delete(dt_iop_liquify_params_t *p, dt_liquify_path_data_t *this)
{
  dt_liquify_path_data_t *n = this;
 
  // clear next
  while(n)
  {
    n->header.type = DT_LIQUIFY_PATH_INVALIDATED;
    n = node_next(p, n);
  }
 
  // clear prev
  n = this;
  while(n)
  {
    n->header.type = DT_LIQUIFY_PATH_INVALIDATED;
    n = node_prev(p, n);
  }
  node_gc(p);
}
 
typedef struct
{
  dt_iop_liquify_params_t params;
  dt_iop_module_t *self;
} dt_iop_liquify_geometry_t;
 
typedef struct
{
  dt_develop_t *develop;
  const dt_dev_pixelpipe_t *pipe;
  /* When set, compose through the geometry service instead of a pixel pipe. Exactly one of
   * `pipe' and `chain' is used: the pipe when a piece is being processed, the chain when the
   * GUI asks the pixel-less service (develop/geometry/geometry.h). The paths are the same
   * fold in both cases -- see _distort_paths(). */
  dt_geometry_chain_t *chain;
  float from_scale;
  float to_scale;
  int transf_direction;
  gboolean from_distort_transform;
} distort_params_t;
 
static void _distort_paths(const struct dt_iop_module_t *module,
                            const distort_params_t *params, const dt_iop_liquify_params_t *p)
{
  int len = 0;
 
  // count nodes
 
  for(int k = 0; k < MAX_NODES; k++)
  {
    dt_liquify_path_data_t *data = (dt_liquify_path_data_t *) &p->nodes[k];
    if(data->header.type == DT_LIQUIFY_PATH_INVALIDATED)
      break;
 
    switch (data->header.type)
    {
    case DT_LIQUIFY_PATH_CURVE_TO_V1:
      len += 2;
      // fall thru
    case DT_LIQUIFY_PATH_MOVE_TO_V1:
    case DT_LIQUIFY_PATH_LINE_TO_V1:
      len += 3;
      break;
    default:
      break;
    }
  }
 
  // create buffer with all points
 
  float *buffer = malloc(sizeof(float) * 2 * len);
  float *b = buffer;
 
  for(int k = 0; k < MAX_NODES; k++)
  {
    dt_liquify_path_data_t *data = (dt_liquify_path_data_t *) &p->nodes[k];
    if(data->header.type == DT_LIQUIFY_PATH_INVALIDATED)
      break;
 
    switch (data->header.type)
    {
    case DT_LIQUIFY_PATH_CURVE_TO_V1:
      *b++ = crealf(data->node.ctrl1) / params->from_scale;
      *b++ = cimagf(data->node.ctrl1) / params->from_scale;
      *b++ = crealf(data->node.ctrl2) / params->from_scale;
      *b++ = cimagf(data->node.ctrl2) / params->from_scale;
      // fall thru
    case DT_LIQUIFY_PATH_MOVE_TO_V1:
    case DT_LIQUIFY_PATH_LINE_TO_V1:
      *b++ = crealf(data->warp.point) / params->from_scale;
      *b++ = cimagf(data->warp.point) / params->from_scale;
      *b++ = crealf(data->warp.strength) / params->from_scale;
      *b++ = cimagf(data->warp.strength) / params->from_scale;
      *b++ = crealf(data->warp.radius) / params->from_scale;
      *b++ = cimagf(data->warp.radius) / params->from_scale;
      break;
    default:
      break;
    }
  }
  if(!IS_NULL_PTR(params->chain))
  {
    /* The geometry service's equivalent of the two branches below. The bound is the caller's
     * own iop_order, exclusive, so this module is not re-entered and the recursion ends. */
    if(params->transf_direction == DT_DEV_TRANSFORM_DIR_ALL)
    {
      dt_geometry_chain_compose(params->chain, module->iop_order, DT_DEV_TRANSFORM_DIR_BACK_EXCL, buffer,
                                len);
      dt_geometry_chain_compose(params->chain, module->iop_order, DT_DEV_TRANSFORM_DIR_FORW_EXCL, buffer,
                                len);
    }
    else
      dt_geometry_chain_compose(params->chain, module->iop_order, params->transf_direction, buffer, len);
  }
  else if(params->from_distort_transform)
  {
    if(params->transf_direction == DT_DEV_TRANSFORM_DIR_ALL)
    {
      dt_dev_distort_transform_locked(params->pipe, module->iop_order,
                                      DT_DEV_TRANSFORM_DIR_BACK_EXCL, buffer, len);
      dt_dev_distort_transform_locked(params->pipe, module->iop_order,
                                      DT_DEV_TRANSFORM_DIR_FORW_EXCL, buffer, len);
    }
    else
      dt_dev_distort_transform_locked(params->pipe, module->iop_order,
                                      params->transf_direction, buffer, len);
  }
  else
  {
    if(params->transf_direction == DT_DEV_TRANSFORM_DIR_ALL)
    {
      dt_dev_distort_transform_plus(params->pipe, module->iop_order, DT_DEV_TRANSFORM_DIR_BACK_EXCL, buffer, len);
      dt_dev_distort_transform_plus(params->pipe, module->iop_order, DT_DEV_TRANSFORM_DIR_FORW_EXCL, buffer, len);
    }
    else
      dt_dev_distort_transform_plus(params->pipe, module->iop_order, params->transf_direction, buffer, len);
  }
 
  // record back the transformed points
 
  b = buffer;
 
  for(int k = 0; k < MAX_NODES; k++)
  {
    dt_liquify_path_data_t *data = (dt_liquify_path_data_t *) &p->nodes[k];
    if(data->header.type == DT_LIQUIFY_PATH_INVALIDATED)
      break;
 
    switch (data->header.type)
    {
       case DT_LIQUIFY_PATH_CURVE_TO_V1:
         data->node.ctrl1 = (b[0] + b[1] * I) * params->to_scale;
         b += 2;
         data->node.ctrl2 = (b[0] + b[1] * I) * params->to_scale;
         b += 2;
         // fall thru
       case DT_LIQUIFY_PATH_MOVE_TO_V1:
       case DT_LIQUIFY_PATH_LINE_TO_V1:
         data->warp.point = (b[0] + b[1] * I) * params->to_scale;
         b += 2;
         data->warp.strength = (b[0] + b[1] * I) * params->to_scale;
         b += 2;
         data->warp.radius = (b[0] + b[1] * I) * params->to_scale;
         b += 2;
         break;
       default:
         break;
    }
  }
 
  dt_free(buffer);
}
 
static void distort_paths_raw_to_piece(const struct dt_iop_module_t *module,
                                        const dt_dev_pixelpipe_t *pipe,
                                        const float roi_in_scale,
                                        dt_iop_liquify_params_t *p,
                                        const gboolean from_distort_transform)
{
  const distort_params_t params = { module->dev, pipe, NULL, 1.f, roi_in_scale, DT_DEV_TRANSFORM_DIR_BACK_EXCL, from_distort_transform };
  _distort_paths(module, &params, p);
}
 
static void distort_paths_raw_to_piece_chain(const struct dt_iop_module_t *module,
                                             dt_geometry_chain_t *chain, const float roi_in_scale,
                                             dt_iop_liquify_params_t *p)
{
  const distort_params_t params = { module->dev, NULL, chain, 1.f, roi_in_scale,
                                    DT_DEV_TRANSFORM_DIR_BACK_EXCL, TRUE };
  _distort_paths(module, &params, p);
}
 
// op-engine code
 
static inline float complex normalize(const float complex v)
{
  if(cabsf(v) < 0.000001f)
    return 1.0f;
  return v / cabsf(v);
}
 
// calculate the linear blend of scalars a and b
 
static inline float mix(const float a, const float b, const float t)
{
  return a + (b - a) * t;
}
 
 
// calculate the linear blend of points p0 and p1
 
static inline float complex cmix(const float complex p0, const float complex p1, const float t)
{
  return p0 + (p1 - p0) * t;
}
 
static void mix_warps(dt_liquify_warp_t *result,
                       const dt_liquify_warp_t *warp1,
                       const dt_liquify_warp_t *warp2,
                       const complex float pt,
                       const float t)
{
  result->type     = warp1->type;
  result->control1 = mix (warp1->control1, warp2->control1, t);
  result->control2 = mix (warp1->control2, warp2->control2, t);
 
  const float radius = mix(cabsf(warp1->radius - warp1->point), cabsf(warp2->radius - warp2->point), t);
  result->radius     = pt + radius;
 
  const complex float p1 = warp1->strength - warp1->point;
  const complex float p2 = warp2->strength - warp2->point;
  float arg1 = cargf(p1);
  float arg2 = cargf(p2);
  gboolean invert = FALSE;
 
  if(arg1 > .0f && arg2 < -(M_PI_F / 2.f))
  {
    invert = TRUE;
    arg1 = M_PI_F - arg1;
    arg2 = -M_PI_F - arg2;
  }
  else if(arg1 < -(M_PI_F / 2.f) && arg2 > .0f)
  {
    invert = TRUE;
    arg1 = -M_PI_F - arg1;
    arg2 = M_PI_F - arg2;
  }
 
  const float r    = mix(cabsf(p1), cabsf(p2), t);
  const float phi  = invert ? M_PI_F - mix(arg1, arg2, t) : mix(arg1, arg2, t);
 
  result->strength = pt + r * cexpf(phi * I);
  result->point    = pt;
}
 
// Interpolate a cubic bezier spline into a series of points.
 
static void interpolate_cubic_bezier(const float complex p0,
                                      const float complex p1,
                                      const float complex p2,
                                      const float complex p3,
                                      float complex buffer[],
                                      const int n)
{
  // convert from bernstein basis to polynomial basis to get faster math
  // See: http://www.tinaja.com/glib/cubemath.pdf
  const float complex A = p3 - 3 * p2 + 3 * p1 -     p0;
  const float complex B =      3 * p2 - 6 * p1 + 3 * p0;
  const float complex C =               3 * p1 - 3 * p0;
  const float complex D =                            p0;
 
  float complex *buf = buffer;
  const float step = 1.0f / n;
  float t = step;
  *buf++ = p0;
 
  for(int i = 1; i < n - 1; ++i)
  {
    *buf++ = ((A * t + B) * t + C) * t + D;
    t += step;
  }
  *buf = p3;
}
 
static GList *interpolate_paths(dt_iop_liquify_params_t *p);
 
/*
  Get approx. arc length of a curve.
 
  Used to approximate the arc length of a bezier curve.
*/
 
static float get_arc_length(const float complex points[], const int n_points)
{
  float length = 0.0f;
  for(int i = 1; i < n_points; i++)
    length += cabsf(points[i-1] - points[i]);
  return length;
}
 
typedef struct
{
  int i;
  float length;
} restart_cookie_t;
 
/*
  Interpolate a point on a curve at a specified arc length.
 
  In a bezier curve the parameter t usually does not correspond to
  the arc length.
*/
 
static float complex point_at_arc_length(const float complex points[], const int n_points,
                                          const float arc_length, restart_cookie_t *restart)
{
  float length = restart ? restart->length : 0.0f;
 
  for(int i = restart ? restart->i : 1; i < n_points; i++)
  {
    const float prev_length = length;
    length += cabsf(points[i-1] - points[i]);
 
    if(length >= arc_length)
    {
      const float t = (arc_length - prev_length) / (length - prev_length);
      if(!IS_NULL_PTR(restart))
      {
        restart->i = i;
        restart->length = prev_length;
      }
      return cmix(points[i - 1], points[i], t);
    }
  }
 
  return points[n_points - 1];
}
 
/*
  Build a lookup table for the warp intensity.
 
  Lookup table for the warp intensity function: f(x). The warp
  intensity function determines how much a pixel is influenced by the
  warp depending from its distance from a central point.
 
  Boundary conditions: f(0) must be 1 and f(@a distance) must be 0.
  f'(0) and f'(@a distance) must both be 0 or we'll get artifacts on
  the picture.
 
  Implementation: a bezier curve with p0 = 0, 1 and p3 = 1, 0. p1 is
  defined by @a control1, 1 and p2 by @a control1, 0.  Because a
  bezier is parameterized on t, we have to reparameterize on x, which
  we do by linear interpolation.
 
  Octave code:
 
  t = linspace(0,1,100);
  grid;
  hold on;
  for steps = 0:0.1:1
    cpoints = [0,1; steps,1; steps,0; 1,0];
    bezier = cbezier2poly(cpoints);
    x = polyval(bezier(1,:), t);
    y = polyval(bezier(2,:), t);
    plot(t, interp1(x, y, t));
  end
  hold off;
*/
 
static float *build_lookup_table(const int distance, const float control1, const float control2)
{
  float complex *clookup = dt_pixelpipe_cache_alloc_align_cache(sizeof(float complex) * (distance + 2), 0);
  if(IS_NULL_PTR(clookup)) return NULL;
 
  interpolate_cubic_bezier(I, control1 + I, control2, 1.0, clookup, distance + 2);
 
  // reparameterize bezier by x and keep only y values
  float *lookup = dt_pixelpipe_cache_alloc_align_float_cache((size_t)(distance + 2), 0);
  if(IS_NULL_PTR(lookup))
  {
    dt_pixelpipe_cache_free_align(clookup);
    return NULL;
  }
  float *ptr = lookup;
  float complex *cptr = clookup + 1;
  const float complex *cptr_end = cptr + distance;
  const float step = 1.0f / (float) distance;
  float x = 0.0f;
 
  *ptr++ = 1.0f;
  for(int i = 1; i < distance && cptr < cptr_end; i++)
  {
    x += step;
    while(crealf(*cptr) < x && cptr < cptr_end)
      cptr++;
    const float dx1 = crealf(cptr[0] - cptr[-1]);
    const float dx2 = x - crealf(cptr[-1]);
    *ptr++ = cimagf(cptr[0]) +(dx2 / dx1) * (cimagf(cptr[0]) - cimagf(cptr[-1]));
  }
  *ptr++ = 0.0f;
 
  dt_pixelpipe_cache_free_align(clookup);
  return lookup;
}
 
static void compute_round_stamp_extent(cairo_rectangle_int_t *const restrict stamp_extent,
                                        const dt_liquify_warp_t *const restrict warp)
{
 
  const int iradius = round(cabsf(warp->radius - warp->point));
  assert(iradius > 0);
 
  stamp_extent->x = stamp_extent->y = -iradius;
  stamp_extent->x += crealf(warp->point);
  stamp_extent->y += cimagf(warp->point);
  stamp_extent->width = stamp_extent->height = 2 * iradius + 1;
}
 
/*
  Compute a round(circular) stamp.
 
  The stamp is a vector field of warp vectors around a center point.
 
  In a linear warp the center point gets a warp of @a strength, while
  points on the circumference of the circle get no warp at all.
  Between center and circumference the warp magnitude tapers off
  following a curve (see: build_lookup_table()).
 
  Note that when applying a linear stamp to a path, we will first rotate its
  vectors into the direction of the path.
 
  In a radial warp the center point and the points on the
  circumference get no warp. Between center and circumference the
  warp magnitude follows a curve with maximum at radius / 0.5
 
  Our stamp is stored in a rectangular region.
*/
 
static int build_round_stamp(float complex **pstamp,
                              cairo_rectangle_int_t *const restrict stamp_extent,
                              const dt_liquify_warp_t *const restrict warp)
{
  const int iradius = round(cabsf(warp->radius - warp->point));
  assert(iradius > 0);
 
  stamp_extent->x = stamp_extent->y = -iradius;
  stamp_extent->width = stamp_extent->height = 2 * iradius + 1;
 
  // 0.5 is factored in so the warp starts to degenerate when the
  // strength arrow crosses the warp radius.
  float complex strength = 0.5f * (warp->strength - warp->point);
  strength = (warp->status & DT_LIQUIFY_STATUS_INTERPOLATED) ?
    (strength * STAMP_RELOCATION) : strength;
  const float abs_strength = cabsf(strength);
 
  float complex *restrict stamp =
    calloc(sizeof(float complex), (size_t)stamp_extent->width * stamp_extent->height);
  if(IS_NULL_PTR(stamp)) return 1;
 
  // lookup table: map of distance from center point => warp
  const int table_size = iradius * LOOKUP_OVERSAMPLE;
  const float *restrict lookup_table = build_lookup_table(table_size, warp->control1, warp->control2);
  if(IS_NULL_PTR(lookup_table))
  {
    dt_free(stamp);
    return 1;
  }
 
  // points into buffer at the center of the circle
  float complex *const center = stamp + 2 * iradius * iradius + 2 * iradius;
 
  // The expensive operation here is hypotf ().  By dividing the
  // circle in quadrants and doing only the inside we have to calculate
  // hypotf only for PI / 16 = 0.196 of the stamp area.
  // We don't do octants to avoid false sharing of cache lines between threads.
  // doesn't work for OSX see issue #7349
  __OMP_PARALLEL_FOR__(firstprivate(LOOKUP_OVERSAMPLE))
  for(int y = 0; y <= iradius; y++)
  {
    for(int x = 0; x <= iradius; x++)
    {
      const float dist = dt_fast_hypotf(x, y);
      const int idist = round(dist * LOOKUP_OVERSAMPLE);
      if(idist >= table_size)
        // idist will only grow bigger in this row
        break;
 
      // pointers into the 4 quadrants of the circle
      // quadrant count is ccw from positive x-axis
      float complex *const q1 = center - y * stamp_extent->width + x;
      float complex *const q2 = center - y * stamp_extent->width - x;
      float complex *const q3 = center + y * stamp_extent->width - x;
      float complex *const q4 = center + y * stamp_extent->width + x;
 
      float abs_lookup = abs_strength * lookup_table[idist] / iradius;
 
      switch (warp->type)
      {
         case DT_LIQUIFY_WARP_TYPE_RADIAL_GROW:
           *q1 = abs_lookup * ( x - y * I);
           *q2 = abs_lookup * (-x - y * I);
           *q3 = abs_lookup * (-x + y * I);
           *q4 = abs_lookup * ( x + y * I);
           break;
 
         case DT_LIQUIFY_WARP_TYPE_RADIAL_SHRINK:
           *q1 = -abs_lookup * ( x - y * I);
           *q2 = -abs_lookup * (-x - y * I);
           *q3 = -abs_lookup * (-x + y * I);
           *q4 = -abs_lookup * ( x + y * I);
           break;
 
         default:
           *q1 = *q2 = *q3 = *q4 = strength * lookup_table[idist];
           break;
      }
    }
  }
 
  dt_pixelpipe_cache_free_align(lookup_table);
  *pstamp = stamp;
  return 0;
}
 
/*
  Applies a stamp at a specified position.
 
  Applies a stamp at the position specified by @a point and adds the
  resulting vector field to the global distortion map @a global_map.
 
  The global distortion map is a map of relative pixel displacements
  encompassing all our paths.
*/
 
static void add_to_global_distortion_map(float complex *global_map,
                                          const cairo_rectangle_int_t *const restrict global_map_extent,
                                          const dt_liquify_warp_t *const restrict warp,
                                          const float complex *const restrict stamp,
                                          const cairo_rectangle_int_t *stamp_extent)
{
  cairo_rectangle_int_t mmext = *stamp_extent;
  mmext.x += (int) round(crealf(warp->point));
  mmext.y += (int) round(cimagf(warp->point));
  cairo_rectangle_int_t cmmext = mmext;
  cairo_region_t *mmreg = cairo_region_create_rectangle(&mmext);
  cairo_region_intersect_rectangle(mmreg, global_map_extent);
  cairo_region_get_extents(mmreg, &cmmext);
  dt_free(mmreg);
 
  #ifdef _OPENMP
  #pragma omp parallel for schedule (static) default (shared)
  #endif
 
  for(int y = cmmext.y; y < cmmext.y + cmmext.height; y++)
  {
    const float complex *const srcrow = stamp + ((y - mmext.y) * mmext.width);
    float complex *const destrow = global_map + ((y - global_map_extent->y) * global_map_extent->width);
 
    for(int x = cmmext.x; x < cmmext.x + cmmext.width; x++)
    {
      destrow[x - global_map_extent->x] -= srcrow[x - mmext.x];
    }
  }
}
 
/*
  Applies the global distortion map to the picture.  The distortion
  map maps points to the position from where the new color of the
  point should be sampled from.  The distortion map is in relative
  device coords.
*/
 
__DT_CLONE_TARGETS__
static void apply_global_distortion_map(struct dt_iop_module_t *module,
                                         const dt_dev_pixelpipe_iop_t *piece,
                                         const float *const restrict in,
                                         float *const restrict out,
                                         const dt_iop_roi_t *const roi_in,
                                         const dt_iop_roi_t *const roi_out,
                                         const int ch,
                                         const float complex *const map,
                                         const cairo_rectangle_int_t *extent)
{
  const int ch_width = ch * roi_in->width;
  const struct dt_interpolation * const interpolation =
    dt_interpolation_new(DT_INTERPOLATION_USERPREF_WARP);
 
  #ifdef _OPENMP
  #pragma omp parallel for schedule (static) default (shared)
  #endif
 
  for(int y = extent->y; y < extent->y + extent->height; y++)
  {
    // point inside roi_out ?
    if(y >= roi_out->y && y < roi_out->y + roi_out->height)
    {
      const float complex *row = map + (y - extent->y) * extent->width;
      float* out_sample = out + ((y - roi_out->y) * roi_out->width +
                               extent->x - roi_out->x) * ch;
      for(int x = extent->x; x < extent->x + extent->width; x++)
      {
        if(
          // point inside roi_out ?
          (x >= roi_out->x && x < roi_out->x + roi_out->width) &&
          // point actually warped ?
          (*row != 0))
        {
          if(ch == 1)
            *out_sample = dt_interpolation_compute_sample(interpolation,
                                                          in,
                                                          x + crealf(*row) - roi_in->x,
                                                          y + cimagf(*row) - roi_in->y,
                                                          roi_in->width,
                                                          roi_in->height,
                                                          ch,
                                                          ch_width);
          else
            dt_interpolation_compute_pixel4c(
              interpolation,
              in,
              out_sample,
              x + crealf(*row) - roi_in->x,
              y + cimagf(*row) - roi_in->y,
              roi_in->width,
              roi_in->height,
              ch_width);
 
        }
        ++row;
        out_sample += ch;
      }
    }
  }
}
 
// calculate the map extent.
 
static GSList *_get_map_extent(const dt_iop_roi_t *roi_out,
                               const GList *interpolated,
                               cairo_rectangle_int_t *map_extent)
{
  const cairo_rectangle_int_t roi_out_rect = { roi_out->x, roi_out->y, roi_out->width, roi_out->height };
  cairo_region_t *roi_out_region = cairo_region_create_rectangle(&roi_out_rect);
  cairo_region_t *map_region = cairo_region_create();
  GSList *in_roi = NULL;
 
  for(const GList *i = interpolated; i; i = g_list_next(i))
  {
    const dt_liquify_warp_t *warp = ((dt_liquify_warp_t *) i->data);
    cairo_rectangle_int_t r;
    compute_round_stamp_extent(&r, warp);
    // add extent if not entirely outside the roi
    if(cairo_region_contains_rectangle(roi_out_region, &r) != CAIRO_REGION_OVERLAP_OUT)
    {
      cairo_region_union_rectangle(map_region, &r);
      in_roi = g_slist_prepend(in_roi, i->data);
    }
  }
 
  // return the paths and the extent of all paths
  cairo_region_get_extents(map_region, map_extent);
  cairo_region_destroy(map_region);
  cairo_region_destroy(roi_out_region);
 
  return g_slist_reverse(in_roi);
}
 
static float complex *create_global_distortion_map(const cairo_rectangle_int_t *map_extent,
                                                   const GSList *interpolated,
                                                   gboolean inverted)
{
  const int mapsize = map_extent->width * map_extent->height;
  if (mapsize == 0)
  {
    // there are no pixels for which we need distortion info, so return right away
    // caller will see the NULL and bypass any further processing of the points it wants to distort
    return NULL;
  }
 
  // allocate distortion map big enough to contain all paths
  float complex *map = dt_pixelpipe_cache_alloc_align_cache(sizeof(float complex) * mapsize, 0);
  if(IS_NULL_PTR(map)) return NULL;
  memset(map, 0, sizeof(float complex) * mapsize);
 
  // build map
  for(const GSList *i = interpolated; i; i = g_slist_next(i))
  {
    const dt_liquify_warp_t *warp = ((dt_liquify_warp_t *) i->data);
    float complex *stamp = NULL;
    cairo_rectangle_int_t r;
    if(build_round_stamp(&stamp, &r, warp) != 0)
    {
      dt_pixelpipe_cache_free_align(map);
      return NULL;
    }
    add_to_global_distortion_map(map, map_extent, warp, stamp, &r);
    dt_free(stamp);
  }
 
  if(inverted)
  {
    float complex * const imap = dt_pixelpipe_cache_alloc_align_cache(sizeof(float complex) * mapsize, 0);
    if(IS_NULL_PTR(imap))
    {
      dt_pixelpipe_cache_free_align(map);
      return NULL;
    }
    memset(imap, 0, sizeof(float complex) * mapsize);
 
    // copy map into imap(inverted map).
    // imap [ n + dx(map[n]) , n + dy(map[n]) ] = -map[n]
 
    #ifdef _OPENMP
    #pragma omp parallel for schedule (static) default (shared)
    #endif
 
    for(int y = 0; y <  map_extent->height; y++)
    {
      const float complex *const row = map + y * map_extent->width;
      for(int x = 0; x < map_extent->width; x++)
      {
        const float complex d = row[x];
        // compute new position (nx,ny) given the displacement d
        const int nx = x + (int)crealf(d);
        const int ny = y + (int)cimagf(d);
 
        // if the point falls into the extent, set it
        if(nx>0 && nx<map_extent->width && ny>0 && ny<map_extent->height)
          imap[nx + ny * map_extent->width] = -d;
      }
    }
 
    dt_pixelpipe_cache_free_align(map);
 
    // now just do a pass to avoid gap with a displacement of zero, note that we do not need high
    // precision here as the inverted distortion mask is only used to compute a final displacement
    // of points.
 
    #ifdef _OPENMP
    #pragma omp parallel for schedule (static) default (shared)
    #endif
 
    for(int y = 0; y <  map_extent->height; y++)
    {
      float complex *const row = imap + y * map_extent->width;
      float complex last[2] = { 0, 0 };
      for(int x = 0; x < map_extent->width / 2 + 1; x++)
      {
        float complex *cl = row + x;
        float complex *cr = row + map_extent->width - x;
        if(x!=0)
        {
          if(*cl == 0) *cl = last[0];
          if(*cr == 0) *cr = last[1];
        }
        last[0] = *cl; last[1] = *cr;
      }
    }
 
    map = imap;
  }
  return map;
}
 
static float complex *build_global_distortion_map(struct dt_iop_module_t *module,
                                                   const dt_dev_pixelpipe_t *pipe,
                                                   const dt_dev_pixelpipe_iop_t *piece,
                                                   const dt_iop_roi_t *roi_in,
                                                   const dt_iop_roi_t *roi_out,
                                                   cairo_rectangle_int_t *map_extent)
{
  // copy params
  dt_iop_liquify_params_t copy_params;
  memcpy(&copy_params, (dt_iop_liquify_params_t *)piece->data, sizeof(dt_iop_liquify_params_t));
 
  distort_paths_raw_to_piece(module, pipe, roi_in->scale, &copy_params, FALSE);
 
  GList *interpolated = interpolate_paths(&copy_params);
  GSList *interpolated_in_roi = _get_map_extent(roi_out, interpolated, map_extent);
 
  float complex *map = create_global_distortion_map(map_extent, interpolated_in_roi, FALSE);
 
  g_slist_free(interpolated_in_roi);
  interpolated_in_roi = NULL;
  g_list_free_full(interpolated, dt_free_gpointer);
  interpolated = NULL;
  return map;
}
 
// 1st pass: how large would the output be, given this input roi?
// this is always called with the full buffer before processing.
void modify_roi_out(struct dt_iop_module_t *module, const struct dt_dev_pixelpipe_t *pipe,
                    struct dt_dev_pixelpipe_iop_t *piece, dt_iop_roi_t *roi_out,
                    const dt_iop_roi_t *roi_in)
{
  // output is same size as input
  *roi_out = *roi_in;
}
 
// 2nd pass: which roi would this operation need as input to fill the given output region?
void modify_roi_in(struct dt_iop_module_t *module, const struct dt_dev_pixelpipe_t *pipe,
                   struct dt_dev_pixelpipe_iop_t *piece, const dt_iop_roi_t *roi_out,
                   dt_iop_roi_t *roi_in)
{
  // Because we move pixels, and we may have to sample a pixel from
  // outside roi_in, we need to expand roi_in to contain all our
  // paths.  But we may ignore paths completely outside of roi_out.
 
  *roi_in = *roi_out;
 
  // copy params
  dt_iop_liquify_params_t copy_params;
  memcpy(&copy_params, (dt_iop_liquify_params_t*)piece->data, sizeof(dt_iop_liquify_params_t));
 
  // Keep ROI expansion on the same processing pipe as piece->data. Preview-only state can be NULL
  // in workers, and mixing pipes here desynchronizes the distortion graph from the current piece.
  distort_paths_raw_to_piece(module, pipe, roi_in->scale, &copy_params, FALSE);
 
  cairo_rectangle_int_t pipe_rect =
    {
      0,
      0,
      lroundf(piece->buf_in.width * roi_in->scale),
      lroundf(piece->buf_in.height * roi_in->scale)
    };
 
  cairo_rectangle_int_t roi_in_rect =
    {
      roi_in->x,
      roi_in->y,
      roi_in->width,
      roi_in->height
    };
  cairo_region_t *roi_in_region = cairo_region_create_rectangle(&roi_in_rect);
 
  // get extent of all paths
  GList *interpolated = interpolate_paths(&copy_params);
  cairo_rectangle_int_t extent;
  GSList *interpolated_in_roi = _get_map_extent(roi_out, interpolated, &extent);
  g_slist_free(interpolated_in_roi);
  interpolated_in_roi = NULL;
  g_list_free_full(interpolated, dt_free_gpointer);
  interpolated = NULL;
 
  // (eventually) extend roi_in
  cairo_region_union_rectangle(roi_in_region, &extent);
  // and clamp to pipe extent
  cairo_region_intersect_rectangle(roi_in_region, &pipe_rect);
 
  // write new extent to roi_in
  cairo_region_get_extents(roi_in_region, &roi_in_rect);
  roi_in->x = roi_in_rect.x;
  roi_in->y = roi_in_rect.y;
  roi_in->width  = roi_in_rect.width;
  roi_in->height = roi_in_rect.height;
 
  // cleanup
  cairo_region_destroy(roi_in_region);
}
 
static int _liquify_warp_points(const dt_iop_liquify_params_t *const params_in,
                                dt_iop_module_t *self, const dt_dev_pixelpipe_t *pipe,
                                dt_geometry_chain_t *chain, float *const restrict points,
                                const size_t points_count, const gboolean inverted)
{
 
  // compute the extent of all points (all computations are done in RAW coordinate)
  float xmin = FLT_MAX, xmax = FLT_MIN, ymin = FLT_MAX, ymax = FLT_MIN;
  __OMP_PARALLEL_FOR__(if(points_count > 100) reduction(min:xmin, ymin) reduction(max:xmax, ymax))
  for(size_t i = 0; i < points_count * 2; i += 2)
  {
    const float x = points[i];
    const float y = points[i + 1];
    xmin = fmin(xmin, x);
    xmax = fmax(xmax, x);
    ymin = fmin(ymin, y);
    ymax = fmax(ymax, y);
  }
 
  cairo_rectangle_int_t extent = { .x = (int)(xmin - .5), .y = (int)(ymin - .5),
                                   .width = (int)(xmax - xmin + 2.5), .height = (int)(ymax - ymin + 2.5) };
 
  if(extent.width > 0 && extent.height > 0)
  {
    // copy params
    dt_iop_liquify_params_t copy_params;
    memcpy(&copy_params, params_in, sizeof(dt_iop_liquify_params_t));
 
    if(!IS_NULL_PTR(chain))
      distort_paths_raw_to_piece_chain(self, chain, 1.f, &copy_params);
    else
      distort_paths_raw_to_piece(self, pipe, 1.f, &copy_params, TRUE);
 
    // create the distortion map for this extent
 
    GList *interpolated = interpolate_paths(&copy_params);
 
    // we need to adjust the extent to be the union enclosing all the points (currently in extent) and
    // the warps that are in (possibly partly) in this same region.
 
    dt_iop_roi_t roi_in = { .x = extent.x, .y = extent.y, .width = extent.width, .height = extent.height };
    GSList *interpolated_in_roi = _get_map_extent(&roi_in, interpolated, &extent);
 
    float complex *map = create_global_distortion_map(&extent, interpolated_in_roi, inverted);
    g_slist_free(interpolated_in_roi);
    interpolated_in_roi = NULL;
    g_list_free_full(interpolated, dt_free_gpointer);
    interpolated = NULL;
 
    if(IS_NULL_PTR(map)) return 0;
 
    const int map_size =  extent.width * extent.height;
    const int x_last = extent.x + extent.width;
    const int y_last = extent.y + extent.height;
 
    // apply distortion to all points (this is a simple displacement given by a vector at this same point in the map)
    __OMP_PARALLEL_FOR__(if(points_count > 100))
    for(size_t i = 0; i < points_count; i++)
    {
      float *px = &points[i*2];
      float *py = &points[i*2+1];
      const float x = *px;
      const float y = *py;
      const int map_offset = ((int)(x - 0.5) - extent.x) + ((int)(y - 0.5) - extent.y) * extent.width;
 
      if(x >= extent.x && x < x_last && y >= extent.y && y < y_last && map_offset >= 0 && map_offset < map_size)
      {
        const float complex dist = map[map_offset];
        *px += crealf(dist);
        *py += cimagf(dist);
      }
    }
 
    dt_pixelpipe_cache_free_align(map);
  }
 
  return 1;
}
 
static void start_drag(dt_iop_liquify_gui_data_t *g, dt_liquify_layer_enum_t layer, dt_liquify_path_data_t *elem)
{
  g->dragging.layer = layer;
  g->dragging.elem = elem;
}
 
static void end_drag(dt_iop_liquify_gui_data_t *g)
{
  g->dragging = NOWHERE;
}
 
static gboolean is_dragging(const dt_iop_liquify_gui_data_t *g)
{
  return !IS_NULL_PTR(g->dragging.elem);
}
 
int distort_transform(dt_iop_module_t *self, const dt_dev_pixelpipe_t *pipe, const dt_dev_pixelpipe_iop_t *piece,
                      float *const restrict points, size_t points_count)
{
  // Recurse on the caller pipe so we reuse the same node list and piece data in preview, export,
  // thumbnail, and mask evaluation paths.
  return _liquify_warp_points((const dt_iop_liquify_params_t *)piece->data, self, pipe, NULL, points,
                              points_count, TRUE);
}
 
int distort_backtransform(dt_iop_module_t *self, const dt_dev_pixelpipe_t *pipe, const dt_dev_pixelpipe_iop_t *piece,
                          float *const restrict points, size_t points_count)
{
  return _liquify_warp_points((const dt_iop_liquify_params_t *)piece->data, self, pipe, NULL, points,
                              points_count, FALSE);
}
 
/* --- the geometry service's view of this module (develop/geometry/geometry.h) ---------
 *
 * modify_roi_out() is the identity, so there is no size map: liquify moves pixels without
 * changing how many there are. What it does need is the one thing no other record needs -- the
 * chain composed around it. Its warps are stored in RAW sensor coordinates, so before it can
 * rasterise anything it has to push its own path nodes through every module upstream of itself,
 * which on the pixel pipe means re-entering the pipe walker mid-walk and here means re-entering
 * the chain (dt_geometry_chain_compose(), bounded BACK_EXCL of this module's own iop_order so
 * the recursion terminates).
 *
 * The cost is inherited, not introduced: each query rasterises a warp map over the extent of
 * the points it was given, so it is O(area), not O(points), exactly as the pipe's own
 * distort_transform() has always been.
 */
 
static int _liquify_geometry_transform(const void *data, const dt_geometry_record_t *const record,
                                       dt_geometry_chain_t *chain, float *points, size_t points_count)
{
  const dt_iop_liquify_geometry_t *const g = (const dt_iop_liquify_geometry_t *)data;
  return _liquify_warp_points(&g->params, g->self, NULL, chain, points, points_count, TRUE);
}
 
static int _liquify_geometry_backtransform(const void *data, const dt_geometry_record_t *const record,
                                           dt_geometry_chain_t *chain, float *points, size_t points_count)
{
  const dt_iop_liquify_geometry_t *const g = (const dt_iop_liquify_geometry_t *)data;
  return _liquify_warp_points(&g->params, g->self, NULL, chain, points, points_count, FALSE);
}
 
static const dt_geometry_vtable_t _liquify_geometry_vtable = {
  .map_size = NULL,   // modify_roi_out() is the identity: liquify changes no dimensions
  .transform = _liquify_geometry_transform,
  .backtransform = _liquify_geometry_backtransform,
};
 
gboolean geometry_record(struct dt_iop_module_t *self, const void *params, dt_geometry_record_t *record)
{
  dt_iop_liquify_geometry_t *g
      = (dt_iop_liquify_geometry_t *)g_malloc0(sizeof(dt_iop_liquify_geometry_t));
  if(IS_NULL_PTR(g)) return FALSE;
 
  /* commit_params() is a straight copy of the parameters, so the record carries that copy. The
   * module pointer is kept because the fold needs its iop_order to bound the recursion and its
   * dev to reach the paths; it is the same module the chain is being built from, so it cannot
   * outlive the record. */
  memcpy(&g->params, params, sizeof(dt_iop_liquify_params_t));
  g->self = self;
 
  record->data = g;
  record->free_data = dt_free_gpointer;
  record->vtable = &_liquify_geometry_vtable;
  return TRUE;
}
 
void distort_mask(struct dt_iop_module_t *self, const struct dt_dev_pixelpipe_t *pipe, struct dt_dev_pixelpipe_iop_t *piece,
                  const float *const in, float *const out, const dt_iop_roi_t *const roi_in,
                  const dt_iop_roi_t *const roi_out)
{
  (void)pipe;
  // 1. copy the whole image (we'll change only a small part of it)
  __OMP_PARALLEL_FOR__()
  for(int i = 0; i < roi_out->height; i++)
  {
    float *destrow = out + (size_t) i * roi_out->width;
    const float *srcrow = in + (size_t) (roi_in->width * (i + roi_out->y - roi_in->y) + roi_out->x - roi_in->x);
 
    memcpy(destrow, srcrow, sizeof(float) * roi_out->width);
  }
 
  // 2. build the distortion map
 
  cairo_rectangle_int_t map_extent;
  float complex *map = build_global_distortion_map(self, pipe, piece, roi_in, roi_out, &map_extent);
  if(IS_NULL_PTR(map))
  {
    if(map_extent.width != 0 && map_extent.height != 0) return;
    return;
  }
 
  // 3. apply the map
 
  if(map_extent.width != 0 && map_extent.height != 0)
    apply_global_distortion_map(self, piece, in, out, roi_in, roi_out, 1, map, &map_extent);
 
  dt_pixelpipe_cache_free_align(map);
 
}
 
__DT_CLONE_TARGETS__
int process(struct dt_iop_module_t *module, const dt_dev_pixelpipe_t *pipe, const dt_dev_pixelpipe_iop_t *piece,
            const void *const in,
             void *const out)
{
  const dt_iop_roi_t *const roi_in = &piece->roi_in;
  const dt_iop_roi_t *const roi_out = &piece->roi_out;
  // 1. copy the whole image (we'll change only a small part of it)
 
  const int ch = 4;
  const int height = MIN(roi_in->height, roi_out->height);
  const int width = MIN(roi_in->width, roi_out->width);
  __OMP_PARALLEL_FOR__()
  for(int i = 0; i < height; i++)
  {
    float *destrow = (float *)out + (size_t)ch * i * roi_out->width;
    const float *srcrow = (float *)in + (size_t)ch * (roi_in->width * (i + roi_out->y - roi_in->y) +
                                                       roi_out->x - roi_in->x);
 
    memcpy(destrow, srcrow, sizeof(float) * ch * width);
  }
 
  // 2. build the distortion map
 
  cairo_rectangle_int_t map_extent;
  float complex *map = build_global_distortion_map(module, pipe, piece, roi_in, roi_out, &map_extent);
  if(IS_NULL_PTR(map))
  {
    if(map_extent.width != 0 && map_extent.height != 0) return 1;
    return 0;
  }
 
  // 3. apply the map
 
  if(map_extent.width != 0 && map_extent.height != 0)
    apply_global_distortion_map(module, piece, in, out, roi_in, roi_out, 1, map, &map_extent);
 
  dt_pixelpipe_cache_free_align(map);
 
  return 0;
}
 
#ifdef HAVE_OPENCL
 
// compute bicubic kernel. See: https://en.wikipedia.org/wiki/Bicubic_interpolation#Bicubic_convolution_algorithm
 
static inline float bicubic(const float a, const float x)
{
  const float absx = fabsf(x);
  if(absx <= 1) return ((a + 2) * absx - (a + 3)) * absx * absx + 1;
  if(absx < 2) return ((a * absx - 5 * a) * absx + 8 * a) * absx - 4 * a;
  return 0.0f;
}
 
// compute Mitchell-Netravali cubic kernel (B=C=1/3): sharp but effectively
// halo-free, unlike Catmull-Rom bicubic and Lanczos which overshoot at edges.
 
static inline float mitchell(const float x)
{
  const float absx = fabsf(x);
  const float x2 = absx * absx;
  const float x3 = x2 * absx;
  if(absx < 1.0f) return (7.0f / 6.0f) * x3 - 2.0f * x2 + 8.0f / 9.0f;
  if(absx < 2.0f) return -(7.0f / 18.0f) * x3 + 2.0f * x2 - (10.0f / 3.0f) * absx + 16.0f / 9.0f;
  return 0.0f;
}
 
typedef struct
{
  int size;
  int resolution;
} dt_liquify_kernel_descriptor_t;
 
typedef cl_mem cl_mem_t;
typedef cl_int cl_int_t;
 
static cl_int_t apply_global_distortion_map_cl(struct dt_iop_module_t *module,
                                                const dt_dev_pixelpipe_t *pipe,
                                                const dt_dev_pixelpipe_iop_t *piece,
                                                const cl_mem_t dev_in,
                                                const cl_mem_t dev_out,
                                                const dt_iop_roi_t *roi_in,
                                                const dt_iop_roi_t *roi_out,
                                                const float complex *map,
                                                const cairo_rectangle_int_t *map_extent)
{
  cl_int_t err = CL_MEM_OBJECT_ALLOCATION_FAILURE;
 
  dt_iop_liquify_global_data_t *gd = (dt_iop_liquify_global_data_t *) module->global_data;
  const int devid = pipe->devid;
 
  const struct dt_interpolation* interpolation = dt_interpolation_new(DT_INTERPOLATION_USERPREF_WARP);
  dt_liquify_kernel_descriptor_t kdesc = { .size = 0, .resolution = 100 };
  float *k = NULL;
 
  switch (interpolation->id)
  {
     case DT_INTERPOLATION_BILINEAR:
       kdesc.size = 1;
       kdesc.resolution = 1;
       k = malloc(sizeof(float) * 2);
       if(k)
       {
        k[0] = 1.0f;
        k[1] = 0.0f;
       }
       break;
     case DT_INTERPOLATION_BICUBIC:
       kdesc.size = 2;
       k = malloc(sizeof(float) * ((size_t)kdesc.size * kdesc.resolution + 1));
       if(k)
        for(int i = 0; i <= kdesc.size * kdesc.resolution; ++i)
          k[i] = bicubic(0.5f, (float) i / kdesc.resolution);
       break;
     case DT_INTERPOLATION_MITCHELL:
       kdesc.size = 2;
       k = malloc(sizeof(float) * ((size_t)kdesc.size * kdesc.resolution + 1));
       if(k)
        for(int i = 0; i <= kdesc.size * kdesc.resolution; ++i)
          k[i] = mitchell((float) i / kdesc.resolution);
       break;
     default:
       return FALSE;
  }
 
  cl_mem_t dev_roi_in = dt_opencl_copy_host_to_device_constant
    (devid, sizeof(dt_iop_roi_t), (void *) roi_in);
 
  cl_mem_t dev_roi_out = dt_opencl_copy_host_to_device_constant
    (devid, sizeof(dt_iop_roi_t), (void *) roi_out);
 
  cl_mem_t dev_map = dt_opencl_copy_host_to_device_constant
    (devid, sizeof(float complex) * map_extent->width * map_extent->height, (void *) map);
 
  cl_mem_t dev_map_extent = dt_opencl_copy_host_to_device_constant
    (devid, sizeof(cairo_rectangle_int_t), (void *) map_extent);
 
  cl_mem_t dev_kdesc = dt_opencl_copy_host_to_device_constant
    (devid, sizeof(dt_liquify_kernel_descriptor_t), (void *) &kdesc);
 
  cl_mem_t dev_kernel = dt_opencl_copy_host_to_device_constant
    (devid, sizeof(float) * (kdesc.size * kdesc.resolution  + 1), (void *) k);
 
  if(IS_NULL_PTR(dev_roi_in) || IS_NULL_PTR(dev_roi_out) || IS_NULL_PTR(dev_map) || IS_NULL_PTR(dev_map_extent)
      || IS_NULL_PTR(dev_kdesc) || IS_NULL_PTR(dev_kernel) || IS_NULL_PTR(k))
    goto error;
 
  dt_opencl_set_kernel_arg(devid, gd->warp_kernel, 0, sizeof(cl_mem), &dev_in);
  dt_opencl_set_kernel_arg(devid, gd->warp_kernel, 1, sizeof(cl_mem), &dev_out);
  dt_opencl_set_kernel_arg(devid, gd->warp_kernel, 2, sizeof(cl_mem), &dev_roi_in);
  dt_opencl_set_kernel_arg(devid, gd->warp_kernel, 3, sizeof(cl_mem), &dev_roi_out);
  dt_opencl_set_kernel_arg(devid, gd->warp_kernel, 4, sizeof(cl_mem), &dev_map);
  dt_opencl_set_kernel_arg(devid, gd->warp_kernel, 5, sizeof(cl_mem), &dev_map_extent);
 
  dt_opencl_set_kernel_arg(devid, gd->warp_kernel, 6, sizeof(cl_mem), &dev_kdesc);
  dt_opencl_set_kernel_arg(devid, gd->warp_kernel, 7, sizeof(cl_mem), &dev_kernel);
 
  const size_t sizes[] = { ROUNDUPDWD(map_extent->width, devid), ROUNDUPDHT(map_extent->height, devid) };
  err = dt_opencl_enqueue_kernel_2d(devid, gd->warp_kernel, sizes);
 
error:
 
  dt_opencl_release_mem_object(dev_kernel);
  dt_opencl_release_mem_object(dev_kdesc);
  dt_opencl_release_mem_object(dev_map_extent);
  dt_opencl_release_mem_object(dev_map);
  dt_opencl_release_mem_object(dev_roi_out);
  dt_opencl_release_mem_object(dev_roi_in);
  dt_free(k);
 
  return err;
}
 
int process_cl(struct dt_iop_module_t *module, const dt_dev_pixelpipe_t *pipe,
                const dt_dev_pixelpipe_iop_t *piece,
                const cl_mem_t dev_in,
                const cl_mem_t dev_out)
{
  const dt_iop_roi_t *const roi_in = &piece->roi_in;
  const dt_iop_roi_t *const roi_out = &piece->roi_out;
  cl_int_t err = -999;
  const int devid = pipe->devid;
  const int height = MIN(roi_in->height, roi_out->height);
  const int width = MIN(roi_in->width, roi_out->width);
 
  // 1. copy the whole image (we'll change only a small part of it)
  {
    size_t src[]    = { roi_out->x - roi_in->x, roi_out->y - roi_in->y, 0 };
    size_t dest[]   = { 0, 0, 0 };
    size_t extent[] = { width, height, 1 };
    err = dt_opencl_enqueue_copy_image(devid, dev_in, dev_out, src, dest, extent);
    if(err != CL_SUCCESS) goto error;
  }
 
  // 2. build the distortion map
  cairo_rectangle_int_t map_extent;
  const float complex *map = build_global_distortion_map(module, pipe, piece, roi_in, roi_out, &map_extent);
  if(IS_NULL_PTR(map))
  {
    if(map_extent.width != 0 && map_extent.height != 0) return FALSE;
    return TRUE;
  }
 
  // 3. apply the map
  if(map_extent.width != 0 && map_extent.height != 0)
    err = apply_global_distortion_map_cl(module, pipe, piece, dev_in, dev_out, roi_in, roi_out, map, &map_extent);
  dt_pixelpipe_cache_free_align(map);
  if(err != CL_SUCCESS) goto error;
 
  return TRUE;
 
error:
  dt_print(DT_DEBUG_OPENCL, "[opencl_liquify] couldn't enqueue kernel! %d\n", err);
  return FALSE;
}
 
#endif
 
void init_global(dt_iop_module_so_t *module)
{
  // called once at startup
  const int program = 17; // from programs.conf
  dt_iop_liquify_global_data_t *gd = (dt_iop_liquify_global_data_t *) malloc(sizeof(dt_iop_liquify_global_data_t));
  module->data = gd;
  gd->warp_kernel = dt_opencl_create_kernel(program, "warp_kernel");
}
 
void cleanup_global(dt_iop_module_so_t *module)
{
  // called once at shutdown
  dt_iop_liquify_global_data_t *gd = (dt_iop_liquify_global_data_t *) module->data;
  dt_opencl_free_kernel(gd->warp_kernel);
  dt_free(module->data);
}
 
void init(dt_iop_module_t *module)
{
  // module is disabled by default
  module->default_enabled = 0;
  module->params_size = sizeof(dt_iop_liquify_params_t);
 
  // all allocated to 0, which is the default
  module->params = calloc(1, module->params_size);
  module->default_params = calloc(1, module->params_size);
}
 
void init_pipe(struct dt_iop_module_t *module, dt_dev_pixelpipe_t *pipe, dt_dev_pixelpipe_iop_t *piece)
{
  piece->data = dt_calloc_align(sizeof(dt_iop_liquify_params_t));
  piece->data_size = sizeof(dt_iop_liquify_params_t);
}
 
void cleanup_pipe(struct dt_iop_module_t *module, dt_dev_pixelpipe_t *pipe, dt_dev_pixelpipe_iop_t *piece)
{
  dt_free_align(piece->data);
  piece->data = NULL;
}
 
/* commit is the synch point between core and gui, so it copies params to pipe data. */
 
void commit_params(struct dt_iop_module_t *module,
                    dt_iop_params_t *params,
                    dt_dev_pixelpipe_t *pipe,
                    dt_dev_pixelpipe_iop_t *piece)
{
  memcpy(piece->data, params, module->params_size);
}
 
// calculate the dot product of 2 vectors.
 
static float cdot(const float complex p0, const float complex p1)
{
#ifdef FP_FAST_FMA
  return fma(crealf(p0), crealf(p1), cimagf(p0) * cimagf(p1));
#else
  return crealf(p0) * crealf(p1) + cimagf(p0) * cimagf(p1);
#endif
}
 
static void draw_rectangle(cairo_t *cr, const float complex pt, const double theta, const double size)
{
  const double x = creal(pt), y = cimag(pt);
  cairo_save(cr);
  cairo_translate(cr, x, y);
  cairo_rotate(cr, theta);
  cairo_rectangle(cr, -size / 2.0, -size / 2.0, size, size);
  cairo_restore(cr);
}
 
static void draw_triangle(cairo_t *cr, const float complex pt, const double theta, const double size)
{
  const double x = creal(pt), y = cimag(pt);
  cairo_save(cr);
  cairo_translate(cr, x, y);
  cairo_rotate(cr, theta);
  cairo_move_to(cr, -size, -size / 2.0);
  cairo_line_to(cr, 0,     0          );
  cairo_line_to(cr, -size, +size / 2.0);
  cairo_close_path(cr);
  cairo_restore(cr);
}
 
static void draw_circle(cairo_t *cr, const float complex pt, const double diameter)
{
  const double x = creal(pt), y = cimag(pt);
  cairo_save(cr);
  cairo_new_sub_path(cr);
  cairo_arc(cr, x, y, diameter / 2.0, 0, 2 * DT_M_PI);
  cairo_restore(cr);
}
 
static void set_source_rgba(cairo_t *cr, dt_liquify_rgba_t rgba)
{
  cairo_set_source_rgba(cr, rgba.red, rgba.green, rgba.blue, rgba.alpha);
}
 
static float get_ui_width(const float scale, const dt_liquify_ui_width_enum_t w)
{
  assert(w >= 0 && w < DT_LIQUIFY_UI_WIDTH_LAST);
  return scale * DT_PIXEL_APPLY_DPI(dt_liquify_ui_widths[w]);
}
 
#define GET_UI_WIDTH(a) (get_ui_width(scale, DT_LIQUIFY_UI_WIDTH_##a))
 
static void set_line_width(cairo_t *cr, double scale, dt_liquify_ui_width_enum_t w)
{
  const double width = get_ui_width(scale, w);
  cairo_set_line_width(cr, width);
}
 
static gboolean detect_drag(const dt_iop_liquify_gui_data_t *g, const double scale, const float complex pt)
{
  // g->last_button1_pressed_pos is valid only while BUTTON1 is down
  return g->last_button1_pressed_pos != -1.0 &&
    cabsf(pt - g->last_button1_pressed_pos) >= (GET_UI_WIDTH(MIN_DRAG) / scale);
}
 
static void update_warp_count(const dt_iop_liquify_gui_data_t *g)
{
  guint warp = 0, node = 0;
  for(int k=0; k<MAX_NODES; k++)
    if(g->params.nodes[k].header.type == DT_LIQUIFY_PATH_INVALIDATED)
      break;
    else
    {
      node++;
      if(g->params.nodes[k].header.type == DT_LIQUIFY_PATH_MOVE_TO_V1)
        warp++;
    }
  char str_warp[10], str_node[20];
  snprintf(str_warp, sizeof(str_warp), "%d", warp);
  snprintf(str_node, sizeof(str_node), "%d", node);
  gtk_label_set_text(g->label_warp, str_warp);
  gtk_label_set_text(g->label_node, str_node);
}
 
static GList *interpolate_paths(dt_iop_liquify_params_t *p)
{
  GList *l = NULL;
  for(int k=0; k<MAX_NODES; k++)
  {
    const dt_liquify_path_data_t *data = &p->nodes[k];
    if(data->header.type == DT_LIQUIFY_PATH_INVALIDATED)
      break;
 
    const float complex *p2 = &data->warp.point;
    const dt_liquify_warp_t *warp2 = &data->warp;
 
    if(data->header.type == DT_LIQUIFY_PATH_MOVE_TO_V1)
    {
      if(data->header.next == -1)
      {
        dt_liquify_warp_t *w = malloc(sizeof(dt_liquify_warp_t));
        *w = *warp2;
        l = g_list_append(l, w);
      }
      continue;
    }
 
    const dt_liquify_path_data_t *prev = node_prev(p, data);
    const dt_liquify_warp_t *warp1 = &prev->warp;
    const float complex *p1 = &prev->warp.point;
    if(data->header.type == DT_LIQUIFY_PATH_LINE_TO_V1)
    {
      const float total_length = cabsf(*p1 - *p2);
      float arc_length = 0.0f;
      while(arc_length < total_length)
      {
        dt_liquify_warp_t *w = malloc(sizeof(dt_liquify_warp_t));
        const float t = arc_length / total_length;
        const float complex pt = cmix(*p1, *p2, t);
        mix_warps(w, warp1, warp2, pt, t);
        w->status = DT_LIQUIFY_STATUS_INTERPOLATED;
        arc_length += cabsf(w->radius - w->point) * STAMP_RELOCATION;
        l = g_list_append(l, w);
      }
      continue;
    }
 
    if(data->header.type == DT_LIQUIFY_PATH_CURVE_TO_V1)
    {
      float complex *buffer = malloc(sizeof(float complex) * INTERPOLATION_POINTS);
      interpolate_cubic_bezier(*p1,
                                data->node.ctrl1,
                                data->node.ctrl2,
                                *p2,
                                buffer,
                                INTERPOLATION_POINTS);
      const float total_length = get_arc_length(buffer, INTERPOLATION_POINTS);
      float arc_length = 0.0f;
      restart_cookie_t restart = { 1, 0.0 };
 
      while(arc_length < total_length)
      {
        dt_liquify_warp_t *w = malloc(sizeof(dt_liquify_warp_t));
        const float t = arc_length / total_length;
        const float complex pt = point_at_arc_length(buffer, INTERPOLATION_POINTS, arc_length, &restart);
        mix_warps(w, warp1, warp2, pt, t);
        w->status = DT_LIQUIFY_STATUS_INTERPOLATED;
        arc_length += cabsf(w->radius - w->point) * STAMP_RELOCATION;
        l = g_list_append(l, w);
      }
      dt_free(buffer);
      continue;
    }
  }
  return l;
}
 
#define FG_COLOR     set_source_rgba(cr, fg_color)
#define BG_COLOR     set_source_rgba(cr, bg_color)
#define VERYTHINLINE set_line_width (cr, scale / 2.0f, DT_LIQUIFY_UI_WIDTH_THINLINE)
#define THINLINE     set_line_width (cr, scale, DT_LIQUIFY_UI_WIDTH_THINLINE)
#define THICKLINE    set_line_width (cr, scale, DT_LIQUIFY_UI_WIDTH_THICKLINE)
 
static void _draw_paths(dt_iop_module_t *module,
                        cairo_t *cr,
                        const float scale,
                        dt_iop_liquify_params_t *p,
                        GList *layers)
{
  const dt_iop_liquify_gui_data_t *g = (dt_iop_liquify_gui_data_t *)dt_iop_gui_data(module);
 
  cairo_set_line_cap(cr, CAIRO_LINE_CAP_ROUND);
 
  // do not display any iterpolated items as slow when:
  //   - we are dragging (pan)
  //   - the button one is pressed
  //   - exception for DT_LIQUIFY_LAYER_STRENGTHPOINT where we want to see the
  //     interpolated strength lines.
  GList *interpolated = (is_dragging(g) || g->last_button1_pressed_pos != -1)
    && (g->last_hit.layer != DT_LIQUIFY_LAYER_STRENGTHPOINT)
    ? NULL
    : interpolate_paths(p);
 
  for(const GList *l = layers; l; l = g_list_next(l))
  {
    const dt_liquify_layer_enum_t layer = (dt_liquify_layer_enum_t) GPOINTER_TO_INT(l->data);
    dt_liquify_rgba_t fg_color = dt_liquify_layers[layer].fg;
    dt_liquify_rgba_t bg_color = dt_liquify_layers[layer].bg;
 
    if(dt_liquify_layers[layer].opacity < 1.0)
      cairo_push_group(cr);
 
    for(int k=0; k<MAX_NODES; k++)
    {
      // this is an empty bin, old invalidated node, nothing more to do
      if(p->nodes[k].header.type == DT_LIQUIFY_PATH_INVALIDATED)
        break;
 
      dt_liquify_path_data_t *data = &p->nodes[k];
      const dt_liquify_path_data_t *prev = node_prev(p, data);
 
      if((dt_liquify_layers[layer].flags & DT_LIQUIFY_LAYER_FLAG_NODE_SELECTED)
         && !data->header.selected)
        continue;
 
      if((dt_liquify_layers[layer].flags & DT_LIQUIFY_LAYER_FLAG_PREV_SELECTED)
         && (IS_NULL_PTR(prev) || !prev->header.selected))
        continue;
 
      fg_color = dt_liquify_layers[layer].fg;
      bg_color = dt_liquify_layers[layer].bg;
 
      if(data->header.selected == layer)
        fg_color = DT_LIQUIFY_COLOR_SELECTED;
 
      if(data->header.hovered == dt_liquify_layers[layer].hover_master)
        fg_color = DT_LIQUIFY_COLOR_HOVER;
 
      cairo_new_path(cr);
 
      const float complex point = data->warp.point;
 
      if(data->header.type == DT_LIQUIFY_PATH_MOVE_TO_V1)
        cairo_move_to(cr, crealf(point), cimagf(point));
 
      if(layer == DT_LIQUIFY_LAYER_RADIUS)
      {
        for(const GList *i = interpolated; i; i = g_list_next(i))
        {
          const dt_liquify_warp_t *pwarp = ((dt_liquify_warp_t *) i->data);
          draw_circle(cr, pwarp->point, 2.0f * cabsf(pwarp->radius - pwarp->point));
        }
        draw_circle(cr, point, 2.0f * cabsf(data->warp.radius - data->warp.point));
        FG_COLOR;
        cairo_fill(cr);
      }
      else if(layer == DT_LIQUIFY_LAYER_HARDNESS1)
      {
        for(const GList *i = interpolated; i; i = g_list_next(i))
        {
          const dt_liquify_warp_t *pwarp = ((dt_liquify_warp_t *) i->data);
          draw_circle(cr, pwarp->point, 2.0f * cabsf(pwarp->radius - pwarp->point) * pwarp->control1);
        }
        FG_COLOR;
        cairo_fill(cr);
      }
      else if(layer == DT_LIQUIFY_LAYER_HARDNESS2)
      {
        for(const GList *i = interpolated; i; i = g_list_next(i))
        {
          const dt_liquify_warp_t *pwarp = ((dt_liquify_warp_t *) i->data);
          draw_circle(cr, pwarp->point, 2.0f * cabsf(pwarp->radius - pwarp->point) * pwarp->control2);
        }
        FG_COLOR;
        cairo_fill(cr);
      }
      else if(layer == DT_LIQUIFY_LAYER_WARPS)
      {
        VERYTHINLINE; FG_COLOR;
        for(const GList *i = interpolated; i; i = g_list_next(i))
        {
          const dt_liquify_warp_t *pwarp = ((dt_liquify_warp_t *) i->data);
          cairo_move_to(cr, crealf(pwarp->point), cimagf(pwarp->point));
          cairo_line_to(cr, crealf(pwarp->strength), cimagf(pwarp->strength));
        }
        cairo_stroke(cr);
 
        for(const GList *i = interpolated; i; i = g_list_next(i))
        {
          const dt_liquify_warp_t *pwarp = ((dt_liquify_warp_t *) i->data);
          const float rot = get_rot(pwarp->type);
          draw_circle(cr, pwarp->point, GET_UI_WIDTH(GIZMO_SMALL));
          draw_triangle(cr, pwarp->strength,
                        cargf(pwarp->strength - pwarp->point) + rot,
                        GET_UI_WIDTH(GIZMO_SMALL) / 3.0);
        }
        BG_COLOR;
        cairo_fill_preserve(cr);
        FG_COLOR;
        cairo_stroke(cr);
      }
      else if(layer == DT_LIQUIFY_LAYER_PATH)
      {
        if((data->header.type == DT_LIQUIFY_PATH_LINE_TO_V1)
            || (data->header.type == DT_LIQUIFY_PATH_CURVE_TO_V1))
        {
          assert(prev);
          cairo_move_to(cr, crealf(prev->warp.point), cimagf(prev->warp.point));
          if(data->header.type == DT_LIQUIFY_PATH_LINE_TO_V1)
            cairo_line_to(cr, crealf(point), cimagf(point));
          if(data->header.type == DT_LIQUIFY_PATH_CURVE_TO_V1)
          {
            cairo_curve_to(cr, crealf(data->node.ctrl1), cimagf(data->node.ctrl1),
                            crealf(data->node.ctrl2), cimagf(data->node.ctrl2),
                            crealf(point), cimagf(point));
          }
          THICKLINE; FG_COLOR;
          cairo_stroke_preserve(cr);
          THINLINE; BG_COLOR;
          cairo_stroke(cr);
        }
      }
      else if(layer == DT_LIQUIFY_LAYER_CENTERPOINT)
      {
        if(data->header.type == DT_LIQUIFY_PATH_MOVE_TO_V1
            || data->header.type == DT_LIQUIFY_PATH_LINE_TO_V1
            || data->header.type == DT_LIQUIFY_PATH_CURVE_TO_V1)
        {
          const float w = GET_UI_WIDTH(GIZMO);
          switch (data->header.node_type)
          {
             case DT_LIQUIFY_NODE_TYPE_CUSP:
               draw_triangle(cr, point - w / 2.0 * I, -DT_M_PI / 2.0, w);
               break;
             case DT_LIQUIFY_NODE_TYPE_SMOOTH:
               draw_rectangle(cr, point, DT_M_PI / 4.0, w);
               break;
             case DT_LIQUIFY_NODE_TYPE_SYMMETRICAL:
               draw_rectangle(cr, point, 0, w);
               break;
             case DT_LIQUIFY_NODE_TYPE_AUTOSMOOTH:
               draw_circle(cr, point, w);
               break;
             default:
               break;
          }
          THINLINE; BG_COLOR;
          cairo_fill_preserve(cr);
          FG_COLOR;
          cairo_stroke(cr);
        }
      }
 
      if(data->header.type == DT_LIQUIFY_PATH_CURVE_TO_V1)
      {
        if(layer == DT_LIQUIFY_LAYER_CTRLPOINT1_HANDLE &&
            !(!IS_NULL_PTR(prev) && prev->header.node_type == DT_LIQUIFY_NODE_TYPE_AUTOSMOOTH))
        {
          THINLINE; FG_COLOR;
          cairo_move_to(cr, crealf(prev->warp.point), cimagf(prev->warp.point));
          cairo_line_to(cr, crealf(data->node.ctrl1), cimagf(data->node.ctrl1));
          cairo_stroke(cr);
        }
        if(layer == DT_LIQUIFY_LAYER_CTRLPOINT2_HANDLE &&
            data->header.node_type != DT_LIQUIFY_NODE_TYPE_AUTOSMOOTH)
        {
          THINLINE; FG_COLOR;
          cairo_move_to(cr, crealf(data->warp.point), cimagf(data->warp.point));
          cairo_line_to(cr, crealf(data->node.ctrl2), cimagf(data->node.ctrl2));
          cairo_stroke(cr);
        }
        if(layer == DT_LIQUIFY_LAYER_CTRLPOINT1 &&
            !(!IS_NULL_PTR(prev) && prev->header.node_type == DT_LIQUIFY_NODE_TYPE_AUTOSMOOTH))
        {
          THINLINE; BG_COLOR;
          draw_circle(cr, data->node.ctrl1, GET_UI_WIDTH(GIZMO_SMALL));
          cairo_fill_preserve(cr);
          FG_COLOR;
          cairo_stroke(cr);
        }
        if(layer == DT_LIQUIFY_LAYER_CTRLPOINT2 &&
            data->header.node_type != DT_LIQUIFY_NODE_TYPE_AUTOSMOOTH)
        {
          THINLINE; BG_COLOR;
          draw_circle(cr, data->node.ctrl2, GET_UI_WIDTH(GIZMO_SMALL));
          cairo_fill_preserve(cr);
          FG_COLOR;
          cairo_stroke(cr);
        }
      }
 
      const dt_liquify_warp_t *warp  = &data->warp;
 
      if(layer == DT_LIQUIFY_LAYER_RADIUSPOINT_HANDLE)
      {
        draw_circle(cr, point, 2.0 * cabsf(warp->radius - point));
        THICKLINE; FG_COLOR;
        cairo_stroke_preserve(cr);
        THINLINE; BG_COLOR;
        cairo_stroke(cr);
      }
 
      if(layer == DT_LIQUIFY_LAYER_RADIUSPOINT)
      {
        THINLINE; BG_COLOR;
        draw_circle(cr, warp->radius, GET_UI_WIDTH(GIZMO_SMALL));
        cairo_fill_preserve(cr);
        FG_COLOR;
        cairo_stroke(cr);
      }
 
      if(layer == DT_LIQUIFY_LAYER_HARDNESSPOINT1_HANDLE)
      {
        draw_circle(cr, point, 2.0 * cabsf(warp->radius - point) * warp->control1);
        THICKLINE; FG_COLOR;
        cairo_stroke_preserve(cr);
        THINLINE; BG_COLOR;
        cairo_stroke(cr);
      }
 
      if(layer == DT_LIQUIFY_LAYER_HARDNESSPOINT2_HANDLE)
      {
        draw_circle(cr, point, 2.0 * cabsf(warp->radius - point) * warp->control2);
        THICKLINE; FG_COLOR;
        cairo_stroke_preserve(cr);
        THINLINE; BG_COLOR;
        cairo_stroke(cr);
      }
 
      if(layer == DT_LIQUIFY_LAYER_HARDNESSPOINT1)
      {
        draw_triangle(cr, cmix(point, warp->radius, warp->control1),
                      cargf(warp->radius - point),
                      GET_UI_WIDTH(GIZMO_SMALL));
        THINLINE; BG_COLOR;
        cairo_fill_preserve(cr);
        FG_COLOR;
        cairo_stroke(cr);
      }
 
      if(layer == DT_LIQUIFY_LAYER_HARDNESSPOINT2)
      {
        draw_triangle(cr, cmix(point, warp->radius, warp->control2),
                      cargf(-(warp->radius - point)),
                      GET_UI_WIDTH(GIZMO_SMALL));
        THINLINE; BG_COLOR;
        cairo_fill_preserve(cr);
        FG_COLOR;
        cairo_stroke(cr);
      }
 
      if(layer == DT_LIQUIFY_LAYER_STRENGTHPOINT_HANDLE)
      {
        cairo_move_to(cr, crealf(point), cimagf(point));
        if(warp->type == DT_LIQUIFY_WARP_TYPE_LINEAR)
        {
          const float complex pt = cmix(point, warp->strength,
                                        1.0 - 0.5
                                        * (GET_UI_WIDTH(GIZMO_SMALL)
                                           / cabsf(warp->strength - point)));
          cairo_line_to(cr, crealf(pt), cimagf(pt));
        }
        else
          draw_circle(cr, point, 2.0 * cabsf(warp->strength - warp->point));
        THICKLINE; FG_COLOR;
        cairo_stroke_preserve(cr);
        THINLINE; BG_COLOR;
        cairo_stroke(cr);
      }
 
      if(layer == DT_LIQUIFY_LAYER_STRENGTHPOINT)
      {
        cairo_move_to(cr, crealf(warp->strength), cimagf(warp->strength));
        const float rot = get_rot(warp->type);
        draw_triangle(cr, warp->strength,
                      cargf(warp->strength - warp->point) + rot,
                      GET_UI_WIDTH(GIZMO_SMALL));
        THINLINE; BG_COLOR;
        cairo_fill_preserve(cr);
        FG_COLOR;
        cairo_stroke(cr);
      }
    }
 
    if(dt_liquify_layers[layer].opacity < 1.0)
    {
      cairo_pop_group_to_source(cr);
      cairo_paint_with_alpha(cr, dt_liquify_layers[layer].opacity);
    }
  }
 
  g_list_free_full(interpolated, dt_free_gpointer);
  interpolated = NULL;
}
 
/*
  Find the nearest point on a cubic bezier curve.
 
  Return the curve parameter t of the point on a cubic bezier curve
  that is nearest to another arbitrary point.  Uses interpolation.
 
  FIXME: Implement a faster method, see:
  http://tog.acm.org/resources/GraphicsGems/gems/NearestPoint.c
*/
 
static float find_nearest_on_curve_t(const float complex p0,
                                      const float complex p1,
                                      const float complex p2,
                                      const float complex p3,
                                      const float complex x,
                                      const int n)
{
  float min_t = 0.0f, min_dist = cabsf(x - p0);
 
  for(int i = 0; i < n; i++)
  {
    const float t = (1.0 * i) / n;
    const float t1 = 1.0 - t;
    const float complex ip =
          t1 * t1 * t1 * p0 +
      3 * t1 * t1 * t  * p1 +
      3 * t1 * t  * t  * p2 +
          t  * t  * t  * p3;
 
    const float dist = cabsf(x - ip);
    if(dist < min_dist)
    {
      min_dist = dist;
      min_t = t;
    }
  }
  return min_t;
}
 
/*
  Find the nearest point on a line.
 
  Return the line parameter t of the point on a line that is nearest
  to another arbitrary point.
*/
 
static float find_nearest_on_line_t(const float complex p0, const float complex p1, const float complex x)
{
  // scalar projection
  const float b     = cabsf(p1 - p0);         // |b|
  const float dotab = cdot(x - p0, p1 - p0);  // |a| * |b| * cos(phi)
  return dotab / (b * b);                     // |a| / |b| * cos(phi)
}
 
// split a cubic bezier at t into two cubic beziers.
 
static void casteljau(const float complex *p0, float complex *p1, float complex *p2, float complex *p3, const float t)
{
  const float complex p01 = *p0 + (*p1 - *p0) * t;
  const float complex p12 = *p1 + (*p2 - *p1) * t;
  const float complex p23 = *p2 + (*p3 - *p2) * t;
 
  const float complex p012 = p01 + (p12 - p01) * t;
  const float complex p123 = p12 + (p23 - p12) * t;
 
  const float complex p0123 = p012 + (p123 - p012) * t;
 
  *p1 = p01;
  *p2 = p012;
  *p3 = p0123;
}
 
#define CHECK_HIT_PT(point)             \
  const float d = cabsf(point - (*pt)); \
  if(d < distance)                      \
  {                                     \
    distance = d;                       \
    hit.layer = layer;                  \
    hit.elem = data;                    \
  }
 
static dt_liquify_hit_t _hit_paths(dt_iop_module_t *module,
                                   dt_iop_liquify_params_t *p,
                                   GList *layers,
                                   const float complex *pt)
{
  dt_liquify_hit_t hit = NOWHERE;
 
  float distance = FLT_MAX;
 
  for(const GList *l = layers; l; l = g_list_next(l))
  {
    const dt_liquify_layer_enum_t layer = (dt_liquify_layer_enum_t) GPOINTER_TO_INT(l->data);
 
    if((dt_liquify_layers[layer].flags & DT_LIQUIFY_LAYER_FLAG_HIT_TEST) == 0)
      continue;
 
    for(int k=0; k<MAX_NODES; k++)
    {
      dt_liquify_path_data_t *data = &p->nodes[k];
      const dt_liquify_path_data_t *prev = node_prev(p, data);
 
      if(p->nodes[k].header.type == DT_LIQUIFY_PATH_INVALIDATED)
        break;
 
      if((dt_liquify_layers[layer].flags & DT_LIQUIFY_LAYER_FLAG_NODE_SELECTED)
          && !data->header.selected)
        continue;
 
      if((dt_liquify_layers[layer].flags & DT_LIQUIFY_LAYER_FLAG_PREV_SELECTED)
          && (IS_NULL_PTR(prev) || !prev->header.selected))
        continue;
 
      const dt_liquify_warp_t *warp  = &data->warp;
      const float complex point = data->warp.point;
 
      if(layer == DT_LIQUIFY_LAYER_PATH)
      {
        if(data->header.type == DT_LIQUIFY_PATH_LINE_TO_V1)
        {
          // remove 5% from start and end of line as non sensible area
          // this is to avoid wrong interaction for center point on both
          // sides.
          const float complex deadzone = (point - prev->warp.point) / 20.0f;
          const float complex lp1 = prev->warp.point + deadzone;
          const float complex lp2 = point - deadzone;
          const float t = find_nearest_on_line_t(lp1, lp2, *pt);
 
          if(t > 0.0f && t < 1.0f)
          {
            const float complex linepoint = cmix(lp1, lp2, t);
            const float d = cabsf(linepoint - *pt);
            if(d < distance)
            {
              distance = d;
              hit.layer = layer;
              hit.elem = data;
            }
          }
        }
        else if(data->header.type == DT_LIQUIFY_PATH_CURVE_TO_V1)
        {
          // remove 5% from start and end of line as non sensible area
          // this is to avoid wrong interaction for center point on both
          // sides.
          const float complex deadzone = (point - prev->warp.point) / 20.0f;
          const float complex lp1 = prev->warp.point + deadzone;
          const float complex lp2 = point - deadzone;
          const float t = find_nearest_on_curve_t(lp1, data->node.ctrl1, data->node.ctrl2, lp2, *pt, INTERPOLATION_POINTS);
 
          if(t > 0.0f && t < 1.0f)
          {
            float complex curvepoint = lp2;
            float complex p1 = data->node.ctrl1;
            float complex p2 = data->node.ctrl2;
            casteljau(&lp1, &p1, &p2, &curvepoint, t);
 
            const float d = cabsf(curvepoint - *pt);
            if(d < distance)
            {
              distance = d;
              hit.layer = layer;
              hit.elem = data;
            }
          }
        }
      }
      else if(layer == DT_LIQUIFY_LAYER_CENTERPOINT)
      {
        if(data->header.type == DT_LIQUIFY_PATH_MOVE_TO_V1
            || data->header.type == DT_LIQUIFY_PATH_LINE_TO_V1
            || data->header.type == DT_LIQUIFY_PATH_CURVE_TO_V1)
        {
          CHECK_HIT_PT(point);
        }
      }
      else if(layer == DT_LIQUIFY_LAYER_RADIUSPOINT)
      {
        CHECK_HIT_PT(warp->radius);
      }
      else if(layer == DT_LIQUIFY_LAYER_HARDNESSPOINT1)
      {
        CHECK_HIT_PT(cmix(point, warp->radius, warp->control1));
      }
      else if(layer == DT_LIQUIFY_LAYER_HARDNESSPOINT2)
      {
        CHECK_HIT_PT(cmix(point, warp->radius, warp->control2));
      }
      else if(layer == DT_LIQUIFY_LAYER_STRENGTHPOINT)
      {
        const float complex p = warp->point - warp->strength;
        CHECK_HIT_PT(warp->strength +  (float)DT_PIXEL_APPLY_DPI(5) * (p / cabsf(p)));
      }
 
      if(data->header.type == DT_LIQUIFY_PATH_CURVE_TO_V1)
      {
        if(layer == DT_LIQUIFY_LAYER_CTRLPOINT1 &&
            !(!IS_NULL_PTR(prev) && prev->header.node_type == DT_LIQUIFY_NODE_TYPE_AUTOSMOOTH))
        {
          CHECK_HIT_PT(data->node.ctrl1);
        }
        if(layer == DT_LIQUIFY_LAYER_CTRLPOINT2 &&
            data->header.node_type != DT_LIQUIFY_NODE_TYPE_AUTOSMOOTH)
        {
          CHECK_HIT_PT(data->node.ctrl2);
        }
      }
    }
  }
 
  if(distance < DT_GUI_MOUSE_EFFECT_RADIUS * 0.5)
    return hit;
  else
    return NOWHERE;
}
 
static void draw_paths(struct dt_iop_module_t *module, cairo_t *cr, const float scale, dt_iop_liquify_params_t *params)
{
  const dt_iop_liquify_gui_data_t *g = (dt_iop_liquify_gui_data_t *)dt_iop_gui_data(module);
  GList *layers = NULL;
 
  for(dt_liquify_layer_enum_t layer = 0; layer < DT_LIQUIFY_LAYER_LAST; ++layer)
  {
    if(gtk_toggle_button_get_active(g->btn_point_tool)
        && (dt_liquify_layers[layer].flags & DT_LIQUIFY_LAYER_FLAG_POINT_TOOL))
      layers = g_list_prepend(layers, GINT_TO_POINTER(layer));
    if(gtk_toggle_button_get_active(g->btn_line_tool)
        && (dt_liquify_layers[layer].flags & DT_LIQUIFY_LAYER_FLAG_LINE_TOOL))
      layers = g_list_prepend(layers, GINT_TO_POINTER(layer));
    if(gtk_toggle_button_get_active(g->btn_curve_tool)
        && (dt_liquify_layers[layer].flags & DT_LIQUIFY_LAYER_FLAG_CURVE_TOOL))
      layers = g_list_prepend(layers, GINT_TO_POINTER(layer));
    if(gtk_toggle_button_get_active(g->btn_node_tool)
        && (dt_liquify_layers[layer].flags & DT_LIQUIFY_LAYER_FLAG_NODE_TOOL))
      layers = g_list_prepend(layers, GINT_TO_POINTER(layer));
  }
  layers = g_list_reverse(layers); // list was built in reverse order, so un-reverse it
 
  _draw_paths(module, cr, scale, params, layers);
 
  g_list_free(layers);
  layers = NULL;
}
 
static dt_liquify_hit_t _hit_test_paths(struct dt_iop_module_t *module,
                                        dt_iop_liquify_params_t *params,
                                        float complex pt)
{
  dt_liquify_hit_t hit = NOWHERE;
  GList *layers = NULL;
 
  for(dt_liquify_layer_enum_t layer = 0; layer < DT_LIQUIFY_LAYER_LAST; ++layer)
  {
    if(dt_liquify_layers[layer].flags & DT_LIQUIFY_LAYER_FLAG_HIT_TEST)
      layers = g_list_prepend(layers, GINT_TO_POINTER(layer));
  }
  layers = g_list_reverse(layers); // list was built in reverse order, so un-reverse it
 
  hit = _hit_paths(module, params, layers, &pt);
  g_list_free(layers);
  layers = NULL;
  return hit;
}
 
static void smooth_path_linsys(size_t n,
                                const float complex *k,
                                float complex *c1,
                                float complex *c2,
                                const int *equation)
{
  --n;
  float *a = malloc(sizeof(float) * n); // subdiagonal
  float *b = malloc(sizeof(float) * n); // main diagonal
  float *c = malloc(sizeof(float) * n); // superdiagonal
  float complex *d = malloc(sizeof(float complex) * n); // right hand side
 
  // Build the tridiagonal matrix.
 
  for(int i = 0; i < n; i++)
  {
    switch (equation[i])
    {
    #define ABCD(A,B,C,D) { { a[i] = A; b[i] = B; c[i] = C; d[i] = D; continue; } }
       case 1:  ABCD(0, 2, 1,       k[i] + 2 * k[i+1]   ); break;
       case 2:  ABCD(1, 4, 1,   4 * k[i] + 2 * k[i+1]   ); break;
       case 3:  ABCD(2, 7, 0,   8 * k[i] +     k[i+1]   ); break;
       case 4:  ABCD(0, 1, 0,                  c1[i]    ); break;
       case 5:  ABCD(0, 1, 0,                  c1[i]    ); break;
       case 6:  ABCD(1, 4, 0,   4 * k[i] +     c2[i]    ); break;
       case 7:  ABCD(0, 1, 0,                  c1[i]    ); break;
       case 8:  ABCD(0, 3, 0,   2 * k[i] +     k[i+1]   ); break;
       case 9:  ABCD(0, 2, 0,       k[i] +     c2[i]    ); break;
    #undef ABCD
    }
  }
 
  // Solve with the Thomas algorithm to compute c1's.  See:
  // http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
 
  for(int i = 1; i < n; i++)
  {
    const float m = a[i] / b[i-1];
    b[i] = b[i] - m * c[i-1];
    d[i] = d[i] - m * d[i-1];
  }
 
  c1[n-1] = d[n-1] / b[n-1];
  for(int i = n - 2; i >= 0; i--)
    c1[i] = (d[i] - c[i] * c1[i+1]) / b[i];
 
  // Now compute the c2's.
 
  for(int i = 0; i < n; i++)
  {
    switch (equation[i])
    {
       // keep end: c2 does not change
       case 5:
       case 6:
       case 9:  break;
 
       // straight end: put c2[i] halfway between c1[i] and k[i+1]
       case 3:
       case 7:
       case 8:  c2[i] = (c1[i] + k[i+1]) / 2;  break;
 
       // smooth end: c2 and c1 are symmetrical around the knot
       default: c2[i] = 2 * k[i+1] - c1[i+1];
    }
  }
 
  dt_free(a);
  dt_free(b);
  dt_free(c);
  dt_free(d);
}
 
static int path_length(dt_iop_liquify_params_t *p, dt_liquify_path_data_t *n)
{
  int count = 1;
  while(n->header.next != -1)
  {
    count++;
    n = &p->nodes[n->header.next];
  }
  return count;
}
 
static void smooth_paths_linsys(dt_iop_liquify_params_t *params)
{
  for(int k = 0; k < MAX_NODES; k++)
  {
    if(params->nodes[k].header.type == DT_LIQUIFY_PATH_INVALIDATED)
      break;
 
    if(params->nodes[k].header.prev != -1)
      continue;
 
    dt_liquify_path_data_t *node = &params->nodes[k];
 
    const size_t n = path_length(params, node);
 
    if(n < 2)
      continue;
 
    float complex *pt   = calloc(n, sizeof(float complex));
    float complex *c1   = calloc(n, sizeof(float complex));
    float complex *c2   = calloc(n, sizeof(float complex));
    int *eqn            = calloc(n, sizeof(int));
    size_t idx = 0;
 
    while(node)
    {
      const dt_liquify_path_data_t *d = (dt_liquify_path_data_t *) node;
      const dt_liquify_path_data_t *p = node_prev(params, node);
      const dt_liquify_path_data_t *n = node_next(params, node);
      const dt_liquify_path_data_t *nn = n ? node_next(params, n) : NULL;
 
      pt[idx] = node->warp.point;
      if(d->header.type == DT_LIQUIFY_PATH_CURVE_TO_V1)
      {
        c1[idx-1] = d->node.ctrl1;
        c2[idx-1] = d->node.ctrl2;
      }
 
      const int autosmooth      = d->header.node_type == DT_LIQUIFY_NODE_TYPE_AUTOSMOOTH;
      const int next_autosmooth = n   &&  n->header.node_type == DT_LIQUIFY_NODE_TYPE_AUTOSMOOTH;
      const int firstseg        = !p  ||  d->header.type != DT_LIQUIFY_PATH_CURVE_TO_V1;
      const int lastseg         = !nn ||  nn->header.type != DT_LIQUIFY_PATH_CURVE_TO_V1;
      const int lineseg         = n   &&  n->header.type == DT_LIQUIFY_PATH_LINE_TO_V1;
 
      // Program the linear system with equations:
      //
      //    START           END
      //    --------------------------
      // 1: straight        smooth
      // 2: smooth          smooth
      // 3: smooth          straight
      // 4: keep            smooth
      // 5: keep            keep
      // 6: smooth          keep
      // 7: keep            straight
      // 8: straight        straight   (== line)
      // 9: straight        keep
 
      if(lineseg)                                                    eqn[idx] = 5;
      else if(!autosmooth && !next_autosmooth)                       eqn[idx] = 5;
      else if(firstseg && lastseg && !autosmooth && next_autosmooth) eqn[idx] = 7;
      else if(firstseg && lastseg && autosmooth && next_autosmooth)  eqn[idx] = 8;
      else if(firstseg && lastseg && autosmooth && !next_autosmooth) eqn[idx] = 9;
      else if(firstseg && autosmooth && !next_autosmooth)            eqn[idx] = 5;
      else if(firstseg && autosmooth)                                eqn[idx] = 1;
      else if(lastseg && autosmooth && next_autosmooth)              eqn[idx] = 3;
      else if(lastseg && !autosmooth && next_autosmooth)             eqn[idx] = 7;
      else if(autosmooth && !next_autosmooth)                        eqn[idx] = 6;
      else if(!autosmooth && next_autosmooth)                        eqn[idx] = 4;
      else                                                            eqn[idx] = 2;
 
      ++idx;
      node = node_next(params, node);
    }
 
    smooth_path_linsys(n, pt, c1, c2, eqn);
 
    // write calculated control points back to list structure
    node = &params->nodes[k];
    node = node_next(params, node);
    idx = 0;
    while(node)
    {
      dt_liquify_path_data_t *d  = (dt_liquify_path_data_t *) node;
      if(d->header.type == DT_LIQUIFY_PATH_CURVE_TO_V1)
      {
        d->node.ctrl1 = c1[idx];
        d->node.ctrl2 = c2[idx];
      }
      ++idx;
      node = node_next(params, node);
    }
 
    dt_free(pt);
    dt_free(c1);
    dt_free(c2);
    dt_free(eqn);
  }
}
 
static dt_liquify_path_data_t *find_hovered(dt_iop_liquify_params_t *p)
{
  for(int k=0; k<MAX_NODES; k++)
    if(p->nodes[k].header.type == DT_LIQUIFY_PATH_INVALIDATED)
      break;
    else if(p->nodes[k].header.hovered)
      return &p->nodes[k];
  return NULL;
}
 
static void init_warp(dt_liquify_warp_t *warp, float complex point)
{
  warp->type     = DT_LIQUIFY_WARP_TYPE_LINEAR;
  warp->point    = point;
  warp->radius   = point;
  warp->strength = point;
  warp->control1 = 0.5;
  warp->control2 = 0.75;
  warp->status   = DT_LIQUIFY_STATUS_NONE;
}
 
static dt_liquify_path_data_t *alloc_move_to(dt_iop_module_t *module, float complex start_point)
{
  dt_iop_liquify_gui_data_t *g = (dt_iop_liquify_gui_data_t *)dt_iop_gui_data(module);
  dt_liquify_path_data_t* m = (dt_liquify_path_data_t*)node_alloc(&g->params, &g->node_index);
  if(m)
  {
    m->header.type = DT_LIQUIFY_PATH_MOVE_TO_V1;
    m->header.node_type = DT_LIQUIFY_NODE_TYPE_AUTOSMOOTH;
    init_warp(&m->warp, start_point);
  }
  return (dt_liquify_path_data_t *)m;
}
 
static dt_liquify_path_data_t *alloc_line_to(dt_iop_module_t *module, float complex end_point)
{
  dt_iop_liquify_gui_data_t *g = (dt_iop_liquify_gui_data_t *)dt_iop_gui_data(module);
  dt_liquify_path_data_t* l = (dt_liquify_path_data_t*)node_alloc(&g->params, &g->node_index);
  if(!IS_NULL_PTR(l))
  {
    l->header.type = DT_LIQUIFY_PATH_LINE_TO_V1;
    l->header.node_type = DT_LIQUIFY_NODE_TYPE_AUTOSMOOTH;
    init_warp(&l->warp, end_point);
  }
  return (dt_liquify_path_data_t *)l;
}
 
static dt_liquify_path_data_t *alloc_curve_to(dt_iop_module_t *module, float complex end_point)
{
  dt_iop_liquify_gui_data_t *g = (dt_iop_liquify_gui_data_t *)dt_iop_gui_data(module);
  dt_liquify_path_data_t* c = (dt_liquify_path_data_t*)node_alloc(&g->params, &g->node_index);
  if(!IS_NULL_PTR(c))
  {
    c->header.type = DT_LIQUIFY_PATH_CURVE_TO_V1;
    c->header.node_type = DT_LIQUIFY_NODE_TYPE_AUTOSMOOTH;
    c->node.ctrl1 = c->node.ctrl2 = 0.0;
    init_warp(&c->warp, end_point);
  }
  return (dt_liquify_path_data_t *)c;
}
 
static void unselect_all(dt_iop_liquify_params_t *p)
{
  for(int k=0; k<MAX_NODES; k++)
    if(p->nodes[k].header.type == DT_LIQUIFY_PATH_INVALIDATED)
      break;
    else
      p->nodes[k].header.selected = 0;
}
 
static float get_zoom_scale(const dt_develop_t *develop)
{
  return dt_dev_get_zoom_scale(develop,  1);
}
 
void gui_post_expose(struct dt_iop_module_t *module,
                      cairo_t *cr,
                      int32_t width,
                      int32_t height,
                      int32_t pointerx,
                      int32_t pointery)
{
  if(IS_NULL_PTR(module))
    return;
  dt_develop_t *develop = module->dev;
  if(IS_NULL_PTR(develop))
    return;
  dt_iop_liquify_gui_data_t *g = (dt_iop_liquify_gui_data_t *)dt_iop_gui_data(module);
  if(IS_NULL_PTR(g))
    return;
 
  const dt_dev_image_geometry_t geometry = dt_dev_geometry_snapshot(develop);
  const float bb_width = geometry.processed_width;
  const float bb_height = geometry.processed_height;
  if(bb_width < 1.0 || bb_height < 1.0)
    return;
 
  // get a copy of all iop params
  dt_iop_gui_enter_critical_section(module);
  update_warp_count(g);
  smooth_paths_linsys(&g->params);
  dt_iop_liquify_params_t copy_params;
  memcpy(&copy_params, &g->params, sizeof(dt_iop_liquify_params_t));
  dt_iop_gui_leave_critical_section(module);
 
  /* Distort all points, through the geometry service. Nothing is resynchronised first: the
   * chain is rebuilt wherever a pipe flag is raised, and this fold excludes the module itself
   * (BACK_EXCL then FORW_EXCL), so the live params being dragged are not what it reads. If the
   * chain cannot answer, its compose leaves the points in RAW coordinates -- which would draw
   * the whole path in the wrong place -- so draw nothing instead. */
  if(!dt_geometry_chain_authoritative(develop->geometry_chain))
    return;
 
  const distort_params_t d_params = { develop, NULL, develop->geometry_chain, 1.0, 1.0,
                                      DT_DEV_TRANSFORM_DIR_ALL, FALSE };
  _distort_paths(module, &d_params, &copy_params);
 
  // You're not supposed to understand this
  const float zoom_scale = get_zoom_scale(develop) * dt_dev_viewport_scaling(develop);
 
  if(dt_dev_rescale_roi_to_input(develop, cr, width, height))
    return;
 
  draw_paths(module, cr, 1.0 / zoom_scale, &copy_params);
}
 
static gboolean btn_make_radio_callback(GtkToggleButton *btn, GdkEventButton *event, dt_iop_module_t *module);
 
void gui_focus(struct dt_iop_module_t *module, gboolean in)
{
  if(!in)
  {
    dt_collection_hint_message(dt_collection_get_global());
    btn_make_radio_callback(NULL, NULL, module);
  }
}
 
static void sync_pipe(struct dt_iop_module_t *module, gboolean history)
{
  if(history)
  {
    const dt_iop_liquify_gui_data_t *g = (dt_iop_liquify_gui_data_t *)dt_iop_gui_data(module);
    // something definitive has happened like button release ... so
    // redraw pipe
    memcpy(module->params, &g->params, sizeof(dt_iop_liquify_params_t));
    dt_dev_add_history_item(module->dev, module, TRUE, TRUE);
  }
  else
  {
    // only moving mouse around, pointing at things or dragging ... so
    // give some cairo feedback, but don't redraw pipe
    dt_control_queue_redraw_center();
  }
}
 
/*
  right-click on node:       Delete node.
  right-click on path:       Delete whole path.
 
  ctrl+click on node:        Cycle symmetrical, smooth, cusp, autosmooth
  ctrl+click on path:        Add node
  ctrl+alt+click on path:    Change line / bezier
 
  ctrl+click on strength:    Cycle linear, grow, shrink
*/
 
static void get_point_scale(struct dt_iop_module_t *module, float x, float y, float complex *pt, float *scale)
{
  float pts[2] = { (float)x, (float)y };
  dt_dev_coordinates_widget_to_image_norm(module->dev, pts, 1);
  dt_dev_coordinates_image_norm_to_image_abs(module->dev, pts, 1);
  dt_dev_distort_backtransform_gui(module->dev,
                                    module->iop_order,DT_DEV_TRANSFORM_DIR_FORW_EXCL, pts, 1);
  dt_dev_distort_backtransform_gui(module->dev,
                                    module->iop_order,DT_DEV_TRANSFORM_DIR_BACK_EXCL, pts, 1);
 
  *scale = get_zoom_scale(module->dev);
  *pt = pts[0] + pts[1] * I;
}
 
static gboolean _is_movable_layer(const dt_liquify_layer_enum_t layer)
{
  switch(layer)
  {
    case DT_LIQUIFY_LAYER_PATH:
    case DT_LIQUIFY_LAYER_CENTERPOINT:
    case DT_LIQUIFY_LAYER_CTRLPOINT1:
    case DT_LIQUIFY_LAYER_CTRLPOINT2:
    case DT_LIQUIFY_LAYER_RADIUSPOINT:
    case DT_LIQUIFY_LAYER_HARDNESSPOINT1:
    case DT_LIQUIFY_LAYER_HARDNESSPOINT2:
    case DT_LIQUIFY_LAYER_STRENGTHPOINT:
      return TRUE;
    default:
      return FALSE;
  }
}
 
int mouse_moved(struct dt_iop_module_t *module,
                 double x,
                 double y,
                 double pressure,
                 int which)
{
  dt_iop_liquify_gui_data_t *g = (dt_iop_liquify_gui_data_t *)dt_iop_gui_data(module);
  gboolean handled = FALSE;
  float complex pt = 0.0f;
  float scale = 0.0f;
 
  get_point_scale(module, x, y, &pt, &scale);
 
  dt_iop_gui_enter_critical_section(module);
 
  const float complex prev_mouse_pos = g->last_mouse_pos;
  g->last_mouse_pos = pt;
 
  // Don't hit test while dragging, you'd only hit the dragged thing
  // anyway.
 
  if(!is_dragging(g))
  {
    dt_liquify_hit_t hit = _hit_test_paths(module, &g->params, pt);
    // Only override the cursor for a movable layer -- otherwise leave darkroom's own default
    // cursor (dot/crosshair/left_ptr, already queued per position by
    // _darkroom_set_default_cursor before this handler runs) instead of forcing "default".
    if(_is_movable_layer(hit.layer))
      dt_control_queue_cursor_by_name("move");
    dt_liquify_path_data_t *last_hovered = find_hovered(&g->params);
    if(hit.elem != last_hovered
       || (!IS_NULL_PTR(last_hovered) && !IS_NULL_PTR(hit.elem)
           && hit.elem->header.hovered != last_hovered->header.hovered))
    {
      if(!IS_NULL_PTR(hit.elem))
        hit.elem->header.hovered = hit.layer;
      if(!IS_NULL_PTR(last_hovered))
        last_hovered->header.hovered = 0;
      // change in hover display
      dt_control_hinter_message(dt_control_get_global(), dt_liquify_layers[hit.layer].hint);
      handled = TRUE;
      goto done;
    }
 
    const gboolean dragged = detect_drag(g, scale, pt);
 
    if(dragged && !IS_NULL_PTR(g->last_hit.elem))
    {
      // start dragging
      start_drag(g, g->last_hit.layer, g->last_hit.elem);
      // nothing more to do, we will refresh on the next call anyway
      // this makes the initial move of a node a bit more fluid.
      handled = TRUE;
      goto done;
    }
 
    if(!IS_NULL_PTR(g->last_hit.elem))
    {
      // an item is selected, so this mouvement is handled and must
      // not trigger any panning.
      handled = TRUE;
    }
  }
  else // we are dragging
  {
    dt_control_queue_cursor_by_name("move");
    dt_liquify_path_data_t *d = g->dragging.elem;
    dt_liquify_path_data_t *n = node_next(&g->params, d);
    dt_liquify_path_data_t *p = node_prev(&g->params, d);
 
    const float complex *start_pt = &d->warp.point;
 
    switch (g->dragging.layer)
    {
       case DT_LIQUIFY_LAYER_PATH:
       {
         if(IS_NULL_PTR(p)) break;
         const float complex delta = prev_mouse_pos == -1 ? 0.0f : pt - prev_mouse_pos;
         if(cabsf(delta) < 1e-6f) break;
 
         dt_liquify_path_data_t *moved_nodes[2] = { p, d };
         for(int i = 0; i < 2; i++)
         {
           dt_liquify_path_data_t *node = moved_nodes[i];
           dt_liquify_path_data_t *node_next_ptr = node_next(&g->params, node);
           dt_liquify_path_data_t *node_prev_ptr = node_prev(&g->params, node);
 
           if(node->header.type == DT_LIQUIFY_PATH_CURVE_TO_V1)
             node->node.ctrl2 += delta;
           if(!IS_NULL_PTR(node_next_ptr) && node_next_ptr->header.type == DT_LIQUIFY_PATH_CURVE_TO_V1)
             node_next_ptr->node.ctrl1 += delta;
           if(!IS_NULL_PTR(node_prev_ptr) && node_prev_ptr->header.type == DT_LIQUIFY_PATH_CURVE_TO_V1)
             node_prev_ptr->node.ctrl2 += delta;
 
           node->warp.radius += delta;
           node->warp.strength += delta;
           node->warp.point += delta;
         }
         break;
       }
 
       case DT_LIQUIFY_LAYER_CENTERPOINT:
         switch (d->header.type)
         {
            case DT_LIQUIFY_PATH_CURVE_TO_V1:
              d->node.ctrl2 += pt - d->warp.point;
              // fall thru
            case DT_LIQUIFY_PATH_MOVE_TO_V1:
            case DT_LIQUIFY_PATH_LINE_TO_V1:
              if(!IS_NULL_PTR(n) && n->header.type == DT_LIQUIFY_PATH_CURVE_TO_V1)
                n->node.ctrl1 += pt - d->warp.point;
              if(!IS_NULL_PTR(p) && p->header.type == DT_LIQUIFY_PATH_CURVE_TO_V1)
                p->node.ctrl2 += pt - d->warp.point;
              d->warp.radius   += pt - d->warp.point;
              d->warp.strength += pt - d->warp.point;
              d->warp.point = pt;
              break;
            default:
              break;
         }
         break;
 
       case DT_LIQUIFY_LAYER_CTRLPOINT1:
         switch (d->header.type)
         {
            case DT_LIQUIFY_PATH_CURVE_TO_V1:
              d->node.ctrl1 = pt;
              if(!IS_NULL_PTR(p) && p->header.type == DT_LIQUIFY_PATH_CURVE_TO_V1)
              {
                switch (p->header.node_type)
                {
                   case DT_LIQUIFY_NODE_TYPE_SMOOTH:
                     p->node.ctrl2 = p->warp.point +
                       cabsf(p->warp.point - p->node.ctrl2) *
                       cexpf(cargf(p->warp.point - pt) * I);
                     break;
                case DT_LIQUIFY_NODE_TYPE_SYMMETRICAL:
                  p->node.ctrl2 = 2 * p->warp.point - pt;
                  break;
                default:
                  break;
                }
              }
              break;
            default:
              break;
         }
         break;
 
       case DT_LIQUIFY_LAYER_CTRLPOINT2:
         switch (d->header.type)
         {
            case DT_LIQUIFY_PATH_CURVE_TO_V1:
              d->node.ctrl2 = pt;
              if(!IS_NULL_PTR(n) && n->header.type == DT_LIQUIFY_PATH_CURVE_TO_V1)
              {
                switch (d->header.node_type)
                {
                   case DT_LIQUIFY_NODE_TYPE_SMOOTH:
                     n->node.ctrl1 = d->warp.point +
                       cabsf(d->warp.point - n->node.ctrl1) *
                       cexpf(cargf(d->warp.point - pt) * I);
                     break;
                   case DT_LIQUIFY_NODE_TYPE_SYMMETRICAL:
                     n->node.ctrl1 = 2 * d->warp.point - pt;
                     break;
                   default:
                     break;
                }
              }
              break;
            default:
              break;
         }
         break;
 
       case DT_LIQUIFY_LAYER_RADIUSPOINT:
         d->warp.radius = pt;
         dt_conf_set_float(CONF_RADIUS, cabsf(d->warp.radius - d->warp.point));
         break;
 
       case DT_LIQUIFY_LAYER_STRENGTHPOINT:
         d->warp.strength = pt;
         dt_conf_set_float(CONF_STRENGTH, cabsf(d->warp.strength - d->warp.point));
         dt_conf_set_float(CONF_ANGLE, cargf(d->warp.strength - d->warp.point));
         break;
 
       case DT_LIQUIFY_LAYER_HARDNESSPOINT1:
         d->warp.control1 = MIN(1.0, cabsf(pt - *start_pt) / cabsf(d->warp.radius - *start_pt));
         break;
 
       case DT_LIQUIFY_LAYER_HARDNESSPOINT2:
         d->warp.control2 = MIN(1.0, cabsf(pt - *start_pt) / cabsf(d->warp.radius - *start_pt));
         break;
 
       default:
         break;
    }
    handled = TRUE;
  }
 
done:
  if(!handled && which != 0 && (!IS_NULL_PTR(g->temp) || !IS_NULL_PTR(g->last_hit.elem)))
  {
    // A drag/edit sequence is in progress, keep consuming motion events
    // so darkroom pan does not steal the interaction.
    handled = TRUE;
  }
  dt_iop_gui_leave_critical_section(module);
  if(handled)
  {
    sync_pipe(module, FALSE);
  }
  return handled;
}
 
static float dt_conf_get_sanitize_float(const char *name, float min, float max, float default_value)
{
  const float value = dt_conf_get_float(name);
  float new_value = CLAMP(value, min, max);
 
  if (default_value != 0.0f && new_value != value) new_value = 0.25f * default_value + 0.75f * value;
 
  dt_conf_set_float(name, new_value);
  return new_value;
}
 
static void get_stamp_params(dt_iop_module_t *module, float *radius, float *r_strength, float *phi)
{
  GtkWidget *widget = dt_gui_main_window();
  GtkAllocation allocation;
  gtk_widget_get_allocation(widget, &allocation);
  const int last_win_min = MIN(allocation.width, allocation.height);
 
  const dt_develop_t *dev = module->dev;
  const float iwd_min = MIN(dt_dev_geometry_raw_width(dev), dt_dev_geometry_raw_height(dev));
  const float proc_wdht_min = MIN(dt_dev_geometry_processed_width(dev), dt_dev_geometry_processed_height(dev));
  const float scale = 1.f / (get_zoom_scale(dev));
  const float im_scale = 0.09f * iwd_min * last_win_min * scale / proc_wdht_min;
 
  *radius = dt_conf_get_sanitize_float(CONF_RADIUS, 0.1f*im_scale, 3.0f*im_scale, im_scale);
  *r_strength = dt_conf_get_sanitize_float(CONF_STRENGTH, 0.5f * *radius, 2.0f * *radius, 1.5f * *radius);
  *phi = dt_conf_get_sanitize_float(CONF_ANGLE, -M_PI, M_PI, 0.0f);
}
/*
  add support for changing the radius and the strength vector for the temp node
 */
int scrolled(struct dt_iop_module_t *module, double x, double y, int up, uint32_t state)
{
  const dt_iop_liquify_gui_data_t *g = (dt_iop_liquify_gui_data_t *)dt_iop_gui_data(module);
 
  // add an option to allow skip mouse events while editing masks
  const gboolean incr = dt_mask_scroll_increases(up);
 
  if(g->temp)
  {
    dt_liquify_warp_t *warp = &g->temp->warp;
    const float complex strength_v = warp->strength - warp->point;
    if(dt_modifier_is(state, 0))
    {
      //  change size
      float radius = 0.0f, r = 0.0f, phi = 0.0f;
      get_stamp_params(module, &radius, &r, &phi);
 
      float factor = 1.0f;
      if(incr)
        factor *= 1.0f / 0.97f;
      else if(!incr && cabsf(warp->radius - warp->point) > 10.0f)
        factor *= 0.97f;
 
      r *= factor;
      radius *= factor;
 
      warp->radius = warp->point + (radius * factor);
      warp->strength = warp->point + r * cexpf(phi * I);
 
      dt_conf_set_float(CONF_RADIUS, radius);
      dt_conf_set_float(CONF_STRENGTH, r);
      return 1;
    }
    else if(dt_modifier_is(state, DT_PRIMARY_MASK))
    {
      //  change the strength direction
      float phi = cargf(strength_v);
      const float r = cabsf(strength_v);
 
      if(incr)
        phi += DT_M_PI_F / 16.0f;
      else
        phi -= DT_M_PI_F / 16.0f;
 
      warp->strength = warp->point + r * cexpf(phi * I);
      dt_conf_set_float(CONF_STRENGTH, r);
      dt_conf_set_float(CONF_ANGLE, phi);
      return 1;
    }
    else if(dt_modifier_is(state, GDK_SHIFT_MASK))
    {
      //  change the strength
      const float phi = cargf(strength_v);
      float r = cabsf(strength_v);
 
      if(incr)
        r *= 1.0f / 0.97f;
      else
        r *= 0.97f;
 
      warp->strength = warp->point + r * cexpf(phi * I);
      dt_conf_set_float(CONF_STRENGTH, r);
      dt_conf_set_float(CONF_ANGLE, phi);
      return 1;
    }
  }
 
  return 0;
}
 
int button_pressed(struct dt_iop_module_t *module,
                    double x,
                    double y,
                    double pressure,
                    int which,
                    int type,
                    uint32_t state)
{
  dt_iop_liquify_gui_data_t *g = (dt_iop_liquify_gui_data_t *)dt_iop_gui_data(module);
  int handled = 0;
  float complex pt = 0.0f;
  float scale = 0.0f;
 
  get_point_scale(module, x, y, &pt, &scale);
 
  dt_iop_gui_enter_critical_section(module);
 
  g->last_mouse_pos = pt;
  g->last_mouse_mods = state;
  if(which == 1)
    g->last_button1_pressed_pos = pt;
 
  if(!is_dragging(g))
    // while dragging you would always hit the dragged thing
    g->last_hit = _hit_test_paths(module, &g->params, pt);
 
  if(which == 2) goto done;
 
  // Point tool
 
  if(which == 1 && gtk_toggle_button_get_active(g->btn_point_tool))
  {
    // always end dragging before manipulating the path list to avoid
    // dangling pointers
    end_drag(g);
 
    if(IS_NULL_PTR(g->temp)) goto done;
    g->status |= DT_LIQUIFY_STATUS_NEW;
    g->status &= ~DT_LIQUIFY_STATUS_PREVIEW;
 
    start_drag(g, DT_LIQUIFY_LAYER_STRENGTHPOINT, g->temp);
    g->last_hit = NOWHERE;
    handled = 1;
    goto done;
  }
 
  // Line tool or curve tool
 
  if(which == 1 && (gtk_toggle_button_get_active(g->btn_line_tool)
                     || gtk_toggle_button_get_active(g->btn_curve_tool)))
  {
    // always end dragging before manipulating the path list to avoid
    // dangling pointers
    end_drag(g);
    if(IS_NULL_PTR(g->temp))
    {
      if(g->last_hit.layer == DT_LIQUIFY_LAYER_CENTERPOINT)
      {
        // continue path
        g->temp = g->last_hit.elem;
      }
      else
      {
        goto done;
      }
    }
 
    g->last_hit = NOWHERE;
    if(gtk_toggle_button_get_active(g->btn_curve_tool))
    {
      start_drag(g, DT_LIQUIFY_LAYER_CTRLPOINT1, g->temp);
    }
    g->status |= DT_LIQUIFY_STATUS_NEW;
    g->status &= ~DT_LIQUIFY_STATUS_PREVIEW;
    handled = 1;
    goto done;
  }
 
  // Node tool
 
  if(gtk_toggle_button_get_active(g->btn_node_tool))
  {
    if(which == 1 && dt_modifier_is(g->last_mouse_mods, DT_PRIMARY_MASK) &&
        (g->last_hit.layer == DT_LIQUIFY_LAYER_CENTERPOINT))
    {
      // cycle node type: smooth -> cusp etc.
      dt_liquify_path_data_t *node = g->last_hit.elem;
      node->header.node_type = (node->header.node_type + 1) % DT_LIQUIFY_NODE_TYPE_LAST;
      handled = 1;
      goto done;
    }
    if(which == 1 && dt_modifier_is(g->last_mouse_mods, DT_PRIMARY_MASK) &&
        (g->last_hit.layer == DT_LIQUIFY_LAYER_STRENGTHPOINT))
    {
      // cycle warp type: linear -> radial etc.
      if(g->last_hit.elem->header.type == DT_LIQUIFY_PATH_MOVE_TO_V1)
      {
        dt_liquify_warp_t *warp = &g->last_hit.elem->warp;
        warp->type = (warp->type + 1) % DT_LIQUIFY_WARP_TYPE_LAST;
      }
      handled = 1;
      goto done;
    }
  }
 
  if(!handled && (which == 1 || which == 3) && (!IS_NULL_PTR(g->temp) || !IS_NULL_PTR(g->last_hit.elem)))
  {
    // Even when the actual edit is finalized on button release, this press
    // starts an interaction sequence and must remain captured by liquify.
    handled = 1;
  }
 
done:
  dt_iop_gui_leave_critical_section(module);
  return handled;
}
 
static void _start_new_shape(dt_iop_module_t *module)
{
  dt_iop_liquify_gui_data_t *g = (dt_iop_liquify_gui_data_t *)dt_iop_gui_data(module);
 
  //  create initial shape at the center
  float complex pt = 0.0f;
  float scale = 1.0f;
  get_point_scale(module, 0.5f * dt_dev_viewport_box_width(module->dev), 0.5f * dt_dev_viewport_box_height(module->dev), &pt, &scale);
  float radius = 0.0f, r = 1.0f, phi = 0.0f;
  get_stamp_params(module, &radius, &r, &phi);
  //  start a new path
  g->temp = alloc_move_to(module, pt);
  g->temp->warp.radius = pt + radius;
  g->temp->warp.strength = pt + r * cexpf(phi * I);
  g->status |= DT_LIQUIFY_STATUS_PREVIEW;
  g->status |= DT_LIQUIFY_STATUS_NEW;
 
  g->just_started = TRUE;
 
 
  start_drag(g, DT_LIQUIFY_LAYER_CENTERPOINT, g->temp);
  g->last_hit = NOWHERE;
}
 
int button_released(struct dt_iop_module_t *module,
                     double x,
                     double y,
                     int which,
                     uint32_t state)
{
  dt_iop_liquify_gui_data_t *g = (dt_iop_liquify_gui_data_t *)dt_iop_gui_data(module);
  int handled = 0;
  float complex pt = 0.0f;
  float scale = 0.0f;
 
  get_point_scale(module, x, y, &pt, &scale);
 
  dt_iop_gui_enter_critical_section(module);
 
  g->last_mouse_pos = pt;
 
  const gboolean dragged = detect_drag(g, scale, pt);
 
  if(which == 1 && g->temp && (g->status & DT_LIQUIFY_STATUS_NEW))
  {
    end_drag(g);
    if(gtk_toggle_button_get_active(g->btn_point_tool))
    {
      g->temp = NULL; // a point is done
      btn_make_radio_callback(g->btn_node_tool, NULL, module);
      handled = 2;
    }
    else if(gtk_toggle_button_get_active(g->btn_line_tool))
    {
      const int prev_index = g->node_index;
      const float complex strength = (g->temp->warp.strength - g->temp->warp.point);
      const float radius = cabsf(g->temp->warp.radius - g->temp->warp.point);
      g->temp = alloc_line_to(module, pt);
      if(IS_NULL_PTR(g->temp)) goto done;
      g->temp->warp.radius = pt + radius;
      g->temp->warp.strength = pt + strength;
      // links
      g->temp->header.prev = prev_index;
      node_get(&g->params, prev_index)->header.next = g->node_index;
      start_drag(g, DT_LIQUIFY_LAYER_CENTERPOINT, g->temp);
      g->just_started = FALSE;
      handled = 1;
    }
    else if(gtk_toggle_button_get_active(g->btn_curve_tool))
    {
      const int prev_index = g->node_index;
      const float complex strength = (g->temp->warp.strength - g->temp->warp.point);
      const float radius = cabsf(g->temp->warp.radius - g->temp->warp.point);
      g->temp = alloc_curve_to(module, pt);
      if(IS_NULL_PTR(g->temp)) goto done;
      g->temp->warp.radius = pt + radius;
      g->temp->warp.strength = pt + strength;
      // links
      g->temp->header.prev = prev_index;
      node_get(&g->params, prev_index)->header.next = g->node_index;
      start_drag(g, DT_LIQUIFY_LAYER_CENTERPOINT, g->temp);
      g->just_started = FALSE;
      handled = 1;
    }
    g->status &= ~DT_LIQUIFY_STATUS_NEW;
    goto done;
  }
 
  if(which == 1 && is_dragging(g))
  {
    end_drag(g);
    handled = 2;
    goto done;
  }
 
  // right click == cancel or delete
  if(which == 3)
  {
    dt_control_hinter_message(dt_control_get_global(), "");
    end_drag(g);
 
    // cancel line or curve creation
    if(g->temp)
    {
      node_delete(&g->params, g->temp);
      g->temp = NULL;
      g->status &= ~DT_LIQUIFY_STATUS_PREVIEW;
      btn_make_radio_callback(g->btn_node_tool, NULL, module);
      handled = 2;
      goto done;
    }
 
    // right click on background toggles node tool
    if(g->last_hit.layer == DT_LIQUIFY_LAYER_BACKGROUND)
    {
      btn_make_radio_callback(g->btn_node_tool, NULL, module);
      handled = 1;
      goto done;
    }
 
    // delete node
    if(g->last_hit.layer == DT_LIQUIFY_LAYER_CENTERPOINT)
    {
      node_delete(&g->params, g->last_hit.elem);
      g->last_hit = NOWHERE;
      handled = 2;
      goto done;
    }
    // delete shape
    if(g->last_hit.layer == DT_LIQUIFY_LAYER_PATH)
    {
      path_delete(&g->params, g->last_hit.elem);
      g->last_hit = NOWHERE;
      handled = 2;
      goto done;
    }
    goto done;
  }
 
  // Node tool
 
  if(gtk_toggle_button_get_active(g->btn_node_tool))
  {
    if(which == 1 && dt_modifier_is(g->last_mouse_mods, 0) && !dragged)
    {
      // select/unselect start/endpoint and clear previous selections
      if(g->last_hit.layer == DT_LIQUIFY_LAYER_CENTERPOINT)
      {
        const int oldsel = !!g->last_hit.elem->header.selected;
        unselect_all(&g->params);
        g->last_hit.elem->header.selected = oldsel ? 0 : g->last_hit.layer;
        handled = 1;
        goto done;
      }
      // unselect all
      if(g->last_hit.layer == DT_LIQUIFY_LAYER_BACKGROUND)
      {
        unselect_all(&g->params);
        handled = 1;
        goto done;
      }
    }
    if(which == 1 && dt_modifier_is(g->last_mouse_mods, GDK_SHIFT_MASK) && !dragged)
    {
      // select/unselect start/endpoint and keep previous selections
      if(g->last_hit.layer == DT_LIQUIFY_LAYER_CENTERPOINT)
      {
        const int oldsel = !!g->last_hit.elem->header.selected;
        g->last_hit.elem->header.selected = oldsel ? 0 : g->last_hit.layer;
        handled = 1;
        goto done;
      }
    }
    if(which == 1 && dt_modifier_is(g->last_mouse_mods, DT_PRIMARY_MASK) && !dragged)
    {
      // add node
      if(g->last_hit.layer == DT_LIQUIFY_LAYER_PATH)
      {
        dt_liquify_path_data_t *e = g->last_hit.elem;
        dt_liquify_path_data_t *prev = node_prev(&g->params, e);
        if(!IS_NULL_PTR(prev) && e->header.type == DT_LIQUIFY_PATH_CURVE_TO_V1)
        {
          // add node to curve
          dt_liquify_path_data_t *curve1 = (dt_liquify_path_data_t *) e;
 
          dt_liquify_path_data_t *curve2 = (dt_liquify_path_data_t *)alloc_curve_to(module, 0);
          if(IS_NULL_PTR(curve2)) goto done;
 
          curve2->node.ctrl1 = curve1->node.ctrl1;
          curve2->node.ctrl2 = curve1->node.ctrl2;
 
          dt_liquify_warp_t *warp1 = &prev->warp;
          dt_liquify_warp_t *warp2 = &curve2->warp;
          dt_liquify_warp_t *warp3 = &e->warp;
 
          const float t = find_nearest_on_curve_t(warp1->point, curve1->node.ctrl1, curve1->node.ctrl2,
                                                   warp3->point, pt, INTERPOLATION_POINTS);
 
          float complex midpoint = warp3->point;
          casteljau(&warp1->point, &curve1->node.ctrl1, &curve1->node.ctrl2, &midpoint, t);
          midpoint = warp1->point;
          casteljau(&warp3->point, &curve2->node.ctrl2, &curve2->node.ctrl1, &midpoint, 1.0 - t);
 
          mix_warps(warp2, warp1, warp3, midpoint, t);
 
          node_insert_before(&g->params, e, (dt_liquify_path_data_t *)curve2);
 
          handled = 2;
          goto done;
        }
        if(!IS_NULL_PTR(prev) && e->header.type == DT_LIQUIFY_PATH_LINE_TO_V1)
        {
          // add node to line
          dt_liquify_warp_t *warp1 = &prev->warp;
          dt_liquify_warp_t *warp3 = &e->warp;
          const float t = find_nearest_on_line_t(warp1->point, warp3->point, pt);
 
          dt_liquify_path_data_t *tmp = alloc_line_to(module, e->warp.point);
          if(IS_NULL_PTR(tmp)) goto done;
 
          dt_liquify_warp_t *warp2 = &tmp->warp;
          const float complex midpoint = cmix(warp1->point, warp3->point, t);
 
          mix_warps(warp2, warp1, warp3, midpoint, t);
          node_insert_before(&g->params, e, tmp);
 
          handled = 2;
          goto done;
        }
      }
    }
    if(which == 1
       && dt_modifier_is(g->last_mouse_mods, GDK_MOD1_MASK | DT_PRIMARY_MASK)
       && !dragged)
    {
      if(g->last_hit.layer == DT_LIQUIFY_LAYER_PATH)
      {
        // change segment
        dt_liquify_path_data_t *e = g->last_hit.elem;
        dt_liquify_path_data_t *prev = node_prev(&g->params, e);
        if(!IS_NULL_PTR(prev) && e->header.type == DT_LIQUIFY_PATH_CURVE_TO_V1)
        {
          // curve -> line
          e->header.type = DT_LIQUIFY_PATH_LINE_TO_V1;
          e->header.node_type = DT_LIQUIFY_NODE_TYPE_AUTOSMOOTH;
          e->header.selected = e->header.hovered = 0;
          handled = 2;
          goto done;
        }
        if(!IS_NULL_PTR(prev) && e->header.type == DT_LIQUIFY_PATH_LINE_TO_V1)
        {
          // line -> curve
          const float complex p0 = prev->warp.point;
          const float complex p1 = e->warp.point;
          dt_liquify_path_data_t *c = (dt_liquify_path_data_t *)e;
          e->header.type = DT_LIQUIFY_PATH_CURVE_TO_V1;
          e->header.node_type = DT_LIQUIFY_NODE_TYPE_AUTOSMOOTH;
          c->node.ctrl1 = (2 * p0 +     p1) / 3.0;
          c->node.ctrl2 = (    p0 + 2 * p1) / 3.0;
 
          handled = 2;
          goto done;
        }
      }
    }
  }
 
done:
  if(which == 1)
    g->last_button1_pressed_pos = -1;
  g->last_hit = NOWHERE;
  dt_iop_gui_leave_critical_section(module);
  if(handled)
  {
    update_warp_count(g);
    sync_pipe(module, handled == 2);
  }
  return handled;
}
 
int key_pressed(struct dt_iop_module_t *self, GdkEventKey *event)
{
  if(IS_NULL_PTR(event)) return 0;
 
  dt_iop_liquify_gui_data_t *g = (dt_iop_liquify_gui_data_t *)dt_iop_gui_data(self);
  if(IS_NULL_PTR(g)) return 0;
  guint key = dt_keys_mainpad_alternatives(event->keyval);
 
 
  const gboolean creating = gtk_toggle_button_get_active(g->btn_point_tool)
                              || gtk_toggle_button_get_active(g->btn_line_tool)
                              || gtk_toggle_button_get_active(g->btn_curve_tool)
                              || (g->status & (DT_LIQUIFY_STATUS_NEW | DT_LIQUIFY_STATUS_PREVIEW));
 
  // Delete last created node while creating a shape.
  if(key == GDK_KEY_BackSpace)
  {
    if(!creating)
      return 0;
 
    const gboolean create_tool_active = gtk_toggle_button_get_active(g->btn_point_tool)
                                        || gtk_toggle_button_get_active(g->btn_line_tool)
                                        || gtk_toggle_button_get_active(g->btn_curve_tool);
    gboolean restart_shape = FALSE;
 
    dt_iop_gui_enter_critical_section(self);
 
    dt_liquify_path_data_t *last = NULL;
    for(int k = 0; k < MAX_NODES; k++)
    {
      if(g->params.nodes[k].header.type == DT_LIQUIFY_PATH_INVALIDATED)
        break;
      last = &g->params.nodes[k];
    }
 
    if(IS_NULL_PTR(last) && IS_NULL_PTR(g->temp))
    {
      restart_shape = create_tool_active;
      dt_iop_gui_leave_critical_section(self);
      if(restart_shape)
      {
        _start_new_shape(self);
        update_warp_count(g);
        sync_pipe(self, TRUE);
      }
      return 1;
    }
 
    dt_liquify_path_data_t *to_delete = !IS_NULL_PTR(g->temp) ? g->temp : last;
 
    end_drag(g);
    const int prev_index = to_delete->header.prev;
    node_delete(&g->params, to_delete);
    g->temp = prev_index >= 0 ? node_get(&g->params, prev_index) : NULL;
    g->node_index = !IS_NULL_PTR(g->temp) ? g->temp->header.idx : 0;
    g->last_hit = NOWHERE;
 
    if(!IS_NULL_PTR(g->temp))
    {
      // Continue creation from the current endpoint by restoring the live drag preview.
      start_drag(g, DT_LIQUIFY_LAYER_CENTERPOINT, g->temp);
      g->status &= ~DT_LIQUIFY_STATUS_NEW;
    }
    else
    {
      restart_shape = create_tool_active;
      g->status &= ~(DT_LIQUIFY_STATUS_PREVIEW | DT_LIQUIFY_STATUS_NEW);
    }
 
    dt_iop_gui_leave_critical_section(self);
 
    if(restart_shape)
      _start_new_shape(self);
 
    update_warp_count(g);
    sync_pipe(self, TRUE);
    return 1;
  }
 
  // Delete selected node outside creation mode.
  if(key == GDK_KEY_Delete)
  {
    if(creating)
      return 0;
 
    dt_iop_gui_enter_critical_section(self);
 
    dt_liquify_path_data_t *selected = NULL;
    for(int k = 0; k < MAX_NODES; k++)
    {
      if(g->params.nodes[k].header.type == DT_LIQUIFY_PATH_INVALIDATED)
        break;
      if(g->params.nodes[k].header.selected == DT_LIQUIFY_LAYER_CENTERPOINT)
      {
        selected = &g->params.nodes[k];
        break;
      }
    }
 
    if(IS_NULL_PTR(selected))
    {
      dt_iop_gui_leave_critical_section(self);
      return 0;
    }
 
    end_drag(g);
    const int deleted_idx = selected->header.idx;
    const int next_idx = selected->header.next;
    const int prev_idx = selected->header.prev;
    node_delete(&g->params, selected);
    g->temp = NULL;
    g->last_hit = NOWHERE;
 
    unselect_all(&g->params);
    int target_idx = (next_idx != -1) ? next_idx : prev_idx;
    if(target_idx > deleted_idx)
      target_idx--;
 
    if(target_idx >= 0)
    {
      dt_liquify_path_data_t *target = node_get(&g->params, target_idx);
      if(!IS_NULL_PTR(target) && target->header.type != DT_LIQUIFY_PATH_INVALIDATED)
        target->header.selected = DT_LIQUIFY_LAYER_CENTERPOINT;
    }
 
    dt_iop_gui_leave_critical_section(self);
 
    update_warp_count(g);
    sync_pipe(self, TRUE);
    return 1;
  }
 
  // Quit current creation or edition on Escape or Enter key
  if(key == GDK_KEY_Escape || key == GDK_KEY_Return)
  {
 
    if(!creating)
      return 0;
 
    dt_iop_gui_enter_critical_section(self);
 
    dt_control_hinter_message(dt_control_get_global(), "");
    end_drag(g);
 
    if(g->temp)
    {
      node_delete(&g->params, g->temp);
      g->temp = NULL;
    }
 
    g->status &= ~(DT_LIQUIFY_STATUS_PREVIEW | DT_LIQUIFY_STATUS_NEW);
    g->last_hit = NOWHERE;
 
    dt_iop_gui_leave_critical_section(self);
 
    gtk_toggle_button_set_active(g->btn_point_tool, FALSE);
    gtk_toggle_button_set_active(g->btn_line_tool, FALSE);
    gtk_toggle_button_set_active(g->btn_curve_tool, FALSE);
    gtk_toggle_button_set_active(g->btn_node_tool, TRUE);
 
    dt_control_hinter_message(dt_control_get_global(), _("click to edit nodes"));
    dt_iop_request_focus(self);
    update_warp_count(g);
    sync_pipe(self, TRUE);
 
    return 1;
  }
 
  return 0;
}
 
// we need this only because darktable has no radiobutton support
 
static gboolean btn_make_radio_callback(GtkToggleButton *btn, GdkEventButton *event, dt_iop_module_t *module)
{
  dt_iop_liquify_gui_data_t *g = (dt_iop_liquify_gui_data_t *)dt_iop_gui_data(module);
 
  // if currently dragging and a form (line or node) has been started, does nothing (expect resetting the toggle button status).
  if(is_dragging(g) && g->temp && node_prev(&g->params, g->temp))
  {
    return TRUE;
  }
 
  dt_control_hinter_message(dt_control_get_global(), "");
 
  // if we are on a preview, it means that a form (point, line, curve) has been started, but no node has yet been placed.
  // in this case we abort the current preview and let the new tool activated.
  if(g->status & DT_LIQUIFY_STATUS_PREVIEW)
  {
    node_delete(&g->params, g->temp);
    g->temp = NULL;
    g->status &= ~DT_LIQUIFY_STATUS_PREVIEW;
  }
 
  // now, let's enable and start a new form safely
  if(IS_NULL_PTR(btn) || !gtk_toggle_button_get_active(btn))
  {
    gtk_toggle_button_set_active(g->btn_point_tool, btn == g->btn_point_tool);
    gtk_toggle_button_set_active(g->btn_line_tool,  btn == g->btn_line_tool);
    gtk_toggle_button_set_active(g->btn_curve_tool, btn == g->btn_curve_tool);
    gtk_toggle_button_set_active(g->btn_node_tool,  btn == g->btn_node_tool);
 
    if(btn == g->btn_point_tool)
      dt_control_hinter_message
        (dt_control_get_global(), _("click and drag to add point\nscroll to change size - "
                              "shift+scroll to change strength - ctrl+scroll to change direction"));
    else if(btn == g->btn_line_tool)
      dt_control_hinter_message
        (dt_control_get_global(), _("click to add line\nscroll to change size - "
                              "shift+scroll to change strength - ctrl+scroll to change direction"));
    else if(btn == g->btn_curve_tool)
      dt_control_hinter_message
        (dt_control_get_global(), _("click to add curve\nscroll to change size - "
                              "shift+scroll to change strength - ctrl+scroll to change direction"));
    else if(btn == g->btn_node_tool)
      dt_control_hinter_message(dt_control_get_global(), _("click to edit nodes"));
 
    //  start the preview mode to show the shape that will be created
 
    if(btn == g->btn_point_tool || btn == g->btn_line_tool || btn == g->btn_curve_tool)
    {
      _start_new_shape(module);
    }
 
    if(!IS_NULL_PTR(btn)) dt_iop_request_focus(module);
  }
  else
  {
    gtk_toggle_button_set_active(btn, FALSE);
  }
 
  sync_pipe(module, FALSE);
 
  return TRUE;
}
 
void gui_update(dt_iop_module_t *module)
{
  dt_iop_liquify_gui_data_t *g = (dt_iop_liquify_gui_data_t *)dt_iop_gui_data(module);
  memcpy(&g->params, module->params, sizeof(dt_iop_liquify_params_t));
  update_warp_count(g);
}
 
void gui_init(dt_iop_module_t *self)
{
  dt_iop_liquify_gui_data_t *g = IOP_GUI_ALLOC(liquify);
 
  // A dummy surface for calculations only, no drawing.
  cairo_surface_t *cs = cairo_image_surface_create(CAIRO_FORMAT_ARGB32, 1, 1);
  cairo_surface_destroy(cs);
 
  g->dragging = NOWHERE;
  g->temp = NULL;
  g->status = 0;
  g->last_mouse_pos =
  g->last_button1_pressed_pos = -1;
  g->last_hit = NOWHERE;
  g->node_index = 0;
 
  self->gui->widget = gtk_box_new(GTK_ORIENTATION_VERTICAL, DT_GUI_BOX_SPACING);
 
  GtkWidget *hbox = gtk_box_new(GTK_ORIENTATION_HORIZONTAL, DT_GUI_BOX_SPACING);
  gtk_widget_set_tooltip_text(hbox, _("use a tool to add warps.\nright-click to remove a warp."));
  gtk_box_pack_start(GTK_BOX(self->gui->widget), hbox, TRUE, TRUE, 0);
 
 
  GtkWidget *lbox = gtk_box_new(GTK_ORIENTATION_HORIZONTAL, DT_GUI_BOX_SPACING);
  gtk_box_pack_start(GTK_BOX(hbox), lbox, FALSE, TRUE, 0);
 
  GtkWidget *labelbox = gtk_box_new(GTK_ORIENTATION_HORIZONTAL, DT_GUI_BOX_SPACING);
  GtkWidget *label = dt_ui_label_new(_("Warps: "));
  gtk_box_pack_start(GTK_BOX(labelbox), label, FALSE, TRUE, 0);
  g->label_warp = GTK_LABEL(dt_ui_label_new("-"));
  gtk_box_pack_start(GTK_BOX(labelbox), GTK_WIDGET(g->label_warp), FALSE, TRUE, 0);
  gtk_box_pack_start(GTK_BOX(lbox), labelbox, FALSE, TRUE, 0);
 
  GtkWidget *labelbox2 = gtk_box_new(GTK_ORIENTATION_HORIZONTAL, DT_GUI_BOX_SPACING);
  GtkWidget *label2 = dt_ui_label_new(_("Nodes: "));
  gtk_box_pack_start(GTK_BOX(labelbox2), label2, FALSE, TRUE, 0);
  g->label_node = GTK_LABEL(dt_ui_label_new("-"));
  gtk_box_pack_start(GTK_BOX(labelbox2), GTK_WIDGET(g->label_node), FALSE, TRUE, 0);
  gtk_box_pack_start(GTK_BOX(lbox), labelbox2, FALSE, TRUE, 0);
 
  hbox = gtk_box_new(GTK_ORIENTATION_HORIZONTAL, DT_GUI_BOX_SPACING);
  gtk_box_pack_start(GTK_BOX(self->gui->widget), hbox, TRUE, TRUE, 0);
 
  g->btn_node_tool = GTK_TOGGLE_BUTTON(dt_iop_togglebutton_new(self, NULL, N_("edit, add and delete nodes"), NULL,
                                       G_CALLBACK(btn_make_radio_callback), TRUE, 0, 0,
                                       dtgtk_cairo_paint_masks_edit, hbox));
 
  g->btn_curve_tool = GTK_TOGGLE_BUTTON(dt_iop_togglebutton_new(self, N_("shapes"), N_("draw curves"), N_("draw multiple curves"),
                                        G_CALLBACK(btn_make_radio_callback), TRUE, 0, 0,
                                        dtgtk_liquify_cairo_paint_curve_tool, hbox));
 
  g->btn_line_tool = GTK_TOGGLE_BUTTON(dt_iop_togglebutton_new(self, N_("shapes"), N_("draw lines"), N_("draw multiple lines"),
                                       G_CALLBACK(btn_make_radio_callback), TRUE, 0, 0,
                                       dtgtk_liquify_cairo_paint_line_tool, hbox));
 
  g->btn_point_tool = GTK_TOGGLE_BUTTON(dt_iop_togglebutton_new(self, N_("shapes"), N_("draw points"), N_("draw multiple points"),
                                         G_CALLBACK(btn_make_radio_callback), TRUE, 0, 0,
                                         dtgtk_liquify_cairo_paint_node_tool, hbox));
 
  dt_liquify_layers[DT_LIQUIFY_LAYER_PATH].hint           = _("drag: move segment - ctrl+click: add node - right click: remove path\n"
                                                              "ctrl+alt+click: toggle line/curve");
  dt_liquify_layers[DT_LIQUIFY_LAYER_CENTERPOINT].hint    = _("click and drag to move - click: show/hide feathering controls\n"
                                                              "ctrl+click: autosmooth, cusp, smooth, symmetrical"
                                                              " - right click to remove");
  dt_liquify_layers[DT_LIQUIFY_LAYER_CTRLPOINT1].hint     = _("drag to change shape of path");
  dt_liquify_layers[DT_LIQUIFY_LAYER_CTRLPOINT2].hint     = _("drag to change shape of path");
  dt_liquify_layers[DT_LIQUIFY_LAYER_RADIUSPOINT].hint    = _("drag to adjust warp radius");
  dt_liquify_layers[DT_LIQUIFY_LAYER_HARDNESSPOINT1].hint = _("drag to adjust hardness (center)");
  dt_liquify_layers[DT_LIQUIFY_LAYER_HARDNESSPOINT2].hint = _("drag to adjust hardness (feather)");
  dt_liquify_layers[DT_LIQUIFY_LAYER_STRENGTHPOINT].hint  = _("drag to adjust warp strength\n"
                                                              "ctrl+click: linear, grow, and shrink");
}
 
void gui_reset(dt_iop_module_t *self)
{
  dt_iop_liquify_gui_data_t *g = (dt_iop_liquify_gui_data_t *)dt_iop_gui_data(self);
  g->dragging = NOWHERE;
  g->temp = NULL;
  g->status = 0;
  btn_make_radio_callback(NULL, NULL, self);
}
 
void gui_cleanup(dt_iop_module_t *self)
{
  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