Branch data Line data Source code
1 : : /*
2 : : * Generate a noise texture for dithering images.
3 : : * Copyright © 2013 Wessel Dankers <wsl@fruit.je>
4 : : *
5 : : * This file is part of libplacebo.
6 : : *
7 : : * libplacebo is free software; you can redistribute it and/or
8 : : * modify it under the terms of the GNU Lesser General Public
9 : : * License as published by the Free Software Foundation; either
10 : : * version 2.1 of the License, or (at your option) any later version.
11 : : *
12 : : * libplacebo is distributed in the hope that it will be useful,
13 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 : : * GNU Lesser General Public License for more details.
16 : : *
17 : : * You should have received a copy of the GNU Lesser General Public
18 : : * License along with libplacebo. If not, see <http://www.gnu.org/licenses/>.
19 : : *
20 : : * The original code is taken from mpv, under the same license.
21 : : */
22 : :
23 : :
24 : : #include <stdint.h>
25 : : #include <stdbool.h>
26 : : #include <stdlib.h>
27 : : #include <inttypes.h>
28 : : #include <string.h>
29 : : #include <assert.h>
30 : : #include <math.h>
31 : :
32 : : #include "common.h"
33 : :
34 : : #include <libplacebo/dither.h>
35 : :
36 : 5 : void pl_generate_bayer_matrix(float *data, int size)
37 : : {
38 [ - + ]: 5 : pl_assert(size >= 0);
39 : :
40 : : // Start with a single entry of 0
41 : 5 : data[0] = 0;
42 : :
43 [ + + ]: 33 : for (int sz = 1; sz < size; sz *= 2) {
44 : : // Make three copies of the current, appropriately shifted and scaled
45 [ + + ]: 295 : for (int y = 0; y < sz; y ++) {
46 [ + + ]: 5812 : for (int x = 0; x < sz; x++) {
47 : 5545 : int offsets[] = {0, sz * size + sz, sz, sz * size};
48 : 5545 : int pos = y * size + x;
49 : :
50 [ + + ]: 22180 : for (int i = 1; i < 4; i++)
51 : 16635 : data[pos + offsets[i]] = data[pos] + i / (4.0 * sz * sz);
52 : : }
53 : : }
54 : : }
55 : 5 : }
56 : :
57 : : #define MAX_SIZEB 8
58 : : #define MAX_SIZE (1 << MAX_SIZEB)
59 : : #define MAX_SIZE2 (MAX_SIZE * MAX_SIZE)
60 : :
61 : : typedef uint_fast32_t index_t;
62 : :
63 : : #define WRAP_SIZE2(k, x) ((index_t)((index_t)(x) & ((k)->size2 - 1)))
64 : : #define XY(k, x, y) ((index_t)(((x) | ((y) << (k)->sizeb))))
65 : :
66 : : struct ctx {
67 : : unsigned int sizeb, size, size2;
68 : : unsigned int gauss_radius;
69 : : unsigned int gauss_middle;
70 : : uint64_t gauss[MAX_SIZE2];
71 : : index_t randomat[MAX_SIZE2];
72 : : bool calcmat[MAX_SIZE2];
73 : : uint64_t gaussmat[MAX_SIZE2];
74 : : index_t unimat[MAX_SIZE2];
75 : : };
76 : :
77 : 9 : static void makegauss(struct ctx *k, unsigned int sizeb)
78 : : {
79 [ - + ]: 9 : pl_assert(sizeb >= 1 && sizeb <= MAX_SIZEB);
80 : :
81 : 9 : k->sizeb = sizeb;
82 : 9 : k->size = 1 << k->sizeb;
83 : 9 : k->size2 = k->size * k->size;
84 : :
85 : 9 : k->gauss_radius = k->size / 2 - 1;
86 : 9 : k->gauss_middle = XY(k, k->gauss_radius, k->gauss_radius);
87 : :
88 : 9 : unsigned int gauss_size = k->gauss_radius * 2 + 1;
89 : 9 : unsigned int gauss_size2 = gauss_size * gauss_size;
90 : :
91 [ + + ]: 32793 : for (index_t c = 0; c < k->size2; c++)
92 : 32784 : k->gauss[c] = 0;
93 : :
94 : 9 : double sigma = -log(1.5 / (double) UINT64_MAX * gauss_size2) / k->gauss_radius;
95 : :
96 [ + + ]: 267 : for (index_t gy = 0; gy <= k->gauss_radius; gy++) {
97 [ + + ]: 4485 : for (index_t gx = 0; gx <= gy; gx++) {
98 : 4227 : int cx = (int)gx - k->gauss_radius;
99 : 4227 : int cy = (int)gy - k->gauss_radius;
100 : 4227 : int sq = cx * cx + cy * cy;
101 : 4227 : double e = exp(-sqrt(sq) * sigma);
102 : 4227 : uint64_t v = e / gauss_size2 * (double) UINT64_MAX;
103 : 4227 : k->gauss[XY(k, gx, gy)] =
104 : 4227 : k->gauss[XY(k, gy, gx)] =
105 : 4227 : k->gauss[XY(k, gx, gauss_size - 1 - gy)] =
106 : 4227 : k->gauss[XY(k, gy, gauss_size - 1 - gx)] =
107 : 4227 : k->gauss[XY(k, gauss_size - 1 - gx, gy)] =
108 : 4227 : k->gauss[XY(k, gauss_size - 1 - gy, gx)] =
109 : 4227 : k->gauss[XY(k, gauss_size - 1 - gx, gauss_size - 1 - gy)] =
110 : 4227 : k->gauss[XY(k, gauss_size - 1 - gy, gauss_size - 1 - gx)] = v;
111 : : }
112 : : }
113 : :
114 : : #ifndef NDEBUG
115 : : uint64_t total = 0;
116 [ + + ]: 32793 : for (index_t c = 0; c < k->size2; c++) {
117 : : uint64_t oldtotal = total;
118 : 32784 : total += k->gauss[c];
119 [ - + ]: 32784 : assert(total >= oldtotal);
120 : : }
121 : : #endif
122 : 9 : }
123 : :
124 : 32784 : static void setbit(struct ctx *k, index_t c)
125 : : {
126 [ + - ]: 32784 : if (k->calcmat[c])
127 : : return;
128 : 32784 : k->calcmat[c] = true;
129 : 32784 : uint64_t *m = k->gaussmat;
130 : 32784 : uint64_t *me = k->gaussmat + k->size2;
131 : 32784 : uint64_t *g = k->gauss + WRAP_SIZE2(k, k->gauss_middle + k->size2 - c);
132 : 32784 : uint64_t *ge = k->gauss + k->size2;
133 [ + + ]: 67158168 : while (g < ge)
134 : 67125384 : *m++ += *g++;
135 : : g = k->gauss;
136 [ + + ]: 67125384 : while (m < me)
137 : 67092600 : *m++ += *g++;
138 : : }
139 : :
140 : 32784 : static index_t getmin(struct ctx *k)
141 : : {
142 : : uint64_t min = UINT64_MAX;
143 : : index_t resnum = 0;
144 : 32784 : unsigned int size2 = k->size2;
145 [ + + ]: 134250768 : for (index_t c = 0; c < size2; c++) {
146 [ + + ]: 134217984 : if (k->calcmat[c])
147 : 67092600 : continue;
148 : 67125384 : uint64_t total = k->gaussmat[c];
149 [ + + ]: 67125384 : if (total <= min) {
150 [ + + ]: 321144 : if (total != min) {
151 : : min = total;
152 : : resnum = 0;
153 : : }
154 : 321144 : k->randomat[resnum++] = c;
155 : : }
156 : : }
157 [ - + ]: 32784 : assert(resnum > 0);
158 [ + + ]: 32784 : if (resnum == 1)
159 : 32727 : return k->randomat[0];
160 [ + + ]: 57 : if (resnum == size2)
161 : 9 : return size2 / 2;
162 : 48 : return k->randomat[rand() % resnum];
163 : : }
164 : :
165 : 9 : static void makeuniform(struct ctx *k)
166 : : {
167 : 9 : unsigned int size2 = k->size2;
168 [ + + ]: 32793 : for (index_t c = 0; c < size2; c++) {
169 : 32784 : index_t r = getmin(k);
170 : 32784 : setbit(k, r);
171 : 32784 : k->unimat[r] = c;
172 : : }
173 : 9 : }
174 : :
175 : 9 : void pl_generate_blue_noise(float *data, int size)
176 : : {
177 [ - + ]: 9 : pl_assert(size > 0);
178 : 9 : int shift = PL_LOG2(size);
179 : :
180 [ - + ]: 9 : pl_assert((1 << shift) == size);
181 : 9 : struct ctx *k = pl_zalloc_ptr(NULL, k);
182 : 9 : makegauss(k, shift);
183 : 9 : makeuniform(k);
184 : 9 : float invscale = k->size2;
185 [ + + ]: 525 : for(index_t y = 0; y < k->size; y++) {
186 [ + + ]: 33300 : for(index_t x = 0; x < k->size; x++)
187 : 32784 : data[x + y * k->size] = k->unimat[XY(k, x, y)] / invscale;
188 : : }
189 : 9 : pl_free(k);
190 : 9 : }
191 : :
192 : : const struct pl_error_diffusion_kernel pl_error_diffusion_simple = {
193 : : .name = "simple",
194 : : .description = "Simple error diffusion",
195 : : .shift = 1,
196 : : .pattern = {{0, 0, 0, 1, 0},
197 : : {0, 0, 1, 0, 0},
198 : : {0, 0, 0, 0, 0}},
199 : : .divisor = 2,
200 : : };
201 : :
202 : : const struct pl_error_diffusion_kernel pl_error_diffusion_false_fs = {
203 : : .name = "false-fs",
204 : : .description = "False Floyd-Steinberg kernel",
205 : : .shift = 1,
206 : : .pattern = {{0, 0, 0, 3, 0},
207 : : {0, 0, 3, 2, 0},
208 : : {0, 0, 0, 0, 0}},
209 : : .divisor = 8,
210 : : };
211 : :
212 : : const struct pl_error_diffusion_kernel pl_error_diffusion_sierra_lite = {
213 : : .name = "sierra-lite",
214 : : .description = "Sierra Lite kernel",
215 : : .shift = 2,
216 : : .pattern = {{0, 0, 0, 2, 0},
217 : : {0, 1, 1, 0, 0},
218 : : {0, 0, 0, 0, 0}},
219 : : .divisor = 4,
220 : : };
221 : :
222 : : const struct pl_error_diffusion_kernel pl_error_diffusion_floyd_steinberg = {
223 : : .name = "floyd-steinberg",
224 : : .description = "Floyd Steinberg kernel",
225 : : .shift = 2,
226 : : .pattern = {{0, 0, 0, 7, 0},
227 : : {0, 3, 5, 1, 0},
228 : : {0, 0, 0, 0, 0}},
229 : : .divisor = 16,
230 : : };
231 : :
232 : : const struct pl_error_diffusion_kernel pl_error_diffusion_atkinson = {
233 : : .name = "atkinson",
234 : : .description = "Atkinson kernel",
235 : : .shift = 2,
236 : : .pattern = {{0, 0, 0, 1, 1},
237 : : {0, 1, 1, 1, 0},
238 : : {0, 0, 1, 0, 0}},
239 : : .divisor = 8,
240 : : };
241 : :
242 : : const struct pl_error_diffusion_kernel pl_error_diffusion_jarvis_judice_ninke = {
243 : : .name = "jarvis-judice-ninke",
244 : : .description = "Jarvis, Judice & Ninke kernel",
245 : : .shift = 3,
246 : : .pattern = {{0, 0, 0, 7, 5},
247 : : {3, 5, 7, 5, 3},
248 : : {1, 3, 5, 3, 1}},
249 : : .divisor = 48,
250 : : };
251 : :
252 : : const struct pl_error_diffusion_kernel pl_error_diffusion_stucki = {
253 : : .name = "stucki",
254 : : .description = "Stucki kernel",
255 : : .shift = 3,
256 : : .pattern = {{0, 0, 0, 8, 4},
257 : : {2, 4, 8, 4, 2},
258 : : {1, 2, 4, 2, 1}},
259 : : .divisor = 42,
260 : : };
261 : :
262 : : const struct pl_error_diffusion_kernel pl_error_diffusion_burkes = {
263 : : .name = "burkes",
264 : : .description = "Burkes kernel",
265 : : .shift = 3,
266 : : .pattern = {{0, 0, 0, 8, 4},
267 : : {2, 4, 8, 4, 2},
268 : : {0, 0, 0, 0, 0}},
269 : : .divisor = 32,
270 : : };
271 : :
272 : : const struct pl_error_diffusion_kernel pl_error_diffusion_sierra2 = {
273 : : .name = "sierra-2",
274 : : .description = "Two-row Sierra",
275 : : .shift = 3,
276 : : .pattern = {{0, 0, 0, 4, 3},
277 : : {1, 2, 3, 2, 1},
278 : : {0, 0, 0, 0, 0}},
279 : : .divisor = 16,
280 : : };
281 : :
282 : : const struct pl_error_diffusion_kernel pl_error_diffusion_sierra3 = {
283 : : .name = "sierra-3",
284 : : .description = "Three-row Sierra",
285 : : .shift = 3,
286 : : .pattern = {{0, 0, 0, 5, 3},
287 : : {2, 4, 5, 4, 2},
288 : : {0, 2, 3, 2, 0}},
289 : : .divisor = 32,
290 : : };
291 : :
292 : : const struct pl_error_diffusion_kernel * const pl_error_diffusion_kernels[] = {
293 : : &pl_error_diffusion_simple,
294 : : &pl_error_diffusion_false_fs,
295 : : &pl_error_diffusion_sierra_lite,
296 : : &pl_error_diffusion_floyd_steinberg,
297 : : &pl_error_diffusion_atkinson,
298 : : &pl_error_diffusion_jarvis_judice_ninke,
299 : : &pl_error_diffusion_stucki,
300 : : &pl_error_diffusion_burkes,
301 : : &pl_error_diffusion_sierra2,
302 : : &pl_error_diffusion_sierra3,
303 : : NULL
304 : : };
305 : :
306 : : const int pl_num_error_diffusion_kernels = PL_ARRAY_SIZE(pl_error_diffusion_kernels) - 1;
307 : :
308 : : // Find the error diffusion kernel with the given name, or NULL on failure.
309 : 0 : const struct pl_error_diffusion_kernel *pl_find_error_diffusion_kernel(const char *name)
310 : : {
311 [ # # ]: 0 : for (int i = 0; i < pl_num_error_diffusion_kernels; i++) {
312 [ # # ]: 0 : if (strcmp(name, pl_error_diffusion_kernels[i]->name) == 0)
313 : 0 : return pl_error_diffusion_kernels[i];
314 : : }
315 : :
316 : : return NULL;
317 : : }
|