When looking at the politics of the American right and the American left, there is no contest for me as to which to support, if these are my only options.
I see the right doing hateful, authoritarian things that seem to violate their own principles. They violate their religious principles as a largely Christian group who say they want to emulate Jesus but then seem to do the exact opposite at every turn. WWJD seems to be about counterexamples. Non-religious principles are also violated, such as personal liberty, small government, and upholding the constitution (all of it), while trying to force their will on everyone else. Republicans are currently pushing us into a dictatorship, thwarting the balance of powers that were set up by colluding across government branches to dismantle the country, and really seeming to be going against everything our constitution is about, cackling in our face the whole time, right out in the open. Most recently this was about the white house openly declaring they wanted to be able to “deport” (sic) citizens to prisons in other countries (source: https://bsky.app/profile/onestpress.bsky.social/post/3lmd2n77ko22g). I am for personal freedom, liberty, small government, and the values of the Constitution. If Republicans stood for what they claim to, I might very well vote Republican, and do so proudly. These people just have a way of twisting these values to be about serving them only, at the cost of others losing these same things, and more.
I see the left (Democrats) trying to do a little bit of good, although they are largely pro corporation (vs the individual), and also support / do shady things militarily to other countries. It isn’t lost on me that the DMCA was signed into law by Bill Clinton which calls for search and seizure without a warrant, or that the over-response to 9/11 hasn’t been curbed at any of the possible points with democrats in power.
There is a lesser of two evils here that I see, and that is Democrats – and it isn’t even close. They aren’t perfect, but they aren’t actively, overtly, trying to do evil and destroy our country. I don’t see any realistic alternatives either.
However, there is a belief that the Democrats did everything they could do to stop Trump and all of this from happening and that other things are to blame – like perhaps not enough voters turned out, and so it’s the fault of the voters.
No matter what the reason is for the loss, for us to consider the Democrats blameless in the current state of affairs is to accept a useless medical diagnosis for the country. It essentially says that we did everything we could, we still lost, and there is no way to win. Where does that leave us other than giving up? What’s the next step besides smugly watching the world burn and saying “I told you so”? Hopelessness is a terminal diagnosis without a treatment plan.
Consider the CEO of a failing company blaming the customers for not buying the product, or the employees for not performing well enough. Even if these things are true, it’s within the CEO’s responsibility to make the company succeed as well as possible, under any circumstances. Just because a failure can be explained doesn’t excuse the responsibility of failure.
Similarly, the democratic leaders we trusted to help preserve our country have failed. They may have failed despite their best efforts, but they still failed.
Consider this: If the Republican misinformation machine of Fox News et al. can mobilize and motivate people into action, why can’t the Democrats?
I’d like to see a stronger party more about individual empowerment, tougher on corporations and their monetary influence in politics, breaking up actual monopolies, and with the energy and cunning to fight literal evil and a stronger conscience to actually do what is right.
Winning is possible, but the Democratic party opposing the Republicans is not up to the task. So the question becomes: What do we do now? I’m not sure, but I’m not giving up, and I hope you aren’t giving up, either. Let’s keep our eyes open and do what we can, even if only small acts of kindness to foster the world we want, until larger ideas come into view that we can all get behind.
If you have other ideas, speak up. We need a bit of a country wide brainstorming session right now to figure out how to attempt to fix this mess.
FYI this blog is proudly woke: Trans rights are human rights, Israel needs to stop the genocide in Palestine, Trump is everything the founding fathers of the USA were trying to protect the country against, and black lives matter.
Still here? Cool. Hello and welcome!
Imagine you sit down to board game night, pop open the box, and the dice are missing. Ack, what do you do? Don’t worry, I got you fam. We can simulate dice with coins, or coin like objects.
I want you to consider a problem: How can you simulate rolling a 6-sided dice using 3 coins? Keep that in your head as you read, and see if you can think of a way to do it before I give away the answer lower down on the page.
If you have one coin, you can flip a 2 sided dice. Tails could be 0, and Heads could be 1. (Add 1 to the numbers if you want to get a 1 or a 2).
If you have two coins, you can flip a 4 sided dice. You multiply one of the coins by 2, and add the other to it. If you know binary, these coins are just binary digits.
Coins
Binary Digits
Numbers
Result
Dice #
Tails, Tails
00
0 + 0
0
1
Tails, Heads
01
0 + 1
1
2
Heads, Tails
10
2 + 0
2
3
Heads, Heads
11
2 + 1
3
4
If you have three coins, you can flip an 8 sided dice in the same way.
Coins
Binary Digits
Numbers
Result
Dice #
Tails, Tails, Tails
000
0 + 0 + 0
0
1
Tails, Tails, Heads
001
0 + 0 + 1
1
2
Tails, Heads, Tails
010
0 + 2 + 0
2
3
Tails, Heads, Heads
011
0 + 2 + 1
3
4
Heads, Tails, Tails
100
4 + 0 + 0
4
5
Heads, Tails, Heads
101
4 + 0 + 1
5
6
Heads, Heads, Tails
110
4 + 2 + 0
6
7
Heads, Heads, Heads
111
4 + 2 + 1
7
8
This all works just fine, and each number that might come up has equal probability of coming up (it is a uniform distribution) so it gives fair rolls. But, it only works if you need dice that have a size that are power of 2. What do you do if you want to roll a 6 sided dice?
So, going back to the puzzle at the beginning: How do you simulate rolling a 6 sided dice using 3 coins?
You might think “each coin has 2 options, and there are three of them, so just add them up!”
If tails are 0 and heads are 1, the lowest number you can get is 0, and the highest number you can get is 3. You could add 1, to make it so this was a dice that gave you random numbers between 1 and 4, but it turns out to not have a uniform distribution!
Coins
Numbers
Result
Dice#
TTT
0 + 0 + 0
0
1
TTH
0 + 0 + 1
1
2
THT
0 + 1 + 0
1
2
THH
0 + 1 + 1
2
3
HTT
1 + 0 + 0
1
2
HTH
1 + 0 + 1
2
3
HHT
1 + 1 + 0
2
3
HHH
1 + 1 + 1
3
4
As you can see, there is only 1 way to get a one, and only 1 way to get a four, but there are 3 ways to get a two, and 3 ways to get a three. That means it’s 3 times as likely to get a 2 or 3, than it is to get a 1 or a 4.
So, adding up 3 coins doesn’t give us a 6 sided dice. It gives us a 4 sided dice, that doesn’t even give fair rolls. Lame!
There is a simple answer though: rejection sampling.
You flip 3 coins and that gives you a number between 1 and 8. If the number is 7 or 8, you ignore it and flip the 3 coins again.
Doing this, the numbers you KEEP (1 through 6) come up with equal probability (1/6th). The only downside is you don’t know how many times you might have to flip the coins to get a valid dice roll. That isn’t the worst thing in the world though. Sometimes when rolling dice, it falls off the table and goes under the couch or something and you have to reroll anyways. You can think of it in the same kind of way.
Advanced Tactic: When you throw out the 7 or 8, you aren’t completely empty handed. That 7 or 8 is actually 1 bit of random information you can keep. It’s essentially a uniform random coin flip, and it has no correlation with the 6 sided dice rolls, or any previous or future bit of info generated in this way. You could stash this bit away somewhere if you wanted – like write down the stream of 7s and 8s on a piece of paper, aka 0s and 1s, and then have a random bit stream for cryptographic purposes later, as a treat. Or you can use that as one of the coin flips for the reroll – flip 2 coins instead of three, and use this bit as the third coin.
So that works for a 6 sided dice, but how about a 20 sided dice?
Well, flipping 5 coins is the same as rolling a 32 sided dice and you can still use rejection sampling. If the number is between 1 and 20, keep it. If it’s bigger than 20 ignore it.
Advanced Tactic: Similarly to flipping 3 coins to either get a 6 sided dice roll, or a coin flip, the information we throw away simulating a 20 sided dice doesn’t have to be a complete waste. Interestingly, if you flip 5 coins to simulate a 20 sided dice, you throw out values between 21 and 32. If you subtract 20 from the number you throw away, you get values between 1 and 12. That’s right! Flipping 5 coins to simulate a d20 will give you a d12 when it fails that you can write down and save for later.
Flipping coins to simulate dice can also help you roll dice that would be very difficult to make in real life, such as 3 sided dice (which you can do with 2 coin flips) or 100 sided dice (which you can do with 7 coin flips).
It’s possible that each time you flip the coins, you only get heads for a very long time, and are there flipping coins for literally years, but it’s very improbable. In the context of computer programming, this “unknown and unbounded” run time is not great, but sometimes, it’s the best strategy there is. Like for instance if you can pick a random uniform point in a square, but you need a random uniform point in a star shape that was drawn by an artist, you can use rejection sampling to make it happen. (for non uniform points in arbitrary shapes, check out blue noise point sets!)
You can actually calculate how much (on average) the rejection sampling will fail though. It’s just the percentage of the values that you throw away.
6 Sided Dice: You use three coins to make 8 values, and throw 2 of them away (7 and 8), which means you should only have to reroll 2/8 or 1/4 of the time.
20 Sided Dice: You use 5 coins to make 32 values, and throw away 12 of them. So you reroll 12/32 or 3/8 of the time.
3 Sided Dice: You use 2 coins to make 4 values and throw away 1 of them. So you reroll 1/4 of the time.
100 Sided Dice: You roll 7 coins to make 128 values and you throw away 28 of them. So you reroll 28/128 or 7/32 of the time.
When you keep the randomness generated from your rejections, you don’t lose anything really, but these formulas should help you understand how much random information is being generated in each stream over time, which in a computerized system could cause latency or thread starvation type scenarios.
One final thing – what do you do if you only have unfair coins (or coin like objects) that don’t have equal odds for coming up heads vs tails? Well, we take a hint from John Von Neumann and flip two coins (or one twice). If they mismatch, you use the first coin as the coin flip, else you ignore it and flip twice again. https://en.wikipedia.org/wiki/Fair_coin#Fair_results_from_a_biased_coin
I stumbled across a very cool YouTube video today that talks about how to look for alternate ways of labeling sides of two dice to give you the same random number distribution as if you added two standard dice together. https://www.youtube.com/watch?v=xVRRykIyll0
Making customized dice give desired probability distributions is an interesting topic (like for game design!), but the math that was used in the video was also pretty neat. I’ll walk you through what they did and then we can talk about why it works, from two different points of view.
(I’m going to refer to a single die as “a dice”, apologies in advance!)
The Problem Statement
If you roll a single dice, the probability of getting each possible value is the same: 16.666…%.
However, if you roll two dice and add them together, the probability of getting each possible value is not equal. You can only get a two by rolling a 1 and a 1, but you can roll a four by rolling a 1 and a 3, a 2 and a 2, or a 3 and a 1. There is only one way to roll a two, but there are three ways to roll a four. So, rolling a four is three times more likely than rolling a two. The probability for each possible value is described by the binomial probability distribution below.
Dice Sum
# Of Ways
Probability (Rounded)
Dice Values
2
1
3%
11
3
2
6%
12, 21
4
3
8%
13, 22, 31
5
4
11%
14, 23, 32, 41
6
5
14%
15, 24, 33, 42, 51
7
6
17%
16, 25, 34, 43, 52, 61
8
5
14%
26, 35, 44, 53, 62
9
4
11%
36, 45, 54, 63
10
3
8%
46, 55, 64
11
2
6%
56, 65
12
1
3%
66
The problem statement is this: Is there a different way to label two dice such that when you add them together, you get the exact same probability distribution as when using standard dice?
It turns out there is exactly one other way to label the dice, and they are called “Sicherman Dice”, named after the person who invented them, in 1978.
Dice A: 1, 2, 2, 3, 3, 4
Dice B: 1, 3, 4, 5, 6, 8
How To Find These Dice
The video explains that you can turn a dice into a polynomial summing up terms of the form . The value is what number is on the dice face, and the value says how many faces have that value.
For instance, a standard dice which has faces 1,2,3,4,5,6 would have this polynomial below, which states that it has one of each number 1 through 6:
If we changed the 2 face to be a 4, we’d have the polynomial below, which says we have faces 1,3,4,4,5,6:
These polynomials are valuable because if you have two of them, you can multiply them together to get a “virtual dice” that is the same as if we added the dice together into a single dice.
If we multiply two standard dice together:
We get the polynomial:
If you add up how many faces there are, that would be a 36 sided dice, and if you look at the number of faces of each value, you can see they follow the same binomial distribution that we showed in the problem statement section, when adding two dice together.
We aren’t done yet though, we need to see if there’s an alternate way to break this 36 sided dice polynomial back into two 6 sided dice, that aren’t just the standard 6 sided dice.
We need to factor that large polynomial into two smaller ones. These smaller polynomials have requirements (constraints) as well, to be valid:
Each needs to be non negative. That is, you can’t have a negative number of face values on a dice. That’d be weird.
The labels need to be positive. 0 is not a valid face value. This means that so that is 0.
The dice should have 6 faces. So if you sum up for all , it should give you 6.
The second bullet point can be described as because it makes all the power terms go away and leaves only the value of .
The third bullet point can be described as because that will make all power terms become 1, and so it will just be the values summed up.
The first step is to decompose the large polynomial into irreducible (prime) factors. In the video, instead of factoring the large polynomial, he factors the standard dice polynomial instead, and duplicates the factors, to get the list below. This video is a good intro to factoring polynomials if you need it: https://www.youtube.com/watch?v=KUMhpKGwpCY
To make new dice, we are going to assign each of those factors above to one of the dice. The factors assigned to a dice are going to be multiplied together to get the final polynomial for the dice.
We start with the second constraint that . The only way that can be true is if each dice has the factor. So, we assign one of those to the first dice, and the other one to the second dice.
Dice A:
Dice B:
Let’s move onto the second constraint that . If we plug 1 in for , the factor we already have for each dice evaluates to 1,and here are the values we get for each of the factors remaining:
Since we are multiplying factors together, and we need both of the dice to evaluate to 6, that means we assign to each dice to multiply by 2, and also assign to each dice to multiply by 3, so that we are multiplying each dice by 6.
Our dice are identical right now, looking like this:
Dice A:
Dice B:
We also have two factors left to assign:
If we gave each dice one of these factors, we’d end up at the standard dice again, so the idea is to give both factors to one of the dice to get this:
Dice A:
Dice B:
We can multiply these out to get the following polynomials:
Dice A:
Dice B:
Listing out the faces of the dice, we find we’ve made the Sicherman dice! Also, all three constraints are satisfied by both dice.
Dice A: 1, 2, 2, 3, 3, 4
Dice B: 1, 3, 4, 5, 6, 8
Explanation As Convolution
When I saw this, it made sense to me in terms of convolution. There are two concepts here to show what I mean.
The first concept is that if you add together two random values, you are also convolving their probability distributions. For discrete random values, such as dice, you are convolving their Probability Mass Function (PMF) which is a normalized histogram. For continuous random values, you are convolving their Probability Density Function (PDF). See this for more information: https://en.wikipedia.org/wiki/Convolution_of_probability_distributions
The second concept is that when you multiply polynomials together, you are actually doing convolution of their coefficients.
For instance, if you wanted to convolve the sequence 1,1 with the sequence 1,1 that would give you the sequence 1,2,1.
We can set this up with polynomials by giving each term in the sequence a unique power of x:
You can see that the coefficients of the result are 1,2,1 and gave us the same answer as convolution.
We got our answer in the coefficients again: 2, 6, 10, 14, 8.
So, if we want to get the probability distribution of the sum of two dice, we know we need to convolve their distribution. By making the dice histogram into polynomials, we can multiply the polynomials to do the convolution, and get our resulting histogram. The polynomial then just needs to be factored into two polynomials that satisfy the dice validity constraints.
Viewed from this perspective, this method seems like a “hacky party trick”.
Interestingly though, the part about factoring a polynomial into smaller polynomials can be seen as deconvolution in this point of view. Deconvolution is a powerful thing, and kind of an advanced topic, but to bring it to something tangible, if you have a blurry image, and you know the kernel that blurred the image, deconvolution lets you get the unblurred source image back.
As we saw in this situation though, factoring a polynomial into two other polynomials multiplied together does not always give a unique solution, and it isn’t an easy process to find and enumerate the solutions. This is also true in deconvolution, unsurprisingly! There are multiple methods for doing deconvolution. This polynomial factoring is equivalent to just one of the methods.
Generating Functions Introduction
The other way I wanted to talk about this problem is from the point of view of generating functions, which generalize this creative use of exponents into a powerful mathematical tool. Generating functions allow you to abstract certain types of problems as polynomial multiplication.
A generating function is a polynomial just like the ones we’ve been looking at: a sum of powers of x multiplied by coefficients.
When looking at a specific term , the value says how many ways there are to select the value .
If we have a single letter A that we can select or not, there is 1 way to select 0 elements, and 1 way to select 1 elements. That gives a generating function which simplifies to:
.
If we have the letters AA, there is 1 way to select 0 elements, 1 way to select 1 element (because the As are identical to each other), and there is 1 way to select 2 elements. That gives the generating function:
If we also have the letters BB, they have the same generating function:
To get the generating function for AABB, we are selecting from AA and from BB and adding them together, so we multiply the generating functions together:
Which equals:
That shows us there is:
1 way to select 0 elements
2 ways to select 1 elements (just A or just B)
3 ways to select 2 elements (AA, BB, AB)
2 ways to select 3 elements (AAB or ABB)
1 way to select 4 elements (AABB)
Example 2 – ABCD
Going back to a single letter choice, we can take the letter or leave it, giving us the generating function:
.
If we have four letters we can choose from (or not) – ABCD – we multiply that generating function by itself four times: once per letter.
That equals:
This shows us there is:
1 way to select 0 elements
4 ways to select 1 element (A, B, C or D)
6 ways to select 2 elements (AB, AC, AD, BC, BD, CD)
4 ways to select 3 elements (BCD, ACD, ABD, ABC)
1 way to select 4 elements (ABCD)
Fun fact: you may notice this is the 5th row of Pascal’s triangle. This relates to binomial expansion and the binomial theorem. This also relates to Bernouli trials, since we are either taking each letter or leaving it. This series of values is also is what happens if you convolve 1,1 against itself, four times!
Problem 1 – Summing Three Numbers
Let’s work through two simple problems to better understand the application of generating functions before talking about the dice problem. These examples come directly from this video https://www.youtube.com/watch?v=-drdeNMoe8w.
For the first problem, we want to know how many ways are there to sum three numbers to get an answer of 12, when each number is between 0 and 6 inclusively.
We want to know how many ways can be true, when .
Choosing a number between 0 and 6 is represented by the generating function . If we want to do that three times – once for A, once for B, once for C – and add them together, we multiply that generating function by itself three times, aka we cube it.
That gives:
The coefficient for is 28, which means there are 28 different ways to add the three numbers between 0 and 6 together, to get 12.
Problem 2 – Picking Jelly Beans
Imagine we have a jar of red, blue and white jellybeans. We want to know how many ways there are to select 20 jelly beans such that:
We have an even number of red jelly beans
We have at least 14 blue jelly beans.
We have less than 5 white jelly beans.
Those three constraints become these three generating functions:
The first and second generating functions could extend beyond into infinity, but we know there are 20 jellybeans total, so we can stop at 20.
Multiplying those three together gives us the polynomial below, which shows us that has a coefficient of 14. There are 14 different ways to select 20 jellybeans while satisfying the constraints.
Explanation As Generating Functions
Looking at the original dice summing problem from the perspective of generating functions, you can see that we have individual dice that have generating functions describing how many ways there are to roll each number.
With a standard 6 sided dice, the generating function only has one way to roll each number:
For our modified dice where we replaced the 2 face with a 4, that meant there was no way to get a 2 anymore, but there were two ways to get a 4:
Multiplying dice polynomials together gave the histogram of how many ways each value could come up if you rolled two dice and added them together.
In the example problems of generating functions in the last section, we only concerned ourselves with a single coefficient, because we wanted to know how many ways a specific value could be selected.
In the dice summing problem though, we are concerned with ALL coefficients. The reason for this is because we are trying to get a specific histogram of possible selected values. A histogram is an array of counts, which is what the coefficients are. So, our dice problem is just a generating function problem which concerns itself with multiple coefficients. That’s still well within the realm of normality for generating functions.
The generating function point of view of this problem helps give a more formal footing to this method, rather than the “hacky party trick” point of view, but I don’t believe it helps you factor the final generating function of the summed dice into the two individual generating functions needed for the alternately labeled dice. It’s just another way of explaining why the polynomial shenanigans work.
Closing
If you’ve done any DSP work, as you learn more about generating functions you may notice they are pretty similar to the z transform (the discrete laplace transform). In generating functions, you can shift the index of the coefficient you need to calculate by multiplying by positive or negative powers of x, just like you can multiply by z to shift a sample in time. I have two write ups on the DSP side of things regarding this if you want to know more:
There is a method to get a single coefficient without having to calculate them all, using the Extended Binomial Theorem. That is explained in the second link above.
In this post I’ll be showing two different ways to scale points along a specific direction (a vector). There isn’t anything novel here, but it gives an example of the type of math encountered in game development, and how you might approach a solution. The first method involves matrices, and the second involves vector projections.
The points in the red circle below are scaled along the dark green vector, with the magnitude of the bright green vector.
This problem came up for me when implementing a Petzval lens effect (https://bbhphoto.org/2017/04/24/the-petzval-lens/) in a bokeh depth of field (https://blog.demofox.org/2018/07/04/pathtraced-depth-of-field-bokeh/) post processing rendering technique. The setup is that I take point samples in the shape of the camera aperture for bokeh and depth of field, and I needed to stretch this sampled shape in a specific direction and magnitude, based on where the pixel is on the screen.
In the end, I just needed to scale points in a specific direction, by a specific amount. If we can do this operation to 1 point, we can do it to N points, so we’ll focus on doing this to a single point.
Method 1 – Matrices
Using matrices to scale a point along a specific direction involves three steps:
Rotate the point so that the direction we want to scale is aligned with the X axis.
Multiply the x axis value by the scaling amount.
Unrotate the point back to the original orientation.
Let’s say we want to scale a point by 4 along the vector [3,1].
First we normalize that vector to [3/sqrt(10), 1/sqrt(10)]. Then we need to get the rotation matrix that will transform [3/sqrt(10), 1/sqrt(10)] to [1,0]. We can do that by making a matrix where the first column is where we want the x axis to point, and the second column to be where we want the y axis to point. For the x axis, we use the normalized vector we already have. For the y axis, we just need a perpendicular vector, which we can get by swapping the x and y components of the vector, and negating one to get [-1/sqrt(10), 3/sqrt(10)]. This operation of swapping x and y and negating one is sometimes called the “2D Cross Product” even though it isn’t really a cross product, and there is no cross product in 2D. That gives us this rotation matrix:
Next we need to make the scaling matrix, which we get by multiplying the vector [4,1] by the identity matrix:
Lastly we need to calculate the unrotation matrix. We need a matrix that rotates by the negative amount of R. We need the inverse matrix of R. The inverse of a rotation matrix is just the transpose matrix, so we can transpose R to make it. Another way to think about it is when we made the matrix before with the first column being the x axis and the second column being the y axis, we are now going to make the first row be the x axis, the second row be the y axis. Rows instead of columns. Whichever explanation makes most sense to you, we end up with this:
Now that we have all the transformations, we can calculate R * S * R’ to get a final matrix that does the transformation we want. I’ll do it in 2 steps in case that helps you follow along, to make sure you get the same numbers.
That is our matrix which scales a point along the vector [3,1], with a magnitude of 4. Let’s put the vector [4,5] through this transformation by multiplying it by the matrix.
For those who are counting instructions, processing a point using this process is 2 multiplies and two adds, to do that 2d vector / matrix product.
Creating the matrix took 16 multiplies and 8 adds (two matrix / matrix multiplies aka eight 2d dot products), but usually, you calculate a matrix like this once and re-use it for many points, which makes the matrix creation basically zero as a percentage of the total amount of calculations done, when amortized across all the points.
Method 2 – Vector Projection
I’m a fan of vector projection techniques. There is a certain intuitiveness in them that is missing from matrix operations, I find.
Using vector projection to scale a point along a specific direction involves these three steps:
Project the point onto the scaling vector and multiply that by the scaling amount.
Project the point onto the perpendicular vector.
Add the scaling vector projection times the scaling vector to the perpendicular vector projection times the perpendicular vector.
We will use the same values from the last section, so we want to scale a point by 4 along the vector [3,1], which we normalize to [3/sqrt(10), 1/sqrt(10)]. We will put the point [4, 5] through this process.
Step 1 is to project our point onto the scaling vector. We do that by doing a dot product between our normalized vector [3/sqrt(10), 1/sqrt(10)], and our point [4, 5]. That gives us the value 17/sqrt(10). We then multiply that by the scaling amount 4 to get 68/sqrt(10).
Step 2 is to project our point onto the perpendicular vector. We can once again use the “2D Cross Product” to get the perpendicular vector. We just flip the x and y component and negate one, to get the vector perpendicular to the scaling vector: [-1/sqrt(10), 3/sqrt(10)]. We can dot product that with our point [4, 5] to get: 11/sqrt(10).
Step 3 is to multiply our projections by the vectors we projected onto, and add the results together. Our scaling vector contribution is 68/sqrt(10) * [3/sqrt(10), 1/sqrt(10)] or [204/10, 68/10]. Our perpendicular vector contribution is 11/sqrt(10)*[-1/sqrt(10), 3/sqrt(10)] or [-11/10, 33/10].
When we add the two values together, we get [193/10, 101/10] or [19.3, 10.1].
That result matches what we got with the matrix operations!
As far as instruction counts, we did two dot products for step 1 and 2 which is 4 multiplies and 2 adds total. Step 3 is 4 multiplies. Step 4 is 2 multiplies and 2 adds. This is a total of 10 multiplies and 4 adds which is a lot more than the matrix version, which was just 2 multiplies and 2 adds.
If you optimized this process to do fewer operations by combining work where you could, you’d eventually end up at the same operations done in the matrix math. Algebra is fun that way.
Higher Dimensions?
Using the matrix method in higher dimensions, making the scaling vector is easy, and making the unrotation matrix is still just taking the transpose of the rotation matrix. It’s more difficult making the rotation matrix though.
In 3D, the scaling direction will be a 3D vector, and you need to come up with two other vectors that are perpendicular to that scaling direction. One way to do this could be to take any vector which is different from the scaling vector, and cross product that with the scaling vector. That will give you a vector perpendicular to both, and you can take that as your second vector. To get the third vector, cross product that vector with the scaling vector. You will then have 3 perpendicular vectors, and an orthonormal basis that you can use to fill out your rotation matrix. The first column is the scaling vector, the second column is the second vector found, and the third column is the third column found.
The cross product only exists in the 3rd and 7th dimension though, so if you are working in a different dimension, or if you don’t want to use the cross product for some reason, another way you can make an orthonormal basis is by using the Gram-Schmidt process. There’s a great video on it here: https://www.youtube.com/watch?v=KOkuTXrv5Gg
For the vector project method, you also need the orthonormal basis vectors to do all the vector projections, before you scale the x axis, and then re-combine the projections, so it boils down to the same issues as the matrix method.
This has methods for generating a high order rotation that moves a point to a particular axis. There is rarely a need need for a Gram Schmidt process and computing the matrix can be made quite cheap 🙂
I think the most efficient way to find a rotation matrix that takes a unit vector A and moves it to another unit vector B (in any dimension) is to find the find the reflection matrix that maps A to C (where C=B with a single component negated – doesn’t matter which one) and then flip the sign of the corresponding row of the matrix to turn it into a rotation.
Finding a reflection matrix that does this requires only a single division in an efficient implementation for any dimension.
if your use case doesn’t need accuracy for scaling amounts near zero, method 2 has a variant that saves you from having to find perpendicular vectors: point + (scaling amount – 1) * dot(normalized scaling vector, point) * normalized scaling vector.
Shuffle Iterator – This means you have an iterator that can walk through the items in a shuffle without actually shuffling the list. It’s seedable too, to have different shuffles.
Low Discrepancy – This means that if you are doing something like numerical integration on the items in the list, you’ll get faster convergence than using a white noise random shuffle.
Random Access – You can ask for any item index in the shuffle and get the result in constant time. If you are shuffling 1,000,000 items and want to know what item will be at the 900,000th place, you don’t have to walk through 900,000 items to find out what’s there. It is just as happy telling you what is at the 900,000th place in the shuffle, as what is at the 0th place.
Inversion – You can also ask the reverse question: At what point in the shuffle does item 1000 come up? It also gives this answer in constant time.
After writing that post, I spent some time trying to find a way to do the same thing in two dimensions. I tried some things with limited success, but didn’t find anything worth sharing.
A couple people suggested I try putting the 1D shuffle iterator output through the Hilbert curve, which is a way of mapping 1d points to 2d points. I finally got around to trying it recently, and yep, it works!
You can also use the Z-order curve (or Morton curve) to do the same thing, and it works, but it doesn’t give as nice results as the Hilbert curve does.
People also probably wonder how Martin Robert’s R2 sequence (from https://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/) would work here, since it generalizes the golden ratio to 2d, and the 1d shuffle iterator is based on the golden ratio, but I couldn’t find a way to make it work. For example, multiplying the R2 sequence by the width and height of the 2D grid and casting to integer causes the same 2d location to be chosen multiple times, which also leaves other 2d locations never chosen. That is fine for many applications, but if you want a shuffle, it really needs to visit each item exactly once during the shuffle. Below is how many duplicates that had at various resolutions. Random is also included (white noise, but not a shuffle) to compare to.
Texture Size
Pixel Count
R2 Duplicate Pixels
R2 Duplicate Percent
Random Duplicate Percent
64×64
4,096
847
20.68%
36.89%
128×128
16,384
3,379
20.62%
37.01%
256×256
65,536
16,375
24.99%
36.78%
512×512
262,144
29,203
11.14%
36.80%
1024×1024
1,048,576
392,921
37.47%
36.78%
2048×2048
4,194,304
1,831,114
43.66%
36.78%
So first up, here’s an RMSE (root mean squared error) graph, integrating an image in the repo cabin.png (aka finding the average pixel color). RMSE is averaged over 1000 tests to reduce noise in the plot, where each test used a different random seed for the shuffles. White uses a standard white noise shuffle. Hilbert and ZOrder use a 1D low discrepancy shuffle, then use their respective curves to turn it into a 2D point.
Here is the 512×512 image being integrated:
So that’s great… Hilbert gives best results generally, but Z Order also does well, compared to a white noise shuffle.
Some things worth noting:
Both Hilbert and Z order curves can be reversed – they can map 1d to 2d, and then 2d back to 1d. That means that this 2d shuffle is reversible as well. To figure out at what point in a shuffle a specific 2D point will appear, you first do the inverse curve transformation to get the 1D index of that point, and then ask the low discrepancy shuffle when that 1D index will show up.
As I’ve written it, this is limited to powers of 2. There are ways of making Hilbert curves that are not powers of 2 in size but I’ll leave that to you implement 😛
I did omit something from the previous test though… What if we didn’t put the 1d shuffle through any curve? What if we treated the image as one long line of pixels (like it is, in memory) and used a 1d shuffle on that? Well, it does pretty well!
However, you should see the actual point sets we are talking about before you make a final judgement.
Below are the White noise shuffled points (left), Z order points (middle) left, Hilbert points (middle right) and 1DShuffler points (right) at different numbers of points, to show how they fill in the space. The last image in each uses a point’s order as a greyscale color to show how the values show up over time.
Click on the images to see them full sized.
64×64
128×128
256×256
512×512
1024×1024
So, while the 1DShuffler may give better convergence, especially near the full sample count, the point set is very regularly ordered and would make noticeable rendering artifacts. This is very similar to the usual trade off of low discrepancy sequences converging faster but having aliasing problems, compared to blue noise which converges more slowly but does not suffer from aliasing, and is isotropic.
You might notice a strange thing in Hilbert at 1024×1024 where it seems to fill in “in clumps”. 2048×2048 seems to do the same, but 4096×4096 goes back to not clumping. It’s strange and my best guess as to what is going on is that Hilbert mapping 1d to 2d means that points nearby in 1d are also nearby in 2d, but the reverse does not hold. Points far away in 1d are not necessarily far away in 2d? I’m not certain though.
Bonus Convergence Tests
Here’s a 512×512 white noise image to integrate, and the convergence graph below.
Here is the same with a blue noise texture. This blue noise made with FAST (https://github.com/electronicarts/fastnoise). This is the same pixel values as the white noise texture above, just re-arranged into a blue noise pattern.
So interestingly, integrating the white noise texture made the “good samplers” do less good. Integrating the blue noise texture made the “good samplers” do even less well and be equal to white noise sampling.
What gives?
This isn’t proven, and isn’t anything that I’ve seen talked about in literature, but here’s my intuition of what’s going on:
“Good sampling” relies on small changes in sampling location giving small changes in sampled value. This way, when it samples a domain in a roughly evenly spaced way, it can be confident that it got fairly representative values of the whole image.
This is true of regular images, where moving a pixel to the left or right is going to usually give you the same color.
This may or may not be true in a white noise texture, randomly.
A blue noise texture however, is made to BREAK this assumption. Small changes in location on a blue noise texture will give big changes in value.
Just a weird observation – good sampling tends to be high frequency, while good integrands tend to be low frequency. Good sampling tends to have negative correlation, while good integrands tend to have positive correlation.
I haven’t made many blog posts this year, due to working on a project at work that I just got open sourced. The project is called Gigi and it is a rapid prototyping and development platform for real time rendering, and GPU programming. We’ve used it to prototype techniques and collaborate with game teams, we’ve used it in published research (FAST noise), and used it to generate code we’ve shared out into the world.
I also used Gigi to do the experiments and make the code that go with this blog post. Both the Gigi source files and the C++ DX12 code generated from those Gigi files are on github at https://github.com/Atrix256/BNOctaves.
Nikita tagged me in this interesting post with a neat looking screenshot. I thought this was interesting and wanted to look deeper so got the details.
The idea is that you start with a blue noise texture. You then add in another blue noise texture, but 2x bigger and multiplied by 1/2. You then add in another blue noise texture again, but 4x bigger, and multiplied by 1/4. Repeat this as many times as desired. While doing this process you sum up the weights (1 + 1/2 + 1/4 + …) and you divide the result by the weight to normalize it.
Results – Blue Noise
Here is a blue noise texture by itself, made using the FAST command line utility (https://github.com/electronicarts/fastnoise). The upper left is the noise texture. The upper right is the DFT to show frequency content, and the bottom half is a histogram. Here we can see that the texture is uniform blue noise. The histogram shows the uniform distribution, and the DFT shows that the noise is blue, because there is a dark circle in the center, showing that the noise has attenuated low frequencies.
With 2 octaves shown below, there are lower frequencies present, and the histogram gets a hump in the middle. The low frequencies increase because when we double the size of the texture, it lowers the frequency content. The reason the histogram is no longer uniform is because of the central limit theorem. Adding random numbers together makes them follow a binomial distribution, and Gaussian at the limit (more info: https://blog.demofox.org/2019/07/30/dice-distributions-noise-colors/).
With 3 octaves shown below, the effects are even more pronounced.
So yeah, adding octaves of blue noise together does produce blue noise! It makes the distribution more Gaussian though, and reduces the frequency cutoff of the blue noise.
Note: It doesn’t seem to matter much if I use the same blue noise texture for each octave, or different ones. You can control this in the demo using the “DifferentNoisePerOctave” checkbox.
Results – Other Noise Types
After going through the work of making this thing, I wanted to see what happened using different noise types too.
White Noise
Here is 1 octave of white noise, then 3 octaves. White noise starts with randomized but roughly equal levels of noise in each frequency. This process looks to make the low frequencies more pronounced, which makes sense since the octaves get bigger and so are lower frequency.
3×3 Binomial Noise
This noise type is something we found to rival “traditional” (gaussian) blue noise while writing the FAST paper. Some intuition here is that while blue noise has a frequency cutoff equal distance from the center (0hz aka DC), this binomial noise frequency cutoff has roughly equal distance from the edges (aka nyquist). More info on FAST here https://www.ea.com/seed/news/spatio-temporal-sampling.
Below is 1 octave, then 3 octaves. It seems to tell the same story, but the 3 octaves seem to have some concentric rings in the DFT which is interesting.
3×3 Box Noise
Here is noise optimized to be filtered away using a 3×3 box filter. 1 octave, then 3 octaves.
5×5 Box Noise
Here is noise optimized to be filtered away using a 5×5 box filter. 1 octave, then 3 octaves.
Results – Low Discrepancy Noise Types
Here are some types of low discrepancy noise.
Interleaved Gradient Noise
Interleaved gradient noise is a low discrepancy object that doesn’t really have a name in formal literature. I call it a low discrepancy grid, which means we can use it as a per pixel source of random numbers, but it has low discrepancy properties spatially. It was invented by Jorge Jimenez at Activision and you can read more about it at https://blog.demofox.org/2022/01/01/interleaved-gradient-noise-a-different-kind-of-low-discrepancy-sequence/.
Here is 1 octave then 3 octaves. The DFT goes from a starfield to a more dense starfield.
// R2 Low discrepancy grid
// A generalization of the golden ratio to 2D
// From https://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/
float R2LDG(uint2 pos)
{
static const float g = 1.32471795724474602596f;
static const float a1 = 1 / g;
static const float a2 = 1 / (g * g);
return frac(float(pos.x) * a1 + float(pos.y) * a2);
}
1 octave, then 3 octaves. The DFT goes from a starfield, to a swirly galaxy. That is pretty cool honestly! The histogram also have a lot less variance than other noise types, which is a cool property.
Wait… Isn’t This Backwards?
When we make Perlin noise, we start with the chunky noise first at full weight, then we add finer detail noise at less weight, like the below with 1 and 8 octaves of perlin noise.
In our noise experiments so far, we’ve been doing it backwards – starting with the smallest detailed noise at full weight, then adding in chunkier noise as lower weight. What if we did it the other way? The demo has a “BlueReverse” noise type that lets you look at this with blue noise. Here is 1, 3 and 5 octaves:
It has an interesting look, but does not result in blue noise!
Use of Gigi – The Editor
I used Gigi to make the demo and do the experiments for this blog post. You can run the demo in the Gigi viewer, or you can run the C++ dx12 code that Gigi generated from the technique I authored.
I want to walk you through how Gigi was used.
Below is the technique in the Gigi editor. The “CS: Display” node makes the noise texture with the desired number of octaves. “Sub: Histogram” makes the histogram graph. “Sub: DFT” makes the DFT. “CS: MakeCombinedOutput” puts the 3 things together into the final images I pasted into this blog post.
If you double click the “Sub: Histogram” subgraph node, it shows you the process for making a histogram for whatever input texture you give it. This is a modular technique that you could grab and use in any of your Gigi projects, like it was used here!
If you double click the “Sub: DFT” subgraph node instead, you would see this, which is also a modular, re-usable technique if ever you want to see the frequency magnitudes of an image. Note: This DFT is not an FFT, it’s a naive calculation so is quite slow.
Use of Gigi – The Viewer
If you click “Open Viewer” in the main technique (BNOctaves.gg) it will open it in the viewer. You can change the settings in the upper left “Variables” tab and see the results instantly. You can also look at the graph at different points in the execution on the right in the “Render Graph” tab, and can inspect the values of pixels etc like you would with a render doc capture. A profiler also shows CPU and GPU time in the lower left. (I told you the DFT was slow, ha!). For each image above in this blog post, i set the parameters like I wanted and then clicked “copy image” to get it to paste into this blog post.
Use of Gigi – The Compiler & Generated Code
I used GigiCompiler.exe to generate “DX12_Application” code for this technique and the results of that are in the “dx12” folder in the repo. You run “MakeSolution.bat” to have cmake create BNOctaves.sln and then can open and run it. Running that, you have the same parameters as you do in the viewer, including a built in profiler.
I did have to make one modification after generating the C++ DX12 code. There are “TODO:”s in the main.cpp to help guide you, but in my case, all I needed to do was copy the “Combined Output” texture to the render target. That is on line 456:
// TODO: Do post execution work here, such as copying an output texture to g_mainRenderTargetResource[backBufferIdx]
CopyTextureToTexture(g_pd3dCommandList, m_BNOctaves->m_output.texture_CombinedOutput, m_BNOctaves->m_output.c_texture_CombinedOutput_endingState, g_mainRenderTargetResource[backBufferIdx], D3D12_RESOURCE_STATE_RENDER_TARGET);
Closing
I sincerely believe that “the age of noise” in graphics has just begun. I more mean noise meant for use in sampling or optimized for filtering, but that type of noise may have use in procedural content generation too, like this post explored.
I am also glad that Gigi is finally open sourced, and that I can use this absolute power tool of research in blog posts like this, while also helping other people understand the value it brings. Hopefully you caught a glimpse of that in this post.
The C++ code that goes with this blog post can be found at https://github.com/Atrix256/GoldenRatioShuffle and includes a header file that contains a shuffler that should be easy to drop into other C++ code.
A shuffle iterator lets you iterate through a shuffle without actually doing the shuffle. It can provide the next item when asked while only storing your current index into the shuffle and the random seed. A friend and I wrote an article showing how to achieve that using mathematics and cryptography: https://www.ea.com/seed/news/constant-time-stateless-shuffling
However, the shuffle made with that technique is a “white noise shuffle”. That is, each item in the shuffle has equal probability of being anywhere in the shuffle and each item in the shuffle an independent random number (uniform random selection without replacement). This is useful anywhere a shuffle is used, but in many situations, the algorithm using the shuffled data could perform better if it could use something better than white noise. I recorded a video that talks about the benefits of going beyond white noise here: https://www.youtube.com/watch?v=tethAU66xaA
In this post, I’ll show how to make a stateless constant time shuffle iterator which:
Takes an irrational number to drive the shuffle sequence
Takes a random seed
Supports random access
Supports inversion
Each briefly explained:
Takes an irrational number to drive the shuffle sequence
The golden ratio is best, but if you need multiple shuffles, having two shuffles based on the same irrational number will cause them to be correlated. Like imagine two sine waves, but starting at different angles (phases / offsets). That correlation will show up in what you are using the shuffles for and be a problem. Using different high quality irrational numbers for each shuffle will make them be more uncorrelated, depending on the numbers used. See the previously mentioned youtube video for more details about high quality vs low quality irrational numbers. Note that this also touches on something I call “co-irrational numbers” that you can read about here: https://blog.demofox.org/2020/07/26/irrational-numbers/
The golden ratio is the best choice, but the square root of 2 is also good, as is the square root of 3. Pi however is a very bad choice!
Takes a random seed
The shuffle in this article is a deterministic sequence that gives every item [0,N) exactly once before repeating, and does so in a quasirandom order.
You can imagine that you could add a constant to every item in that sequence and apply modulo N to bring it back to being in [0,N). That constant is the random seed and that is how we apply the seed. So, for a list of N items long, the random seed can be any value between 0 and N-1.
Supports Random Access
You can either ask for the next item in the list, or you can ask for the Mth item in the list.
That is, you can have a shuffle with 4 billion items in it and ask “what is the 1 billionth item in the shuffle?” and it will tell you in constant time – taking just as much computing power / time as if you asked what the 1st item in the shuffle was.
Supports Inversion
You can ask when a specific item will show up in the shuffle.
Again, if you are shuffling 4 billion items, you can ask it “At what point in the shuffle does the 5th item appear?” and it will give you the answer, in constant time again.
Let’s dive into the details.
Low Discrepancy Sequence to Low Discrepancy Shuffle
A well known and powerful one dimensional low discrepancy sequence (LDS) is the golden ratio LDS, which is calculated like this in C++:
The first time you call that function, you can pass any value in between 0 and 1 as the “initial state” or “random seed”. Every time you call it after that, it returns a value between 0 and 1 that you can use as a “random number” in your code. This works for numerical integration (quasi monte carlo integration), but also improves things like calculating treasure drops when a player kills a monster in a video game – see the previously mentioned youtube video for more info.
If you wanted to use these values in [0,1) to shuffle N items, you might think to multiply those values by N and clamp to [0,N-1), or something similar.
That works pretty well, but there are repeats! With an initial starting state of 0.55, here are 10 items shuffled:
5 1 7 4 0 6 2 8 4 1
That sequence has 1 and 4 twice, and misses 3 and 9.
There is a neat trick from Marc Reynolds to convert the golden ratio into an integer for when you want integer sequences (http://marc-b-reynolds.github.io/distribution/2020/01/24/Rank1Pre.html). He does it by multiplying the golden ratio by the number of integers there are and rounding to odd (round up or down, to whichever is the odd number). This treats the [0,N) integer range as if it were the real number range [0,1). It’s basically changing from working in floating point, to fixed point, and gives better results than using the golden ratio with conversion to integers.
That is nice, but doing that, it can still have repeats. Rounding to odd helps it have a longer period sequence, but it is not always a maximal period sequence. To make a maximal period sequence, we would need our integer approximation of the golden ratio to be coprime to the number of items in our shuffle.
We can search for such a number by checking if the number is coprime. If not, look at the number -1, then +1, then -2, then +2 and so on. We will eventually find the number coprime to N that most closely approximates the golden ratio. The farther we get from the golden ratio, the more our low discrepancy properties are degraded, so we want to be as close as possible.
Here’s the graph of the real number GoldenRatio*N (blue) and the integer coprime found (red) up to N=1000, you can see the coprime found is very close.
Here is a graph of the difference between the two. It never goes above 3.5.
So to make a low discrepancy shuffle iterator happens in two steps.
If we have an irrational number A, and want to shuffle N items, we find the integer coprime to N, which is closes to the real number (A*N). We’ll call this P.
Once we have P, we can call this function iteratively to get the next item in the shuffle:
unsigned int Shuffle_NextValue(unsigned int lastValue, unsigned int P, unsigned int N)
{
return (lastValue + P) % N;
}
To give this shuffle a random seed, the first time you call this function, give it a lastValue which is a random integer in [0,N).
Let’s look again at the golden ratio sequence before of 10 items which started at index 5 (a starting value of 0.548813522 for the LDS) and had some duplicates and missing items. We shuffled N=10 items, which gives us a coprime P=7, and we’ll start at index 5.
Sequence
Duplicates
Missing
LDS
5 1 7 4 0 6 2 8 4 1
1 4
3 9
Coprime
5 2 9 6 3 0 7 4 1 8
–
–
The coprime approximation of the golden ratio gets rid of duplicates / missing values, while still being a nice low discrepancy sequence. Neighboring values are “very different” from each other.
Here is a 64 item shuffle with coprime P=41. The index started at 26, and the LDS started at 0.417021990.
The code that goes with this post lets you see the error rate of LDS at different item counts and with different seeds, but the number of missing items in the LDS version of the shuffle is usually around 10-20%. The coprime version never misses any items, and is always pretty close to the integer (fixed point) approximation of the golden ratio.
Random Access
The code for the random access shuffle iterator is shown as an iterative process: you start at the first item and call Shuffle_NextValue() to get the next item. Rinse and repeat until you call it N times and the shuffle starts over.
The shuffle iterator does support random access though. If you want to know what item is at a specific shuffle index, the result is (index*P)%N. if index is a large enough number, you can start to hit problems with overflow though. To help this a bit, you can first modulus index by N, since the shuffle repeats itself every N steps.
That gives you this following code, which does give random access, but can still have integer overflow problems if index and N are large enough. So, be careful and take extra precautions when needed!
unsigned int Shuffle_RandomAccess(
unsigned int index,
unsigned int shuffleSeed,
unsigned int P,
unsigned int N)
{
return (((index % N) * P) + shuffleSeed) % N;
}
Inversion
If we want to know what location a specific item will show up in the shuffle, we can do that. Ignoring the integer overflow problems, and random seed for a moment, the formula for our random shuffle is this:
If we know item, P and N, and we want to know index, we are trying to calculate the “Modular Multiplicative Inverse”, which you can do using the “Extended Euclidean Algorithm”. It’s an iterative algorithm, but is not complicated or super expensive to calculate.
The extended euclidean algorithm tells you what you need to multiply P by to get 1. it gives you the x below.
Once you know how many steps it takes to get to 1, you can double it to get to 2, or triple it to get to 3 and so on. So if you want to know how many times it takes adding P to 0 to get to a specific number, you just multiply x by that number. That number may be much larger than N but you can just put it through modulo N since the sequence repeats every N times.
To handle the seed, you need to subtract out how many times you add P to 0 to get to the seed value.
In the end, the inversion algorithm looks like this:
Notice on the last step we are adding an extra numItems before subtracting stepsToSeed. This is to handle the case where stepsToSeed is larger than stepsToValue and prevent integer underflow.
If you needed to do inversion often, you could store off that t value to make inversion lightning fast, without having to run the extended Euclidean algorithm each time.
To help show why this shuffle is meaningful, we take 10,000 items valued 0 to 9,999, divide them by 9,999 to make them be between 0 and 1, and we iteratively average them from index 0 to index 1, index 2, index 3 and on up to index 9,999.
We do this for both a white noise shuffled array, and also this post’s low discrepancy shuffle iterator. We do that 1000 times and collect the average error.
Below is the graph of that data, which shows the low discrepancy shuffle iterator (blue) far outperforms the white noise shuffled array (orange). It gives a much better average much sooner which shows that it’s low discrepancy and has the desirable numerical properties we would expect from a golden ratio sequence, while still being a shuffle.
Hopefully you enjoyed this post. If you end up using this shuffle iterator for anything and see benefits, I’d love to hear about it here, or on mastodon, at https://mastodon.gamedev.place/@demofox
The C++ code that goes with this blog post can be found at https://github.com/Atrix256/GoldenRatioShuffle and includes a header file that contains a shuffler that should be easy to drop into other C++ code.
Here it is as shadercode (hlsl) as well:
// Example calculations
// 128x128 is 16384 pixels to shuffle, so numItems = 16384
// Coprime is 10127, calculated by code at that repo to be the coprime value closest to the golden ratio.
// StepsToUnity is 5999, calculated by code at that repo, using the chinese remainder theorem.
// Seed is the random seed of the shuffle, such as the constant "435".
// Return shuffle[index]
uint LDSShuffle1D_GetValueAtIndex(uint index, uint numItems, uint coprime, uint seed)
{
return ((index % numItems) * coprime + seed) % numItems;
}
// return the index where shuffle[index] == value
uint LDSShuffle1D_GetIndexOfValue(uint value, uint numItems, uint stepsToUnity, uint seed)
{
uint stepsToValue = (stepsToUnity * value) % numItems;
uint stepsToSeed = (stepsToUnity * seed) % numItems;
return (stepsToValue + numItems - stepsToSeed) % numItems;
}
Disclaimer: This is based on information I believe to be true, in 2024, in the United States, and in the state of California specifically. Laws change over time and there are subtleties to all these things, so please talk to a financial professional before taking specific actions. I do hope this helps illuminates how to get the most out of these benefits though.
If you are a programmer (and mhigherany other fields), you are likely to encounter an Employee Stock Purchase Plan (ESPP) at your work. If you work in an in demand area, you are also likely to encounter Restricted Stock Units (RSUs) as well. This article aims to help illuminate some things I wish I knew / realized earlier in my career.
Let’s talk about how these work.
How Do ESPPs Work?
ESPPs let you deduct money out of your pay check automatically, into an account, where periodically (like once a quarter) it is used to buy company stock. There is commonly a “look back” window of time where the price you pay per share of the stock is the lowest price of the stock in that window of time. Also commonly, the stock is offered at a discount on top of that, say 5%, 10% or 15% off.
The money taken out of your paycheck is post-tax, meaning you pay income tax on the money you use to buy the stock.
When you sell the stock, the amount of discount you had on buying the shares become income, which follows regular tax income rules for both federal and state taxes. If the stock has gone up, capital gains applies to the increased value.
Example:
Let’s say you had 1700$ in your account when it was time for the stock buy to happen. Let’s say the current price of the stock was 130$, but that there was a 6 month “look back” window, and the lowest price of the stock in the last 6 months was 100$. You get to purchase stock at 100$, instead of at the 130$ it currently costs. Buying the stock already gives you a profit of 30$ a share! On top of this, let’s say you get a 15$ discount. That 15% discount applies to the 100$ price, meaning you get to buy shares at 85$ instead of the current market price of 130$. That allows you to buy 20 shares of stock for that 1700$, instead of the 13 you could buy at the current market price.
Let’s say you hold onto the stock for a while and it goes up to 250$ a share, and then you decide to sell.
Since you bought 20 shares of stock at 85$, when the stock was actually 130$ a share, that discount of 45$ per share becomes taxable income. 45$ * 20 shares = 900$ of increased taxable income on your tax return, to perhaps push you into a different tax bracket.
However, the stock also went up 120$ a share. 120$ * 20 shares = 2400$ of “capital gains”. In the united states, for federal taxes, if you’ve held this stock for less than 1 year, you pay “short term capital gains” tax, which means the capital gains are treated as regular income, that can push you into a different tax bracket and all of that. If you’ve held onto it for longer than a year, you pay “long term capital gains” which is taxed based on your income, but most folks will be in the 15% bracket.
In the state of California, there is no distinction between long and short term capital gains. The gains all become taxable income.
RSU stands for “Restricted Stock Unit” and they are stock that is promised to you that vests (becomes yours) over a specified schedule. For instance, you may be given 480 RSUs when you are hired, with 120 vesting (being given to you) at your 1 year anniversary, and 30 more every quarter for the 3 years after that until all are 480 shares are vested.
When the stock vests, it becomes yours, but the value of the stock becomes taxable income that you pay that year. Whenever you sell the shares and if the value has gone up, you then have to pay short term gains if you’ve held the stock less than a year, or long term if you’ve held it more than a year. Same as ESPPs.
Lets look at some motivations of these benefits. Reality is complicated per usual, with a mix of both optimistic and cynical truths.
ESPP
An optimistic viewpoint of ESPPs is that it’s a way to invest in the company because you believe it’s going to do well. The company thinks so too and is helping you out by giving these shares to you at a discount. How awesome!
A more cynical viewpoint is that any money a company pays an employee that goes into the ESPP program never actually left the company. It can effectively reduce payroll costs, turning some of those guaranteed payroll funds into not guaranteed stock holder returns, where the employee takes on the risk. If the company does well, yeah, the employees with ESPP will cash out some, but there will be extra money around to fund that. If the company does poorly, the employees with ESPP will lose, and get their money back at a reduced rate, which is great for the company, effectively having to pay less payroll if things are going poorly.
RSUs
An optimistic viewpoint here is that a company is going to do well, so giving you stock is going to benefit you in the future when the stock goes up.
A reality of the vesting schedule is that you are less likely to leave if you have significant amount of stock that hasn’t vested yet. It helps retention.
A more cynical viewpoint is that RSUs are better for the company than a cash bonus, because the money is still invested in the company (they didn’t actually have to pay anything out). Also, you are given shares, not dollar amounts, so the risk of how much that is actually worth goes to you, the employee.
Bottom Line
The bottom line is that giving stock is a way of not actually giving anything in the short term, and putting any longer term risk on you, the employee. It’s a pretty clever business move frankly. Despite this, these benefits are still beneficial to you – you will get that money eventually!
How Can I Take Most Advantage of These Benefits?
So, companies definitely have an angle on these benefits and are certainly going to do everything in their power to get the most out of it from their side. They are likely to offer educational material that lists the benefits that help them, but not say the things I’m about to say below. So how can you make the most of it for yourself?
ESPPs
If you believe the stock is going to go up, and you want the stock, go ahead and buy the stock! Max out the amount you contribute if you can afford it and sell when you are ready.
If you don’t believe the stock is going to go up, another way of thinking about ESPPs is that they are allowing you to buy money at a discount. If your ESPP program gives you a 15% discount on the stock, that means you only pay 85 cents for every dollar of stock you buy. This stock is paid for with money you are already paying income tax on so buying this stock doesn’t do anything special for you as far as taxes are concerned. You are welcome to then IMMEDIATELY sell the stock, right after purchase and get your dollar back for every 85 cents you paid. It’s true that you then have to take the 15 cents on the dollar as income and pay income tax on it, but that is the same as if you got a raise. If you immediately sell your ESPP stock after buying it, you are essentially giving yourself a raise. Why wouldn’t you do that? And why wouldn’t you max out how much you were contributing to ESPP to max out your raise?
Something to watch out for though is that if the stock is sold to you a lot lower because of the “look back” window, the difference between what you bought it for and what you sold it for becomes “capital gains”, with the short term vs long term implications I explained above. It might be that you prefer to hold onto that stock for a year before selling, to pay a smaller amount of tax on it. If this issue comes up, it’s a good problem to have though. The problem is basically “What do I do with all this money?!”
In short, if you don’t want the stock, and the stock is flat, or decreasing in value, you can max out your ESPP, and sell immediately after the buying event, to give yourself a raise.
RSUs
Again, if you believe the stock is going to go up, and you want the stock, go ahead and keep these RSUs.
If you don’t believe the stock is going to go up, what should you do?
When your RSUs vest, you pay income tax on the value they are at the time of vesting. This is equivalent to the company paying you more, and you go out and buy stock with that extra pay. If this wasn’t your plan for the money, sell it right when it vests! Convert it back into cash. The stock isn’t going to have changed price (except normal market fluctuations), so you aren’t going to have any gains. Take the cash and walk away, as if the stock purchase never happened. This is fine and doesn’t affect your taxes at all.
Furthermore, let’s say you are working at a company where you either don’t have any more RSUs waiting to vest, or you only have a small amount remaining. If you are in a position where you were given RSUs, you can likely go to a competitor, and get more RSUs from them as a signing bonus. Even if you have the same title, job function, and pay, go get some RSUs from the competition and look out for #1 (yourself). When a company is interested in retaining employees, they will offer them more RSUs to keep this vesting train going. If your company is not doing that for you, they are either not very smart, or they are not interested in keeping you. Either way, it might be a good time to jump elsewhere 🙂
Watching The Stock Market
When you have a lot of money in the stock market, it can be nerve wracking. Having “decent fractions of your pay” in the stock market is a lot of money IMO. You might feel like you need to be constantly watching the market to protect your money, in case there is a big spike up or a big spike down.
A nice thing to help this is there are things called “stop orders” and “limit orders” that can help automate some simple tasks to help you not be constantly glued to stock tickers. There are also more complicated things such as “stop limit” orders.
For instance, you can say “Sell this much stock if the price drops below 400$ a share” with a stop order. Sometimes stock can be volatile, and it could drop down to 200$ a share for a moment, then jump back up to 450$ a share. With a stop order, your stock might sell during that moment when it dropped down to 200, and you get sort of screwed, selling your stock too cheaply and losing money. A stop limit order lets you set a minimum price so you can say “Sell this much stock if the price drops below 400$ a share, but don’t sell if it goes below 350$ a share”.
With something like that in place, you can happily go off and have lunch without having to check your phone every few minutes just in case LOL.
ESPP and RSUs can be nice perks when working for companies, but it’s important to know how these mechanisms work, to help make the best decisions for yourself and your own situation.
If you want the stock, by all means, keep it! Otherwise…
ESPPs can be a way to give yourself a raise. You deserve it, give yourself that raise!
RSUs are as if you were given extra pay, but the company bought stock with it before they gave it to you. Sell that stock if you weren’t going to buy stock with that money anyways, and take that cash, for the same tax implications.
I recently published a set 3 of articles on machine learning in plain C++ for game developers, as well as a video version of the same content. This article explains an improvement to plain gradient descent called Adam and can be read as a follow up to those articles.
Let’s say we wanted to find the values of the x parameters in the function below that gives the lowest y value:
This is an optimization problem and functions of this form are called cost functions. Gradient descent is a method for doing this optimization, where you start with a random set of x values, and “roll down hill” by making small adjustments to those x values that lead to lower y values. You can stop either after a specific number of steps, or when the y value is below some threshold.
The gradient of a function is a vector of the partial derivatives of the x parameters, and the vector as a whole points the steepest direction up hill on the surface of the function. If you want to go the steepest direction down hill, you just go in the opposite direction.
A challenge is that the steepest path down hill is most often a curved path, not a straight line. To deal with this, you can take small steps down hill and re-evaluate the gradient at each step to figure out the direction to step in next. In machine learning, this step size is called the learning rate and is often called α (“alpha”) and is tuned by hand.
Smaller steps give you a more accurate result, but may be computationally wasteful, taking many small steps when you could take fewer, larger steps towards the goal. This makes for longer than necessary computation (or training) times.
Larger steps might be able to get you nearer to the goal faster, but if the minimum value of the function lies in a deep, narrow canyon, large steps won’t be able to get you into the deepest part of the function.
Having to adjust this learning rate by hand, and having the appropriate learning rate change as you do gradient descent is one of the two problems Adam addresses, by turning down the learning rate when the terrain gets more complicated / bumpy.
Another challenge with gradient descent is that it’s possible to get stuck in a shallow valley when there is a much deeper canyon waiting to be found. This is called a local minima, because while it’s the deepest point locally, it isn’t the globally deepest part of the function, or the global minimum.
A common way of dealing with this is to add momentum to our imaginary ball rolling down hill, so that it can escape local minima and hopefully find deeper, more global minima. If not, it should just roll back into the minima it did find. Momentum is the second benefit Adam gives over plain gradient descent.
Adam is a popular, and simple, method of improving gradient descent and is heavily used in machine learning, which uses gradient descent to train neural networks.
Adam
Adam has 4 parameters, but only alpha usually needs to be tuned.
Alpha – This is essentially the learning rate.
Beta1 – The rate of exponential moving average of the derivative (First Moment). Default is 0.9.
Beta2 – The same for the derivative squared (Second Moment). Default is 0.999.
Epsilon – a value to prevent division by zero. Default is 1e-8.
Adam works by keeping an exponential moving average of the derivative (the first moment of the derivative), as well as a moving average of the derivative squared (the second moment of the derivative), each with their own smoothing factor. In this way, this gives Adam an imperfect memory of the derivatives it has seen, tending to forget older derivative values, and weight newer derivative values more highly. The squared derivative moving average is also more sensitive to large outliers.
The code for doing this is straight forward:
static const float c_beta1 = 0.9f;
static const float c_beta2 = 0.999f;
m = c_beta1 * m + (1.0f - c_beta1) * derivative;
v = c_beta2 * v + (1.0f - c_beta2) * derivative * derivative;
Interestingly, the m calculation is the same as doing a lerp from m to derivative with a t value of 0.1. If you work in real time rendering, that should look familiar, as that is how TAA (Temporal Anti Aliasing) works as well, also often using a value of 0.1. v works similarly, but uses a t value in the lerp of 0.001.
So, m and v are a moving average of recent derivative and derivative squared values respectively.
The learning rate alpha is then multiplied by m and divided by the square root of v, to calculate an adjusted derivative to use for gradient descent.
float adjustedDerivative = alpha * m / sqrt(v);
if the derivative was a constant, m would be that constant value, and so would sqrt(v), which means alpha would be the amount to adjust the parameter by. if m was bigger than sqrt(v), it would increase the learning rate, else if m is smaller than sqrt(v), it would decrease the learning rate. Through this simple operation, Adam gives both momentum to escape local minima, as well as the ability to slow down when there are more details to explore on the surface of a function.
We aren’t quite done though. If you work in graphics, you likely have noticed that when using TAA, the first several frames rendered are darker than they should be, before enough frames have been processed. This same problem occurs in Adam, meaning that when it’s first starting out, both m and v are biased towards zero. Adam does bias correction to account for this, essentially making the numbers larger (or for TAA, it’d make the pixels brighter). It calculates an “m hat” and a “v hat”, and that is what is used to calculate the learning rate instead of the raw m and v. Since v hat could possibly be zero, we have to account for that in the division.
Putting it all together, here is a small struct you can use to put a derivative through Adam to get an adjusted derivative. You would have one of these per parameter to your function, so that Adam was able to adjust the learning rate of each parameter individually.
// Note: alpha needs to be tuned, but beta1, beta2, epsilon as usually fine as is
struct Adam
{
// Internal state
float m = 0.0f;
float v = 0.0f;
// Internal state calculated for convenience
// If you have a bunch of derivatives, you would probably want to store / calculate these once
// for the entire gradient, instead of each derivative like this is doing.
float beta1Decayed = 1.0f;
float beta2Decayed = 1.0f;
float GetAdjustedDerivative(float derivative, float alpha)
{
// Adam parameters
static const float c_beta1 = 0.9f;
static const float c_beta2 = 0.999f;
static const float c_epsilon = 1e-8f;
// exponential moving average of first and second moment
m = c_beta1 * m + (1.0f - c_beta1) * derivative;
v = c_beta2 * v + (1.0f - c_beta2) * derivative * derivative;
// bias correction
beta1Decayed *= c_beta1;
beta2Decayed *= c_beta2;
float mhat = m / (1.0f - beta1Decayed);
float vhat = v / (1.0f - beta2Decayed);
// Adam adjusted derivative
return alpha * mhat / (std::sqrt(vhat) + c_epsilon);
}
};
An interesting aspect of Adam is that you can see that it adjusts the step size (alpha) by saying “be bigger” or “be smaller” but it never gives a specific value that is an optimal step size, and you still have to tune the learning rate value by hand. To me, that looks like something we’ll likely have a better answer to in the future.
I also find it interesting that adam calculates an adjusted learning rate (alpha) and doesn’t multiply the derivative by that, but just uses that for the derivative instead. That makes it so while the length of a function’s gradient may change from location to location, this actual step size is more fixed, and is just affected by how the gradient changes over time.
Another thing people do with gradient descent is start with a larger learning rate, and decrease it over time. This helps the gradient descent move quickly in the beginning to get closer to the global minimum, and then as it slows down, it helps it refine the search. Interestingly, some people do this ALONG SIDE of Adam, which to me makes some sense. Since Adam can only slow down or speed up a learning rate, but not set it to a specific value, it seems like it could be a good idea to decrease the learning rate over time as well. But, then you have another hyper parameter to tune – how much should alpha decrease each iteration?
Results
As a test, I’ve compared Adam to plain gradient descent, on finding the minimum height location of a 2D terrain made of randomized Gaussians. The code that made these diagrams and data can be found at https://github.com/Atrix256/Adam
The difference between the two methods is not as big as I’d hoped, but I think Adam does better in higher dimensions (more than the 2 parameters of x and y), where the per parameter (per dimension) adjustable training weight means Adam becomes more maneuverable. In this case for only 2 dimensions, I’m not real impressed. Here’s a relevant quote from Andrej Karpathy (https://karpathy.github.io/2019/04/25/recipe/) though:
Below are graphs to show how gradient descent (GD) with different learning rates compares to Adam with different alpha parameters, and each to themselves. Also below are animated gifs showing several points of each type searching for the minimum. Gradient descent is on the left, with Adam on the right. The colored bars in the center show the relative average heights of each particle type.
In the animated images, I used only 16 greys to show the Gaussians, to make it easier to read as a topo map, but the actual heights being optimized against are continuous.
Seed = 3684993765
Seed =181763002
Seed = 1730885985
Links
The paper that introduced Adam is surprisingly easy to read. It’s called “Adam: A Method for Stochastic Optimization” https://arxiv.org/abs/1412.6980
Here are some other resources relating to Adam, and other gradient descent methods.
This post explains how to use sliced optimal transport to make blue noise point sets. The plain, well commented C++ code that goes along with this post, which made the point sets and diagrams, is at https://github.com/Atrix256/SOTPointSets.
Let’s say we wanted a set of points randomly placed in 2D, but with even density – aka we wanted 2D blue noise points in a square.
We could start by generating 1000 white noise points, then using sliced optimal transport to make them evenly spaced:
Project all points onto a random unit vector. Calculate how much to move each point to make them be evenly spaced on that 1D line, in the same order.
Calculate the average of that process done 64 times, then move the points that amount.
Repeat this process 1000 times.
Doing that, we start with the points on the left and end up with the points on the right.
That looks like good blue noise and only took 0.56 seconds to generate, but confusingly, the points are on a disk and are denser near the edges. We’ll make the density uniform in this section, and will make the points be in a square in the next section.
The reason we are getting a disk is because we are projecting the points onto random directions and making the points be on that line, up to a maximum distance away from the origin. That is forcing them to be inside of a disk.
The reason the points are denser near the edges is because we are making the points evenly spaced on the 1D lines, but there is less room vertically on the disk near the edges. That means the points have less space to occupy at the edges, causing them to bunch up.
To account for this, we need to make there be more points in the center, and fewer at the edges. More specifically, we need each point to claim an even amount of disk area, not an even amount of line width.
Below shows a circle divided into evenly spaced slices on the x axis (top), and slices of equal area (bottom). (diagram made at https://www.geogebra.org/geometry)
We can do this in a couple of steps:
Get the formula for the height of the disk has at a given x position. We’ll call this y=f(x).
Integrate that function to get a y=g(x) function that tells you how much area the disk has at, and to the left of, a given x position.
Divide g(x) by the area of the disk to get a function h(x) that gives a value between 0 and 1.
Invert h(x) by changing it from y=h(x) to x=h(y) and solving for y. That gives us a function y=i(x).
Plug the evenly spaced values we were using on the line into the function i, to get positions on the line that give us equal values of disk area.
This might look more familiar from a statistics point of view. f(x) is an unnormalized PDF, g(x) is an unnormalized CDF, h(x) is a CDF (normalized), and i(x) is an inverse CDF.
If we have a disk centered at (0,0) with radius 0.5, the height of the disk at location x can be found using the Pythagorean theorem. We know x, and the hypotenuse is the radius 0.5, so . That only gives us the top half of the disk height, but we can double it to get the bottom half as well: .
We can take the integral of that to get . Dividing that by the area of the disk, we get
Instead of inverting that function, I made a table with 1000 (x,y) pairs of x and y=h(x). In the table, x is between -0.5 and +0.5, but I also made y be between -0.5 and +0.5. That is a bit different than a CDF where y goes from 0 to 1, but this is just a convention I’m choosing; a standard CDF would work as well. To evaluate i(x), i find the location where the given value shows up in the table as a y value, and give back the x that gave that y value, using linear interpolation.
The sliced optimal transport sampling paper also used a numerical inversion of the CDF, but they did it differently, involving gradient descent. I am happy with the results I got, so think the table is good enough.
Doing that, we get blue noise points nicely distributed on a disk, and it only took 0.73 seconds to generate.
Blue noise points on a disk are useful for a couple things. If you want to select points to sample on a disk light source (such as the sun, perhaps), you can use blue noise points to get good coverage over the sampling domain, and without aliasing problems from a low sample count. You can also take the (x,y) values of these 2D points and add a Z component with the positive value that makes it a normalized vector () and that gives you blue noise distributed cosine weighted hemispherical samples, useful for importance sampling the term in lighting calculations.
Fun fact: this projection of a density function (PDF) onto a line is actually a partial Radon transform, which is from Tomography, and relates to how we are able to make images from xrays.
Generating Points in a Square
While points in a disk are useful, we started by trying to make points in a square. To do that, we’ll need to project the area of the square down onto each 1D line, like we did the circle, and use that area to drive where the points go, to make each point have equal area of the square. This is more challenging than the circle case, because a circle is rotationally symmetric, but a square isn’t, so the function to give us evenly spaced points also has to take the projection direction into account.
We ultimately want to project the area of a rotated square onto a 1d line, and use that to give an equal amount of area to each point we distribute along that line.
My mathematically talented friend William Donnelly came up with this function (which is in squaredcdf.h in the repo). The value u should be between -0.5 and +0.5, and will return the x value of where the point should go on the line defined by direction d. The square has width of 1, and is centered at the origin, so the unrotated square corners are at (+/- 0.5, +/- 0.5). A sketch of the derivation can be found at https://blog.demofox.org/2023/12/24/deriving-the-inverse-cdf-of-a-rotated-square-projected-onto-a-line/.
inline float InverseCDF(float u, float2 d)
{
float c = std::max(std::abs(d.x), std::abs(d.y));
float s = std::min(std::abs(d.x), std::abs(d.y));
if (2 * c * std::abs(u) < c - s)
{
return c * u;
}
else
{
float t = 0.5f * (c + s) - sqrtf(2.0f * s * c * (0.5f - std::abs(u)));
return (u < 0) ? -t : t;
}
}
If we use that to space the points along each line, again doing a batch size of 64, and doing 1000 batches, we end up with this, which took 0.89 seconds to generate.
The first image is the starting state, the 2nd and 3rd are the ending states, and last is the DFT to show the characteristic looking blue noise spectrum with a dark circle in the center of suppressed low frequencies.
Regarding the third image, if you think through what’s happening with random projections, and distributing points according to those 1D projections, there is nothing explicitly keeping the points in the square. It’s surprising to me they stay in there so well. The third image shows how well they are constrained to that square.
The blue noise points don’t look that great though. There are some bent line type patterns near the edges, and especially the corners of the square. Unfortunately, these points seem to have reached a point of “overconvergence” where they start losing their randomization. You can see the same thing in the sliced optimal transport paper. One solution to this problem is to just not do as many iterations. Below is what the point set looks like after one tenth of the iterations, or 100 batches. The DFT shows stronger lower frequencies, but the visible patterns in the point sets are gone. As we’ll see further down, adjusting the batch size may also help.
Comparison To Other Blue Noise Point Sets
Below are shown three types of blue noise point sets, each having 1000 points. Points are on the left, and DFT is on the right to show frequency characteristics.
Dart Throwing – Points are placed in a square going from [0,1]. Each point is generated using white noise and is accepted if the distance to any existing point is greater than 0.022. I had to tune that distance threshold by hand to make it as large as possible, but not so large that it couldn’t be solved for 1000 points. Wrap around / toroidal distance is used: (https://blog.demofox.org/2017/10/01/calculating-the-distance-between-points-in-wrap-around-toroidal-space/)
Sliced Optimal Transport aka SOT – The method in this post, using a batch size of 64, and 100 iterations.
I’d say dart throwing is the worst, with a less pronounced dark ring in the frequency domain. For the best, I think it’s a toss up between SOT and MBC. MBC suppresses more of the low frequencies, but SOT more strongly suppresses the frequencies it does suppress. We saw earlier that doing more iterations can increase the size of the dark circle in the DFT, but that it makes the point set too ordered. It may be situational which of these you want. There is a big difference between using blue noise point sets for organic looking object placement in a game, and using blue noise point sets for sparse sampling meant to then be filtered by a matching gaussian filter. The first is an artistic choice and the second is a mathematical one.
There other methods for generating blue noise point sets. “Gaussian Blue Noise” by Ahmed et al. is the state of the art, I believe: https://arxiv.org/abs/2206.07798
Algorithm parameters
A batch size of 64 is used in generating all the point sets so far. Here’s what other batch sizes look like, at 1000 iterations.
it looks like increasing batch may also be a way to get rid of the regularity we saw before at a batch size of 64, after 1000 iterations. More investigation needed here, but being able to run an algorithm to completion, instead of having to stop it early at some unknown point, is a much better way to run an algorithm.
Here is a log/log convergence graph that shows average pixel movement each batch. The x axis is the total number of random lines projected, which is constant across all batch sizes. The number of iterations is increased for smaller batch sizes and decreased for larger batch sizes to give an “equal computation done” comparison. This isn’t necessarily a graph of quality though, it just shows each batch size reaching the final state, whatever quality that may be. More investigation is needed.
We’ve been using white noise for generating the random directions to project onto. What if we instead used the golden ratio low discrepancy sequence to generate pseudorandom rotations that were more evenly spaced over the circle? Below is a graph that shows that for both square and circle, 1000 iterations of a batch size of 64. The golden ratio variants move slightly less, and move less erratically, but the difference seems pretty small.
I can’t tell much of a difference in the disk points:
And the square points look real similar too. The DFT does seem smoother for the golden ratio points though, and the light circle around the dark center seems to be less lumpy. This might be worth looking into more later on, but it seems like a minor improvement.
I think it could also be interesting to try starting with a stratified point set, or a low discrepancy point set, instead of white noise, before running the SOT algorithm. That might help find a better minimum, or find it faster, by starting with a point set that was already pretty well evenly spaced. Perhaps more generally, this algorithm could start with any other blue noise point set generation method and refine the results, assuming that the points created were lower quality than this algorithm can provide.
Using Density Maps
Generating these sample points using sliced optimal transport involves projecting the area of where we want our points to be down onto a 1D line, and using that projected area as a PDF for where points should go to get equal amounts of that projected area.
What if instead of projecting down a shape, which projects a boolean “in shape” or “not in shape” onto the line, what if we instead used a greyscale image as a density map and projected the pixel values down onto the line to make a PDF?
I did that, projecting each pixel of a 512×512 image down onto 1000 buckets along the line, adding the pixel value into each bucket it touches, but multiplying the pixel value by how much of the bucket it takes up. The pixel values are also inverted so that bright pixels are places where dots should be less dense, because I’m doing black dots on a white background.
Below right is the image used as the density map. Below left is the starting point and below middle is the final result, using 20,000 points, 1,000 iterations and a batch size of 64. It took 50 seconds to generate, and I think it looks pretty nice! I’m amazed that doing random 1d projections of an image onto a line results in such a clean image with such sharp details.
The reason this takes so much longer to generate that the other point sets seen so far is the looping through all 512×512=262,144 pixels and projecting them down onto the 1d line for each random projection vector. I’ve heard that there are ways to make this faster by working in frequency space, but don’t know the details of that.
In game development, perhaps this could be used to place trees in a forest, where a greyscale density texture controlled the density of the trees on the map.
I haven’t ever seen it before, but you could probably use a density map with both dart throwing and Mitchell’s best candidate as well. Both of those algorithms calculate distance between points. You could get the density value for each point, convert that density to a radius value, and subtract the two radius values from the distance between the points. It would be interesting to compare their quality vs these results.
Generating Multiclass Samples
Multiclass blue noise point sets are blue noise point sets where each point belongs to a class. Each class of points should independently be blue noise, but combined combinations of classes should also be blue noise.
The “Sliced Optimal Transport Sampling” paper that this post is exploring has support for this. Every other iteration, after projecting the points onto the line and sorting them, they then make sure the classes are interleaved on that line. They show it for 2 and 3 classes with equal point counts.
If using this for placing objects on a map, you might have one class for trees, another class for bushes, and another class for boulders. However, you might not want an equal number of these objects, or equivalently, you may want more space between trees, than you want space between bushes. To do that, you’d need to have different point counts per class.
Luckily that ends up being pretty easy to do. Let’s say we have three classes with weights 1, 4 and 16. Those sum up to 21. When generating your random points, you can use the point index to calculate which class it is:
int class = 2;
if (index % 21 == 0)
class = 0;
else if (index % 21 < 5)
class = 1;
Then, when doing the “interleave” step that is done every other iteration, after sorting the points, we make sure that there is a class 0, then four class 1s, then sixteen class 2s, repeating that pattern over and over.
Doing that in a square, with 1000 points, 1000 iterations, and a batch size of 64 gives us this, which took 0.57 seconds to generate. Click the image to see it full sized in another window. The top row shows the classes as colors RGB. The middle shows the combined point set, and bottom shows the DFT.
This works just fine with density maps too, generating this in 49 seconds:
Optimal Transport
Generating multi class and uneven density blue noise point sets works well using sliced optimal transport, but where is the actual optimal transport of this, and what exactly is being transported?
Going back to the last post’s explanation of optimal transport being about finding which bakeries should deliver to which cafe’s to minimize the total distance traveled, the initial point sets are the bakeries, and instead of there being discrete cafe’s, the PDF / density maps define a fuzzy “general desire for baked goods in an area”.
The optimal transport here is finding where to deliver baked goods to most evenly match demand for baked goods, and doing so, such that the sum of the distance traveled by all baked goods is minimized.
When we run the sliced optimal transport algorithm, we do eventually end up with the points being in the “optimal transport” position at the end, but the path that the points traveled to get there are not the optimal transport solution.
The optimal transport solution is the line from the initial point locations to their ending locations. That optimal transport solution will be a straight line, while the path the points took in SEARCH of the optimal transport solution may not be straight lines.
Below is the evolution of 20,000 points as they search for the optimal transport solution, over 100 iterations, with a batch size of 64. The points are not taking the path of optimal transport, the points are searching for their final position for the optimal transport solution.
Below is an animation that interpolates from the starting point set, to the ending point set, over the same amount of time (100 frames). This is the points taking the path of the optimal transport solution!
Closing
The paper which introduced dart throwing was “Stochastic sampling in computer graphics”, Cook 1986. That paper explains how in human eyes and some animals, the fovea (the place where vision is the sharpest) places photoreceptors in a hexagonal packing. This is interesting because a hexagonal grid is more evenly distributed than a regular grid, which has diagonal distances 50% longer than cardinal distances. Outside of the fovea, where the photoreceptors are sparser, a blue noise pattern is used, which is known to be roughly evenly spaced, but randomized which avoids aliasing.
To me, this is nature giving us a hint for how to do sampling or rendering. For best quality and higher sample counts, a hex grid (or more generally, a uniform grid) is best, and for lower / sparser sample counts, blue noise is best.
I hope you enjoyed this. Next up will be a post looking more deeply at an interesting 2022 paper “Scalable multi-class sampling via filtered sliced optimal transport” (https://www.iliyan.com/publications/ScalableMultiClassSampling). That paper is what convinced me I needed to learn optimal transport, and is what led to this series of blog posts.
. The value 
is what number is on the dice face, and the value 
says how many faces have that value.
so that 
is 0.
because it makes all the 
power terms go away and leaves only the value of 
.
because that will make all 
. The only way that can be true is if each dice has the 
. If we plug 1 in for 
, the 
which simplifies to:
.
can be true, when ![A,B,C \in [0,6]](https://s0.wp.com/latex.php?latex=A%2CB%2CC+%5Cin+%5B0%2C6%5D&bg=000000&fg=666666&s=0&c=20201002)
.
. If we want to do that three times – once for A, once for B, once for C – and add them together, we multiply that generating function by itself three times, aka we cube it.
is 28, which means there are 28 different ways to add the three numbers between 0 and 6 together, to get 12.
into infinity, but we know there are 20 jellybeans total, so we can stop at 20.
. That only gives us the top half of the disk height, but we can double it to get the bottom half as well: 
.
. Dividing that by the area of the disk, we get 
) and that gives you blue noise distributed cosine weighted hemispherical samples, useful for importance sampling the 
term in lighting calculations.
