blind-gauss-blur.c (10542B)
1 /* See LICENSE file for copyright and license details. */ 2 #include "common.h" 3 4 USAGE("[-j jobs] [-s spread | -s 'auto'] [-acghvy] sd-stream") 5 6 static int chroma = 0; 7 static int noalpha = 0; 8 static int glow = 0; 9 static int vertical = 0; 10 static int horizontal = 0; 11 static int measure_y_only = 0; 12 static int auto_spread = 0; 13 static size_t jobs = 1; 14 static size_t spread = 0; 15 static void *original = NULL; 16 17 /* 18 * This is not a regular simple gaussian blur implementation. 19 * This implementation is able to apply different levels of 20 * blur on different pixels. It's therefore written a bit 21 * oddly. Instead of going through each pixel and calculate 22 * the new value for each pixel, it goes through each pixel 23 * and smears it out to the other pixels. 24 */ 25 26 #define BLUR_PIXEL_PROLOGUE(TYPE, DIR)\ 27 if (sig[i1][3] == 0)\ 28 goto no_blur_##DIR;\ 29 if (chroma || measure_y_only) {\ 30 k[0] = sig[i1][1] * sig[i1][3];\ 31 if (auto_spread)\ 32 spread = k[0] > 0 ? (size_t)(k[0] * 3 + (TYPE)0.5) : 0;\ 33 blur[2] = blur[0] = k[0] > 0;\ 34 c[0] = k[0] *= k[0] * 2, c[0] = sqrt(c[0] * (TYPE)M_PI);\ 35 k[0] = 1 / -k[0], c[0] = 1 / c[0];\ 36 if (chroma) {\ 37 k[2] = k[0];\ 38 c[2] = c[0];\ 39 c[1] = k[1] = 0;\ 40 blur[1] = 0;\ 41 } else {\ 42 k[2] = k[1] = k[0];\ 43 c[2] = c[1] = c[0];\ 44 blur[1] = blur[0];\ 45 }\ 46 } else {\ 47 if (auto_spread)\ 48 spread = 0;\ 49 for (i = 0; i < 3; i++) {\ 50 k[i] = sig[i1][i] * sig[i1][3];\ 51 if (auto_spread && k[i] > 0 && spread < (size_t)(k[i] * 3 + (TYPE)0.5))\ 52 spread = (size_t)(k[i] * 3 + (TYPE)0.5);\ 53 blur[i] = k[i] > 0;\ 54 c[i] = k[i] *= k[i] * 2, c[i] = sqrt(c[i] * (TYPE)M_PI);\ 55 k[i] = 1 / -k[i], c[i] = 1 / c[i];\ 56 }\ 57 }\ 58 if (blur[0] + blur[1] + blur[2] == 0)\ 59 goto no_blur_##DIR;\ 60 if (auto_spread && spread < 1)\ 61 spread = 1; 62 63 #define BLUR_PIXEL(TYPE, START, LOOP, DISTANCE)\ 64 if (k[0] == k[1] && k[1] == k[2]) {\ 65 START;\ 66 for (LOOP) {\ 67 d = (TYPE)(DISTANCE);\ 68 d *= d;\ 69 m = c[0] * exp(d * k[0]);\ 70 img[i2][0] += clr[i1][0] * m;\ 71 img[i2][1] += clr[i1][1] * m;\ 72 img[i2][2] += clr[i1][2] * m;\ 73 img[i2][3] += clr[i1][3] * m;\ 74 }\ 75 } else {\ 76 blurred = 0;\ 77 for (i = 0; i < 3; i++) {\ 78 if (blur[i])\ 79 blurred += 1;\ 80 else\ 81 img[i1][i] += clr[i1][i];\ 82 }\ 83 for (i = 0; i < 3; i++) {\ 84 if (!blur[i])\ 85 continue;\ 86 START;\ 87 for (LOOP) {\ 88 d = (TYPE)(DISTANCE);\ 89 d *= d;\ 90 m = c[i] * exp(d * k[i]);\ 91 img[i2][i] += clr[i1][i] * m;\ 92 img[i2][3] += clr[i1][3] * m / (TYPE)blurred;\ 93 }\ 94 }\ 95 } 96 97 #define BLUR_PIXEL_EPILOGUE(DIR)\ 98 continue;\ 99 no_blur_##DIR:\ 100 img[i1][0] = clr[i1][0];\ 101 img[i1][1] = clr[i1][1];\ 102 img[i1][2] = clr[i1][2];\ 103 img[i1][3] = clr[i1][3]; 104 105 #define BLUR(TYPE, DIR, SETSTART, SETEND, START, LOOP, DISTANCE)\ 106 do {\ 107 memset(img, 0, colour->frame_size);\ 108 start = 0, end = colour->height;\ 109 is_master = efork_jobs(&start, &end, jobs, &children);\ 110 i1 = start * colour->width;\ 111 for (y1 = start; y1 < end; y1++) {\ 112 for (x1 = 0; x1 < colour->width; x1++, i1++) {\ 113 BLUR_PIXEL_PROLOGUE(TYPE, DIR);\ 114 if (spread) {\ 115 SETSTART;\ 116 SETEND;\ 117 }\ 118 BLUR_PIXEL(TYPE, START, LOOP, DISTANCE);\ 119 BLUR_PIXEL_EPILOGUE(DIR);\ 120 }\ 121 }\ 122 ejoin_jobs(is_master, children);\ 123 } while (0) 124 125 #define PROCESS(TYPE)\ 126 do {\ 127 typedef TYPE pixel_t[4];\ 128 \ 129 pixel_t *restrict clr = (pixel_t *)cbuf;\ 130 pixel_t *restrict sig = (pixel_t *)sbuf;\ 131 pixel_t *restrict orig = original;\ 132 pixel_t *img = (pixel_t *)output;\ 133 pixel_t c, k;\ 134 size_t x1, y1, i1, x2, y2, i2;\ 135 TYPE d, m;\ 136 int i, blurred, blur[3] = {0, 0, 0};\ 137 size_t start, end, x2start, x2end, y2start, y2end;\ 138 int is_master;\ 139 pid_t *children;\ 140 \ 141 y2start = x2start = 0;\ 142 x2end = colour->width;\ 143 y2end = colour->height;\ 144 \ 145 if (glow)\ 146 memcpy(orig, clr, colour->frame_size);\ 147 if (chroma || !noalpha) {\ 148 start = 0, end = colour->height;\ 149 is_master = efork_jobs(&start, &end, jobs, &children);\ 150 \ 151 /* premultiply alpha channel */\ 152 if (!noalpha) {\ 153 i1 = start * colour->width;\ 154 for (y1 = start; y1 < end; y1++) {\ 155 for (x1 = 0; x1 < colour->width; x1++, i1++) {\ 156 clr[i1][0] *= clr[i1][3];\ 157 clr[i1][1] *= clr[i1][3];\ 158 clr[i1][2] *= clr[i1][3];\ 159 }\ 160 }\ 161 }\ 162 \ 163 /* store original image */\ 164 if (glow) {\ 165 i1 = start * colour->width;\ 166 memcpy(orig + i1, clr + i1, (end - start) * colour->width * sizeof(*clr));\ 167 }\ 168 \ 169 /* convert colour model */\ 170 if (chroma) {\ 171 i1 = start * colour->width;\ 172 for (y1 = start; y1 < end; y1++) {\ 173 for (x1 = 0; x1 < colour->width; x1++, i1++) {\ 174 clr[i1][0] = clr[i1][0] / (TYPE)D65_XYZ_X - clr[i1][1];\ 175 clr[i1][2] = clr[i1][2] / (TYPE)D65_XYZ_Z - clr[i1][1];\ 176 /* 177 * Explaination: 178 * Y is the luma and ((X / Xn - Y / Yn), (Z / Zn - Y / Yn)) 179 * is the chroma (according to CIELAB), where (Xn, Yn, Zn) 180 * is the white point. 181 */\ 182 }\ 183 }\ 184 }\ 185 /* Conversion makes no difference if blur is applied to all 186 * parameters: 187 * 188 * Gaussian blur: 189 * 190 * ∞ ∞ 191 * ⌠ ⌠ V(x,y) −((x−x₀)² + (y−y₀)²)/(2σ²) 192 * V′ (x₀,y₀) = │ │ ────── e dxdy 193 * σ ⌡ ⌡ 2πσ² 194 * −∞ −∞ 195 * 196 * With linear transformation, F: 197 * 198 * ∞ ∞ 199 * ⌠ ⌠ F(V(x,y)) −((x−x₀)² + (y−y₀)²)/(2σ²) 200 * V′ (x₀,y₀) = F⁻¹ │ │ ───────── e dxdy 201 * σ ⌡ ⌡ 2πσ² 202 * −∞ −∞ 203 * 204 * ∞ ∞ 205 * ⌠ ⌠ ⎛V(x,y) −((x−x₀)² + (y−y₀)²)/(2σ²)⎞ 206 * V′ (x₀,y₀) = F⁻¹ │ │ F⎜────── e ⎟ dxdy 207 * σ ⌡ ⌡ ⎝ 2πσ² ⎠ 208 * −∞ −∞ 209 * 210 * ∞ ∞ 211 * ⌠ ⌠ V(x,y) −((x−x₀)² + (y−y₀)²)/(2σ²) 212 * V′ (x₀,y₀) = (F⁻¹ ∘ F) │ │ ────── e dxdy 213 * σ ⌡ ⌡ 2πσ² 214 * −∞ −∞ 215 * 216 * ∞ ∞ 217 * ⌠ ⌠ V(x,y) −((x−x₀)² + (y−y₀)²)/(2σ²) 218 * V′ (x₀,y₀) = │ │ ────── e dxdy 219 * σ ⌡ ⌡ 2πσ² 220 * −∞ −∞ 221 * 222 * Just like expected, the colour space should not affect the 223 * result of guassian blur as long as it is linear. 224 */\ 225 \ 226 ejoin_jobs(is_master, children);\ 227 }\ 228 \ 229 /* blur */\ 230 if (horizontal)\ 231 BLUR(TYPE, horizontal,\ 232 x2start = spread > x1 ? 0 : x1 - spread,\ 233 x2end = spread + 1 > colour->width - x1 ? colour->width : x1 + spread + 1,\ 234 i2 = y1 * colour->width + x2start,\ 235 x2 = x2start; x2 < x2end; (x2++, i2++),\ 236 (ssize_t)x1 - (ssize_t)x2);\ 237 if (horizontal && vertical)\ 238 memcpy(clr, img, colour->frame_size);\ 239 if (vertical)\ 240 BLUR(TYPE, vertical,\ 241 y2start = spread > y1 ? 0 : y1 - spread,\ 242 y2end = spread + 1 > colour->height - y1 ? colour->height : y1 + spread + 1,\ 243 i2 = y2start * colour->width + x1,\ 244 y2 = y2start; y2 < y2end; (y2++, i2 += colour->width),\ 245 (ssize_t)y1 - (ssize_t)y2);\ 246 \ 247 start = 0, end = colour->height;\ 248 is_master = efork_jobs(&start, &end, jobs, &children);\ 249 \ 250 /* convert back to CIE XYZ */\ 251 if (chroma) {\ 252 i1 = start * colour->width;\ 253 for (y1 = start; y1 < end; y1++) {\ 254 for (x1 = 0; x1 < colour->width; x1++, i1++) {\ 255 img[i1][0] = (img[i1][0] + img[i1][1]) * (TYPE)D65_XYZ_X;\ 256 img[i1][2] = (img[i1][2] + img[i1][1]) * (TYPE)D65_XYZ_Z;\ 257 }\ 258 }\ 259 }\ 260 \ 261 /* apply glow */\ 262 if (glow) {\ 263 i1 = start * colour->width;\ 264 for (y1 = start; y1 < end; y1++) {\ 265 for (x1 = 0; x1 < colour->width; x1++, i1++) {\ 266 img[i1][0] += orig[i1][0];\ 267 img[i1][1] += orig[i1][1];\ 268 img[i1][2] += orig[i1][2];\ 269 img[i1][3] += orig[i1][3];\ 270 }\ 271 }\ 272 }\ 273 \ 274 /* unpremultiply alpha channel */\ 275 i1 = start * colour->width;\ 276 for (y1 = start; y1 < end; y1++) {\ 277 for (x1 = 0; x1 < colour->width; x1++, i1++) {\ 278 if (img[i1][3] != 0)\ 279 continue;\ 280 img[i1][0] /= img[i1][3];\ 281 img[i1][1] /= img[i1][3];\ 282 img[i1][2] /= img[i1][3];\ 283 }\ 284 }\ 285 \ 286 /* ensure the video if opaque if -a was used */\ 287 if (noalpha) {\ 288 i1 = start * colour->width;\ 289 for (y1 = start; y1 < end; y1++)\ 290 for (x1 = 0; x1 < colour->width; x1++, i1++)\ 291 img[i1][3] = 1;\ 292 }\ 293 \ 294 ejoin_jobs(is_master, children);\ 295 \ 296 (void) sigma;\ 297 } while (0) 298 299 static void 300 process_lf(char *restrict output, char *restrict cbuf, char *restrict sbuf, 301 struct stream *colour, struct stream *sigma) 302 { 303 PROCESS(double); 304 } 305 306 static void 307 process_f(char *restrict output, char *restrict cbuf, char *restrict sbuf, 308 struct stream *colour, struct stream *sigma) 309 { 310 PROCESS(float); 311 } 312 313 int 314 main(int argc, char *argv[]) 315 { 316 struct stream colour, sigma; 317 char *arg; 318 void (*process)(char *restrict output, char *restrict cbuf, char *restrict sbuf, 319 struct stream *colour, struct stream *sigma); 320 321 ARGBEGIN { 322 case 'a': 323 noalpha = 1; 324 break; 325 case 'c': 326 chroma = 1; 327 break; 328 case 'g': 329 glow = 1; 330 break; 331 case 'h': 332 horizontal = 1; 333 break; 334 case 'v': 335 vertical = 1; 336 break; 337 case 'y': 338 measure_y_only = 1; 339 break; 340 case 'j': 341 jobs = etozu_flag('j', UARGF(), 1, SHRT_MAX); 342 break; 343 case 's': 344 arg = UARGF(); 345 if (!strcmp(arg, "auto")) 346 auto_spread = 1; 347 else 348 spread = etozu_flag('s', arg, 1, SIZE_MAX); 349 break; 350 default: 351 usage(); 352 } ARGEND; 353 354 if (argc != 1) 355 usage(); 356 357 if (!vertical && !horizontal) 358 vertical = horizontal = 1; 359 360 eopen_stream(&colour, NULL); 361 eopen_stream(&sigma, argv[0]); 362 363 SELECT_PROCESS_FUNCTION(&colour); 364 CHECK_CHANS(&colour, == 3, == (measure_y_only ? 1 : colour.luma_chan)); 365 CHECK_N_CHAN(&colour, 4, 4); 366 367 echeck_compat(&colour, &sigma); 368 369 if (jobs > colour.height) 370 jobs = colour.height; 371 372 if (glow) 373 original = emalloc(colour.frame_size); 374 375 fprint_stream_head(stdout, &colour); 376 efflush(stdout, "<stdout>"); 377 process_each_frame_two_streams(&colour, &sigma, STDOUT_FILENO, "<stdout>", process); 378 free(original); 379 return 0; 380 }