Ansel 0.0
A darktable fork - bloat + design vision
Loading...
Searching...
No Matches
diffuse.c
Go to the documentation of this file.
1/*
2 This file is part of the Ansel project.
3 Copyright (C) 2021-2023, 2025-2026 Aurélien PIERRE.
4 Copyright (C) 2021 Chris Elston.
5 Copyright (C) 2021 Hubert Kowalski.
6 Copyright (C) 2021 luzpaz.
7 Copyright (C) 2021-2022 Pascal Obry.
8 Copyright (C) 2021-2022 quovadit.
9 Copyright (C) 2021 Ralf Brown.
10 Copyright (C) 2021-2022 Sakari Kapanen.
11 Copyright (C) 2021 Victor Forsiuk.
12 Copyright (C) 2022 Diederik Ter Rahe.
13 Copyright (C) 2022 Hanno Schwalm.
14 Copyright (C) 2022 Martin Bařinka.
15 Copyright (C) 2022 Philipp Lutz.
16 Copyright (C) 2023, 2025 Guillaume Stutin.
17 Copyright (C) 2023 Luca Zulberti.
18 Copyright (C) 2024 Alynx Zhou.
19
20 Ansel is free software: you can redistribute it and/or modify
21 it under the terms of the GNU General Public License as published by
22 the Free Software Foundation, either version 3 of the License, or
23 (at your option) any later version.
24
25 Ansel is distributed in the hope that it will be useful,
26 but WITHOUT ANY WARRANTY; without even the implied warranty of
27 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
28 GNU General Public License for more details.
29
30 You should have received a copy of the GNU General Public License
31 along with Ansel. If not, see <http://www.gnu.org/licenses/>.
32*/
33
34#ifdef HAVE_CONFIG_H
35#include "config.h"
36#endif
37#include "widgets/bauhaus.h"
38#include "pixel/bspline.h"
39#include "system/macros.h"
40#include "system/mem_alloc.h"
42#include "common/logging.h"
43#include "system/openmp.h"
44#include "system/simd.h"
47#include "pixel/dwt.h"
48#include "develop/iop_profile.h"
49#include "common/opencl.h"
50#include "develop/develop.h"
51#include "develop/imageop_gui.h"
52#include "iop/noise_generator.h"
53#include "math/openmp_maths.h"
54#include "develop/tiling.h"
55
56#include "gui/presets.h"
57#include "iop/iop_api.h"
58#include "widgets/label.h"
59
60// Set to one to output intermediate image steps as PFM in /tmp
61#define DEBUG_DUMP_PFM 0
62
63// Diffuse v3 adds a new parameter that allows more aggressive "sharpening"
64// (and mathematically-accurate multiscale handling...)
65// on coarse scales, ensuring each HF details band is normalized to the same
66// energy. This makes the `radius span` parameter much more impactful.
67#define DIFFUSE_V3 0
68
69#if DIFFUSE_V3
71#else
73#endif
74
75#define MAX_NUM_SCALES 10
77{
78 // global parameters
79 int iterations; // $MIN: 0 $MAX: 500 $DEFAULT: 1 $DESCRIPTION: "iterations"
80 float sharpness; // $MIN: -1. $MAX: 1. $DEFAULT: 0. $DESCRIPTION: "sharpness"
81 int radius; // $MIN: 0 $MAX: 2048 $DEFAULT: 8 $DESCRIPTION: "radius span"
82 float regularization; // $MIN: 0. $MAX: 6. $DEFAULT: 0. $DESCRIPTION: "edge sensitivity"
83 float variance_threshold; // $MIN: -2. $MAX: 2. $DEFAULT: 0. $DESCRIPTION: "edge threshold"
84
85 float anisotropy_first; // $MIN: -10. $MAX: 10. $DEFAULT: 0. $DESCRIPTION: "1st order anisotropy"
86 float anisotropy_second; // $MIN: -10. $MAX: 10. $DEFAULT: 0. $DESCRIPTION: "2nd order anisotropy"
87 float anisotropy_third; // $MIN: -10. $MAX: 10. $DEFAULT: 0. $DESCRIPTION: "3rd order anisotropy"
88 float anisotropy_fourth; // $MIN: -10. $MAX: 10. $DEFAULT: 0. $DESCRIPTION: "4th order anisotropy"
89
90 float threshold; // $MIN: 0. $MAX: 8. $DEFAULT: 0. $DESCRIPTION: "luminance masking threshold"
91
92 float first; // $MIN: -1. $MAX: 1. $DEFAULT: 0. $DESCRIPTION: "1st order speed"
93 float second; // $MIN: -1. $MAX: 1. $DEFAULT: 0. $DESCRIPTION: "2nd order speed"
94 float third; // $MIN: -1. $MAX: 1. $DEFAULT: 0. $DESCRIPTION: "3rd order speed"
95 float fourth; // $MIN: -1. $MAX: 1. $DEFAULT: 0. $DESCRIPTION: "4th order speed"
96
97 // v2
98 int radius_center; // $MIN: 0 $MAX: 1024 $DEFAULT: 0 $DESCRIPTION: "central radius"
99
100 // new versions add params mandatorily at the end, so we can memcpy old parameters at the beginning
101
102 // v3 : Ansel 1.0
103 // bool normalize_band_energy; // $DEFAULT: FALSE $DESCRIPTION: "normalize coarse scales"
104 // When disabled, this will boost coarse scale sharpening by a lot.
105 // There is no reason to enable it for new edits,
106 // it's there to keep compatiblity with old edits.
107
109
110
116
129
130
131// only copy params struct to avoid a commit_params()
133
136{
137 memcpy(piece->data, params, self->params_size);
138 piece->cache_output_on_ram = TRUE;
139}
140
141
142typedef enum dt_isotropy_t
143{
144 DT_ISOTROPY_ISOTROPE = 0, // diffuse in all directions with same intensity
145 DT_ISOTROPY_ISOPHOTE = 1, // diffuse more in the isophote direction (orthogonal to gradient)
146 DT_ISOTROPY_GRADIENT = 2 // diffuse more in the gradient direction
148
149
151static inline dt_isotropy_t check_isotropy_mode(const float anisotropy)
152{
153 // user param is negative, positive or zero. The sign encodes the direction of diffusion, the magnitude encodes the ratio of anisotropy
154 // ultimately, the anisotropy factor needs to be positive before going into the exponential
155 if(anisotropy == 0.f)
157 else if(anisotropy > 0.f)
159 else
160 return DT_ISOTROPY_GRADIENT; // if(anisotropy > 0.f)
161}
162
163
164const char *name()
165{
166 return _("diffuse or _sharpen");
167}
168
169const char *aliases()
170{
171 return _("diffusion|deconvolution|blur|sharpening");
172}
173
174const char **description(struct dt_iop_module_t *self)
175{
176 return dt_iop_set_description(self,
177 _("simulate directional diffusion of light with heat transfer model\n"
178 "to apply an iterative edge-oriented blur,\n"
179 "inpaint damaged parts of the image,"
180 "or to remove blur with blind deconvolution."),
181 _("corrective and creative"),
182 _("linear, RGB, scene-referred"),
183 _("linear, RGB"),
184 _("linear, RGB, scene-referred"));
185}
186
188{
189 return IOP_GROUP_SHARPNESS;
190}
191
196
198{
199 return IOP_CS_RGB;
200}
201
202int legacy_params(dt_iop_module_t *self, const void *const old_params, const int old_version, void *new_params,
203 const int new_version)
204{
205 if(old_version == 1 && new_version == 2)
206 {
207 typedef struct dt_iop_diffuse_params_v1_t
208 {
209 // global parameters
210 int iterations;
211 float sharpness;
212 int radius;
213 float regularization;
214 float variance_threshold;
215
216 float anisotropy_first;
217 float anisotropy_second;
218 float anisotropy_third;
219 float anisotropy_fourth;
220
221 float threshold;
222
223 float first;
224 float second;
225 float third;
226 float fourth;
227 } dt_iop_diffuse_params_v1_t;
228
229 dt_iop_diffuse_params_v1_t *o = (dt_iop_diffuse_params_v1_t *)old_params;
232
233 *n = *d; // start with a fresh copy of default parameters
234
235 // copy common parameters
236 memcpy(n, o, sizeof(dt_iop_diffuse_params_v1_t));
237
238 // init only new parameters
239 n->radius_center = 0;
240
241#if !DIFFUSE_V3
242 // When version 3 will be out, we need to handle v1 -> v2 -> v3 conversion,
243 // so don't return just yet.
244 return 0;
245#endif
246 }
247
248#if DIFFUSE_V3
249 if(old_version == 2 && new_version == 3)
250 {
251 typedef struct dt_iop_diffuse_params_v2_t
252 {
253 // global parameters
254 int iterations; // $MIN: 0 $MAX: 500 $DEFAULT: 1 $DESCRIPTION: "iterations"
255 float sharpness; // $MIN: -1. $MAX: 1. $DEFAULT: 0. $DESCRIPTION: "sharpness"
256 int radius; // $MIN: 0 $MAX: 2048 $DEFAULT: 8 $DESCRIPTION: "radius span"
257 float regularization; // $MIN: 0. $MAX: 8. $DEFAULT: 0. $DESCRIPTION: "edge sensitivity"
258 float variance_threshold; // $MIN: -3. $MAX: 3. $DEFAULT: 0. $DESCRIPTION: "edge threshold"
259
260 float anisotropy_first; // $MIN: -100. $MAX: 100. $DEFAULT: 0. $DESCRIPTION: "1st order anisotropy"
261 float anisotropy_second; // $MIN: -100. $MAX: 100. $DEFAULT: 0. $DESCRIPTION: "2nd order anisotropy"
262 float anisotropy_third; // $MIN: -100. $MAX: 100. $DEFAULT: 0. $DESCRIPTION: "3rd order anisotropy"
263 float anisotropy_fourth; // $MIN: -100. $MAX: 100. $DEFAULT: 0. $DESCRIPTION: "4th order anisotropy"
264
265 float threshold; // $MIN: 0. $MAX: 8. $DEFAULT: 0. $DESCRIPTION: "luminance masking threshold"
266
267 float first; // $MIN: -1. $MAX: 1. $DEFAULT: 0. $DESCRIPTION: "1st order speed"
268 float second; // $MIN: -1. $MAX: 1. $DEFAULT: 0. $DESCRIPTION: "2nd order speed"
269 float third; // $MIN: -1. $MAX: 1. $DEFAULT: 0. $DESCRIPTION: "3rd order speed"
270 float fourth; // $MIN: -1. $MAX: 1. $DEFAULT: 0. $DESCRIPTION: "4th order speed"
271
272 // v2
273 int radius_center; // $MIN: 0 $MAX: 1024 $DEFAULT: 0 $DESCRIPTION: "central radius"
274
275 } dt_iop_diffuse_params_v2_t;
276
277 dt_iop_diffuse_params_v2_t *o = (dt_iop_diffuse_params_v2_t *)old_params;
280
281 *n = *d; // start with a fresh copy of default parameters
282
283 // copy common parameters
284 memcpy(n, o, sizeof(dt_iop_diffuse_params_v2_t));
285
286 // init only new parameters
287 n->normalize_band_energy = 1; // legacy compatiblity
288
289 return 0;
290 }
291#endif
292
293 return 1;
294}
295
297{
299 memset(&p, 0, sizeof(p));
300 p.radius_center = 0;
301
302 // deblurring presets
303 p.sharpness = 0.0f;
304 p.threshold = 0.0f;
305 p.variance_threshold = +0.f;
306 p.regularization = 1.f;
307
308 p.anisotropy_first = +2.f;
309 p.anisotropy_second = 0.f;
310 p.anisotropy_third = +2.f;
311 p.anisotropy_fourth = 0.f;
312
313 p.first = -0.25f;
314 p.second = +0.125f;
315 p.third = -0.125f;
316 p.fourth = +0.0625f;
317
318 p.radius = 8;
319 p.iterations = 8;
320 dt_gui_presets_add_generic(_("lens deblur: soft"), self->op, self->version(), &p, sizeof(p), 1,
322
323 p.radius = 10;
324 p.iterations = 16;
325 dt_gui_presets_add_generic(_("lens deblur: medium"), self->op, self->version(), &p, sizeof(p), 1,
327
328 p.radius = 12;
329 p.iterations = 24;
330 dt_gui_presets_add_generic(_("lens deblur: hard"), self->op, self->version(), &p, sizeof(p), 1,
332
333 p.iterations = 10;
334 p.radius = 512;
335 p.sharpness = 0.f;
336 p.variance_threshold = 0.f;
337 p.regularization = 2.5f;
338
339 p.first = -0.20f;
340 p.second = +0.10f;
341 p.third = -0.20f;
342 p.fourth = +0.10f;
343
344 p.anisotropy_first = 2.f;
345 p.anisotropy_second = 0.f;
346 p.anisotropy_third = 2.f;
347 p.anisotropy_fourth = 0.f;
348
349 dt_gui_presets_add_generic(_("dehaze"), self->op, self->version(), &p, sizeof(p), 1,
351
352 /* Denoise presets, retuned for issue #879 (they were too faint to see).
353 * Values are the RMSE optimum of a two-stage parametric sweep over
354 * synthetic Poisson-Gaussian noise, calibrated per ISO bucket from the
355 * median of 216 modern-camera profiles in noiseprofiles.json and validated
356 * across three pictures: fine ISO 400-800, medium 800-1600, coarse
357 * 1600-3200. The faintness had two causes: speeds far below useful, and
358 * edge sensitivity 4.0 throttling diffusion image-wide (2.5 wins at every
359 * level). The optima all sit at 0.20-0.25 of the 1.0 band-speed stability
360 * budget, and RMSE stays flat up to 1.0 — no degeneration nearby. Fine and
361 * medium share speeds on purpose: their differentiation is the radius,
362 * which is the physically right axis for grain size; coarse grain wants
363 * pure 3rd-order diffusion. */
364 p.iterations = 32;
365 p.sharpness = 0.f;
366 p.threshold = 0.f;
367 p.variance_threshold = -0.f;
368 p.regularization = 2.5f;
369
370 p.anisotropy_first = +2.f;
371 p.anisotropy_second = 0.f;
372 p.anisotropy_third = +2.f;
373 p.anisotropy_fourth = 0.f;
374
375 p.radius = 1;
376 p.radius_center = 2;
377
378 p.first = +0.10f;
379 p.second = 0.f;
380 p.third = +0.10f;
381 p.fourth = 0.f;
382 dt_gui_presets_add_generic(_("denoise: fine"), self->op, self->version(), &p, sizeof(p), 1, DEVELOP_BLEND_CS_RGB_SCENE);
383
384 p.radius = 3;
385 p.radius_center = 4;
386
387 p.first = +0.10f;
388 p.second = 0.f;
389 p.third = +0.10f;
390 p.fourth = 0.f;
391 dt_gui_presets_add_generic(_("denoise: medium"), self->op, self->version(), &p, sizeof(p), 1, DEVELOP_BLEND_CS_RGB_SCENE);
392
393 p.radius = 6;
394 p.radius_center = 8;
395
396 p.first = 0.f;
397 p.second = 0.f;
398 p.third = +0.25f;
399 p.fourth = 0.f;
400 dt_gui_presets_add_generic(_("denoise: coarse"), self->op, self->version(), &p, sizeof(p), 1, DEVELOP_BLEND_CS_RGB_SCENE);
401
402 p.radius_center = 0;
403
404 p.iterations = 2;
405 p.radius = 32;
406 p.sharpness = 0.0f;
407 p.threshold = 0.0f;
408 p.variance_threshold = 0.f;
409 p.regularization = 4.f;
410
411 p.anisotropy_first = +4.f;
412 p.anisotropy_second = +4.f;
413 p.anisotropy_third = +4.f;
414 p.anisotropy_fourth = +4.f;
415
416 p.first = +1.f;
417 p.second = +1.f;
418 p.third = +1.f;
419 p.fourth = +1.f;
420 dt_gui_presets_add_generic(_("surface blur"), self->op, self->version(), &p, sizeof(p), 1, DEVELOP_BLEND_CS_RGB_SCENE);
421
422 p.iterations = 1;
423 p.radius = 32;
424 p.sharpness = 0.0f;
425 p.threshold = 0.0f;
426 p.variance_threshold = 0.f;
427 p.regularization = 0.f;
428
429 p.anisotropy_first = 0.f;
430 p.anisotropy_second = 0.f;
431 p.anisotropy_third = 0.f;
432 p.anisotropy_fourth = 0.f;
433
434 p.first = +0.5f;
435 p.second = +0.5f;
436 p.third = +0.5f;
437 p.fourth = +0.5f;
438 dt_gui_presets_add_generic(_("bloom"), self->op, self->version(), &p, sizeof(p), 1, DEVELOP_BLEND_CS_RGB_SCENE);
439
440 p.iterations = 1;
441 p.radius = 4;
442 p.sharpness = 0.0f;
443 p.threshold = 0.0f;
444 p.variance_threshold = 0.f;
445 p.regularization = 1.f;
446
447 p.anisotropy_first = +1.f;
448 p.anisotropy_second = +1.f;
449 p.anisotropy_third = +1.f;
450 p.anisotropy_fourth = +1.f;
451
452 p.first = -0.25f;
453 p.second = -0.25f;
454 p.third = -0.25f;
455 p.fourth = -0.25f;
456 dt_gui_presets_add_generic(_("sharpen demosaicing (no AA filter)"), self->op, self->version(), &p, sizeof(p), 1,
458
459 p.radius = 8;
460 dt_gui_presets_add_generic(_("sharpen demosaicing (AA filter)"), self->op, self->version(), &p, sizeof(p), 1,
462
463 p.iterations = 4;
464 p.radius = 64;
465 p.sharpness = 0.0f;
466 p.threshold = 0.0f;
467 p.variance_threshold = 0.f;
468 p.regularization = 2.f;
469
470 p.anisotropy_first = 0.f;
471 p.anisotropy_second = 0.f;
472 p.anisotropy_third = +4.f;
473 p.anisotropy_fourth = +4.f;
474
475 p.first = 0.f;
476 p.second = 0.f;
477 p.third = +0.5f;
478 p.fourth = +0.5f;
479 dt_gui_presets_add_generic(_("simulate watercolor"), self->op, self->version(), &p, sizeof(p), 1, DEVELOP_BLEND_CS_RGB_SCENE);
480
481 p.iterations = 50;
482 p.radius = 64;
483 p.sharpness = 0.0f;
484 p.threshold = 0.0f;
485 p.variance_threshold = 0.f;
486 p.regularization = 4.f;
487
488 p.anisotropy_first = -5.f;
489 p.anisotropy_second = -5.f;
490 p.anisotropy_third = -5.f;
491 p.anisotropy_fourth = -5.f;
492
493 p.first = -1.f;
494 p.second = -1.f;
495 p.third = -1.f;
496 p.fourth = -1.f;
497 dt_gui_presets_add_generic(_("simulate line drawing"), self->op, self->version(), &p, sizeof(p), 1, DEVELOP_BLEND_CS_RGB_SCENE);
498
499 // local contrast
500 p.sharpness = 0.0f;
501 p.threshold = 0.0f;
502 p.variance_threshold = 0.f;
503
504 p.anisotropy_first = -2.5f;
505 p.anisotropy_second = 0.f;
506 p.anisotropy_third = 0.f;
507 p.anisotropy_fourth = -2.5f;
508
509 p.first = -0.50f;
510 p.second = 0.f;
511 p.third = 0.f;
512 p.fourth = -0.50f;
513
514 p.iterations = 10;
515 p.radius = 333;
516 p.radius_center = 512;
517 p.regularization = 0.1f;
518 dt_gui_presets_add_generic(_("add local contrast"), self->op, self->version(), &p, sizeof(p), 1,
520
521 p.iterations = 32;
522 p.radius = 4;
523 p.radius_center = 0;
524 p.sharpness = 0.0f;
525 p.threshold = 1.41f;
526 p.variance_threshold = 0.f;
527 p.regularization = 0.f;
528
529 p.anisotropy_first = +0.f;
530 p.anisotropy_second = +0.f;
531 p.anisotropy_third = +0.f;
532 p.anisotropy_fourth = +2.f;
533
534 p.first = +0.0f;
535 p.second = +0.0f;
536 p.third = +0.0f;
537 p.fourth = +0.5f;
538 dt_gui_presets_add_generic(_("inpaint highlights"), self->op, self->version(), &p, sizeof(p), 1, DEVELOP_BLEND_CS_RGB_SCENE);
539
540 // fast presets for slow hardware
541 p.radius_center = 0;
542 p.radius = 128;
543 p.sharpness = 0.0f;
544 p.threshold = 0.0f;
545 p.variance_threshold = 0.f;
546 p.regularization = 0.f;
547
548 p.anisotropy_first = 0.f;
549 p.anisotropy_second = 0.f;
550 p.anisotropy_third = 5.f;
551 p.anisotropy_fourth = 0.f;
552
553 p.first = 0.f;
554 p.second = 0.f;
555 p.third = -0.50f;
556 p.fourth = 0.f;
557
558 p.iterations = 1;
559 dt_gui_presets_add_generic(_("fast sharpness"), self->op, self->version(), &p, sizeof(p), 1,
561
562 p.radius_center = 512;
563 p.radius = 512;
564 p.sharpness = 0.0f;
565 p.threshold = 0.0f;
566 p.variance_threshold = 0.f;
567 p.regularization = 0.f;
568
569
570 p.anisotropy_first = 0.f;
571 p.anisotropy_second = 0.f;
572 p.anisotropy_third = 5.f;
573 p.anisotropy_fourth = 0.f;
574
575 p.first = 0.f;
576 p.second = 0.f;
577 p.third = -0.50f;
578 p.fourth = 0.f;
579
580 p.iterations = 1;
581 dt_gui_presets_add_generic(_("fast local contrast"), self->op, self->version(), &p, sizeof(p), 1,
583}
584
585void tiling_callback(struct dt_iop_module_t *self, const struct dt_dev_pixelpipe_t *pipe, const struct dt_dev_pixelpipe_iop_t *piece, struct dt_develop_tiling_t *tiling)
586{
587 const dt_iop_roi_t *const roi_in = &piece->roi_in;
589
590 const float zoom = dt_dev_get_module_scale(pipe, roi_in);
591 const float final_radius = (data->radius + data->radius_center) * 2.f / zoom;
592 const int diffusion_scales = num_steps_to_reach_equivalent_sigma(B_SPLINE_SIGMA, final_radius);
593 const int scales = CLAMP(diffusion_scales, 1, MAX_NUM_SCALES);
594 const int max_filter_radius = (1 << scales);
595
596 // Account for the exact full-frame buffers kept alive by the CPU/OpenCL paths:
597 // one borrowed input, one output, two temp ping-pong buffers, two low-pass ping-pong
598 // buffers, one stored detail buffer per wavelet scale, and one 8-bit mask.
599 tiling->factor = 6.0625f + scales;
600 tiling->factor_cl = 6.0625f + scales;
601
602 tiling->maxbuf = 1.0f;
603 tiling->maxbuf_cl = 1.0f;
604 tiling->overhead = 0;
605 tiling->overlap = max_filter_radius;
606 tiling->xalign = 1;
607 tiling->yalign = 1;
608 return;
609}
610
612static inline void init_reconstruct(float *const restrict reconstructed, const size_t width,
613 const size_t height)
614{
615// init the reconstructed buffer with non-clipped and partially clipped pixels
616 __OMP_PARALLEL_FOR_SIMD__(aligned(reconstructed:64))
617 for(size_t k = 0; k < height * width * 4; k++) reconstructed[k] = 0.f;
618
619}
620
621
622// Discretization parameters for the Partial Derivative Equation solver
623#define H 1 // spatial step
624#define KAPPA 0.25f // 0.25 if h = 1, 1 if h = 2
625
626
627static inline __attribute__((always_inline)) void find_gradients(const dt_aligned_pixel_simd_t pixels[9],
628 dt_aligned_pixel_simd_t xy[2])
629{
630 // Compute the gradient with centered finite differences in a 3x3 stencil
631 // warning : x is vertical, y is horizontal
632 const dt_aligned_pixel_simd_t half = dt_simd_set1(0.5f);
633 xy[0] = (pixels[7] - pixels[1]) * half;
634 xy[1] = (pixels[5] - pixels[3]) * half;
635}
636
637static inline __attribute__((always_inline)) void find_laplacians(const dt_aligned_pixel_simd_t pixels[9],
638 dt_aligned_pixel_simd_t xy[2])
639{
640 // Compute the laplacian with centered finite differences in a 3x3 stencil
641 // warning : x is vertical, y is horizontal
642 const dt_aligned_pixel_simd_t two = dt_simd_set1(2.f);
643 xy[0] = (pixels[7] + pixels[1]) - two * pixels[4];
644 xy[1] = (pixels[5] + pixels[3]) - two * pixels[4];
645}
646
647
648static inline __attribute__((always_inline)) void rotation_matrix_isophote(
649 const dt_aligned_pixel_simd_t c2, const dt_aligned_pixel_simd_t cos_theta_sin_theta,
650 const dt_aligned_pixel_simd_t cos_theta2, const dt_aligned_pixel_simd_t sin_theta2,
651 dt_aligned_pixel_simd_t a[2][2])
652{
653 // Write the coefficients of a square symmetrical matrice of rotation of the gradient :
654 // [[ a11, a12 ],
655 // [ a12, a22 ]]
656 // taken from https://www.researchgate.net/publication/220663968
657 // c dampens the gradient direction
658 a[0][0] = cos_theta2 + c2 * sin_theta2;
659 a[1][1] = c2 * cos_theta2 + sin_theta2;
660 a[0][1] = a[1][0] = (c2 - dt_simd_set1(1.f)) * cos_theta_sin_theta;
661}
662
663static inline __attribute__((always_inline)) void rotation_matrix_gradient(
664 const dt_aligned_pixel_simd_t c2, const dt_aligned_pixel_simd_t cos_theta_sin_theta,
665 const dt_aligned_pixel_simd_t cos_theta2, const dt_aligned_pixel_simd_t sin_theta2,
666 dt_aligned_pixel_simd_t a[2][2])
667{
668 // Write the coefficients of a square symmetrical matrice of rotation of the gradient :
669 // [[ a11, a12 ],
670 // [ a12, a22 ]]
671 // based on https://www.researchgate.net/publication/220663968 and inverted
672 // c dampens the isophote direction
673 a[0][0] = c2 * cos_theta2 + sin_theta2;
674 a[1][1] = cos_theta2 + c2 * sin_theta2;
675 a[0][1] = a[1][0] = (dt_simd_set1(1.f) - c2) * cos_theta_sin_theta;
676}
677
678
679static inline __attribute__((always_inline)) void build_matrix(const dt_aligned_pixel_simd_t a[2][2],
680 dt_aligned_pixel_simd_t kernel[9])
681{
682 const dt_aligned_pixel_simd_t half = dt_simd_set1(0.5f);
683 const dt_aligned_pixel_simd_t minus_two = dt_simd_set1(-2.f);
684 const dt_aligned_pixel_simd_t b11 = a[0][1] * half;
685 const dt_aligned_pixel_simd_t b13 = -b11;
686 const dt_aligned_pixel_simd_t b22 = minus_two * (a[0][0] + a[1][1]);
687
688 // build the kernel of rotated anisotropic laplacian
689 // from https://www.researchgate.net/publication/220663968 :
690 // [ [ a12 / 2, a22, -a12 / 2 ],
691 // [ a11, -2 (a11 + a22), a11 ],
692 // [ -a12 / 2, a22, a12 / 2 ] ]
693 // N.B. we have flipped the signs of the a12 terms
694 // compared to the paper. There's probably a mismatch
695 // of coordinate convention between the paper and the
696 // original derivation of this convolution mask
697 // (Witkin 1991, https://doi.org/10.1145/127719.122750).
698 kernel[0] = b11;
699 kernel[1] = a[1][1];
700 kernel[2] = b13;
701 kernel[3] = a[0][0];
702 kernel[4] = b22;
703 kernel[5] = a[0][0];
704 kernel[6] = b13;
705 kernel[7] = a[1][1];
706 kernel[8] = b11;
707}
708
709static inline __attribute__((always_inline)) void isotrope_laplacian(dt_aligned_pixel_simd_t kernel[9])
710{
711 // see in https://eng.aurelienpierre.com/2021/03/rotation-invariant-laplacian-for-2d-grids/#Second-order-isotropic-finite-differences
712 // for references (Oono & Puri)
713 const dt_aligned_pixel_simd_t corner = dt_simd_set1(0.25f);
714 const dt_aligned_pixel_simd_t edge = dt_simd_set1(0.5f);
715 const dt_aligned_pixel_simd_t center = dt_simd_set1(-3.f);
716 kernel[0] = corner;
717 kernel[1] = edge;
718 kernel[2] = corner;
719 kernel[3] = edge;
720 kernel[4] = center;
721 kernel[5] = edge;
722 kernel[6] = corner;
723 kernel[7] = edge;
724 kernel[8] = corner;
725}
726
727static inline __attribute__((always_inline)) void compute_kernel(
728 const dt_aligned_pixel_simd_t c2, const dt_aligned_pixel_simd_t cos_theta_sin_theta,
729 const dt_aligned_pixel_simd_t cos_theta2, const dt_aligned_pixel_simd_t sin_theta2,
730 const dt_isotropy_t isotropy_type, dt_aligned_pixel_simd_t kernel[9])
731{
732 // Build the matrix of rotation with anisotropy
733
734 switch(isotropy_type)
735 {
737 default:
738 {
739 isotrope_laplacian(kernel);
740 break;
741 }
743 {
744 dt_aligned_pixel_simd_t a[2][2] = { { dt_simd_set1(0.f) } };
745 rotation_matrix_isophote(c2, cos_theta_sin_theta, cos_theta2, sin_theta2, a);
746 build_matrix(a, kernel);
747 break;
748 }
750 {
751 dt_aligned_pixel_simd_t a[2][2] = { { dt_simd_set1(0.f) } };
752 rotation_matrix_gradient(c2, cos_theta_sin_theta, cos_theta2, sin_theta2, a);
753 build_matrix(a, kernel);
754 break;
755 }
756 }
757}
758
760static inline void heat_PDE_diffusion(const float *const restrict high_freq, const float *const restrict low_freq,
761 const uint8_t *const restrict mask, const int has_mask,
762 float *const restrict output, const size_t width,
763 const size_t height, const dt_aligned_pixel_simd_t anisotropy,
764 const dt_isotropy_t isotropy_type[4],
765 const float variance_threshold, const int mult,
766 const float normalized_regularization,
767 const dt_aligned_pixel_simd_t ABCD, const float strength,
768 const int use_nontemporal)
769{
770 // Simultaneous inpainting for image structure and texture using anisotropic heat transfer model
771 // https://www.researchgate.net/publication/220663968
772 // modified as follow :
773 // * apply it in a multi-scale wavelet setup : we basically solve it twice, on the wavelets LF and HF layers.
774 // * replace the manual texture direction/distance selection by an automatic detection similar to the structure one,
775 // * generalize the framework for isotropic diffusion and anisotropic weighted on the isophote direction
776 // * add an HF-band energy regularization to better avoid edges.
777 // The sharpness setting mimics the contrast equalizer effect by simply multiplying the HF by some gain.
778
779 float *const restrict out = DT_IS_ALIGNED(output);
780 const float *const restrict LF = DT_IS_ALIGNED(low_freq);
781 const float *const restrict HF = DT_IS_ALIGNED(high_freq);
782 const dt_aligned_pixel_simd_t zero = dt_simd_set1(0.f);
783 const dt_aligned_pixel_simd_t flt_min = dt_simd_set1(1e-8f);
784 const dt_aligned_pixel_simd_t variance_threshold_v = dt_simd_set1(variance_threshold);
785 const dt_aligned_pixel_simd_t normalized_regularization_v = dt_simd_set1(normalized_regularization);
786 const dt_aligned_pixel_simd_t strength_v = dt_simd_set1(strength);
787
789 for(size_t row = 0; row < height; ++row)
790 {
791 // interleave the order in which we process the rows so that we minimize cache misses
792 const size_t i = dwt_interleave_rows(row, height, mult);
793 // compute the 'above' and 'below' coordinates, clamping them to the image, once for the entire row
794 const size_t i_neighbours[3]
795 = { MAX((int)(i - mult * H), (int)0) * width, // x - mult
796 i * width, // x
797 MIN((int)(i + mult * H), (int)height - 1) * width }; // x + mult
798 for(size_t j = 0; j < width; ++j)
799 {
800 const size_t idx = (i * width + j);
801 const size_t index = idx * 4;
802 const uint8_t opacity = (has_mask) ? mask[idx] : 1;
803
804 if(opacity)
805 {
806 // non-local neighbours coordinates
807 const size_t j_neighbours[3]
808 = { MAX((int)(j - mult * H), (int)0), // y - mult
809 j, // y
810 MIN((int)(j + mult * H), (int)width - 1) }; // y + mult
811
812 // fetch non-local pixels and store them locally and contiguously
813 dt_aligned_pixel_simd_t neighbour_pixel_HF[9];
814 dt_aligned_pixel_simd_t neighbour_pixel_LF[9];
815 dt_aligned_pixel_simd_t energy = zero;
816
817 for(size_t ii = 0; ii < 3; ii++)
818 for(size_t jj = 0; jj < 3; jj++)
819 {
820 const size_t neighbor = 4 * (i_neighbours[ii] + j_neighbours[jj]);
821 const dt_aligned_pixel_simd_t hf_value = dt_load_simd_aligned(HF + neighbor);
822 const dt_aligned_pixel_simd_t lf_value = dt_load_simd_aligned(LF + neighbor);
823 neighbour_pixel_HF[3 * ii + jj] = hf_value;
824 neighbour_pixel_LF[3 * ii + jj] = lf_value;
825 // Clamp LF to a strictly positive floor to avoid divide-by-zero in
826 // the HF/LF energy estimate without branching per channel.
827 const dt_aligned_pixel_simd_t safe_lf = dt_simd_max_zero(lf_value - flt_min) + flt_min;
828 const dt_aligned_pixel_simd_t ratio = hf_value / safe_lf;
829 energy += ratio * ratio;
830 }
831
832 // normalized_regularization already folds together the user
833 // regularization, the 3x3-support averaging factor, the physical blur
834 // radius carried by the current wavelet band and its scale normalization.
835 energy = dt_simd_max_zero(variance_threshold_v + energy * normalized_regularization_v - flt_min) + flt_min;
836
837 // build the local anisotropic convolution filters for gradients and laplacians
838 dt_aligned_pixel_simd_t lf_gradient[2], hf_gradient[2]; // x, y for each channel
839 find_gradients(neighbour_pixel_LF, lf_gradient);
840 find_gradients(neighbour_pixel_HF, hf_gradient);
841
842 // c² in https://www.researchgate.net/publication/220663968
843 dt_aligned_pixel_simd_t c2[4];
844 dt_aligned_pixel_simd_t grad_x = lf_gradient[0];
845 dt_aligned_pixel_simd_t grad_y = lf_gradient[1];
846 dt_aligned_pixel_simd_t c2_first = zero;
847 dt_aligned_pixel_simd_t c2_third = zero;
848 dt_aligned_pixel_simd_t cos_theta_grad_sq = zero;
849 dt_aligned_pixel_simd_t sin_theta_grad_sq = zero;
850 dt_aligned_pixel_simd_t cos_theta_sin_theta_grad = zero;
852 {
853 const float magnitude_grad = dt_fast_hypotf(grad_x[c], grad_y[c]);
854 c2_first[c] = -magnitude_grad * anisotropy[0];
855 c2_third[c] = -magnitude_grad * anisotropy[2];
856 // Compute cos/sin(arg(grad)) with a branchless normalization, forcing
857 // arg(grad)=0 when magnitude is zero.
858 const float nonzero = (magnitude_grad != 0.f);
859 const float inv_mag = 1.f / (magnitude_grad + (1.f - nonzero));
860 grad_x[c] = grad_x[c] * inv_mag + (1.f - nonzero); // cos(0)
861 grad_y[c] = grad_y[c] * inv_mag; // sin(0)
862 // Warning : now gradient = { cos(arg(grad)) , sin(arg(grad)) }
863 cos_theta_grad_sq[c] = sqf(grad_x[c]);
864 sin_theta_grad_sq[c] = sqf(grad_y[c]);
865 cos_theta_sin_theta_grad[c] = grad_x[c] * grad_y[c];
866 }
867
868 c2[0] = c2_first;
869 c2[2] = c2_third;
870 dt_aligned_pixel_simd_t lapl_x = hf_gradient[0];
871 dt_aligned_pixel_simd_t lapl_y = hf_gradient[1];
872 dt_aligned_pixel_simd_t c2_second = zero;
873 dt_aligned_pixel_simd_t c2_fourth = zero;
874 dt_aligned_pixel_simd_t cos_theta_lapl_sq = zero;
875 dt_aligned_pixel_simd_t sin_theta_lapl_sq = zero;
876 dt_aligned_pixel_simd_t cos_theta_sin_theta_lapl = zero;
878 {
879 const float magnitude_lapl = dt_fast_hypotf(lapl_x[c], lapl_y[c]);
880 c2_second[c] = -magnitude_lapl * anisotropy[1];
881 c2_fourth[c] = -magnitude_lapl * anisotropy[3];
882 // Compute cos/sin(arg(lapl)) with a branchless normalization, forcing
883 // arg(lapl)=0 when magnitude is zero.
884 const float nonzero = (magnitude_lapl != 0.f);
885 const float inv_mag = 1.f / (magnitude_lapl + (1.f - nonzero));
886 lapl_x[c] = lapl_x[c] * inv_mag + (1.f - nonzero); // cos(0)
887 lapl_y[c] = lapl_y[c] * inv_mag; // sin(0)
888 // Warning : now laplacian = { cos(arg(lapl)) , sin(arg(lapl)) }
889 cos_theta_lapl_sq[c] = sqf(lapl_x[c]);
890 sin_theta_lapl_sq[c] = sqf(lapl_y[c]);
891 cos_theta_sin_theta_lapl[c] = lapl_x[c] * lapl_y[c];
892 }
893 c2[1] = c2_second;
894 c2[3] = c2_fourth;
895
896 // elements of c2 need to be expf(mag*anistropy), but we haven't applied the expf() yet. Do that now.
897 for(size_t k = 0; k < 4; k++)
898 for_each_channel(c) c2[k][c] = dt_fast_expf(c2[k][c]);
899
900 dt_aligned_pixel_simd_t kern_first[9], kern_second[9], kern_third[9], kern_fourth[9];
901 compute_kernel(c2[0], cos_theta_sin_theta_grad, cos_theta_grad_sq, sin_theta_grad_sq, isotropy_type[0],
902 kern_first);
903 compute_kernel(c2[1], cos_theta_sin_theta_lapl, cos_theta_lapl_sq, sin_theta_lapl_sq, isotropy_type[1],
904 kern_second);
905 compute_kernel(c2[2], cos_theta_sin_theta_grad, cos_theta_grad_sq, sin_theta_grad_sq, isotropy_type[2],
906 kern_third);
907 compute_kernel(c2[3], cos_theta_sin_theta_lapl, cos_theta_lapl_sq, sin_theta_lapl_sq, isotropy_type[3],
908 kern_fourth);
909
910 dt_aligned_pixel_simd_t derivatives[4] = { zero, zero, zero, zero };
911 // Convolve filters and accumulate the local HF band energy over the
912 // current 3x3 support. This is not a statistical variance estimator:
913 // HF is a band-pass residual, so we normalize each sample by the
914 // corresponding LF value before squaring it, then normalize the summed
915 // ratio by the physical kernel-variance increment of the current
916 // wavelet band.
917 for(size_t k = 0; k < 9; k++)
918 {
919 derivatives[0] = kern_first[k] * neighbour_pixel_LF[k] + derivatives[0];
920 derivatives[1] = kern_second[k] * neighbour_pixel_LF[k] + derivatives[1];
921 derivatives[2] = kern_third[k] * neighbour_pixel_HF[k] + derivatives[2];
922 derivatives[3] = kern_fourth[k] * neighbour_pixel_HF[k] + derivatives[3];
923 }
924
925 // compute the update
926 dt_aligned_pixel_simd_t update = derivatives[0] * ABCD[0];
927 update = derivatives[1] * ABCD[1] + update;
928 update = derivatives[2] * ABCD[2] + update;
929 update = derivatives[3] * ABCD[3] + update;
930 const dt_aligned_pixel_simd_t acc = neighbour_pixel_HF[4] * strength_v + update / energy;
931
932 if(use_nontemporal)
933 dt_store_simd_nontemporal(out + index, dt_simd_max_zero(acc + neighbour_pixel_LF[4]));
934 else
935 dt_store_simd_aligned(out + index, dt_simd_max_zero(acc + neighbour_pixel_LF[4]));
936 }
937 else
938 {
939 // only copy input to output, do nothing
940 if(use_nontemporal)
941 dt_store_simd_nontemporal(out + index, dt_simd_max_zero(dt_load_simd_aligned(HF + index)
942 + dt_load_simd_aligned(LF + index)));
943 else
944 dt_store_simd_aligned(out + index, dt_simd_max_zero(dt_load_simd_aligned(HF + index)
945 + dt_load_simd_aligned(LF + index)));
946 }
947 }
948 }
949
950
951 if(use_nontemporal)
952 dt_omploop_sfence(); // ensure the final nontemporal writeback completes before the caller reads out
953}
954
955static inline float compute_anisotropy_factor(const float user_param)
956{
957 // compute the inverse of the K param in c evaluation from
958 // https://www.researchgate.net/publication/220663968
959 // but in a perceptually-even way, for better GUI interaction
960 return sqf(user_param);
961}
962
963#if DEBUG_DUMP_PFM
965static void dump_PFM(const char *filename, const float* out, const uint32_t w, const uint32_t h)
966{
967 FILE *f = g_fopen(filename, "wb");
968 fprintf(f, "PF\n%d %d\n-1.0\n", w, h);
969 for(int j = h - 1 ; j >= 0 ; j--)
970 for(int i = 0 ; i < w ; i++)
971 for(int c = 0 ; c < 3 ; c++)
972 fwrite(out + (j * w + i) * 4 + c, 1, sizeof(float), f);
973 fclose(f);
974}
975#endif
976
978static inline int wavelets_process(const float *const restrict in, float *const restrict reconstructed,
979 const uint8_t *const restrict mask, const size_t width,
980 const size_t height, const dt_iop_diffuse_data_t *const data,
981 const float zoom, const int scales,
982 const int has_mask,
983 float *const restrict HF[MAX_NUM_SCALES],
984 float *const restrict LF_odd,
985 float *const restrict LF_even)
986{
987 const dt_aligned_pixel_simd_t anisotropy
992
993 const dt_isotropy_t DT_ALIGNED_PIXEL isotropy_type[4]
998
999 const float regularization = powf(10.f, data->regularization) - 1.f;
1000 const float variance_threshold = powf(10.f, data->variance_threshold);
1001
1002 // À trous decimated wavelet decompose
1003 // there is a paper from a guy we know that explains it : https://jo.dreggn.org/home/2010_atrous.pdf
1004 // the wavelets decomposition here is the same as the equalizer/atrous module,
1005 float *restrict residual; // will store the temp buffer containing the last step of blur
1006 // allocate a one-row temporary buffer for the decomposition
1007 size_t padded_size;
1008 float *const tempbuf = dt_pixelpipe_cache_alloc_perthread_float(4 * width, &padded_size); //TODO: alloc in caller
1009 if(IS_NULL_PTR(tempbuf)) return 1;
1010
1011 for(int s = 0; s < scales; ++s)
1012 {
1013 const int mult = 1 << s;
1014
1015 const float *restrict buffer_in;
1016 float *restrict buffer_out;
1017
1018 if(s == 0)
1019 {
1020 buffer_in = in;
1021 buffer_out = LF_odd;
1022 }
1023 else if(s % 2 != 0)
1024 {
1025 buffer_in = LF_odd;
1026 buffer_out = LF_even;
1027 }
1028 else
1029 {
1030 buffer_in = LF_even;
1031 buffer_out = LF_odd;
1032 }
1033
1034 decompose_2D_Bspline(buffer_in, HF[s], buffer_out, width, height, mult, tempbuf, padded_size);
1035
1036 residual = buffer_out;
1037
1038#if DEBUG_DUMP_PFM
1039 char name[64];
1040 sprintf(name, "/tmp/scale-input-%i.pfm", s);
1041 dump_PFM(name, buffer_in, width, height);
1042
1043 sprintf(name, "/tmp/scale-blur-%i.pfm", s);
1044 dump_PFM(name, buffer_out, width, height);
1045#endif
1046 }
1048
1049 // will store the temp buffer NOT containing the last step of blur
1050 float *restrict temp = (residual == LF_even) ? LF_odd : LF_even;
1051 int count = 0;
1052
1053 for(int s = scales - 1; s > -1; --s)
1054 {
1055 const int mult = 1 << s;
1056 const float current_radius = equivalent_sigma_at_step(B_SPLINE_SIGMA, s);
1057 const float real_radius = current_radius * zoom;
1058
1059#if DIFFUSE_V3
1060 const float normalized_regularization =
1061 (data->normalize_band_energy)
1062 ? regularization * sqf(real_radius) / 9.f
1063 : regularization / 9.f;
1064#else
1065 const float normalized_regularization = regularization / 9.f * sqf(real_radius);
1066#endif
1067
1068 const float norm = expf(-sqf(real_radius - (float)data->radius_center) / sqf(data->radius));
1069
1070 const dt_aligned_pixel_simd_t ABCD = { data->first * KAPPA * norm,
1071 data->second * KAPPA * norm,
1072 data->third * KAPPA * norm,
1073 data->fourth * KAPPA * norm };
1074 const float strength = data->sharpness * norm + 1.f;
1075
1076 const float *restrict buffer_in;
1077 float *restrict buffer_out;
1078
1079 if(count == 0)
1080 {
1081 buffer_in = residual;
1082 buffer_out = temp;
1083 }
1084 else if(count % 2 != 0)
1085 {
1086 buffer_in = temp;
1087 buffer_out = residual;
1088 }
1089 else
1090 {
1091 buffer_in = residual;
1092 buffer_out = temp;
1093 }
1094
1095 if(s == 0) buffer_out = reconstructed;
1096
1097 heat_PDE_diffusion(HF[s], buffer_in, mask, has_mask, buffer_out, width, height,
1098 anisotropy, isotropy_type, variance_threshold, mult,
1099 normalized_regularization, ABCD, strength, (s == 0));
1100
1101 count++;
1102 }
1103
1104 return 0;
1105}
1106
1107
1109static inline void build_mask(const float *const restrict input, uint8_t *const restrict mask,
1110 const float threshold, const size_t width, const size_t height)
1111{
1112 __OMP_PARALLEL_FOR_SIMD__(aligned(mask, input : 64))
1113 for(size_t k = 0; k < height * width * 4; k += 4)
1114 {
1115 // TRUE if any channel is above threshold
1116 mask[k / 4] = (input[k] > threshold || input[k + 1] > threshold || input[k + 2] > threshold);
1117 }
1118
1119}
1120
1122static inline void inpaint_mask(float *const restrict inpainted, const float *const restrict original,
1123 const uint8_t *const restrict mask, const size_t width,
1124 const size_t height)
1125{
1126 // init the reconstruction with noise inside the masked areas
1128 for(size_t k = 0; k < height * width * 4; k += 4)
1129 {
1130 if(mask[k / 4])
1131 {
1132 const uint32_t i = k / width;
1133 const uint32_t j = k - i;
1134 uint32_t DT_ALIGNED_ARRAY state[4]
1135 = { splitmix32(j + 1), splitmix32((uint64_t)(j + 1) * (i + 3)),
1136 splitmix32(1337), splitmix32(666) };
1141
1142 for_four_channels(c, aligned(inpainted, original, state:64))
1143 inpainted[k + c] = fabsf(gaussian_noise(original[k + c], original[k + c], i % 2 || j % 2, state));
1144 }
1145 else
1146 {
1147 for_four_channels(c, aligned(original, inpainted:64))
1148 inpainted[k + c] = original[k + c];
1149 }
1150 }
1151
1152}
1153
1156 const void *const restrict ivoid, void *const restrict ovoid)
1157{
1158 const dt_iop_roi_t *const roi_in = &piece->roi_in;
1159 const dt_iop_roi_t *const roi_out = &piece->roi_out;
1160 const dt_iop_diffuse_data_t *const data = (dt_iop_diffuse_data_t *)piece->data;
1161
1162 float *restrict in = DT_IS_ALIGNED((float *const restrict)ivoid);
1163 float *const restrict out = DT_IS_ALIGNED((float *const restrict)ovoid);
1164
1165 float *const restrict temp1 = dt_pixelpipe_cache_alloc_align_float((size_t)roi_out->width * roi_out->height * 4, pipe);
1166 float *const restrict temp2 = dt_pixelpipe_cache_alloc_align_float((size_t)roi_out->width * roi_out->height * 4, pipe);
1167
1168 float *restrict temp_in = NULL;
1169 float *restrict temp_out = NULL;
1170 int err = 0;
1171
1172 uint8_t *const restrict mask = dt_pixelpipe_cache_alloc_align(
1173 sizeof(uint8_t) * roi_out->width * roi_out->height,
1174 pipe);
1175
1176 const float zoom = dt_dev_get_module_scale(pipe, roi_in);
1177 const float final_radius = (data->radius + data->radius_center) * 2.f / zoom;
1178 // No legacy iteration remap is applied here anymore. The current solver uses
1179 // the historical a-trous band order and kernel-variance increments exactly,
1180 // so any extra factor would be content-dependent and belong to pixel math,
1181 // not to the user parameter itself.
1182 const int iterations = MAX((int)ceilf((float)data->iterations), 1);
1183 const int diffusion_scales = num_steps_to_reach_equivalent_sigma(B_SPLINE_SIGMA, final_radius);
1184 const int scales = CLAMP(diffusion_scales, 1, MAX_NUM_SCALES);
1185
1186 gboolean out_of_memory = (IS_NULL_PTR(temp1)) || (IS_NULL_PTR(temp2));
1187 // One full-resolution buffer per stored wavelet band.
1188 float *restrict HF[MAX_NUM_SCALES] = { NULL };
1189 for(int s = 0; s < scales; s++)
1190 {
1191 HF[s] = dt_pixelpipe_cache_alloc_align_float(roi_out->width * roi_out->height * 4, pipe);
1192 if(!HF[s]) out_of_memory = TRUE;
1193 }
1194 // Two ping-pong low-pass buffers reused by the decomposition/synthesis.
1195 float *const restrict LF_odd = dt_pixelpipe_cache_alloc_align_float(roi_out->width * roi_out->height * 4, pipe);
1196 float *const restrict LF_even = dt_pixelpipe_cache_alloc_align_float(roi_out->width * roi_out->height * 4, pipe);
1197
1198 // PAUSE !
1199 // check that all buffers exist before processing,
1200 // because we use a lot of memory here.
1201 if(IS_NULL_PTR(mask) || IS_NULL_PTR(temp1) || IS_NULL_PTR(temp2) || IS_NULL_PTR(LF_odd) || IS_NULL_PTR(LF_even) || out_of_memory)
1202 {
1203 err = 1;
1204 goto error;
1205 }
1206
1207 const int has_mask = (data->threshold > 0.f);
1208
1209 if(has_mask)
1210 {
1211 // build a boolean mask, TRUE where image is above threshold, FALSE otherwise
1212 build_mask(in, mask, data->threshold, roi_out->width, roi_out->height);
1213
1214 // init the inpainting area with noise
1215 inpaint_mask(temp1, in, mask, roi_out->width, roi_out->height);
1216
1217 in = temp1;
1218 }
1219
1220 for(int it = 0; it < iterations; it++)
1221 {
1222 if(it == 0)
1223 {
1224 temp_in = in;
1225 temp_out = temp2;
1226 }
1227 else if(it % 2 == 0)
1228 {
1229 temp_in = temp1;
1230 temp_out = temp2;
1231 }
1232 else
1233 {
1234 temp_in = temp2;
1235 temp_out = temp1;
1236 }
1237
1238 if(it == (int)iterations - 1)
1239 temp_out = out;
1240
1241 if(wavelets_process(temp_in, temp_out, mask, roi_out->width, roi_out->height,
1242 data, zoom, scales, has_mask, HF, LF_odd, LF_even))
1243 {
1244 err = 1;
1245 goto error;
1246 }
1247 }
1248
1249error:
1255 for(int s = 0; s < scales; s++)
1256 if(HF[s]) dt_pixelpipe_cache_free_align(HF[s]);
1257 return err;
1258}
1259
1260#if HAVE_OPENCL
1261static inline cl_int wavelets_process_cl(const int devid, cl_mem in, cl_mem reconstructed, cl_mem mask,
1262 const size_t sizes[3], const int width, const int height,
1263 const dt_iop_diffuse_data_t *const data,
1265 const float zoom, const int scales,
1266 const int has_mask,
1267 cl_mem HF[MAX_NUM_SCALES],
1268 cl_mem LF_odd,
1269 cl_mem LF_even)
1270{
1271 cl_int err = -999;
1272
1273 const dt_aligned_pixel_simd_t anisotropy
1278
1279 /*
1280 fprintf(stdout, "anisotropy : %f ; %f ; %f ; %f \n",
1281 anisotropy[0], anisotropy[1], anisotropy[2], anisotropy[3]);
1282 */
1283
1284 const dt_isotropy_t DT_ALIGNED_PIXEL isotropy_type[4]
1289
1290 /*
1291 fprintf(stdout, "type : %d ; %d ; %d ; %d \n",
1292 isotropy_type[0], isotropy_type[1], isotropy_type[2], isotropy_type[3]);
1293 */
1294
1295 const float regularization = powf(10.f, data->regularization) - 1.f;
1296 const float variance_threshold = powf(10.f, data->variance_threshold);
1297 // Same a-trous decomposition as the CPU path, mirrored in OpenCL.
1298 cl_mem residual;
1299
1300 for(int s = 0; s < scales; ++s)
1301 {
1302 const int mult = 1 << s;
1303
1304 cl_mem buffer_in;
1305 cl_mem buffer_out;
1306
1307 if(s == 0)
1308 {
1309 buffer_in = in;
1310 buffer_out = LF_odd;
1311 }
1312 else if(s % 2 != 0)
1313 {
1314 buffer_in = LF_odd;
1315 buffer_out = LF_even;
1316 }
1317 else
1318 {
1319 buffer_in = LF_even;
1320 buffer_out = LF_odd;
1321 }
1322
1323 // Compute wavelets low-frequency scales
1324 const int clamp_lf = 1;
1325 int hblocksize;
1326 dt_opencl_local_buffer_t hlocopt = (dt_opencl_local_buffer_t){ .xoffset = 2 * mult, .xfactor = 1,
1327 .yoffset = 0, .yfactor = 1,
1328 .cellsize = 4 * sizeof(float), .overhead = 0,
1329 .sizex = 1 << 16, .sizey = 1 };
1331 hblocksize = hlocopt.sizex;
1332 else
1333 hblocksize = 1;
1334
1335 // Keep the same separable order as the CPU path: vertical pass first,
1336 // store its intermediate into HF[s], then horizontal pass builds LF.
1337 int vblocksize;
1338 dt_opencl_local_buffer_t vlocopt = (dt_opencl_local_buffer_t){ .xoffset = 0, .xfactor = 1,
1339 .yoffset = 2 * mult, .yfactor = 1,
1340 .cellsize = 4 * sizeof(float), .overhead = 0,
1341 .sizex = 1, .sizey = 1 << 16 };
1343 vblocksize = vlocopt.sizey;
1344 else
1345 vblocksize = 1;
1346
1347 if(vblocksize > 1)
1348 {
1349 const size_t vertical_sizes[3] = { ROUNDUPDWD(width, devid), ROUNDUP(height, vblocksize), 1 };
1350 const size_t vertical_local[3] = { 1, vblocksize, 1 };
1351 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_vertical_local, 0, sizeof(cl_mem), (void *)&buffer_in);
1352 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_vertical_local, 1, sizeof(cl_mem), (void *)&HF[s]);
1353 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_vertical_local, 2, sizeof(int), (void *)&width);
1354 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_vertical_local, 3, sizeof(int), (void *)&height);
1355 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_vertical_local, 4, sizeof(int), (void *)&mult);
1356 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_vertical_local, 5, sizeof(int), (void *)&clamp_lf);
1358 (vblocksize + 4 * mult) * 4 * sizeof(float), NULL);
1360 vertical_sizes, vertical_local);
1361 }
1362 else
1363 {
1364 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_vertical, 0, sizeof(cl_mem), (void *)&buffer_in);
1365 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_vertical, 1, sizeof(cl_mem), (void *)&HF[s]);
1366 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_vertical, 2, sizeof(int), (void *)&width);
1367 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_vertical, 3, sizeof(int), (void *)&height);
1368 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_vertical, 4, sizeof(int), (void *)&mult);
1369 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_vertical, 5, sizeof(int), (void *)&clamp_lf);
1371 }
1372 if(err != CL_SUCCESS) return err;
1373
1374 if(hblocksize > 1)
1375 {
1376 const size_t horizontal_sizes[3] = { ROUNDUP(width, hblocksize), ROUNDUPDHT(height, devid), 1 };
1377 const size_t horizontal_local[3] = { hblocksize, 1, 1 };
1378 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_horizontal_local, 0, sizeof(cl_mem), (void *)&HF[s]);
1379 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_horizontal_local, 1, sizeof(cl_mem), (void *)&buffer_out);
1380 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_horizontal_local, 2, sizeof(int), (void *)&width);
1381 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_horizontal_local, 3, sizeof(int), (void *)&height);
1382 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_horizontal_local, 4, sizeof(int), (void *)&mult);
1383 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_horizontal_local, 5, sizeof(int), (void *)&clamp_lf);
1385 (hblocksize + 4 * mult) * 4 * sizeof(float), NULL);
1387 horizontal_sizes, horizontal_local);
1388 }
1389 else
1390 {
1391 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_horizontal, 0, sizeof(cl_mem), (void *)&HF[s]);
1392 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_horizontal, 1, sizeof(cl_mem), (void *)&buffer_out);
1393 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_horizontal, 2, sizeof(int), (void *)&width);
1394 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_horizontal, 3, sizeof(int), (void *)&height);
1395 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_horizontal, 4, sizeof(int), (void *)&mult);
1396 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_bspline_horizontal, 5, sizeof(int), (void *)&clamp_lf);
1398 }
1399 if(err != CL_SUCCESS) return err;
1400
1401 // Compute wavelets high-frequency scales and backup the maximum of texture over the RGB channels
1402 // Note : HF = detail - LF
1403 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_wavelets_detail, 0, sizeof(cl_mem), (void *)&buffer_in);
1404 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_wavelets_detail, 1, sizeof(cl_mem), (void *)&buffer_out);
1405 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_wavelets_detail, 2, sizeof(cl_mem), (void *)&HF[s]);
1406 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_wavelets_detail, 3, sizeof(int), (void *)&width);
1407 dt_opencl_set_kernel_arg(devid, gd->kernel_filmic_wavelets_detail, 4, sizeof(int), (void *)&height);
1409 if(err != CL_SUCCESS) return err;
1410
1411 residual = buffer_out;
1412 }
1413
1414 // Ping-pong low-pass buffer not currently holding the coarsest residual.
1415 cl_mem temp = (residual == LF_even) ? LF_odd : LF_even;
1416 int count = 0;
1417
1418 for(int s = scales - 1; s > -1; --s)
1419 {
1420 const int mult = 1 << s;
1421 const float current_radius = equivalent_sigma_at_step(B_SPLINE_SIGMA, s);
1422 const float real_radius = current_radius * zoom;
1423
1424#if DIFFUSE_V3
1425 const float normalized_regularization =
1426 (data->normalize_band_energy)
1427 ? regularization * sqf(real_radius) / 9.f
1428 : regularization / 9.f;
1429#else
1430 const float normalized_regularization = regularization / 9.f * sqf(real_radius);
1431#endif
1432
1433 const float norm = expf(-sqf(real_radius - (float)data->radius_center) / sqf(data->radius));
1434
1435 const dt_aligned_pixel_simd_t ABCD = { data->first * KAPPA * norm,
1436 data->second * KAPPA * norm,
1437 data->third * KAPPA * norm,
1438 data->fourth * KAPPA * norm };
1439 const float strength = data->sharpness * norm + 1.f;
1440
1441 cl_mem buffer_in;
1442 cl_mem buffer_out;
1443
1444 if(count == 0)
1445 {
1446 buffer_in = residual;
1447 buffer_out = temp;
1448 }
1449 else if(count % 2 != 0)
1450 {
1451 buffer_in = temp;
1452 buffer_out = residual;
1453 }
1454 else
1455 {
1456 buffer_in = residual;
1457 buffer_out = temp;
1458 }
1459
1460 if(s == 0) buffer_out = reconstructed;
1461
1462 // Compute wavelets low-frequency scales
1463 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_pde, 0, sizeof(cl_mem), (void *)&HF[s]);
1464 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_pde, 1, sizeof(cl_mem), (void *)&buffer_in);
1465 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_pde, 2, sizeof(cl_mem), (void *)&mask);
1466 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_pde, 3, sizeof(int), (void *)&has_mask);
1467 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_pde, 4, sizeof(cl_mem), (void *)&buffer_out);
1468 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_pde, 5, sizeof(int), (void *)&width);
1469 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_pde, 6, sizeof(int), (void *)&height);
1470 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_pde, 7, 4 * sizeof(float), (void *)&anisotropy);
1471 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_pde, 8, 4 * sizeof(dt_isotropy_t), (void *)&isotropy_type);
1472 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_pde, 9, sizeof(float), (void *)&normalized_regularization);
1473 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_pde, 10, sizeof(float), (void *)&variance_threshold);
1474 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_pde, 11, sizeof(int), (void *)&mult);
1475 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_pde, 12, 4 * sizeof(float), (void *)&ABCD);
1476 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_pde, 13, sizeof(float), (void *)&strength);
1477 err = dt_opencl_enqueue_kernel_2d(devid, gd->kernel_diffuse_pde, sizes);
1478 if(err != CL_SUCCESS) return err;
1479
1480 count++;
1481 }
1482
1483 return err;
1484}
1485
1486int process_cl(struct dt_iop_module_t *self, const dt_dev_pixelpipe_t *pipe, const dt_dev_pixelpipe_iop_t *piece, cl_mem dev_in, cl_mem dev_out)
1487{
1488 const dt_iop_roi_t *const roi_in = &piece->roi_in;
1489 const dt_iop_roi_t *const roi_out = &piece->roi_out;
1490 const dt_iop_diffuse_data_t *const data = (dt_iop_diffuse_data_t *)piece->data;
1492
1493 int out_of_memory = FALSE;
1494
1495 cl_int err = -999;
1496
1497 const int devid = pipe->devid;
1498 const int width = roi_in->width;
1499 const int height = roi_in->height;
1500
1501 size_t sizes[] = { ROUNDUPDWD(width, devid), ROUNDUPDHT(height, devid), 1 };
1502
1503 cl_mem in = dev_in;
1504
1505 cl_mem temp1 = dt_opencl_alloc_device(devid, sizes[0], sizes[1], sizeof(float) * 4);
1506 cl_mem temp2 = dt_opencl_alloc_device(devid, sizes[0], sizes[1], sizeof(float) * 4);
1507
1508 cl_mem temp_in = NULL;
1509 cl_mem temp_out = NULL;
1510
1511 cl_mem mask = dt_opencl_alloc_device(devid, sizes[0], sizes[1], sizeof(uint8_t));
1512
1513 const float zoom = dt_dev_get_module_scale(pipe, roi_in);
1514 const float final_radius = (data->radius + data->radius_center) * 2.f / zoom;
1515 // See the CPU path above: iterations stay in user space because the current
1516 // solver already matches the historical a-trous band ordering and kernel
1517 // variance increments. There is no content-independent remap left to apply.
1518 const int iterations = MAX((int)ceilf((float)data->iterations), 1);
1519 const int diffusion_scales = num_steps_to_reach_equivalent_sigma(B_SPLINE_SIGMA, final_radius);
1520 const int scales = CLAMP(diffusion_scales, 1, MAX_NUM_SCALES);
1521 // One device buffer per stored wavelet band.
1522 cl_mem HF[MAX_NUM_SCALES] = { NULL };
1523 for(int s = 0; s < scales; s++)
1524 {
1525 HF[s] = dt_opencl_alloc_device(devid, sizes[0], sizes[1], sizeof(float) * 4);
1526 if(!HF[s]) out_of_memory = TRUE;
1527 }
1528 // Two low-pass ping-pong buffers reused across all scales.
1529 cl_mem LF_even = dt_opencl_alloc_device(devid, sizes[0], sizes[1], sizeof(float) * 4);
1530 cl_mem LF_odd = dt_opencl_alloc_device(devid, sizes[0], sizes[1], sizeof(float) * 4);
1531
1532 // PAUSE !
1533 // check that all buffers exist before processing,
1534 // because we use a lot of memory here.
1535 if(IS_NULL_PTR(mask) || IS_NULL_PTR(temp1) || IS_NULL_PTR(temp2) || IS_NULL_PTR(LF_odd) || IS_NULL_PTR(LF_even) || out_of_memory)
1536 {
1537 err = CL_MEM_OBJECT_ALLOCATION_FAILURE;
1538 goto error;
1539 }
1540
1541 const int has_mask = (data->threshold > 0.f);
1542
1543 if(has_mask)
1544 {
1545 // build a boolean mask, TRUE where image is above threshold, FALSE otherwise
1546 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_build_mask, 0, sizeof(cl_mem), (void *)&in);
1547 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_build_mask, 1, sizeof(cl_mem), (void *)&mask);
1548 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_build_mask, 2, sizeof(float), (void *)&data->threshold);
1549 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_build_mask, 3, sizeof(int), (void *)&roi_out->width);
1550 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_build_mask, 4, sizeof(int), (void *)&roi_out->height);
1552 if(err != CL_SUCCESS) goto error;
1553
1554 // init the inpainting area with noise
1555 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_inpaint_mask, 0, sizeof(cl_mem), (void *)&temp1);
1556 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_inpaint_mask, 1, sizeof(cl_mem), (void *)&in);
1557 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_inpaint_mask, 2, sizeof(cl_mem), (void *)&mask);
1558 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_inpaint_mask, 3, sizeof(int), (void *)&roi_out->width);
1559 dt_opencl_set_kernel_arg(devid, gd->kernel_diffuse_inpaint_mask, 4, sizeof(int), (void *)&roi_out->height);
1561 if(err != CL_SUCCESS) goto error;
1562
1563 in = temp1;
1564 }
1565
1566 for(int it = 0; it < iterations; it++)
1567 {
1568 if(it == 0)
1569 {
1570 temp_in = in;
1571 temp_out = temp2;
1572 }
1573 else if(it % 2 == 0)
1574 {
1575 temp_in = temp1;
1576 temp_out = temp2;
1577 }
1578 else
1579 {
1580 temp_in = temp2;
1581 temp_out = temp1;
1582 }
1583
1584 if(it == (int)iterations - 1) temp_out = dev_out;
1585 err = wavelets_process_cl(devid, temp_in, temp_out, mask, sizes, width, height,
1586 data, gd, zoom, scales, has_mask, HF, LF_odd, LF_even);
1587 if(err != CL_SUCCESS) goto error;
1588 }
1589
1590 // cleanup and exit on success
1596 for(int s = 0; s < scales; s++) dt_opencl_release_mem_object(HF[s]);
1597 return TRUE;
1598
1599error:
1605 for(int s = 0; s < scales; s++) dt_opencl_release_mem_object(HF[s]);
1606
1607 dt_print(DT_DEBUG_OPENCL, "[opencl_diffuse] couldn't enqueue kernel! %d\n", err);
1608 return FALSE;
1609}
1610
1612{
1613 const int program = 33; // diffuse.cl in programs.conf
1615
1616 module->data = gd;
1617 gd->kernel_diffuse_build_mask = dt_opencl_create_kernel(program, "build_mask");
1618 gd->kernel_diffuse_inpaint_mask = dt_opencl_create_kernel(program, "inpaint_mask");
1619 gd->kernel_diffuse_pde = dt_opencl_create_kernel(program, "diffuse_pde");
1620
1621 const int wavelets = 35; // bspline.cl, from programs.conf
1622 gd->kernel_filmic_bspline_horizontal = dt_opencl_create_kernel(wavelets, "blur_2D_Bspline_horizontal");
1623 gd->kernel_filmic_bspline_vertical = dt_opencl_create_kernel(wavelets, "blur_2D_Bspline_vertical");
1624 gd->kernel_filmic_bspline_horizontal_local = dt_opencl_create_kernel(wavelets, "blur_2D_Bspline_horizontal_local");
1625 gd->kernel_filmic_bspline_vertical_local = dt_opencl_create_kernel(wavelets, "blur_2D_Bspline_vertical_local");
1626 gd->kernel_filmic_wavelets_detail = dt_opencl_create_kernel(wavelets, "wavelets_detail_level");
1627}
1628
1629
1644#endif
1645
1646
1647void gui_init(struct dt_iop_module_t *self)
1648{
1650 self->gui->widget = gtk_box_new(GTK_ORIENTATION_VERTICAL, DT_GUI_BOX_SPACING);
1651
1652 gtk_box_pack_start(GTK_BOX(self->gui->widget), dt_ui_section_label_new(_("properties")), FALSE, FALSE, 0);
1653
1654 g->iterations = dt_bauhaus_slider_from_params(self, "iterations");
1655 dt_bauhaus_slider_set_soft_range(g->iterations, 1., 128);
1656 gtk_widget_set_tooltip_text(g->iterations,
1657 _("more iterations make the effect stronger but the module slower.\n"
1658 "this is analogous to giving more time to the diffusion reaction.\n"
1659 "if you plan on sharpening or inpainting, \n"
1660 "more iterations help reconstruction."));
1661
1662 g->radius_center = dt_bauhaus_slider_from_params(self, "radius_center");
1663 dt_bauhaus_slider_set_soft_range(g->radius_center, 0., 512.);
1664 dt_bauhaus_slider_set_format(g->radius_center, " px");
1665 gtk_widget_set_tooltip_text(
1666 g->radius_center, _("main scale of the diffusion.\n"
1667 "zero makes diffusion act on the finest details more heavily.\n"
1668 "non-zero defines the size of the details to diffuse heavily.\n"
1669 "for deblurring and denoising, set to zero.\n"
1670 "increase to act on local contrast instead."));
1671
1672 g->radius = dt_bauhaus_slider_from_params(self, "radius");
1673 dt_bauhaus_slider_set_soft_range(g->radius, 1., 512.);
1674 dt_bauhaus_slider_set_format(g->radius, " px");
1675 gtk_widget_set_tooltip_text(
1676 g->radius, _("width of the diffusion around the central radius.\n"
1677 "high values diffuse on a large band of radii.\n"
1678 "low values diffuse closer to the central radius.\n"
1679 "if you plan on deblurring, \n"
1680 "the radius should be around the width of your lens blur."));
1681
1682 GtkWidget *label_speed = dt_ui_section_label_new(_("speed (sharpen \342\206\224 diffuse)"));
1683 gtk_box_pack_start(GTK_BOX(self->gui->widget), label_speed, FALSE, FALSE, 0);
1684
1685 g->first = dt_bauhaus_slider_from_params(self, "first");
1687 dt_bauhaus_slider_set_format(g->first, "%");
1688 gtk_widget_set_tooltip_text(g->first, _("diffusion speed of low-frequency wavelet layers\n"
1689 "in the direction of 1st order anisotropy (set below).\n\n"
1690 "negative values sharpen, \n"
1691 "positive values diffuse and blur, \n"
1692 "zero does nothing."));
1693
1694 g->second = dt_bauhaus_slider_from_params(self, "second");
1695 dt_bauhaus_slider_set_digits(g->second, 4);
1696 dt_bauhaus_slider_set_format(g->second, "%");
1697 gtk_widget_set_tooltip_text(g->second, _("diffusion speed of low-frequency wavelet layers\n"
1698 "in the direction of 2nd order anisotropy (set below).\n\n"
1699 "negative values sharpen, \n"
1700 "positive values diffuse and blur, \n"
1701 "zero does nothing."));
1702
1703 g->third = dt_bauhaus_slider_from_params(self, "third");
1705 dt_bauhaus_slider_set_format(g->third, "%");
1706 gtk_widget_set_tooltip_text(g->third, _("diffusion speed of high-frequency wavelet layers\n"
1707 "in the direction of 3rd order anisotropy (set below).\n\n"
1708 "negative values sharpen, \n"
1709 "positive values diffuse and blur, \n"
1710 "zero does nothing."));
1711
1712 g->fourth = dt_bauhaus_slider_from_params(self, "fourth");
1713 dt_bauhaus_slider_set_digits(g->fourth, 4);
1714 dt_bauhaus_slider_set_format(g->fourth, "%");
1715 gtk_widget_set_tooltip_text(g->fourth, _("diffusion speed of high-frequency wavelet layers\n"
1716 "in the direction of 4th order anisotropy (set below).\n\n"
1717 "negative values sharpen, \n"
1718 "positive values diffuse and blur, \n"
1719 "zero does nothing."));
1720
1721 GtkWidget *label_direction = dt_ui_section_label_new(_("direction"));
1722 gtk_box_pack_start(GTK_BOX(self->gui->widget), label_direction, FALSE, FALSE, 0);
1723
1724 g->anisotropy_first = dt_bauhaus_slider_from_params(self, "anisotropy_first");
1725 dt_bauhaus_slider_set_digits(g->anisotropy_first, 4);
1726 dt_bauhaus_slider_set_format(g->anisotropy_first, "%");
1727 gtk_widget_set_tooltip_text(g->anisotropy_first, _("direction of 1st order speed (set above).\n\n"
1728 "negative values follow gradients more closely, \n"
1729 "positive values rather avoid edges (isophotes), \n"
1730 "zero affects both equally (isotropic)."));
1731
1732 g->anisotropy_second = dt_bauhaus_slider_from_params(self, "anisotropy_second");
1733 dt_bauhaus_slider_set_digits(g->anisotropy_second, 4);
1734 dt_bauhaus_slider_set_format(g->anisotropy_second, "%");
1735 gtk_widget_set_tooltip_text(g->anisotropy_second,_("direction of 2nd order speed (set above).\n\n"
1736 "negative values follow gradients more closely, \n"
1737 "positive values rather avoid edges (isophotes), \n"
1738 "zero affects both equally (isotropic)."));
1739
1740 g->anisotropy_third = dt_bauhaus_slider_from_params(self, "anisotropy_third");
1741 dt_bauhaus_slider_set_digits(g->anisotropy_third, 4);
1742 dt_bauhaus_slider_set_format(g->anisotropy_third, "%");
1743 gtk_widget_set_tooltip_text(g->anisotropy_third,_("direction of 3rd order speed (set above).\n\n"
1744 "negative values follow gradients more closely, \n"
1745 "positive values rather avoid edges (isophotes), \n"
1746 "zero affects both equally (isotropic)."));
1747
1748 g->anisotropy_fourth = dt_bauhaus_slider_from_params(self, "anisotropy_fourth");
1749 dt_bauhaus_slider_set_digits(g->anisotropy_fourth, 4);
1750 dt_bauhaus_slider_set_format(g->anisotropy_fourth, "%");
1751 gtk_widget_set_tooltip_text(g->anisotropy_fourth,_("direction of 4th order speed (set above).\n\n"
1752 "negative values follow gradients more closely, \n"
1753 "positive values rather avoid edges (isophotes), \n"
1754 "zero affects both equally (isotropic)."));
1755
1756 gtk_box_pack_start(GTK_BOX(self->gui->widget), dt_ui_section_label_new(_("edge management")), FALSE, FALSE, 0);
1757
1758 g->sharpness = dt_bauhaus_slider_from_params(self, "sharpness");
1759 dt_bauhaus_slider_set_format(g->sharpness, "%");
1760 gtk_widget_set_tooltip_text(g->sharpness,
1761 _("increase or decrease the sharpness of the highest frequencies.\n"
1762 "can be used to keep details after blooming,\n"
1763 "for standalone sharpening set speed to negative values."));
1764
1765 g->regularization = dt_bauhaus_slider_from_params(self, "regularization");
1766 gtk_widget_set_tooltip_text(g->regularization,
1767 _("define the sensitivity of the variance penalty for edges.\n"
1768 "increase to exclude more edges from diffusion,\n"
1769 "if fringes or halos appear."));
1770
1771 g->variance_threshold = dt_bauhaus_slider_from_params(self, "variance_threshold");
1772 gtk_widget_set_tooltip_text(g->variance_threshold,
1773 _("define the variance threshold between edge amplification and penalty.\n"
1774 "decrease if you want pixels on smooth surfaces get a boost,\n"
1775 "increase if you see noise appear on smooth surfaces or\n"
1776 "if dark areas seem oversharpened compared to bright areas."));
1777
1778
1779 gtk_box_pack_start(GTK_BOX(self->gui->widget), dt_ui_section_label_new(_("diffusion spatiality")), FALSE, FALSE, 0);
1780
1781 g->threshold = dt_bauhaus_slider_from_params(self, "threshold");
1782 dt_bauhaus_slider_set_format(g->threshold, "%");
1783 dt_bauhaus_slider_set_digits(g->threshold, 2);
1784 gtk_widget_set_tooltip_text(g->threshold,
1785 _("luminance threshold for the mask.\n"
1786 "0. disables the luminance masking and applies the module on the whole image.\n"
1787 "any higher value excludes pixels with luminance lower than the threshold.\n"
1788 "this can be used to inpaint highlights."));
1789}
1790// clang-format off
1791// modelines: These editor modelines have been set for all relevant files by tools/update_modelines.py
1792// vim: shiftwidth=2 expandtab tabstop=2 cindent
1793// kate: tab-indents: off; indent-width 2; replace-tabs on; indent-mode cstyle; remove-trailing-spaces modified;
1794// clang-format on
static void error(char *msg)
Definition ashift_lsd.c:202
#define TRUE
Definition ashift_lsd.c:162
#define FALSE
Definition ashift_lsd.c:158
void dt_bauhaus_slider_set_soft_range(GtkWidget *widget, float soft_min, float soft_max)
Definition bauhaus.c:1498
void dt_bauhaus_slider_set_digits(GtkWidget *widget, int val)
Definition bauhaus.c:3339
void dt_bauhaus_slider_set_format(GtkWidget *widget, const char *format)
Definition bauhaus.c:3403
@ DEVELOP_BLEND_CS_RGB_SCENE
Definition blend.h:58
#define B_SPLINE_SIGMA
Definition bspline.h:38
static unsigned int num_steps_to_reach_equivalent_sigma(const float sigma_filter, const float sigma_final)
Definition bspline.h:65
static float equivalent_sigma_at_step(const float sigma, const unsigned int s)
Definition bspline.h:52
static void decompose_2D_Bspline(const float *const restrict in, float *const restrict HF, float *const restrict LF, const size_t width, const size_t height, const int mult, float *const tempbuf, size_t padded_size)
Definition bspline.h:351
return vector dt_simd_set1(valid ?(scaling+NORM_MIN) :NORM_MIN)
@ IOP_CS_RGB
const float f
struct _GtkWidget GtkWidget
GtkWidget, opaque, spelled exactly as GTK spells it.
Definition colorspaces.h:98
const float threshold
const dt_colormatrix_t dt_aligned_pixel_t out
dt_store_simd_aligned(out, dt_mat3x4_mul_vec4(vin, dt_colormatrix_row_to_simd(matrix, 0), dt_colormatrix_row_to_simd(matrix, 1), dt_colormatrix_row_to_simd(matrix, 2)))
static const int row
static float strength(float value, float strength)
Definition colorzones.c:431
static unsigned int splitmix32(const unsigned long seed)
static float xoshiro128plus(uint state[4])
void dt_iop_params_t
Definition dev_history.h:43
const char ** description(struct dt_iop_module_t *self)
Definition diffuse.c:174
int default_group()
Definition diffuse.c:187
static __DT_CLONE_TARGETS__ int wavelets_process(const float *const restrict in, float *const restrict reconstructed, const uint8_t *const restrict mask, const size_t width, const size_t height, const dt_iop_diffuse_data_t *const data, const float zoom, const int scales, const int has_mask, float *const restrict HF[10], float *const restrict LF_odd, float *const restrict LF_even)
Definition diffuse.c:978
__DT_CLONE_TARGETS__ int process(dt_iop_module_t *self, const dt_dev_pixelpipe_t *pipe, const dt_dev_pixelpipe_iop_t *piece, const void *const restrict ivoid, void *const restrict ovoid)
Definition diffuse.c:1155
void commit_params(struct dt_iop_module_t *self, dt_iop_params_t *params, dt_dev_pixelpipe_t *pipe, dt_dev_pixelpipe_iop_t *piece)
Definition diffuse.c:134
dt_isotropy_t
Definition diffuse.c:143
@ DT_ISOTROPY_ISOPHOTE
Definition diffuse.c:145
@ DT_ISOTROPY_ISOTROPE
Definition diffuse.c:144
@ DT_ISOTROPY_GRADIENT
Definition diffuse.c:146
static float compute_anisotropy_factor(const float user_param)
Definition diffuse.c:955
const char * aliases()
Definition diffuse.c:169
static __DT_CLONE_TARGETS__ void build_mask(const float *const restrict input, uint8_t *const restrict mask, const float threshold, const size_t width, const size_t height)
Definition diffuse.c:1109
static __DT_CLONE_TARGETS__ void heat_PDE_diffusion(const float *const restrict high_freq, const float *const restrict low_freq, const uint8_t *const restrict mask, const int has_mask, float *const restrict output, const size_t width, const size_t height, const dt_aligned_pixel_simd_t anisotropy, const dt_isotropy_t isotropy_type[4], const float variance_threshold, const int mult, const float normalized_regularization, const dt_aligned_pixel_simd_t ABCD, const float strength, const int use_nontemporal)
Definition diffuse.c:760
const char * name()
Definition diffuse.c:164
void gui_init(struct dt_iop_module_t *self)
Definition diffuse.c:1647
void tiling_callback(struct dt_iop_module_t *self, const struct dt_dev_pixelpipe_t *pipe, const struct dt_dev_pixelpipe_iop_t *piece, struct dt_develop_tiling_t *tiling)
Definition diffuse.c:585
void cleanup_global(dt_iop_module_so_t *module)
Definition diffuse.c:1630
static cl_int wavelets_process_cl(const int devid, cl_mem in, cl_mem reconstructed, cl_mem mask, const size_t sizes[3], const int width, const int height, const dt_iop_diffuse_data_t *const data, dt_iop_diffuse_global_data_t *const gd, const float zoom, const int scales, const int has_mask, cl_mem HF[10], cl_mem LF_odd, cl_mem LF_even)
Definition diffuse.c:1261
int default_colorspace(dt_iop_module_t *self, dt_dev_pixelpipe_t *pipe, const dt_dev_pixelpipe_iop_t *piece)
Definition diffuse.c:197
int flags()
Definition diffuse.c:192
void init_presets(dt_iop_module_so_t *self)
Definition diffuse.c:296
#define H
Definition diffuse.c:623
#define MAX_NUM_SCALES
Definition diffuse.c:75
static __DT_CLONE_TARGETS__ void init_reconstruct(float *const restrict reconstructed, const size_t width, const size_t height)
Definition diffuse.c:612
#define KAPPA
Definition diffuse.c:624
static dt_isotropy_t check_isotropy_mode(const float anisotropy)
Definition diffuse.c:151
void init_global(dt_iop_module_so_t *module)
Definition diffuse.c:1611
int process_cl(struct dt_iop_module_t *self, const dt_dev_pixelpipe_t *pipe, const dt_dev_pixelpipe_iop_t *piece, cl_mem dev_in, cl_mem dev_out)
Definition diffuse.c:1486
int legacy_params(dt_iop_module_t *self, const void *const old_params, const int old_version, void *new_params, const int new_version)
Definition diffuse.c:202
static __DT_CLONE_TARGETS__ void inpaint_mask(float *const restrict inpainted, const float *const restrict original, const uint8_t *const restrict mask, const size_t width, const size_t height)
Definition diffuse.c:1122
static int dwt_interleave_rows(const int rowid, const int height, const int stride)
Definition dwt.h:93
void dt_gui_presets_add_generic(const char *name, dt_dev_operation_t op, const int32_t version, const void *params, const int32_t params_size, const int32_t enabled, const dt_develop_blend_colorspace_t blend_cst)
const char ** dt_iop_set_description(dt_iop_module_t *module, const char *main_text, const char *purpose, const char *input, const char *process, const char *output)
Definition imageop.c:1787
float dt_dev_get_module_scale(const dt_dev_pixelpipe_t *const pipe, const dt_iop_roi_t *const roi_in)
Definition imageop.c:133
@ IOP_FLAGS_INCLUDE_IN_STYLES
Definition imageop.h:185
@ IOP_FLAGS_SUPPORTS_BLENDING
Definition imageop.h:186
@ IOP_FLAGS_ALLOW_TILING
Definition imageop.h:188
@ IOP_GROUP_SHARPNESS
Definition imageop.h:160
GtkWidget * dt_bauhaus_slider_from_params(dt_iop_module_t *self, const char *param)
#define IOP_GUI_ALLOC(module)
Definition imageop_gui.h:93
void *const ovoid
static float kernel(const float *x, const float *y)
GtkWidget * dt_ui_section_label_new(const gchar *str)
Definition label.c:114
#define ABCD(A, B, C, D)
@ DT_DEBUG_OPENCL
Definition logging.h:57
void dt_print(dt_debug_thread_t thread, const char *msg,...) __attribute__((format(printf
Print to stdout when thread is enabled, prefixed with seconds since startup.
float *const restrict const size_t k
#define IS_NULL_PTR(p)
C is way too permissive with !=, == and if(var) checks, which can mean too many things depending on w...
Definition macros.h:96
#define DT_ALIGNED_PIXEL
Align a 4-float pixel on 16 bytes, enough for SSE. Same struct-member caveat as DT_ALIGNED_ARRAY,...
Definition mem_alloc.h:85
#define DT_ALIGNED_ARRAY
Align an object on a cacheline boundary, so AVX2 can load it whole.
Definition mem_alloc.h:80
#define DT_IS_ALIGNED(x)
Promise the compiler that x is already cacheline-aligned, and return it.
Definition mem_alloc.h:51
#define dt_free(ptr)
g_free() ptr and set it to NULL, skipping both if it is already NULL.
Definition mem_alloc.h:171
uint32_t width
Definition mipmap_cache.c:0
uint32_t height
Definition mipmap_cache.c:1
#define DT_MODULE_INTROSPECTION(MODVER, PARAMSTYPE)
DT_MODULE() for a module whose params struct is introspected.
int dt_opencl_local_buffer_opt(const int devid, const int kernel, dt_opencl_local_buffer_t *factors)
Definition opencl.c:3713
int dt_opencl_enqueue_kernel_2d(const int dev, const int kernel, const size_t *sizes)
Definition opencl.c:2554
void * dt_opencl_alloc_device(const int devid, const int width, const int height, const int bpp)
Definition opencl.c:2894
int dt_opencl_create_kernel(const int prog, const char *name)
Definition opencl.c:2448
void dt_opencl_free_kernel(const int kernel)
Definition opencl.c:2491
int dt_opencl_set_kernel_arg(const int dev, const int kernel, const int num, const size_t size, const void *arg)
Definition opencl.c:2545
int dt_opencl_enqueue_kernel_2d_with_local(const int dev, const int kernel, const size_t *sizes, const size_t *local)
Definition opencl.c:2560
void dt_opencl_release_mem_object(cl_mem mem)
Definition opencl.c:2805
#define ROUNDUP(a, n)
Definition opencl.h:82
#define ROUNDUPDHT(a, b)
Definition opencl.h:86
#define ROUNDUPDWD(a, b)
Definition opencl.h:85
#define dt_omploop_sfence()
Definition openmp.h:164
#define __OMP_DECLARE_SIMD__(...)
Definition openmp.h:100
#define __OMP_PARALLEL_FOR__(...)
Definition openmp.h:95
#define __OMP_PARALLEL_FOR_SIMD__(...)
Definition openmp.h:96
#define dt_pixelpipe_cache_alloc_align(size, pipe)
#define dt_pixelpipe_cache_free_align(mem)
#define dt_pixelpipe_cache_alloc_align_float(pixels, pipe)
#define dt_pixelpipe_cache_alloc_perthread_float(n, padded_size)
#define for_each_channel(_var,...)
Definition simd.h:87
float dt_aligned_pixel_simd_t __attribute__((vector_size(16), aligned(16)))
Apply one channel's tone curve to each of the three colour channels, or pass the channel through unto...
Definition simd.h:55
#define for_four_channels(_var,...)
Definition simd.h:89
const float uint32_t state[4]
unsigned __int64 uint64_t
Definition strptime.c:75
struct dt_iop_module_t *void * data
GtkWidget * anisotropy_third
Definition diffuse.c:114
GtkWidget * anisotropy_first
Definition diffuse.c:114
GtkWidget * regularization
Definition diffuse.c:113
GtkWidget * regularization_first
Definition diffuse.c:114
GtkWidget * anisotropy_second
Definition diffuse.c:114
GtkWidget * radius_center
Definition diffuse.c:113
GtkWidget * variance_threshold
Definition diffuse.c:114
GtkWidget * anisotropy_fourth
Definition diffuse.c:114
GtkWidget * widget
Definition imageop_gui.h:47
dt_dev_operation_t op
Definition imageop.h:235
dt_iop_global_data_t * data
Definition imageop.h:238
dt_iop_params_t * default_params
Definition imageop.h:333
struct dt_iop_module_gui_t * gui
Definition imageop.h:346
dt_iop_global_data_t * global_data
Definition imageop.h:337
int32_t params_size
Definition imageop.h:335
Region of interest passed through the pixelpipe.
Definition format.h:49
int width
Definition format.h:50
int height
Definition format.h:50
#define __DT_CLONE_TARGETS__
#define MIN(a, b)
Definition thinplate.c:32
#define MAX(a, b)
Definition thinplate.c:29
#define DT_GUI_BOX_SPACING