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
pl_generate_bayer_matrix(float * data,int size)34 void pl_generate_bayer_matrix(float *data, int size)
35 {
36 pl_assert(size >= 0);
37
38 // Start with a single entry of 0
39 data[0] = 0;
40
41 for (int sz = 1; sz < size; sz *= 2) {
42 // Make three copies of the current, appropriately shifted and scaled
43 for (int y = 0; y < sz; y ++) {
44 for (int x = 0; x < sz; x++) {
45 int offsets[] = {0, sz * size + sz, sz, sz * size};
46 int pos = y * size + x;
47
48 for (int i = 1; i < 4; i++)
49 data[pos + offsets[i]] = data[pos] + i / (4.0 * sz * sz);
50 }
51 }
52 }
53 }
54
55 #define MAX_SIZEB 8
56 #define MAX_SIZE (1 << MAX_SIZEB)
57 #define MAX_SIZE2 (MAX_SIZE * MAX_SIZE)
58
59 typedef uint_fast32_t index_t;
60
61 #define WRAP_SIZE2(k, x) ((index_t)((index_t)(x) & ((k)->size2 - 1)))
62 #define XY(k, x, y) ((index_t)(((x) | ((y) << (k)->sizeb))))
63
64 struct ctx {
65 unsigned int sizeb, size, size2;
66 unsigned int gauss_radius;
67 unsigned int gauss_middle;
68 uint64_t gauss[MAX_SIZE2];
69 index_t randomat[MAX_SIZE2];
70 bool calcmat[MAX_SIZE2];
71 uint64_t gaussmat[MAX_SIZE2];
72 index_t unimat[MAX_SIZE2];
73 };
74
makegauss(struct ctx * k,unsigned int sizeb)75 static void makegauss(struct ctx *k, unsigned int sizeb)
76 {
77 pl_assert(sizeb >= 1 && sizeb <= MAX_SIZEB);
78
79 k->sizeb = sizeb;
80 k->size = 1 << k->sizeb;
81 k->size2 = k->size * k->size;
82
83 k->gauss_radius = k->size / 2 - 1;
84 k->gauss_middle = XY(k, k->gauss_radius, k->gauss_radius);
85
86 unsigned int gauss_size = k->gauss_radius * 2 + 1;
87 unsigned int gauss_size2 = gauss_size * gauss_size;
88
89 for (index_t c = 0; c < k->size2; c++)
90 k->gauss[c] = 0;
91
92 double sigma = -log(1.5 / (double) UINT64_MAX * gauss_size2) / k->gauss_radius;
93
94 for (index_t gy = 0; gy <= k->gauss_radius; gy++) {
95 for (index_t gx = 0; gx <= gy; gx++) {
96 int cx = (int)gx - k->gauss_radius;
97 int cy = (int)gy - k->gauss_radius;
98 int sq = cx * cx + cy * cy;
99 double e = exp(-sqrt(sq) * sigma);
100 uint64_t v = e / gauss_size2 * (double) UINT64_MAX;
101 k->gauss[XY(k, gx, gy)] =
102 k->gauss[XY(k, gy, gx)] =
103 k->gauss[XY(k, gx, gauss_size - 1 - gy)] =
104 k->gauss[XY(k, gy, gauss_size - 1 - gx)] =
105 k->gauss[XY(k, gauss_size - 1 - gx, gy)] =
106 k->gauss[XY(k, gauss_size - 1 - gy, gx)] =
107 k->gauss[XY(k, gauss_size - 1 - gx, gauss_size - 1 - gy)] =
108 k->gauss[XY(k, gauss_size - 1 - gy, gauss_size - 1 - gx)] = v;
109 }
110 }
111
112 #ifndef NDEBUG
113 uint64_t total = 0;
114 for (index_t c = 0; c < k->size2; c++) {
115 uint64_t oldtotal = total;
116 total += k->gauss[c];
117 assert(total >= oldtotal);
118 }
119 #endif
120 }
121
setbit(struct ctx * k,index_t c)122 static void setbit(struct ctx *k, index_t c)
123 {
124 if (k->calcmat[c])
125 return;
126 k->calcmat[c] = true;
127 uint64_t *m = k->gaussmat;
128 uint64_t *me = k->gaussmat + k->size2;
129 uint64_t *g = k->gauss + WRAP_SIZE2(k, k->gauss_middle + k->size2 - c);
130 uint64_t *ge = k->gauss + k->size2;
131 while (g < ge)
132 *m++ += *g++;
133 g = k->gauss;
134 while (m < me)
135 *m++ += *g++;
136 }
137
getmin(struct ctx * k)138 static index_t getmin(struct ctx *k)
139 {
140 uint64_t min = UINT64_MAX;
141 index_t resnum = 0;
142 unsigned int size2 = k->size2;
143 for (index_t c = 0; c < size2; c++) {
144 if (k->calcmat[c])
145 continue;
146 uint64_t total = k->gaussmat[c];
147 if (total <= min) {
148 if (total != min) {
149 min = total;
150 resnum = 0;
151 }
152 k->randomat[resnum++] = c;
153 }
154 }
155 assert(resnum > 0);
156 if (resnum == 1)
157 return k->randomat[0];
158 if (resnum == size2)
159 return size2 / 2;
160 return k->randomat[rand() % resnum];
161 }
162
makeuniform(struct ctx * k)163 static void makeuniform(struct ctx *k)
164 {
165 unsigned int size2 = k->size2;
166 for (index_t c = 0; c < size2; c++) {
167 index_t r = getmin(k);
168 setbit(k, r);
169 k->unimat[r] = c;
170 }
171 }
172
pl_generate_blue_noise(float * data,int size)173 void pl_generate_blue_noise(float *data, int size)
174 {
175 pl_assert(size > 0);
176 int shift = PL_LOG2(size);
177
178 pl_assert((1 << shift) == size);
179 struct ctx *k = pl_zalloc_ptr(NULL, k);
180 makegauss(k, shift);
181 makeuniform(k);
182 float invscale = k->size2;
183 for(index_t y = 0; y < k->size; y++) {
184 for(index_t x = 0; x < k->size; x++)
185 data[x + y * k->size] = k->unimat[XY(k, x, y)] / invscale;
186 }
187 pl_free(k);
188 }
189