The code that goes with this post can be found at https://github.com/Atrix256/ToUniform

My last blog post was on filtering a stream of uniform white noise to make a stream of red or blue noise, and then using the CDF (inverted inverse CDF) to transform the non uniform red or blue noise back into a uniform distribution while preserving the noise color. It is at https://blog.demofox.org/2023/02/20/uniform-1d-red-noise-blue-noise/

After that post, Bart commented:
https://mastodon.gamedev.place/@demofox/109909020853454095

That had me thinking about how you could make a CDF of any distribution, numerically. You could then store it in a LUT, or perhaps use least squares to fit a polynomial to the CDF.

I also got this interesting comment from Graeme.
https://mastodon.gamedev.place/@graeme/109916215756255970

I had never seen that before, only ever seeing that you can subtract a low pass filter from the identity filter to get a high pass filter (https://blog.demofox.org/2022/02/26/image-sharpening-convolution-kernels/). It made me want to check out that kind of a filter to see what the differences could be.

So, I started by generating 10 million values of white noise and then filtering it. I made a histogram of the results and divided the counts by 10 million to normalize it into a PDF. From there I made a CDF by making a table where each CDF bucket value was the sum of all PDF values with the same or lower index. From there I made a 1024 entry CDF table, as well as a smaller 64 entry one, and I also fit the CDF with piecewise polynomials (https://blog.demofox.org/2022/06/29/piecewise-least-squares-curve-fitting/).

With any of those 3 CDF approximations, I could filter white noise in the same way I did when generating the CDF, put the colored noise through the CDF function, and get out colored noise that was uniform distribution.

I wasn’t very happy with the results though. They were decent, but I felt like I was leaving some quality on the table. Going from real values to indices and back was a bit fiddly, even when getting the nuances correct. I mentioned that it was surprisingly hard to make a usable CDF from a histogram, and Nikita gave this great comment.

It took me a minute to understand but this is quite an amazing insight: If you generate a lot of numbers from an unknown PDF and sort them you now have a LUT of the inverse CDF! If you are wondering why that is such an exciting thing, give this a read https://blog.demofox.org/2017/08/05/generating-random-numbers-from-a-specific-distribution-by-inverting-the-cdf/.

The intuition here is that regions of values that are more probable will take up more indices and so if you randomly select an index, you are more likely to select from those regions.

So, if you have a black box that is giving you random numbers, you can generate random numbers from the same distribution just by storing a bunch of the numbers it gives you, sorting them, and using uniform random numbers to select items from that list. Better yet, select a random position between 0 and 1, multiply by the size of the list, and use linear interpolation between the indices when you get a fractional index.

Andrew Helmer let me know that this is called the “Empirical Distribution Function” (thanks!) https://en.wikipedia.org/wiki/Empirical_distribution_function

In our case where we are trying to make non uniform numbers be uniform again, we actually want the inverted inverse CDF, or in other words, we want the CDF. We can get a CDF value from the inverse CDF table by taking a value between 0 and 1 and finding where it occurs in the inverse CDF table. We’ll get a fractional index out if we use linear interpolation, and we can divide by the size of the array to get a value between 0 and 1. This is the CDF value for that input value. If we do this lookup for N evenly spaced values between 0 and 1, we can get a CDF look up table out.

I did that with a 1024 sized look up table, a 64 sized look up table, and I also fit a piecewise polynomial to the CDF, with piecewise least squares.

The results were improved! Fewer moving parts and fewer conversions. Nice.

I’ll show the results in a minute but here are the tests I did:

• Box 3 Red Noise: convolved the white noise with the kernel [1, 1, 1]
• Box 3 Blue Noise: convolved the white noise with the kernel [-1, 1, -1]. This is an unsharp mask because I subtracted the red noise kernel from [0,2,0].
• Box 5 Red Noise: [1,1,1,1,1]
• Box 5 Blue Noise 1: [-1, -1, 1, -1, -1]
• Box 5 Blue Noise 2: [1, -1, 1, -1, 1]
• Gauss 1.0 Blue Noise: A sigma 1 gaussian kernel, with alternating signs to make it a high pass filter [0.0002, -0.0060, 0.0606, -0.2417, 0.3829, -0.2417, 0.0606, -0.0060, 0.0002]. Kernel calculated from here http://demofox.org/gauss.html
• FIRHPF – a finite impulse response high pass filter. It’s just convolved with the kernel [0.5, -1.0, 0.5]. Kernel calculated from here http://demofox.org/DSPFIR/FIR.html/
• IIRHPF – an infinite impulse response high pass filter. It has the same “x coefficients” as the previous filter, but has a y coefficient of 0.9. Filter calculated from here http://demofox.org/DSPIIR/IIR.html.
• FIRLPF – the low pass filter version of FIRHPF. kernel [0.25, 0.5, 0.25]
• IIRLPF – the low pass filter version of IIRHPF. x coefficients are the same as the previous filter, but has a y coefficient of -0.9.
• Void And Cluster – I used the void and cluster algorithm to make 100,000 blue noise values, and I did this 100 times to get 10 million total blue noise values. There is a “seam” in the blue noise every 100,000 values but that shouldn’t be a big deal.
• Final BN LUT – I took the sized 64 CDF LUT from the FIRHPF test and used it to turn filtered white noise back to uniform.
• Final BN Polynomial – I took the 4 piecewise quadratic functions that fit the CDF of FIRHPF and used it to turn filtered white noise back to uniform.
• Final RN Polynomial – The filter from FIRLPF is used, along with the polynomial from the previous filter to make it uniform again. Luckily FIRLPF and IIRLPF have the same histogram and CDF, so can use the same polynomials!

CDFs

In the below, the blue line is the 4 piecewise quadratic functions fit to the CDF. The green dotted line is the size 1024 CDF table. The orange line is the size 64 CDF table.

For void and cluster, and the “Final BN” tests, the source values are uniform, so it doesn’t make much sense to make a CDF and convert them to uniform etc, but it wasn’t trivial to disable that testing for them.

Here is the error of the CDF approximations, using the size 1024 CDF LUT as the ground truth to calculate the error.

Frequencies (DFTs)

Here are the DFTs to show the frequencies in the streams of noise. Converting filtered white noise to uniform does change the frequency a bit, but looking at the low frequency magnitudes of void and cluster, that might just be part of the nature of uniform distributed colored noise. I’m not sure.

The blue line is the original values from E.g. filtering white noise. The other colors are when converting them to uniform.

Histograms

Here are the histograms, to see if they are appropriately uniform. The histograms come in sets of 4: the filtered noise, and then 3 methods for converting the filtered noise back to uniform. The graphs should be red from top to bottom, then left to right.

Streaming Blue Noise

Looking at the results, the FIRHPF results look pretty decent, with the polynomial fit having a better histogram than the size 64 CDF LUT. Only being a 3 tap filter means that it only needs to remember the last 2 “x values” that it has seen too, as a “delay line” for the FIR filter. Decent quality results and efficient to obtain. nice. The results of that are in the diagrams above as “Final BN Polynomial”.

I tried taking only the first half of the fits – either 32 of the 64 LUT entries, or only the first two of the four polynomials. For the second half of the curve, when x > 0.5, you set x to1-x, and the resulting y value gets flipped as 1-y too. I wasn’t getting great results though, even when adding a least squares constraint that CDF(0.5) = 0.5. I’m saving that battle for a future day. If you solve it though, leave a comment!

Below is the C++ code for streaming scalar valued blue noise – aka a blue noise PRNG. Also the same for red noise. The full source code can be found at https://github.com/Atrix256/ToUniform.

If you have a need of other colors of noise, of other histograms, hopefully, this post can help you achieve that.

class BlueNoiseStreamPolynomial
{
public:
	BlueNoiseStreamPolynomial(pcg32_random_t rng)
		: m_rng(rng)
	{
		m_lastValues[0] = RandomFloat01();
		m_lastValues[1] = RandomFloat01();
	}
	
	float Next()
	{
		// Filter uniform white noise to remove low frequencies and make it blue.
		// A side effect is the noise becomes non uniform.
		static const float xCoefficients[3] = {0.5f, -1.0f, 0.5f};

		float value = RandomFloat01();

		float y =
			value * xCoefficients[0] +
			m_lastValues[0] * xCoefficients[1] +
			m_lastValues[1] * xCoefficients[2];

		m_lastValues[1] = m_lastValues[0];
		m_lastValues[0] = value;

		// the noise is also [-1,1] now, normalize to [0,1]
		float x = y * 0.5f + 0.5f;

		// Make the noise uniform again by putting it through a piecewise cubic polynomial approximation of the CDF
		// Switched to Horner's method polynomials, and a polynomial array to avoid branching, per Marc Reynolds. Thanks!
		float polynomialCoefficients[16] = {
			5.25964f, 0.039474f, 0.000708779f, 0.0f,
			-5.20987f, 7.82905f, -1.93105f, 0.159677f,
			-5.22644f, 7.8272f, -1.91677f, 0.15507f,
			5.23882f, -15.761f, 15.8054f, -4.28323f
		};
		int first = std::min(int(x * 4.0f), 3) * 4;
		return polynomialCoefficients[first + 3] + x * (polynomialCoefficients[first + 2] + x * (polynomialCoefficients[first + 1] + x * polynomialCoefficients[first + 0]));
	}

private:
	float RandomFloat01()
	{
		// return a uniform white noise random float between 0 and 1.
		// Can use whatever RNG you want, such as std::mt19937.
		return ldexpf((float)pcg32_random_r(&m_rng), -32);
	}

	pcg32_random_t m_rng;
	float m_lastValues[2] = {};
};

class RedNoiseStreamPolynomial
{
public:
	RedNoiseStreamPolynomial(pcg32_random_t rng)
		: m_rng(rng)
	{
		m_lastValues[0] = RandomFloat01();
		m_lastValues[1] = RandomFloat01();
	}

	float Next()
	{
		// Filter uniform white noise to remove high frequencies and make it red.
		// A side effect is the noise becomes non uniform.
		static const float xCoefficients[3] = { 0.25f, 0.5f, 0.25f };

		float value = RandomFloat01();

		float y =
			value * xCoefficients[0] +
			m_lastValues[0] * xCoefficients[1] +
			m_lastValues[1] * xCoefficients[2];

		m_lastValues[1] = m_lastValues[0];
		m_lastValues[0] = value;

		float x = y;

		// Make the noise uniform again by putting it through a piecewise cubic polynomial approximation of the CDF
		// Switched to Horner's method polynomials, and a polynomial array to avoid branching, per Marc Reynolds. Thanks!
		float polynomialCoefficients[16] = {
			5.25964f, 0.039474f, 0.000708779f, 0.0f,
			-5.20987f, 7.82905f, -1.93105f, 0.159677f,
			-5.22644f, 7.8272f, -1.91677f, 0.15507f,
			5.23882f, -15.761f, 15.8054f, -4.28323f
		};
		int first = std::min(int(x * 4.0f), 3) * 4;
		return polynomialCoefficients[first + 3] + x * (polynomialCoefficients[first + 2] + x * (polynomialCoefficients[first + 1] + x * polynomialCoefficients[first + 0]));
	}

private:
	float RandomFloat01()
	{
		// return a uniform white noise random float between 0 and 1.
		// Can use whatever RNG you want, such as std::mt19937.
		return ldexpf((float)pcg32_random_r(&m_rng), -32);
	}

	pcg32_random_t m_rng;
	float m_lastValues[2] = {};
};

Update 3/12/23

Nick Appleton has another way to make streaming blue noise (https://mastodon.gamedev.place/@nickappleton/110012203119490715).

I coded it up and put it through my analysis so you could compare apples to apples. It looks pretty good!

Histograms

DFTs

Code

// From Nick Appleton:
// https://mastodon.gamedev.place/@nickappleton/110009300197779505
// But I'm using this for the single bit random value needed per number:
// https://blog.demofox.org/2013/07/07/a-super-tiny-random-number-generator/
// Which comes from:
// http://www.woodmann.com/forum/showthread.php?3100-super-tiny-PRNG
class BlueNoiseStreamAppleton
{
public:
	BlueNoiseStreamAppleton(unsigned int seed)
		: m_seed(seed)
		, m_p(0.0f)
	{
	}

	float Next()
	{
		float ret = (GenerateRandomBit() ? 1.0f : -1.0f) / 2.0f - m_p;
		m_p = ret / 2.0f;
		return ret;
	}

private:
	bool GenerateRandomBit()
	{
		m_seed += (m_seed * m_seed) | 5;
		return (m_seed & 0x80000000) != 0;
	}

	unsigned int m_seed;
	float m_p;
};