1 /*
2 This file is part of darktable,
3 Copyright (C) 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 #ifdef HAVE_CONFIG_H
19 #include "config.h"
20 #endif
21
22 #include "bauhaus/bauhaus.h"
23 #include "develop/imageop.h"
24 #include "develop/imageop_gui.h"
25 #include "gui/color_picker_proxy.h"
26 #include "gui/gtk.h"
27 #include "iop/iop_api.h"
28 #include "common/gaussian.h"
29 #include "common/fast_guided_filter.h"
30
31 #include <gtk/gtk.h>
32 #include <stdlib.h>
33
34 DT_MODULE_INTROSPECTION(1, dt_iop_cacorrectrgb_params_t)
35
36 /**
37 * Description of the approach
38 *
39 ** The problem
40 * chromatic aberration appear when:
41 * (1) channels are misaligned
42 * (2) or if some channel is more blurry than another.
43 *
44 * example case (1):
45 * _________
46 * _________| first channel
47 * _______
48 * ___________| second channel
49 * ^^ chromatic aberration
50 *
51 * other example case (1):
52 * _________
53 * _________| first channel
54 * ___________
55 * |_______ second channel
56 * ^^ chromatic aberration
57 *
58 * example case (2):
59 * _________
60 * | first channel
61 * _________|
62 * ________
63 * / second channel
64 * /
65 * ________/
66 * ^^^ chromatic aberration
67 *
68 * note that case (1) can already be partially corrected using the lens
69 * correction module.
70 *
71 ** Requirements for the solution
72 * - handle both cases
73 * - preserve borders as much as possible
74 * - be fast to compute
75 *
76 ** The solution
77 * The main idea is to represent 2 channels as a function of the third one.
78 *
79 * a very simple function is: guided = a * guide
80 * where a = blur(guided) / blur(guide)
81 * But this function is too simple to cope with borders.
82 *
83 * We stick with the idea of having guided channel as a factor of
84 * the guide channel, but instead of having a locally constant factor
85 * a, we use a factor that depends on the value of the guide pixel:
86 * guided = a(guide) * guide
87 *
88 * Our function a(guide) is pretty simple, it is a weighted average
89 * between 2 values (one high and one low), where the weights are
90 * dependent on the guide pixel value.
91 *
92 * Now, how do we determine these high and low value.
93 *
94 * We compute 2 manifolds.
95 * manifolds are partial local averages:
96 * some pixels are not used in the averages.
97 *
98 * for the lower manifold, we average only pixels whose guide value are below
99 * a local average of the guide.
100 * for the higher manifold, we average only pixels whose guide value are above
101 * a local average of the guide.
102 *
103 * for example here:
104 * _________
105 * _ _ _ _ _| _ _ _ _ _ _ _ average
106 * _________|
107 *
108 * ^^^^^^^^^ pixels below average (will be used to compute lower manifold)
109 *
110 * ^^^^^^^^^ pixels above average (will be used to compute higher manifold)
111 *
112 * As we want to write the guided channel as a ratio of the guide channel,
113 * we compute the manifolds on:
114 * - the guide channel
115 * - log difference between guide and guided
116 *
117 * using the log difference gives much better result than using directly the
118 * guided channel in the manifolds computation and computing the ratio after
119 * that, because averaging in linear makes lower manifolds harder to estimate
120 * accurately.
121 * Note that the repartition of pixels into higher and lower manifold
122 * computation is done by taking into account ONLY the guide channel.
123 *
124 * Once we have our 2 manifolds, with an average log difference for each of them
125 * (i.e. an average ratio), we can do a weighted mean to get the result.
126 * We weight more one ratio or the other depending to how close the guide pixel
127 * is from one manifold or another.
128 *
129 **/
130
131 typedef enum dt_iop_cacorrectrgb_guide_channel_t
132 {
133 DT_CACORRECT_RGB_R = 0, // $DESCRIPTION: "red"
134 DT_CACORRECT_RGB_G = 1, // $DESCRIPTION: "green"
135 DT_CACORRECT_RGB_B = 2 // $DESCRIPTION: "blue"
136 } dt_iop_cacorrectrgb_guide_channel_t;
137
138 typedef enum dt_iop_cacorrectrgb_mode_t
139 {
140 DT_CACORRECT_MODE_STANDARD = 0, // $DESCRIPTION: "standard"
141 DT_CACORRECT_MODE_DARKEN = 1, // $DESCRIPTION: "darken only"
142 DT_CACORRECT_MODE_BRIGHTEN = 2 // $DESCRIPTION: "brighten only"
143 } dt_iop_cacorrectrgb_mode_t;
144
145 typedef struct dt_iop_cacorrectrgb_params_t
146 {
147 dt_iop_cacorrectrgb_guide_channel_t guide_channel; // $DEFAULT: DT_CACORRECT_RGB_G $DESCRIPTION: "guide"
148 float radius; // $MIN: 1 $MAX: 500 $DEFAULT: 5 $DESCRIPTION: "radius"
149 float strength; // $MIN: 0 $MAX: 4 $DEFAULT: 0.5 $DESCRIPTION: "strength"
150 dt_iop_cacorrectrgb_mode_t mode; // $DEFAULT: DT_CACORRECT_MODE_STANDARD $DESCRIPTION: "correction mode"
151 gboolean refine_manifolds; // $MIN: FALSE $MAX: TRUE $DEFAULT: FALSE $DESCRIPTION: "very large chromatic aberration"
152 } dt_iop_cacorrectrgb_params_t;
153
154 typedef struct dt_iop_cacorrectrgb_gui_data_t
155 {
156 GtkWidget *guide_channel, *radius, *strength, *mode, *refine_manifolds;
157 } dt_iop_cacorrectrgb_gui_data_t;
158
name()159 const char *name()
160 {
161 return _("chromatic aberrations");
162 }
163
description(struct dt_iop_module_t * self)164 const char *description(struct dt_iop_module_t *self)
165 {
166 return dt_iop_set_description(self, _("correct chromatic aberrations"),
167 _("corrective"),
168 _("linear, raw, scene-referred"),
169 _("linear, raw"),
170 _("linear, raw, scene-referred"));
171 }
172
flags()173 int flags()
174 {
175 return IOP_FLAGS_INCLUDE_IN_STYLES | IOP_FLAGS_SUPPORTS_BLENDING;
176 }
177
default_group()178 int default_group()
179 {
180 return IOP_GROUP_CORRECT | IOP_GROUP_TECHNICAL;
181 }
182
default_colorspace(dt_iop_module_t * self,dt_dev_pixelpipe_t * pipe,dt_dev_pixelpipe_iop_t * piece)183 int default_colorspace(dt_iop_module_t *self, dt_dev_pixelpipe_t *pipe, dt_dev_pixelpipe_iop_t *piece)
184 {
185 return iop_cs_rgb;
186 }
187
commit_params(dt_iop_module_t * self,dt_iop_params_t * p1,dt_dev_pixelpipe_t * pipe,dt_dev_pixelpipe_iop_t * piece)188 void commit_params(dt_iop_module_t *self, dt_iop_params_t *p1, dt_dev_pixelpipe_t *pipe, dt_dev_pixelpipe_iop_t *piece)
189 {
190 memcpy(piece->data, p1, self->params_size);
191 }
192
normalize_manifolds(const float * const restrict blurred_in,float * const restrict blurred_manifold_lower,float * const restrict blurred_manifold_higher,const size_t width,const size_t height,const dt_iop_cacorrectrgb_guide_channel_t guide)193 static void normalize_manifolds(const float *const restrict blurred_in, float *const restrict blurred_manifold_lower, float *const restrict blurred_manifold_higher, const size_t width, const size_t height, const dt_iop_cacorrectrgb_guide_channel_t guide)
194 {
195 #ifdef _OPENMP
196 #pragma omp parallel for default(none) \
197 dt_omp_firstprivate(blurred_in, blurred_manifold_lower, blurred_manifold_higher, width, height, guide) \
198 schedule(simd:static)
199 #endif
200 for(size_t k = 0; k < width * height; k++)
201 {
202 const float weighth = fmaxf(blurred_manifold_higher[k * 4 + 3], 1E-2f);
203 const float weightl = fmaxf(blurred_manifold_lower[k * 4 + 3], 1E-2f);
204
205 // normalize guide
206 const float highg = blurred_manifold_higher[k * 4 + guide] / weighth;
207 const float lowg = blurred_manifold_lower[k * 4 + guide] / weightl;
208
209 blurred_manifold_higher[k * 4 + guide] = highg;
210 blurred_manifold_lower[k * 4 + guide] = lowg;
211
212 // normalize and unlog other channels
213 for(size_t kc = 0; kc <= 1; kc++)
214 {
215 const size_t c = (kc + guide + 1) % 3;
216 const float highc = blurred_manifold_higher[k * 4 + c] / weighth;
217 const float lowc = blurred_manifold_lower[k * 4 + c] / weightl;
218 blurred_manifold_higher[k * 4 + c] = exp2f(highc) * highg;
219 blurred_manifold_lower[k * 4 + c] = exp2f(lowc) * lowg;
220 }
221
222 // replace by average if weight is too small
223 if(weighth < 0.05f)
224 {
225 // we make a smooth transition between full manifold at
226 // weighth = 0.05f to full average at weighth = 0.01f
227 const float w = (weighth - 0.01f) / (0.05f - 0.01f);
228 for_each_channel(c,aligned(blurred_manifold_higher,blurred_in))
229 {
230 blurred_manifold_higher[k * 4 + c] = w * blurred_manifold_higher[k * 4 + c]
231 + (1.0f - w) * blurred_in[k * 4 + c];
232 }
233 }
234 if(weightl < 0.05f)
235 {
236 // we make a smooth transition between full manifold at
237 // weightl = 0.05f to full average at weightl = 0.01f
238 const float w = (weightl - 0.01f) / (0.05f - 0.01f);
239 for_each_channel(c,aligned(blurred_manifold_lower,blurred_in))
240 {
241 blurred_manifold_lower[k * 4 + c] = w * blurred_manifold_lower[k * 4 + c]
242 + (1.0f - w) * blurred_in[k * 4 + c];
243 }
244 }
245 }
246 }
247
248 #define DT_CACORRECTRGB_MAX_EV_DIFF 2.0f
get_manifolds(const float * const restrict in,const size_t width,const size_t height,const float sigma,const float sigma2,const dt_iop_cacorrectrgb_guide_channel_t guide,float * const restrict manifolds,gboolean refine_manifolds)249 static void get_manifolds(const float* const restrict in, const size_t width, const size_t height,
250 const float sigma, const float sigma2,
251 const dt_iop_cacorrectrgb_guide_channel_t guide,
252 float* const restrict manifolds, gboolean refine_manifolds)
253 {
254 float *const restrict blurred_in = dt_alloc_align_float(width * height * 4);
255 float *const restrict manifold_higher = dt_alloc_align_float(width * height * 4);
256 float *const restrict manifold_lower = dt_alloc_align_float(width * height * 4);
257 float *const restrict blurred_manifold_higher = dt_alloc_align_float(width * height * 4);
258 float *const restrict blurred_manifold_lower = dt_alloc_align_float(width * height * 4);
259
260 dt_aligned_pixel_t max = {INFINITY, INFINITY, INFINITY, INFINITY};
261 dt_aligned_pixel_t min = {-INFINITY, -INFINITY, -INFINITY, 0.0f};
262 // start with a larger blur to estimate the manifolds if we refine them
263 // later on
264 const float blur_size = refine_manifolds ? sigma2 : sigma;
265 dt_gaussian_t *g = dt_gaussian_init(width, height, 4, max, min, blur_size, 0);
266 if(!g) return;
267 dt_gaussian_blur_4c(g, in, blurred_in);
268
269 // construct the manifolds
270 // higher manifold is the blur of all pixels that are above average,
271 // lower manifold is the blur of all pixels that are below average
272 // we use the guide channel to categorize the pixels as above or below average
273 #ifdef _OPENMP
274 #pragma omp parallel for default(none) \
275 dt_omp_firstprivate(in, blurred_in, manifold_lower, manifold_higher, width, height, guide) \
276 schedule(simd:static)
277 #endif
278 for(size_t k = 0; k < width * height; k++)
279 {
280 const float pixelg = fmaxf(in[k * 4 + guide], 1E-6f);
281 const float avg = blurred_in[k * 4 + guide];
282 float weighth = (pixelg >= avg);
283 float weightl = (pixelg <= avg);
284 float logdiffs[2];
285 for(size_t kc = 0; kc <= 1; kc++)
286 {
287 const size_t c = (kc + guide + 1) % 3;
288 const float pixel = fmaxf(in[k * 4 + c], 1E-6f);
289 const float log_diff = log2f(pixel / pixelg);
290 logdiffs[kc] = log_diff;
291 }
292 // regularization of logdiff to avoid too many problems with noise:
293 // we lower the weights of pixels with too high logdiff
294 const float maxlogdiff = fmaxf(fabsf(logdiffs[0]), fabsf(logdiffs[1]));
295 if(maxlogdiff > DT_CACORRECTRGB_MAX_EV_DIFF)
296 {
297 const float correction_weight = DT_CACORRECTRGB_MAX_EV_DIFF / maxlogdiff;
298 weightl *= correction_weight;
299 weighth *= correction_weight;
300 }
301 for(size_t kc = 0; kc <= 1; kc++)
302 {
303 const size_t c = (kc + guide + 1) % 3;
304 manifold_higher[k * 4 + c] = logdiffs[kc] * weighth;
305 manifold_lower[k * 4 + c] = logdiffs[kc] * weightl;
306 }
307 manifold_higher[k * 4 + guide] = pixelg * weighth;
308 manifold_lower[k * 4 + guide] = pixelg * weightl;
309 manifold_higher[k * 4 + 3] = weighth;
310 manifold_lower[k * 4 + 3] = weightl;
311 }
312
313 dt_gaussian_blur_4c(g, manifold_higher, blurred_manifold_higher);
314 dt_gaussian_blur_4c(g, manifold_lower, blurred_manifold_lower);
315 dt_gaussian_free(g);
316
317 normalize_manifolds(blurred_in, blurred_manifold_lower, blurred_manifold_higher, width, height, guide);
318
319 // note that manifolds were constructed based on the value and average
320 // of the guide channel ONLY.
321 // this implies that the "higher" manifold in the channel c may be
322 // actually lower than the "lower" manifold of that channel.
323 // This happens in the following example:
324 // guide: 1_____
325 // |_____0
326 // guided: _____1
327 // 0_____|
328 // here the higher manifold of guide is equal to 1, its lower manifold is
329 // equal to 0. The higher manifold of the guided channel is equal to 0
330 // as it is the average of the values where the guide is higher than its
331 // average, and the lower manifold of the guided channel is equal to 1.
332
333 if(refine_manifolds)
334 {
335 g = dt_gaussian_init(width, height, 4, max, min, sigma, 0);
336 if(!g) return;
337 dt_gaussian_blur_4c(g, in, blurred_in);
338
339 // refine the manifolds
340 // improve result especially on very degraded images
341 // we use a blur of normal size for this step
342 #ifdef _OPENMP
343 #pragma omp parallel for default(none) \
344 dt_omp_firstprivate(in, blurred_in, manifold_lower, manifold_higher, blurred_manifold_lower, blurred_manifold_higher, width, height, guide) \
345 schedule(simd:static)
346 #endif
347 for(size_t k = 0; k < width * height; k++)
348 {
349 // in order to refine the manifolds, we will compute weights
350 // for which all channels will have a contribution.
351 // this will allow to avoid taking too much into account pixels
352 // that have wrong values due to the chromatic aberration
353 //
354 // for example, here:
355 // guide: 1_____
356 // |_____0
357 // guided: 1______
358 // |____0
359 // ^ this pixel makes the estimated lower manifold erroneous
360 // here, the higher and lower manifolds values computed are:
361 // _______|_higher_|________lower_________|
362 // guide | 1 | 0 |
363 // guided | 1 |(1 + 4 * 0) / 5 = 0.2 |
364 //
365 // the lower manifold of the guided is 0.2 if we consider only the guide
366 //
367 // at this step of the algorithm, we know estimates of manifolds
368 //
369 // we can refine the manifolds by computing weights that reduce the influence
370 // of pixels that are probably suffering from chromatic aberrations
371 const float pixelg = log2f(fmaxf(in[k * 4 + guide], 1E-6f));
372 const float highg = log2f(fmaxf(blurred_manifold_higher[k * 4 + guide], 1E-6f));
373 const float lowg = log2f(fmaxf(blurred_manifold_lower[k * 4 + guide], 1E-6f));
374 const float avgg = log2f(fmaxf(blurred_in[k * 4 + guide], 1E-6f));
375
376 float w = 1.0f;
377 for(size_t kc = 0; kc <= 1; kc++)
378 {
379 const size_t c = (guide + kc + 1) % 3;
380 // weight by considering how close pixel is for a manifold,
381 // and how close the log difference between the channels is
382 // close to the wrong log difference between the channels.
383
384 const float pixel = log2f(fmaxf(in[k * 4 + c], 1E-6f));
385 const float highc = log2f(fmaxf(blurred_manifold_higher[k * 4 + c], 1E-6f));
386 const float lowc = log2f(fmaxf(blurred_manifold_lower[k * 4 + c], 1E-6f));
387
388 // find how likely the pixel is part of a chromatic aberration
389 // (lowc, lowg) and (highc, highg) are valid points
390 // (lowc, highg) and (highc, lowg) are chromatic aberrations
391 const float dist_to_ll = fabsf(pixelg - lowg - pixel + lowc);
392 const float dist_to_hh = fabsf(pixelg - highg - pixel + highc);
393 const float dist_to_lh = fabsf((pixelg - pixel) - (highg - lowc));
394 const float dist_to_hl = fabsf((pixelg - pixel) - (lowg - highc));
395
396 float dist_to_good = 1.0f;
397 if(fabsf(pixelg - lowg) < fabsf(pixelg - highg))
398 dist_to_good = dist_to_ll;
399 else
400 dist_to_good = dist_to_hh;
401
402 float dist_to_bad = 1.0f;
403 if(fabsf(pixelg - lowg) < fabsf(pixelg - highg))
404 dist_to_bad = dist_to_hl;
405 else
406 dist_to_bad = dist_to_lh;
407
408 // make w higher if close to good, and smaller if close to bad.
409 w *= 1.0f * (0.2f + 1.0f / fmaxf(dist_to_good, 0.1f)) / (0.2f + 1.0f / fmaxf(dist_to_bad, 0.1f));
410 }
411
412 if(pixelg > avgg)
413 {
414 float logdiffs[2];
415 for(size_t kc = 0; kc <= 1; kc++)
416 {
417 const size_t c = (guide + kc + 1) % 3;
418 const float pixel = fmaxf(in[k * 4 + c], 1E-6f);
419 const float log_diff = log2f(pixel) - pixelg;
420 logdiffs[kc] = log_diff;
421 }
422 // regularization of logdiff to avoid too many problems with noise:
423 // we lower the weights of pixels with too high logdiff
424 const float maxlogdiff = fmaxf(fabsf(logdiffs[0]), fabsf(logdiffs[1]));
425 if(maxlogdiff > DT_CACORRECTRGB_MAX_EV_DIFF)
426 {
427 const float correction_weight = DT_CACORRECTRGB_MAX_EV_DIFF / maxlogdiff;
428 w *= correction_weight;
429 }
430 for(size_t kc = 0; kc <= 1; kc++)
431 {
432 const size_t c = (kc + guide + 1) % 3;
433 manifold_higher[k * 4 + c] = logdiffs[kc] * w;
434 }
435 manifold_higher[k * 4 + guide] = fmaxf(in[k * 4 + guide], 0.0f) * w;
436 manifold_higher[k * 4 + 3] = w;
437 // manifold_lower still contains the values from first iteration
438 // -> reset it.
439 for_four_channels(c)
440 {
441 manifold_lower[k * 4 + c] = 0.0f;
442 }
443 }
444 else
445 {
446 float logdiffs[2];
447 for(size_t kc = 0; kc <= 1; kc++)
448 {
449 const size_t c = (guide + kc + 1) % 3;
450 const float pixel = fmaxf(in[k * 4 + c], 1E-6f);
451 const float log_diff = log2f(pixel) - pixelg;
452 logdiffs[kc] = log_diff;
453 }
454 // regularization of logdiff to avoid too many problems with noise:
455 // we lower the weights of pixels with too high logdiff
456 const float maxlogdiff = fmaxf(fabsf(logdiffs[0]), fabsf(logdiffs[1]));
457 if(maxlogdiff > DT_CACORRECTRGB_MAX_EV_DIFF)
458 {
459 const float correction_weight = DT_CACORRECTRGB_MAX_EV_DIFF / maxlogdiff;
460 w *= correction_weight;
461 }
462 for(size_t kc = 0; kc <= 1; kc++)
463 {
464 const size_t c = (kc + guide + 1) % 3;
465 manifold_lower[k * 4 + c] = logdiffs[kc] * w;
466 }
467 manifold_lower[k * 4 + guide] = fmaxf(in[k * 4 + guide], 0.0f) * w;
468 manifold_lower[k * 4 + 3] = w;
469 // manifold_higher still contains the values from first iteration
470 // -> reset it.
471 for(size_t c = 0; c < 4; c++)
472 {
473 manifold_higher[k * 4 + c] = 0.0f;
474 }
475 }
476 }
477
478 dt_gaussian_blur_4c(g, manifold_higher, blurred_manifold_higher);
479 dt_gaussian_blur_4c(g, manifold_lower, blurred_manifold_lower);
480 normalize_manifolds(blurred_in, blurred_manifold_lower, blurred_manifold_higher, width, height, guide);
481 dt_gaussian_free(g);
482 }
483
484 dt_free_align(manifold_lower);
485 dt_free_align(manifold_higher);
486
487 // store all manifolds in the same structure to make upscaling faster
488 #ifdef _OPENMP
489 #pragma omp parallel for simd default(none) \
490 dt_omp_firstprivate(manifolds, blurred_manifold_lower, blurred_manifold_higher, width, height, guide) \
491 schedule(simd:static) aligned(manifolds, blurred_manifold_lower, blurred_manifold_higher:64)
492 #endif
493 for(size_t k = 0; k < width * height; k++)
494 {
495 for(size_t c = 0; c < 3; c++)
496 {
497 manifolds[k * 6 + c] = blurred_manifold_higher[k * 4 + c];
498 manifolds[k * 6 + 3 + c] = blurred_manifold_lower[k * 4 + c];
499 }
500 }
501 dt_free_align(blurred_in);
502 dt_free_align(blurred_manifold_lower);
503 dt_free_align(blurred_manifold_higher);
504 }
505 #undef DT_CACORRECTRGB_MAX_EV_DIFF
506
apply_correction(const float * const restrict in,const float * const restrict manifolds,const size_t width,const size_t height,const float sigma,const dt_iop_cacorrectrgb_guide_channel_t guide,const dt_iop_cacorrectrgb_mode_t mode,float * const restrict out)507 static void apply_correction(const float* const restrict in,
508 const float* const restrict manifolds,
509 const size_t width, const size_t height, const float sigma,
510 const dt_iop_cacorrectrgb_guide_channel_t guide,
511 const dt_iop_cacorrectrgb_mode_t mode,
512 float* const restrict out)
513
514 {
515 #ifdef _OPENMP
516 #pragma omp parallel for default(none) \
517 dt_omp_firstprivate(in, width, height, guide, manifolds, out, sigma, mode) \
518 schedule(simd:static)
519 #endif
520 for(size_t k = 0; k < width * height; k++)
521 {
522 const float high_guide = fmaxf(manifolds[k * 6 + guide], 1E-6f);
523 const float low_guide = fmaxf(manifolds[k * 6 + 3 + guide], 1E-6f);
524 const float log_high = log2f(high_guide);
525 const float log_low = log2f(low_guide);
526 const float dist_low_high = log_high - log_low;
527 const float pixelg = fmaxf(in[k * 4 + guide], 0.0f);
528 const float log_pixg = log2f(fminf(fmaxf(pixelg, low_guide), high_guide));
529
530 // determine how close our pixel is from the low manifold compared to the
531 // high manifold.
532 // if pixel value is lower or equal to the low manifold, weight_low = 1.0f
533 // if pixel value is higher or equal to the high manifold, weight_low = 0.0f
534 float weight_low = fabsf(log_high - log_pixg) / fmaxf(dist_low_high, 1E-6f);
535 // if the manifolds are very close, we are likely to introduce discontinuities
536 // and to have a meaningless "weight_low".
537 // thus in these cases make dist closer to 0.5.
538 // we set a threshold of 0.25f EV min.
539 const float threshold_dist_low_high = 0.25f;
540 if(dist_low_high < threshold_dist_low_high)
541 {
542 const float weight = dist_low_high / threshold_dist_low_high;
543 // dist_low_high = threshold_dist_low_high => dist
544 // dist_low_high = 0.0 => 0.5f
545 weight_low = weight_low * weight + 0.5f * (1.0f - weight);
546 }
547 const float weight_high = fmaxf(1.0f - weight_low, 0.0f);
548
549 for(size_t kc = 0; kc <= 1; kc++)
550 {
551 const size_t c = (guide + kc + 1) % 3;
552 const float pixelc = fmaxf(in[k * 4 + c], 0.0f);
553
554 const float ratio_high_manifolds = manifolds[k * 6 + c] / high_guide;
555 const float ratio_low_manifolds = manifolds[k * 6 + 3 + c] / low_guide;
556 // weighted geometric mean between the ratios.
557 const float ratio = powf(ratio_low_manifolds, weight_low) * powf(ratio_high_manifolds, weight_high);
558
559 const float outp = pixelg * ratio;
560
561 switch(mode)
562 {
563 case DT_CACORRECT_MODE_STANDARD:
564 out[k * 4 + c] = outp;
565 break;
566 case DT_CACORRECT_MODE_DARKEN:
567 out[k * 4 + c] = fminf(outp, pixelc);
568 break;
569 case DT_CACORRECT_MODE_BRIGHTEN:
570 out[k * 4 + c] = fmaxf(outp, pixelc);
571 break;
572 }
573 }
574
575 out[k * 4 + guide] = pixelg;
576 out[k * 4 + 3] = in[k * 4 + 3];
577 }
578 }
579
reduce_artifacts(const float * const restrict in,const size_t width,const size_t height,const float sigma,const dt_iop_cacorrectrgb_guide_channel_t guide,const float safety,float * const restrict out)580 static void reduce_artifacts(const float* const restrict in,
581 const size_t width, const size_t height, const float sigma,
582 const dt_iop_cacorrectrgb_guide_channel_t guide,
583 const float safety,
584 float* const restrict out)
585
586 {
587 // in_out contains the 2 guided channels of in, and the 2 guided channels of out
588 // it allows to blur all channels in one 4-channel gaussian blur instead of 2
589 float *const restrict DT_ALIGNED_PIXEL in_out = dt_alloc_align_float(width * height * 4);
590 #ifdef _OPENMP
591 #pragma omp parallel for default(none) \
592 dt_omp_firstprivate(in, out, in_out, width, height, guide) \
593 schedule(simd:static)
594 #endif
595 for(size_t k = 0; k < width * height; k++)
596 {
597 for(size_t kc = 0; kc <= 1; kc++)
598 {
599 const size_t c = (guide + kc + 1) % 3;
600 in_out[k * 4 + kc * 2 + 0] = in[k * 4 + c];
601 in_out[k * 4 + kc * 2 + 1] = out[k * 4 + c];
602 }
603 }
604
605 float *const restrict blurred_in_out = dt_alloc_align_float(width * height * 4);
606 dt_aligned_pixel_t max = {INFINITY, INFINITY, INFINITY, INFINITY};
607 dt_aligned_pixel_t min = {0.0f, 0.0f, 0.0f, 0.0f};
608 dt_gaussian_t *g = dt_gaussian_init(width, height, 4, max, min, sigma, 0);
609 if(!g) return;
610 dt_gaussian_blur_4c(g, in_out, blurred_in_out);
611 dt_gaussian_free(g);
612 dt_free_align(in_out);
613
614 // we consider that even with chromatic aberration, local average should
615 // be close to be accurate.
616 // thus, the local average of output should be similar to the one of the input
617 // if they are not, the algorithm probably washed out colors too much or
618 // may have produced artifacts.
619 // we do a weighted average between input and output, keeping more input if
620 // the local averages are very different.
621 // we use the same weight for all channels, as using different weights
622 // introduces artifacts in practice.
623 #ifdef _OPENMP
624 #pragma omp parallel for default(none) \
625 dt_omp_firstprivate(in, out, blurred_in_out, width, height, guide, safety) \
626 schedule(simd:static)
627 #endif
628 for(size_t k = 0; k < width * height; k++)
629 {
630 float w = 1.0f;
631 for(size_t kc = 0; kc <= 1; kc++)
632 {
633 const float avg_in = log2f(fmaxf(blurred_in_out[k * 4 + kc * 2 + 0], 1E-6f));
634 const float avg_out = log2f(fmaxf(blurred_in_out[k * 4 + kc * 2 + 1], 1E-6f));
635 w *= expf(-fmaxf(fabsf(avg_out - avg_in), 0.01f) * safety);
636 }
637 for(size_t kc = 0; kc <= 1; kc++)
638 {
639 const size_t c = (guide + kc + 1) % 3;
640 out[k * 4 + c] = fmaxf(1.0f - w, 0.0f) * fmaxf(in[k * 4 + c], 0.0f) + w * fmaxf(out[k * 4 + c], 0.0f);
641 }
642 }
643 dt_free_align(blurred_in_out);
644 }
645
reduce_chromatic_aberrations(const float * const restrict in,const size_t width,const size_t height,const size_t ch,const float sigma,const float sigma2,const dt_iop_cacorrectrgb_guide_channel_t guide,const dt_iop_cacorrectrgb_mode_t mode,const gboolean refine_manifolds,const float safety,float * const restrict out)646 static void reduce_chromatic_aberrations(const float* const restrict in,
647 const size_t width, const size_t height,
648 const size_t ch, const float sigma, const float sigma2,
649 const dt_iop_cacorrectrgb_guide_channel_t guide,
650 const dt_iop_cacorrectrgb_mode_t mode,
651 const gboolean refine_manifolds,
652 const float safety,
653 float* const restrict out)
654
655 {
656 const float downsize = fminf(3.0f, sigma);
657 const size_t ds_width = width / downsize;
658 const size_t ds_height = height / downsize;
659 float *const restrict ds_in = dt_alloc_align_float(ds_width * ds_height * 4);
660 // we use only one variable for both higher and lower manifolds in order
661 // to save time by doing only one bilinear interpolation instead of 2.
662 float *const restrict ds_manifolds = dt_alloc_align_float(ds_width * ds_height * 6);
663 // Downsample the image for speed-up
664 interpolate_bilinear(in, width, height, ds_in, ds_width, ds_height, 4);
665
666 // Compute manifolds
667 get_manifolds(ds_in, ds_width, ds_height, sigma / downsize, sigma2 / downsize, guide, ds_manifolds, refine_manifolds);
668 dt_free_align(ds_in);
669
670 // upscale manifolds
671 float *const restrict manifolds = dt_alloc_align_float(width * height * 6);
672 interpolate_bilinear(ds_manifolds, ds_width, ds_height, manifolds, width, height, 6);
673 dt_free_align(ds_manifolds);
674
675 apply_correction(in, manifolds, width, height, sigma, guide, mode, out);
676 dt_free_align(manifolds);
677
678 reduce_artifacts(in, width, height, sigma, guide, safety, out);
679 }
680
process(struct dt_iop_module_t * self,dt_dev_pixelpipe_iop_t * piece,const void * const ivoid,void * const ovoid,const dt_iop_roi_t * const roi_in,const dt_iop_roi_t * const roi_out)681 void process(struct dt_iop_module_t *self, dt_dev_pixelpipe_iop_t *piece, const void *const ivoid, void *const ovoid,
682 const dt_iop_roi_t *const roi_in, const dt_iop_roi_t *const roi_out)
683 {
684 if (!dt_iop_have_required_input_format(4 /*we need full-color pixels*/, self, piece->colors,
685 ivoid, ovoid, roi_in, roi_out))
686 return; // ivoid has been copied to ovoid and the module's trouble flag has been set
687
688 dt_iop_cacorrectrgb_params_t *d = (dt_iop_cacorrectrgb_params_t *)piece->data;
689 // used to adjuste blur level depending on size. Don't amplify noise if magnified > 100%
690 const float scale = fmaxf(piece->iscale / roi_in->scale, 1.f);
691 const int ch = piece->colors;
692 const size_t width = roi_out->width;
693 const size_t height = roi_out->height;
694 const float* in = (float*)ivoid;
695 float* out = (float*)ovoid;
696 const float sigma = fmaxf(d->radius / scale, 1.0f);
697 const float sigma2 = fmaxf(d->radius * d->radius / scale, 1.0f);
698
699 // whether to be very conservative in preserving the original image, or to
700 // keep algorithm result even if it overshoots
701 const float safety = powf(20.0f, 1.0f - d->strength);
702 reduce_chromatic_aberrations(in, width, height, ch, sigma, sigma2, d->guide_channel, d->mode, d->refine_manifolds, safety, out);
703 }
704
gui_update(dt_iop_module_t * self)705 void gui_update(dt_iop_module_t *self)
706 {
707 dt_iop_cacorrectrgb_gui_data_t *g = (dt_iop_cacorrectrgb_gui_data_t *)self->gui_data;
708 dt_iop_cacorrectrgb_params_t *p = (dt_iop_cacorrectrgb_params_t *)self->params;
709
710 dt_bauhaus_combobox_set_from_value(g->guide_channel, p->guide_channel);
711 dt_bauhaus_slider_set_soft(g->radius, p->radius);
712 dt_bauhaus_slider_set_soft(g->strength, p->strength);
713 dt_bauhaus_combobox_set_from_value(g->mode, p->mode);
714 gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(g->refine_manifolds), p->refine_manifolds);
715 }
716
reload_defaults(dt_iop_module_t * module)717 void reload_defaults(dt_iop_module_t *module)
718 {
719 dt_iop_cacorrectrgb_params_t *d = (dt_iop_cacorrectrgb_params_t *)module->default_params;
720
721 d->guide_channel = DT_CACORRECT_RGB_G;
722 d->radius = 5.0f;
723 d->strength = 0.5f;
724 d->mode = DT_CACORRECT_MODE_STANDARD;
725 d->refine_manifolds = FALSE;
726
727 dt_iop_cacorrectrgb_gui_data_t *g = (dt_iop_cacorrectrgb_gui_data_t *)module->gui_data;
728 if(g)
729 {
730 dt_bauhaus_combobox_set_default(g->guide_channel, d->guide_channel);
731 dt_bauhaus_slider_set_default(g->radius, d->radius);
732 dt_bauhaus_slider_set_soft_range(g->radius, 1.0, 20.0);
733 dt_bauhaus_slider_set_default(g->strength, d->strength);
734 dt_bauhaus_combobox_set_default(g->mode, d->mode);
735 gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(g->refine_manifolds), d->refine_manifolds);
736 }
737 }
738
gui_init(dt_iop_module_t * self)739 void gui_init(dt_iop_module_t *self)
740 {
741 dt_iop_cacorrectrgb_gui_data_t *g = IOP_GUI_ALLOC(cacorrectrgb);
742 self->widget = gtk_box_new(GTK_ORIENTATION_VERTICAL, DT_BAUHAUS_SPACE);
743 g->guide_channel = dt_bauhaus_combobox_from_params(self, "guide_channel");
744 gtk_widget_set_tooltip_text(g->guide_channel, _("channel used as a reference to\n"
745 "correct the other channels.\n"
746 "use sharpest channel if some\n"
747 "channels are blurry.\n"
748 "try changing guide channel if you\n"
749 "have artifacts."));
750 g->radius = dt_bauhaus_slider_from_params(self, "radius");
751 gtk_widget_set_tooltip_text(g->radius, _("increase for stronger correction"));
752 g->strength = dt_bauhaus_slider_from_params(self, "strength");
753 gtk_widget_set_tooltip_text(g->strength, _("balance between smoothing colors\n"
754 "and preserving them.\n"
755 "high values can lead to overshooting\n"
756 "and edge bleeding."));
757
758 gtk_box_pack_start(GTK_BOX(self->widget), dt_ui_section_label_new(_("advanced parameters")), TRUE, TRUE, 0);
759 g->mode = dt_bauhaus_combobox_from_params(self, "mode");
760 gtk_widget_set_tooltip_text(g->mode, _("correction mode to use.\n"
761 "can help with multiple\n"
762 "instances for very damaged\n"
763 "images.\n"
764 "darken only is particularly\n"
765 "efficient to correct blue\n"
766 "chromatic aberration."));
767 g->refine_manifolds = dt_bauhaus_toggle_from_params(self, "refine_manifolds");
768 gtk_widget_set_tooltip_text(g->refine_manifolds, _("runs an iterative approach\n"
769 "with several radii.\n"
770 "improves result on images\n"
771 "with very large chromatic\n"
772 "aberrations, but can smooth\n"
773 "colors too much on other\n"
774 "images."));
775 }
776