Skip to content

Instantly share code, notes, and snippets.

@kajott
Created June 26, 2019 10:49
Show Gist options
  • Save kajott/d9f9bb93043040bfe2f48f4f499903d8 to your computer and use it in GitHub Desktop.
Save kajott/d9f9bb93043040bfe2f48f4f499903d8 to your computer and use it in GitHub Desktop.
Blue Noise texture generation using the void-and-cluster algorithm
#if 0 // self-compiling code
gcc -std=c99 -Wall -Wextra -pedantic -Werror -g -O4 -march=native $0 -lm || exit 1
exec time ./a.out
#endif
// Blue Noise texture generation using the void-and-cluster algorithm
// implemented by Martin Fiedler <keyj@emphy.de>
// using an algorithm description written by Alan Wolfe:
// https://blog.demofox.org/2019/06/25/generating-blue-noise-textures-with-void-and-cluster/
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
// <parameter> log2 image size; 8 = 256x256 pixels
#define LOG2_SIZE 8
// <parameter> energy function falloff (Gaussian sigma value)
#define SIGMA 1.9
// <parameter> fraction of the image area to "seed" with initial random points
#define INITIAL_POINT_DENSITY 10
#define SIZE (1 << LOG2_SIZE) // image size (width = height)
#define AREA (SIZE * SIZE) // image area (wifth * height)
#define POINT_SHIFT ((2 * LOG2_SIZE) - 8) // shift amount between point ID and grayscale value
#define XY_MASK ((1 << LOG2_SIZE) - 1) // coordinate mask
#define I2X(i) ((i) & XY_MASK) // extract X coordinate from point index
#define I2Y(i) ((i) >> LOG2_SIZE) // extract Y coordinate from point index
#define XY2I(x,y) (((y) << LOG2_SIZE) + (x)) // build point index from coordinates
// structure representing the current processing state
struct map_state {
uint32_t bitmap[AREA / 32]; // dot occupancy bitmap
float emap[AREA]; // energy map
int imin, imax; // index of minimum and maximum energy spots
int points; // number of active points in the bitmap
} state1, state2, *state = &state1;
// get/set/clear individual bits in the current state bitmap
#define BMP_GET(i) ((state->bitmap[((i) >> 5)] >> ((i) & 31)) & 1u)
#define BMP_SET(i) do { state->bitmap[(i) >> 5] |= 1u << ((i) & 31); } while (0)
#define BMP_CLR(i) do { state->bitmap[(i) >> 5] &= ~(1u << ((i) & 31)); } while (0)
// energy distribution LUT
// (twice the image dimensions so that we can sample from a window)
float etab[(SIZE * 2) * (SIZE * 2)];
// TODO: Maybe it's faster to shrink the LUT to image size and handle wrap-
// around in update_map()? The large LUT might be cache-unfriendly ...
// resulting noise image
uint8_t noise[AREA];
// update the current energy map (and imax/imin) after adding (sign=+1)
// or removing (sign=-1) a point at index
void update_map(int index, float sign) {
int px = I2X(index);
int py = I2Y(index);
const float* pt = &etab[((SIZE - py) * (SIZE * 2)) + (SIZE - px)];
float* pm = state->emap;
float emin = 1e9f, emax = 0.0f;
int i = 0;
for (int y = 0; y < SIZE; ++y) {
for (int x = 0; x < SIZE; ++x) {
float e = (*pm++ += *pt++ * sign);
if (BMP_GET(i)) {
if (e > emax) { emax = e; state->imax = i; }
} else {
if (e < emin) { emin = e; state->imin = i; }
}
++i;
}
pt += SIZE;
}
}
int main(void) {
srand(time(NULL));
FILE *f;
// initialize the noise image to white (saves us 1/256th of the work, yay!)
memset(noise, 0xFF, AREA);
// prepare energy look-up table, using symmetries and repetitions:
// only create the upper-left 1/16th explicitly, then mirror it to get the
// upper-left quarter, then repeat that again
float *row = etab;
for (int y = 0; y < (SIZE / 2); ++y) {
for (int x = 0; x < (SIZE / 2); ++x) {
row[x] = row[SIZE - x] = (float)exp((x * x + y * y) * (-0.5 / (SIGMA * SIGMA)));
}
memcpy(&row[SIZE], row, SIZE * sizeof(float));
row += SIZE * 2;
}
for (int y = 0; y < (SIZE / 2); ++y) {
memcpy(&etab[(SIZE - y) * (SIZE * 2)], &etab[y * (SIZE * 2)], (SIZE * 2) * sizeof(float));
}
memcpy(&etab[SIZE * (SIZE * 2)], etab, SIZE * (SIZE * 2) * sizeof(float));
// debug: dump energy look-up table
#if 0
{
uint8_t conv[(SIZE * 2) * (SIZE * 2)];
for (int i = 0; i < ((SIZE * 2) * (SIZE * 2)); ++i) {
float f = etab[i];
for (int contrast = 10; contrast; contrast--) {
f = 1.0f - (1.0f - f) * (1.0f - f);
}
int v = (int)(f * 255.0 + 0.5);
conv[i] = (v > 255) ? 255 : (v < 0) ? 0 : (uint8_t) v;
}
f = fopen("bluenoise_energy_lut.pgm", "wb");
fprintf(f, "P5\n%d %d\n255\n", SIZE*2, SIZE*2);
fwrite(conv, SIZE*2, SIZE*2, f);
fclose(f);
}
#endif
// set initial points
while (state->points < (AREA / INITIAL_POINT_DENSITY)) {
int i = XY2I(rand() & XY_MASK, rand() & XY_MASK);
if (!BMP_GET(i)) {
BMP_SET(i); update_map(i, +1.0f);
++state->points;
}
}
printf("set %d initial points\n", state->points);
// re-distribute initial points
int redist = 0;
int last_point = AREA;
while (state->imax != last_point) {
BMP_CLR(state->imax); update_map(state->imax, -1.0f);
last_point = state->imin;
BMP_SET(state->imin); update_map(state->imin, +1.0f);
++redist;
}
printf("re-distributed initial points %d times\n", redist);
// serialize initial points
memcpy(&state2, &state1, sizeof(struct map_state));
while (state->points) {
noise[state->imax] = (--state->points) >> POINT_SHIFT;
BMP_CLR(state->imax); update_map(state->imax, -1.0f);
}
state = &state2;
// create points until half full
while (state->points < (AREA / 2)) {
noise[state->imin] = (state->points++) >> POINT_SHIFT;
BMP_SET(state->imin); update_map(state->imin, +1.0f);
}
printf("filled to %d points\n", state->points);
// invert energy map
for (int i = 0; i < (AREA / 32); ++i) {
state->bitmap[i] ^= (uint32_t)(-1);
}
memset(state->emap, 0, sizeof(state->emap));
for (int i = 0; i < AREA; ++i) {
if (BMP_GET(i)) { update_map(i, +1.0f); }
}
printf("inverted energy map\n");
// create points until full
while (state->points < (255 << POINT_SHIFT)) {
noise[state->imax] = (state->points++) >> POINT_SHIFT;
BMP_CLR(state->imax); update_map(state->imax, -1.0f);
}
printf("filled to %d points\n", state->points);
// debug: dump energy map
#if 0
{
uint8_t conv[(SIZE * 2) * SIZE];
float scale = 0.0f;
for (int i = 0; i < AREA; ++i) {
if (state->emap[i] > scale) { scale = state->emap[i]; }
}
printf("max energy in map: %g\n", scale);
scale = 255.0f / scale;
for (int i = 0; i < AREA; ++i) {
int v = (int)(state->emap[i] * scale);
conv[i] = (v > 255) ? 255 : (v < 0) ? 0 : (uint8_t) v;
}
f = fopen("bluenoise_energy_map.pgm", "wb");
fprintf(f, "P5\n%d %d\n255\n", SIZE, SIZE);
fwrite(conv, SIZE, SIZE, f);
fclose(f);
}
#endif
// write noise image
f = fopen("bluenoise.pgm", "wb");
fprintf(f, "P5\n%d %d\n255\n", SIZE, SIZE);
fwrite(noise, SIZE, SIZE, f);
fclose(f);
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment