Ansel 0.0
A darktable fork - bloat + design vision
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
fast_guided_filter.h
Go to the documentation of this file.
1/*
2 This file is part of darktable,
3 Copyright (C) 2019-2021 darktable developers.
4
5 darktable is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9
10 darktable is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with darktable. If not, see <http://www.gnu.org/licenses/>.
17*/
18
19#pragma once
20
21#include <assert.h>
22#include <math.h>
23#include <stdlib.h>
24#include <stdio.h>
25#include <string.h>
26#include <time.h>
27
28#include "common/box_filters.h"
29#include "common/darktable.h"
30#include "common/imagebuf.h"
31#include "control/control.h"
32
33
34/* NOTE: this code complies with the optimizations in "common/extra_optimizations.h".
35 * Consider including that at the beginning of a *.c file where you use this
36 * header (provided the rest of the code complies).
37 **/
38
39
40#define MIN_FLOAT exp2f(-16.0f)
41
42
48
49
50/***
51 * DOCUMENTATION
52 *
53 * Fast Iterative Guided filter for surface blur
54 *
55 * This is a fast vectorized implementation of guided filter for grey images optimized for
56 * the special case where the guiding and the guided image are the same, which is useful
57 * for edge-aware surface blur.
58 *
59 * Since the guided filter is a linear application, we can safely downscale
60 * the guiding and the guided image by a factor of 4, using a bilinear interpolation,
61 * compute the guidance at this scale, then upscale back to the original size
62 * and get a free 10x speed-up.
63 *
64 * Then, the vectorization adds another substantial speed-up. Overall, it brings a x50 to x200
65 * speed-up compared to the guided_filter.h lib. Of course, it requires every buffer to be
66 * 64-bits aligned.
67 *
68 * On top of the default guided filter, several pre- and post-processing options are provided :
69 *
70 * - mask quantization : perform a posterization of the guiding image in log2 space to
71 * help the guiding to produce smoother areas,
72 *
73 * - blending : perform a regular (linear) blending of a and b parameters after the
74 * variance analysis (aka the by-the-book guided filter), or a geometric mean of the filter output (by-the-book)
75 * and the original image, which produces a pleasing trade-off.
76 *
77 * - iterations : apply the guided filtering recursively, with kernel size increasing by sqrt(2)
78 * between each iteration, to diffuse the filter and soften edges transitions.
79 *
80 * Reference :
81 * Kaiming He, Jian Sun, Microsoft : https://arxiv.org/abs/1505.00996
82 **/
83
84
85 #ifdef _OPENMP
86#pragma omp declare simd
87#endif
89static inline float fast_clamp(const float value, const float bottom, const float top)
90{
91 // vectorizable clamping between bottom and top values
92 return fmaxf(fminf(value, top), bottom);
93}
94
95
97static inline void interpolate_bilinear(const float *const restrict in, const size_t width_in, const size_t height_in,
98 float *const restrict out, const size_t width_out, const size_t height_out,
99 const size_t ch)
100{
101 // Fast vectorized bilinear interpolation on ch channels
102#ifdef _OPENMP
103#pragma omp parallel for collapse(2) default(none) \
104 dt_omp_firstprivate(in, out, width_out, height_out, width_in, height_in, ch) \
105 schedule(simd:static)
106#endif
107 for(size_t i = 0; i < height_out; i++)
108 {
109 for(size_t j = 0; j < width_out; j++)
110 {
111 // Relative coordinates of the pixel in output space
112 const float x_out = (float)j /(float)width_out;
113 const float y_out = (float)i /(float)height_out;
114
115 // Corresponding absolute coordinates of the pixel in input space
116 const float x_in = x_out * (float)width_in;
117 const float y_in = y_out * (float)height_in;
118
119 // Nearest neighbours coordinates in input space
120 size_t x_prev = (size_t)floorf(x_in);
121 size_t x_next = x_prev + 1;
122 size_t y_prev = (size_t)floorf(y_in);
123 size_t y_next = y_prev + 1;
124
125 x_prev = (x_prev < width_in) ? x_prev : width_in - 1;
126 x_next = (x_next < width_in) ? x_next : width_in - 1;
127 y_prev = (y_prev < height_in) ? y_prev : height_in - 1;
128 y_next = (y_next < height_in) ? y_next : height_in - 1;
129
130 // Nearest pixels in input array (nodes in grid)
131 const size_t Y_prev = y_prev * width_in;
132 const size_t Y_next = y_next * width_in;
133 const float *const Q_NW = (float *)in + (Y_prev + x_prev) * ch;
134 const float *const Q_NE = (float *)in + (Y_prev + x_next) * ch;
135 const float *const Q_SE = (float *)in + (Y_next + x_next) * ch;
136 const float *const Q_SW = (float *)in + (Y_next + x_prev) * ch;
137
138 // Spatial differences between nodes
139 const float Dy_next = (float)y_next - y_in;
140 const float Dy_prev = 1.f - Dy_next; // because next - prev = 1
141 const float Dx_next = (float)x_next - x_in;
142 const float Dx_prev = 1.f - Dx_next; // because next - prev = 1
143
144 // Interpolate over ch layers
145 float *const pixel_out = (float *)out + (i * width_out + j) * ch;
146
147//#pragma unroll //LLVM warns it can't unroll -- presumably because 'ch' is not a constant
148 for(size_t c = 0; c < ch; c++)
149 {
150 pixel_out[c] = Dy_prev * (Q_SW[c] * Dx_next + Q_SE[c] * Dx_prev) +
151 Dy_next * (Q_NW[c] * Dx_next + Q_NE[c] * Dx_prev);
152 }
153 }
154 }
155}
156
157
159static inline void variance_analyse(const float *const restrict guide, // I
160 const float *const restrict mask, //p
161 float *const restrict ab,
162 const size_t width, const size_t height,
163 const int radius, const float feathering)
164{
165 // Compute a box average (filter) on a grey image over a window of size 2*radius + 1
166 // then get the variance of the guide and covariance with its mask
167 // output a and b, the linear blending params
168 // p, the mask is the quantised guide I
169
170 const size_t Ndim = width * height;
171 const size_t Ndimch = Ndim * 4;
172
173 /*
174 * input is array of struct : { { guide , mask, guide * guide, guide * mask } }
175 */
176 float *const restrict input = dt_alloc_align_float(Ndimch);
177 if(input == NULL) goto error;
178
179 // Pre-multiply guide and mask and pack all inputs into an array of 4x1 SIMD struct
180#ifdef _OPENMP
181#pragma omp parallel for default(none) \
182 dt_omp_firstprivate(guide, mask, Ndim, radius, input) \
183 schedule(simd:static)
184#endif
185 for(size_t k = 0; k < Ndim; k++)
186 {
187 const size_t index = k * 4;
188 input[index] = guide[k];
189 input[index + 1] = mask[k];
190 input[index + 2] = guide[k] * guide[k];
191 input[index + 3] = guide[k] * mask[k];
192 }
193
194 // blur the guide and mask as a four-channel image to exploit data locality and SIMD
195 dt_box_mean(input, height, width, 4, radius, 1);
196
197 // blend the result and store in output buffer
198#ifdef _OPENMP
199#pragma omp parallel for default(none) \
200 dt_omp_firstprivate(ab, input, width, height, feathering) \
201 schedule(static)
202#endif
203 for(size_t idx = 0; idx < width*height; idx++)
204 {
205 const float d = fmaxf((input[4*idx+2] - input[4*idx+0] * input[4*idx+0]) + feathering, 1e-15f); // avoid division by 0.
206 const float a = (input[4*idx+3] - input[4*idx+0] * input[4*idx+1]) / d;
207 const float b = input[4*idx+1] - a * input[4*idx+0];
208 ab[2*idx] = a;
209 ab[2*idx+1] = b;
210 }
211
212error:;
213 if(input) dt_free_align(input);
214}
215
216
218static inline void apply_linear_blending(float *const restrict image,
219 const float *const restrict ab,
220 const size_t num_elem)
221{
222#ifdef _OPENMP
223#pragma omp parallel for simd default(none) \
224dt_omp_firstprivate(image, ab, num_elem) \
225schedule(simd:static) aligned(image, ab:64)
226#endif
227 for(size_t k = 0; k < num_elem; k++)
228 {
229 // Note : image[k] is positive at the outside of the luminance mask
230 image[k] = fmaxf(image[k] * ab[k * 2] + ab[k * 2 + 1], MIN_FLOAT);
231 }
232}
233
234
236static inline void apply_linear_blending_w_geomean(float *const restrict image,
237 const float *const restrict ab,
238 const size_t num_elem)
239{
240#ifdef _OPENMP
241#pragma omp parallel for simd default(none) \
242dt_omp_firstprivate(image, ab, num_elem) \
243schedule(simd:static) aligned(image, ab:64)
244#endif
245 for(size_t k = 0; k < num_elem; k++)
246 {
247 // Note : image[k] is positive at the outside of the luminance mask
248 image[k] = sqrtf(image[k] * fmaxf(image[k] * ab[k * 2] + ab[k * 2 + 1], MIN_FLOAT));
249 }
250}
251
252
254static inline void quantize(const float *const restrict image,
255 float *const restrict out,
256 const size_t num_elem,
257 const float sampling, const float clip_min, const float clip_max)
258{
259 // Quantize in exposure levels evenly spaced in log by sampling
260
261 if(sampling == 0.0f)
262 {
263 // No-op
264 dt_iop_image_copy(out, image, num_elem);
265 }
266 else if(sampling == 1.0f)
267 {
268 // fast track
269#ifdef _OPENMP
270#pragma omp parallel for simd default(none) \
271dt_omp_firstprivate(image, out, num_elem, sampling, clip_min, clip_max) \
272schedule(simd:static) aligned(image, out:64)
273#endif
274 for(size_t k = 0; k < num_elem; k++)
275 out[k] = fast_clamp(exp2f(floorf(log2f(image[k]))), clip_min, clip_max);
276 }
277
278 else
279 {
280 // slow track
281#ifdef _OPENMP
282#pragma omp parallel for simd default(none) \
283dt_omp_firstprivate(image, out, num_elem, sampling, clip_min, clip_max) \
284schedule(simd:static) aligned(image, out:64)
285#endif
286 for(size_t k = 0; k < num_elem; k++)
287 out[k] = fast_clamp(exp2f(floorf(log2f(image[k]) / sampling) * sampling), clip_min, clip_max);
288 }
289}
290
291
293static inline void fast_surface_blur(float *const restrict image,
294 const size_t width, const size_t height,
295 const int radius, float feathering, const int iterations,
296 const dt_iop_guided_filter_blending_t filter, const float scale,
297 const float quantization, const float quantize_min, const float quantize_max)
298{
299 // Works in-place on a grey image
300
301 // A down-scaling of 4 seems empirically safe and consistent no matter the image zoom level
302 // see reference paper above for proof.
303 const float scaling = 4.0f;
304 const int ds_radius = (radius < 4) ? 1 : radius / scaling;
305
306 const size_t ds_height = height / scaling;
307 const size_t ds_width = width / scaling;
308
309 const size_t num_elem_ds = ds_width * ds_height;
310 const size_t num_elem = width * height;
311
312 float *const restrict ds_image = dt_alloc_sse_ps(dt_round_size_sse(num_elem_ds));
313 float *const restrict ds_mask = dt_alloc_sse_ps(dt_round_size_sse(num_elem_ds));
314 float *const restrict ds_ab = dt_alloc_sse_ps(dt_round_size_sse(num_elem_ds * 2));
315 float *const restrict ab = dt_alloc_sse_ps(dt_round_size_sse(num_elem * 2));
316
317 if(!ds_image || !ds_mask || !ds_ab || !ab)
318 {
319 dt_control_log(_("fast guided filter failed to allocate memory, check your RAM settings"));
320 goto clean;
321 }
322
323 // Downsample the image for speed-up
324 interpolate_bilinear(image, width, height, ds_image, ds_width, ds_height, 1);
325
326 // Iterations of filter models the diffusion, sort of
327 for(int i = 0; i < iterations; ++i)
328 {
329 // (Re)build the mask from the quantized image to help guiding
330 quantize(ds_image, ds_mask, ds_width * ds_height, quantization, quantize_min, quantize_max);
331
332 // Perform the patch-wise variance analyse to get
333 // the a and b parameters for the linear blending s.t. mask = a * I + b
334 variance_analyse(ds_mask, ds_image, ds_ab, ds_width, ds_height, ds_radius, feathering);
335
336 // Compute the patch-wise average of parameters a and b
337 dt_box_mean(ds_ab, ds_height, ds_width, 2, ds_radius, 1);
338
339 if(i != iterations - 1)
340 {
341 // Process the intermediate filtered image
342 apply_linear_blending(ds_image, ds_ab, num_elem_ds);
343 }
344 }
345
346 // Upsample the blending parameters a and b
347 interpolate_bilinear(ds_ab, ds_width, ds_height, ab, width, height, 2);
348
349 // Finally, blend the guided image
350 if(filter == DT_GF_BLENDING_LINEAR)
351 apply_linear_blending(image, ab, num_elem);
352 else if(filter == DT_GF_BLENDING_GEOMEAN)
353 apply_linear_blending_w_geomean(image, ab, num_elem);
354
355clean:
356 if(ab) dt_free_align(ab);
357 if(ds_ab) dt_free_align(ds_ab);
358 if(ds_mask) dt_free_align(ds_mask);
359 if(ds_image) dt_free_align(ds_image);
360}
361
362// clang-format off
363// modelines: These editor modelines have been set for all relevant files by tools/update_modelines.py
364// vim: shiftwidth=2 expandtab tabstop=2 cindent
365// kate: tab-indents: off; indent-width 2; replace-tabs on; indent-mode cstyle; remove-trailing-spaces modified;
366// clang-format on
static void error(char *msg)
Definition ashift_lsd.c:191
int width
Definition bilateral.h:1
int height
Definition bilateral.h:1
void dt_box_mean(float *const buf, const size_t height, const size_t width, const int ch, const int radius, const unsigned iterations)
Definition box_filters.c:1225
void dt_control_log(const char *msg,...)
Definition control.c:424
static float * dt_alloc_align_float(size_t pixels)
Definition darktable.h:345
static size_t dt_round_size_sse(const size_t size)
Definition darktable.h:285
#define __DT_CLONE_TARGETS__
Definition darktable.h:249
#define dt_free_align(A)
Definition darktable.h:334
static void * dt_alloc_sse_ps(size_t pixels)
Definition darktable.h:356
static __DT_CLONE_TARGETS__ void variance_analyse(const float *const restrict guide, const float *const restrict mask, float *const restrict ab, const size_t width, const size_t height, const int radius, const float feathering)
Definition fast_guided_filter.h:159
dt_iop_guided_filter_blending_t
Definition fast_guided_filter.h:44
@ DT_GF_BLENDING_LINEAR
Definition fast_guided_filter.h:45
@ DT_GF_BLENDING_GEOMEAN
Definition fast_guided_filter.h:46
static __DT_CLONE_TARGETS__ void quantize(const float *const restrict image, float *const restrict out, const size_t num_elem, const float sampling, const float clip_min, const float clip_max)
Definition fast_guided_filter.h:254
#define MIN_FLOAT
Definition fast_guided_filter.h:40
static __DT_CLONE_TARGETS__ void fast_surface_blur(float *const restrict image, const size_t width, const size_t height, const int radius, float feathering, const int iterations, const dt_iop_guided_filter_blending_t filter, const float scale, const float quantization, const float quantize_min, const float quantize_max)
Definition fast_guided_filter.h:293
static __DT_CLONE_TARGETS__ void apply_linear_blending_w_geomean(float *const restrict image, const float *const restrict ab, const size_t num_elem)
Definition fast_guided_filter.h:236
static __DT_CLONE_TARGETS__ void apply_linear_blending(float *const restrict image, const float *const restrict ab, const size_t num_elem)
Definition fast_guided_filter.h:218
static __DT_CLONE_TARGETS__ float fast_clamp(const float value, const float bottom, const float top)
Definition fast_guided_filter.h:89
static __DT_CLONE_TARGETS__ void interpolate_bilinear(const float *const restrict in, const size_t width_in, const size_t height_in, float *const restrict out, const size_t width_out, const size_t height_out, const size_t ch)
Definition fast_guided_filter.h:97
void dt_iop_image_copy(float *const __restrict__ out, const float *const __restrict__ in, const size_t nfloats)
Definition imagebuf.c:134