Quelle  make-blue-noise.c

/* 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)
 *
 *     http://cv.ulichney.com/papers/1993-void-cluster.pdf
 *
 * 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;
}