Source code
Revision control
Copy as Markdown
Other Tools
/* Blue noise generation using the void-and-cluster method as described in
*
* The void-and-cluster method for dither array generation
* Ulichney, Robert A (1993)
*
*
* Note that running with openmp (-DUSE_OPENMP) will trigger additional
* randomness due to computing reductions in parallel, and is not recommended
* unless generating very large dither arrays.
*/
#include <assert.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include <stdio.h>
/* Booleans and utility functions */
#ifndef TRUE
# define TRUE 1
#endif
#ifndef FALSE
# define FALSE 0
#endif
typedef int bool_t;
int
imin (int x, int y)
{
return x < y ? x : y;
}
/* Memory allocation */
void *
malloc_abc (unsigned int a, unsigned int b, unsigned int c)
{
if (a >= INT32_MAX / b)
return NULL;
else if (a * b >= INT32_MAX / c)
return NULL;
else
return malloc (a * b * c);
}
/* Random number generation */
typedef uint32_t xorwow_state_t[5];
uint32_t
xorwow_next (xorwow_state_t *state)
{
uint32_t s = (*state)[0],
t = (*state)[3];
(*state)[3] = (*state)[2];
(*state)[2] = (*state)[1];
(*state)[1] = s;
t ^= t >> 2;
t ^= t << 1;
t ^= s ^ (s << 4);
(*state)[0] = t;
(*state)[4] += 362437;
return t + (*state)[4];
}
float
xorwow_float (xorwow_state_t *s)
{
return (xorwow_next (s) >> 9) / (float)((1 << 23) - 1);
}
/* Floating point matrices
*
* Used to cache the cluster sizes.
*/
typedef struct matrix_t {
int width;
int height;
float *buffer;
} matrix_t;
bool_t
matrix_init (matrix_t *matrix, int width, int height)
{
float *buffer;
if (!matrix)
return FALSE;
buffer = malloc_abc (width, height, sizeof (float));
if (!buffer)
return FALSE;
matrix->buffer = buffer;
matrix->width = width;
matrix->height = height;
return TRUE;
}
bool_t
matrix_copy (matrix_t *dst, matrix_t const *src)
{
float *srcbuf = src->buffer,
*srcend = src->buffer + src->width * src->height,
*dstbuf = dst->buffer;
if (dst->width != src->width || dst->height != src->height)
return FALSE;
while (srcbuf < srcend)
*dstbuf++ = *srcbuf++;
return TRUE;
}
float *
matrix_get (matrix_t *matrix, int x, int y)
{
return &matrix->buffer[y * matrix->width + x];
}
void
matrix_destroy (matrix_t *matrix)
{
free (matrix->buffer);
}
/* Binary patterns */
typedef struct pattern_t {
int width;
int height;
bool_t *buffer;
} pattern_t;
bool_t
pattern_init (pattern_t *pattern, int width, int height)
{
bool_t *buffer;
if (!pattern)
return FALSE;
buffer = malloc_abc (width, height, sizeof (bool_t));
if (!buffer)
return FALSE;
pattern->buffer = buffer;
pattern->width = width;
pattern->height = height;
return TRUE;
}
bool_t
pattern_copy (pattern_t *dst, pattern_t const *src)
{
bool_t *srcbuf = src->buffer,
*srcend = src->buffer + src->width * src->height,
*dstbuf = dst->buffer;
if (dst->width != src->width || dst->height != src->height)
return FALSE;
while (srcbuf < srcend)
*dstbuf++ = *srcbuf++;
return TRUE;
}
bool_t *
pattern_get (pattern_t *pattern, int x, int y)
{
return &pattern->buffer[y * pattern->width + x];
}
void
pattern_fill_white_noise (pattern_t *pattern, float fraction,
xorwow_state_t *s)
{
bool_t *buffer = pattern->buffer;
bool_t *end = buffer + (pattern->width * pattern->height);
while (buffer < end)
*buffer++ = xorwow_float (s) < fraction;
}
void
pattern_destroy (pattern_t *pattern)
{
free (pattern->buffer);
}
/* Dither arrays */
typedef struct array_t {
int width;
int height;
uint32_t *buffer;
} array_t;
bool_t
array_init (array_t *array, int width, int height)
{
uint32_t *buffer;
if (!array)
return FALSE;
buffer = malloc_abc (width, height, sizeof (uint32_t));
if (!buffer)
return FALSE;
array->buffer = buffer;
array->width = width;
array->height = height;
return TRUE;
}
uint32_t *
array_get (array_t *array, int x, int y)
{
return &array->buffer[y * array->width + x];
}
bool_t
array_save_ppm (array_t *array, const char *filename)
{
FILE *f = fopen(filename, "wb");
int i = 0;
int bpp = 2;
uint8_t buffer[1024];
if (!f)
return FALSE;
if (array->width * array->height - 1 < 256)
bpp = 1;
fprintf(f, "P5 %d %d %d\n", array->width, array->height,
array->width * array->height - 1);
while (i < array->width * array->height)
{
int j = 0;
for (; j < 1024 / bpp && j < array->width * array->height; ++j)
{
uint32_t v = array->buffer[i + j];
if (bpp == 2)
{
buffer[2 * j] = v & 0xff;
buffer[2 * j + 1] = (v & 0xff00) >> 8;
} else {
buffer[j] = v;
}
}
fwrite((void *)buffer, bpp, j, f);
i += j;
}
if (fclose(f) != 0)
return FALSE;
return TRUE;
}
bool_t
array_save (array_t *array, const char *filename)
{
int x, y;
FILE *f = fopen(filename, "wb");
if (!f)
return FALSE;
fprintf (f,
"/* WARNING: This file is generated by make-blue-noise.c\n"
" * Please edit that file instead of this one.\n"
" */\n"
"\n"
"#ifndef BLUE_NOISE_%dX%d_H\n"
"#define BLUE_NOISE_%dX%d_H\n"
"\n"
"#include <stdint.h>\n"
"\n", array->width, array->height, array->width, array->height);
fprintf (f, "static const uint16_t dither_blue_noise_%dx%d[%d] = {\n",
array->width, array->height, array->width * array->height);
for (y = 0; y < array->height; ++y)
{
fprintf (f, " ");
for (x = 0; x < array->width; ++x)
{
if (x != 0)
fprintf (f, ", ");
fprintf (f, "%d", *array_get (array, x, y));
}
fprintf (f, ",\n");
}
fprintf (f, "};\n");
fprintf (f, "\n#endif /* BLUE_NOISE_%dX%d_H */\n",
array->width, array->height);
if (fclose(f) != 0)
return FALSE;
return TRUE;
}
void
array_destroy (array_t *array)
{
free (array->buffer);
}
/* Dither array generation */
bool_t
compute_cluster_sizes (pattern_t *pattern, matrix_t *matrix)
{
int width = pattern->width,
height = pattern->height;
if (matrix->width != width || matrix->height != height)
return FALSE;
int px, py, qx, qy, dx, dy;
float tsqsi = 2.f * 1.5f * 1.5f;
#ifdef USE_OPENMP
#pragma omp parallel for default (none) \
private (py, px, qy, qx, dx, dy) \
shared (height, width, pattern, matrix, tsqsi)
#endif
for (py = 0; py < height; ++py)
{
for (px = 0; px < width; ++px)
{
bool_t pixel = *pattern_get (pattern, px, py);
float dist = 0.f;
for (qx = 0; qx < width; ++qx)
{
dx = imin (abs (qx - px), width - abs (qx - px));
dx = dx * dx;
for (qy = 0; qy < height; ++qy)
{
dy = imin (abs (qy - py), height - abs (qy - py));
dy = dy * dy;
dist += (pixel == *pattern_get (pattern, qx, qy))
* expf (- (dx + dy) / tsqsi);
}
}
*matrix_get (matrix, px, py) = dist;
}
}
return TRUE;
}
bool_t
swap_pixel (pattern_t *pattern, matrix_t *matrix, int x, int y)
{
int width = pattern->width,
height = pattern->height;
bool_t new;
float f,
dist = 0.f,
tsqsi = 2.f * 1.5f * 1.5f;
int px, py, dx, dy;
bool_t b;
new = !*pattern_get (pattern, x, y);
*pattern_get (pattern, x, y) = new;
if (matrix->width != width || matrix->height != height)
return FALSE;
#ifdef USE_OPENMP
#pragma omp parallel for reduction (+:dist) default (none) \
private (px, py, dx, dy, b, f) \
shared (x, y, width, height, pattern, matrix, new, tsqsi)
#endif
for (py = 0; py < height; ++py)
{
dy = imin (abs (py - y), height - abs (py - y));
dy = dy * dy;
for (px = 0; px < width; ++px)
{
dx = imin (abs (px - x), width - abs (px - x));
dx = dx * dx;
b = (*pattern_get (pattern, px, py) == new);
f = expf (- (dx + dy) / tsqsi);
*matrix_get (matrix, px, py) += (2 * b - 1) * f;
dist += b * f;
}
}
*matrix_get (matrix, x, y) = dist;
return TRUE;
}
void
largest_cluster (pattern_t *pattern, matrix_t *matrix,
bool_t pixel, int *xmax, int *ymax)
{
int width = pattern->width,
height = pattern->height;
int x, y;
float vmax = -INFINITY;
#ifdef USE_OPENMP
#pragma omp parallel default (none) \
private (x, y) \
shared (height, width, pattern, matrix, pixel, xmax, ymax, vmax)
#endif
{
int xbest = -1,
ybest = -1;
#ifdef USE_OPENMP
float vbest = -INFINITY;
#pragma omp for reduction (max: vmax) collapse (2)
#endif
for (y = 0; y < height; ++y)
{
for (x = 0; x < width; ++x)
{
if (*pattern_get (pattern, x, y) != pixel)
continue;
if (*matrix_get (matrix, x, y) > vmax)
{
vmax = *matrix_get (matrix, x, y);
#ifdef USE_OPENMP
vbest = vmax;
#endif
xbest = x;
ybest = y;
}
}
}
#ifdef USE_OPENMP
#pragma omp barrier
#pragma omp critical
{
if (vmax == vbest)
{
*xmax = xbest;
*ymax = ybest;
}
}
#else
*xmax = xbest;
*ymax = ybest;
#endif
}
assert (vmax > -INFINITY);
}
void
generate_initial_binary_pattern (pattern_t *pattern, matrix_t *matrix)
{
int xcluster = 0,
ycluster = 0,
xvoid = 0,
yvoid = 0;
for (;;)
{
largest_cluster (pattern, matrix, TRUE, &xcluster, &ycluster);
assert (*pattern_get (pattern, xcluster, ycluster) == TRUE);
swap_pixel (pattern, matrix, xcluster, ycluster);
largest_cluster (pattern, matrix, FALSE, &xvoid, &yvoid);
assert (*pattern_get (pattern, xvoid, yvoid) == FALSE);
swap_pixel (pattern, matrix, xvoid, yvoid);
if (xcluster == xvoid && ycluster == yvoid)
return;
}
}
bool_t
generate_dither_array (array_t *array,
pattern_t const *prototype, matrix_t const *matrix,
pattern_t *temp_pattern, matrix_t *temp_matrix)
{
int width = prototype->width,
height = prototype->height;
int x, y, rank;
int initial_rank = 0;
if (array->width != width || array->height != height)
return FALSE;
// Make copies of the prototype and associated sizes matrix since we will
// trash them
if (!pattern_copy (temp_pattern, prototype))
return FALSE;
if (!matrix_copy (temp_matrix, matrix))
return FALSE;
// Compute initial rank
for (y = 0; y < height; ++y)
{
for (x = 0; x < width; ++x)
{
if (*pattern_get (temp_pattern, x, y))
initial_rank += 1;
*array_get (array, x, y) = 0;
}
}
// Phase 1
for (rank = initial_rank; rank > 0; --rank)
{
largest_cluster (temp_pattern, temp_matrix, TRUE, &x, &y);
swap_pixel (temp_pattern, temp_matrix, x, y);
*array_get (array, x, y) = rank - 1;
}
// Make copies again for phases 2 & 3
if (!pattern_copy (temp_pattern, prototype))
return FALSE;
if (!matrix_copy (temp_matrix, matrix))
return FALSE;
// Phase 2 & 3
for (rank = initial_rank; rank < width * height; ++rank)
{
largest_cluster (temp_pattern, temp_matrix, FALSE, &x, &y);
swap_pixel (temp_pattern, temp_matrix, x, y);
*array_get (array, x, y) = rank;
}
return TRUE;
}
bool_t
generate (int size, xorwow_state_t *s,
char const *c_filename, char const *ppm_filename)
{
bool_t ok = TRUE;
pattern_t prototype, temp_pattern;
array_t array;
matrix_t matrix, temp_matrix;
printf ("Generating %dx%d blue noise...\n", size, size);
if (!pattern_init (&prototype, size, size))
return FALSE;
if (!pattern_init (&temp_pattern, size, size))
{
pattern_destroy (&prototype);
return FALSE;
}
if (!matrix_init (&matrix, size, size))
{
pattern_destroy (&temp_pattern);
pattern_destroy (&prototype);
return FALSE;
}
if (!matrix_init (&temp_matrix, size, size))
{
matrix_destroy (&matrix);
pattern_destroy (&temp_pattern);
pattern_destroy (&prototype);
return FALSE;
}
if (!array_init (&array, size, size))
{
matrix_destroy (&temp_matrix);
matrix_destroy (&matrix);
pattern_destroy (&temp_pattern);
pattern_destroy (&prototype);
return FALSE;
}
printf("Filling initial binary pattern with white noise...\n");
pattern_fill_white_noise (&prototype, .1, s);
printf("Initializing cluster sizes...\n");
if (!compute_cluster_sizes (&prototype, &matrix))
{
fprintf (stderr, "Error while computing cluster sizes\n");
ok = FALSE;
goto out;
}
printf("Generating initial binary pattern...\n");
generate_initial_binary_pattern (&prototype, &matrix);
printf("Generating dither array...\n");
if (!generate_dither_array (&array, &prototype, &matrix,
&temp_pattern, &temp_matrix))
{
fprintf (stderr, "Error while generating dither array\n");
ok = FALSE;
goto out;
}
printf("Saving dither array...\n");
if (!array_save (&array, c_filename))
{
fprintf (stderr, "Error saving dither array\n");
ok = FALSE;
goto out;
}
#if SAVE_PPM
if (!array_save_ppm (&array, ppm_filename))
{
fprintf (stderr, "Error saving dither array PPM\n");
ok = FALSE;
goto out;
}
#else
(void)ppm_filename;
#endif
printf("All done!\n");
out:
array_destroy (&array);
matrix_destroy (&temp_matrix);
matrix_destroy (&matrix);
pattern_destroy (&temp_pattern);
pattern_destroy (&prototype);
return ok;
}
int
main (void)
{
xorwow_state_t s = {1185956906, 12385940, 983948, 349208051, 901842};
if (!generate (64, &s, "blue-noise-64x64.h", "blue-noise-64x64.ppm"))
return -1;
return 0;
}