make-blue-noise.c (13556B)
1 /* Blue noise generation using the void-and-cluster method as described in 2 * 3 * The void-and-cluster method for dither array generation 4 * Ulichney, Robert A (1993) 5 * 6 * http://cv.ulichney.com/papers/1993-void-cluster.pdf 7 * 8 * Note that running with openmp (-DUSE_OPENMP) will trigger additional 9 * randomness due to computing reductions in parallel, and is not recommended 10 * unless generating very large dither arrays. 11 */ 12 13 #include <assert.h> 14 #include <stdlib.h> 15 #include <stdint.h> 16 #include <math.h> 17 #include <stdio.h> 18 19 /* Booleans and utility functions */ 20 21 #ifndef TRUE 22 # define TRUE 1 23 #endif 24 25 #ifndef FALSE 26 # define FALSE 0 27 #endif 28 29 typedef int bool_t; 30 31 int 32 imin (int x, int y) 33 { 34 return x < y ? x : y; 35 } 36 37 /* Memory allocation */ 38 void * 39 malloc_abc (unsigned int a, unsigned int b, unsigned int c) 40 { 41 if (a >= INT32_MAX / b) 42 return NULL; 43 else if (a * b >= INT32_MAX / c) 44 return NULL; 45 else 46 return malloc (a * b * c); 47 } 48 49 /* Random number generation */ 50 typedef uint32_t xorwow_state_t[5]; 51 52 uint32_t 53 xorwow_next (xorwow_state_t *state) 54 { 55 uint32_t s = (*state)[0], 56 t = (*state)[3]; 57 (*state)[3] = (*state)[2]; 58 (*state)[2] = (*state)[1]; 59 (*state)[1] = s; 60 61 t ^= t >> 2; 62 t ^= t << 1; 63 t ^= s ^ (s << 4); 64 65 (*state)[0] = t; 66 (*state)[4] += 362437; 67 68 return t + (*state)[4]; 69 } 70 71 float 72 xorwow_float (xorwow_state_t *s) 73 { 74 return (xorwow_next (s) >> 9) / (float)((1 << 23) - 1); 75 } 76 77 /* Floating point matrices 78 * 79 * Used to cache the cluster sizes. 80 */ 81 typedef struct matrix_t { 82 int width; 83 int height; 84 float *buffer; 85 } matrix_t; 86 87 bool_t 88 matrix_init (matrix_t *matrix, int width, int height) 89 { 90 float *buffer; 91 92 if (!matrix) 93 return FALSE; 94 95 buffer = malloc_abc (width, height, sizeof (float)); 96 97 if (!buffer) 98 return FALSE; 99 100 matrix->buffer = buffer; 101 matrix->width = width; 102 matrix->height = height; 103 104 return TRUE; 105 } 106 107 bool_t 108 matrix_copy (matrix_t *dst, matrix_t const *src) 109 { 110 float *srcbuf = src->buffer, 111 *srcend = src->buffer + src->width * src->height, 112 *dstbuf = dst->buffer; 113 114 if (dst->width != src->width || dst->height != src->height) 115 return FALSE; 116 117 while (srcbuf < srcend) 118 *dstbuf++ = *srcbuf++; 119 120 return TRUE; 121 } 122 123 float * 124 matrix_get (matrix_t *matrix, int x, int y) 125 { 126 return &matrix->buffer[y * matrix->width + x]; 127 } 128 129 void 130 matrix_destroy (matrix_t *matrix) 131 { 132 free (matrix->buffer); 133 } 134 135 /* Binary patterns */ 136 typedef struct pattern_t { 137 int width; 138 int height; 139 bool_t *buffer; 140 } pattern_t; 141 142 bool_t 143 pattern_init (pattern_t *pattern, int width, int height) 144 { 145 bool_t *buffer; 146 147 if (!pattern) 148 return FALSE; 149 150 buffer = malloc_abc (width, height, sizeof (bool_t)); 151 152 if (!buffer) 153 return FALSE; 154 155 pattern->buffer = buffer; 156 pattern->width = width; 157 pattern->height = height; 158 159 return TRUE; 160 } 161 162 bool_t 163 pattern_copy (pattern_t *dst, pattern_t const *src) 164 { 165 bool_t *srcbuf = src->buffer, 166 *srcend = src->buffer + src->width * src->height, 167 *dstbuf = dst->buffer; 168 169 if (dst->width != src->width || dst->height != src->height) 170 return FALSE; 171 172 while (srcbuf < srcend) 173 *dstbuf++ = *srcbuf++; 174 175 return TRUE; 176 } 177 178 bool_t * 179 pattern_get (pattern_t *pattern, int x, int y) 180 { 181 return &pattern->buffer[y * pattern->width + x]; 182 } 183 184 void 185 pattern_fill_white_noise (pattern_t *pattern, float fraction, 186 xorwow_state_t *s) 187 { 188 bool_t *buffer = pattern->buffer; 189 bool_t *end = buffer + (pattern->width * pattern->height); 190 191 while (buffer < end) 192 *buffer++ = xorwow_float (s) < fraction; 193 } 194 195 void 196 pattern_destroy (pattern_t *pattern) 197 { 198 free (pattern->buffer); 199 } 200 201 /* Dither arrays */ 202 typedef struct array_t { 203 int width; 204 int height; 205 uint32_t *buffer; 206 } array_t; 207 208 bool_t 209 array_init (array_t *array, int width, int height) 210 { 211 uint32_t *buffer; 212 213 if (!array) 214 return FALSE; 215 216 buffer = malloc_abc (width, height, sizeof (uint32_t)); 217 218 if (!buffer) 219 return FALSE; 220 221 array->buffer = buffer; 222 array->width = width; 223 array->height = height; 224 225 return TRUE; 226 } 227 228 uint32_t * 229 array_get (array_t *array, int x, int y) 230 { 231 return &array->buffer[y * array->width + x]; 232 } 233 234 bool_t 235 array_save_ppm (array_t *array, const char *filename) 236 { 237 FILE *f = fopen(filename, "wb"); 238 239 int i = 0; 240 int bpp = 2; 241 uint8_t buffer[1024]; 242 243 if (!f) 244 return FALSE; 245 246 if (array->width * array->height - 1 < 256) 247 bpp = 1; 248 249 fprintf(f, "P5 %d %d %d\n", array->width, array->height, 250 array->width * array->height - 1); 251 while (i < array->width * array->height) 252 { 253 int j = 0; 254 for (; j < 1024 / bpp && j < array->width * array->height; ++j) 255 { 256 uint32_t v = array->buffer[i + j]; 257 if (bpp == 2) 258 { 259 buffer[2 * j] = v & 0xff; 260 buffer[2 * j + 1] = (v & 0xff00) >> 8; 261 } else { 262 buffer[j] = v; 263 } 264 } 265 266 fwrite((void *)buffer, bpp, j, f); 267 i += j; 268 } 269 270 if (fclose(f) != 0) 271 return FALSE; 272 273 return TRUE; 274 } 275 276 bool_t 277 array_save (array_t *array, const char *filename) 278 { 279 int x, y; 280 FILE *f = fopen(filename, "wb"); 281 282 if (!f) 283 return FALSE; 284 285 fprintf (f, 286 "/* WARNING: This file is generated by make-blue-noise.c\n" 287 " * Please edit that file instead of this one.\n" 288 " */\n" 289 "\n" 290 "#ifndef BLUE_NOISE_%dX%d_H\n" 291 "#define BLUE_NOISE_%dX%d_H\n" 292 "\n" 293 "#include <stdint.h>\n" 294 "\n", array->width, array->height, array->width, array->height); 295 296 fprintf (f, "static const uint16_t dither_blue_noise_%dx%d[%d] = {\n", 297 array->width, array->height, array->width * array->height); 298 299 for (y = 0; y < array->height; ++y) 300 { 301 fprintf (f, " "); 302 for (x = 0; x < array->width; ++x) 303 { 304 if (x != 0) 305 fprintf (f, ", "); 306 307 fprintf (f, "%d", *array_get (array, x, y)); 308 } 309 310 fprintf (f, ",\n"); 311 } 312 fprintf (f, "};\n"); 313 314 fprintf (f, "\n#endif /* BLUE_NOISE_%dX%d_H */\n", 315 array->width, array->height); 316 317 if (fclose(f) != 0) 318 return FALSE; 319 320 return TRUE; 321 } 322 323 void 324 array_destroy (array_t *array) 325 { 326 free (array->buffer); 327 } 328 329 /* Dither array generation */ 330 bool_t 331 compute_cluster_sizes (pattern_t *pattern, matrix_t *matrix) 332 { 333 int width = pattern->width, 334 height = pattern->height; 335 336 if (matrix->width != width || matrix->height != height) 337 return FALSE; 338 339 int px, py, qx, qy, dx, dy; 340 float tsqsi = 2.f * 1.5f * 1.5f; 341 342 #ifdef USE_OPENMP 343 #pragma omp parallel for default (none) \ 344 private (py, px, qy, qx, dx, dy) \ 345 shared (height, width, pattern, matrix, tsqsi) 346 #endif 347 for (py = 0; py < height; ++py) 348 { 349 for (px = 0; px < width; ++px) 350 { 351 bool_t pixel = *pattern_get (pattern, px, py); 352 float dist = 0.f; 353 354 for (qx = 0; qx < width; ++qx) 355 { 356 dx = imin (abs (qx - px), width - abs (qx - px)); 357 dx = dx * dx; 358 359 for (qy = 0; qy < height; ++qy) 360 { 361 dy = imin (abs (qy - py), height - abs (qy - py)); 362 dy = dy * dy; 363 364 dist += (pixel == *pattern_get (pattern, qx, qy)) 365 * expf (- (dx + dy) / tsqsi); 366 } 367 } 368 369 *matrix_get (matrix, px, py) = dist; 370 } 371 } 372 373 return TRUE; 374 } 375 376 bool_t 377 swap_pixel (pattern_t *pattern, matrix_t *matrix, int x, int y) 378 { 379 int width = pattern->width, 380 height = pattern->height; 381 382 bool_t new; 383 384 float f, 385 dist = 0.f, 386 tsqsi = 2.f * 1.5f * 1.5f; 387 388 int px, py, dx, dy; 389 bool_t b; 390 391 new = !*pattern_get (pattern, x, y); 392 *pattern_get (pattern, x, y) = new; 393 394 if (matrix->width != width || matrix->height != height) 395 return FALSE; 396 397 398 #ifdef USE_OPENMP 399 #pragma omp parallel for reduction (+:dist) default (none) \ 400 private (px, py, dx, dy, b, f) \ 401 shared (x, y, width, height, pattern, matrix, new, tsqsi) 402 #endif 403 for (py = 0; py < height; ++py) 404 { 405 dy = imin (abs (py - y), height - abs (py - y)); 406 dy = dy * dy; 407 408 for (px = 0; px < width; ++px) 409 { 410 dx = imin (abs (px - x), width - abs (px - x)); 411 dx = dx * dx; 412 413 b = (*pattern_get (pattern, px, py) == new); 414 f = expf (- (dx + dy) / tsqsi); 415 *matrix_get (matrix, px, py) += (2 * b - 1) * f; 416 417 dist += b * f; 418 } 419 } 420 421 *matrix_get (matrix, x, y) = dist; 422 return TRUE; 423 } 424 425 void 426 largest_cluster (pattern_t *pattern, matrix_t *matrix, 427 bool_t pixel, int *xmax, int *ymax) 428 { 429 int width = pattern->width, 430 height = pattern->height; 431 432 int x, y; 433 434 float vmax = -INFINITY; 435 436 #ifdef USE_OPENMP 437 #pragma omp parallel default (none) \ 438 private (x, y) \ 439 shared (height, width, pattern, matrix, pixel, xmax, ymax, vmax) 440 #endif 441 { 442 int xbest = -1, 443 ybest = -1; 444 445 #ifdef USE_OPENMP 446 float vbest = -INFINITY; 447 448 #pragma omp for reduction (max: vmax) collapse (2) 449 #endif 450 for (y = 0; y < height; ++y) 451 { 452 for (x = 0; x < width; ++x) 453 { 454 if (*pattern_get (pattern, x, y) != pixel) 455 continue; 456 457 if (*matrix_get (matrix, x, y) > vmax) 458 { 459 vmax = *matrix_get (matrix, x, y); 460 #ifdef USE_OPENMP 461 vbest = vmax; 462 #endif 463 xbest = x; 464 ybest = y; 465 } 466 } 467 } 468 469 #ifdef USE_OPENMP 470 #pragma omp barrier 471 #pragma omp critical 472 { 473 if (vmax == vbest) 474 { 475 *xmax = xbest; 476 *ymax = ybest; 477 } 478 } 479 #else 480 *xmax = xbest; 481 *ymax = ybest; 482 #endif 483 } 484 485 assert (vmax > -INFINITY); 486 } 487 488 void 489 generate_initial_binary_pattern (pattern_t *pattern, matrix_t *matrix) 490 { 491 int xcluster = 0, 492 ycluster = 0, 493 xvoid = 0, 494 yvoid = 0; 495 496 for (;;) 497 { 498 largest_cluster (pattern, matrix, TRUE, &xcluster, &ycluster); 499 assert (*pattern_get (pattern, xcluster, ycluster) == TRUE); 500 swap_pixel (pattern, matrix, xcluster, ycluster); 501 502 largest_cluster (pattern, matrix, FALSE, &xvoid, &yvoid); 503 assert (*pattern_get (pattern, xvoid, yvoid) == FALSE); 504 swap_pixel (pattern, matrix, xvoid, yvoid); 505 506 if (xcluster == xvoid && ycluster == yvoid) 507 return; 508 } 509 } 510 511 bool_t 512 generate_dither_array (array_t *array, 513 pattern_t const *prototype, matrix_t const *matrix, 514 pattern_t *temp_pattern, matrix_t *temp_matrix) 515 { 516 int width = prototype->width, 517 height = prototype->height; 518 519 int x, y, rank; 520 521 int initial_rank = 0; 522 523 if (array->width != width || array->height != height) 524 return FALSE; 525 526 // Make copies of the prototype and associated sizes matrix since we will 527 // trash them 528 if (!pattern_copy (temp_pattern, prototype)) 529 return FALSE; 530 531 if (!matrix_copy (temp_matrix, matrix)) 532 return FALSE; 533 534 // Compute initial rank 535 for (y = 0; y < height; ++y) 536 { 537 for (x = 0; x < width; ++x) 538 { 539 if (*pattern_get (temp_pattern, x, y)) 540 initial_rank += 1; 541 542 *array_get (array, x, y) = 0; 543 } 544 } 545 546 // Phase 1 547 for (rank = initial_rank; rank > 0; --rank) 548 { 549 largest_cluster (temp_pattern, temp_matrix, TRUE, &x, &y); 550 swap_pixel (temp_pattern, temp_matrix, x, y); 551 *array_get (array, x, y) = rank - 1; 552 } 553 554 // Make copies again for phases 2 & 3 555 if (!pattern_copy (temp_pattern, prototype)) 556 return FALSE; 557 558 if (!matrix_copy (temp_matrix, matrix)) 559 return FALSE; 560 561 // Phase 2 & 3 562 for (rank = initial_rank; rank < width * height; ++rank) 563 { 564 largest_cluster (temp_pattern, temp_matrix, FALSE, &x, &y); 565 swap_pixel (temp_pattern, temp_matrix, x, y); 566 *array_get (array, x, y) = rank; 567 } 568 569 return TRUE; 570 } 571 572 bool_t 573 generate (int size, xorwow_state_t *s, 574 char const *c_filename, char const *ppm_filename) 575 { 576 bool_t ok = TRUE; 577 578 pattern_t prototype, temp_pattern; 579 array_t array; 580 matrix_t matrix, temp_matrix; 581 582 printf ("Generating %dx%d blue noise...\n", size, size); 583 584 if (!pattern_init (&prototype, size, size)) 585 return FALSE; 586 587 if (!pattern_init (&temp_pattern, size, size)) 588 { 589 pattern_destroy (&prototype); 590 return FALSE; 591 } 592 593 if (!matrix_init (&matrix, size, size)) 594 { 595 pattern_destroy (&temp_pattern); 596 pattern_destroy (&prototype); 597 return FALSE; 598 } 599 600 if (!matrix_init (&temp_matrix, size, size)) 601 { 602 matrix_destroy (&matrix); 603 pattern_destroy (&temp_pattern); 604 pattern_destroy (&prototype); 605 return FALSE; 606 } 607 608 if (!array_init (&array, size, size)) 609 { 610 matrix_destroy (&temp_matrix); 611 matrix_destroy (&matrix); 612 pattern_destroy (&temp_pattern); 613 pattern_destroy (&prototype); 614 return FALSE; 615 } 616 617 printf("Filling initial binary pattern with white noise...\n"); 618 pattern_fill_white_noise (&prototype, .1, s); 619 620 printf("Initializing cluster sizes...\n"); 621 if (!compute_cluster_sizes (&prototype, &matrix)) 622 { 623 fprintf (stderr, "Error while computing cluster sizes\n"); 624 ok = FALSE; 625 goto out; 626 } 627 628 printf("Generating initial binary pattern...\n"); 629 generate_initial_binary_pattern (&prototype, &matrix); 630 631 printf("Generating dither array...\n"); 632 if (!generate_dither_array (&array, &prototype, &matrix, 633 &temp_pattern, &temp_matrix)) 634 { 635 fprintf (stderr, "Error while generating dither array\n"); 636 ok = FALSE; 637 goto out; 638 } 639 640 printf("Saving dither array...\n"); 641 if (!array_save (&array, c_filename)) 642 { 643 fprintf (stderr, "Error saving dither array\n"); 644 ok = FALSE; 645 goto out; 646 } 647 648 #if SAVE_PPM 649 if (!array_save_ppm (&array, ppm_filename)) 650 { 651 fprintf (stderr, "Error saving dither array PPM\n"); 652 ok = FALSE; 653 goto out; 654 } 655 #else 656 (void)ppm_filename; 657 #endif 658 659 printf("All done!\n"); 660 661 out: 662 array_destroy (&array); 663 matrix_destroy (&temp_matrix); 664 matrix_destroy (&matrix); 665 pattern_destroy (&temp_pattern); 666 pattern_destroy (&prototype); 667 return ok; 668 } 669 670 int 671 main (void) 672 { 673 xorwow_state_t s = {1185956906, 12385940, 983948, 349208051, 901842}; 674 675 if (!generate (64, &s, "blue-noise-64x64.h", "blue-noise-64x64.ppm")) 676 return -1; 677 678 return 0; 679 }