Building a quadtree filter in GLSL, using a probabilistic approach

In this article, we will try to build a quadtree filter running in GLSL. This article does not intend to provide a way to build a quadtree of exploitable data, but rather an artistic filter that can be applied over an image, grouping together sections where the colors are similar.

1. How do we decide if a quad of our tree will be divided into 4 smaller quads

First step is to find a way to determine when we will divide a quad of the image into smaller quad and when we don’t. Let’s consider a really basic example:

Figure 1. An example image divided into 4 quads

In our study case, we have 4 quads, the first 4 quads of our tree. With a graphical analysis, we know by instinct that quad 1 & 2 should be divided into 4 smaller quads because the level of details (cf: the variation of the colors) is too important, while quad 4 can be left as it is since we don’t have any detail in it.

Quad 3 is a bit trickier. There are some color variations in it, however not enough for it to be obvious. Because we are working with GLSL, we will need to sample the texture within each quad n times, therefore we will need to pick random points in each quad to determine the color variations based on these samples. The next figure demonstrates the result of the sampling on our example:

Figure 2. Sampling 9 points in each of the 4 quads

If for instance, we take 9 samples at the same position of each of the quads, we immediatly notice the issue if details are located in a small area of a quad: details are missing from the samples we picked. This will become a limitation of our algorithm. We will trade computational power with precision, and rely on the samples we take from the quad. If we land on enough of the pixels where there is detail, we will be able to extrapolate whether or not the active quad needs to be divided again, otherwise we will consider that the quad does not have enough detail to be splitted again and we will leave it. Thus the probabilistic approach.

2. The color variation of a serie of samples: variance

In probability theory and statistics, variance is the expectation of the squared deviation of a random variable from its mean. Informally, it measures how far a set of (random) numbers are spread out from their average value. That’s what we are looking for in our case: for how much does our samples spread out from their average. Using the variance allows us to take a small number of samples within a quad and find, with a great certitude, by how much the color of the samples are appart from each other.

In our case, and for a generic purpose, we need to compute the variance for the red, green and blue channels and then average the 3 results. If we compute the variance over the grayscale, it will be faster but we could miss some informations. For instance, if we have 3 samples ([1, 0, 0], [0, 1, 0], [0, 0, 1]), the variance on the grayscale will be 0.

The following formulas will describe how to compute the variance over 1 color component, which can be extended to the 3 components. If we have $n$ samples, we can compute the average $\overline{x}$ likeso:

$\overline{x} = \frac{1}{n} \cdot \sum_{i=0}^{n} x_i$

We can then compute the variance $V$ using the following formula:

$V = (\frac{1}{n} \cdot \sum_{i=0}^{n} {x_i}^2) - {\overline{x}}^2$

Once the variance is computed, we can compare it to a defined variance threshold, and if the variance of a particular quad is above, we need to divide this quad into 4 smaller ones, and repeat the operation.

Great, we now have a way to determine if a quad needs to be divided. Let’s port this on GLSL.

3. Thinking our shader

For those familiar with it, when working with GLSL you need to think differently. Because our shader will be executed on each texel in parallel, there are a few criteria required for our algorithm to work:

we need to know in which quad the texel we’re working with belongs
for every texel within a quad, we need to find the same variance over the quad

Until the variance of the active quad of the texel is greater than the threshold, we’ll keep dividing the active square into a smaller one. Let’s first meet our first criteria. The following illustration demonstrates how we proceed to subdivide the quads in which our texel belongs, iteration by iteration:

Figure 3. Iterating into smaller and smaller quads

To compute the coordinates of the quads, we can use a simple trick. Let $N_b$ be the active number of divisions. If $N_b = 1$ , we don’t have any subdivisions, if $N_b = 2$ , we have 2 vertical and 2 horizontal subdivisions thus 4 squares, if $N_b = 4$ , we have 8 squares… so on, so forth. Each time we multiply $N_b$ by 2, we divide the quads into 4 smaller quads. Given the value of $N_b$ and the coordinates $P$ of the texel, we can compute the coordinates $P$ of the center of the quad in which the texel belongs using the following formula (reminder: P coordinates belongs to a unit square):

$C = \frac{floor(P * {N_b}) + 0.5}{N_b}$

And the length $L$ of a side of the quad:

$L = \frac{1}{N_b}$

Knowing the center and the size of a side is enough.

Let’s now meet the 2nd criteria: find a way to pick the same samples

To find “random” points into a quad, that are the same for every texel of the quad, we need to find the common property they share. We have 2: the center of the quad and the size of a side. We also need a deterministic function which, given the same input, will output the same value, and small variations in the input results in chaotic variations in the output: hash functions. Hash functions are a very common technique used in GLSL to create noise. I will let you dive into the topic if you’re interested. Let me show you our hash function:

vec2 hash22(vec2 p) { 
    float n = sin(dot(p, vec2(41, 289)));
    return fract(vec2(262144, 32768)*n);    
}

It takes any vec2 as an input and returns a vec2 as an output. To pick $n$ samples of a quad, we feed to that hash function the center of the quad + $({i_n}, 0)$ where $i_n$ is the index of the sample. Our hash function outputs values between 0 and 1. We’ll need to subtract 0.5 to both coordinates to pick points arround the center of the quad. Also, we need to downscale the output vector so that it fits into the quad we want to sample. Which, translated into GLSL code, gives:

for (int i = 0; i < NB_SAMPLES; i++) {
    float fi = float(i); // cast index to float for computations
    // pick a random 2d point using the center of the active quad as input
    // this ensures that for every point belonging to the active quad, 
    // we pick the same samples
    vec2 rnd = hash22(center.xy + vec2(fi, 0.0)) - 0.5;
    vec2 sampleCoord = quadCenter + rnd * quadSize;
}

4. Writing our shader

This part is pretty straightforward since everything is defined. I will not go into the details but rather provide the whole code on shadertoy. In this example, the mouse X drives the variance threshold: