"Fossies" - the Fresh Open Source Software Archive

Member "darktable-3.6.1/src/iop/ashift.c" (10 Sep 2021, 169035 Bytes) of package /linux/misc/darktable-3.6.1.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "ashift.c" see the Fossies "Dox" file reference documentation and the last Fossies "Diffs" side-by-side code changes report: 3.4.1.1_vs_3.6.0.

    1 /*
    2   This file is part of darktable,
    3   Copyright (C) 2016-2021 darktable developers.
    4 
    5   darktable is free software: you can redistribute it and/or modify
    6   it under the terms of the GNU General Public License as published by
    7   the Free Software Foundation, either version 3 of the License, or
    8   (at your option) any later version.
    9 
   10   darktable is distributed in the hope that it will be useful,
   11   but WITHOUT ANY WARRANTY; without even the implied warranty of
   12   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   13   GNU General Public License for more details.
   14 
   15   You should have received a copy of the GNU General Public License
   16   along with darktable.  If not, see <http://www.gnu.org/licenses/>.
   17 */
   18 
   19 #ifdef HAVE_CONFIG_H
   20 #include "config.h"
   21 #endif
   22 #include "bauhaus/bauhaus.h"
   23 #include "common/bilateral.h"
   24 #include "common/colorspaces_inline_conversions.h"
   25 #include "common/debug.h"
   26 #include "common/imagebuf.h"
   27 #include "common/interpolation.h"
   28 #include "common/math.h"
   29 #include "common/opencl.h"
   30 #include "control/control.h"
   31 #include "develop/develop.h"
   32 #include "develop/imageop.h"
   33 #include "develop/imageop_gui.h"
   34 #include "develop/tiling.h"
   35 #include "dtgtk/button.h"
   36 #include "dtgtk/resetlabel.h"
   37 #include "gui/accelerators.h"
   38 #include "gui/draw.h"
   39 #include "gui/gtk.h"
   40 #include "gui/guides.h"
   41 #include "gui/presets.h"
   42 #include "iop/iop_api.h"
   43 #include "libs/modulegroups.h"
   44 
   45 #include <assert.h>
   46 #include <gtk/gtk.h>
   47 #include <inttypes.h>
   48 #include <math.h>
   49 #include <stdlib.h>
   50 #include <string.h>
   51 
   52 // Inspiration to this module comes from the program ShiftN (http://www.shiftn.de) by
   53 // Marcus Hebel.
   54 
   55 // Thanks to Marcus for his support when implementing part of the ShiftN functionality
   56 // to darktable.
   57 
   58 #define ROTATION_RANGE 10                   // allowed min/max default range for rotation parameter
   59 #define ROTATION_RANGE_SOFT 20              // allowed min/max range for rotation parameter with manual adjustment
   60 #define LENSSHIFT_RANGE 1.0                 // allowed min/max default range for lensshift parameters
   61 #define LENSSHIFT_RANGE_SOFT 2.0            // allowed min/max range for lensshift parameters with manual adjustment
   62 #define SHEAR_RANGE 0.2                     // allowed min/max range for shear parameter
   63 #define SHEAR_RANGE_SOFT 0.5                // allowed min/max range for shear parameter with manual adjustment
   64 #define MIN_LINE_LENGTH 5                   // the minimum length of a line in pixels to be regarded as relevant
   65 #define MAX_TANGENTIAL_DEVIATION 30         // by how many degrees a line may deviate from the +/-180 and +/-90 to be regarded as relevant
   66 #define LSD_SCALE 0.99                      // LSD: scaling factor for line detection
   67 #define LSD_SIGMA_SCALE 0.6                 // LSD: sigma for Gaussian filter is computed as sigma = sigma_scale/scale
   68 #define LSD_QUANT 2.0                       // LSD: bound to the quantization error on the gradient norm
   69 #define LSD_ANG_TH 22.5                     // LSD: gradient angle tolerance in degrees
   70 #define LSD_LOG_EPS 0.0                     // LSD: detection threshold: -log10(NFA) > log_eps
   71 #define LSD_DENSITY_TH 0.7                  // LSD: minimal density of region points in rectangle
   72 #define LSD_N_BINS 1024                     // LSD: number of bins in pseudo-ordering of gradient modulus
   73 #define LSD_GAMMA 0.45                      // gamma correction to apply on raw images prior to line detection
   74 #define RANSAC_RUNS 400                     // how many iterations to run in ransac
   75 #define RANSAC_EPSILON 2                    // starting value for ransac epsilon (in -log10 units)
   76 #define RANSAC_EPSILON_STEP 1               // step size of epsilon optimization (log10 units)
   77 #define RANSAC_ELIMINATION_RATIO 60         // percentage of lines we try to eliminate as outliers
   78 #define RANSAC_OPTIMIZATION_STEPS 5         // home many steps to optimize epsilon
   79 #define RANSAC_OPTIMIZATION_DRY_RUNS 50     // how man runs per optimization steps
   80 #define RANSAC_HURDLE 5                     // hurdle rate: the number of lines below which we do a complete permutation instead of random sampling
   81 #define MINIMUM_FITLINES 4                  // minimum number of lines needed for automatic parameter fit
   82 #define NMS_EPSILON 1e-3                    // break criterion for Nelder-Mead simplex
   83 #define NMS_SCALE 1.0                       // scaling factor for Nelder-Mead simplex
   84 #define NMS_ITERATIONS 400                  // number of iterations for Nelder-Mead simplex
   85 #define NMS_CROP_EPSILON 100.0              // break criterion for Nelder-Mead simplex on crop fitting
   86 #define NMS_CROP_SCALE 0.5                  // scaling factor for Nelder-Mead simplex on crop fitting
   87 #define NMS_CROP_ITERATIONS 100             // number of iterations for Nelder-Mead simplex on crop fitting
   88 #define NMS_ALPHA 1.0                       // reflection coefficient for Nelder-Mead simplex
   89 #define NMS_BETA 0.5                        // contraction coefficient for Nelder-Mead simplex
   90 #define NMS_GAMMA 2.0                       // expansion coefficient for Nelder-Mead simplex
   91 #define DEFAULT_F_LENGTH 28.0               // focal length we assume if no exif data are available
   92 
   93 // define to get debugging output
   94 #undef ASHIFT_DEBUG
   95 
   96 #define SQR(a) ((a) * (a))
   97 
   98 // For line detection we use the LSD algorithm as published by Rafael Grompone:
   99 //
  100 //  "LSD: a Line Segment Detector" by Rafael Grompone von Gioi,
  101 //  Jeremie Jakubowicz, Jean-Michel Morel, and Gregory Randall,
  102 //  Image Processing On Line, 2012. DOI:10.5201/ipol.2012.gjmr-lsd
  103 //  http://dx.doi.org/10.5201/ipol.2012.gjmr-lsd
  104 #include "ashift_lsd.c"
  105 
  106 // For parameter optimization we are using the Nelder-Mead simplex method
  107 // implemented by Michael F. Hutt.
  108 #include "ashift_nmsimplex.c"
  109 
  110 
  111 DT_MODULE_INTROSPECTION(4, dt_iop_ashift_params_t)
  112 
  113 
  114 const char *name()
  115 {
  116   return _("perspective correction");
  117 }
  118 
  119 const char *aliases()
  120 {
  121   return _("keystone|distortion");
  122 }
  123 
  124 const char *description(struct dt_iop_module_t *self)
  125 {
  126   return dt_iop_set_description(self, _("distort perspective automatically"),
  127                                       _("corrective"),
  128                                       _("linear, RGB, scene-referred"),
  129                                       _("geometric, RGB"),
  130                                       _("linear, RGB, scene-referred"));
  131 }
  132 
  133 int flags()
  134 {
  135   return IOP_FLAGS_ALLOW_TILING | IOP_FLAGS_TILING_FULL_ROI | IOP_FLAGS_ONE_INSTANCE | IOP_FLAGS_ALLOW_FAST_PIPE;
  136 }
  137 
  138 int default_group()
  139 {
  140   return IOP_GROUP_CORRECT | IOP_GROUP_TECHNICAL;
  141 }
  142 
  143 int operation_tags()
  144 {
  145   return IOP_TAG_DISTORT;
  146 }
  147 
  148 int operation_tags_filter()
  149 {
  150   // switch off clipping and decoration, we want to see the full image.
  151   return IOP_TAG_DECORATION | IOP_TAG_CLIPPING;
  152 }
  153 
  154 int default_colorspace(dt_iop_module_t *self, dt_dev_pixelpipe_t *pipe, dt_dev_pixelpipe_iop_t *piece)
  155 {
  156   return iop_cs_rgb;
  157 }
  158 
  159 typedef enum dt_iop_ashift_homodir_t
  160 {
  161   ASHIFT_HOMOGRAPH_FORWARD,
  162   ASHIFT_HOMOGRAPH_INVERTED
  163 } dt_iop_ashift_homodir_t;
  164 
  165 typedef enum dt_iop_ashift_linetype_t
  166 {
  167   ASHIFT_LINE_IRRELEVANT   = 0,       // the line is found to be not interesting
  168                                       // eg. too short, or not horizontal or vertical
  169   ASHIFT_LINE_RELEVANT     = 1 << 0,  // the line is relevant for us
  170   ASHIFT_LINE_DIRVERT      = 1 << 1,  // the line is (mostly) vertical, else (mostly) horizontal
  171   ASHIFT_LINE_SELECTED     = 1 << 2,  // the line is selected for fitting
  172   ASHIFT_LINE_VERTICAL_NOT_SELECTED   = ASHIFT_LINE_RELEVANT | ASHIFT_LINE_DIRVERT,
  173   ASHIFT_LINE_HORIZONTAL_NOT_SELECTED = ASHIFT_LINE_RELEVANT,
  174   ASHIFT_LINE_VERTICAL_SELECTED = ASHIFT_LINE_RELEVANT | ASHIFT_LINE_DIRVERT | ASHIFT_LINE_SELECTED,
  175   ASHIFT_LINE_HORIZONTAL_SELECTED = ASHIFT_LINE_RELEVANT | ASHIFT_LINE_SELECTED,
  176   ASHIFT_LINE_MASK = ASHIFT_LINE_RELEVANT | ASHIFT_LINE_DIRVERT | ASHIFT_LINE_SELECTED
  177 } dt_iop_ashift_linetype_t;
  178 
  179 typedef enum dt_iop_ashift_linecolor_t
  180 {
  181   ASHIFT_LINECOLOR_GREY    = 0,
  182   ASHIFT_LINECOLOR_GREEN   = 1,
  183   ASHIFT_LINECOLOR_RED     = 2,
  184   ASHIFT_LINECOLOR_BLUE    = 3,
  185   ASHIFT_LINECOLOR_YELLOW  = 4
  186 } dt_iop_ashift_linecolor_t;
  187 
  188 typedef enum dt_iop_ashift_fitaxis_t
  189 {
  190   ASHIFT_FIT_NONE          = 0,       // none
  191   ASHIFT_FIT_ROTATION      = 1 << 0,  // flag indicates to fit rotation angle
  192   ASHIFT_FIT_LENS_VERT     = 1 << 1,  // flag indicates to fit vertical lens shift
  193   ASHIFT_FIT_LENS_HOR      = 1 << 2,  // flag indicates to fit horizontal lens shift
  194   ASHIFT_FIT_SHEAR         = 1 << 3,  // flag indicates to fit shear parameter
  195   ASHIFT_FIT_LINES_VERT    = 1 << 4,  // use vertical lines for fitting
  196   ASHIFT_FIT_LINES_HOR     = 1 << 5,  // use horizontal lines for fitting
  197   ASHIFT_FIT_LENS_BOTH = ASHIFT_FIT_LENS_VERT | ASHIFT_FIT_LENS_HOR,
  198   ASHIFT_FIT_LINES_BOTH = ASHIFT_FIT_LINES_VERT | ASHIFT_FIT_LINES_HOR,
  199   ASHIFT_FIT_VERTICALLY = ASHIFT_FIT_ROTATION | ASHIFT_FIT_LENS_VERT | ASHIFT_FIT_LINES_VERT,
  200   ASHIFT_FIT_HORIZONTALLY = ASHIFT_FIT_ROTATION | ASHIFT_FIT_LENS_HOR | ASHIFT_FIT_LINES_HOR,
  201   ASHIFT_FIT_BOTH = ASHIFT_FIT_ROTATION | ASHIFT_FIT_LENS_VERT | ASHIFT_FIT_LENS_HOR |
  202                     ASHIFT_FIT_LINES_VERT | ASHIFT_FIT_LINES_HOR,
  203   ASHIFT_FIT_VERTICALLY_NO_ROTATION = ASHIFT_FIT_LENS_VERT | ASHIFT_FIT_LINES_VERT,
  204   ASHIFT_FIT_HORIZONTALLY_NO_ROTATION = ASHIFT_FIT_LENS_HOR | ASHIFT_FIT_LINES_HOR,
  205   ASHIFT_FIT_BOTH_NO_ROTATION = ASHIFT_FIT_LENS_VERT | ASHIFT_FIT_LENS_HOR |
  206                                 ASHIFT_FIT_LINES_VERT | ASHIFT_FIT_LINES_HOR,
  207   ASHIFT_FIT_BOTH_SHEAR = ASHIFT_FIT_ROTATION | ASHIFT_FIT_LENS_VERT | ASHIFT_FIT_LENS_HOR |
  208                     ASHIFT_FIT_SHEAR | ASHIFT_FIT_LINES_VERT | ASHIFT_FIT_LINES_HOR,
  209   ASHIFT_FIT_ROTATION_VERTICAL_LINES = ASHIFT_FIT_ROTATION | ASHIFT_FIT_LINES_VERT,
  210   ASHIFT_FIT_ROTATION_HORIZONTAL_LINES = ASHIFT_FIT_ROTATION | ASHIFT_FIT_LINES_HOR,
  211   ASHIFT_FIT_ROTATION_BOTH_LINES = ASHIFT_FIT_ROTATION | ASHIFT_FIT_LINES_VERT | ASHIFT_FIT_LINES_HOR,
  212   ASHIFT_FIT_FLIP = ASHIFT_FIT_LENS_VERT | ASHIFT_FIT_LENS_HOR | ASHIFT_FIT_LINES_VERT | ASHIFT_FIT_LINES_HOR
  213 } dt_iop_ashift_fitaxis_t;
  214 
  215 typedef enum dt_iop_ashift_nmsresult_t
  216 {
  217   NMS_SUCCESS = 0,
  218   NMS_NOT_ENOUGH_LINES = 1,
  219   NMS_DID_NOT_CONVERGE = 2,
  220   NMS_INSANE = 3
  221 } dt_iop_ashift_nmsresult_t;
  222 
  223 typedef enum dt_iop_ashift_enhance_t
  224 {
  225   ASHIFT_ENHANCE_NONE       = 0,
  226   ASHIFT_ENHANCE_EDGES      = 1 << 0,
  227   ASHIFT_ENHANCE_DETAIL     = 1 << 1,
  228   ASHIFT_ENHANCE_HORIZONTAL = 0x100,
  229   ASHIFT_ENHANCE_VERTICAL   = 0x200
  230 } dt_iop_ashift_enhance_t;
  231 
  232 typedef enum dt_iop_ashift_mode_t
  233 {
  234   ASHIFT_MODE_GENERIC = 0, // $DESCRIPTION: "generic"
  235   ASHIFT_MODE_SPECIFIC = 1 // $DESCRIPTION: "specific"
  236 } dt_iop_ashift_mode_t;
  237 
  238 typedef enum dt_iop_ashift_crop_t
  239 {
  240   ASHIFT_CROP_OFF = 0,    // $DESCRIPTION: "off"
  241   ASHIFT_CROP_LARGEST = 1,// $DESCRIPTION: "largest area"
  242   ASHIFT_CROP_ASPECT = 2  // $DESCRIPTION: "original format"
  243 } dt_iop_ashift_crop_t;
  244 
  245 typedef enum dt_iop_ashift_bounding_t
  246 {
  247   ASHIFT_BOUNDING_OFF = 0,
  248   ASHIFT_BOUNDING_SELECT = 1,
  249   ASHIFT_BOUNDING_DESELECT = 2
  250 } dt_iop_ashift_bounding_t;
  251 
  252 typedef enum dt_iop_ashift_jobcode_t
  253 {
  254   ASHIFT_JOBCODE_NONE = 0,
  255   ASHIFT_JOBCODE_GET_STRUCTURE = 1,
  256   ASHIFT_JOBCODE_FIT = 2
  257 } dt_iop_ashift_jobcode_t;
  258 
  259 typedef struct dt_iop_ashift_params1_t
  260 {
  261   float rotation;
  262   float lensshift_v;
  263   float lensshift_h;
  264   int toggle;
  265 } dt_iop_ashift_params1_t;
  266 
  267 typedef struct dt_iop_ashift_params2_t
  268 {
  269   float rotation;
  270   float lensshift_v;
  271   float lensshift_h;
  272   float f_length;
  273   float crop_factor;
  274   float orthocorr;
  275   float aspect;
  276   dt_iop_ashift_mode_t mode;
  277   int toggle;
  278 } dt_iop_ashift_params2_t;
  279 
  280 typedef struct dt_iop_ashift_params3_t
  281 {
  282   float rotation;
  283   float lensshift_v;
  284   float lensshift_h;
  285   float f_length;
  286   float crop_factor;
  287   float orthocorr;
  288   float aspect;
  289   dt_iop_ashift_mode_t mode;
  290   int toggle;
  291   dt_iop_ashift_crop_t cropmode;
  292   float cl;
  293   float cr;
  294   float ct;
  295   float cb;
  296 } dt_iop_ashift_params3_t;
  297 
  298 typedef struct dt_iop_ashift_params_t
  299 {
  300   float rotation;    // $MIN: -ROTATION_RANGE_SOFT $MAX: ROTATION_RANGE_SOFT $DEFAULT: 0.0
  301   float lensshift_v; // $MIN: -LENSSHIFT_RANGE_SOFT $MAX: LENSSHIFT_RANGE_SOFT $DEFAULT: 0.0 $DESCRIPTION: "lens shift (vertical)"
  302   float lensshift_h; // $MIN: -LENSSHIFT_RANGE_SOFT $MAX: LENSSHIFT_RANGE_SOFT $DEFAULT: 0.0 $DESCRIPTION: "lens shift (horizontal)"
  303   float shear;       // $MIN: -SHEAR_RANGE_SOFT $MAX: SHEAR_RANGE_SOFT $DEFAULT: 0.0 $DESCRIPTION: "shear"
  304   float f_length;    // $MIN: 1.0 $MAX: 2000.0 $DEFAULT: DEFAULT_F_LENGTH $DESCRIPTION: "focal length"
  305   float crop_factor; // $MIN: 0.5 $MAX: 10.0 $DEFAULT: 1.0 $DESCRIPTION: "crop factor"
  306   float orthocorr;   // $MIN: 0.0 $MAX: 100.0 $DEFAULT: 100.0 $DESCRIPTION: "lens dependence"
  307   float aspect;      // $MIN: 0.5 $MAX: 2.0 $DEFAULT: 1.0 $DESCRIPTION: "aspect adjust"
  308   dt_iop_ashift_mode_t mode;     // $DEFAULT: ASHIFT_MODE_GENERIC $DESCRIPTION: "lens model"
  309   int toggle;                    // $DEFAULT: 0
  310   dt_iop_ashift_crop_t cropmode; // $DEFAULT: ASHIFT_CROP_OFF $DESCRIPTION: "automatic cropping"
  311   float cl;          // $DEFAULT: 0.0
  312   float cr;          // $DEFAULT: 1.0
  313   float ct;          // $DEFAULT: 0.0
  314   float cb;          // $DEFAULT: 1.0
  315 } dt_iop_ashift_params_t;
  316 
  317 typedef struct dt_iop_ashift_line_t
  318 {
  319   float p1[3];
  320   float p2[3];
  321   float length;
  322   float width;
  323   float weight;
  324   dt_iop_ashift_linetype_t type;
  325   // homogeneous coordinates:
  326   float L[3];
  327 } dt_iop_ashift_line_t;
  328 
  329 typedef struct dt_iop_ashift_points_idx_t
  330 {
  331   size_t offset;
  332   int length;
  333   int near;
  334   int bounded;
  335   dt_iop_ashift_linetype_t type;
  336   dt_iop_ashift_linecolor_t color;
  337   // bounding box:
  338   float bbx, bby, bbX, bbY;
  339 } dt_iop_ashift_points_idx_t;
  340 
  341 typedef struct dt_iop_ashift_fit_params_t
  342 {
  343   int params_count;
  344   dt_iop_ashift_linetype_t linetype;
  345   dt_iop_ashift_linetype_t linemask;
  346   dt_iop_ashift_line_t *lines;
  347   int lines_count;
  348   int width;
  349   int height;
  350   float weight;
  351   float f_length_kb;
  352   float orthocorr;
  353   float aspect;
  354   float rotation;
  355   float lensshift_v;
  356   float lensshift_h;
  357   float shear;
  358   float rotation_range;
  359   float lensshift_v_range;
  360   float lensshift_h_range;
  361   float shear_range;
  362 } dt_iop_ashift_fit_params_t;
  363 
  364 typedef struct dt_iop_ashift_cropfit_params_t
  365 {
  366   int width;
  367   int height;
  368   float x;
  369   float y;
  370   float alpha;
  371   float homograph[3][3];
  372   float edges[4][3];
  373 } dt_iop_ashift_cropfit_params_t;
  374 
  375 typedef struct dt_iop_ashift_gui_data_t
  376 {
  377   GtkWidget *rotation;
  378   GtkWidget *lensshift_v;
  379   GtkWidget *lensshift_h;
  380   GtkWidget *shear;
  381   GtkWidget *guide_lines;
  382   GtkWidget *cropmode;
  383   GtkWidget *mode;
  384   GtkWidget *specifics;
  385   GtkWidget *f_length;
  386   GtkWidget *crop_factor;
  387   GtkWidget *orthocorr;
  388   GtkWidget *aspect;
  389   GtkWidget *fit_v;
  390   GtkWidget *fit_h;
  391   GtkWidget *fit_both;
  392   GtkWidget *structure;
  393   GtkWidget *clean;
  394   GtkWidget *eye;
  395   int lines_suppressed;
  396   int fitting;
  397   int isflipped;
  398   int show_guides;
  399   int isselecting;
  400   int isdeselecting;
  401   dt_iop_ashift_bounding_t isbounding;
  402   float near_delta;
  403   int selecting_lines_version;
  404   float rotation_range;
  405   float lensshift_v_range;
  406   float lensshift_h_range;
  407   float shear_range;
  408   dt_iop_ashift_line_t *lines;
  409   int lines_in_width;
  410   int lines_in_height;
  411   int lines_x_off;
  412   int lines_y_off;
  413   int lines_count;
  414   int vertical_count;
  415   int horizontal_count;
  416   int lines_version;
  417   float vertical_weight;
  418   float horizontal_weight;
  419   float *points;
  420   dt_iop_ashift_points_idx_t *points_idx;
  421   int points_lines_count;
  422   int points_version;
  423   float *buf;
  424   int buf_width;
  425   int buf_height;
  426   int buf_x_off;
  427   int buf_y_off;
  428   float buf_scale;
  429   uint64_t lines_hash;
  430   uint64_t grid_hash;
  431   uint64_t buf_hash;
  432   dt_iop_ashift_fitaxis_t lastfit;
  433   float lastx;
  434   float lasty;
  435   float crop_cx;
  436   float crop_cy;
  437   dt_iop_ashift_jobcode_t jobcode;
  438   int jobparams;
  439   gboolean adjust_crop;
  440   float cl; // shadow copy of dt_iop_ashift_data_t.cl
  441   float cr; // shadow copy of dt_iop_ashift_data_t.cr
  442   float ct; // shadow copy of dt_iop_ashift_data_t.ct
  443   float cb; // shadow copy of dt_iop_ashift_data_t.cb
  444 } dt_iop_ashift_gui_data_t;
  445 
  446 typedef struct dt_iop_ashift_data_t
  447 {
  448   float rotation;
  449   float lensshift_v;
  450   float lensshift_h;
  451   float shear;
  452   float f_length_kb;
  453   float orthocorr;
  454   float aspect;
  455   float cl;
  456   float cr;
  457   float ct;
  458   float cb;
  459 } dt_iop_ashift_data_t;
  460 
  461 typedef struct dt_iop_ashift_global_data_t
  462 {
  463   int kernel_ashift_bilinear;
  464   int kernel_ashift_bicubic;
  465   int kernel_ashift_lanczos2;
  466   int kernel_ashift_lanczos3;
  467 } dt_iop_ashift_global_data_t;
  468 
  469 int legacy_params(dt_iop_module_t *self, const void *const old_params, const int old_version,
  470                   void *new_params, const int new_version)
  471 {
  472   if(old_version == 1 && new_version == 4)
  473   {
  474     const dt_iop_ashift_params1_t *old = old_params;
  475     dt_iop_ashift_params_t *new = new_params;
  476     new->rotation = old->rotation;
  477     new->lensshift_v = old->lensshift_v;
  478     new->lensshift_h = old->lensshift_h;
  479     new->shear = 0.0f;
  480     new->toggle = old->toggle;
  481     new->f_length = DEFAULT_F_LENGTH;
  482     new->crop_factor = 1.0f;
  483     new->orthocorr = 100.0f;
  484     new->aspect = 1.0f;
  485     new->mode = ASHIFT_MODE_GENERIC;
  486     new->cropmode = ASHIFT_CROP_OFF;
  487     new->cl = 0.0f;
  488     new->cr = 1.0f;
  489     new->ct = 0.0f;
  490     new->cb = 1.0f;
  491     return 0;
  492   }
  493   if(old_version == 2 && new_version == 4)
  494   {
  495     const dt_iop_ashift_params2_t *old = old_params;
  496     dt_iop_ashift_params_t *new = new_params;
  497     new->rotation = old->rotation;
  498     new->lensshift_v = old->lensshift_v;
  499     new->lensshift_h = old->lensshift_h;
  500     new->shear = 0.0f;
  501     new->toggle = old->toggle;
  502     new->f_length = old->f_length;
  503     new->crop_factor = old->crop_factor;
  504     new->orthocorr = old->orthocorr;
  505     new->aspect = old->aspect;
  506     new->mode = old->mode;
  507     new->cropmode = ASHIFT_CROP_OFF;
  508     new->cl = 0.0f;
  509     new->cr = 1.0f;
  510     new->ct = 0.0f;
  511     new->cb = 1.0f;
  512     return 0;
  513   }
  514   if(old_version == 3 && new_version == 4)
  515   {
  516     const dt_iop_ashift_params3_t *old = old_params;
  517     dt_iop_ashift_params_t *new = new_params;
  518     new->rotation = old->rotation;
  519     new->lensshift_v = old->lensshift_v;
  520     new->lensshift_h = old->lensshift_h;
  521     new->shear = 0.0f;
  522     new->toggle = old->toggle;
  523     new->f_length = old->f_length;
  524     new->crop_factor = old->crop_factor;
  525     new->orthocorr = old->orthocorr;
  526     new->aspect = old->aspect;
  527     new->mode = old->mode;
  528     new->cropmode = old->cropmode;
  529     new->cl = old->cl;
  530     new->cr = old->cr;
  531     new->ct = old->ct;
  532     new->cb = old->cb;
  533     return 0;
  534   }
  535 
  536   return 1;
  537 }
  538 
  539 // normalized product of two 3x1 vectors
  540 // dst needs to be different from v1 and v2
  541 static inline void vec3prodn(float *dst, const float *const v1, const float *const v2)
  542 {
  543   const float l1 = v1[1] * v2[2] - v1[2] * v2[1];
  544   const float l2 = v1[2] * v2[0] - v1[0] * v2[2];
  545   const float l3 = v1[0] * v2[1] - v1[1] * v2[0];
  546 
  547   // normalize so that l1^2 + l2^2 + l3^3 = 1
  548   const float sq = sqrtf(l1 * l1 + l2 * l2 + l3 * l3);
  549 
  550   const float f = sq > 0.0f ? 1.0f / sq : 1.0f;
  551 
  552   dst[0] = l1 * f;
  553   dst[1] = l2 * f;
  554   dst[2] = l3 * f;
  555 }
  556 
  557 // normalize a 3x1 vector so that x^2 + y^2 + z^2 = 1
  558 // dst and v may be the same
  559 static inline void vec3norm(float *dst, const float *const v)
  560 {
  561   const float sq = sqrtf(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
  562 
  563   // special handling for an all-zero vector
  564   const float f = sq > 0.0f ? 1.0f / sq : 1.0f;
  565 
  566   dst[0] = v[0] * f;
  567   dst[1] = v[1] * f;
  568   dst[2] = v[2] * f;
  569 }
  570 
  571 // normalize a 3x1 vector so that x^2 + y^2 = 1; a useful normalization for
  572 // lines in homogeneous coordinates
  573 // dst and v may be the same
  574 static inline void vec3lnorm(float *dst, const float *const v)
  575 {
  576   const float sq = sqrtf(v[0] * v[0] + v[1] * v[1]);
  577 
  578   // special handling for a point vector of the image center
  579   const float f = sq > 0.0f ? 1.0f / sq : 1.0f;
  580 
  581   dst[0] = v[0] * f;
  582   dst[1] = v[1] * f;
  583   dst[2] = v[2] * f;
  584 }
  585 
  586 
  587 // scalar product of two 3x1 vectors
  588 static inline float vec3scalar(const float *const v1, const float *const v2)
  589 {
  590   return (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]);
  591 }
  592 
  593 // check if 3x1 vector is (very close to) null
  594 static inline int vec3isnull(const float *const v)
  595 {
  596   const float eps = 1e-10f;
  597   return (fabsf(v[0]) < eps && fabsf(v[1]) < eps && fabsf(v[2]) < eps);
  598 }
  599 
  600 #ifdef ASHIFT_DEBUG
  601 static void print_roi(const dt_iop_roi_t *roi, const char *label)
  602 {
  603   printf("{ %5d  %5d  %5d  %5d  %.6f } %s\n", roi->x, roi->y, roi->width, roi->height, roi->scale, label);
  604 }
  605 #endif
  606 
  607 static inline void shadow_crop_box(dt_iop_ashift_params_t *p, dt_iop_ashift_gui_data_t *g)
  608 {
  609   // copy actual crop box values into shadow variables
  610   g->cl = p->cl;
  611   g->cr = p->cr;
  612   g->ct = p->ct;
  613   g->cb = p->cb;
  614 }
  615 
  616 static void clear_shadow_crop_box(dt_iop_ashift_gui_data_t *g)
  617 {
  618   // reset the crop to the full image
  619   g->cl = 0.0f;
  620   g->cr = 1.0f;
  621   g->ct = 0.0f;
  622   g->cb = 1.0f;
  623 }
  624 
  625 static inline void commit_crop_box(dt_iop_ashift_params_t *p, dt_iop_ashift_gui_data_t *g)
  626 {
  627   // copy shadow values for crop box into actual parameters
  628   p->cl = g->cl;
  629   p->cr = g->cr;
  630   p->ct = g->ct;
  631   p->cb = g->cb;
  632 }
  633 
  634 static inline void swap_shadow_crop_box(dt_iop_ashift_params_t *p, dt_iop_ashift_gui_data_t *g)
  635 {
  636   // exchange shadow values and actual crop values
  637   // this is needed for a temporary commit to be able to properly update the undo history
  638   float tmp;
  639   tmp = p->cl; p->cl = g->cl; g->cl = tmp;
  640   tmp = p->cr; p->cr = g->cr; g->cr = tmp;
  641   tmp = p->ct; p->ct = g->ct; g->ct = tmp;
  642   tmp = p->cb; p->cb = g->cb; g->cb = tmp;
  643 }
  644 
  645 #define MAT3SWAP(a, b) { float (*tmp)[3] = (a); (a) = (b); (b) = tmp; }
  646 
  647 static void homography(float *homograph, const float angle, const float shift_v, const float shift_h,
  648                        const float shear, const float f_length_kb, const float orthocorr, const float aspect,
  649                        const int width, const int height, dt_iop_ashift_homodir_t dir)
  650 {
  651   // calculate homograph that combines all translations, rotations
  652   // and warping into one single matrix operation.
  653   // this is heavily leaning on ShiftN where the homographic matrix expects
  654   // input in (y : x : 1) format. in the darktable world we want to keep the
  655   // (x : y : 1) convention. therefore we need to flip coordinates first and
  656   // make sure that output is in correct format after corrections are applied.
  657 
  658   const float u = width;
  659   const float v = height;
  660 
  661   const float phi = M_PI * angle / 180.0f;
  662   const float cosi = cosf(phi);
  663   const float sini = sinf(phi);
  664   const float ascale = sqrtf(aspect);
  665 
  666   // most of this comes from ShiftN
  667   const float f_global = f_length_kb;
  668   const float horifac = 1.0f - orthocorr / 100.0f;
  669   const float exppa_v = expf(shift_v);
  670   const float fdb_v = f_global / (14.4f + (v / u - 1) * 7.2f);
  671   const float rad_v = fdb_v * (exppa_v - 1.0f) / (exppa_v + 1.0f);
  672   const float alpha_v = CLAMP(atanf(rad_v), -1.5f, 1.5f);
  673   const float rt_v = sinf(0.5f * alpha_v);
  674   const float r_v = fmaxf(0.1f, 2.0f * (horifac - 1.0f) * rt_v * rt_v + 1.0f);
  675 
  676   const float vertifac = 1.0f - orthocorr / 100.0f;
  677   const float exppa_h = expf(shift_h);
  678   const float fdb_h = f_global / (14.4f + (u / v - 1) * 7.2f);
  679   const float rad_h = fdb_h * (exppa_h - 1.0f) / (exppa_h + 1.0f);
  680   const float alpha_h = CLAMP(atanf(rad_h), -1.5f, 1.5f);
  681   const float rt_h = sinf(0.5f * alpha_h);
  682   const float r_h = fmaxf(0.1f, 2.0f * (vertifac - 1.0f) * rt_h * rt_h + 1.0f);
  683 
  684 
  685   // three intermediate buffers for matrix calculation ...
  686   float m1[3][3], m2[3][3], m3[3][3];
  687 
  688   // ... and some pointers to handle them more intuitively
  689   float (*mwork)[3] = m1;
  690   float (*minput)[3] = m2;
  691   float (*moutput)[3] = m3;
  692 
  693   // Step 1: flip x and y coordinates (see above)
  694   memset(minput, 0, sizeof(float) * 9);
  695   minput[0][1] = 1.0f;
  696   minput[1][0] = 1.0f;
  697   minput[2][2] = 1.0f;
  698 
  699 
  700   // Step 2: rotation of image around its center
  701   memset(mwork, 0, sizeof(float) * 9);
  702   mwork[0][0] = cosi;
  703   mwork[0][1] = -sini;
  704   mwork[1][0] = sini;
  705   mwork[1][1] = cosi;
  706   mwork[0][2] = -0.5f * v * cosi + 0.5f * u * sini + 0.5f * v;
  707   mwork[1][2] = -0.5f * v * sini - 0.5f * u * cosi + 0.5f * u;
  708   mwork[2][2] = 1.0f;
  709 
  710   // multiply mwork * minput -> moutput
  711   mat3mul((float *)moutput, (float *)mwork, (float *)minput);
  712 
  713 
  714   // Step 3: apply shearing
  715   memset(mwork, 0, sizeof(float) * 9);
  716   mwork[0][0] = 1.0f;
  717   mwork[0][1] = shear;
  718   mwork[1][1] = 1.0f;
  719   mwork[1][0] = shear;
  720   mwork[2][2] = 1.0f;
  721 
  722   // moutput (of last calculation) -> minput
  723   MAT3SWAP(minput, moutput);
  724   // multiply mwork * minput -> moutput
  725   mat3mul((float *)moutput, (float *)mwork, (float *)minput);
  726 
  727 
  728   // Step 4: apply vertical lens shift effect
  729   memset(mwork, 0, sizeof(float) * 9);
  730   mwork[0][0] = exppa_v;
  731   mwork[1][0] = 0.5f * ((exppa_v - 1.0f) * u) / v;
  732   mwork[1][1] = 2.0f * exppa_v / (exppa_v + 1.0f);
  733   mwork[1][2] = -0.5f * ((exppa_v - 1.0f) * u) / (exppa_v + 1.0f);
  734   mwork[2][0] = (exppa_v - 1.0f) / v;
  735   mwork[2][2] = 1.0f;
  736 
  737   // moutput (of last calculation) -> minput
  738   MAT3SWAP(minput, moutput);
  739   // multiply mwork * minput -> moutput
  740   mat3mul((float *)moutput, (float *)mwork, (float *)minput);
  741 
  742 
  743   // Step 5: horizontal compression
  744   memset(mwork, 0, sizeof(float) * 9);
  745   mwork[0][0] = 1.0f;
  746   mwork[1][1] = r_v;
  747   mwork[1][2] = 0.5f * u * (1.0f - r_v);
  748   mwork[2][2] = 1.0f;
  749 
  750   // moutput (of last calculation) -> minput
  751   MAT3SWAP(minput, moutput);
  752   // multiply mwork * minput -> moutput
  753   mat3mul((float *)moutput, (float *)mwork, (float *)minput);
  754 
  755 
  756   // Step 6: flip x and y back again
  757   memset(mwork, 0, sizeof(float) * 9);
  758   mwork[0][1] = 1.0f;
  759   mwork[1][0] = 1.0f;
  760   mwork[2][2] = 1.0f;
  761 
  762   // moutput (of last calculation) -> minput
  763   MAT3SWAP(minput, moutput);
  764   // multiply mwork * minput -> moutput
  765   mat3mul((float *)moutput, (float *)mwork, (float *)minput);
  766 
  767 
  768   // from here output vectors would be in (x : y : 1) format
  769 
  770   // Step 7: now we can apply horizontal lens shift with the same matrix format as above
  771   memset(mwork, 0, sizeof(float) * 9);
  772   mwork[0][0] = exppa_h;
  773   mwork[1][0] = 0.5f * ((exppa_h - 1.0f) * v) / u;
  774   mwork[1][1] = 2.0f * exppa_h / (exppa_h + 1.0f);
  775   mwork[1][2] = -0.5f * ((exppa_h - 1.0f) * v) / (exppa_h + 1.0f);
  776   mwork[2][0] = (exppa_h - 1.0f) / u;
  777   mwork[2][2] = 1.0f;
  778 
  779   // moutput (of last calculation) -> minput
  780   MAT3SWAP(minput, moutput);
  781   // multiply mwork * minput -> moutput
  782   mat3mul((float *)moutput, (float *)mwork, (float *)minput);
  783 
  784 
  785   // Step 8: vertical compression
  786   memset(mwork, 0, sizeof(float) * 9);
  787   mwork[0][0] = 1.0f;
  788   mwork[1][1] = r_h;
  789   mwork[1][2] = 0.5f * v * (1.0f - r_h);
  790   mwork[2][2] = 1.0f;
  791 
  792   // moutput (of last calculation) -> minput
  793   MAT3SWAP(minput, moutput);
  794   // multiply mwork * minput -> moutput
  795   mat3mul((float *)moutput, (float *)mwork, (float *)minput);
  796 
  797 
  798   // Step 9: apply aspect ratio scaling
  799   memset(mwork, 0, sizeof(float) * 9);
  800   mwork[0][0] = 1.0f * ascale;
  801   mwork[1][1] = 1.0f / ascale;
  802   mwork[2][2] = 1.0f;
  803 
  804   // moutput (of last calculation) -> minput
  805   MAT3SWAP(minput, moutput);
  806   // multiply mwork * minput -> moutput
  807   mat3mul((float *)moutput, (float *)mwork, (float *)minput);
  808 
  809 
  810   // Step 10: find x/y offsets and apply according correction so that
  811   // no negative coordinates occur in output vector
  812   float umin = FLT_MAX, vmin = FLT_MAX;
  813   // visit all four corners
  814   for(int y = 0; y < height; y += height - 1)
  815     for(int x = 0; x < width; x += width - 1)
  816     {
  817       float pi[3], po[3];
  818       pi[0] = x;
  819       pi[1] = y;
  820       pi[2] = 1.0f;
  821       // moutput expects input in (x:y:1) format and gives output as (x:y:1)
  822       mat3mulv(po, (float *)moutput, pi);
  823       umin = fmin(umin, po[0] / po[2]);
  824       vmin = fmin(vmin, po[1] / po[2]);
  825     }
  826 
  827   memset(mwork, 0, sizeof(float) * 9);
  828   mwork[0][0] = 1.0f;
  829   mwork[1][1] = 1.0f;
  830   mwork[2][2] = 1.0f;
  831   mwork[0][2] = -umin;
  832   mwork[1][2] = -vmin;
  833 
  834   // moutput (of last calculation) -> minput
  835   MAT3SWAP(minput, moutput);
  836   // multiply mwork * minput -> moutput
  837   mat3mul((float *)moutput, (float *)mwork, (float *)minput);
  838 
  839 
  840   // on request we either keep the final matrix for forward conversions
  841   // or produce an inverted matrix for backward conversions
  842   if(dir == ASHIFT_HOMOGRAPH_FORWARD)
  843   {
  844     // we have what we need -> copy it to the right place
  845     memcpy(homograph, moutput, sizeof(float) * 9);
  846   }
  847   else
  848   {
  849     // generate inverted homograph (mat3inv function defined in colorspaces.c)
  850     if(mat3inv((float *)homograph, (float *)moutput))
  851     {
  852       // in case of error we set to unity matrix
  853       memset(mwork, 0, sizeof(float) * 9);
  854       mwork[0][0] = 1.0f;
  855       mwork[1][1] = 1.0f;
  856       mwork[2][2] = 1.0f;
  857       memcpy(homograph, mwork, sizeof(float) * 9);
  858     }
  859   }
  860 }
  861 #undef MAT3SWAP
  862 
  863 
  864 // check if module parameters are set to all neutral values in which case the module's
  865 // output is identical to its input
  866 static inline int isneutral(const dt_iop_ashift_data_t *data)
  867 {
  868   // values lower than this have no visible effect
  869   const float eps = 1.0e-4f;
  870 
  871   return(fabs(data->rotation) < eps &&
  872          fabs(data->lensshift_v) < eps &&
  873          fabs(data->lensshift_h) < eps &&
  874          fabs(data->shear) < eps &&
  875          fabs(data->aspect - 1.0f) < eps &&
  876          data->cl < eps &&
  877          1.0f - data->cr < eps &&
  878          data->ct < eps &&
  879          1.0f - data->cb < eps);
  880 }
  881 
  882 
  883 int distort_transform(dt_iop_module_t *self, dt_dev_pixelpipe_iop_t *piece, float *const restrict points, size_t points_count)
  884 {
  885   const dt_iop_ashift_data_t *const data = (dt_iop_ashift_data_t *)piece->data;
  886 
  887   // nothing to be done if parameters are set to neutral values
  888   if(isneutral(data)) return 1;
  889 
  890   float DT_ALIGNED_ARRAY homograph[3][3];
  891   homography((float *)homograph, data->rotation, data->lensshift_v, data->lensshift_h, data->shear, data->f_length_kb,
  892              data->orthocorr, data->aspect, piece->buf_in.width, piece->buf_in.height, ASHIFT_HOMOGRAPH_FORWARD);
  893 
  894   // clipping offset
  895   const float fullwidth = (float)piece->buf_out.width / (data->cr - data->cl);
  896   const float fullheight = (float)piece->buf_out.height / (data->cb - data->ct);
  897   const float cx = fullwidth * data->cl;
  898   const float cy = fullheight * data->ct;
  899 
  900 #ifdef _OPENMP
  901 #pragma omp parallel for simd default(none) \
  902   dt_omp_firstprivate(cx, cy, points_count, points, homograph) \
  903   schedule(static) if(points_count > 100) aligned(points, homograph:64)
  904 #endif
  905   for(size_t i = 0; i < points_count * 2; i += 2)
  906   {
  907     float DT_ALIGNED_PIXEL pi[3] = { points[i], points[i + 1], 1.0f };
  908     float DT_ALIGNED_PIXEL po[3];
  909     mat3mulv(po, (float *)homograph, pi);
  910     points[i] = po[0] / po[2] - cx;
  911     points[i + 1] = po[1] / po[2] - cy;
  912   }
  913 
  914   return 1;
  915 }
  916 
  917 
  918 int distort_backtransform(dt_iop_module_t *self, dt_dev_pixelpipe_iop_t *piece, float *points,
  919                           size_t points_count)
  920 {
  921   const dt_iop_ashift_data_t *const data = (dt_iop_ashift_data_t *)piece->data;
  922 
  923   // nothing to be done if parameters are set to neutral values
  924   if(isneutral(data)) return 1;
  925 
  926   float DT_ALIGNED_ARRAY ihomograph[3][3];
  927   homography((float *)ihomograph, data->rotation, data->lensshift_v, data->lensshift_h, data->shear, data->f_length_kb,
  928              data->orthocorr, data->aspect, piece->buf_in.width, piece->buf_in.height, ASHIFT_HOMOGRAPH_INVERTED);
  929 
  930   // clipping offset
  931   const float fullwidth = (float)piece->buf_out.width / (data->cr - data->cl);
  932   const float fullheight = (float)piece->buf_out.height / (data->cb - data->ct);
  933   const float cx = fullwidth * data->cl;
  934   const float cy = fullheight * data->ct;
  935 
  936 #ifdef _OPENMP
  937 #pragma omp parallel for simd default(none) \
  938   dt_omp_firstprivate(cx, cy, points, points_count, ihomograph) \
  939   schedule(static) if(points_count > 100) aligned(ihomograph, points:64)
  940 #endif
  941   for(size_t i = 0; i < points_count * 2; i += 2)
  942   {
  943     float DT_ALIGNED_PIXEL pi[3] = { points[i] + cx, points[i + 1] + cy, 1.0f };
  944     float DT_ALIGNED_PIXEL po[3];
  945     mat3mulv(po, (float *)ihomograph, pi);
  946     points[i] = po[0] / po[2];
  947     points[i + 1] = po[1] / po[2];
  948   }
  949 
  950   return 1;
  951 }
  952 
  953 void distort_mask(struct dt_iop_module_t *self, struct dt_dev_pixelpipe_iop_t *piece, const float *const in,
  954                   float *const out, const dt_iop_roi_t *const roi_in, const dt_iop_roi_t *const roi_out)
  955 {
  956   const dt_iop_ashift_data_t *const data = (dt_iop_ashift_data_t *)piece->data;
  957 
  958   // if module is set to neutral parameters we just copy input->output and are done
  959   if(isneutral(data))
  960   {
  961     dt_iop_image_copy_by_size(out, in, roi_out->width, roi_out->height, 1);
  962     return;
  963   }
  964 
  965   const struct dt_interpolation *interpolation = dt_interpolation_new(DT_INTERPOLATION_USERPREF_WARP);
  966 
  967   float ihomograph[3][3];
  968   homography((float *)ihomograph, data->rotation, data->lensshift_v, data->lensshift_h, data->shear, data->f_length_kb,
  969              data->orthocorr, data->aspect, piece->buf_in.width, piece->buf_in.height, ASHIFT_HOMOGRAPH_INVERTED);
  970 
  971   // clipping offset
  972   const float fullwidth = (float)piece->buf_out.width / (data->cr - data->cl);
  973   const float fullheight = (float)piece->buf_out.height / (data->cb - data->ct);
  974   const float cx = roi_out->scale * fullwidth * data->cl;
  975   const float cy = roi_out->scale * fullheight * data->ct;
  976 
  977 #ifdef _OPENMP
  978 #pragma omp parallel for default(none) \
  979   dt_omp_firstprivate(cx, cy, in, out, roi_in, roi_out) \
  980   shared(ihomograph, interpolation) \
  981   schedule(static)
  982 #endif
  983   // go over all pixels of output image
  984   for(int j = 0; j < roi_out->height; j++)
  985   {
  986     float *const restrict _out = out + (size_t)j * roi_out->width;
  987     for(int i = 0; i < roi_out->width; i++)
  988     {
  989       float pin[3], pout[3];
  990 
  991       // convert output pixel coordinates to original image coordinates
  992       pout[0] = roi_out->x + i + cx;
  993       pout[1] = roi_out->y + j + cy;
  994       pout[0] /= roi_out->scale;
  995       pout[1] /= roi_out->scale;
  996       pout[2] = 1.0f;
  997 
  998       // apply homograph
  999       mat3mulv(pin, (float *)ihomograph, pout);
 1000 
 1001       // convert to input pixel coordinates
 1002       pin[0] /= pin[2];
 1003       pin[1] /= pin[2];
 1004       pin[0] *= roi_in->scale;
 1005       pin[1] *= roi_in->scale;
 1006       pin[0] -= roi_in->x;
 1007       pin[1] -= roi_in->y;
 1008 
 1009       // get output values by interpolation from input image
 1010       dt_interpolation_compute_pixel1c(interpolation, in, _out + i, pin[0], pin[1], roi_in->width,
 1011                                        roi_in->height, roi_in->width);
 1012     }
 1013   }
 1014 }
 1015 
 1016 void modify_roi_out(struct dt_iop_module_t *self, struct dt_dev_pixelpipe_iop_t *piece, dt_iop_roi_t *roi_out,
 1017                     const dt_iop_roi_t *roi_in)
 1018 {
 1019   dt_iop_ashift_data_t *data = (dt_iop_ashift_data_t *)piece->data;
 1020   *roi_out = *roi_in;
 1021 
 1022   // nothing more to be done if parameters are set to neutral values
 1023   if(isneutral(data)) return;
 1024 
 1025   float homograph[3][3];
 1026   homography((float *)homograph, data->rotation, data->lensshift_v, data->lensshift_h, data->shear, data->f_length_kb,
 1027              data->orthocorr, data->aspect, piece->buf_in.width, piece->buf_in.height, ASHIFT_HOMOGRAPH_FORWARD);
 1028 
 1029   float xm = FLT_MAX, xM = -FLT_MAX, ym = FLT_MAX, yM = -FLT_MAX;
 1030 
 1031   // go through all four vertices of input roi and convert coordinates to output
 1032   for(int y = 0; y < roi_in->height; y += roi_in->height - 1)
 1033   {
 1034     for(int x = 0; x < roi_in->width; x += roi_in->width - 1)
 1035     {
 1036       float pin[3], pout[3];
 1037 
 1038       // convert from input coordinates to original image coordinates
 1039       pin[0] = roi_in->x + x;
 1040       pin[1] = roi_in->y + y;
 1041       pin[0] /= roi_in->scale;
 1042       pin[1] /= roi_in->scale;
 1043       pin[2] = 1.0f;
 1044 
 1045       // apply homograph
 1046       mat3mulv(pout, (float *)homograph, pin);
 1047 
 1048       // convert to output image coordinates
 1049       pout[0] /= pout[2];
 1050       pout[1] /= pout[2];
 1051       pout[0] *= roi_out->scale;
 1052       pout[1] *= roi_out->scale;
 1053       xm = MIN(xm, pout[0]);
 1054       xM = MAX(xM, pout[0]);
 1055       ym = MIN(ym, pout[1]);
 1056       yM = MAX(yM, pout[1]);
 1057     }
 1058   }
 1059 
 1060   float width = xM - xm + 1;
 1061   float height = yM - ym + 1;
 1062 
 1063   // clipping adjustments
 1064   width *= data->cr - data->cl;
 1065   height *= data->cb - data->ct;
 1066 
 1067   roi_out->width = floorf(width);
 1068   roi_out->height = floorf(height);
 1069 
 1070 #ifdef ASHIFT_DEBUG
 1071   print_roi(roi_in, "roi_in (going into modify_roi_out)");
 1072   print_roi(roi_out, "roi_out (after modify_roi_out)");
 1073 #endif
 1074 }
 1075 
 1076 void modify_roi_in(struct dt_iop_module_t *self, struct dt_dev_pixelpipe_iop_t *piece,
 1077                    const dt_iop_roi_t *const roi_out, dt_iop_roi_t *roi_in)
 1078 {
 1079   dt_iop_ashift_data_t *data = (dt_iop_ashift_data_t *)piece->data;
 1080   *roi_in = *roi_out;
 1081 
 1082   // nothing more to be done if parameters are set to neutral values
 1083   if(isneutral(data)) return;
 1084 
 1085   float ihomograph[3][3];
 1086   homography((float *)ihomograph, data->rotation, data->lensshift_v, data->lensshift_h, data->shear, data->f_length_kb,
 1087              data->orthocorr, data->aspect, piece->buf_in.width, piece->buf_in.height, ASHIFT_HOMOGRAPH_INVERTED);
 1088 
 1089   const float orig_w = roi_in->scale * piece->buf_in.width;
 1090   const float orig_h = roi_in->scale * piece->buf_in.height;
 1091 
 1092   // clipping offset
 1093   const float fullwidth = (float)piece->buf_out.width / (data->cr - data->cl);
 1094   const float fullheight = (float)piece->buf_out.height / (data->cb - data->ct);
 1095   const float cx = roi_out->scale * fullwidth * data->cl;
 1096   const float cy = roi_out->scale * fullheight * data->ct;
 1097 
 1098   float xm = FLT_MAX, xM = -FLT_MAX, ym = FLT_MAX, yM = -FLT_MAX;
 1099 
 1100   // go through all four vertices of output roi and convert coordinates to input
 1101   for(int y = 0; y < roi_out->height; y += roi_out->height - 1)
 1102   {
 1103     for(int x = 0; x < roi_out->width; x += roi_out->width - 1)
 1104     {
 1105       float pin[3], pout[3];
 1106 
 1107       // convert from output image coordinates to original image coordinates
 1108       pout[0] = roi_out->x + x + cx;
 1109       pout[1] = roi_out->y + y + cy;
 1110       pout[0] /= roi_out->scale;
 1111       pout[1] /= roi_out->scale;
 1112       pout[2] = 1.0f;
 1113 
 1114       // apply homograph
 1115       mat3mulv(pin, (float *)ihomograph, pout);
 1116 
 1117       // convert to input image coordinates
 1118       pin[0] /= pin[2];
 1119       pin[1] /= pin[2];
 1120       pin[0] *= roi_in->scale;
 1121       pin[1] *= roi_in->scale;
 1122       xm = MIN(xm, pin[0]);
 1123       xM = MAX(xM, pin[0]);
 1124       ym = MIN(ym, pin[1]);
 1125       yM = MAX(yM, pin[1]);
 1126     }
 1127   }
 1128 
 1129   const struct dt_interpolation *interpolation = dt_interpolation_new(DT_INTERPOLATION_USERPREF_WARP);
 1130   roi_in->x = fmaxf(0.0f, xm - interpolation->width);
 1131   roi_in->y = fmaxf(0.0f, ym - interpolation->width);
 1132   roi_in->width = fminf(ceilf(orig_w) - roi_in->x, xM - roi_in->x + 1 + interpolation->width);
 1133   roi_in->height = fminf(ceilf(orig_h) - roi_in->y, yM - roi_in->y + 1 + interpolation->width);
 1134 
 1135   // sanity check.
 1136   roi_in->x = CLAMP(roi_in->x, 0, (int)floorf(orig_w));
 1137   roi_in->y = CLAMP(roi_in->y, 0, (int)floorf(orig_h));
 1138   roi_in->width = CLAMP(roi_in->width, 1, (int)floorf(orig_w) - roi_in->x);
 1139   roi_in->height = CLAMP(roi_in->height, 1, (int)floorf(orig_h) - roi_in->y);
 1140 #ifdef ASHIFT_DEBUG
 1141   print_roi(roi_out, "roi_out (going into modify_roi_in)");
 1142   print_roi(roi_in, "roi_in (after modify_roi_in)");
 1143 #endif
 1144 }
 1145 
 1146 // simple conversion of rgb image into greyscale variant suitable for line segment detection
 1147 // the lsd routines expect input as *double, roughly in the range [0.0; 256.0]
 1148 static void rgb2grey256(const float *const in, double *const out, const int width, const int height)
 1149 {
 1150   const size_t npixels = (size_t)width * height;
 1151 
 1152 #ifdef _OPENMP
 1153 #pragma omp parallel for default(none) \
 1154   dt_omp_firstprivate(npixels) \
 1155   dt_omp_sharedconst(in, out) \
 1156   schedule(static)
 1157 #endif
 1158   for(int index = 0; index < npixels; index++)
 1159   {
 1160     out[index] = (0.3f * in[4*index+0] + 0.59f * in[4*index+1] + 0.11f * in[4*index+2]) * 256.0;
 1161   }
 1162 }
 1163 
 1164 // sobel edge enhancement in one direction
 1165 static void edge_enhance_1d(const double *in, double *out, const int width, const int height,
 1166                             dt_iop_ashift_enhance_t dir)
 1167 {
 1168   // Sobel kernels for both directions
 1169   const double hkernel[3][3] = { { 1.0, 0.0, -1.0 }, { 2.0, 0.0, -2.0 }, { 1.0, 0.0, -1.0 } };
 1170   const double vkernel[3][3] = { { 1.0, 2.0, 1.0 }, { 0.0, 0.0, 0.0 }, { -1.0, -2.0, -1.0 } };
 1171   const int kwidth = 3;
 1172   const int khwidth = kwidth / 2;
 1173 
 1174   // select kernel
 1175   const double *kernel = (dir == ASHIFT_ENHANCE_HORIZONTAL) ? (const double *)hkernel : (const double *)vkernel;
 1176 
 1177 #ifdef _OPENMP
 1178 #pragma omp parallel for default(none) \
 1179   dt_omp_firstprivate(height, width, khwidth, kwidth) \
 1180   shared(in, out, kernel) \
 1181   schedule(static)
 1182 #endif
 1183   // loop over image pixels and perform sobel convolution
 1184   for(int j = khwidth; j < height - khwidth; j++)
 1185   {
 1186     const double *inp = in + (size_t)j * width + khwidth;
 1187     double *outp = out + (size_t)j * width + khwidth;
 1188     for(int i = khwidth; i < width - khwidth; i++, inp++, outp++)
 1189     {
 1190       double sum = 0.0f;
 1191       for(int jj = 0; jj < kwidth; jj++)
 1192       {
 1193         const int k = jj * kwidth;
 1194         const int l = (jj - khwidth) * width;
 1195         for(int ii = 0; ii < kwidth; ii++)
 1196         {
 1197           sum += inp[l + ii - khwidth] * kernel[k + ii];
 1198         }
 1199       }
 1200       *outp = sum;
 1201     }
 1202   }
 1203 
 1204 #ifdef _OPENMP
 1205 #pragma omp parallel for default(none) \
 1206   dt_omp_firstprivate(height, width, khwidth) \
 1207   shared(out) \
 1208   schedule(static)
 1209 #endif
 1210   // border fill in output buffer, so we don't get pseudo lines at image frame
 1211   for(int j = 0; j < height; j++)
 1212     for(int i = 0; i < width; i++)
 1213     {
 1214       double val = out[j * width + i];
 1215 
 1216       if(j < khwidth)
 1217         val = out[(khwidth - j) * width + i];
 1218       else if(j >= height - khwidth)
 1219         val = out[(j - khwidth) * width + i];
 1220       else if(i < khwidth)
 1221         val = out[j * width + (khwidth - i)];
 1222       else if(i >= width - khwidth)
 1223         val = out[j * width + (i - khwidth)];
 1224 
 1225       out[j * width + i] = val;
 1226 
 1227       // jump over center of image
 1228       if(i == khwidth && j >= khwidth && j < height - khwidth) i = width - khwidth;
 1229     }
 1230 }
 1231 
 1232 // edge enhancement in both directions
 1233 static int edge_enhance(const double *in, double *out, const int width, const int height)
 1234 {
 1235   double *Gx = NULL;
 1236   double *Gy = NULL;
 1237 
 1238   Gx = malloc(sizeof(double) * width * height);
 1239   if(Gx == NULL) goto error;
 1240 
 1241   Gy = malloc(sizeof(double) * width * height);
 1242   if(Gy == NULL) goto error;
 1243 
 1244   // perform edge enhancement in both directions
 1245   edge_enhance_1d(in, Gx, width, height, ASHIFT_ENHANCE_HORIZONTAL);
 1246   edge_enhance_1d(in, Gy, width, height, ASHIFT_ENHANCE_VERTICAL);
 1247 
 1248 // calculate absolute values
 1249 #ifdef _OPENMP
 1250 #pragma omp parallel for default(none) \
 1251   dt_omp_firstprivate(height, width) \
 1252   shared(Gx, Gy, out) \
 1253   schedule(static)
 1254 #endif
 1255   for(size_t k = 0; k < (size_t)width * height; k++)
 1256   {
 1257     out[k] = sqrt(Gx[k] * Gx[k] + Gy[k] * Gy[k]);
 1258   }
 1259 
 1260   free(Gx);
 1261   free(Gy);
 1262   return TRUE;
 1263 
 1264 error:
 1265   if(Gx) free(Gx);
 1266   if(Gy) free(Gy);
 1267   return FALSE;
 1268 }
 1269 
 1270 // XYZ -> sRGB matrix
 1271 static void XYZ_to_sRGB(const float *XYZ, float *sRGB)
 1272 {
 1273   sRGB[0] =  3.1338561f * XYZ[0] - 1.6168667f * XYZ[1] - 0.4906146f * XYZ[2];
 1274   sRGB[1] = -0.9787684f * XYZ[0] + 1.9161415f * XYZ[1] + 0.0334540f * XYZ[2];
 1275   sRGB[2] =  0.0719453f * XYZ[0] - 0.2289914f * XYZ[1] + 1.4052427f * XYZ[2];
 1276 }
 1277 
 1278 // sRGB -> XYZ matrix
 1279 static void sRGB_to_XYZ(const float *sRGB, float *XYZ)
 1280 {
 1281   XYZ[0] = 0.4360747f * sRGB[0] + 0.3850649f * sRGB[1] + 0.1430804f * sRGB[2];
 1282   XYZ[1] = 0.2225045f * sRGB[0] + 0.7168786f * sRGB[1] + 0.0606169f * sRGB[2];
 1283   XYZ[2] = 0.0139322f * sRGB[0] + 0.0971045f * sRGB[1] + 0.7141733f * sRGB[2];
 1284 }
 1285 
 1286 // detail enhancement via bilateral grid (function arguments in and out may represent identical buffers)
 1287 static int detail_enhance(const float *const in, float *const out, const int width, const int height)
 1288 {
 1289   const float sigma_r = 5.0f;
 1290   const float sigma_s = fminf(width, height) * 0.02f;
 1291   const float detail = 10.0f;
 1292   const size_t npixels = (size_t)width * height;
 1293   int success = TRUE;
 1294 
 1295   // we need to convert from RGB to Lab first;
 1296   // as colors don't matter we are safe to assume data to be sRGB
 1297 
 1298   // convert RGB input to Lab, use output buffer for intermediate storage
 1299 #ifdef _OPENMP
 1300 #pragma omp parallel for default(none) \
 1301   dt_omp_firstprivate(npixels) \
 1302   dt_omp_sharedconst(in, out) \
 1303   schedule(static)
 1304 #endif
 1305   for(size_t index = 0; index < 4*npixels; index += 4)
 1306   {
 1307     float DT_ALIGNED_PIXEL XYZ[4];
 1308     sRGB_to_XYZ(in + index, XYZ);
 1309     dt_XYZ_to_Lab(XYZ, out + index);
 1310   }
 1311 
 1312   // bilateral grid detail enhancement
 1313   dt_bilateral_t *b = dt_bilateral_init(width, height, sigma_s, sigma_r);
 1314 
 1315   if(b != NULL)
 1316   {
 1317     dt_bilateral_splat(b, out);
 1318     dt_bilateral_blur(b);
 1319     dt_bilateral_slice_to_output(b, out, out, detail);
 1320     dt_bilateral_free(b);
 1321   }
 1322   else
 1323     success = FALSE;
 1324 
 1325   // convert resulting Lab to RGB output
 1326 #ifdef _OPENMP
 1327 #pragma omp parallel for default(none) \
 1328   dt_omp_firstprivate(npixels) \
 1329   dt_omp_sharedconst(out) \
 1330   schedule(static)
 1331 #endif
 1332   for(size_t index = 0; index < 4*npixels; index += 4)
 1333   {
 1334     float XYZ[3];
 1335     dt_Lab_to_XYZ(out + index, XYZ);
 1336     XYZ_to_sRGB(XYZ, out + index);
 1337   }
 1338 
 1339   return success;
 1340 }
 1341 
 1342 // apply gamma correction to RGB buffer (function arguments in and out may represent identical buffers)
 1343 static void gamma_correct(const float *const in, float *const out, const int width, const int height)
 1344 {
 1345   const size_t npixels = (size_t)width * height;
 1346 #ifdef _OPENMP
 1347 #pragma omp parallel for default(none) \
 1348   dt_omp_firstprivate(npixels) \
 1349   dt_omp_sharedconst(in, out) \
 1350   schedule(static)
 1351 #endif
 1352   for(int index = 0; index < 4*npixels; index += 4)
 1353   {
 1354     for(int c = 0; c < 3; c++)
 1355       out[index+c] = powf(in[index+c], LSD_GAMMA);
 1356   }
 1357 }
 1358 
 1359 // do actual line_detection based on LSD algorithm and return results according
 1360 // to this module's conventions
 1361 static int line_detect(float *in, const int width, const int height, const int x_off, const int y_off,
 1362                        const float scale, dt_iop_ashift_line_t **alines, int *lcount, int *vcount, int *hcount,
 1363                        float *vweight, float *hweight, dt_iop_ashift_enhance_t enhance, const int is_raw)
 1364 {
 1365   double *greyscale = NULL;
 1366   double *lsd_lines = NULL;
 1367   dt_iop_ashift_line_t *ashift_lines = NULL;
 1368 
 1369   int vertical_count = 0;
 1370   int horizontal_count = 0;
 1371   float vertical_weight = 0.0f;
 1372   float horizontal_weight = 0.0f;
 1373 
 1374   // apply gamma correction if image is raw
 1375   if(is_raw)
 1376   {
 1377     gamma_correct(in, in, width, height);
 1378   }
 1379 
 1380   // if requested perform an additional detail enhancement step
 1381   if(enhance & ASHIFT_ENHANCE_DETAIL)
 1382   {
 1383     (void)detail_enhance(in, in, width, height);
 1384   }
 1385 
 1386   // allocate intermediate buffers
 1387   greyscale = malloc(sizeof(double) * width * height);
 1388   if(greyscale == NULL) goto error;
 1389 
 1390   // convert to greyscale image
 1391   rgb2grey256(in, greyscale, width, height);
 1392 
 1393   // if requested perform an additional edge enhancement step
 1394   if(enhance & ASHIFT_ENHANCE_EDGES)
 1395   {
 1396     (void)edge_enhance(greyscale, greyscale, width, height);
 1397   }
 1398 
 1399   // call the line segment detector LSD;
 1400   // LSD stores the number of found lines in lines_count.
 1401   // it returns structural details as vector 'double lines[7 * lines_count]'
 1402   int lines_count;
 1403 
 1404   lsd_lines = LineSegmentDetection(&lines_count, greyscale, width, height,
 1405                                    LSD_SCALE, LSD_SIGMA_SCALE, LSD_QUANT,
 1406                                    LSD_ANG_TH, LSD_LOG_EPS, LSD_DENSITY_TH,
 1407                                    LSD_N_BINS, NULL, NULL, NULL);
 1408 
 1409   // we count the lines that we really want to use
 1410   int lct = 0;
 1411   if(lines_count > 0)
 1412   {
 1413     // aggregate lines data into our own structures
 1414     ashift_lines = (dt_iop_ashift_line_t *)malloc(sizeof(dt_iop_ashift_line_t) * lines_count);
 1415     if(ashift_lines == NULL) goto error;
 1416 
 1417     for(int n = 0; n < lines_count; n++)
 1418     {
 1419       float x1 = lsd_lines[n * 7 + 0];
 1420       float y1 = lsd_lines[n * 7 + 1];
 1421       float x2 = lsd_lines[n * 7 + 2];
 1422       float y2 = lsd_lines[n * 7 + 3];
 1423 
 1424       // check for lines running along image borders and skip them.
 1425       // these would likely be false-positives which could result
 1426       // from any kind of processing artifacts
 1427       if((fabsf(x1 - x2) < 1 && fmaxf(x1, x2) < 2) ||
 1428          (fabsf(x1 - x2) < 1 && fminf(x1, x2) > width - 3) ||
 1429          (fabsf(y1 - y2) < 1 && fmaxf(y1, y2) < 2) ||
 1430          (fabsf(y1 - y2) < 1 && fminf(y1, y2) > height - 3))
 1431         continue;
 1432 
 1433       // line position in absolute coordinates
 1434       float px1 = x_off + x1;
 1435       float py1 = y_off + y1;
 1436       float px2 = x_off + x2;
 1437       float py2 = y_off + y2;
 1438 
 1439       // scale back to input buffer
 1440       px1 /= scale;
 1441       py1 /= scale;
 1442       px2 /= scale;
 1443       py2 /= scale;
 1444 
 1445       // store as homogeneous coordinates
 1446       ashift_lines[lct].p1[0] = px1;
 1447       ashift_lines[lct].p1[1] = py1;
 1448       ashift_lines[lct].p1[2] = 1.0f;
 1449       ashift_lines[lct].p2[0] = px2;
 1450       ashift_lines[lct].p2[1] = py2;
 1451       ashift_lines[lct].p2[2] = 1.0f;
 1452 
 1453       // calculate homogeneous coordinates of connecting line (defined by the two points)
 1454       vec3prodn(ashift_lines[lct].L, ashift_lines[lct].p1, ashift_lines[lct].p2);
 1455 
 1456       // normalaze line coordinates so that x^2 + y^2 = 1
 1457       // (this will always succeed as L is a real line connecting two real points)
 1458       vec3lnorm(ashift_lines[lct].L, ashift_lines[lct].L);
 1459 
 1460       // length and width of rectangle (see LSD)
 1461       ashift_lines[lct].length = sqrt((px2 - px1) * (px2 - px1) + (py2 - py1) * (py2 - py1));
 1462       ashift_lines[lct].width = lsd_lines[n * 7 + 4] / scale;
 1463 
 1464       // ...  and weight (= length * width * angle precision)
 1465       const float weight = ashift_lines[lct].length * ashift_lines[lct].width * lsd_lines[n * 7 + 5];
 1466       ashift_lines[lct].weight = weight;
 1467 
 1468 
 1469       const float angle = atan2f(py2 - py1, px2 - px1) / M_PI * 180.0f;
 1470       const int vertical = fabsf(fabsf(angle) - 90.0f) < MAX_TANGENTIAL_DEVIATION ? 1 : 0;
 1471       const int horizontal = fabsf(fabsf(fabsf(angle) - 90.0f) - 90.0f) < MAX_TANGENTIAL_DEVIATION ? 1 : 0;
 1472 
 1473       const int relevant = ashift_lines[lct].length > MIN_LINE_LENGTH ? 1 : 0;
 1474 
 1475       // register type of line
 1476       dt_iop_ashift_linetype_t type = ASHIFT_LINE_IRRELEVANT;
 1477       if(vertical && relevant)
 1478       {
 1479         type = ASHIFT_LINE_VERTICAL_SELECTED;
 1480         vertical_count++;
 1481         vertical_weight += weight;
 1482       }
 1483       else if(horizontal && relevant)
 1484       {
 1485         type = ASHIFT_LINE_HORIZONTAL_SELECTED;
 1486         horizontal_count++;
 1487         horizontal_weight += weight;
 1488       }
 1489       ashift_lines[lct].type = type;
 1490 
 1491       // the next valid line
 1492       lct++;
 1493     }
 1494   }
 1495 #ifdef ASHIFT_DEBUG
 1496     printf("%d lines (vertical %d, horizontal %d, not relevant %d)\n", lines_count, vertical_count,
 1497            horizontal_count, lct - vertical_count - horizontal_count);
 1498     float xmin = FLT_MAX, xmax = FLT_MIN, ymin = FLT_MAX, ymax = FLT_MIN;
 1499     for(int n = 0; n < lct; n++)
 1500     {
 1501       xmin = fmin(xmin, fmin(ashift_lines[n].p1[0], ashift_lines[n].p2[0]));
 1502       xmax = fmax(xmax, fmax(ashift_lines[n].p1[0], ashift_lines[n].p2[0]));
 1503       ymin = fmin(ymin, fmin(ashift_lines[n].p1[1], ashift_lines[n].p2[1]));
 1504       ymax = fmax(ymax, fmax(ashift_lines[n].p1[1], ashift_lines[n].p2[1]));
 1505       printf("x1 %.0f, y1 %.0f, x2 %.0f, y2 %.0f, length %.0f, width %f, X %f, Y %f, Z %f, type %d, scalars %f %f\n",
 1506              ashift_lines[n].p1[0], ashift_lines[n].p1[1], ashift_lines[n].p2[0], ashift_lines[n].p2[1],
 1507              ashift_lines[n].length, ashift_lines[n].width,
 1508              ashift_lines[n].L[0], ashift_lines[n].L[1], ashift_lines[n].L[2], ashift_lines[n].type,
 1509              vec3scalar(ashift_lines[n].p1, ashift_lines[n].L),
 1510              vec3scalar(ashift_lines[n].p2, ashift_lines[n].L));
 1511     }
 1512     printf("xmin %.0f, xmax %.0f, ymin %.0f, ymax %.0f\n", xmin, xmax, ymin, ymax);
 1513 #endif
 1514 
 1515   // store results in provided locations
 1516   *lcount = lct;
 1517   *vcount = vertical_count;
 1518   *vweight = vertical_weight;
 1519   *hcount = horizontal_count;
 1520   *hweight = horizontal_weight;
 1521   *alines = ashift_lines;
 1522 
 1523   // free intermediate buffers
 1524   free(lsd_lines);
 1525   free(greyscale);
 1526   return lct > 0 ? TRUE : FALSE;
 1527 
 1528 error:
 1529   free(lsd_lines);
 1530   free(greyscale);
 1531   return FALSE;
 1532 }
 1533 
 1534 // get image from buffer, analyze for structure and save results
 1535 static int get_structure(dt_iop_module_t *module, dt_iop_ashift_enhance_t enhance)
 1536 {
 1537   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)module->gui_data;
 1538 
 1539   float *buffer = NULL;
 1540   int width = 0;
 1541   int height = 0;
 1542   int x_off = 0;
 1543   int y_off = 0;
 1544   float scale = 0.0f;
 1545 
 1546   dt_iop_gui_enter_critical_section(module);
 1547   // read buffer data if they are available
 1548   if(g->buf != NULL)
 1549   {
 1550     width = g->buf_width;
 1551     height = g->buf_height;
 1552     x_off = g->buf_x_off;
 1553     y_off = g->buf_y_off;
 1554     scale = g->buf_scale;
 1555 
 1556     // create a temporary buffer to hold image data
 1557     buffer = malloc(sizeof(float) * 4 * (size_t)width * height);
 1558     if(buffer != NULL)
 1559       dt_iop_image_copy_by_size(buffer, g->buf, width, height, 4);
 1560   }
 1561   dt_iop_gui_leave_critical_section(module);
 1562 
 1563   if(buffer == NULL) goto error;
 1564 
 1565   // get rid of old structural data
 1566   g->lines_count = 0;
 1567   g->vertical_count = 0;
 1568   g->horizontal_count = 0;
 1569   free(g->lines);
 1570   g->lines = NULL;
 1571 
 1572   dt_iop_ashift_line_t *lines;
 1573   int lines_count;
 1574   int vertical_count;
 1575   int horizontal_count;
 1576   float vertical_weight;
 1577   float horizontal_weight;
 1578 
 1579   // get new structural data
 1580   if(!line_detect(buffer, width, height, x_off, y_off, scale, &lines, &lines_count,
 1581                   &vertical_count, &horizontal_count, &vertical_weight, &horizontal_weight,
 1582                   enhance, dt_image_is_raw(&module->dev->image_storage)))
 1583     goto error;
 1584 
 1585   // save new structural data
 1586   g->lines_in_width = width;
 1587   g->lines_in_height = height;
 1588   g->lines_x_off = x_off;
 1589   g->lines_y_off = y_off;
 1590   g->lines_count = lines_count;
 1591   g->vertical_count = vertical_count;
 1592   g->horizontal_count = horizontal_count;
 1593   g->vertical_weight = vertical_weight;
 1594   g->horizontal_weight = horizontal_weight;
 1595   g->lines_version++;
 1596   g->lines_suppressed = 0;
 1597   g->lines = lines;
 1598 
 1599   free(buffer);
 1600   return TRUE;
 1601 
 1602 error:
 1603   free(buffer);
 1604   return FALSE;
 1605 }
 1606 
 1607 
 1608 // swap two integer values
 1609 static inline void swap(int *a, int *b)
 1610 {
 1611   int tmp = *a;
 1612   *a = *b;
 1613   *b = tmp;
 1614 }
 1615 
 1616 // do complete permutations
 1617 static int quickperm(int *a, int *p, const int N, int *i)
 1618 {
 1619   if(*i >= N) return FALSE;
 1620 
 1621   p[*i]--;
 1622   int j = (*i % 2 == 1) ? p[*i] : 0;
 1623   swap(&a[j], &a[*i]);
 1624   *i = 1;
 1625   while(p[*i] == 0)
 1626   {
 1627     p[*i] = *i;
 1628     (*i)++;
 1629   }
 1630   return TRUE;
 1631 }
 1632 
 1633 // Fisher-Yates shuffle
 1634 static void shuffle(int *a, const int N)
 1635 {
 1636   for(int i = 0; i < N; i++)
 1637   {
 1638     int j = i + rand() % (N - i);
 1639     swap(&a[j], &a[i]);
 1640   }
 1641 }
 1642 
 1643 // factorial function
 1644 static int fact(const int n)
 1645 {
 1646   return (n == 1 ? 1 : n * fact(n - 1));
 1647 }
 1648 
 1649 // We use a pseudo-RANSAC algorithm to elminiate ouliers from our set of lines. The
 1650 // original RANSAC works on linear optimization problems. Our model is nonlinear. We
 1651 // take advantage of the fact that lines interesting for our model are vantage lines
 1652 // that meet in one vantage point for each subset of lines (vertical/horizontal).
 1653 // Strategy: we construct a model by (random) sampling within the subset of lines and
 1654 // calculate the vantage point. Then we check the "distance" of all other lines to the
 1655 // vantage point. The model that gives highest number of lines combined with the highest
 1656 // total weight and lowest overall "distance" wins.
 1657 // Disadvantage: compared to the original RANSAC we don't get any model parameters that
 1658 // we could use for the following NMS fit.
 1659 // Self-tuning: we optimize "epsilon", the hurdle rate to reject a line as an outlier,
 1660 // by a number of dry runs first. The target average percentage value of lines to eliminate as
 1661 // outliers (without judging on the quality of the model) is given by RANSAC_ELIMINATION_RATIO,
 1662 // note: the actual percentage of outliers removed in the final run will be lower because we
 1663 // will finally look for the best quality model with the optimized epsilon and that quality value also
 1664 // encloses the number of good lines
 1665 static void ransac(const dt_iop_ashift_line_t *lines, int *index_set, int *inout_set,
 1666                   const int set_count, const float total_weight, const int xmin, const int xmax,
 1667                   const int ymin, const int ymax)
 1668 {
 1669   if(set_count < 3) return;
 1670 
 1671   const size_t set_size = set_count * sizeof(int);
 1672   int *best_set = malloc(set_size);
 1673   memcpy(best_set, index_set, set_size);
 1674   int *best_inout = calloc(1, set_size);
 1675 
 1676   float best_quality = 0.0f;
 1677 
 1678   // hurdle value epsilon for rejecting a line as an outlier will be self-tuning
 1679   // in a number of dry runs
 1680   float epsilon = powf(10.0f, -RANSAC_EPSILON);
 1681   float epsilon_step = RANSAC_EPSILON_STEP;
 1682   // some accounting variables for self-tuning
 1683   int lines_eliminated = 0;
 1684   int valid_runs = 0;
 1685 
 1686   // number of runs to optimize epsilon
 1687   const int optiruns = RANSAC_OPTIMIZATION_STEPS * RANSAC_OPTIMIZATION_DRY_RUNS;
 1688   // go for complete permutations on small set sizes, else for random sample consensus
 1689   const int riter = (set_count > RANSAC_HURDLE) ? RANSAC_RUNS : fact(set_count);
 1690 
 1691   // some data needed for quickperm
 1692   int *perm = malloc(sizeof(int) * (set_count + 1));
 1693   for(int n = 0; n < set_count + 1; n++) perm[n] = n;
 1694   int piter = 1;
 1695 
 1696   // inout holds good/bad qualification for each line
 1697   int *inout = malloc(set_size);
 1698 
 1699   for(int r = 0; r < optiruns + riter; r++)
 1700   {
 1701     // get random or systematic variation of index set
 1702     if(set_count > RANSAC_HURDLE || r < optiruns)
 1703       shuffle(index_set, set_count);
 1704     else
 1705       (void)quickperm(index_set, perm, set_count, &piter);
 1706 
 1707     // summed quality evaluation of this run
 1708     float quality = 0.0f;
 1709 
 1710     // we build a model ouf of the first two lines
 1711     const float *L1 = lines[index_set[0]].L;
 1712     const float *L2 = lines[index_set[1]].L;
 1713 
 1714     // get intersection point (ideally a vantage point)
 1715     float V[3];
 1716     vec3prodn(V, L1, L2);
 1717 
 1718     // catch special cases:
 1719     // a) L1 and L2 are identical -> V is NULL -> no valid vantage point
 1720     // b) vantage point lies inside image frame (no chance to correct for this case)
 1721     if(vec3isnull(V) ||
 1722        (fabsf(V[2]) > 0.0f &&
 1723         V[0]/V[2] >= xmin &&
 1724         V[1]/V[2] >= ymin &&
 1725         V[0]/V[2] <= xmax &&
 1726         V[1]/V[2] <= ymax))
 1727     {
 1728       // no valid model
 1729       quality = 0.0f;
 1730     }
 1731     else
 1732     {
 1733       // valid model
 1734 
 1735       // normalize V so that x^2 + y^2 + z^2 = 1
 1736       vec3norm(V, V);
 1737 
 1738       // the two lines constituting the model are part of the set
 1739       inout[0] = 1;
 1740       inout[1] = 1;
 1741 
 1742       // go through all remaining lines, check if they are within the model, and
 1743       // mark that fact in inout[].
 1744       // summarize a quality parameter for all lines within the model
 1745       for(int n = 2; n < set_count; n++)
 1746       {
 1747         // L is normalized so that x^2 + y^2 = 1
 1748         const float *L3 = lines[index_set[n]].L;
 1749 
 1750         // we take the absolute value of the dot product of V and L as a measure
 1751         // of the "distance" between point and line. Note that this is not the real euclidean
 1752         // distance but - with the given normalization - just a pragmatically selected number
 1753         // that goes to zero if V lies on L and increases the more V and L are apart
 1754         const float d = fabsf(vec3scalar(V, L3));
 1755 
 1756         // depending on d we either include or exclude the point from the set
 1757         inout[n] = (d < epsilon) ? 1 : 0;
 1758 
 1759         float q;
 1760 
 1761         if(inout[n] == 1)
 1762         {
 1763           // a quality parameter that depends 1/3 on the number of lines within the model,
 1764           // 1/3 on their weight, and 1/3 on their weighted distance d to the vantage point
 1765           q = 0.33f / (float)set_count
 1766               + 0.33f * lines[index_set[n]].weight / total_weight
 1767               + 0.33f * (1.0f - d / epsilon) * (float)set_count * lines[index_set[n]].weight / total_weight;
 1768         }
 1769         else
 1770         {
 1771           q = 0.0f;
 1772           lines_eliminated++;
 1773         }
 1774 
 1775         quality += q;
 1776       }
 1777       valid_runs++;
 1778     }
 1779 
 1780     if(r < optiruns)
 1781     {
 1782       // on last run of each self-tuning step
 1783       if((r % RANSAC_OPTIMIZATION_DRY_RUNS) == (RANSAC_OPTIMIZATION_DRY_RUNS - 1) && (valid_runs > 0))
 1784       {
 1785 #ifdef ASHIFT_DEBUG
 1786         printf("ransac self-tuning (run %d): epsilon %f", r, epsilon);
 1787 #endif
 1788         // average ratio of lines that we eliminated with the given epsilon
 1789         float ratio = 100.0f * (float)lines_eliminated / ((float)set_count * valid_runs);
 1790         // adjust epsilon accordingly
 1791         if(ratio < RANSAC_ELIMINATION_RATIO)
 1792           epsilon = powf(10.0f, log10(epsilon) - epsilon_step);
 1793         else if(ratio > RANSAC_ELIMINATION_RATIO)
 1794           epsilon = powf(10.0f, log10(epsilon) + epsilon_step);
 1795 #ifdef ASHIFT_DEBUG
 1796         printf(" (elimination ratio %f) -> %f\n", ratio, epsilon);
 1797 #endif
 1798         // reduce step-size for next optimization round
 1799         epsilon_step /= 2.0f;
 1800         lines_eliminated = 0;
 1801         valid_runs = 0;
 1802       }
 1803     }
 1804     else
 1805     {
 1806       // in the "real" runs check against the best model found so far
 1807       if(quality > best_quality)
 1808       {
 1809         memcpy(best_set, index_set, set_size);
 1810         memcpy(best_inout, inout, set_size);
 1811         best_quality = quality;
 1812       }
 1813     }
 1814 
 1815 #ifdef ASHIFT_DEBUG
 1816     // report some statistics
 1817     int count = 0, lastcount = 0;
 1818     for(int n = 0; n < set_count; n++) count += best_inout[n];
 1819     for(int n = 0; n < set_count; n++) lastcount += inout[n];
 1820     printf("ransac run %d: best qual %.6f, eps %.6f, line count %d of %d (this run: qual %.5f, count %d (%2f%%))\n", r,
 1821            best_quality, epsilon, count, set_count, quality, lastcount, 100.0f * lastcount / (float)set_count);
 1822 #endif
 1823   }
 1824 
 1825   // store back best set
 1826   memcpy(index_set, best_set, set_size);
 1827   memcpy(inout_set, best_inout, set_size);
 1828 
 1829   free(inout);
 1830   free(perm);
 1831   free(best_inout);
 1832   free(best_set);
 1833 }
 1834 
 1835 
 1836 // try to clean up structural data by eliminating outliers and thereby increasing
 1837 // the chance of a convergent fitting
 1838 static int remove_outliers(dt_iop_module_t *module)
 1839 {
 1840   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)module->gui_data;
 1841 
 1842   const int width = g->lines_in_width;
 1843   const int height = g->lines_in_height;
 1844   const int xmin = g->lines_x_off;
 1845   const int ymin = g->lines_y_off;
 1846   const int xmax = xmin + width;
 1847   const int ymax = ymin + height;
 1848 
 1849   // holds the index set of lines we want to work on
 1850   int *lines_set = malloc(sizeof(int) * g->lines_count);
 1851   // holds the result of ransac
 1852   int *inout_set = malloc(sizeof(int) * g->lines_count);
 1853 
 1854   // some accounting variables
 1855   int vnb = 0, vcount = 0;
 1856   int hnb = 0, hcount = 0;
 1857 
 1858   // just to be on the safe side
 1859   if(g->lines == NULL) goto error;
 1860 
 1861   // generate index list for the vertical lines
 1862   for(int n = 0; n < g->lines_count; n++)
 1863   {
 1864     // is this a selected vertical line?
 1865     if((g->lines[n].type & ASHIFT_LINE_MASK) != ASHIFT_LINE_VERTICAL_SELECTED)
 1866       continue;
 1867 
 1868     lines_set[vnb] = n;
 1869     inout_set[vnb] = 0;
 1870     vnb++;
 1871   }
 1872 
 1873   // it only makes sense to call ransac if we have more than two lines
 1874   if(vnb > 2)
 1875     ransac(g->lines, lines_set, inout_set, vnb, g->vertical_weight,
 1876            xmin, xmax, ymin, ymax);
 1877 
 1878   // adjust line selected flag according to the ransac results
 1879   for(int n = 0; n < vnb; n++)
 1880   {
 1881     const int m = lines_set[n];
 1882     if(inout_set[n] == 1)
 1883     {
 1884       g->lines[m].type |= ASHIFT_LINE_SELECTED;
 1885       vcount++;
 1886     }
 1887     else
 1888       g->lines[m].type &= ~ASHIFT_LINE_SELECTED;
 1889   }
 1890   // update number of vertical lines
 1891   g->vertical_count = vcount;
 1892   g->lines_version++;
 1893 
 1894   // now generate index list for the horizontal lines
 1895   for(int n = 0; n < g->lines_count; n++)
 1896   {
 1897     // is this a selected horizontal line?
 1898     if((g->lines[n].type & ASHIFT_LINE_MASK) != ASHIFT_LINE_HORIZONTAL_SELECTED)
 1899       continue;
 1900 
 1901     lines_set[hnb] = n;
 1902     inout_set[hnb] = 0;
 1903     hnb++;
 1904   }
 1905 
 1906   // it only makes sense to call ransac if we have more than two lines
 1907   if(hnb > 2)
 1908     ransac(g->lines, lines_set, inout_set, hnb, g->horizontal_weight,
 1909            xmin, xmax, ymin, ymax);
 1910 
 1911   // adjust line selected flag according to the ransac results
 1912   for(int n = 0; n < hnb; n++)
 1913   {
 1914     const int m = lines_set[n];
 1915     if(inout_set[n] == 1)
 1916     {
 1917       g->lines[m].type |= ASHIFT_LINE_SELECTED;
 1918       hcount++;
 1919     }
 1920     else
 1921       g->lines[m].type &= ~ASHIFT_LINE_SELECTED;
 1922   }
 1923   // update number of horizontal lines
 1924   g->horizontal_count = hcount;
 1925   g->lines_version++;
 1926 
 1927   free(inout_set);
 1928   free(lines_set);
 1929 
 1930   return TRUE;
 1931 
 1932 error:
 1933   free(inout_set);
 1934   free(lines_set);
 1935   return FALSE;
 1936 }
 1937 
 1938 // utility function to map a variable in [min; max] to [-INF; + INF]
 1939 static inline double logit(double x, double min, double max)
 1940 {
 1941   const double eps = 1.0e-6;
 1942   // make sure p does not touch the borders of its definition area,
 1943   // not critical for data accuracy as logit() is only used on initial fit parameters
 1944   double p = CLAMP((x - min) / (max - min), eps, 1.0 - eps);
 1945 
 1946   return (2.0 * atanh(2.0 * p - 1.0));
 1947 }
 1948 
 1949 // inverted function to logit()
 1950 static inline double ilogit(double L, double min, double max)
 1951 {
 1952   double p = 0.5 * (1.0 + tanh(0.5 * L));
 1953 
 1954   return (p * (max - min) + min);
 1955 }
 1956 
 1957 // helper function for simplex() return quality parameter for the given model
 1958 // strategy:
 1959 //    * generate homography matrix out of fixed parameters and fitting parameters
 1960 //    * apply homography to all end points of affected lines
 1961 //    * generate new line out of transformed end points
 1962 //    * calculate scalar product s of line with perpendicular axis
 1963 //    * sum over weighted s^2 values
 1964 static double model_fitness(double *params, void *data)
 1965 {
 1966   dt_iop_ashift_fit_params_t *fit = (dt_iop_ashift_fit_params_t *)data;
 1967 
 1968   // just for convenience: get shorter names
 1969   dt_iop_ashift_line_t *lines = fit->lines;
 1970   const int lines_count = fit->lines_count;
 1971   const int width = fit->width;
 1972   const int height = fit->height;
 1973   const float f_length_kb = fit->f_length_kb;
 1974   const float orthocorr = fit->orthocorr;
 1975   const float aspect = fit->aspect;
 1976 
 1977   float rotation = fit->rotation;
 1978   float lensshift_v = fit->lensshift_v;
 1979   float lensshift_h = fit->lensshift_h;
 1980   float shear = fit->shear;
 1981   float rotation_range = fit->rotation_range;
 1982   float lensshift_v_range = fit->lensshift_v_range;
 1983   float lensshift_h_range = fit->lensshift_h_range;
 1984   float shear_range = fit->shear_range;
 1985 
 1986   int pcount = 0;
 1987 
 1988   // fill in fit parameters from params[]. Attention: order matters!!!
 1989   if(isnan(rotation))
 1990   {
 1991     rotation = ilogit(params[pcount], -rotation_range, rotation_range);
 1992     pcount++;
 1993   }
 1994 
 1995   if(isnan(lensshift_v))
 1996   {
 1997     lensshift_v = ilogit(params[pcount], -lensshift_v_range, lensshift_v_range);
 1998     pcount++;
 1999   }
 2000 
 2001   if(isnan(lensshift_h))
 2002   {
 2003     lensshift_h = ilogit(params[pcount], -lensshift_h_range, lensshift_h_range);
 2004     pcount++;
 2005   }
 2006 
 2007   if(isnan(shear))
 2008   {
 2009     shear = ilogit(params[pcount], -shear_range, shear_range);
 2010     pcount++;
 2011   }
 2012 
 2013   assert(pcount == fit->params_count);
 2014 
 2015   // the possible reference axes
 2016   const float Av[3] = { 1.0f, 0.0f, 0.0f };
 2017   const float Ah[3] = { 0.0f, 1.0f, 0.0f };
 2018 
 2019   // generate homograph out of the parameters
 2020   float homograph[3][3];
 2021   homography((float *)homograph, rotation, lensshift_v, lensshift_h, shear, f_length_kb,
 2022              orthocorr, aspect, width, height, ASHIFT_HOMOGRAPH_FORWARD);
 2023 
 2024   // accounting variables
 2025   double sumsq_v = 0.0;
 2026   double sumsq_h = 0.0;
 2027   double weight_v = 0.0;
 2028   double weight_h = 0.0;
 2029   int count_v = 0;
 2030   int count_h = 0;
 2031   int count = 0;
 2032 
 2033   // iterate over all lines
 2034   for(int n = 0; n < lines_count; n++)
 2035   {
 2036     // check if this is a line which we must skip
 2037     if((lines[n].type & fit->linemask) != fit->linetype)
 2038       continue;
 2039 
 2040     // the direction of this line (vertical?)
 2041     const int isvertical = lines[n].type & ASHIFT_LINE_DIRVERT;
 2042 
 2043     // select the perpendicular reference axis
 2044     const float *A = isvertical ? Ah : Av;
 2045 
 2046     // apply homographic transformation to the end points
 2047     float P1[3], P2[3];
 2048     mat3mulv(P1, (float *)homograph, lines[n].p1);
 2049     mat3mulv(P2, (float *)homograph, lines[n].p2);
 2050 
 2051     // get line connecting the two points
 2052     float L[3];
 2053     vec3prodn(L, P1, P2);
 2054 
 2055     // normalize L so that x^2 + y^2 = 1; makes sure that
 2056     // y^2 = 1 / (1 + m^2) and x^2 = m^2 / (1 + m^2) with m defining the slope of the line
 2057     vec3lnorm(L, L);
 2058 
 2059     // get scalar product of line L with orthogonal axis A -> gives 0 if line is perpendicular
 2060     float s = vec3scalar(L, A);
 2061 
 2062     // sum up weighted s^2 for both directions individually
 2063     sumsq_v += isvertical ? s * s * lines[n].weight : 0.0;
 2064     weight_v  += isvertical ? lines[n].weight : 0.0;
 2065     count_v += isvertical ? 1 : 0;
 2066     sumsq_h += !isvertical ? s * s * lines[n].weight : 0.0;
 2067     weight_h  += !isvertical ? lines[n].weight : 0.0;
 2068     count_h += !isvertical ? 1 : 0;
 2069     count++;
 2070   }
 2071 
 2072   const double v = weight_v > 0.0f && count > 0 ? sumsq_v / weight_v * (float)count_v / count : 0.0;
 2073   const double h = weight_h > 0.0f && count > 0 ? sumsq_h / weight_h * (float)count_h / count : 0.0;
 2074 
 2075   double sum = sqrt(1.0 - (1.0 - v) * (1.0 - h)) * 1.0e6;
 2076   //double sum = sqrt(v + h) * 1.0e6;
 2077 
 2078 #ifdef ASHIFT_DEBUG
 2079   printf("fitness with rotation %f, lensshift_v %f, lensshift_h %f, shear %f -> lines %d, quality %10f\n",
 2080          rotation, lensshift_v, lensshift_h, shear, count, sum);
 2081 #endif
 2082 
 2083   return sum;
 2084 }
 2085 
 2086 // setup all data structures for fitting and call NM simplex
 2087 static dt_iop_ashift_nmsresult_t nmsfit(dt_iop_module_t *module, dt_iop_ashift_params_t *p, dt_iop_ashift_fitaxis_t dir)
 2088 {
 2089   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)module->gui_data;
 2090 
 2091   if(!g->lines) return NMS_NOT_ENOUGH_LINES;
 2092   if(dir == ASHIFT_FIT_NONE) return NMS_SUCCESS;
 2093 
 2094   double params[4];
 2095   int pcount = 0;
 2096   int enough_lines = TRUE;
 2097 
 2098   // initialize fit parameters
 2099   dt_iop_ashift_fit_params_t fit;
 2100   fit.lines = g->lines;
 2101   fit.lines_count = g->lines_count;
 2102   fit.width = g->lines_in_width;
 2103   fit.height = g->lines_in_height;
 2104   fit.f_length_kb = (p->mode == ASHIFT_MODE_GENERIC) ? DEFAULT_F_LENGTH : p->f_length * p->crop_factor;
 2105   fit.orthocorr = (p->mode == ASHIFT_MODE_GENERIC) ? 0.0f : p->orthocorr;
 2106   fit.aspect = (p->mode == ASHIFT_MODE_GENERIC) ? 1.0f : p->aspect;
 2107   fit.rotation = p->rotation;
 2108   fit.lensshift_v = p->lensshift_v;
 2109   fit.lensshift_h = p->lensshift_h;
 2110   fit.shear = p->shear;
 2111   fit.rotation_range = g->rotation_range;
 2112   fit.lensshift_v_range = g->lensshift_v_range;
 2113   fit.lensshift_h_range = g->lensshift_h_range;
 2114   fit.shear_range = g->shear_range;
 2115   fit.linetype = ASHIFT_LINE_RELEVANT | ASHIFT_LINE_SELECTED;
 2116   fit.linemask = ASHIFT_LINE_MASK;
 2117   fit.params_count = 0;
 2118   fit.weight = 0.0f;
 2119 
 2120   // if the image is flipped and if we do not want to fit both lens shift
 2121   // directions or none at all, then we need to change direction
 2122   dt_iop_ashift_fitaxis_t mdir = dir;
 2123   if((mdir & ASHIFT_FIT_LENS_BOTH) != ASHIFT_FIT_LENS_BOTH &&
 2124      (mdir & ASHIFT_FIT_LENS_BOTH) != 0)
 2125   {
 2126     // flip all directions
 2127     mdir ^= g->isflipped ? ASHIFT_FIT_FLIP : 0;
 2128     // special case that needs to be corrected
 2129     mdir |= (mdir & ASHIFT_FIT_LINES_BOTH) == 0 ? ASHIFT_FIT_LINES_BOTH : 0;
 2130   }
 2131 
 2132 
 2133   // prepare fit structure and starting parameters for simplex fit.
 2134   // note: the sequence of parameters in params[] needs to match the
 2135   // respective order in dt_iop_ashift_fit_params_t. Parameters which are
 2136   // to be fittet are marked with NAN in the fit structure. Non-NAN
 2137   // parameters are assumed to be constant.
 2138   if(mdir & ASHIFT_FIT_ROTATION)
 2139   {
 2140     // we fit rotation
 2141     fit.params_count++;
 2142     params[pcount] = logit(fit.rotation, -fit.rotation_range, fit.rotation_range);
 2143     pcount++;
 2144     fit.rotation = NAN;
 2145   }
 2146 
 2147   if(mdir & ASHIFT_FIT_LENS_VERT)
 2148   {
 2149     // we fit vertical lens shift
 2150     fit.params_count++;
 2151     params[pcount] = logit(fit.lensshift_v, -fit.lensshift_v_range, fit.lensshift_v_range);
 2152     pcount++;
 2153     fit.lensshift_v = NAN;
 2154   }
 2155 
 2156   if(mdir & ASHIFT_FIT_LENS_HOR)
 2157   {
 2158     // we fit horizontal lens shift
 2159     fit.params_count++;
 2160     params[pcount] = logit(fit.lensshift_h, -fit.lensshift_h_range, fit.lensshift_h_range);
 2161     pcount++;
 2162     fit.lensshift_h = NAN;
 2163   }
 2164 
 2165   if(mdir & ASHIFT_FIT_SHEAR)
 2166   {
 2167     // we fit the shear parameter
 2168     fit.params_count++;
 2169     params[pcount] = logit(fit.shear, -fit.shear_range, fit.shear_range);
 2170     pcount++;
 2171     fit.shear = NAN;
 2172   }
 2173 
 2174   if(mdir & ASHIFT_FIT_LINES_VERT)
 2175   {
 2176     // we use vertical lines for fitting
 2177     fit.linetype |= ASHIFT_LINE_DIRVERT;
 2178     fit.weight += g->vertical_weight;
 2179     enough_lines = enough_lines && (g->vertical_count >= MINIMUM_FITLINES);
 2180   }
 2181 
 2182   if(mdir & ASHIFT_FIT_LINES_HOR)
 2183   {
 2184     // we use horizontal lines for fitting
 2185     fit.linetype |= 0;
 2186     fit.weight += g->horizontal_weight;
 2187     enough_lines = enough_lines && (g->horizontal_count >= MINIMUM_FITLINES);
 2188   }
 2189 
 2190   // this needs to come after ASHIFT_FIT_LINES_VERT and ASHIFT_FIT_LINES_HOR
 2191   if((mdir & ASHIFT_FIT_LINES_BOTH) == ASHIFT_FIT_LINES_BOTH)
 2192   {
 2193     // if we use fitting in both directions we need to
 2194     // adjust fit.linetype and fit.linemask to match all selected lines
 2195     fit.linetype = ASHIFT_LINE_RELEVANT | ASHIFT_LINE_SELECTED;
 2196     fit.linemask = ASHIFT_LINE_RELEVANT | ASHIFT_LINE_SELECTED;
 2197   }
 2198 
 2199   // error case: we do not run simplex if there are not enough lines
 2200   if(!enough_lines)
 2201   {
 2202 #ifdef ASHIFT_DEBUG
 2203     printf("optimization not possible: insufficient number of lines\n");
 2204 #endif
 2205     return NMS_NOT_ENOUGH_LINES;
 2206   }
 2207 
 2208   // start the simplex fit
 2209   int iter = simplex(model_fitness, params, fit.params_count, NMS_EPSILON, NMS_SCALE, NMS_ITERATIONS, NULL, (void*)&fit);
 2210 
 2211   // error case: the fit did not converge
 2212   if(iter >= NMS_ITERATIONS)
 2213   {
 2214 #ifdef ASHIFT_DEBUG
 2215     printf("optimization not successful: maximum number of iterations reached (%d)\n", iter);
 2216 #endif
 2217     return NMS_DID_NOT_CONVERGE;
 2218   }
 2219 
 2220   // fit was successful: now consolidate the results (order matters!!!)
 2221   pcount = 0;
 2222   fit.rotation = isnan(fit.rotation) ? ilogit(params[pcount++], -fit.rotation_range, fit.rotation_range) : fit.rotation;
 2223   fit.lensshift_v = isnan(fit.lensshift_v) ? ilogit(params[pcount++], -fit.lensshift_v_range, fit.lensshift_v_range) : fit.lensshift_v;
 2224   fit.lensshift_h = isnan(fit.lensshift_h) ? ilogit(params[pcount++], -fit.lensshift_h_range, fit.lensshift_h_range) : fit.lensshift_h;
 2225   fit.shear = isnan(fit.shear) ? ilogit(params[pcount++], -fit.shear_range, fit.shear_range) : fit.shear;
 2226 #ifdef ASHIFT_DEBUG
 2227   printf("params after optimization (%d iterations): rotation %f, lensshift_v %f, lensshift_h %f, shear %f\n",
 2228          iter, fit.rotation, fit.lensshift_v, fit.lensshift_h, fit.shear);
 2229 #endif
 2230 
 2231   // sanity check: in case of extreme values the image gets distorted so strongly that it spans an insanely huge area. we check that
 2232   // case and assume values that increase the image area by more than a factor of 4 as being insane.
 2233   float homograph[3][3];
 2234   homography((float *)homograph, fit.rotation, fit.lensshift_v, fit.lensshift_h, fit.shear, fit.f_length_kb,
 2235              fit.orthocorr, fit.aspect, fit.width, fit.height, ASHIFT_HOMOGRAPH_FORWARD);
 2236 
 2237   // visit all four corners and find maximum span
 2238   float xm = FLT_MAX, xM = -FLT_MAX, ym = FLT_MAX, yM = -FLT_MAX;
 2239   for(int y = 0; y < fit.height; y += fit.height - 1)
 2240     for(int x = 0; x < fit.width; x += fit.width - 1)
 2241     {
 2242       float pi[3], po[3];
 2243       pi[0] = x;
 2244       pi[1] = y;
 2245       pi[2] = 1.0f;
 2246       mat3mulv(po, (float *)homograph, pi);
 2247       po[0] /= po[2];
 2248       po[1] /= po[2];
 2249       xm = fmin(xm, po[0]);
 2250       ym = fmin(ym, po[1]);
 2251       xM = fmax(xM, po[0]);
 2252       yM = fmax(yM, po[1]);
 2253     }
 2254 
 2255   if((xM - xm) * (yM - ym) > 4.0f * fit.width * fit.height)
 2256   {
 2257 #ifdef ASHIFT_DEBUG
 2258     printf("optimization not successful: degenerate case with area growth factor (%f) exceeding limits\n",
 2259            (xM - xm) * (yM - ym) / (fit.width * fit.height));
 2260 #endif
 2261     return NMS_INSANE;
 2262   }
 2263 
 2264   // now write the results into structure p
 2265   p->rotation = fit.rotation;
 2266   p->lensshift_v = fit.lensshift_v;
 2267   p->lensshift_h = fit.lensshift_h;
 2268   p->shear = fit.shear;
 2269   return NMS_SUCCESS;
 2270 }
 2271 
 2272 #ifdef ASHIFT_DEBUG
 2273 // only used in development phase. call model_fitness() with current parameters and
 2274 // print some useful information
 2275 static void model_probe(dt_iop_module_t *module, dt_iop_ashift_params_t *p, dt_iop_ashift_fitaxis_t dir)
 2276 {
 2277   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)module->gui_data;
 2278 
 2279   if(!g->lines) return;
 2280   if(dir == ASHIFT_FIT_NONE) return;
 2281 
 2282   double params[4];
 2283   int enough_lines = TRUE;
 2284 
 2285   // initialize fit parameters
 2286   dt_iop_ashift_fit_params_t fit;
 2287   fit.lines = g->lines;
 2288   fit.lines_count = g->lines_count;
 2289   fit.width = g->lines_in_width;
 2290   fit.height = g->lines_in_height;
 2291   fit.f_length_kb = (p->mode == ASHIFT_MODE_GENERIC) ? DEFAULT_F_LENGTH : p->f_length * p->crop_factor;
 2292   fit.orthocorr = (p->mode == ASHIFT_MODE_GENERIC) ? 0.0f : p->orthocorr;
 2293   fit.aspect = (p->mode == ASHIFT_MODE_GENERIC) ? 1.0f : p->aspect;
 2294   fit.rotation = p->rotation;
 2295   fit.lensshift_v = p->lensshift_v;
 2296   fit.lensshift_h = p->lensshift_h;
 2297   fit.shear = p->shear;
 2298   fit.linetype = ASHIFT_LINE_RELEVANT | ASHIFT_LINE_SELECTED;
 2299   fit.linemask = ASHIFT_LINE_MASK;
 2300   fit.params_count = 0;
 2301   fit.weight = 0.0f;
 2302 
 2303   // if the image is flipped and if we do not want to fit both lens shift
 2304   // directions or none at all, then we need to change direction
 2305   dt_iop_ashift_fitaxis_t mdir = dir;
 2306   if((mdir & ASHIFT_FIT_LENS_BOTH) != ASHIFT_FIT_LENS_BOTH &&
 2307      (mdir & ASHIFT_FIT_LENS_BOTH) != 0)
 2308   {
 2309     // flip all directions
 2310     mdir ^= g->isflipped ? ASHIFT_FIT_FLIP : 0;
 2311     // special case that needs to be corrected
 2312     mdir |= (mdir & ASHIFT_FIT_LINES_BOTH) == 0 ? ASHIFT_FIT_LINES_BOTH : 0;
 2313   }
 2314 
 2315   if(mdir & ASHIFT_FIT_LINES_VERT)
 2316   {
 2317     // we use vertical lines for fitting
 2318     fit.linetype |= ASHIFT_LINE_DIRVERT;
 2319     fit.weight += g->vertical_weight;
 2320     enough_lines = enough_lines && (g->vertical_count >= MINIMUM_FITLINES);
 2321   }
 2322 
 2323   if(mdir & ASHIFT_FIT_LINES_HOR)
 2324   {
 2325     // we use horizontal lines for fitting
 2326     fit.linetype |= 0;
 2327     fit.weight += g->horizontal_weight;
 2328     enough_lines = enough_lines && (g->horizontal_count >= MINIMUM_FITLINES);
 2329   }
 2330 
 2331   // this needs to come after ASHIFT_FIT_LINES_VERT and ASHIFT_FIT_LINES_HOR
 2332   if((mdir & ASHIFT_FIT_LINES_BOTH) == ASHIFT_FIT_LINES_BOTH)
 2333   {
 2334     // if we use fitting in both directions we need to
 2335     // adjust fit.linetype and fit.linemask to match all selected lines
 2336     fit.linetype = ASHIFT_LINE_RELEVANT | ASHIFT_LINE_SELECTED;
 2337     fit.linemask = ASHIFT_LINE_RELEVANT | ASHIFT_LINE_SELECTED;
 2338   }
 2339 
 2340   double quality = model_fitness(params, (void *)&fit);
 2341 
 2342   printf("model fitness: %.8f (rotation %f, lensshift_v %f, lensshift_h %f, shear %f)\n",
 2343          quality, p->rotation, p->lensshift_v, p->lensshift_h, p->shear);
 2344 }
 2345 #endif
 2346 
 2347 // function to keep crop fitting parameters within constraints
 2348 static void crop_constraint(double *params, int pcount)
 2349 {
 2350   if(pcount > 0) params[0] = fabs(params[0]);
 2351   if(pcount > 1) params[1] = fabs(params[1]);
 2352   if(pcount > 2) params[2] = fabs(params[2]);
 2353 
 2354   if(pcount > 0 && params[0] > 1.0) params[0] = 1.0 - params[0];
 2355   if(pcount > 1 && params[1] > 1.0) params[1] = 1.0 - params[1];
 2356   if(pcount > 2 && params[2] > 0.5*M_PI) params[2] = 0.5*M_PI - params[2];
 2357 }
 2358 
 2359 // helper function for getting the best fitting crop area;
 2360 // returns the negative area of the largest rectangle that fits within the
 2361 // defined image with a given rectangle's center and its aspect angle;
 2362 // the trick: the rectangle center coordinates are given in the input
 2363 // image coordinates so we know for sure that it also lies within the image after
 2364 // conversion to the output coordinates
 2365 static double crop_fitness(double *params, void *data)
 2366 {
 2367   dt_iop_ashift_cropfit_params_t *cropfit = (dt_iop_ashift_cropfit_params_t *)data;
 2368 
 2369   const float wd = cropfit->width;
 2370   const float ht = cropfit->height;
 2371 
 2372   // get variable and constant parameters, respectively
 2373   const float x = isnan(cropfit->x) ? params[0] : cropfit->x;
 2374   const float y = isnan(cropfit->y) ? params[1] : cropfit->y;
 2375   const float alpha = isnan(cropfit->alpha) ? params[2] : cropfit->alpha;
 2376 
 2377   // the center of the rectangle in input image coordinates
 2378   const float Pc[3] = { x * wd, y * ht, 1.0f };
 2379 
 2380   // convert to the output image coordinates and normalize
 2381   float P[3];
 2382   mat3mulv(P, (float *)cropfit->homograph, Pc);
 2383   P[0] /= P[2];
 2384   P[1] /= P[2];
 2385   P[2] = 1.0f;
 2386 
 2387   // two auxiliary points (some arbitrary distance away from P) to construct the diagonals
 2388   const float Pa[2][3] = { { P[0] + 10.0f * cosf(alpha), P[1] + 10.0f * sinf(alpha), 1.0f },
 2389                            { P[0] + 10.0f * cosf(alpha), P[1] - 10.0f * sinf(alpha), 1.0f } };
 2390 
 2391   // the two diagonals: D = P x Pa
 2392   float D[2][3];
 2393   vec3prodn(D[0], P, Pa[0]);
 2394   vec3prodn(D[1], P, Pa[1]);
 2395 
 2396   // find all intersection points of all four edges with both diagonals (I = E x D);
 2397   // the shortest distance d2min of the intersection point I to the crop area center P determines
 2398   // the size of the crop area that still fits into the image (for the given center and aspect angle)
 2399   float d2min = FLT_MAX;
 2400   for(int k = 0; k < 4; k++)
 2401     for(int l = 0; l < 2; l++)
 2402     {
 2403       // the intersection point
 2404       float I[3];
 2405       vec3prodn(I, cropfit->edges[k], D[l]);
 2406 
 2407       // special case: I is all null -> E and D are identical -> P lies on E -> d2min = 0
 2408       if(vec3isnull(I))
 2409       {
 2410         d2min = 0.0f;
 2411         break;
 2412       }
 2413 
 2414       // special case: I[2] is 0.0f -> E and D are parallel and intersect at infinity -> no relevant point
 2415       if(I[2] == 0.0f)
 2416         continue;
 2417 
 2418       // the default case -> normalize I
 2419       I[0] /= I[2];
 2420       I[1] /= I[2];
 2421 
 2422       // calculate distance from I to P
 2423       const float d2 = SQR(P[0] - I[0]) + SQR(P[1] - I[1]);
 2424 
 2425       // the minimum distance over all intersection points
 2426       d2min = MIN(d2min, d2);
 2427     }
 2428 
 2429   // calculate the area of the rectangle
 2430   const float A = 2.0f * d2min * sinf(2.0f * alpha);
 2431 
 2432 #ifdef ASHIFT_DEBUG
 2433   printf("crop fitness with x %f, y %f, angle %f -> distance %f, area %f\n",
 2434          x, y, alpha, d2min, A);
 2435 #endif
 2436   // and return -A to allow Nelder-Mead simplex to search for the minimum
 2437   return -A;
 2438 }
 2439 
 2440 // strategy: for a given center of the crop area and a specific aspect angle
 2441 // we calculate the largest crop area that still lies within the output image;
 2442 // now we allow a Nelder-Mead simplex to search for the center coordinates
 2443 // (and optionally the aspect angle) that delivers the largest overall crop area.
 2444 static void do_crop(dt_iop_module_t *module, dt_iop_ashift_params_t *p)
 2445 {
 2446   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)module->gui_data;
 2447 
 2448   // skip if fitting is still running
 2449   if(g->fitting) return;
 2450 
 2451   // reset fit margins if auto-cropping is off
 2452   if(p->cropmode == ASHIFT_CROP_OFF)
 2453   {
 2454     clear_shadow_crop_box(g);
 2455     commit_crop_box(p,g);
 2456     return;
 2457   }
 2458 
 2459   g->fitting = 1;
 2460 
 2461   double params[3];
 2462   int pcount;
 2463 
 2464   // get parameters for the homograph
 2465   const float f_length_kb = (p->mode == ASHIFT_MODE_GENERIC) ? DEFAULT_F_LENGTH : p->f_length * p->crop_factor;
 2466   const float orthocorr = (p->mode == ASHIFT_MODE_GENERIC) ? 0.0f : p->orthocorr;
 2467   const float aspect = (p->mode == ASHIFT_MODE_GENERIC) ? 1.0f : p->aspect;
 2468   const float rotation = p->rotation;
 2469   const float lensshift_v = p->lensshift_v;
 2470   const float lensshift_h = p->lensshift_h;
 2471   const float shear = p->shear;
 2472 
 2473   // prepare structure of constant parameters
 2474   dt_iop_ashift_cropfit_params_t cropfit;
 2475   cropfit.width = g->buf_width;
 2476   cropfit.height = g->buf_height;
 2477   homography((float *)cropfit.homograph, rotation, lensshift_v, lensshift_h, shear, f_length_kb,
 2478              orthocorr, aspect, cropfit.width, cropfit.height, ASHIFT_HOMOGRAPH_FORWARD);
 2479 
 2480   const float wd = cropfit.width;
 2481   const float ht = cropfit.height;
 2482 
 2483   // the four vertices of the image in input image coordinates
 2484   const float Vc[4][3] = { { 0.0f, 0.0f, 1.0f },
 2485                            { 0.0f,   ht, 1.0f },
 2486                            {   wd,   ht, 1.0f },
 2487                            {   wd, 0.0f, 1.0f } };
 2488 
 2489   // convert the vertices to output image coordinates
 2490   float V[4][3];
 2491   for(int n = 0; n < 4; n++)
 2492     mat3mulv(V[n], (float *)cropfit.homograph, Vc[n]);
 2493 
 2494   // get width and height of output image for later use
 2495   float xmin = FLT_MAX, ymin = FLT_MAX, xmax = FLT_MIN, ymax = FLT_MIN;
 2496   for(int n = 0; n < 4; n++)
 2497   {
 2498     // normalize V
 2499     V[n][0] /= V[n][2];
 2500     V[n][1] /= V[n][2];
 2501     V[n][2] = 1.0f;
 2502     xmin = MIN(xmin, V[n][0]);
 2503     xmax = MAX(xmax, V[n][0]);
 2504     ymin = MIN(ymin, V[n][1]);
 2505     ymax = MAX(ymax, V[n][1]);
 2506   }
 2507   const float owd = xmax - xmin;
 2508   const float oht = ymax - ymin;
 2509 
 2510   // calculate the lines defining the four edges of the image area: E = V[n] x V[n+1]
 2511   for(int n = 0; n < 4; n++)
 2512     vec3prodn(cropfit.edges[n], V[n], V[(n + 1) % 4]);
 2513 
 2514   // initial fit parameters: crop area is centered and aspect angle is that of the original image
 2515   // number of parameters: fit only crop center coordinates with a fixed aspect ratio, or fit all three variables
 2516   if(p->cropmode == ASHIFT_CROP_LARGEST)
 2517   {
 2518     params[0] = 0.5;
 2519     params[1] = 0.5;
 2520     params[2] = atan2f((float)cropfit.height, (float)cropfit.width);
 2521     cropfit.x = NAN;
 2522     cropfit.y = NAN;
 2523     cropfit.alpha = NAN;
 2524     pcount = 3;
 2525   }
 2526   else //(p->cropmode == ASHIFT_CROP_ASPECT)
 2527   {
 2528     params[0] = 0.5;
 2529     params[1] = 0.5;
 2530     cropfit.x = NAN;
 2531     cropfit.y = NAN;
 2532     cropfit.alpha = atan2f((float)cropfit.height, (float)cropfit.width);
 2533     pcount = 2;
 2534   }
 2535 
 2536   // start the simplex fit
 2537   const int iter = simplex(crop_fitness, params, pcount, NMS_CROP_EPSILON, NMS_CROP_SCALE, NMS_CROP_ITERATIONS,
 2538                            crop_constraint, (void*)&cropfit);
 2539   // in case the fit did not converge -> failed
 2540   if(iter >= NMS_CROP_ITERATIONS) goto failed;
 2541 
 2542   // the fit did converge -> get clipping margins out of params:
 2543   cropfit.x = isnan(cropfit.x) ? params[0] : cropfit.x;
 2544   cropfit.y = isnan(cropfit.y) ? params[1] : cropfit.y;
 2545   cropfit.alpha = isnan(cropfit.alpha) ? params[2] : cropfit.alpha;
 2546 
 2547   // the area of the best fitting rectangle
 2548   const float A = fabs(crop_fitness(params, (void*)&cropfit));
 2549 
 2550   // unlikely to happen but we need to catch this case
 2551   if(A == 0.0f) goto failed;
 2552 
 2553   // we need the half diagonal of that rectangle (this is in output image dimensions);
 2554   // no need to check for division by zero here as this case implies A == 0.0f, caught above
 2555   const float d = sqrtf(A / (2.0f * sinf(2.0f * cropfit.alpha)));
 2556 
 2557   // the rectangle's center in input image (homogeneous) coordinates
 2558   const float Pc[3] = { cropfit.x * wd, cropfit.y * ht, 1.0f };
 2559 
 2560   // convert rectangle center to output image coordinates and normalize
 2561   float P[3];
 2562   mat3mulv(P, (float *)cropfit.homograph, Pc);
 2563   P[0] /= P[2];
 2564   P[1] /= P[2];
 2565 
 2566   // calculate clipping margins relative to output image dimensions
 2567   g->cl = CLAMP((P[0] - d * cosf(cropfit.alpha)) / owd, 0.0f, 1.0f);
 2568   g->cr = CLAMP((P[0] + d * cosf(cropfit.alpha)) / owd, 0.0f, 1.0f);
 2569   g->ct = CLAMP((P[1] - d * sinf(cropfit.alpha)) / oht, 0.0f, 1.0f);
 2570   g->cb = CLAMP((P[1] + d * sinf(cropfit.alpha)) / oht, 0.0f, 1.0f);
 2571 
 2572   // final sanity check
 2573   if(g->cr - g->cl <= 0.0f || g->cb - g->ct <= 0.0f) goto failed;
 2574 
 2575   g->fitting = 0;
 2576 
 2577 #ifdef ASHIFT_DEBUG
 2578   printf("margins after crop fitting: iter %d, x %f, y %f, angle %f, crop area (%f %f %f %f), width %f, height %f\n",
 2579          iter, cropfit.x, cropfit.y, cropfit.alpha, g->cl, g->cr, g->ct, g->cb, wd, ht);
 2580 #endif
 2581   dt_control_queue_redraw_center();
 2582   return;
 2583 
 2584 failed:
 2585   // in case of failure: reset clipping margins, set "automatic cropping" parameter
 2586   // to "off" state, and display warning message
 2587   clear_shadow_crop_box(g);
 2588   commit_crop_box(p,g);
 2589   p->cropmode = ASHIFT_CROP_OFF;
 2590   dt_bauhaus_combobox_set(g->cropmode, p->cropmode);
 2591   g->fitting = 0;
 2592   dt_control_log(_("automatic cropping failed"));
 2593   return;
 2594 }
 2595 
 2596 // manually adjust crop area by shifting its center
 2597 static void crop_adjust(dt_iop_module_t *module, const dt_iop_ashift_params_t *const p,
 2598                         const float newx, const float newy)
 2599 {
 2600   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)module->gui_data;
 2601 
 2602   // skip if fitting is still running
 2603   if(g->fitting) return;
 2604 
 2605   // get parameters for the homograph
 2606   const float f_length_kb = (p->mode == ASHIFT_MODE_GENERIC) ? DEFAULT_F_LENGTH : p->f_length * p->crop_factor;
 2607   const float orthocorr = (p->mode == ASHIFT_MODE_GENERIC) ? 0.0f : p->orthocorr;
 2608   const float aspect = (p->mode == ASHIFT_MODE_GENERIC) ? 1.0f : p->aspect;
 2609   const float rotation = p->rotation;
 2610   const float lensshift_v = p->lensshift_v;
 2611   const float lensshift_h = p->lensshift_h;
 2612   const float shear = p->shear;
 2613   const float wd = g->buf_width;
 2614   const float ht = g->buf_height;
 2615 
 2616   const float alpha = atan2f(ht, wd);
 2617 
 2618   float homograph[3][3];
 2619   homography((float *)homograph, rotation, lensshift_v, lensshift_h, shear, f_length_kb,
 2620              orthocorr, aspect, wd, ht, ASHIFT_HOMOGRAPH_FORWARD);
 2621 
 2622   // the four vertices of the image in input image coordinates
 2623   const float Vc[4][3] = { { 0.0f, 0.0f, 1.0f },
 2624                            { 0.0f,   ht, 1.0f },
 2625                            {   wd,   ht, 1.0f },
 2626                            {   wd, 0.0f, 1.0f } };
 2627 
 2628   // convert the vertices to output image coordinates
 2629   float V[4][3];
 2630   for(int n = 0; n < 4; n++)
 2631     mat3mulv(V[n], (float *)homograph, Vc[n]);
 2632 
 2633   // get width and height of output image
 2634   float xmin = FLT_MAX, ymin = FLT_MAX, xmax = FLT_MIN, ymax = FLT_MIN;
 2635   for(int n = 0; n < 4; n++)
 2636   {
 2637     // normalize V
 2638     V[n][0] /= V[n][2];
 2639     V[n][1] /= V[n][2];
 2640     V[n][2] = 1.0f;
 2641     xmin = MIN(xmin, V[n][0]);
 2642     xmax = MAX(xmax, V[n][0]);
 2643     ymin = MIN(ymin, V[n][1]);
 2644     ymax = MAX(ymax, V[n][1]);
 2645   }
 2646   const float owd = xmax - xmin;
 2647   const float oht = ymax - ymin;
 2648 
 2649   // calculate the lines defining the four edges of the image area: E = V[n] x V[n+1]
 2650   float E[4][3];
 2651   for(int n = 0; n < 4; n++)
 2652     vec3prodn(E[n], V[n], V[(n + 1) % 4]);
 2653 
 2654   // the center of the rectangle in output image coordinates
 2655   const float P[3] = { newx * owd, newy * oht, 1.0f };
 2656 
 2657   // two auxiliary points (some arbitrary distance away from P) to construct the diagonals
 2658   const float Pa[2][3] = { { P[0] + 10.0f * cosf(alpha), P[1] + 10.0f * sinf(alpha), 1.0f },
 2659                            { P[0] + 10.0f * cosf(alpha), P[1] - 10.0f * sinf(alpha), 1.0f } };
 2660 
 2661   // the two diagonals: D = P x Pa
 2662   float D[2][3];
 2663   vec3prodn(D[0], P, Pa[0]);
 2664   vec3prodn(D[1], P, Pa[1]);
 2665 
 2666   // find all intersection points of all four edges with both diagonals (I = E x D);
 2667   // the shortest distance d2min of the intersection point I to the crop area center P determines
 2668   // the size of the crop area that still fits into the image (for the given center and aspect angle)
 2669   float d2min = FLT_MAX;
 2670   for(int k = 0; k < 4; k++)
 2671     for(int l = 0; l < 2; l++)
 2672     {
 2673       // the intersection point
 2674       float I[3];
 2675       vec3prodn(I, E[k], D[l]);
 2676 
 2677       // special case: I is all null -> E and D are identical -> P lies on E -> d2min = 0
 2678       if(vec3isnull(I))
 2679       {
 2680         d2min = 0.0f;
 2681         break;
 2682       }
 2683 
 2684       // special case: I[2] is 0.0f -> E and D are parallel and intersect at infinity -> no relevant point
 2685       if(I[2] == 0.0f)
 2686         continue;
 2687 
 2688       // the default case -> normalize I
 2689       I[0] /= I[2];
 2690       I[1] /= I[2];
 2691 
 2692       // calculate distance from I to P
 2693       const float d2 = SQR(P[0] - I[0]) + SQR(P[1] - I[1]);
 2694 
 2695       // the minimum distance over all intersection points
 2696       d2min = MIN(d2min, d2);
 2697     }
 2698 
 2699   const float d = sqrtf(d2min);
 2700 
 2701   // do not allow crop area to drop below 1% of input image area
 2702   const float A = 2.0f * d * d * sinf(2.0f * alpha);
 2703   if(A < 0.01f * wd * ht) return;
 2704 
 2705   // calculate clipping margins relative to output image dimensions
 2706   g->cl = CLAMP((P[0] - d * cosf(alpha)) / owd, 0.0f, 1.0f);
 2707   g->cr = CLAMP((P[0] + d * cosf(alpha)) / owd, 0.0f, 1.0f);
 2708   g->ct = CLAMP((P[1] - d * sinf(alpha)) / oht, 0.0f, 1.0f);
 2709   g->cb = CLAMP((P[1] + d * sinf(alpha)) / oht, 0.0f, 1.0f);
 2710 
 2711 #ifdef ASHIFT_DEBUG
 2712   printf("margins after crop adjustment: x %f, y %f, angle %f, crop area (%f %f %f %f), width %f, height %f\n",
 2713          0.5f * (g->cl + g->cr), 0.5f * (g->ct + g->cb), alpha, g->cl, g->cr, g->ct, g->cb, wd, ht);
 2714 #endif
 2715   return;
 2716 }
 2717 
 2718 // helper function to start analysis for structural data and report about errors
 2719 static int do_get_structure(dt_iop_module_t *module, dt_iop_ashift_params_t *p,
 2720                             dt_iop_ashift_enhance_t enhance)
 2721 {
 2722   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)module->gui_data;
 2723 
 2724   if(g->fitting) return FALSE;
 2725 
 2726   g->fitting = 1;
 2727 
 2728   dt_iop_gui_enter_critical_section(module);
 2729   float *b = g->buf;
 2730   dt_iop_gui_leave_critical_section(module);
 2731 
 2732   if(b == NULL)
 2733   {
 2734     dt_control_log(_("data pending - please repeat"));
 2735     goto error;
 2736   }
 2737 
 2738   if(!get_structure(module, enhance))
 2739   {
 2740     dt_control_log(_("could not detect structural data in image"));
 2741 #ifdef ASHIFT_DEBUG
 2742     // find out more
 2743     printf("do_get_structure: buf %p, buf_hash %lu, buf_width %d, buf_height %d, lines %p, lines_count %d\n",
 2744            g->buf, g->buf_hash, g->buf_width, g->buf_height, g->lines, g->lines_count);
 2745 #endif
 2746     goto error;
 2747   }
 2748 
 2749   if(!remove_outliers(module))
 2750   {
 2751     dt_control_log(_("could not run outlier removal"));
 2752 #ifdef ASHIFT_DEBUG
 2753     // find out more
 2754     printf("remove_outliers: buf %p, buf_hash %lu, buf_width %d, buf_height %d, lines %p, lines_count %d\n",
 2755            g->buf, g->buf_hash, g->buf_width, g->buf_height, g->lines, g->lines_count);
 2756 #endif
 2757     goto error;
 2758   }
 2759 
 2760   g->fitting = 0;
 2761   return TRUE;
 2762 
 2763 error:
 2764   g->fitting = 0;
 2765   return FALSE;
 2766 }
 2767 
 2768 // helper function to clean structural data
 2769 static int do_clean_structure(dt_iop_module_t *module, dt_iop_ashift_params_t *p)
 2770 {
 2771   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)module->gui_data;
 2772 
 2773   if(g->fitting) return FALSE;
 2774 
 2775   g->fitting = 1;
 2776   g->lines_count = 0;
 2777   g->vertical_count = 0;
 2778   g->horizontal_count = 0;
 2779   free(g->lines);
 2780   g->lines = NULL;
 2781   g->lines_version++;
 2782   g->lines_suppressed = 0;
 2783   g->fitting = 0;
 2784   return TRUE;
 2785 }
 2786 
 2787 // helper function to start parameter fit and report about errors
 2788 static int do_fit(dt_iop_module_t *module, dt_iop_ashift_params_t *p, dt_iop_ashift_fitaxis_t dir)
 2789 {
 2790   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)module->gui_data;
 2791 
 2792   if(g->fitting) return FALSE;
 2793 
 2794   // if no structure available get it
 2795   if(g->lines == NULL)
 2796     if(!do_get_structure(module, p, ASHIFT_ENHANCE_NONE)) goto error;
 2797 
 2798   g->fitting = 1;
 2799 
 2800   dt_iop_ashift_nmsresult_t res = nmsfit(module, p, dir);
 2801 
 2802   switch(res)
 2803   {
 2804     case NMS_NOT_ENOUGH_LINES:
 2805       dt_control_log(_("not enough structure for automatic correction"));
 2806       goto error;
 2807       break;
 2808     case NMS_DID_NOT_CONVERGE:
 2809     case NMS_INSANE:
 2810       dt_control_log(_("automatic correction failed, please correct manually"));
 2811       goto error;
 2812       break;
 2813     case NMS_SUCCESS:
 2814     default:
 2815       break;
 2816   }
 2817 
 2818   g->fitting = 0;
 2819 
 2820   // finally apply cropping
 2821   do_crop(module, p);
 2822   return TRUE;
 2823 
 2824 error:
 2825   g->fitting = 0;
 2826   return FALSE;
 2827 }
 2828 
 2829 void process(struct dt_iop_module_t *self, dt_dev_pixelpipe_iop_t *piece, const void *const ivoid,
 2830              void *const ovoid, const dt_iop_roi_t *const roi_in, const dt_iop_roi_t *const roi_out)
 2831 {
 2832   dt_iop_ashift_data_t *data = (dt_iop_ashift_data_t *)piece->data;
 2833   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 2834 
 2835   const int ch = piece->colors;
 2836   const int ch_width = ch * roi_in->width;
 2837 
 2838   // only for preview pipe: collect input buffer data and do some other evaluations
 2839   if(self->dev->gui_attached && g && (piece->pipe->type & DT_DEV_PIXELPIPE_PREVIEW) == DT_DEV_PIXELPIPE_PREVIEW)
 2840   {
 2841     // we want to find out if the final output image is flipped in relation to this iop
 2842     // so we can adjust the gui labels accordingly
 2843     const float pr_d = self->dev->preview_downsampling;
 2844     const int width = roi_in->width;
 2845     const int height = roi_in->height;
 2846     const int x_off = roi_in->x;
 2847     const int y_off = roi_in->y;
 2848     const float scale = roi_in->scale / pr_d;
 2849 
 2850     // origin of image and opposite corner as reference points
 2851     float points[4] = { 0.0f, 0.0f, (float)piece->buf_in.width, (float)piece->buf_in.height };
 2852     float ivec[2] = { points[2] - points[0], points[3] - points[1] };
 2853     float ivecl = sqrtf(ivec[0] * ivec[0] + ivec[1] * ivec[1]);
 2854 
 2855     // where do they go?
 2856     dt_dev_distort_backtransform_plus(self->dev, self->dev->preview_pipe, self->iop_order,
 2857                                       DT_DEV_TRANSFORM_DIR_FORW_EXCL, points, 2);
 2858 
 2859     float ovec[2] = { points[2] - points[0], points[3] - points[1] };
 2860     float ovecl = sqrtf(ovec[0] * ovec[0] + ovec[1] * ovec[1]);
 2861 
 2862     // angle between input vector and output vector
 2863     float alpha = acos(CLAMP((ivec[0] * ovec[0] + ivec[1] * ovec[1]) / (ivecl * ovecl), -1.0f, 1.0f));
 2864 
 2865     // we are interested if |alpha| is in the range of 90° +/- 45° -> we assume the image is flipped
 2866     int isflipped = fabs(fmod(alpha + M_PI, M_PI) - M_PI / 2.0f) < M_PI / 4.0f ? 1 : 0;
 2867 
 2868     // did modules prior to this one in pixelpipe have changed? -> check via hash value
 2869     uint64_t hash = dt_dev_hash_plus(self->dev, self->dev->preview_pipe, self->iop_order, DT_DEV_TRANSFORM_DIR_BACK_EXCL);
 2870 
 2871     dt_iop_gui_enter_critical_section(self);
 2872     g->isflipped = isflipped;
 2873 
 2874     // save a copy of preview input buffer for parameter fitting
 2875     if(g->buf == NULL || (size_t)g->buf_width * g->buf_height < (size_t)width * height)
 2876     {
 2877       // if needed allocate buffer
 2878       free(g->buf); // a no-op if g->buf is NULL
 2879       // only get new buffer if no old buffer available or old buffer does not fit in terms of size
 2880       g->buf = malloc(sizeof(float) * 4 * width * height);
 2881     }
 2882 
 2883     if(g->buf /* && hash != g->buf_hash */)
 2884     {
 2885       // copy data
 2886       dt_iop_image_copy_by_size(g->buf, ivoid, width, height, ch);
 2887 
 2888       g->buf_width = width;
 2889       g->buf_height = height;
 2890       g->buf_x_off = x_off;
 2891       g->buf_y_off = y_off;
 2892       g->buf_scale = scale;
 2893       g->buf_hash = hash;
 2894     }
 2895 
 2896     dt_iop_gui_leave_critical_section(self);
 2897   }
 2898 
 2899   // if module is set to neutral parameters we just copy input->output and are done
 2900   if(isneutral(data))
 2901   {
 2902     dt_iop_image_copy_by_size(ovoid, ivoid, roi_out->width, roi_out->height, ch);
 2903     return;
 2904   }
 2905 
 2906   const struct dt_interpolation *interpolation = dt_interpolation_new(DT_INTERPOLATION_USERPREF_WARP);
 2907 
 2908   float ihomograph[3][3];
 2909   homography((float *)ihomograph, data->rotation, data->lensshift_v, data->lensshift_h, data->shear, data->f_length_kb,
 2910              data->orthocorr, data->aspect, piece->buf_in.width, piece->buf_in.height, ASHIFT_HOMOGRAPH_INVERTED);
 2911 
 2912   // clipping offset
 2913   const float fullwidth = (float)piece->buf_out.width / (data->cr - data->cl);
 2914   const float fullheight = (float)piece->buf_out.height / (data->cb - data->ct);
 2915   const float cx = roi_out->scale * fullwidth * data->cl;
 2916   const float cy = roi_out->scale * fullheight * data->ct;
 2917 
 2918 #ifdef _OPENMP
 2919 #pragma omp parallel for default(none) \
 2920   dt_omp_firstprivate(ch, ch_width, cx, cy, ivoid, ovoid, roi_in, roi_out) \
 2921   shared(ihomograph, interpolation) \
 2922   schedule(static)
 2923 #endif
 2924   // go over all pixels of output image
 2925   for(int j = 0; j < roi_out->height; j++)
 2926   {
 2927     float *const restrict out = ((float *)ovoid) + (size_t)ch * j * roi_out->width;
 2928     for(int i = 0; i < roi_out->width; i++)
 2929     {
 2930       float pin[3], pout[3];
 2931 
 2932       // convert output pixel coordinates to original image coordinates
 2933       pout[0] = roi_out->x + i + cx;
 2934       pout[1] = roi_out->y + j + cy;
 2935       pout[0] /= roi_out->scale;
 2936       pout[1] /= roi_out->scale;
 2937       pout[2] = 1.0f;
 2938 
 2939       // apply homograph
 2940       mat3mulv(pin, (float *)ihomograph, pout);
 2941 
 2942       // convert to input pixel coordinates
 2943       pin[0] /= pin[2];
 2944       pin[1] /= pin[2];
 2945       pin[0] *= roi_in->scale;
 2946       pin[1] *= roi_in->scale;
 2947       pin[0] -= roi_in->x;
 2948       pin[1] -= roi_in->y;
 2949 
 2950       // get output values by interpolation from input image
 2951       dt_interpolation_compute_pixel4c(interpolation, (float *)ivoid, out + ch*i, pin[0], pin[1], roi_in->width,
 2952                                        roi_in->height, ch_width);
 2953     }
 2954   }
 2955 }
 2956 
 2957 #ifdef HAVE_OPENCL
 2958 int process_cl(struct dt_iop_module_t *self, dt_dev_pixelpipe_iop_t *piece, cl_mem dev_in, cl_mem dev_out,
 2959                const dt_iop_roi_t *const roi_in, const dt_iop_roi_t *const roi_out)
 2960 {
 2961   dt_iop_ashift_data_t *d = (dt_iop_ashift_data_t *)piece->data;
 2962   dt_iop_ashift_global_data_t *gd = (dt_iop_ashift_global_data_t *)self->global_data;
 2963   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 2964 
 2965   const int devid = piece->pipe->devid;
 2966   const int iwidth = roi_in->width;
 2967   const int iheight = roi_in->height;
 2968   const int width = roi_out->width;
 2969   const int height = roi_out->height;
 2970 
 2971   cl_int err = -999;
 2972   cl_mem dev_homo = NULL;
 2973 
 2974   // only for preview pipe: collect input buffer data and do some other evaluations
 2975   if(self->dev->gui_attached && g && (piece->pipe->type & DT_DEV_PIXELPIPE_PREVIEW) == DT_DEV_PIXELPIPE_PREVIEW)
 2976   {
 2977     // we want to find out if the final output image is flipped in relation to this iop
 2978     // so we can adjust the gui labels accordingly
 2979     const float pr_d = self->dev->preview_downsampling;
 2980     const int x_off = roi_in->x;
 2981     const int y_off = roi_in->y;
 2982     const float scale = roi_in->scale / pr_d;
 2983 
 2984     // origin of image and opposite corner as reference points
 2985     float points[4] = { 0.0f, 0.0f, (float)piece->buf_in.width, (float)piece->buf_in.height };
 2986     float ivec[2] = { points[2] - points[0], points[3] - points[1] };
 2987     float ivecl = sqrtf(ivec[0] * ivec[0] + ivec[1] * ivec[1]);
 2988 
 2989     // where do they go?
 2990     dt_dev_distort_backtransform_plus(self->dev, self->dev->preview_pipe, self->iop_order,
 2991                                       DT_DEV_TRANSFORM_DIR_FORW_EXCL, points, 2);
 2992 
 2993     float ovec[2] = { points[2] - points[0], points[3] - points[1] };
 2994     float ovecl = sqrtf(ovec[0] * ovec[0] + ovec[1] * ovec[1]);
 2995 
 2996     // angle between input vector and output vector
 2997     float alpha = acos(CLAMP((ivec[0] * ovec[0] + ivec[1] * ovec[1]) / (ivecl * ovecl), -1.0f, 1.0f));
 2998 
 2999     // we are interested if |alpha| is in the range of 90° +/- 45° -> we assume the image is flipped
 3000     int isflipped = fabs(fmod(alpha + M_PI, M_PI) - M_PI / 2.0f) < M_PI / 4.0f ? 1 : 0;
 3001 
 3002     // do modules coming before this one in pixelpipe have changed? -> check via hash value
 3003     uint64_t hash = dt_dev_hash_plus(self->dev, self->dev->preview_pipe, self->iop_order, DT_DEV_TRANSFORM_DIR_BACK_EXCL);
 3004 
 3005     dt_iop_gui_enter_critical_section(self);
 3006     g->isflipped = isflipped;
 3007 
 3008     // save a copy of preview input buffer for parameter fitting
 3009     if(g->buf == NULL || (size_t)g->buf_width * g->buf_height < (size_t)iwidth * iheight)
 3010     {
 3011       // if needed allocate buffer
 3012       free(g->buf); // a no-op if g->buf is NULL
 3013       // only get new buffer if no old buffer or old buffer does not fit in terms of size
 3014       g->buf = malloc(sizeof(float) * 4 * iwidth * iheight);
 3015     }
 3016 
 3017     if(g->buf /* && hash != g->buf_hash */)
 3018     {
 3019       // copy data
 3020       err = dt_opencl_copy_device_to_host(devid, g->buf, dev_in, iwidth, iheight, sizeof(float) * 4);
 3021 
 3022       g->buf_width = iwidth;
 3023       g->buf_height = iheight;
 3024       g->buf_x_off = x_off;
 3025       g->buf_y_off = y_off;
 3026       g->buf_scale = scale;
 3027       g->buf_hash = hash;
 3028     }
 3029     dt_iop_gui_leave_critical_section(self);
 3030     if(err != CL_SUCCESS) goto error;
 3031   }
 3032 
 3033   // if module is set to neutral parameters we just copy input->output and are done
 3034   if(isneutral(d))
 3035   {
 3036     size_t origin[] = { 0, 0, 0 };
 3037     size_t region[] = { width, height, 1 };
 3038     err = dt_opencl_enqueue_copy_image(devid, dev_in, dev_out, origin, origin, region);
 3039     if(err != CL_SUCCESS) goto error;
 3040     return TRUE;
 3041   }
 3042 
 3043   float ihomograph[3][3];
 3044   homography((float *)ihomograph, d->rotation, d->lensshift_v, d->lensshift_h, d->shear, d->f_length_kb,
 3045              d->orthocorr, d->aspect, piece->buf_in.width, piece->buf_in.height, ASHIFT_HOMOGRAPH_INVERTED);
 3046 
 3047   // clipping offset
 3048   const float fullwidth = (float)piece->buf_out.width / (d->cr - d->cl);
 3049   const float fullheight = (float)piece->buf_out.height / (d->cb - d->ct);
 3050   const float cx = roi_out->scale * fullwidth * d->cl;
 3051   const float cy = roi_out->scale * fullheight * d->ct;
 3052 
 3053   dev_homo = dt_opencl_copy_host_to_device_constant(devid, sizeof(float) * 9, ihomograph);
 3054   if(dev_homo == NULL) goto error;
 3055 
 3056   const int iroi[2] = { roi_in->x, roi_in->y };
 3057   const int oroi[2] = { roi_out->x, roi_out->y };
 3058   const float in_scale = roi_in->scale;
 3059   const float out_scale = roi_out->scale;
 3060   const float clip[2] = { cx, cy };
 3061 
 3062   size_t sizes[] = { ROUNDUPWD(width), ROUNDUPHT(height), 1 };
 3063 
 3064   const struct dt_interpolation *interpolation = dt_interpolation_new(DT_INTERPOLATION_USERPREF_WARP);
 3065 
 3066   int ldkernel = -1;
 3067 
 3068   switch(interpolation->id)
 3069   {
 3070     case DT_INTERPOLATION_BILINEAR:
 3071       ldkernel = gd->kernel_ashift_bilinear;
 3072       break;
 3073     case DT_INTERPOLATION_BICUBIC:
 3074       ldkernel = gd->kernel_ashift_bicubic;
 3075       break;
 3076     case DT_INTERPOLATION_LANCZOS2:
 3077       ldkernel = gd->kernel_ashift_lanczos2;
 3078       break;
 3079     case DT_INTERPOLATION_LANCZOS3:
 3080       ldkernel = gd->kernel_ashift_lanczos3;
 3081       break;
 3082     default:
 3083       goto error;
 3084   }
 3085 
 3086   dt_opencl_set_kernel_arg(devid, ldkernel, 0, sizeof(cl_mem), (void *)&dev_in);
 3087   dt_opencl_set_kernel_arg(devid, ldkernel, 1, sizeof(cl_mem), (void *)&dev_out);
 3088   dt_opencl_set_kernel_arg(devid, ldkernel, 2, sizeof(int), (void *)&width);
 3089   dt_opencl_set_kernel_arg(devid, ldkernel, 3, sizeof(int), (void *)&height);
 3090   dt_opencl_set_kernel_arg(devid, ldkernel, 4, sizeof(int), (void *)&iwidth);
 3091   dt_opencl_set_kernel_arg(devid, ldkernel, 5, sizeof(int), (void *)&iheight);
 3092   dt_opencl_set_kernel_arg(devid, ldkernel, 6, 2 * sizeof(int), (void *)iroi);
 3093   dt_opencl_set_kernel_arg(devid, ldkernel, 7, 2 * sizeof(int), (void *)oroi);
 3094   dt_opencl_set_kernel_arg(devid, ldkernel, 8, sizeof(float), (void *)&in_scale);
 3095   dt_opencl_set_kernel_arg(devid, ldkernel, 9, sizeof(float), (void *)&out_scale);
 3096   dt_opencl_set_kernel_arg(devid, ldkernel, 10, 2 * sizeof(float), (void *)clip);
 3097   dt_opencl_set_kernel_arg(devid, ldkernel, 11, sizeof(cl_mem), (void *)&dev_homo);
 3098   err = dt_opencl_enqueue_kernel_2d(devid, ldkernel, sizes);
 3099   if(err != CL_SUCCESS) goto error;
 3100 
 3101   dt_opencl_release_mem_object(dev_homo);
 3102   return TRUE;
 3103 
 3104 error:
 3105   dt_opencl_release_mem_object(dev_homo);
 3106   dt_print(DT_DEBUG_OPENCL, "[opencl_ashift] couldn't enqueue kernel! %d\n", err);
 3107   return FALSE;
 3108 }
 3109 #endif
 3110 
 3111 // gather information about "near"-ness in g->points_idx
 3112 static void get_near(const float *points, dt_iop_ashift_points_idx_t *points_idx,
 3113                      const int lines_count, float pzx, float pzy, float delta)
 3114 {
 3115   const float delta2 = delta * delta;
 3116 
 3117   for(int n = 0; n < lines_count; n++)
 3118   {
 3119     points_idx[n].near = 0;
 3120 
 3121     // skip irrelevant lines
 3122     if(points_idx[n].type == ASHIFT_LINE_IRRELEVANT)
 3123       continue;
 3124 
 3125     // first check if the mouse pointer is outside the bounding box of the line -> skip this line
 3126     if(pzx < points_idx[n].bbx - delta &&
 3127        pzx > points_idx[n].bbX + delta &&
 3128        pzy < points_idx[n].bby - delta &&
 3129        pzy > points_idx[n].bbY + delta)
 3130       continue;
 3131 
 3132     // pointer is inside bounding box
 3133     size_t offset = points_idx[n].offset;
 3134     const int length = points_idx[n].length;
 3135 
 3136     // sanity check (this should not happen)
 3137     if(length < 2) continue;
 3138 
 3139     // check line point by point
 3140     for(int l = 0; l < length; l++, offset++)
 3141     {
 3142       float dx = pzx - points[offset * 2];
 3143       float dy = pzy - points[offset * 2 + 1];
 3144 
 3145       if(dx * dx + dy * dy < delta2)
 3146       {
 3147         points_idx[n].near = 1;
 3148         break;
 3149       }
 3150     }
 3151   }
 3152 }
 3153 
 3154 // mark lines which are inside a rectangular area in isbounding mode
 3155 static void get_bounded_inside(const float *points, dt_iop_ashift_points_idx_t *points_idx,
 3156                                const int points_lines_count, float pzx, float pzy,
 3157                                float pzx2, float pzy2, dt_iop_ashift_bounding_t mode)
 3158 {
 3159   // get bounding box coordinates
 3160   float ax = pzx;
 3161   float ay = pzy;
 3162   float bx = pzx2;
 3163   float by = pzy2;
 3164   if(pzx > pzx2)
 3165   {
 3166     ax = pzx2;
 3167     bx = pzx;
 3168   }
 3169   if(pzy > pzy2)
 3170   {
 3171     ay = pzy2;
 3172     by = pzy;
 3173   }
 3174 
 3175   // we either look for the selected or the deselected lines
 3176   dt_iop_ashift_linetype_t mask = ASHIFT_LINE_SELECTED;
 3177   dt_iop_ashift_linetype_t state = (mode == ASHIFT_BOUNDING_DESELECT) ? ASHIFT_LINE_SELECTED : 0;
 3178 
 3179   for(int n = 0; n < points_lines_count; n++)
 3180   {
 3181     // mark line as "not near" and "not bounded"
 3182     points_idx[n].near = 0;
 3183     points_idx[n].bounded = 0;
 3184 
 3185     // skip irrelevant lines
 3186     if(points_idx[n].type == ASHIFT_LINE_IRRELEVANT)
 3187       continue;
 3188 
 3189     // is the line inside the box ?
 3190     if(points_idx[n].bbx >= ax && points_idx[n].bbx <= bx && points_idx[n].bbX >= ax
 3191        && points_idx[n].bbX <= bx && points_idx[n].bby >= ay && points_idx[n].bby <= by
 3192        && points_idx[n].bbY >= ay && points_idx[n].bbY <= by)
 3193     {
 3194       points_idx[n].bounded = 1;
 3195       // only mark "near"-ness of those lines we are interested in
 3196       points_idx[n].near = ((points_idx[n].type & mask) != state) ? 0 : 1;
 3197     }
 3198   }
 3199 }
 3200 
 3201 // generate hash value for lines taking into account only the end point coordinates
 3202 static uint64_t get_lines_hash(const dt_iop_ashift_line_t *lines, const int lines_count)
 3203 {
 3204   uint64_t hash = 5381;
 3205   for(int n = 0; n < lines_count; n++)
 3206   {
 3207     float v[4] = { lines[n].p1[0], lines[n].p1[1], lines[n].p2[0], lines[n].p2[1] };
 3208     union {
 3209         float f;
 3210         uint32_t u;
 3211     } x;
 3212 
 3213     for(size_t i = 0; i < 4; i++) {
 3214       x.f = v[i];
 3215       hash = ((hash << 5) + hash) ^ x.u;
 3216     }
 3217   }
 3218   return hash;
 3219 }
 3220 
 3221 // update color information in points_idx if lines have changed in terms of type (but not in terms
 3222 // of number or position)
 3223 static int update_colors(struct dt_iop_module_t *self, dt_iop_ashift_points_idx_t *points_idx,
 3224                          int points_lines_count)
 3225 {
 3226   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 3227 
 3228   // is the display flipped relative to the original image?
 3229   const int isflipped = g->isflipped;
 3230 
 3231   // go through all lines
 3232   for(int n = 0; n < points_lines_count; n++)
 3233   {
 3234     const dt_iop_ashift_linetype_t type = points_idx[n].type;
 3235 
 3236     // set line color according to line type/orientation
 3237     // note: if the screen display is flipped versus the original image we need
 3238     // to respect that fact in the color selection
 3239     if((type & ASHIFT_LINE_MASK) == ASHIFT_LINE_VERTICAL_SELECTED)
 3240       points_idx[n].color = isflipped ? ASHIFT_LINECOLOR_BLUE : ASHIFT_LINECOLOR_GREEN;
 3241     else if((type & ASHIFT_LINE_MASK) == ASHIFT_LINE_VERTICAL_NOT_SELECTED)
 3242       points_idx[n].color = isflipped ? ASHIFT_LINECOLOR_YELLOW : ASHIFT_LINECOLOR_RED;
 3243     else if((type & ASHIFT_LINE_MASK) == ASHIFT_LINE_HORIZONTAL_SELECTED)
 3244       points_idx[n].color = isflipped ? ASHIFT_LINECOLOR_GREEN : ASHIFT_LINECOLOR_BLUE;
 3245     else if((type & ASHIFT_LINE_MASK) == ASHIFT_LINE_HORIZONTAL_NOT_SELECTED)
 3246       points_idx[n].color = isflipped ? ASHIFT_LINECOLOR_RED : ASHIFT_LINECOLOR_YELLOW;
 3247     else
 3248       points_idx[n].color = ASHIFT_LINECOLOR_GREY;
 3249   }
 3250 
 3251   return TRUE;
 3252 }
 3253 
 3254 // get all the points to display lines in the gui
 3255 static int get_points(struct dt_iop_module_t *self, const dt_iop_ashift_line_t *lines, const int lines_count,
 3256                       const int lines_version, float **points, dt_iop_ashift_points_idx_t **points_idx,
 3257                       int *points_lines_count, float scale)
 3258 {
 3259   dt_develop_t *dev = self->dev;
 3260   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 3261 
 3262   dt_iop_ashift_points_idx_t *my_points_idx = NULL;
 3263   float *my_points = NULL;
 3264 
 3265   // is the display flipped relative to the original image?
 3266   const int isflipped = g->isflipped;
 3267 
 3268   // allocate new index array
 3269   my_points_idx = (dt_iop_ashift_points_idx_t *)malloc(sizeof(dt_iop_ashift_points_idx_t) * lines_count);
 3270   if(my_points_idx == NULL) goto error;
 3271 
 3272   // account for total number of points
 3273   size_t total_points = 0;
 3274 
 3275   // first step: basic initialization of my_points_idx and counting of total_points
 3276   for(int n = 0; n < lines_count; n++)
 3277   {
 3278     const int length = lines[n].length;
 3279 
 3280     total_points += length;
 3281 
 3282     my_points_idx[n].length = length;
 3283     my_points_idx[n].near = 0;
 3284     my_points_idx[n].bounded = 0;
 3285 
 3286     const dt_iop_ashift_linetype_t type = lines[n].type;
 3287     my_points_idx[n].type = type;
 3288 
 3289     // set line color according to line type/orientation
 3290     // note: if the screen display is flipped versus the original image we need
 3291     // to respect that fact in the color selection
 3292     if((type & ASHIFT_LINE_MASK) == ASHIFT_LINE_VERTICAL_SELECTED)
 3293       my_points_idx[n].color = isflipped ? ASHIFT_LINECOLOR_BLUE : ASHIFT_LINECOLOR_GREEN;
 3294     else if((type & ASHIFT_LINE_MASK) == ASHIFT_LINE_VERTICAL_NOT_SELECTED)
 3295       my_points_idx[n].color = isflipped ? ASHIFT_LINECOLOR_YELLOW : ASHIFT_LINECOLOR_RED;
 3296     else if((type & ASHIFT_LINE_MASK) == ASHIFT_LINE_HORIZONTAL_SELECTED)
 3297       my_points_idx[n].color = isflipped ? ASHIFT_LINECOLOR_GREEN : ASHIFT_LINECOLOR_BLUE;
 3298     else if((type & ASHIFT_LINE_MASK) == ASHIFT_LINE_HORIZONTAL_NOT_SELECTED)
 3299       my_points_idx[n].color = isflipped ? ASHIFT_LINECOLOR_RED : ASHIFT_LINECOLOR_YELLOW;
 3300     else
 3301       my_points_idx[n].color = ASHIFT_LINECOLOR_GREY;
 3302   }
 3303 
 3304   // now allocate new points buffer
 3305   my_points = (float *)malloc(sizeof(float) * 2 * total_points);
 3306   if(my_points == NULL) goto error;
 3307 
 3308   // second step: generate points for each line
 3309   for(int n = 0, offset = 0; n < lines_count; n++)
 3310   {
 3311     my_points_idx[n].offset = offset;
 3312 
 3313     float x = lines[n].p1[0] / scale;
 3314     float y = lines[n].p1[1] / scale;
 3315     const int length = lines[n].length;
 3316 
 3317     const float dx = (lines[n].p2[0] / scale - x) / (float)(length - 1);
 3318     const float dy = (lines[n].p2[1] / scale - y) / (float)(length - 1);
 3319 
 3320     for(int l = 0; l < length && offset < total_points; l++, offset++)
 3321     {
 3322       my_points[2 * offset] = x;
 3323       my_points[2 * offset + 1] = y;
 3324 
 3325       x += dx;
 3326       y += dy;
 3327     }
 3328   }
 3329 
 3330   // third step: transform all points
 3331   if(!dt_dev_distort_transform_plus(dev, dev->preview_pipe, self->iop_order, DT_DEV_TRANSFORM_DIR_FORW_INCL, my_points, total_points))
 3332     goto error;
 3333 
 3334   // fourth step: get bounding box in final coordinates (used later for checking "near"-ness to mouse pointer)
 3335   for(int n = 0; n < lines_count; n++)
 3336   {
 3337     float xmin = FLT_MAX, xmax = FLT_MIN, ymin = FLT_MAX, ymax = FLT_MIN;
 3338 
 3339     size_t offset = my_points_idx[n].offset;
 3340     int length = my_points_idx[n].length;
 3341 
 3342     for(int l = 0; l < length; l++)
 3343     {
 3344       xmin = fmin(xmin, my_points[2 * offset]);
 3345       xmax = fmax(xmax, my_points[2 * offset]);
 3346       ymin = fmin(ymin, my_points[2 * offset + 1]);
 3347       ymax = fmax(ymax, my_points[2 * offset + 1]);
 3348     }
 3349 
 3350     my_points_idx[n].bbx = xmin;
 3351     my_points_idx[n].bbX = xmax;
 3352     my_points_idx[n].bby = ymin;
 3353     my_points_idx[n].bbY = ymax;
 3354   }
 3355 
 3356   // check if lines_version has changed in-between -> too bad: we can forget about all we did :(
 3357   if(g->lines_version > lines_version)
 3358     goto error;
 3359 
 3360   *points = my_points;
 3361   *points_idx = my_points_idx;
 3362   *points_lines_count = lines_count;
 3363 
 3364   return TRUE;
 3365 
 3366 error:
 3367   if(my_points_idx != NULL) free(my_points_idx);
 3368   if(my_points != NULL) free(my_points);
 3369   return FALSE;
 3370 }
 3371 
 3372 // does this gui have focus?
 3373 static int gui_has_focus(struct dt_iop_module_t *self)
 3374 {
 3375   return (self->dev->gui_module == self
 3376           && dt_dev_modulegroups_get_activated(darktable.develop) != DT_MODULEGROUP_BASICS);
 3377 }
 3378 
 3379 /* this function replaces this sentence, it calls distort_transform() for this module on the pipe
 3380 if(!dt_dev_distort_transform_plus(self->dev, self->dev->preview_pipe, self->priority, self->priority + 1,
 3381       (float *)V, 4))
 3382 */
 3383 static int call_distort_transform(dt_develop_t *dev, dt_dev_pixelpipe_t *pipe, struct dt_iop_module_t *self,
 3384                                   float *points, size_t points_count)
 3385 {
 3386   int ret = 0;
 3387   dt_dev_pixelpipe_iop_t *piece = dt_dev_distort_get_iop_pipe(self->dev, self->dev->preview_pipe, self);
 3388   if(!piece) return ret;
 3389   if(piece->module == self && /*piece->enabled && */  //see note below
 3390      !(dev->gui_module && dev->gui_module->operation_tags_filter() & piece->module->operation_tags()))
 3391   {
 3392     ret = piece->module->distort_transform(piece->module, piece, points, points_count);
 3393   }
 3394   return ret;
 3395   //NOTE: piece->enabled is FALSE for exactly the first mouse_moved event following a button_pressed event
 3396   //  when ASHIFT_CROP_ASPECT is active, which causes the first gui_post_expose call on starting to resize
 3397   //  the crop box to draw the center image without the crop overlay, resulting in an annoying visual glitch.
 3398   //  Removing the check appears to have no adverse effects and eliminates the glitch.
 3399 }
 3400 
 3401 void gui_post_expose(struct dt_iop_module_t *self, cairo_t *cr, int32_t width, int32_t height,
 3402                      int32_t pointerx, int32_t pointery)
 3403 {
 3404   dt_develop_t *dev = self->dev;
 3405   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 3406   dt_iop_ashift_params_t *p = (dt_iop_ashift_params_t *)self->params;
 3407 
 3408   // the usual rescaling stuff
 3409   const float wd = dev->preview_pipe->backbuf_width;
 3410   const float ht = dev->preview_pipe->backbuf_height;
 3411   if(wd < 1.0 || ht < 1.0) return;
 3412   const float pr_d = dev->preview_downsampling;
 3413   const float zoom_y = dt_control_get_dev_zoom_y();
 3414   const float zoom_x = dt_control_get_dev_zoom_x();
 3415   const dt_dev_zoom_t zoom = dt_control_get_dev_zoom();
 3416   const int closeup = dt_control_get_dev_closeup();
 3417   const float zoom_scale = dt_dev_get_zoom_scale(dev, zoom, 1<<closeup, 1);
 3418 
 3419   // we draw the cropping area; we need x_off/y_off/width/height which is only available
 3420   // after g->buf has been processed
 3421   if(g->buf && (p->cropmode != ASHIFT_CROP_OFF) && self->enabled)
 3422   {
 3423     // roi data of the preview pipe input buffer
 3424 
 3425     const float iwd = g->buf_width / pr_d;
 3426     const float iht = g->buf_height / pr_d;
 3427     const float ixo = g->buf_x_off / pr_d;
 3428     const float iyo = g->buf_y_off / pr_d;
 3429 
 3430     // the four corners of the input buffer of this module
 3431     float V[4][2] = { { ixo,        iyo       },
 3432                       { ixo,        iyo + iht },
 3433                       { ixo + iwd,  iyo + iht },
 3434                       { ixo + iwd,  iyo       } };
 3435 
 3436     // convert coordinates of corners to coordinates of this module's output
 3437     if(!call_distort_transform(self->dev, self->dev->preview_pipe, self, (float *)V, 4))
 3438       return;
 3439 
 3440     // get x/y-offset as well as width and height of output buffer
 3441     float xmin = FLT_MAX, ymin = FLT_MAX, xmax = FLT_MIN, ymax = FLT_MIN;
 3442     for(int n = 0; n < 4; n++)
 3443     {
 3444       xmin = MIN(xmin, V[n][0]);
 3445       xmax = MAX(xmax, V[n][0]);
 3446       ymin = MIN(ymin, V[n][1]);
 3447       ymax = MAX(ymax, V[n][1]);
 3448     }
 3449     const float owd = xmax - xmin;
 3450     const float oht = ymax - ymin;
 3451 
 3452     // the four clipping corners
 3453     float C[4][2] = { { xmin + g->cl * owd, ymin + g->ct * oht },
 3454                       { xmin + g->cl * owd, ymin + g->cb * oht },
 3455                       { xmin + g->cr * owd, ymin + g->cb * oht },
 3456                       { xmin + g->cr * owd, ymin + g->ct * oht } };
 3457 
 3458     // convert clipping corners to final output image
 3459     if(!dt_dev_distort_transform_plus(self->dev, self->dev->preview_pipe, self->iop_order,
 3460                                      DT_DEV_TRANSFORM_DIR_FORW_EXCL, (float *)C, 4))
 3461       return;
 3462 
 3463     cairo_save(cr);
 3464 
 3465     double dashes = DT_PIXEL_APPLY_DPI(5.0) / zoom_scale;
 3466     cairo_set_dash(cr, &dashes, 0, 0);
 3467 
 3468     cairo_rectangle(cr, 0, 0, width, height);
 3469     cairo_clip(cr);
 3470 
 3471     // mask parts of image outside of clipping area in dark grey
 3472     cairo_set_source_rgba(cr, .2, .2, .2, .8);
 3473     cairo_set_fill_rule(cr, CAIRO_FILL_RULE_EVEN_ODD);
 3474     cairo_rectangle(cr, 0, 0, width, height);
 3475     cairo_translate(cr, width / 2.0, height / 2.0);
 3476     cairo_scale(cr, zoom_scale, zoom_scale);
 3477     cairo_translate(cr, -.5f * wd - zoom_x * wd, -.5f * ht - zoom_y * ht);
 3478     cairo_move_to(cr, C[0][0], C[0][1]);
 3479     cairo_line_to(cr, C[1][0], C[1][1]);
 3480     cairo_line_to(cr, C[2][0], C[2][1]);
 3481     cairo_line_to(cr, C[3][0], C[3][1]);
 3482     cairo_close_path(cr);
 3483     cairo_fill(cr);
 3484 
 3485     // draw white outline around clipping area
 3486     dt_draw_set_color_overlay(cr, 0.7, 1.0);
 3487     cairo_set_line_width(cr, 2.0 / zoom_scale);
 3488     cairo_move_to(cr, C[0][0], C[0][1]);
 3489     cairo_line_to(cr, C[1][0], C[1][1]);
 3490     cairo_line_to(cr, C[2][0], C[2][1]);
 3491     cairo_line_to(cr, C[3][0], C[3][1]);
 3492     cairo_close_path(cr);
 3493     cairo_stroke(cr);
 3494 
 3495     // if adjusting crop, draw indicator
 3496     if (g->adjust_crop && p->cropmode == ASHIFT_CROP_ASPECT)
 3497     {
 3498       const double x1 = C[0][0];
 3499       const double x2 = fabs(x1 - C[1][0]) < 0.001f ? C[2][0] : C[1][0];
 3500       const double y1 = C[0][1];
 3501       const double y2 = fabs(y1 - C[1][1]) < 0.001f ? C[2][1] : C[1][1];
 3502 
 3503       const double xpos = (x1 + x2) / 2.0f;
 3504       const double ypos = (y1 + y2) / 2.0f;
 3505       const double base_size = fabs(x1 - x2);
 3506       const double size_circle = base_size / 30.0f;
 3507       const double size_line = base_size / 5.0f;
 3508       const double size_arrow = base_size / 25.0f;
 3509 
 3510       cairo_set_line_width(cr, 2.0 / zoom_scale);
 3511       dt_draw_set_color_overlay(cr, 0.7, 1.0);
 3512       cairo_arc (cr, xpos, ypos, size_circle, 0, 2.0 * M_PI);
 3513       cairo_stroke(cr);
 3514       cairo_fill(cr);
 3515 
 3516       cairo_set_line_width(cr, 2.0 / zoom_scale);
 3517       dt_draw_set_color_overlay(cr, 0.7, 1.0);
 3518 
 3519       // horizontal line
 3520       cairo_move_to(cr, xpos - size_line, ypos);
 3521       cairo_line_to(cr, xpos + size_line, ypos);
 3522 
 3523       cairo_move_to(cr, xpos - size_line, ypos);
 3524       cairo_rel_line_to(cr, size_arrow, size_arrow);
 3525       cairo_move_to(cr, xpos - size_line, ypos);
 3526       cairo_rel_line_to(cr, size_arrow, -size_arrow);
 3527 
 3528       cairo_move_to(cr, xpos + size_line, ypos);
 3529       cairo_rel_line_to(cr, -size_arrow, size_arrow);
 3530       cairo_move_to(cr, xpos + size_line, ypos);
 3531       cairo_rel_line_to(cr, -size_arrow, -size_arrow);
 3532 
 3533       // vertical line
 3534       cairo_move_to(cr, xpos, ypos - size_line);
 3535       cairo_line_to(cr, xpos, ypos + size_line);
 3536 
 3537       cairo_move_to(cr, xpos, ypos - size_line);
 3538       cairo_rel_line_to(cr, -size_arrow, size_arrow);
 3539       cairo_move_to(cr, xpos, ypos - size_line);
 3540       cairo_rel_line_to(cr, size_arrow, size_arrow);
 3541 
 3542       cairo_move_to(cr, xpos, ypos + size_line);
 3543       cairo_rel_line_to(cr, -size_arrow, -size_arrow);
 3544       cairo_move_to(cr, xpos, ypos + size_line);
 3545       cairo_rel_line_to(cr, size_arrow, -size_arrow);
 3546 
 3547       cairo_stroke(cr);
 3548     }
 3549 
 3550     cairo_restore(cr);
 3551   }
 3552 
 3553   // show guide lines on request
 3554   if(g->show_guides)
 3555   {
 3556     dt_guides_t *guide = (dt_guides_t *)darktable.guides->data;
 3557     double dashes = DT_PIXEL_APPLY_DPI(5.0);
 3558     cairo_save(cr);
 3559     cairo_rectangle(cr, 0, 0, width, height);
 3560     cairo_clip(cr);
 3561     cairo_set_line_width(cr, DT_PIXEL_APPLY_DPI(1.0));
 3562     cairo_set_source_rgb(cr, .8, .8, .8);
 3563     cairo_set_dash(cr, &dashes, 1, 0);
 3564     guide->draw(cr, 0, 0, width, height, 1.0, guide->user_data);
 3565     cairo_stroke_preserve(cr);
 3566     cairo_set_dash(cr, &dashes, 0, 0);
 3567     cairo_set_source_rgba(cr, 0.3, .3, .3, .8);
 3568     cairo_stroke(cr);
 3569     cairo_restore(cr);
 3570   }
 3571 
 3572   // structural data are currently being collected or fit procedure is running? -> skip
 3573   if(g->fitting) return;
 3574 
 3575   // no structural data or visibility switched off? -> stop here
 3576   if(g->lines == NULL || g->lines_suppressed || !gui_has_focus(self)) return;
 3577 
 3578   // get hash value that changes if distortions from here to the end of the pixelpipe changed
 3579   uint64_t hash = dt_dev_hash_distort(dev);
 3580   // get hash value that changes if coordinates of lines have changed
 3581   uint64_t lines_hash = get_lines_hash(g->lines, g->lines_count);
 3582 
 3583   // points data are missing or outdated, or distortion has changed?
 3584   if(g->points == NULL || g->points_idx == NULL || hash != g->grid_hash ||
 3585     (g->lines_version > g->points_version && g->lines_hash != lines_hash))
 3586   {
 3587     // we need to reprocess points
 3588     free(g->points);
 3589     g->points = NULL;
 3590     free(g->points_idx);
 3591     g->points_idx = NULL;
 3592     g->points_lines_count = 0;
 3593 
 3594     if(!get_points(self, g->lines, g->lines_count, g->lines_version, &g->points, &g->points_idx,
 3595                    &g->points_lines_count, pr_d))
 3596       return;
 3597 
 3598     g->points_version = g->lines_version;
 3599     g->grid_hash = hash;
 3600     g->lines_hash = lines_hash;
 3601   }
 3602   else if(g->lines_hash == lines_hash)
 3603   {
 3604     // update line type information in points_idx
 3605     for(int n = 0; n < g->points_lines_count; n++)
 3606       g->points_idx[n].type = g->lines[n].type;
 3607 
 3608     // coordinates of lines are unchanged -> we only need to update colors
 3609     if(!update_colors(self, g->points_idx, g->points_lines_count))
 3610       return;
 3611 
 3612     g->points_version = g->lines_version;
 3613   }
 3614 
 3615   // a final check
 3616   if(g->points == NULL || g->points_idx == NULL) return;
 3617 
 3618   cairo_save(cr);
 3619   cairo_rectangle(cr, 0, 0, width, height);
 3620   cairo_clip(cr);
 3621   cairo_translate(cr, width / 2.0, height / 2.0);
 3622   cairo_scale(cr, zoom_scale, zoom_scale);
 3623   cairo_translate(cr, -.5f * wd - zoom_x * wd, -.5f * ht - zoom_y * ht);
 3624 
 3625   // this must match the sequence of enum dt_iop_ashift_linecolor_t!
 3626   const float line_colors[5][4] =
 3627   { { 0.3f, 0.3f, 0.3f, 0.8f },                    // grey (misc. lines)
 3628     { 0.0f, 1.0f, 0.0f, 0.8f },                    // green (selected vertical lines)
 3629     { 0.8f, 0.0f, 0.0f, 0.8f },                    // red (de-selected vertical lines)
 3630     { 0.0f, 0.0f, 1.0f, 0.8f },                    // blue (selected horizontal lines)
 3631     { 0.8f, 0.8f, 0.0f, 0.8f } };                  // yellow (de-selected horizontal lines)
 3632 
 3633   cairo_set_line_cap(cr, CAIRO_LINE_CAP_ROUND);
 3634 
 3635   // now draw all lines
 3636   for(int n = 0; n < g->points_lines_count; n++)
 3637   {
 3638     // is the near flag set? -> draw line a bit thicker
 3639     if(g->points_idx[n].near)
 3640       cairo_set_line_width(cr, DT_PIXEL_APPLY_DPI(3.0) / zoom_scale);
 3641     else
 3642       cairo_set_line_width(cr, DT_PIXEL_APPLY_DPI(1.5) / zoom_scale);
 3643 
 3644     // the color of this line
 3645     const float *color = line_colors[g->points_idx[n].color];
 3646     cairo_set_source_rgba(cr, color[0], color[1], color[2], color[3]);
 3647 
 3648     size_t offset = g->points_idx[n].offset;
 3649     const int length = g->points_idx[n].length;
 3650 
 3651     // sanity check (this should not happen)
 3652     if(length < 2) continue;
 3653 
 3654     // set starting point of multi-segment line
 3655     cairo_move_to(cr, g->points[offset * 2], g->points[offset * 2 + 1]);
 3656 
 3657     offset++;
 3658     // draw individual line segments
 3659     for(int l = 1; l < length; l++, offset++)
 3660     {
 3661       cairo_line_to(cr, g->points[offset * 2], g->points[offset * 2 + 1]);
 3662     }
 3663 
 3664     // finally stroke the line
 3665     cairo_stroke(cr);
 3666   }
 3667 
 3668   // and we draw the selection box if any
 3669   if(g->isbounding != ASHIFT_BOUNDING_OFF)
 3670   {
 3671     float pzx = 0.0f, pzy = 0.0f;
 3672     dt_dev_get_pointer_zoom_pos(dev, pointerx, pointery, &pzx, &pzy);
 3673     pzx += 0.5f;
 3674     pzy += 0.5f;
 3675 
 3676     double dashed[] = { 4.0, 4.0 };
 3677     dashed[0] /= zoom_scale;
 3678     dashed[1] /= zoom_scale;
 3679     int len = sizeof(dashed) / sizeof(dashed[0]);
 3680 
 3681     cairo_rectangle(cr, g->lastx * wd, g->lasty * ht, (pzx - g->lastx) * wd,
 3682                    (pzy - g->lasty) * ht);
 3683     cairo_set_source_rgba(cr, .3, .3, .3, .8);
 3684     cairo_set_line_width(cr, 1.0 / zoom_scale);
 3685     cairo_set_dash(cr, dashed, len, 0);
 3686     cairo_stroke_preserve(cr);
 3687     cairo_set_source_rgba(cr, .8, .8, .8, .8);
 3688     cairo_set_dash(cr, dashed, len, 4);
 3689     cairo_stroke(cr);
 3690   }
 3691 
 3692   // indicate which area is used for "near"-ness detection when selecting/deselecting lines
 3693   if(g->near_delta > 0)
 3694   {
 3695     float pzx = 0.0f, pzy = 0.0f;
 3696     dt_dev_get_pointer_zoom_pos(dev, pointerx, pointery, &pzx, &pzy);
 3697     pzx += 0.5f;
 3698     pzy += 0.5f;
 3699 
 3700     double dashed[] = { 4.0, 4.0 };
 3701     dashed[0] /= zoom_scale;
 3702     dashed[1] /= zoom_scale;
 3703     int len = sizeof(dashed) / sizeof(dashed[0]);
 3704 
 3705     cairo_arc(cr, pzx * wd, pzy * ht, g->near_delta, 0, 2.0 * M_PI);
 3706 
 3707     cairo_set_source_rgba(cr, .3, .3, .3, .8);
 3708     cairo_set_line_width(cr, 1.0 / zoom_scale);
 3709     cairo_set_dash(cr, dashed, len, 0);
 3710     cairo_stroke_preserve(cr);
 3711     cairo_set_source_rgba(cr, .8, .8, .8, .8);
 3712     cairo_set_dash(cr, dashed, len, 4);
 3713     cairo_stroke(cr);
 3714   }
 3715 
 3716   cairo_restore(cr);
 3717 }
 3718 
 3719 // update the number of selected vertical and horizontal lines
 3720 static void update_lines_count(const dt_iop_ashift_line_t *lines, const int lines_count,
 3721                                int *vertical_count, int *horizontal_count)
 3722 {
 3723   int vlines = 0;
 3724   int hlines = 0;
 3725 
 3726   for(int n = 0; n < lines_count; n++)
 3727   {
 3728     if((lines[n].type & ASHIFT_LINE_MASK) == ASHIFT_LINE_VERTICAL_SELECTED)
 3729       vlines++;
 3730     else if((lines[n].type & ASHIFT_LINE_MASK) == ASHIFT_LINE_HORIZONTAL_SELECTED)
 3731       hlines++;
 3732   }
 3733 
 3734   *vertical_count = vlines;
 3735   *horizontal_count = hlines;
 3736 }
 3737 
 3738 int mouse_moved(struct dt_iop_module_t *self, double x, double y, double pressure, int which)
 3739 {
 3740   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 3741   int handled = 0;
 3742 
 3743   const float wd = self->dev->preview_pipe->backbuf_width;
 3744   const float ht = self->dev->preview_pipe->backbuf_height;
 3745   if(wd < 1.0 || ht < 1.0) return 1;
 3746 
 3747   float pzx = 0.0f, pzy = 0.0f;
 3748   dt_dev_get_pointer_zoom_pos(self->dev, x, y, &pzx, &pzy);
 3749   pzx += 0.5f;
 3750   pzy += 0.5f;
 3751 
 3752   if (g->adjust_crop)
 3753   {
 3754     dt_iop_ashift_params_t *p = (dt_iop_ashift_params_t *)self->params;
 3755 
 3756     float pts[4] = { pzx, pzy, 1.0f, 1.0f };
 3757     dt_dev_distort_backtransform_plus(self->dev, self->dev->preview_pipe, self->iop_order,
 3758                                       DT_DEV_TRANSFORM_DIR_FORW_INCL, pts, 2);
 3759 
 3760     const float newx = g->crop_cx + (pts[0] - pts[2]) - g->lastx;
 3761     const float newy = g->crop_cy + (pts[1] - pts[3]) - g->lasty;
 3762 
 3763     crop_adjust(self, p, newx, newy);
 3764     dt_control_queue_redraw_center();
 3765     return TRUE;
 3766   }
 3767 
 3768   // if visibility of lines is switched off or no lines available, we would normally adjust the crop box
 3769   // but since g->adjust_crop was FALSE, we have nothing to do
 3770   if(g->lines_suppressed || g->lines == NULL)
 3771     return TRUE;
 3772 
 3773   // if in rectangle selecting mode adjust "near"-ness of lines according to
 3774   // the rectangular selection
 3775   if(g->isbounding != ASHIFT_BOUNDING_OFF)
 3776   {
 3777     if(wd >= 1.0 && ht >= 1.0)
 3778     {
 3779       // mark lines inside the rectangle
 3780       get_bounded_inside(g->points, g->points_idx, g->points_lines_count, pzx * wd, pzy * ht, g->lastx * wd,
 3781                          g->lasty * ht, g->isbounding);
 3782     }
 3783 
 3784     dt_control_queue_redraw_center();
 3785     return FALSE;
 3786   }
 3787 
 3788   // gather information about "near"-ness in g->points_idx
 3789   get_near(g->points, g->points_idx, g->points_lines_count, pzx * wd, pzy * ht, g->near_delta);
 3790 
 3791   // if we are in sweeping mode iterate over lines as we move the pointer and change "selected" state.
 3792   if(g->isdeselecting || g->isselecting)
 3793   {
 3794     for(int n = 0; g->selecting_lines_version == g->lines_version && n < g->points_lines_count; n++)
 3795     {
 3796       if(g->points_idx[n].near == 0)
 3797         continue;
 3798 
 3799       if(g->isdeselecting)
 3800         g->lines[n].type &= ~ASHIFT_LINE_SELECTED;
 3801       else if(g->isselecting)
 3802         g->lines[n].type |= ASHIFT_LINE_SELECTED;
 3803 
 3804       handled = 1;
 3805     }
 3806   }
 3807 
 3808   if(handled)
 3809   {
 3810     update_lines_count(g->lines, g->lines_count, &g->vertical_count, &g->horizontal_count);
 3811     g->lines_version++;
 3812     g->selecting_lines_version++;
 3813   }
 3814 
 3815   dt_control_queue_redraw_center();
 3816 
 3817   // if not in sweeping mode we need to pass the event
 3818   return (g->isdeselecting || g->isselecting);
 3819 }
 3820 
 3821 int button_pressed(struct dt_iop_module_t *self, double x, double y, double pressure, int which, int type,
 3822                    uint32_t state)
 3823 {
 3824   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 3825   int handled = 0;
 3826 
 3827   float pzx = 0.0f, pzy = 0.0f;
 3828   dt_dev_get_pointer_zoom_pos(self->dev, x, y, &pzx, &pzy);
 3829   pzx += 0.5f;
 3830   pzy += 0.5f;
 3831 
 3832   const float wd = self->dev->preview_pipe->backbuf_width;
 3833   const float ht = self->dev->preview_pipe->backbuf_height;
 3834   if(wd < 1.0 || ht < 1.0) return 1;
 3835 
 3836 
 3837   // if visibility of lines is switched off or no lines available -> potentially adjust crop area
 3838   if(g->lines_suppressed || g->lines == NULL)
 3839   {
 3840     dt_iop_ashift_params_t *p = (dt_iop_ashift_params_t *)self->params;
 3841     if (p->cropmode == ASHIFT_CROP_ASPECT)
 3842     {
 3843       dt_control_change_cursor(GDK_HAND1);
 3844       g->adjust_crop = TRUE;
 3845 
 3846       float pts[4] = { pzx, pzy, 1.0f, 1.0f };
 3847       dt_dev_distort_backtransform_plus(self->dev, self->dev->preview_pipe, self->iop_order,
 3848                                         DT_DEV_TRANSFORM_DIR_FORW_INCL, pts, 2);
 3849 
 3850       g->lastx = pts[0] - pts[2];
 3851       g->lasty = pts[1] - pts[3];
 3852       g->crop_cx = 0.5f * (g->cl + g->cr);
 3853       g->crop_cy = 0.5f * (g->ct + g->cb);
 3854       return TRUE;
 3855     }
 3856     else
 3857       return FALSE;
 3858   }
 3859 
 3860   // remember lines version at this stage so we can continuously monitor if the
 3861   // lines have changed in-between
 3862   g->selecting_lines_version = g->lines_version;
 3863 
 3864   // if shift button is pressed go into bounding mode (selecting or deselecting
 3865   // in a rectangle area)
 3866   if(dt_modifier_is(state, GDK_SHIFT_MASK))
 3867   {
 3868     g->lastx = pzx;
 3869     g->lasty = pzy;
 3870 
 3871     g->isbounding = (which == 3) ? ASHIFT_BOUNDING_DESELECT : ASHIFT_BOUNDING_SELECT;
 3872     dt_control_change_cursor(GDK_CROSS);
 3873 
 3874     return TRUE;
 3875   }
 3876 
 3877   dt_dev_zoom_t zoom = dt_control_get_dev_zoom();
 3878   const int closeup = dt_control_get_dev_closeup();
 3879   const float min_scale = dt_dev_get_zoom_scale(self->dev, DT_ZOOM_FIT, 1<<closeup, 0);
 3880   const float cur_scale = dt_dev_get_zoom_scale(self->dev, zoom, 1<<closeup, 0);
 3881 
 3882   // if we are zoomed out (no panning possible) and we have lines to display we take control
 3883   const int take_control = (cur_scale == min_scale) && (g->points_lines_count > 0);
 3884 
 3885   g->near_delta = dt_conf_get_float("plugins/darkroom/ashift/near_delta");
 3886 
 3887   // gather information about "near"-ness in g->points_idx
 3888   get_near(g->points, g->points_idx, g->points_lines_count, pzx * wd, pzy * ht, g->near_delta);
 3889 
 3890   // iterate over all lines close to the pointer and change "selected" state.
 3891   // left-click selects and right-click deselects the line
 3892   for(int n = 0; g->selecting_lines_version == g->lines_version && n < g->points_lines_count; n++)
 3893   {
 3894     if(g->points_idx[n].near == 0)
 3895       continue;
 3896 
 3897     if(which == 3)
 3898       g->lines[n].type &= ~ASHIFT_LINE_SELECTED;
 3899     else
 3900       g->lines[n].type |= ASHIFT_LINE_SELECTED;
 3901 
 3902     handled = 1;
 3903   }
 3904 
 3905   // we switch into sweeping mode either if we anyhow take control
 3906   // or if cursor was close to a line when button was pressed. in other
 3907   // cases we hand over the event (for image panning)
 3908   if((take_control || handled) && which == 3)
 3909   {
 3910     dt_control_change_cursor(GDK_PIRATE);
 3911     g->isdeselecting = 1;
 3912   }
 3913   else if(take_control || handled)
 3914   {
 3915     dt_control_change_cursor(GDK_PLUS);
 3916     g->isselecting = 1;
 3917   }
 3918 
 3919   if(handled)
 3920   {
 3921     update_lines_count(g->lines, g->lines_count, &g->vertical_count, &g->horizontal_count);
 3922     g->lines_version++;
 3923     g->selecting_lines_version++;
 3924   }
 3925 
 3926   return (take_control || handled);
 3927 }
 3928 
 3929 int button_released(struct dt_iop_module_t *self, double x, double y, int which, uint32_t state)
 3930 {
 3931   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 3932 
 3933   dt_control_change_cursor(GDK_LEFT_PTR);
 3934   if (g->adjust_crop)
 3935   {
 3936     // stop adjust crop
 3937     g->adjust_crop = FALSE;
 3938     dt_iop_ashift_params_t *p = (dt_iop_ashift_params_t *)self->params;
 3939     swap_shadow_crop_box(p,g);  // temporarily update the crop box in p
 3940     dt_dev_add_history_item(darktable.develop, self, TRUE);
 3941     swap_shadow_crop_box(p,g);  // restore p
 3942   }
 3943 
 3944   // finalize the isbounding mode
 3945   // if user has released the shift button in-between -> do nothing
 3946   if(g->isbounding != ASHIFT_BOUNDING_OFF && dt_modifier_is(state, GDK_SHIFT_MASK))
 3947   {
 3948     int handled = 0;
 3949 
 3950     // we compute the rectangle selection
 3951     float pzx = 0.0f, pzy = 0.0f;
 3952     dt_dev_get_pointer_zoom_pos(self->dev, x, y, &pzx, &pzy);
 3953 
 3954     pzx += 0.5f;
 3955     pzy += 0.5f;
 3956 
 3957     const float wd = self->dev->preview_pipe->backbuf_width;
 3958     const float ht = self->dev->preview_pipe->backbuf_height;
 3959 
 3960     if(wd >= 1.0 && ht >= 1.0)
 3961     {
 3962       // mark lines inside the rectangle
 3963       get_bounded_inside(g->points, g->points_idx, g->points_lines_count, pzx * wd, pzy * ht, g->lastx * wd,
 3964                          g->lasty * ht, g->isbounding);
 3965 
 3966       // select or deselect lines within the rectangle according to isbounding state
 3967       for(int n = 0; g->selecting_lines_version == g->lines_version && n < g->points_lines_count; n++)
 3968       {
 3969         if(g->points_idx[n].bounded == 0) continue;
 3970 
 3971         if(g->isbounding == ASHIFT_BOUNDING_DESELECT)
 3972           g->lines[n].type &= ~ASHIFT_LINE_SELECTED;
 3973         else
 3974           g->lines[n].type |= ASHIFT_LINE_SELECTED;
 3975 
 3976         handled = 1;
 3977       }
 3978 
 3979       if(handled)
 3980       {
 3981         update_lines_count(g->lines, g->lines_count, &g->vertical_count, &g->horizontal_count);
 3982         g->lines_version++;
 3983         g->selecting_lines_version++;
 3984       }
 3985 
 3986     dt_control_queue_redraw_center();
 3987     }
 3988   }
 3989 
 3990   // end of sweeping/isbounding mode
 3991   g->isselecting = g->isdeselecting = 0;
 3992   g->isbounding = ASHIFT_BOUNDING_OFF;
 3993   g->near_delta = 0;
 3994   g->lastx = g->lasty = -1.0f;
 3995   g->crop_cx = g->crop_cy = -1.0f;
 3996 
 3997   return 0;
 3998 }
 3999 
 4000 int scrolled(struct dt_iop_module_t *self, double x, double y, int up, uint32_t state)
 4001 {
 4002   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 4003 
 4004   // do nothing if visibility of lines is switched off or no lines available
 4005   if(g->lines_suppressed || g->lines == NULL)
 4006     return FALSE;
 4007 
 4008   if(g->near_delta > 0 && (g->isdeselecting || g->isselecting))
 4009   {
 4010     int handled = 0;
 4011 
 4012     float pzx = 0.0f, pzy = 0.0f;
 4013     dt_dev_get_pointer_zoom_pos(self->dev, x, y, &pzx, &pzy);
 4014     pzx += 0.5f;
 4015     pzy += 0.5f;
 4016 
 4017     const float wd = self->dev->preview_pipe->backbuf_width;
 4018     const float ht = self->dev->preview_pipe->backbuf_height;
 4019 
 4020     float near_delta = dt_conf_get_float("plugins/darkroom/ashift/near_delta");
 4021     const float amount = up ? 0.8f : 1.25f;
 4022     near_delta = MAX(4.0f, MIN(near_delta * amount, 100.0f));
 4023     dt_conf_set_float("plugins/darkroom/ashift/near_delta", near_delta);
 4024     g->near_delta = near_delta;
 4025 
 4026     // gather information about "near"-ness in g->points_idx
 4027     get_near(g->points, g->points_idx, g->points_lines_count, pzx * wd, pzy * ht, g->near_delta);
 4028 
 4029     // iterate over all lines close to the pointer and change "selected" state.
 4030     for(int n = 0; g->selecting_lines_version == g->lines_version && n < g->points_lines_count; n++)
 4031     {
 4032       if(g->points_idx[n].near == 0)
 4033         continue;
 4034 
 4035       if(g->isdeselecting)
 4036         g->lines[n].type &= ~ASHIFT_LINE_SELECTED;
 4037       else if(g->isselecting)
 4038         g->lines[n].type |= ASHIFT_LINE_SELECTED;
 4039 
 4040       handled = 1;
 4041     }
 4042 
 4043     if(handled)
 4044     {
 4045       update_lines_count(g->lines, g->lines_count, &g->vertical_count, &g->horizontal_count);
 4046       g->lines_version++;
 4047       g->selecting_lines_version++;
 4048     }
 4049 
 4050     dt_control_queue_redraw_center();
 4051     return TRUE;
 4052   }
 4053 
 4054   return FALSE;
 4055 }
 4056 
 4057 
 4058 void gui_changed(dt_iop_module_t *self, GtkWidget *w, void *previous)
 4059 {
 4060   dt_iop_ashift_params_t *p = (dt_iop_ashift_params_t *)self->params;
 4061   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 4062 
 4063 #ifdef ASHIFT_DEBUG
 4064   model_probe(self, p, g->lastfit);
 4065 #endif
 4066   do_crop(self, p);
 4067   commit_crop_box(p,g);
 4068 
 4069   if(w == g->mode)
 4070   {
 4071     gtk_widget_set_visible(g->specifics, p->mode == ASHIFT_MODE_SPECIFIC);
 4072   }
 4073 }
 4074 
 4075 static void cropmode_callback(GtkWidget *widget, gpointer user_data)
 4076 {
 4077   if(darktable.gui->reset) return;
 4078 
 4079   dt_iop_module_t *self = (dt_iop_module_t *)user_data;
 4080   dt_iop_ashift_params_t *p = (dt_iop_ashift_params_t *)self->params;
 4081   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 4082 
 4083   if(g->lines != NULL && !g->lines_suppressed)
 4084   {
 4085     g->lines_suppressed = 1;
 4086     gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(g->eye), g->lines_suppressed);
 4087   }
 4088 
 4089   swap_shadow_crop_box(p,g);    //temporarily update real crop box
 4090   dt_dev_add_history_item(darktable.develop, self, TRUE);
 4091   swap_shadow_crop_box(p,g);
 4092 }
 4093 
 4094 static void guide_lines_callback(GtkWidget *widget, gpointer user_data)
 4095 {
 4096   dt_iop_module_t *self = (dt_iop_module_t *)user_data;
 4097   if(darktable.gui->reset) return;
 4098   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 4099   g->show_guides = dt_bauhaus_combobox_get(widget);
 4100   dt_iop_request_focus(self);
 4101   dt_control_queue_redraw_center();
 4102 }
 4103 
 4104 static int fit_v_button_clicked(GtkWidget *widget, GdkEventButton *event, gpointer user_data)
 4105 {
 4106   dt_iop_module_t *self = (dt_iop_module_t *)user_data;
 4107   if(darktable.gui->reset) return FALSE;
 4108 
 4109   if(event->button == 1)
 4110   {
 4111     dt_iop_ashift_params_t *p = (dt_iop_ashift_params_t *)self->params;
 4112     dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 4113 
 4114     const int control = dt_modifiers_include(event->state, GDK_CONTROL_MASK);
 4115     const int shift = dt_modifiers_include(event->state, GDK_SHIFT_MASK);
 4116 
 4117     dt_iop_ashift_fitaxis_t fitaxis = ASHIFT_FIT_NONE;
 4118 
 4119     if(control)
 4120       g->lastfit = fitaxis = ASHIFT_FIT_ROTATION_VERTICAL_LINES;
 4121     else if(shift)
 4122       g->lastfit = fitaxis = ASHIFT_FIT_VERTICALLY_NO_ROTATION;
 4123     else
 4124       g->lastfit = fitaxis = ASHIFT_FIT_VERTICALLY;
 4125 
 4126     dt_iop_request_focus(self);
 4127 
 4128     if(self->enabled)
 4129     {
 4130       // module is enable -> we process directly
 4131       if(do_fit(self, p, fitaxis))
 4132       {
 4133         ++darktable.gui->reset;
 4134         dt_bauhaus_slider_set_soft(g->rotation, p->rotation);
 4135         dt_bauhaus_slider_set_soft(g->lensshift_v, p->lensshift_v);
 4136         dt_bauhaus_slider_set_soft(g->lensshift_h, p->lensshift_h);
 4137         dt_bauhaus_slider_set_soft(g->shear, p->shear);
 4138         --darktable.gui->reset;
 4139       }
 4140     }
 4141     else
 4142     {
 4143       // module is not enabled -> invoke it and queue the job to be processed once
 4144       // the preview image is ready
 4145       g->jobcode = ASHIFT_JOBCODE_FIT;
 4146       g->jobparams = g->lastfit = fitaxis;
 4147       p->toggle ^= 1;
 4148     }
 4149 
 4150     dt_dev_add_history_item(darktable.develop, self, TRUE); //also calls dt_control_queue_redraw_center
 4151     return TRUE;
 4152   }
 4153   return FALSE;
 4154 }
 4155 
 4156 static int fit_h_button_clicked(GtkWidget *widget, GdkEventButton *event, gpointer user_data)
 4157 {
 4158   dt_iop_module_t *self = (dt_iop_module_t *)user_data;
 4159   if(darktable.gui->reset) return FALSE;
 4160 
 4161   if(event->button == 1)
 4162   {
 4163     dt_iop_ashift_params_t *p = (dt_iop_ashift_params_t *)self->params;
 4164     dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 4165 
 4166     const int control = dt_modifiers_include(event->state, GDK_CONTROL_MASK);
 4167     const int shift = dt_modifiers_include(event->state, GDK_SHIFT_MASK);
 4168 
 4169     dt_iop_ashift_fitaxis_t fitaxis = ASHIFT_FIT_NONE;
 4170 
 4171     if(control)
 4172       g->lastfit = fitaxis = ASHIFT_FIT_ROTATION_HORIZONTAL_LINES;
 4173     else if(shift)
 4174       g->lastfit = fitaxis = ASHIFT_FIT_HORIZONTALLY_NO_ROTATION;
 4175     else
 4176       g->lastfit = fitaxis = ASHIFT_FIT_HORIZONTALLY;
 4177 
 4178     dt_iop_request_focus(self);
 4179 
 4180     if(self->enabled)
 4181     {
 4182       // module is enable -> we process directly
 4183       if(do_fit(self, p, fitaxis))
 4184       {
 4185         ++darktable.gui->reset;
 4186         dt_bauhaus_slider_set_soft(g->rotation, p->rotation);
 4187         dt_bauhaus_slider_set_soft(g->lensshift_v, p->lensshift_v);
 4188         dt_bauhaus_slider_set_soft(g->lensshift_h, p->lensshift_h);
 4189         dt_bauhaus_slider_set_soft(g->shear, p->shear);
 4190         --darktable.gui->reset;
 4191       }
 4192     }
 4193     else
 4194     {
 4195       // module is not enabled -> invoke it and queue the job to be processed once
 4196       // the preview image is ready
 4197       g->jobcode = ASHIFT_JOBCODE_FIT;
 4198       g->jobparams = g->lastfit = fitaxis;
 4199       p->toggle ^= 1;
 4200     }
 4201 
 4202     dt_dev_add_history_item(darktable.develop, self, TRUE); //also calls dt_control_queue_redraw_center
 4203     return TRUE;
 4204   }
 4205   return FALSE;
 4206 }
 4207 
 4208 static int fit_both_button_clicked(GtkWidget *widget, GdkEventButton *event, gpointer user_data)
 4209 {
 4210   dt_iop_module_t *self = (dt_iop_module_t *)user_data;
 4211   if(darktable.gui->reset) return FALSE;
 4212 
 4213   if(event->button == 1)
 4214   {
 4215     dt_iop_ashift_params_t *p = (dt_iop_ashift_params_t *)self->params;
 4216     dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 4217 
 4218     const int control = dt_modifiers_include(event->state, GDK_CONTROL_MASK);
 4219     const int shift = dt_modifiers_include(event->state, GDK_SHIFT_MASK);
 4220 
 4221     dt_iop_ashift_fitaxis_t fitaxis = ASHIFT_FIT_NONE;
 4222 
 4223     if(control && shift)
 4224       fitaxis = ASHIFT_FIT_BOTH;
 4225     else if(control)
 4226       fitaxis = ASHIFT_FIT_ROTATION_BOTH_LINES;
 4227     else if(shift)
 4228       fitaxis = ASHIFT_FIT_BOTH_NO_ROTATION;
 4229     else
 4230       fitaxis = ASHIFT_FIT_BOTH_SHEAR;
 4231 
 4232     dt_iop_request_focus(self);
 4233 
 4234     if(self->enabled)
 4235     {
 4236       // module is enable -> we process directly
 4237       if(do_fit(self, p, fitaxis))
 4238       {
 4239         ++darktable.gui->reset;
 4240         dt_bauhaus_slider_set_soft(g->rotation, p->rotation);
 4241         dt_bauhaus_slider_set_soft(g->lensshift_v, p->lensshift_v);
 4242         dt_bauhaus_slider_set_soft(g->lensshift_h, p->lensshift_h);
 4243         dt_bauhaus_slider_set_soft(g->shear, p->shear);
 4244         --darktable.gui->reset;
 4245       }
 4246     }
 4247     else
 4248     {
 4249       // module is not enabled -> invoke it and queue the job to be processed once
 4250       // the preview image is ready
 4251       g->jobcode = ASHIFT_JOBCODE_FIT;
 4252       g->jobparams = g->lastfit = fitaxis;
 4253       p->toggle ^= 1;
 4254     }
 4255 
 4256     dt_dev_add_history_item(darktable.develop, self, TRUE); //also calls dt_control_queue_redraw_center
 4257     return TRUE;
 4258   }
 4259   return FALSE;
 4260 }
 4261 
 4262 static int structure_button_clicked(GtkWidget *widget, GdkEventButton *event, gpointer user_data)
 4263 {
 4264   dt_iop_module_t *self = (dt_iop_module_t *)user_data;
 4265   if(darktable.gui->reset) return FALSE;
 4266 
 4267   if(event->button == 1)
 4268   {
 4269     dt_iop_ashift_params_t *p = (dt_iop_ashift_params_t *)self->params;
 4270     dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 4271 
 4272     const int control = dt_modifiers_include(event->state, GDK_CONTROL_MASK);
 4273     const int shift = dt_modifiers_include(event->state, GDK_SHIFT_MASK);
 4274 
 4275     dt_iop_ashift_enhance_t enhance;
 4276 
 4277     if(control && shift)
 4278       enhance = ASHIFT_ENHANCE_EDGES | ASHIFT_ENHANCE_DETAIL;
 4279     else if(shift)
 4280       enhance = ASHIFT_ENHANCE_DETAIL;
 4281     else if(control)
 4282       enhance = ASHIFT_ENHANCE_EDGES;
 4283     else
 4284       enhance = ASHIFT_ENHANCE_NONE;
 4285 
 4286     dt_iop_request_focus(self);
 4287 
 4288     if(self->enabled)
 4289     {
 4290       // module is enabled -> process directly
 4291       (void)do_get_structure(self, p, enhance);
 4292     }
 4293     else
 4294     {
 4295       // module is not enabled -> invoke it and queue the job to be processed once
 4296       // the preview image is ready
 4297       g->jobcode = ASHIFT_JOBCODE_GET_STRUCTURE;
 4298       g->jobparams = enhance;
 4299       p->toggle ^= 1;
 4300     }
 4301 
 4302     dt_dev_add_history_item(darktable.develop, self, TRUE); // also calls dt_control_queue_redraw_center
 4303     return TRUE;
 4304   }
 4305   return FALSE;
 4306 }
 4307 
 4308 static void clean_button_clicked(GtkButton *button, gpointer user_data)
 4309 {
 4310   dt_iop_module_t *self = (dt_iop_module_t *)user_data;
 4311   if(darktable.gui->reset) return;
 4312   dt_iop_ashift_params_t *p = (dt_iop_ashift_params_t *)self->params;
 4313   (void)do_clean_structure(self, p);
 4314   dt_iop_request_focus(self);
 4315   dt_control_queue_redraw_center();
 4316 }
 4317 
 4318 static void eye_button_toggled(GtkToggleButton *togglebutton, gpointer user_data)
 4319 {
 4320   dt_iop_module_t *self = (dt_iop_module_t *)user_data;
 4321   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 4322   if(darktable.gui->reset) return;
 4323   if(g->lines == NULL)
 4324   {
 4325     g->lines_suppressed = 0;
 4326     gtk_toggle_button_set_active(togglebutton, 0);
 4327   }
 4328   else
 4329   {
 4330     g->lines_suppressed = gtk_toggle_button_get_active(togglebutton);
 4331   }
 4332   dt_iop_request_focus(self);
 4333   dt_control_queue_redraw_center();
 4334 }
 4335 
 4336 // routine that is called after preview image has been processed. we use it
 4337 // to perform structure collection or fitting in case those have been triggered while
 4338 // the module had not yet been enabled
 4339 static void process_after_preview_callback(gpointer instance, gpointer user_data)
 4340 {
 4341   dt_iop_module_t *self = (dt_iop_module_t *)user_data;
 4342   dt_iop_ashift_params_t *p = (dt_iop_ashift_params_t *)self->params;
 4343   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 4344 
 4345   dt_iop_ashift_jobcode_t jobcode = g->jobcode;
 4346   int jobparams = g->jobparams;
 4347 
 4348   // purge
 4349   g->jobcode = ASHIFT_JOBCODE_NONE;
 4350   g->jobparams = 0;
 4351 
 4352   if(darktable.gui->reset) return;
 4353 
 4354   switch(jobcode)
 4355   {
 4356     case ASHIFT_JOBCODE_GET_STRUCTURE:
 4357       (void)do_get_structure(self, p, (dt_iop_ashift_enhance_t)jobparams);
 4358       break;
 4359 
 4360     case ASHIFT_JOBCODE_FIT:
 4361       if(do_fit(self, p, (dt_iop_ashift_fitaxis_t)jobparams))
 4362       {
 4363         ++darktable.gui->reset;
 4364         dt_bauhaus_slider_set_soft(g->rotation, p->rotation);
 4365         dt_bauhaus_slider_set_soft(g->lensshift_v, p->lensshift_v);
 4366         dt_bauhaus_slider_set_soft(g->lensshift_h, p->lensshift_h);
 4367         dt_bauhaus_slider_set_soft(g->shear, p->shear);
 4368         --darktable.gui->reset;
 4369       }
 4370       dt_dev_add_history_item(darktable.develop, self, TRUE);
 4371       break;
 4372 
 4373     case ASHIFT_JOBCODE_NONE:
 4374     default:
 4375       break;
 4376   }
 4377 
 4378   dt_control_queue_redraw_center();
 4379 }
 4380 
 4381 void commit_params(struct dt_iop_module_t *self, dt_iop_params_t *p1, dt_dev_pixelpipe_t *pipe,
 4382                    dt_dev_pixelpipe_iop_t *piece)
 4383 {
 4384   dt_iop_ashift_params_t *p = (dt_iop_ashift_params_t *)p1;
 4385   dt_iop_ashift_data_t *d = (dt_iop_ashift_data_t *)piece->data;
 4386 
 4387   d->rotation = p->rotation;
 4388   d->lensshift_v = p->lensshift_v;
 4389   d->lensshift_h = p->lensshift_h;
 4390   d->shear = p->shear;
 4391   d->f_length_kb = (p->mode == ASHIFT_MODE_GENERIC) ? DEFAULT_F_LENGTH : p->f_length * p->crop_factor;
 4392   d->orthocorr = (p->mode == ASHIFT_MODE_GENERIC) ? 0.0f : p->orthocorr;
 4393   d->aspect = (p->mode == ASHIFT_MODE_GENERIC) ? 1.0f : p->aspect;
 4394 
 4395   if(gui_has_focus(self))
 4396   {
 4397     // if gui has focus we want to see the full uncropped image
 4398     d->cl = 0.0f;
 4399     d->cr = 1.0f;
 4400     d->ct = 0.0f;
 4401     d->cb = 1.0f;
 4402   }
 4403   else
 4404   {
 4405     d->cl = p->cl;
 4406     d->cr = p->cr;
 4407     d->ct = p->ct;
 4408     d->cb = p->cb;
 4409   }
 4410 }
 4411 
 4412 void init_pipe(struct dt_iop_module_t *self, dt_dev_pixelpipe_t *pipe, dt_dev_pixelpipe_iop_t *piece)
 4413 {
 4414   dt_iop_ashift_data_t *d = (dt_iop_ashift_data_t *)calloc(1, sizeof(dt_iop_ashift_data_t));
 4415   piece->data = (void *)d;
 4416 }
 4417 
 4418 void cleanup_pipe(struct dt_iop_module_t *self, dt_dev_pixelpipe_t *pipe, dt_dev_pixelpipe_iop_t *piece)
 4419 {
 4420   free(piece->data);
 4421   piece->data = NULL;
 4422 }
 4423 
 4424 void gui_update(struct dt_iop_module_t *self)
 4425 {
 4426   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 4427   dt_iop_ashift_params_t *p = (dt_iop_ashift_params_t *)self->params;
 4428 
 4429   dt_bauhaus_slider_set_soft(g->rotation, p->rotation);
 4430   dt_bauhaus_slider_set_soft(g->lensshift_v, p->lensshift_v);
 4431   dt_bauhaus_slider_set_soft(g->lensshift_h, p->lensshift_h);
 4432   dt_bauhaus_slider_set_soft(g->shear, p->shear);
 4433   dt_bauhaus_slider_set_soft(g->f_length, p->f_length);
 4434   dt_bauhaus_slider_set_soft(g->crop_factor, p->crop_factor);
 4435   dt_bauhaus_slider_set(g->orthocorr, p->orthocorr);
 4436   dt_bauhaus_slider_set(g->aspect, p->aspect);
 4437   dt_bauhaus_combobox_set(g->mode, p->mode);
 4438   dt_bauhaus_combobox_set(g->guide_lines, g->show_guides);
 4439   dt_bauhaus_combobox_set(g->cropmode, p->cropmode);
 4440   gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(g->eye), 0);
 4441 
 4442   gtk_widget_set_visible(g->specifics, p->mode == ASHIFT_MODE_SPECIFIC);
 4443 
 4444   // copy crop box into shadow variables
 4445   shadow_crop_box(p,g);
 4446 }
 4447 
 4448 void reload_defaults(dt_iop_module_t *module)
 4449 {
 4450   // our module is disabled by default
 4451   module->default_enabled = 0;
 4452 
 4453   int isflipped = 0;
 4454   float f_length = DEFAULT_F_LENGTH;
 4455   float crop_factor = 1.0f;
 4456 
 4457   // try to get information on orientation, focal length and crop factor from image data
 4458   if(module->dev)
 4459   {
 4460     const dt_image_t *img = &module->dev->image_storage;
 4461     // orientation only needed as a-priori information to correctly label some sliders
 4462     // before pixelpipe has been set up. later we will get a definite result by
 4463     // assessing the pixelpipe
 4464     isflipped = (img->orientation == ORIENTATION_ROTATE_CCW_90_DEG
 4465                  || img->orientation == ORIENTATION_ROTATE_CW_90_DEG) ? 1 : 0;
 4466 
 4467     // focal length should be available in exif data if lens is electronically coupled to the camera
 4468     f_length = isfinite(img->exif_focal_length) && img->exif_focal_length > 0.0f ? img->exif_focal_length : f_length;
 4469     // crop factor of the camera is often not available and user will need to set it manually in the gui
 4470     crop_factor = isfinite(img->exif_crop) && img->exif_crop > 0.0f ? img->exif_crop : crop_factor;
 4471   }
 4472 
 4473   // init defaults:
 4474   ((dt_iop_ashift_params_t *)module->default_params)->f_length = f_length;
 4475   ((dt_iop_ashift_params_t *)module->default_params)->crop_factor = crop_factor;
 4476 
 4477   // reset gui elements
 4478   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)module->gui_data;
 4479   if(g)
 4480   {
 4481 
 4482     char string_v[256];
 4483     char string_h[256];
 4484 
 4485     snprintf(string_v, sizeof(string_v), _("lens shift (%s)"), isflipped ? _("horizontal") : _("vertical"));
 4486     snprintf(string_h, sizeof(string_h), _("lens shift (%s)"), isflipped ? _("vertical") : _("horizontal"));
 4487 
 4488     dt_bauhaus_widget_set_label(g->lensshift_v, NULL, string_v);
 4489     dt_bauhaus_widget_set_label(g->lensshift_h, NULL, string_h);
 4490 
 4491     dt_bauhaus_slider_set_default(g->f_length, f_length);
 4492     dt_bauhaus_slider_set_default(g->crop_factor, crop_factor);
 4493 
 4494     dt_iop_gui_enter_critical_section(module);
 4495     free(g->buf);
 4496     g->buf = NULL;
 4497     g->buf_width = 0;
 4498     g->buf_height = 0;
 4499     g->buf_x_off = 0;
 4500     g->buf_y_off = 0;
 4501     g->buf_scale = 1.0f;
 4502     g->buf_hash = 0;
 4503     g->isflipped = -1;
 4504     g->lastfit = ASHIFT_FIT_NONE;
 4505     dt_iop_gui_leave_critical_section(module);
 4506 
 4507     g->fitting = 0;
 4508     free(g->lines);
 4509     g->lines = NULL;
 4510     g->lines_count =0;
 4511     g->horizontal_count = 0;
 4512     g->vertical_count = 0;
 4513     g->grid_hash = 0;
 4514     g->lines_hash = 0;
 4515     g->rotation_range = ROTATION_RANGE_SOFT;
 4516     g->lensshift_v_range = LENSSHIFT_RANGE_SOFT;
 4517     g->lensshift_h_range = LENSSHIFT_RANGE_SOFT;
 4518     g->shear_range = SHEAR_RANGE_SOFT;
 4519     g->lines_suppressed = 0;
 4520     g->lines_version = 0;
 4521     g->show_guides = 0;
 4522     g->isselecting = 0;
 4523     g->isdeselecting = 0;
 4524     g->isbounding = ASHIFT_BOUNDING_OFF;
 4525     g->near_delta = 0;
 4526     g->selecting_lines_version = 0;
 4527 
 4528     free(g->points);
 4529     g->points = NULL;
 4530     free(g->points_idx);
 4531     g->points_idx = NULL;
 4532     g->points_lines_count = 0;
 4533     g->points_version = 0;
 4534 
 4535     g->jobcode = ASHIFT_JOBCODE_NONE;
 4536     g->jobparams = 0;
 4537     g->adjust_crop = FALSE;
 4538     g->lastx = g->lasty = -1.0f;
 4539     g->crop_cx = g->crop_cy = 1.0f;
 4540   }
 4541 }
 4542 
 4543 
 4544 void init_global(dt_iop_module_so_t *module)
 4545 {
 4546   dt_iop_ashift_global_data_t *gd
 4547       = (dt_iop_ashift_global_data_t *)malloc(sizeof(dt_iop_ashift_global_data_t));
 4548   module->data = gd;
 4549 
 4550   const int program = 2; // basic.cl, from programs.conf
 4551   gd->kernel_ashift_bilinear = dt_opencl_create_kernel(program, "ashift_bilinear");
 4552   gd->kernel_ashift_bicubic = dt_opencl_create_kernel(program, "ashift_bicubic");
 4553   gd->kernel_ashift_lanczos2 = dt_opencl_create_kernel(program, "ashift_lanczos2");
 4554   gd->kernel_ashift_lanczos3 = dt_opencl_create_kernel(program, "ashift_lanczos3");
 4555 }
 4556 
 4557 void cleanup_global(dt_iop_module_so_t *module)
 4558 {
 4559   dt_iop_ashift_global_data_t *gd = (dt_iop_ashift_global_data_t *)module->data;
 4560   dt_opencl_free_kernel(gd->kernel_ashift_bilinear);
 4561   dt_opencl_free_kernel(gd->kernel_ashift_bicubic);
 4562   dt_opencl_free_kernel(gd->kernel_ashift_lanczos2);
 4563   dt_opencl_free_kernel(gd->kernel_ashift_lanczos3);
 4564   free(module->data);
 4565   module->data = NULL;
 4566 }
 4567 
 4568 // adjust labels of lens shift parameters according to flip status of image
 4569 static gboolean draw(GtkWidget *widget, cairo_t *cr, dt_iop_module_t *self)
 4570 {
 4571   dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 4572   if(darktable.gui->reset) return FALSE;
 4573 
 4574   dt_iop_gui_enter_critical_section(self);
 4575   const int isflipped = g->isflipped;
 4576   dt_iop_gui_leave_critical_section(self);
 4577 
 4578   if(isflipped == -1) return FALSE;
 4579 
 4580   char string_v[256];
 4581   char string_h[256];
 4582 
 4583   snprintf(string_v, sizeof(string_v), _("lens shift (%s)"), isflipped ? _("horizontal") : _("vertical"));
 4584   snprintf(string_h, sizeof(string_h), _("lens shift (%s)"), isflipped ? _("vertical") : _("horizontal"));
 4585 
 4586   ++darktable.gui->reset;
 4587   dt_bauhaus_widget_set_label(g->lensshift_v, NULL, string_v);
 4588   dt_bauhaus_widget_set_label(g->lensshift_h, NULL, string_h);
 4589   gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(g->eye), g->lines_suppressed);
 4590   --darktable.gui->reset;
 4591 
 4592   return FALSE;
 4593 }
 4594 
 4595 static void _event_preview_updated_callback(gpointer instance, dt_iop_module_t *self)
 4596 {
 4597   if(self->dev->gui_module != self)
 4598   {
 4599     dt_image_update_final_size(self->dev->preview_pipe->output_imgid);
 4600   }
 4601   DT_DEBUG_CONTROL_SIGNAL_DISCONNECT(darktable.signals, G_CALLBACK(_event_preview_updated_callback), self);
 4602 }
 4603 
 4604 void gui_focus(struct dt_iop_module_t *self, gboolean in)
 4605 {
 4606   if(self->enabled)
 4607   {
 4608     dt_iop_ashift_params_t *p = (dt_iop_ashift_params_t *)self->params;
 4609     dt_iop_ashift_gui_data_t *g = (dt_iop_ashift_gui_data_t *)self->gui_data;
 4610     if (in)
 4611     {
 4612       shadow_crop_box(p,g);
 4613       dt_control_queue_redraw_center();
 4614     }
 4615     else
 4616     {
 4617       // once the pipe is recomputed, we want to update final sizes
 4618       DT_DEBUG_CONTROL_SIGNAL_CONNECT(darktable.signals, DT_SIGNAL_DEVELOP_PREVIEW_PIPE_FINISHED,
 4619                                       G_CALLBACK(_event_preview_updated_callback), self);
 4620       commit_crop_box(p,g);
 4621     }
 4622   }
 4623 }
 4624 
 4625 static float log10_curve(GtkWidget *self, float inval, dt_bauhaus_curve_t dir)
 4626 {
 4627   float outval;
 4628   if(dir == DT_BAUHAUS_SET)
 4629   {
 4630     outval = log10f(inval * 999.0f + 1.0f) / 3.0f;
 4631   }
 4632   else
 4633   {
 4634     outval = (expf(M_LN10 * inval * 3.0f) - 1.0f) / 999.0f;
 4635   }
 4636   return outval;
 4637 }
 4638 
 4639 static float log2_curve(GtkWidget *self, float inval, dt_bauhaus_curve_t dir)
 4640 {
 4641   float outval;
 4642   if(dir == DT_BAUHAUS_SET)
 4643   {
 4644       outval = log2f(inval * 1.5f + 0.5f) / 2.0f + 0.5f;
 4645   }
 4646   else
 4647   {
 4648     outval = (exp2f(inval * 2.0 - 1.0) - 0.5f) / 1.5f;
 4649   }
 4650   return outval;
 4651 }
 4652 
 4653 void gui_init(struct dt_iop_module_t *self)
 4654 {
 4655   dt_iop_ashift_gui_data_t *g = IOP_GUI_ALLOC(ashift);
 4656 
 4657   dt_iop_gui_enter_critical_section(self); //not actually needed, we're the only one with a pointer to this instance
 4658   g->buf = NULL;
 4659   g->buf_width = 0;
 4660   g->buf_height = 0;
 4661   g->buf_x_off = 0;
 4662   g->buf_y_off = 0;
 4663   g->buf_scale = 1.0f;
 4664   g->buf_hash = 0;
 4665   g->isflipped = -1;
 4666   g->lastfit = ASHIFT_FIT_NONE;
 4667   dt_iop_gui_leave_critical_section(self);
 4668 
 4669   g->fitting = 0;
 4670   g->lines = NULL;
 4671   g->lines_count = 0;
 4672   g->vertical_count = 0;
 4673   g->horizontal_count = 0;
 4674   g->lines_version = 0;
 4675   g->lines_suppressed = 0;
 4676   g->points = NULL;
 4677   g->points_idx = NULL;
 4678   g->points_lines_count = 0;
 4679   g->points_version = 0;
 4680   g->grid_hash = 0;
 4681   g->lines_hash = 0;
 4682   g->rotation_range = ROTATION_RANGE_SOFT;
 4683   g->lensshift_v_range = LENSSHIFT_RANGE_SOFT;
 4684   g->lensshift_h_range = LENSSHIFT_RANGE_SOFT;
 4685   g->shear_range = SHEAR_RANGE_SOFT;
 4686   g->show_guides = 0;
 4687   g->isselecting = 0;
 4688   g->isdeselecting = 0;
 4689   g->isbounding = ASHIFT_BOUNDING_OFF;
 4690   g->near_delta = 0;
 4691   g->selecting_lines_version = 0;
 4692 
 4693   g->jobcode = ASHIFT_JOBCODE_NONE;
 4694   g->jobparams = 0;
 4695   g->adjust_crop = FALSE;
 4696   g->lastx = g->lasty = -1.0f;
 4697   g->crop_cx = g->crop_cy = 1.0f;
 4698 
 4699   g->rotation = dt_bauhaus_slider_from_params(self, N_("rotation"));
 4700   dt_bauhaus_slider_set_format(g->rotation, "%.2f°");
 4701   dt_bauhaus_slider_set_soft_range(g->rotation, -ROTATION_RANGE, ROTATION_RANGE);
 4702 
 4703   g->lensshift_v = dt_bauhaus_slider_from_params(self, "lensshift_v");
 4704   dt_bauhaus_slider_set_soft_range(g->lensshift_v, -LENSSHIFT_RANGE, LENSSHIFT_RANGE);
 4705   dt_bauhaus_slider_set_digits(g->lensshift_v, 3);
 4706 
 4707   g->lensshift_h = dt_bauhaus_slider_from_params(self, "lensshift_h");
 4708   dt_bauhaus_slider_set_soft_range(g->lensshift_h, -LENSSHIFT_RANGE, LENSSHIFT_RANGE);
 4709   dt_bauhaus_slider_set_digits(g->lensshift_h, 3);
 4710 
 4711   g->shear = dt_bauhaus_slider_from_params(self, "shear");
 4712   dt_bauhaus_slider_set_soft_range(g->shear, -SHEAR_RANGE, SHEAR_RANGE);
 4713 
 4714   g->guide_lines = dt_bauhaus_combobox_new(self);
 4715   dt_bauhaus_widget_set_label(g->guide_lines, NULL, N_("guides"));
 4716   dt_bauhaus_combobox_add(g->guide_lines, _("off"));
 4717   dt_bauhaus_combobox_add(g->guide_lines, _("on"));
 4718   gtk_box_pack_start(GTK_BOX(self->widget), g->guide_lines, TRUE, TRUE, 0);
 4719 
 4720   g->cropmode = dt_bauhaus_combobox_from_params(self, "cropmode");
 4721   g_signal_connect(G_OBJECT(g->cropmode), "value-changed", G_CALLBACK(cropmode_callback), self);
 4722 
 4723   g->mode = dt_bauhaus_combobox_from_params(self, "mode");
 4724 
 4725   GtkWidget *saved_widget = self->widget;
 4726   self->widget = g->specifics = gtk_box_new(GTK_ORIENTATION_VERTICAL, DT_BAUHAUS_SPACE);
 4727 
 4728   g->f_length = dt_bauhaus_slider_from_params(self, "f_length");
 4729   dt_bauhaus_slider_set_soft_range(g->f_length, 10.0f, 1000.0f);
 4730   dt_bauhaus_slider_set_curve(g->f_length, log10_curve);
 4731   dt_bauhaus_slider_set_format(g->f_length, "%.0fmm");
 4732   dt_bauhaus_slider_set_step(g->f_length, 1.0);
 4733 
 4734   g->crop_factor = dt_bauhaus_slider_from_params(self, "crop_factor");
 4735   dt_bauhaus_slider_set_soft_range(g->crop_factor, 1.0f, 2.0f);
 4736 
 4737   g->orthocorr = dt_bauhaus_slider_from_params(self, "orthocorr");
 4738   dt_bauhaus_slider_set_format(g->orthocorr, "%.0f%%");
 4739   // this parameter could serve to finetune between generic model (0%) and specific model (100%).
 4740   // however, users can more easily get the same effect with the aspect adjust parameter so we keep
 4741   // this one hidden.
 4742   gtk_widget_set_no_show_all(g->orthocorr, TRUE);
 4743   gtk_widget_set_visible(g->orthocorr, FALSE);
 4744 
 4745   g->aspect = dt_bauhaus_slider_from_params(self, "aspect");
 4746   dt_bauhaus_slider_set_curve(g->aspect, log2_curve);
 4747 
 4748   self->widget = saved_widget;
 4749   gtk_box_pack_start(GTK_BOX(self->widget), g->specifics, TRUE, TRUE, 0);
 4750 
 4751   GtkGrid *grid = GTK_GRID(gtk_grid_new());
 4752   gtk_grid_set_row_spacing(grid, 2 * DT_BAUHAUS_SPACE);
 4753   gtk_grid_set_column_spacing(grid, DT_PIXEL_APPLY_DPI(10));
 4754 
 4755   gtk_grid_attach(grid, dt_ui_label_new(_("automatic fit")), 0, 0, 1, 1);
 4756 
 4757   g->fit_v = dtgtk_button_new(dtgtk_cairo_paint_perspective, CPF_STYLE_FLAT | 1, NULL);
 4758   gtk_widget_set_hexpand(GTK_WIDGET(g->fit_v), TRUE);
 4759   gtk_grid_attach(grid, g->fit_v, 1, 0, 1, 1);
 4760 
 4761   g->fit_h = dtgtk_button_new(dtgtk_cairo_paint_perspective, CPF_STYLE_FLAT | 2, NULL);
 4762   gtk_widget_set_hexpand(GTK_WIDGET(g->fit_h), TRUE);
 4763   gtk_grid_attach(grid, g->fit_h, 2, 0, 1, 1);
 4764 
 4765   g->fit_both = dtgtk_button_new(dtgtk_cairo_paint_perspective, CPF_STYLE_FLAT | 3, NULL);
 4766   gtk_widget_set_hexpand(GTK_WIDGET(g->fit_both), TRUE);
 4767   gtk_grid_attach(grid, g->fit_both, 3, 0, 1, 1);
 4768 
 4769   gtk_grid_attach(grid, dt_ui_label_new(_("get structure")), 0, 1, 1, 1);
 4770 
 4771   g->structure = dtgtk_button_new(dtgtk_cairo_paint_structure, CPF_STYLE_FLAT, NULL);
 4772   gtk_widget_set_hexpand(GTK_WIDGET(g->structure), TRUE);
 4773   gtk_grid_attach(grid, g->structure, 1, 1, 1, 1);
 4774 
 4775   g->clean = dtgtk_button_new(dtgtk_cairo_paint_cancel, CPF_STYLE_FLAT, NULL);
 4776   gtk_widget_set_hexpand(GTK_WIDGET(g->clean), TRUE);
 4777   gtk_grid_attach(grid, g->clean, 2, 1, 1, 1);
 4778 
 4779   g->eye = dtgtk_togglebutton_new(dtgtk_cairo_paint_eye_toggle, CPF_STYLE_FLAT, NULL);
 4780   gtk_widget_set_hexpand(GTK_WIDGET(g->eye), TRUE);
 4781   gtk_grid_attach(grid, g->eye, 3, 1, 1, 1);
 4782 
 4783   gtk_box_pack_start(GTK_BOX(self->widget), GTK_WIDGET(grid), TRUE, TRUE, 0);
 4784 
 4785   gtk_widget_set_tooltip_text(g->rotation, _("rotate image"));