Generating Noise
In my previous post I talked about generating voronoi noise, but when the subject of noise comes up, most people tend to think of perlin noise or simplex noise rather than anything else.
However, one doesn’t have to go that fancy to get relatively decent results. But, without any further ado, let’s dive into it and explore the world of noise.
To start things off on the right foot, let’s generate some good old static noise.
static void process(
rgba_t *restrict output,
const int w,
const int h
)
{
int i, n, c;
for(i = 0, n = w * h; i < n; i++)
{
c = rand() & 0xFF;
output[i] = (rgba_t) {
.r = c,
.g = c,
.b = c,
.a = 0xFF
};
}
}
Do not worry, there’s no poltergeist coming through if you stare at it for too long. I pinky promise.
Let’s see what happens if we blur this by simply taking the average of the surrounding pixels.
static void randomize(
rgba_t *restrict output,
const int w,
const int h,
const unsigned int seed
)
{
int i, n, c;
srand(seed);
for(i = 0, n = w * h; i < n; i++)
{
c = rand() & 0xFF;
output[i] = (rgba_t) {
.r = c,
.g = c,
.b = c,
.a = 0xFF
};
}
}
static void process(
rgba_t *restrict output,
const int w,
const int h
)
{
int x, y, yo, you, yod, ww, hh;
uint32_t c;
rgba_t cu, cd, cl, cr;
ww = w - 1;
hh = h - 1;
for(y = 0; y < h; y++)
{
yo = y * w;
you = max(y - 1, 0) * w;
yod = min(y + 1, hh) * w;
for(x = 0; x < w; x++)
{
cu = output[x + you];
cd = output[x + yod];
cl = output[max(x - 1, 0) + yo];
cr = output[min(x + 1, ww) + yo];
c = (cu.r + cd.r + cl.r + cr.r) >> 2;
output[x + yo] = (rgba_t) {
.r = c,
.g = c,
.b = c,
.a = 0xFF
};
}
}
}
This works best with a relatively small size and at least 4 passes of blur.
w = 128;
h = 128;
randomize(output, w, h, seed);
for(i = 0; i < 4; i++)
process(output, w, h);
Now that is not too shabby, but we can do better than that, but first we’ll divorce the randomization from the output image and fill out a separate array with values in the range of [0, 255]
.
static void randomize(
uint8_t *restrict noise,
const int w,
const int h,
const unsigned int seed
)
{
int i, n;
srand(seed);
for(i = 0, n = w * h; i < n; i++)
noise[i] = rand() & 0xFF;
}
Then we modify the processing portion to sample N pixels around the current pixel in the range of [-N/2, +N/2]
.
static void process(
const uint8_t *restrict noise,
rgba_t *restrict output,
const int w,
const int h,
const int num_samples
)
{
int i, x, y, xx, yy, yo, ww, hh, num_samples;
uint32_t c, cc;
float num_samples;
ww = w - 1;
hh = h - 1;
half_num_samples = num_samples >> 1;
inv_num_samples = 1.0f / (num_samples << 1);
for(y = 0; y < h; y++)
{
yo = y * w;
for(x = 0; x < w; x++)
{
cc = 0;
for(i = -half_num_samples; i <= half_num_samples; i++)
{
xx = clamp(x + i, 0, ww);
yy = clamp(y + i, 0, hh);
cc += noise[xx + y * w];
cc += noise[x + yy * w];
}
c = cc * inv_num_samples;
output[x + yo] = (rgba_t) {
.r = c,
.g = c,
.b = c,
.a = 0xFF
};
}
}
}
Not bad, not terrible! But, sadly it’s not the noise we are looking for.
Let’s combine the two techniques into one and see if we can get closer to the kind of noise that we are looking for. Namely, something more akin to what one gets with those “noise cloud” type filters in amazing programs like photoshop or gimp. I am sure that you know what I am talking about.
#define NOISE_MAX (1 << 14)
#define NOISE_MAX_MASK (NOISE_MAX - 1)
#define NOISE_MAX_MASK_INV (1.0f / (float) NOISE_MAX_MASK)
static void randomize(
float *restrict noise,
const int w,
const int h,
const unsigned int seed
)
{
int i, n;
srand(seed);
for(i = 0, n = w * h; i < n; i++)
noise[i] = (rand() & NOISE_MAX_MASK) * NOISE_MAX_MASK_INV;
}
#define fract(a) ((a) - ((int)(a)))
#define lerp(a, b, t) ((1 - (t)) * (a) + (t) * (b))
static void process(
const float *restrict noise,
rgba_t *restrict output,
const int w,
const int h,
const int num_samples
)
{
int x, y, i, xx0, yy0, xx1, yy1, yyo0, yyo1, yo, ww, hh;
float inv_num_samples, inv_scale, xx, yy, fx, fy;
float ac, v0, v1, c, c0, c1, c2, c3;
ww = w - 1;
hh = h - 1;
inv_num_samples = 127.0f / (float) num_samples;
for(y = 0; y < h; y++)
{
yo = y * w;
for(x = 0; x < w; x++)
{
ac = 0;
for(i = 1; i <= num_samples; i <<= 1)
{
inv_scale = 1.0f / (float) i;
xx = x * inv_scale;
yy = y * inv_scale;
xx0 = ((int) xx) & ww;
yy0 = ((int) yy) & hh;
xx1 = (xx0 + 1) & ww;
yy1 = (yy0 + 1) & hh;
yyo0 = yy0 * w;
yyo1 = yy1 * w;
c0 = noise[xx0 + yyo0];
c1 = noise[xx1 + yyo0];
c2 = noise[xx0 + yyo1];
c3 = noise[xx1 + yyo1];
fx = fract(xx);
fy = fract(yy);
v0 = lerp(c0, c1, fx);
v1 = lerp(c2, c3, fx);
ac += lerp(v0, v1, fy) * i;
}
c = ac * inv_num_samples;
output[x + yo] = (rgba_t) {
.r = c,
.g = c,
.b = c,
.a = 0xFF
};
}
}
}
First, we changed the randomize
function to generate random floating point values in the range of [0,1]
in order to increase accuracy and simplify things down the line.
Second, we now sample from [1, num_samples]
in increments of powers of two. In other words, if num_samples = 128
, then we end up sampling [1, 2, 4, 8, 16, 32, 64, 128]
.
Third, we scale the x
and y
coordinates by the current sample size. We could avoid the lonely division there, but it’s not end of the world since everything else is using multiplications and additions.
I am talking about this portion of the code of course, in case you’ve been scratching your head.
inv_scale = 1.0f / (float) i;
xx = x * inv_scale;
yy = y * inv_scale;
Forth, instead of simply taking the mean average of the surrounding pixels we perform bilinear interpolation, which is just a fancy word for linearly interpolating on the x
and y
axis at the same time.
Let’s take a deeper dive into this because this is really the most interesting part and it’s also where the so called magic happens.
xx0 = ((int) xx) & ww;
yy0 = ((int) yy) & hh;
xx1 = (xx0 + 1) & ww;
yy1 = (yy0 + 1) & hh;
yyo0 = yy0 * w;
yyo1 = yy1 * w;
c0 = noise[xx0 + yyo0];
c1 = noise[xx1 + yyo0];
c2 = noise[xx0 + yyo1];
c3 = noise[xx1 + yyo1];
fx = fract(xx);
fy = fract(yy);
v0 = lerp(c0, c1, fx);
v1 = lerp(c2, c3, fx);
ac += lerp(v0, v1, fy) * i;
I know that it looks a bit intimidating, but it’s simpler than it looks as we’ll see in a second and basically it boils down to the following operations:
- clamp
x0
andy0
- offset
x0
andy0
by+1
and clamp each - sample the noise at the following coordinates:
[x0, y0]
[x1, y0]
[x0, y1]
[x1, y1]
- interpolate
[x0, y0]
and[x1, y0]
by usingt = fract(xx)
- interpolate
[x0, y1]
and[x1, y1]
by usingt = fract(xx)
- interpolate the result of the two previous interpolations by using
t = fract(yy)
- scale the result of the interpolation by the current sample size;
lerp * i
- accumulate the result into
ac
Last by not least, we normalize the accumulated result and bring it into the [0, 127]
range.
c = ac * inv_num_samples;
If we expand this, we get the expression seen below; which of course has an unnecessary division.
c = 127 * ac / num_samples;
That’s all. Anyway, let’s see what we get if we execute this by using the following values.
variable | value |
---|---|
w | 1024 |
h | 1024 |
num_samples | 128 |
seed | 4213666 |
Voilà! Not as nice as good old perlin noise or its friends, but it’s something that one can throw together without too much fuss in a couple of minutes and it’s reasonably cheap in terms of speed.
Benchmark
$ time noise noise.png 1024 128 4213666
wrote output image to 'noise.png' ...
real 0m0.279s
user 0m0.268s
sys 0m0.010s
$ time noise noise.png 1024 128 4213666
wrote output image to 'noise.png' ...
real 0m0.271s
user 0m0.263s
sys 0m0.007s
$ time noise noise.png 1024 128 4213666
wrote output image to 'noise.png' ...
real 0m0.271s
user 0m0.266s
sys 0m0.005s
Source Code
You can find the full source code below as per usual.
/*
MIT LICENSE
Copyright (c) 2023, Mihail Szabolcs
*/
#include <stdbool.h>
#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#define STB_IMAGE_WRITE_IMPLEMENTATION
#include "stb_image_write.h"
#ifndef UNUSED
#define UNUSED(x) (void)(x)
#endif
#ifndef DEFAULT_SIZE
#define DEFAULT_SIZE 1024
#endif
#ifndef DEFAULT_NUM_SAMPLES
#define DEFAULT_NUM_SAMPLES 128
#endif
#define NOISE_MAX (1 << 14)
#define NOISE_MAX_MASK (NOISE_MAX - 1)
#define NOISE_MAX_MASK_INV (1.0f / (float) NOISE_MAX_MASK)
#define min(a, b) ((a) < (b) ? (a) : (b))
#define max(a, b) ((a) > (b) ? (a) : (b))
#define clamp(x, a, b) max(a, min(b, x))
#define fract(a) ((a) - ((int)(a)))
#define lerp(a, b, t) ((1 - (t)) * (a) + (t) * (b))
typedef union
{
uint32_t value;
struct
{
uint32_t r : 8;
uint32_t g : 8;
uint32_t b : 8;
uint32_t a : 8;
};
} rgba_t;
static void randomize(
float *restrict output,
const int w,
const int h,
const unsigned int seed
);
static void process(
const float *restrict noise,
rgba_t *restrict output,
const int w,
const int h,
const int num_samples
);
int main(int argc, char *argv[])
{
int w, h, num_samples, ret;
unsigned int seed;
float *noise;
rgba_t *output;
if(argc < 2)
{
fprintf(
stderr,
"usage: %s output.png [[size] [num_samples] [seed]]\n",
argv[0]
);
return EXIT_FAILURE;
}
if(argc > 2)
{
w = clamp(atoi(argv[2]), 16, 8192);
h = w;
}
else
{
w = DEFAULT_SIZE;
h = DEFAULT_SIZE;
}
if(argc > 3)
num_samples = clamp(atoi(argv[3]), 1, 4096);
else
num_samples = DEFAULT_NUM_SAMPLES;
if(argc > 4)
seed = atoi(argv[4]);
else
seed = time(NULL);
noise = malloc(w * h * sizeof(float));
if(noise == NULL)
{
fprintf(stderr, "failed to allocate memory for output image\n");
return EXIT_FAILURE;
}
output = malloc(w * h * sizeof(rgba_t));
if(output == NULL)
{
free(noise);
fprintf(stderr, "failed to allocate memory for output image\n");
return EXIT_FAILURE;
}
randomize(noise, w, h, seed);
process(noise, output, w, h, num_samples);
ret = stbi_write_png(
argv[1],
w,
h,
sizeof(rgba_t),
output,
w * sizeof(rgba_t)
);
if(ret == 0)
{
fprintf(stderr, "failed to write output image: '%s'\n", argv[1]);
ret = EXIT_FAILURE;
}
else
{
fprintf(stdout, "wrote output image to '%s' ...\n", argv[1]);
ret = EXIT_SUCCESS;
}
free(output);
free(noise);
return ret;
}
static void randomize(
float *restrict noise,
const int w,
const int h,
const unsigned int seed
)
{
int i, n;
srand(seed);
for(i = 0, n = w * h; i < n; i++)
noise[i] = (rand() & NOISE_MAX_MASK) * NOISE_MAX_MASK_INV;
}
static void process(
const float *restrict noise,
rgba_t *restrict output,
const int w,
const int h,
const int num_samples
)
{
int x, y, i, xx0, yy0, xx1, yy1, yyo0, yyo1, yo, ww, hh;
float inv_num_samples, inv_scale, xx, yy, fx, fy;
float ac, v0, v1, c, c0, c1, c2, c3;
ww = w - 1;
hh = h - 1;
inv_num_samples = 127.0f / (float) num_samples;
for(y = 0; y < h; y++)
{
yo = y * w;
for(x = 0; x < w; x++)
{
ac = 0;
for(i = 1; i <= num_samples; i <<= 1)
{
inv_scale = 1.0f / (float) i;
xx = x * inv_scale;
yy = y * inv_scale;
xx0 = ((int) xx) & ww;
yy0 = ((int) yy) & hh;
xx1 = (xx0 + 1) & ww;
yy1 = (yy0 + 1) & hh;
yyo0 = yy0 * w;
yyo1 = yy1 * w;
c0 = noise[xx0 + yyo0];
c1 = noise[xx1 + yyo0];
c2 = noise[xx0 + yyo1];
c3 = noise[xx1 + yyo1];
fx = fract(xx);
fy = fract(yy);
v0 = lerp(c0, c1, fx);
v1 = lerp(c2, c3, fx);
ac += lerp(v0, v1, fy) * i;
}
c = ac * inv_num_samples;
output[x + yo] = (rgba_t) {
.r = c,
.g = c,
.b = c,
.a = 0xFF
};
}
}
}
/* vim: set ts=4 sw=4 sts=4 noet: */
2023-12-17