1/*2 This file is part of darktable,3 Copyright (C) 2016-2026 darktable developers.45 darktable is free software: you can redistribute it and/or modify6 it under the terms of the GNU General Public License as published by7 the Free Software Foundation, either version 3 of the License, or8 (at your option) any later version.910 darktable is distributed in the hope that it will be useful,11 but WITHOUT ANY WARRANTY; without even the implied warranty of12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the13 GNU General Public License for more details.1415 You should have received a copy of the GNU General Public License16 along with darktable. If not, see <http://www.gnu.org/licenses/>.17*/1819// Extracted and adapted for Upright, 2026-09-21. See NOTICE.md.2021#include "compat.h"2223#define ROTATION_RANGE 10 // allowed min/max default range for rotation parameter24#define ROTATION_RANGE_SOFT 180 // allowed min/max range for rotation parameter with manual adjustment25#define LENSSHIFT_RANGE 1.0 // allowed min/max default range for lensshift parameters26#define LENSSHIFT_RANGE_SOFT 2.0 // allowed min/max range for lensshift parameters with manual adjustment27#define SHEAR_RANGE 0.2 // allowed min/max range for shear parameter28#define SHEAR_RANGE_SOFT 0.5 // allowed min/max range for shear parameter with manual adjustment29#define MIN_LINE_LENGTH 5 // the minimum length of a line in pixels to be regarded as relevant30#define MAX_TANGENTIAL_DEVIATION 30 // by how many degrees a line may deviate from the +/-180 and +/-90 to be regarded as relevant31#define LSD_SCALE 0.99 // LSD: scaling factor for line detection32#define LSD_SIGMA_SCALE 0.6 // LSD: sigma for Gaussian filter is computed as sigma = sigma_scale/scale33#define LSD_QUANT 2.0 // LSD: bound to the quantization error on the gradient norm34#define LSD_ANG_TH 22.5 // LSD: gradient angle tolerance in degrees35#define LSD_LOG_EPS 0.0 // LSD: detection threshold: -log10(NFA) > log_eps36#define LSD_DENSITY_TH 0.7 // LSD: minimal density of region points in rectangle37#define LSD_N_BINS 1024 // LSD: number of bins in pseudo-ordering of gradient modulus38#define LSD_GAMMA 0.45 // gamma correction to apply on raw images prior to line detection39#define RANSAC_RUNS 400 // how many iterations to run in ransac40#define RANSAC_EPSILON 2 // starting value for ransac epsilon (in -log10 units)41#define RANSAC_EPSILON_STEP 1 // step size of epsilon optimization (log10 units)42#define RANSAC_ELIMINATION_RATIO 60 // percentage of lines we try to eliminate as outliers43#define RANSAC_OPTIMIZATION_STEPS 5 // home many steps to optimize epsilon44#define RANSAC_OPTIMIZATION_DRY_RUNS 50 // how man runs per optimization steps45#define RANSAC_HURDLE 5 // hurdle rate: the number of lines below which we do a complete permutation instead of random sampling46#define MINIMUM_FITLINES 2 // minimum number of lines needed for automatic parameter fit47#define NMS_EPSILON 1e-3 // break criterion for Nelder-Mead simplex48#define NMS_SCALE 1.0 // scaling factor for Nelder-Mead simplex49#define NMS_ITERATIONS 400 // number of iterations for Nelder-Mead simplex50#define NMS_CROP_EPSILON 100.0 // break criterion for Nelder-Mead simplex on crop fitting51#define NMS_CROP_SCALE 0.5 // scaling factor for Nelder-Mead simplex on crop fitting52#define NMS_CROP_ITERATIONS 100 // number of iterations for Nelder-Mead simplex on crop fitting53#define NMS_ALPHA 1.0 // reflection coefficient for Nelder-Mead simplex54#define NMS_BETA 0.5 // contraction coefficient for Nelder-Mead simplex55#define NMS_GAMMA 2.0 // expansion coefficient for Nelder-Mead simplex56#define DEFAULT_F_LENGTH 28.0 // focal length we assume if no exif data are available5758// define to get debugging output59#undef ASHIFT_DEBUG6061// maximum number of drawn lines that can be saved in parameters62// any change in this value needs to upgrade parameters version !63#define MAX_SAVED_LINES 50646566#include "lsd.inc"67#include "upstream/ashift_nmsimplex.c"6869typedef enum dt_iop_ashift_method_t70{71 ASHIFT_METHOD_NONE = 0,72 ASHIFT_METHOD_AUTO,73 ASHIFT_METHOD_QUAD,74 ASHIFT_METHOD_LINES75} dt_iop_ashift_method_t;7677typedef enum dt_iop_ashift_homodir_t78{79 ASHIFT_HOMOGRAPH_FORWARD,80 ASHIFT_HOMOGRAPH_INVERTED81} dt_iop_ashift_homodir_t;8283typedef enum dt_iop_ashift_linetype_t84{85 ASHIFT_LINE_IRRELEVANT = 0, // the line is found to be not interesting86 // eg. too short, or not horizontal or vertical87 ASHIFT_LINE_RELEVANT = 1 << 0, // the line is relevant for us88 ASHIFT_LINE_DIRVERT = 1 << 1, // the line is (mostly) vertical, else (mostly) horizontal89 ASHIFT_LINE_SELECTED = 1 << 2, // the line is selected for fitting90 ASHIFT_LINE_VERTICAL_NOT_SELECTED = ASHIFT_LINE_RELEVANT | ASHIFT_LINE_DIRVERT,91 ASHIFT_LINE_HORIZONTAL_NOT_SELECTED = ASHIFT_LINE_RELEVANT,92 ASHIFT_LINE_VERTICAL_SELECTED = ASHIFT_LINE_RELEVANT | ASHIFT_LINE_DIRVERT | ASHIFT_LINE_SELECTED,93 ASHIFT_LINE_HORIZONTAL_SELECTED = ASHIFT_LINE_RELEVANT | ASHIFT_LINE_SELECTED,94 ASHIFT_LINE_MASK = ASHIFT_LINE_RELEVANT | ASHIFT_LINE_DIRVERT | ASHIFT_LINE_SELECTED95} dt_iop_ashift_linetype_t;9697typedef enum dt_iop_ashift_linecolor_t98{99 ASHIFT_LINECOLOR_GREY = 0,100 ASHIFT_LINECOLOR_GREEN = 1,101 ASHIFT_LINECOLOR_RED = 2,102 ASHIFT_LINECOLOR_BLUE = 3,103 ASHIFT_LINECOLOR_YELLOW = 4104} dt_iop_ashift_linecolor_t;105106typedef enum dt_iop_ashift_fitaxis_t107{108 ASHIFT_FIT_NONE = 0, // none109 ASHIFT_FIT_ROTATION = 1 << 0, // flag indicates to fit rotation angle110 ASHIFT_FIT_LENS_VERT = 1 << 1, // flag indicates to fit vertical lens shift111 ASHIFT_FIT_LENS_HOR = 1 << 2, // flag indicates to fit horizontal lens shift112 ASHIFT_FIT_SHEAR = 1 << 3, // flag indicates to fit shear parameter113 ASHIFT_FIT_LINES_VERT = 1 << 4, // use vertical lines for fitting114 ASHIFT_FIT_LINES_HOR = 1 << 5, // use horizontal lines for fitting115 ASHIFT_FIT_LENS_BOTH = ASHIFT_FIT_LENS_VERT | ASHIFT_FIT_LENS_HOR,116 ASHIFT_FIT_LINES_BOTH = ASHIFT_FIT_LINES_VERT | ASHIFT_FIT_LINES_HOR,117 ASHIFT_FIT_VERTICALLY = ASHIFT_FIT_ROTATION | ASHIFT_FIT_LENS_VERT | ASHIFT_FIT_LINES_VERT,118 ASHIFT_FIT_HORIZONTALLY = ASHIFT_FIT_ROTATION | ASHIFT_FIT_LENS_HOR | ASHIFT_FIT_LINES_HOR,119 ASHIFT_FIT_BOTH = ASHIFT_FIT_ROTATION | ASHIFT_FIT_LENS_VERT | ASHIFT_FIT_LENS_HOR |120 ASHIFT_FIT_LINES_VERT | ASHIFT_FIT_LINES_HOR,121 ASHIFT_FIT_VERTICALLY_NO_ROTATION = ASHIFT_FIT_LENS_VERT | ASHIFT_FIT_LINES_VERT,122 ASHIFT_FIT_HORIZONTALLY_NO_ROTATION = ASHIFT_FIT_LENS_HOR | ASHIFT_FIT_LINES_HOR,123 ASHIFT_FIT_BOTH_NO_ROTATION = ASHIFT_FIT_LENS_VERT | ASHIFT_FIT_LENS_HOR |124 ASHIFT_FIT_LINES_VERT | ASHIFT_FIT_LINES_HOR,125 ASHIFT_FIT_BOTH_SHEAR = ASHIFT_FIT_ROTATION | ASHIFT_FIT_LENS_VERT | ASHIFT_FIT_LENS_HOR |126 ASHIFT_FIT_SHEAR | ASHIFT_FIT_LINES_VERT | ASHIFT_FIT_LINES_HOR,127 ASHIFT_FIT_ROTATION_VERTICAL_LINES = ASHIFT_FIT_ROTATION | ASHIFT_FIT_LINES_VERT,128 ASHIFT_FIT_ROTATION_HORIZONTAL_LINES = ASHIFT_FIT_ROTATION | ASHIFT_FIT_LINES_HOR,129 ASHIFT_FIT_ROTATION_BOTH_LINES = ASHIFT_FIT_ROTATION | ASHIFT_FIT_LINES_VERT | ASHIFT_FIT_LINES_HOR,130 ASHIFT_FIT_FLIP = ASHIFT_FIT_LENS_VERT | ASHIFT_FIT_LENS_HOR | ASHIFT_FIT_LINES_VERT | ASHIFT_FIT_LINES_HOR131} dt_iop_ashift_fitaxis_t;132133typedef enum dt_iop_ashift_nmsresult_t134{135 NMS_SUCCESS = 0,136 NMS_NOT_ENOUGH_LINES = 1,137 NMS_DID_NOT_CONVERGE = 2,138 NMS_INSANE = 3139} dt_iop_ashift_nmsresult_t;140141typedef enum dt_iop_ashift_enhance_t142{143 ASHIFT_ENHANCE_NONE = 0,144 ASHIFT_ENHANCE_EDGES = 1 << 0,145 ASHIFT_ENHANCE_DETAIL = 1 << 1,146 ASHIFT_ENHANCE_HORIZONTAL = 0x100,147 ASHIFT_ENHANCE_VERTICAL = 0x200148} dt_iop_ashift_enhance_t;149150typedef enum dt_iop_ashift_mode_t151{152 ASHIFT_MODE_GENERIC = 0, // $DESCRIPTION: "generic"153 ASHIFT_MODE_SPECIFIC = 1 // $DESCRIPTION: "specific"154} dt_iop_ashift_mode_t;155156typedef enum dt_iop_ashift_crop_t157{158 ASHIFT_CROP_OFF = 0, // $DESCRIPTION: "off"159 ASHIFT_CROP_LARGEST = 1,// $DESCRIPTION: "largest area"160 ASHIFT_CROP_ASPECT = 2 // $DESCRIPTION: "original format"161} dt_iop_ashift_crop_t;162163typedef enum dt_iop_ashift_bounding_t164{165 ASHIFT_BOUNDING_OFF = 0,166 ASHIFT_BOUNDING_SELECT = 1,167 ASHIFT_BOUNDING_DESELECT = 2168} dt_iop_ashift_bounding_t;169170typedef enum dt_iop_ashift_jobcode_t171{172 ASHIFT_JOBCODE_NONE = 0,173 ASHIFT_JOBCODE_GET_STRUCTURE = 1,174 ASHIFT_JOBCODE_FIT = 2,175 ASHIFT_JOBCODE_GET_STRUCTURE_LINES = 3,176 ASHIFT_JOBCODE_GET_STRUCTURE_QUAD = 4,177 ASHIFT_JOBCODE_DO_CROP = 5178} dt_iop_ashift_jobcode_t;179180typedef struct dt_iop_ashift_params_t181{182 float rotation; // $MIN: -ROTATION_RANGE_SOFT $MAX: ROTATION_RANGE_SOFT $DEFAULT: 0.0183 float lensshift_v; // $MIN: -LENSSHIFT_RANGE_SOFT $MAX: LENSSHIFT_RANGE_SOFT $DEFAULT: 0.0 $DESCRIPTION: "lens shift (vertical)"184 float lensshift_h; // $MIN: -LENSSHIFT_RANGE_SOFT $MAX: LENSSHIFT_RANGE_SOFT $DEFAULT: 0.0 $DESCRIPTION: "lens shift (horizontal)"185 float shear; // $MIN: -SHEAR_RANGE_SOFT $MAX: SHEAR_RANGE_SOFT $DEFAULT: 0.0 $DESCRIPTION: "shear"186 float f_length; // $MIN: 1.0 $MAX: 2000.0 $DEFAULT: DEFAULT_F_LENGTH $DESCRIPTION: "focal length"187 float crop_factor; // $MIN: 0.5 $MAX: 10.0 $DEFAULT: 1.0 $DESCRIPTION: "crop factor"188 float orthocorr; // $MIN: 0.0 $MAX: 100.0 $DEFAULT: 100.0 $DESCRIPTION: "lens dependence"189 float aspect; // $MIN: 0.5 $MAX: 2.0 $DEFAULT: 1.0 $DESCRIPTION: "aspect adjust"190 dt_iop_ashift_mode_t mode; // $DEFAULT: ASHIFT_MODE_GENERIC $DESCRIPTION: "lens model"191 dt_iop_ashift_crop_t cropmode; // $DEFAULT: ASHIFT_CROP_LARGEST $DESCRIPTION: "automatic cropping"192 float cl; // $DEFAULT: 0.0 $MIN: 0.0 $MAX: 1.0193 float cr; // $DEFAULT: 1.0 $MIN: 0.0 $MAX: 1.0194 float ct; // $DEFAULT: 0.0 $MIN: 0.0 $MAX: 1.0195 float cb; // $DEFAULT: 1.0 $MIN: 0.0 $MAX: 1.0196 float last_drawn_lines[MAX_SAVED_LINES * 4];197 int last_drawn_lines_count;198 float last_quad_lines[8];199} dt_iop_ashift_params_t;200201typedef struct dt_iop_ashift_line_t202{203 float p1[3];204 float p2[3];205 float length;206 float width;207 float weight;208 dt_iop_ashift_linetype_t type;209 // homogeneous coordinates:210 float DT_ALIGNED_PIXEL L[3];211} dt_iop_ashift_line_t;212213typedef struct dt_iop_ashift_points_idx_t214{215 size_t offset;216 int length;217 int near;218 int bounded;219 dt_iop_ashift_linetype_t type;220 dt_iop_ashift_linecolor_t color;221 // bounding box:222 float bbx, bby, bbX, bbY;223} dt_iop_ashift_points_idx_t;224225typedef struct dt_iop_ashift_fit_params_t226{227 int params_count;228 dt_iop_ashift_linetype_t linetype;229 dt_iop_ashift_linetype_t linemask;230 dt_iop_ashift_line_t *lines;231 int lines_count;232 int width;233 int height;234 float weight;235 float f_length_kb;236 float orthocorr;237 float aspect;238 float rotation;239 float lensshift_v;240 float lensshift_h;241 float shear;242 float rotation_range;243 float lensshift_v_range;244 float lensshift_h_range;245 float shear_range;246} dt_iop_ashift_fit_params_t;247248249typedef struct {250 dt_iop_ashift_line_t *lines;251 int lines_count, lines_in_width, lines_in_height, lines_x_off, lines_y_off;252 int vertical_count, horizontal_count, lines_version, isflipped;253 float vertical_weight, horizontal_weight;254 float rotation_range, lensshift_v_range, lensshift_h_range, shear_range;255} dt_iop_ashift_gui_data_t;256typedef struct { dt_iop_ashift_gui_data_t *gui_data; } dt_iop_module_t;257258// normalized product of two 3x1 vectors259// dst needs to be different from v1 and v2260static inline void vec3prodn(float *dst, const float *const v1, const float *const v2)261{262 const float l1 = v1[1] * v2[2] - v1[2] * v2[1];263 const float l2 = v1[2] * v2[0] - v1[0] * v2[2];264 const float l3 = v1[0] * v2[1] - v1[1] * v2[0];265266 // normalize so that l1^2 + l2^2 + l3^3 = 1267 const float sq = sqrtf(l1 * l1 + l2 * l2 + l3 * l3);268269 const float f = sq > 0.0f ? 1.0f / sq : 1.0f;270271 dst[0] = l1 * f;272 dst[1] = l2 * f;273 dst[2] = l3 * f;274}275276// normalize a 3x1 vector so that x^2 + y^2 + z^2 = 1277// dst and v may be the same278static inline void vec3norm(float *dst, const float *const v)279{280 const float sq = sqrtf(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);281282 // special handling for an all-zero vector283 const float f = sq > 0.0f ? 1.0f / sq : 1.0f;284285 dst[0] = v[0] * f;286 dst[1] = v[1] * f;287 dst[2] = v[2] * f;288}289290// normalize a 3x1 vector so that x^2 + y^2 = 1; a useful normalization for291// lines in homogeneous coordinates292// dst and v may be the same293static inline void vec3lnorm(float *dst, const float *const v)294{295 const float sq = dt_fast_hypotf(v[0], v[1]);296297 // special handling for a point vector of the image center298 const float f = sq > 0.0f ? 1.0f / sq : 1.0f;299300 dst[0] = v[0] * f;301 dst[1] = v[1] * f;302 dst[2] = v[2] * f;303}304305306// scalar product of two 3x1 vectors307static inline float vec3scalar(const float *const v1, const float *const v2)308{309 return (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]);310}311312// check if 3x1 vector is (very close to) null313static inline int vec3isnull(const float *const v)314{315 const float eps = 1e-10f;316 return (fabsf(v[0]) < eps && fabsf(v[1]) < eps && fabsf(v[2]) < eps);317}318319320#define MAT3SWAP(a, b) { float (*tmp)[3] = (a); (a) = (b); (b) = tmp; }321322static void _homography(float *homograph,323 const float angle,324 const float shift_v,325 const float shift_h,326 const float shear,327 const float f_length_kb,328 const float orthocorr,329 const float aspect,330 const int width,331 const int height,332 const dt_iop_ashift_homodir_t dir)333{334 // calculate homograph that combines all translations, rotations335 // and warping into one single matrix operation.336 // this is heavily leaning on ShiftN where the homographic matrix expects337 // input in (y : x : 1) format. in the darktable world we want to keep the338 // (x : y : 1) convention. therefore we need to flip coordinates first and339 // make sure that output is in correct format after corrections are applied.340341 const float u = width;342 const float v = height;343344 const float phi = deg2radf(angle);345 const float cosi = cosf(phi);346 const float sini = sinf(phi);347 const float ascale = sqrtf(aspect);348349 // most of this comes from ShiftN350 const float f_global = f_length_kb;351 const float horifac = 1.0f - orthocorr / 100.0f;352 const float exppa_v = expf(shift_v);353 const float fdb_v = f_global / (14.4f + (v / u - 1) * 7.2f);354 const float rad_v = fdb_v * (exppa_v - 1.0f) / (exppa_v + 1.0f);355 const float alpha_v = CLAMP(atanf(rad_v), -1.5f, 1.5f);356 const float rt_v = sinf(0.5f * alpha_v);357 const float r_v = fmaxf(0.1f, 2.0f * (horifac - 1.0f) * rt_v * rt_v + 1.0f);358359 const float vertifac = 1.0f - orthocorr / 100.0f;360 const float exppa_h = expf(shift_h);361 const float fdb_h = f_global / (14.4f + (u / v - 1) * 7.2f);362 const float rad_h = fdb_h * (exppa_h - 1.0f) / (exppa_h + 1.0f);363 const float alpha_h = CLAMP(atanf(rad_h), -1.5f, 1.5f);364 const float rt_h = sinf(0.5f * alpha_h);365 const float r_h = fmaxf(0.1f, 2.0f * (vertifac - 1.0f) * rt_h * rt_h + 1.0f);366367368 // three intermediate buffers for matrix calculation ...369 float DT_ALIGNED_ARRAY m1[3][3];370 float DT_ALIGNED_ARRAY m2[3][3];371 float DT_ALIGNED_ARRAY m3[3][3];372373 // ... and some pointers to handle them more intuitively374 float (*mwork)[3] = m1;375 float (*minput)[3] = m2;376 float (*moutput)[3] = m3;377378 // Step 1: flip x and y coordinates (see above)379 memset(minput, 0, sizeof(float) * 9);380 minput[0][1] = 1.0f;381 minput[1][0] = 1.0f;382 minput[2][2] = 1.0f;383384385 // Step 2: rotation of image around its center386 memset(mwork, 0, sizeof(float) * 9);387 mwork[0][0] = cosi;388 mwork[0][1] = -sini;389 mwork[1][0] = sini;390 mwork[1][1] = cosi;391 mwork[0][2] = -0.5f * v * cosi + 0.5f * u * sini + 0.5f * v;392 mwork[1][2] = -0.5f * v * sini - 0.5f * u * cosi + 0.5f * u;393 mwork[2][2] = 1.0f;394395 // multiply mwork * minput -> moutput396 mat3mul((float *)moutput, (float *)mwork, (float *)minput);397398399 // Step 3: apply shearing400 memset(mwork, 0, sizeof(float) * 9);401 mwork[0][0] = 1.0f;402 mwork[0][1] = shear;403 mwork[1][1] = 1.0f;404 mwork[1][0] = shear;405 mwork[2][2] = 1.0f;406407 // moutput (of last calculation) -> minput408 MAT3SWAP(minput, moutput);409 // multiply mwork * minput -> moutput410 mat3mul((float *)moutput, (float *)mwork, (float *)minput);411412413 // Step 4: apply vertical lens shift effect414 memset(mwork, 0, sizeof(float) * 9);415 mwork[0][0] = exppa_v;416 mwork[1][0] = 0.5f * ((exppa_v - 1.0f) * u) / v;417 mwork[1][1] = 2.0f * exppa_v / (exppa_v + 1.0f);418 mwork[1][2] = -0.5f * ((exppa_v - 1.0f) * u) / (exppa_v + 1.0f);419 mwork[2][0] = (exppa_v - 1.0f) / v;420 mwork[2][2] = 1.0f;421422 // moutput (of last calculation) -> minput423 MAT3SWAP(minput, moutput);424 // multiply mwork * minput -> moutput425 mat3mul((float *)moutput, (float *)mwork, (float *)minput);426427428 // Step 5: horizontal compression429 memset(mwork, 0, sizeof(float) * 9);430 mwork[0][0] = 1.0f;431 mwork[1][1] = r_v;432 mwork[1][2] = 0.5f * u * (1.0f - r_v);433 mwork[2][2] = 1.0f;434435 // moutput (of last calculation) -> minput436 MAT3SWAP(minput, moutput);437 // multiply mwork * minput -> moutput438 mat3mul((float *)moutput, (float *)mwork, (float *)minput);439440441 // Step 6: flip x and y back again442 memset(mwork, 0, sizeof(float) * 9);443 mwork[0][1] = 1.0f;444 mwork[1][0] = 1.0f;445 mwork[2][2] = 1.0f;446447 // moutput (of last calculation) -> minput448 MAT3SWAP(minput, moutput);449 // multiply mwork * minput -> moutput450 mat3mul((float *)moutput, (float *)mwork, (float *)minput);451452453 // from here output vectors would be in (x : y : 1) format454455 // Step 7: now we can apply horizontal lens shift with the same matrix format as above456 memset(mwork, 0, sizeof(float) * 9);457 mwork[0][0] = exppa_h;458 mwork[1][0] = 0.5f * ((exppa_h - 1.0f) * v) / u;459 mwork[1][1] = 2.0f * exppa_h / (exppa_h + 1.0f);460 mwork[1][2] = -0.5f * ((exppa_h - 1.0f) * v) / (exppa_h + 1.0f);461 mwork[2][0] = (exppa_h - 1.0f) / u;462 mwork[2][2] = 1.0f;463464 // moutput (of last calculation) -> minput465 MAT3SWAP(minput, moutput);466 // multiply mwork * minput -> moutput467 mat3mul((float *)moutput, (float *)mwork, (float *)minput);468469470 // Step 8: vertical compression471 memset(mwork, 0, sizeof(float) * 9);472 mwork[0][0] = 1.0f;473 mwork[1][1] = r_h;474 mwork[1][2] = 0.5f * v * (1.0f - r_h);475 mwork[2][2] = 1.0f;476477 // moutput (of last calculation) -> minput478 MAT3SWAP(minput, moutput);479 // multiply mwork * minput -> moutput480 mat3mul((float *)moutput, (float *)mwork, (float *)minput);481482483 // Step 9: apply aspect ratio scaling484 memset(mwork, 0, sizeof(float) * 9);485 mwork[0][0] = 1.0f * ascale;486 mwork[1][1] = 1.0f / ascale;487 mwork[2][2] = 1.0f;488489 // moutput (of last calculation) -> minput490 MAT3SWAP(minput, moutput);491 // multiply mwork * minput -> moutput492 mat3mul((float *)moutput, (float *)mwork, (float *)minput);493494495 // Step 10: find x/y offsets and apply according correction so that496 // no negative coordinates occur in output vector497 float umin = FLT_MAX, vmin = FLT_MAX;498 // visit all four corners499 for(int y = 0; y < height; y += height - 1)500 for(int x = 0; x < width; x += width - 1)501 {502 float pi[3], po[3];503 pi[0] = x;504 pi[1] = y;505 pi[2] = 1.0f;506 // moutput expects input in (x:y:1) format and gives output as (x:y:1)507 mat3mulv(po, (float *)moutput, pi);508 umin = MIN(umin, po[0] / po[2]);509 vmin = MIN(vmin, po[1] / po[2]);510 }511512 memset(mwork, 0, sizeof(float) * 9);513 mwork[0][0] = 1.0f;514 mwork[1][1] = 1.0f;515 mwork[2][2] = 1.0f;516 mwork[0][2] = -umin;517 mwork[1][2] = -vmin;518519 // moutput (of last calculation) -> minput520 MAT3SWAP(minput, moutput);521 // multiply mwork * minput -> moutput522 mat3mul((float *)moutput, (float *)mwork, (float *)minput);523524525 // on request we either keep the final matrix for forward conversions526 // or produce an inverted matrix for backward conversions527 if(dir == ASHIFT_HOMOGRAPH_FORWARD)528 {529 // we have what we need -> copy it to the right place530 memcpy(homograph, moutput, sizeof(float) * 9);531 }532 else533 {534 // generate inverted homograph (mat3inv function defined in colorspaces.c)535 if(mat3inv((float *)homograph, (float *)moutput))536 {537 // in case of error we set to unity matrix538 memset(mwork, 0, sizeof(float) * 9);539 mwork[0][0] = 1.0f;540 mwork[1][1] = 1.0f;541 mwork[2][2] = 1.0f;542 memcpy(homograph, mwork, sizeof(float) * 9);543 }544 }545}546#undef MAT3SWAP547548549550static void rgb2grey256(const float *const in,551 double *const out,552 const int width,553 const int height)554{555 const size_t npixels = (size_t)width * height;556557 DT_OMP_FOR()558 for(size_t index = 0; index < npixels; index++)559 {560 out[index] = (0.3f * in[4*index+0]561 + 0.59f * in[4*index+1]562 + 0.11f * in[4*index+2]) * 256.0;563 }564}565566567static gboolean line_detect(float *in,568 const int width,569 const int height,570 const int x_off,571 const int y_off,572 const float scale,573 dt_iop_ashift_line_t **alines,574 int *lcount,575 int *vcount,576 int *hcount,577 float *vweight,578 float *hweight,579 const dt_iop_ashift_enhance_t enhance,580 const gboolean is_raw)581{582 double *greyscale = NULL;583 double *lsd_lines = NULL;584 dt_iop_ashift_line_t *ashift_lines = NULL;585586 int vertical_count = 0;587 int horizontal_count = 0;588 float vertical_weight = 0.0f;589 float horizontal_weight = 0.0f;590591 // Browser input is already display-referred; no RAW/detail preprocessing.592 // allocate intermediate buffers593 greyscale = malloc(sizeof(double) * width * height);594 if(greyscale == NULL) goto error;595596 // convert to greyscale image597 rgb2grey256(in, greyscale, width, height);598599 // call the line segment detector LSD;600 // LSD stores the number of found lines in lines_count.601 // it returns structural details as vector 'double lines[7 * lines_count]'602 int lines_count;603604 lsd_lines = LineSegmentDetection(&lines_count, greyscale, width, height,605 LSD_SCALE, LSD_SIGMA_SCALE, LSD_QUANT,606 LSD_ANG_TH, LSD_LOG_EPS, LSD_DENSITY_TH,607 LSD_N_BINS, NULL, NULL, NULL);608609 // we count the lines that we really want to use610 int lct = 0;611 if(lines_count > 0)612 {613 // aggregate lines data into our own structures614 ashift_lines = malloc(sizeof(dt_iop_ashift_line_t) * lines_count);615616 if(ashift_lines == NULL) goto error;617618 for(int n = 0; n < lines_count; n++)619 {620 const float x1 = lsd_lines[n * 7 + 0];621 const float y1 = lsd_lines[n * 7 + 1];622 const float x2 = lsd_lines[n * 7 + 2];623 const float y2 = lsd_lines[n * 7 + 3];624625 // check for lines running along image borders and skip them.626 // these would likely be false-positives which could result627 // from any kind of processing artifacts628 if((fabsf(x1 - x2) < 1 && fmaxf(x1, x2) < 2)629 || (fabsf(x1 - x2) < 1 && fminf(x1, x2) > width - 3)630 || (fabsf(y1 - y2) < 1 && fmaxf(y1, y2) < 2)631 || (fabsf(y1 - y2) < 1 && fminf(y1, y2) > height - 3))632 continue;633634 // line position in absolute coordinates635 float px1 = x_off + x1;636 float py1 = y_off + y1;637 float px2 = x_off + x2;638 float py2 = y_off + y2;639640 // scale back to input buffer641 px1 /= scale;642 py1 /= scale;643 px2 /= scale;644 py2 /= scale;645646 // store as homogeneous coordinates647 ashift_lines[lct].p1[0] = px1;648 ashift_lines[lct].p1[1] = py1;649 ashift_lines[lct].p1[2] = 1.0f;650 ashift_lines[lct].p2[0] = px2;651 ashift_lines[lct].p2[1] = py2;652 ashift_lines[lct].p2[2] = 1.0f;653654 // calculate homogeneous coordinates of connecting line (defined by the two points)655 vec3prodn(ashift_lines[lct].L, ashift_lines[lct].p1, ashift_lines[lct].p2);656657 // normalaze line coordinates so that x^2 + y^2 = 1658 // (this will always succeed as L is a real line connecting two real points)659 vec3lnorm(ashift_lines[lct].L, ashift_lines[lct].L);660661 // length and width of rectangle (see LSD)662 ashift_lines[lct].length = sqrt((px2 - px1) * (px2 - px1)663 + (py2 - py1) * (py2 - py1));664 ashift_lines[lct].width = lsd_lines[n * 7 + 4] / scale;665666 // ... and weight (= length * width * angle precision)667 const float weight = ashift_lines[lct].length668 * ashift_lines[lct].width669 * lsd_lines[n * 7 + 5];670671 ashift_lines[lct].weight = weight;672673674 const float angle = rad2degf(atan2f(py2 - py1, px2 - px1));675 const gboolean vertical = fabsf(fabsf(angle) - 90.0f) < MAX_TANGENTIAL_DEVIATION;676 const gboolean horizontal = fabsf(fabsf(fabsf(angle) - 90.0f) - 90.0f) < MAX_TANGENTIAL_DEVIATION;677 const gboolean relevant = ashift_lines[lct].length > MIN_LINE_LENGTH;678679 // register type of line680 dt_iop_ashift_linetype_t type = ASHIFT_LINE_IRRELEVANT;681682 if(vertical && relevant)683 {684 type = ASHIFT_LINE_VERTICAL_SELECTED;685 vertical_count++;686 vertical_weight += weight;687 }688 else if(horizontal && relevant)689 {690 type = ASHIFT_LINE_HORIZONTAL_SELECTED;691 horizontal_count++;692 horizontal_weight += weight;693 }694 ashift_lines[lct].type = type;695696 // the next valid line697 lct++;698 }699 }700#ifdef ASHIFT_DEBUG701 printf("%d lines (vertical %d, horizontal %d, not relevant %d)\n",702 lines_count, vertical_count,703 horizontal_count, lct - vertical_count - horizontal_count);704 float xmin = FLT_MAX, xmax = FLT_MIN, ymin = FLT_MAX, ymax = FLT_MIN;705 for(int n = 0; n < lct; n++)706 {707 xmin = fmin(xmin, fmin(ashift_lines[n].p1[0], ashift_lines[n].p2[0]));708 xmax = fmax(xmax, fmax(ashift_lines[n].p1[0], ashift_lines[n].p2[0]));709 ymin = fmin(ymin, fmin(ashift_lines[n].p1[1], ashift_lines[n].p2[1]));710 ymax = fmax(ymax, fmax(ashift_lines[n].p1[1], ashift_lines[n].p2[1]));711 printf("x1 %.0f, y1 %.0f, x2 %.0f, y2 %.0f,"712 " length %.0f, width %f, X %f, Y %f, Z %f, type %d, scalars %f %f\n",713 ashift_lines[n].p1[0], ashift_lines[n].p1[1], ashift_lines[n].p2[0],714 ashift_lines[n].p2[1],715 ashift_lines[n].length, ashift_lines[n].width,716 ashift_lines[n].L[0], ashift_lines[n].L[1], ashift_lines[n].L[2],717 ashift_lines[n].type,718 vec3scalar(ashift_lines[n].p1, ashift_lines[n].L),719 vec3scalar(ashift_lines[n].p2, ashift_lines[n].L));720 }721 printf("xmin %.0f, xmax %.0f, ymin %.0f, ymax %.0f\n", xmin, xmax, ymin, ymax);722#endif723724 // store results in provided locations725 *lcount = lct;726 *vcount = vertical_count;727 *vweight = vertical_weight;728 *hcount = horizontal_count;729 *hweight = horizontal_weight;730 *alines = ashift_lines;731732 // free intermediate buffers733 free(lsd_lines);734 free(greyscale);735 return lct > 0;736737error:738 free(lsd_lines);739 free(greyscale);740 return FALSE;741}742743744static inline void swap(int *a, int *b)745{746 const int tmp = *a;747 *a = *b;748 *b = tmp;749}750751// do complete permutations752static gboolean quickperm(int *a, int *p, const int N, int *i)753{754 if(*i >= N) return FALSE;755756 p[*i]--;757 const int j = (*i % 2 == 1) ? p[*i] : 0;758 swap(&a[j], &a[*i]);759 *i = 1;760 while(p[*i] == 0)761 {762 p[*i] = *i;763 (*i)++;764 }765 return TRUE;766}767768// Fisher-Yates shuffle769static void shuffle(int *a, const int N)770{771 for(int i = 0; i < N; i++)772 {773 const int j = i + rand() % (N - i);774 swap(&a[j], &a[i]);775 }776}777778// factorial function779static int fact(const int n)780{781 return (n == 1 ? 1 : n * fact(n - 1));782}783784// We use a pseudo-RANSAC algorithm to elminiate ouliers from our set of lines. The785// original RANSAC works on linear optimization problems. Our model is nonlinear. We786// take advantage of the fact that lines interesting for our model are vantage lines787// that meet in one vantage point for each subset of lines (vertical/horizontal).788// Strategy: we construct a model by (random) sampling within the subset of lines and789// calculate the vantage point. Then we check the "distance" of all other lines to the790// vantage point. The model that gives highest number of lines combined with the highest791// total weight and lowest overall "distance" wins.792// Disadvantage: compared to the original RANSAC we don't get any model parameters that793// we could use for the following NMS fit.794//795// Self-tuning: we optimize "epsilon", the hurdle rate to reject a796// line as an outlier, by a number of dry runs first. The target797// average percentage value of lines to eliminate as outliers (without798// judging on the quality of the model) is given by799// RANSAC_ELIMINATION_RATIO,800//801// note: the actual percentage of outliers removed in the final run802// will be lower because we will finally look for the best quality803// model with the optimized epsilon and that quality value also804// encloses the number of good lines.805806static void ransac(const dt_iop_ashift_line_t *lines,807 int *index_set,808 int *inout_set,809 const int set_count,810 const float total_weight,811 const int xmin,812 const int xmax,813 const int ymin,814 const int ymax)815{816 if(set_count < 3) return;817818 const size_t set_size = set_count * sizeof(int);819 int *best_set = malloc(set_size);820 memcpy(best_set, index_set, set_size);821 int *best_inout = calloc(1, set_size);822823 float best_quality = 0.0f;824825 // hurdle value epsilon for rejecting a line as an outlier will be self-tuning826 // in a number of dry runs827 float epsilon = powf(10.0f, -RANSAC_EPSILON);828 float epsilon_step = RANSAC_EPSILON_STEP;829 // some accounting variables for self-tuning830 int lines_eliminated = 0;831 int valid_runs = 0;832833 // number of runs to optimize epsilon834 const int optiruns = RANSAC_OPTIMIZATION_STEPS * RANSAC_OPTIMIZATION_DRY_RUNS;835 // go for complete permutations on small set sizes, else for random sample consensus836 const int riter = (set_count > RANSAC_HURDLE) ? RANSAC_RUNS : fact(set_count);837838 // some data needed for quickperm839 int *perm = malloc(sizeof(int) * (set_count + 1));840 for(int n = 0; n < set_count + 1; n++) perm[n] = n;841 int piter = 1;842843 // inout holds good/bad qualification for each line844 int *inout = malloc(set_size);845846 for(int r = 0; r < optiruns + riter; r++)847 {848 // get random or systematic variation of index set849 if(set_count > RANSAC_HURDLE || r < optiruns)850 shuffle(index_set, set_count);851 else852 (void)quickperm(index_set, perm, set_count, &piter);853854 // summed quality evaluation of this run855 float quality = 0.0f;856857 // we build a model ouf of the first two lines858 const float *L1 = lines[index_set[0]].L;859 const float *L2 = lines[index_set[1]].L;860861 // get intersection point (ideally a vantage point)862 float DT_ALIGNED_PIXEL V[3];863 vec3prodn(V, L1, L2);864865 // catch special cases:866 // a) L1 and L2 are identical -> V is NULL -> no valid vantage point867 // b) vantage point lies inside image frame (no chance to correct for this case)868 if(vec3isnull(V) ||869 (fabsf(V[2]) > 0.0f &&870 V[0]/V[2] >= xmin &&871 V[1]/V[2] >= ymin &&872 V[0]/V[2] <= xmax &&873 V[1]/V[2] <= ymax))874 {875 // no valid model876 quality = 0.0f;877 }878 else879 {880 // valid model881882 // normalize V so that x^2 + y^2 + z^2 = 1883 vec3norm(V, V);884885 // the two lines constituting the model are part of the set886 inout[0] = 1;887 inout[1] = 1;888889 // go through all remaining lines, check if they are within the model, and890 // mark that fact in inout[].891 // summarize a quality parameter for all lines within the model892 for(int n = 2; n < set_count; n++)893 {894 // L is normalized so that x^2 + y^2 = 1895 const float *L3 = lines[index_set[n]].L;896897 // we take the absolute value of the dot product of V and L as898 // a measure of the "distance" between point and line. Note899 // that this is not the real euclidean distance but - with the900 // given normalization - just a pragmatically selected number901 // that goes to zero if V lies on L and increases the more V902 // and L are apart903 const float d = fabsf(vec3scalar(V, L3));904905 // depending on d we either include or exclude the point from the set906 inout[n] = (d < epsilon) ? 1 : 0;907908 float q;909910 if(inout[n] == 1)911 {912 // a quality parameter that depends 1/3 on the number of913 // lines within the model, 1/3 on their weight, and 1/3 on914 // their weighted distance d to the vantage point915 q = 0.33f / (float)set_count916 + 0.33f * lines[index_set[n]].weight / total_weight917 + 0.33f * (1.0f - d / epsilon)918 * (float)set_count * lines[index_set[n]].weight / total_weight;919 }920 else921 {922 q = 0.0f;923 lines_eliminated++;924 }925926 quality += q;927 }928 valid_runs++;929 }930931 if(r < optiruns)932 {933 // on last run of each self-tuning step934 if((r % RANSAC_OPTIMIZATION_DRY_RUNS) == (RANSAC_OPTIMIZATION_DRY_RUNS - 1)935 && (valid_runs > 0))936 {937#ifdef ASHIFT_DEBUG938 printf("ransac self-tuning (run %d): epsilon %f", r, epsilon);939#endif940 // average ratio of lines that we eliminated with the given epsilon941 float ratio = 100.0f * (float)lines_eliminated / ((float)set_count * valid_runs);942 // adjust epsilon accordingly943 if(ratio < RANSAC_ELIMINATION_RATIO)944 epsilon = powf(10.0f, log10(epsilon) - epsilon_step);945 else if(ratio > RANSAC_ELIMINATION_RATIO)946 epsilon = powf(10.0f, log10(epsilon) + epsilon_step);947#ifdef ASHIFT_DEBUG948 printf(" (elimination ratio %f) -> %f\n", ratio, epsilon);949#endif950 // reduce step-size for next optimization round951 epsilon_step /= 2.0f;952 lines_eliminated = 0;953 valid_runs = 0;954 }955 }956 else957 {958 // in the "real" runs check against the best model found so far959 if(quality > best_quality)960 {961 memcpy(best_set, index_set, set_size);962 memcpy(best_inout, inout, set_size);963 best_quality = quality;964 }965 }966967#ifdef ASHIFT_DEBUG968 // report some statistics969 int count = 0, lastcount = 0;970 for(int n = 0; n < set_count; n++) count += best_inout[n];971 for(int n = 0; n < set_count; n++) lastcount += inout[n];972 printf("ransac run %d: best qual %.6f, eps %.6f,"973 " line count %d of %d (this run: qual %.5f, count %d (%2f%%))\n", r,974 best_quality, epsilon, count, set_count, quality, lastcount,975 100.0f * lastcount / (float)set_count);976#endif977 }978979 // store back best set980 memcpy(index_set, best_set, set_size);981 memcpy(inout_set, best_inout, set_size);982983 free(inout);984 free(perm);985 free(best_inout);986 free(best_set);987}988989990// try to clean up structural data by eliminating outliers and thereby increasing991// the chance of a convergent fitting992static gboolean _remove_outliers(const dt_iop_module_t *self)993{994 dt_iop_ashift_gui_data_t *g = self->gui_data;995996 const int width = g->lines_in_width;997 const int height = g->lines_in_height;998 const int xmin = g->lines_x_off;999 const int ymin = g->lines_y_off;1000 const int xmax = xmin + width;1001 const int ymax = ymin + height;10021003 if(g->lines_count < 1)1004 return TRUE;1005 // holds the index set of lines we want to work on1006 int *lines_set = malloc(sizeof(int) * g->lines_count);1007 // holds the result of ransac1008 int *inout_set = malloc(sizeof(int) * g->lines_count);10091010 if(!lines_set || !inout_set)1011 goto error;1012 // some accounting variables1013 int vnb = 0, vcount = 0;1014 int hnb = 0, hcount = 0;10151016 // just to be on the safe side1017 if(g->lines == NULL) goto error;10181019 // generate index list for the vertical lines1020 for(int n = 0; n < g->lines_count; n++)1021 {1022 // is this a selected vertical line?1023 if((g->lines[n].type & ASHIFT_LINE_MASK) != ASHIFT_LINE_VERTICAL_SELECTED)1024 continue;10251026 lines_set[vnb] = n;1027 inout_set[vnb] = 1;1028 vnb++;1029 }10301031 // it only makes sense to call ransac if we have more than two lines1032 if(vnb > 2)1033 ransac(g->lines, lines_set, inout_set, vnb, g->vertical_weight,1034 xmin, xmax, ymin, ymax);10351036 // adjust line selected flag according to the ransac results1037 for(int n = 0; n < vnb; n++)1038 {1039 const int m = lines_set[n];1040 if(inout_set[n] == 1)1041 {1042 g->lines[m].type |= ASHIFT_LINE_SELECTED;1043 vcount++;1044 }1045 else1046 g->lines[m].type &= ~ASHIFT_LINE_SELECTED;1047 }1048 // update number of vertical lines1049 g->vertical_count = vcount;1050 g->lines_version++;10511052 // now generate index list for the horizontal lines1053 for(int n = 0; n < g->lines_count; n++)1054 {1055 // is this a selected horizontal line?1056 if((g->lines[n].type & ASHIFT_LINE_MASK) != ASHIFT_LINE_HORIZONTAL_SELECTED)1057 continue;10581059 lines_set[hnb] = n;1060 inout_set[hnb] = 1;1061 hnb++;1062 }10631064 // it only makes sense to call ransac if we have more than two lines1065 if(hnb > 2)1066 ransac(g->lines, lines_set, inout_set, hnb, g->horizontal_weight,1067 xmin, xmax, ymin, ymax);10681069 // adjust line selected flag according to the ransac results1070 for(int n = 0; n < hnb; n++)1071 {1072 const int m = lines_set[n];1073 if(inout_set[n] == 1)1074 {1075 g->lines[m].type |= ASHIFT_LINE_SELECTED;1076 hcount++;1077 }1078 else1079 g->lines[m].type &= ~ASHIFT_LINE_SELECTED;1080 }1081 // update number of horizontal lines1082 g->horizontal_count = hcount;1083 g->lines_version++;10841085 free(inout_set);1086 free(lines_set);10871088 return TRUE;10891090error:1091 free(inout_set);1092 free(lines_set);1093 return FALSE;1094}10951096// utility function to map a variable in [min; max] to [-INF; + INF]1097static inline double logit(const double x, const double min, const double max)1098{1099 const double eps = 1.0e-6;1100 // make sure p does not touch the borders of its definition area,1101 // not critical for data accuracy as logit() is only used on initial fit parameters1102 const double p = CLAMP((x - min) / (max - min), eps, 1.0 - eps);11031104 return (2.0 * atanh(2.0 * p - 1.0));1105}11061107// inverted function to logit()1108static inline double ilogit(const double L, const double min, const double max)1109{1110 const double p = 0.5 * (1.0 + tanh(0.5 * L));11111112 return (p * (max - min) + min);1113}11141115// helper function for simplex() return quality parameter for the given model1116// strategy:1117// * generate homography matrix out of fixed parameters and fitting parameters1118// * apply homography to all end points of affected lines1119// * generate new line out of transformed end points1120// * calculate scalar product s of line with perpendicular axis1121// * sum over weighted s^2 values11221123static double model_fitness(double *params, void *data)1124{1125 const dt_iop_ashift_fit_params_t *fit = (dt_iop_ashift_fit_params_t *)data;11261127 // just for convenience: get shorter names1128 const dt_iop_ashift_line_t *lines = fit->lines;1129 const int lines_count = fit->lines_count;1130 const int width = fit->width;1131 const int height = fit->height;1132 const float f_length_kb = fit->f_length_kb;1133 const float orthocorr = fit->orthocorr;1134 const float aspect = fit->aspect;11351136 float rotation = fit->rotation;1137 float lensshift_v = fit->lensshift_v;1138 float lensshift_h = fit->lensshift_h;1139 float shear = fit->shear;1140 const float rotation_range = fit->rotation_range;1141 const float lensshift_v_range = fit->lensshift_v_range;1142 const float lensshift_h_range = fit->lensshift_h_range;1143 const float shear_range = fit->shear_range;11441145 int pcount = 0;11461147 // fill in fit parameters from params[]. Attention: order matters!!!1148 if(dt_isnan(rotation))1149 {1150 rotation = ilogit(params[pcount], -rotation_range, rotation_range);1151 pcount++;1152 }11531154 if(dt_isnan(lensshift_v))1155 {1156 lensshift_v = ilogit(params[pcount], -lensshift_v_range, lensshift_v_range);1157 pcount++;1158 }11591160 if(dt_isnan(lensshift_h))1161 {1162 lensshift_h = ilogit(params[pcount], -lensshift_h_range, lensshift_h_range);1163 pcount++;1164 }11651166 if(dt_isnan(shear))1167 {1168 shear = ilogit(params[pcount], -shear_range, shear_range);1169 pcount++;1170 }11711172 assert(pcount == fit->params_count);11731174 // the possible reference axes1175 const float DT_ALIGNED_PIXEL Av[3] = { 1.0f, 0.0f, 0.0f };1176 const float DT_ALIGNED_PIXEL Ah[3] = { 0.0f, 1.0f, 0.0f };11771178 // generate homograph out of the parameters1179 float DT_ALIGNED_ARRAY homograph[3][3];1180 _homography((float *)homograph, rotation, lensshift_v, lensshift_h, shear, f_length_kb,1181 orthocorr, aspect, width, height, ASHIFT_HOMOGRAPH_FORWARD);11821183 // accounting variables1184 double sumsq_v = 0.0;1185 double sumsq_h = 0.0;1186 double weight_v = 0.0;1187 double weight_h = 0.0;1188 int count_v = 0;1189 int count_h = 0;1190 int count = 0;11911192 // iterate over all lines1193 for(int n = 0; n < lines_count; n++)1194 {1195 // check if this is a line which we must skip1196 if((lines[n].type & fit->linemask) != fit->linetype)1197 continue;11981199 // the direction of this line (vertical?)1200 const int isvertical = lines[n].type & ASHIFT_LINE_DIRVERT;12011202 // select the perpendicular reference axis1203 const float *A = isvertical ? Ah : Av;12041205 // apply homographic transformation to the end points1206 float DT_ALIGNED_PIXEL P1[3];1207 float DT_ALIGNED_PIXEL P2[3];1208 mat3mulv(P1, (float *)homograph, lines[n].p1);1209 mat3mulv(P2, (float *)homograph, lines[n].p2);12101211 // get line connecting the two points1212 float DT_ALIGNED_PIXEL L[3];1213 vec3prodn(L, P1, P2);12141215 // normalize L so that x^2 + y^2 = 1; makes sure that1216 // y^2 = 1 / (1 + m^2) and x^2 = m^2 / (1 + m^2) with m defining the slope of the line1217 vec3lnorm(L, L);12181219 // get scalar product of line L with orthogonal axis A -> gives 01220 // if line is perpendicular1221 const float s = vec3scalar(L, A);12221223 // sum up weighted s^2 for both directions individually1224 sumsq_v += isvertical ? s * s * lines[n].weight : 0.0;1225 weight_v += isvertical ? lines[n].weight : 0.0;1226 count_v += isvertical ? 1 : 0;1227 sumsq_h += !isvertical ? s * s * lines[n].weight : 0.0;1228 weight_h += !isvertical ? lines[n].weight : 0.0;1229 count_h += !isvertical ? 1 : 0;1230 count++;1231 }12321233 const double v = weight_v > 0.0f && count > 01234 ? sumsq_v / weight_v * (float)count_v / count1235 : 0.0;12361237 const double h = weight_h > 0.0f && count > 01238 ? sumsq_h / weight_h * (float)count_h / count1239 : 0.0;12401241 double sum = sqrt(1.0 - (1.0 - v) * (1.0 - h)) * 1.0e6;1242 //double sum = sqrt(v + h) * 1.0e6;12431244#ifdef ASHIFT_DEBUG1245 printf("fitness with rotation %f, lensshift_v %f,"1246 " lensshift_h %f, shear %f -> lines %d, quality %10f\n",1247 rotation, lensshift_v, lensshift_h, shear, count, sum);1248#endif12491250 return sum;1251}12521253// setup all data structures for fitting and call NM simplex1254static dt_iop_ashift_nmsresult_t nmsfit(dt_iop_module_t *self,1255 dt_iop_ashift_params_t *p,1256 const dt_iop_ashift_fitaxis_t dir)1257{1258 const dt_iop_ashift_gui_data_t *g = self->gui_data;12591260 if(!g->lines) return NMS_NOT_ENOUGH_LINES;1261 if(dir == ASHIFT_FIT_NONE) return NMS_SUCCESS;12621263 double params[4];1264 int pcount = 0;1265 int enough_lines = TRUE;12661267 // initialize fit parameters1268 dt_iop_ashift_fit_params_t fit;1269 fit.lines = g->lines;1270 fit.lines_count = g->lines_count;1271 fit.width = g->lines_in_width;1272 fit.height = g->lines_in_height;1273 fit.f_length_kb = (p->mode == ASHIFT_MODE_GENERIC)1274 ? DEFAULT_F_LENGTH1275 : p->f_length * p->crop_factor;1276 fit.orthocorr = (p->mode == ASHIFT_MODE_GENERIC) ? 0.0f : p->orthocorr;1277 fit.aspect = (p->mode == ASHIFT_MODE_GENERIC) ? 1.0f : p->aspect;1278 fit.rotation = p->rotation;1279 fit.lensshift_v = p->lensshift_v;1280 fit.lensshift_h = p->lensshift_h;1281 fit.shear = p->shear;1282 fit.rotation_range = g->rotation_range;1283 fit.lensshift_v_range = g->lensshift_v_range;1284 fit.lensshift_h_range = g->lensshift_h_range;1285 fit.shear_range = g->shear_range;1286 fit.linetype = ASHIFT_LINE_RELEVANT | ASHIFT_LINE_SELECTED;1287 fit.linemask = ASHIFT_LINE_MASK;1288 fit.params_count = 0;1289 fit.weight = 0.0f;12901291 // if the image is flipped and if we do not want to fit both lens shift1292 // directions or none at all, then we need to change direction1293 dt_iop_ashift_fitaxis_t mdir = dir;1294 if((mdir & ASHIFT_FIT_LENS_BOTH) != ASHIFT_FIT_LENS_BOTH &&1295 (mdir & ASHIFT_FIT_LENS_BOTH) != 0)1296 {1297 // flip all directions1298 mdir ^= g->isflipped ? ASHIFT_FIT_FLIP : 0;1299 // special case that needs to be corrected1300 mdir |= (mdir & ASHIFT_FIT_LINES_BOTH) == 0 ? ASHIFT_FIT_LINES_BOTH : 0;1301 }130213031304 // prepare fit structure and starting parameters for simplex fit.1305 // note: the sequence of parameters in params[] needs to match the1306 // respective order in dt_iop_ashift_fit_params_t. Parameters which are1307 // to be fittet are marked with NAN in the fit structure. Non-NAN1308 // parameters are assumed to be constant.1309 if(mdir & ASHIFT_FIT_ROTATION)1310 {1311 // we fit rotation1312 fit.params_count++;1313 params[pcount] = logit(fit.rotation, -fit.rotation_range, fit.rotation_range);1314 pcount++;1315 fit.rotation = NAN;1316 }13171318 if(mdir & ASHIFT_FIT_LENS_VERT)1319 {1320 // we fit vertical lens shift1321 fit.params_count++;1322 params[pcount] = logit(fit.lensshift_v, -fit.lensshift_v_range, fit.lensshift_v_range);1323 pcount++;1324 fit.lensshift_v = NAN;1325 }13261327 if(mdir & ASHIFT_FIT_LENS_HOR)1328 {1329 // we fit horizontal lens shift1330 fit.params_count++;1331 params[pcount] = logit(fit.lensshift_h, -fit.lensshift_h_range, fit.lensshift_h_range);1332 pcount++;1333 fit.lensshift_h = NAN;1334 }13351336 if(mdir & ASHIFT_FIT_SHEAR)1337 {1338 // we fit the shear parameter1339 fit.params_count++;1340 params[pcount] = logit(fit.shear, -fit.shear_range, fit.shear_range);1341 pcount++;1342 fit.shear = NAN;1343 }13441345 if(mdir & ASHIFT_FIT_LINES_VERT)1346 {1347 // we use vertical lines for fitting1348 fit.linetype |= ASHIFT_LINE_DIRVERT;1349 fit.weight += g->vertical_weight;1350 enough_lines = enough_lines && (g->vertical_count >= MINIMUM_FITLINES);1351 }13521353 if(mdir & ASHIFT_FIT_LINES_HOR)1354 {1355 // we use horizontal lines for fitting1356 fit.linetype |= 0;1357 fit.weight += g->horizontal_weight;1358 enough_lines = enough_lines && (g->horizontal_count >= MINIMUM_FITLINES);1359 }13601361 // this needs to come after ASHIFT_FIT_LINES_VERT and ASHIFT_FIT_LINES_HOR1362 if((mdir & ASHIFT_FIT_LINES_BOTH) == ASHIFT_FIT_LINES_BOTH)1363 {1364 // if we use fitting in both directions we need to1365 // adjust fit.linetype and fit.linemask to match all selected lines1366 fit.linetype = ASHIFT_LINE_RELEVANT | ASHIFT_LINE_SELECTED;1367 fit.linemask = ASHIFT_LINE_RELEVANT | ASHIFT_LINE_SELECTED;1368 }13691370 // error case: we do not run simplex if there are not enough lines1371 if(!enough_lines)1372 {1373#ifdef ASHIFT_DEBUG1374 printf("optimization not possible: insufficient number of lines\n");1375#endif1376 return NMS_NOT_ENOUGH_LINES;1377 }13781379 // start the simplex fit1380 const int iter = simplex(model_fitness, params, fit.params_count,1381 NMS_EPSILON, NMS_SCALE, NMS_ITERATIONS, NULL, (void*)&fit);13821383 // error case: the fit did not converge1384 if(iter >= NMS_ITERATIONS)1385 {1386#ifdef ASHIFT_DEBUG1387 printf("optimization not successful: maximum number of iterations reached (%d)\n", iter);1388#endif1389 return NMS_DID_NOT_CONVERGE;1390 }13911392 // fit was successful: now consolidate the results (order matters!!!)1393 pcount = 0;1394 fit.rotation = dt_isnan(fit.rotation)1395 ? ilogit(params[pcount++], -fit.rotation_range, fit.rotation_range)1396 : fit.rotation;13971398 fit.lensshift_v = dt_isnan(fit.lensshift_v)1399 ? ilogit(params[pcount++], -fit.lensshift_v_range, fit.lensshift_v_range)1400 : fit.lensshift_v;14011402 fit.lensshift_h = dt_isnan(fit.lensshift_h)1403 ? ilogit(params[pcount++], -fit.lensshift_h_range, fit.lensshift_h_range)1404 : fit.lensshift_h;14051406 fit.shear = dt_isnan(fit.shear)1407 ? ilogit(params[pcount++], -fit.shear_range, fit.shear_range)1408 : fit.shear;14091410#ifdef ASHIFT_DEBUG1411 printf("params after optimization (%d iterations): rotation %f,"1412 " lensshift_v %f, lensshift_h %f, shear %f\n",1413 iter, fit.rotation, fit.lensshift_v, fit.lensshift_h, fit.shear);1414#endif14151416 // sanity check: in case of extreme values the image gets distorted1417 // so strongly that it spans an insanely huge area. we check that1418 // case and assume values that increase the image area by more than1419 // a factor of 4 as being insane.1420 float DT_ALIGNED_ARRAY homograph[3][3];1421 _homography((float *)homograph, fit.rotation, fit.lensshift_v, fit.lensshift_h,1422 fit.shear, fit.f_length_kb,1423 fit.orthocorr, fit.aspect, fit.width, fit.height, ASHIFT_HOMOGRAPH_FORWARD);14241425 // visit all four corners and find maximum span1426 float xm = FLT_MAX, xM = -FLT_MAX, ym = FLT_MAX, yM = -FLT_MAX;1427 for(int y = 0; y < fit.height; y += fit.height - 1)1428 for(int x = 0; x < fit.width; x += fit.width - 1)1429 {1430 float DT_ALIGNED_PIXEL pi[3], DT_ALIGNED_PIXEL po[3];1431 pi[0] = x;1432 pi[1] = y;1433 pi[2] = 1.0f;1434 mat3mulv(po, (float *)homograph, pi);1435 po[0] /= po[2];1436 po[1] /= po[2];1437 xm = MIN(xm, po[0]);1438 ym = MIN(ym, po[1]);1439 xM = MAX(xM, po[0]);1440 yM = MAX(yM, po[1]);1441 }14421443 if((xM - xm) * (yM - ym) > 4.0f * fit.width * fit.height)1444 {1445#ifdef ASHIFT_DEBUG1446 printf("optimization not successful: degenerate case with"1447 " area growth factor (%f) exceeding limits\n",1448 (xM - xm) * (yM - ym) / (fit.width * fit.height));1449#endif1450 return NMS_INSANE;1451 }14521453 // now write the results into structure p1454 p->rotation = fit.rotation;1455 p->lensshift_v = fit.lensshift_v;1456 p->lensshift_h = fit.lensshift_h;1457 p->shear = fit.shear;1458 return NMS_SUCCESS;1459}1460