@MaximeHeckel

On Crafting Painterly Shaders

October 29, 2024 / 35 min read

Last Updated: October 29, 2024

Writing a shader that can reproduce the look and feel of aquarelle, watercolor, or gouache to obtain a more painterly output for my WebGL scenes has always been a long-term goal of mine. Inspired by the work of very talented 3D artists such as @simonxxoo or @arpeegee, the contrast between paintings and the added dimension allowed by 3D renderers was always very appealing to me. On top of that, my recent work with stylized shaders to mimic the hand-drawn Moebius art style emphasized not only that obtaining such stylized output was possible but also that post-processing was more likely than not the key to emulating any artistic style.

After several months of on/off research that frankly was not leading to much, I stumbled upon a smoothing filter I never heard about before: the Kuwahara filter. Lucky for me, it turned out that many (very smart) people published papers about it and its ability to transform any image input into a painting-like work of art. By implementing it into a custom post-processing pass and going through many ups and downs in the process, I got very close to my original objective by building what is essentially a real-time 3D painting running in the browser:

This article is the culmination of months of work, trial and error, and research to craft the perfect painterly shader for your next WebGL project. In it, we will deep dive into the characteristics of the Kuwahara filter that make it so great at transforming our work into paintings, as well as look into several techniques highlighted in the many papers I read to improve the original implementation and have it return a sharper, and more anisotropic output. Finally, we'll go through a bit of color correction and the use of textures to achieve a satisfying and accurate painting effect.

Painterly post-processing with the Kuwahara filter

Up until then, my only attempts at painterly shader consisted of several implementations of custom materials focused on emulating specific aspects of paint that I instinctively listed as essential to get a convincing output:

  • ArrowAn icon representing an arrow
    The absence of texture/normal detail
  • ArrowAn icon representing an arrow
    Preserved edges
  • ArrowAn icon representing an arrow
    Quantized colors

Below is an example of one of my "best" iterations from that exploratory work. It notably uses a Voronoi noise to smudge normals and a paintbrush texture for a more realistic look and feel.

However, despite my best efforts, this work was dead in the water as my implementation was severely flawed and only worked with a subset of meshes. From then on, I knew I had to go the post-processing route to solve most of those issues, but I didn't know how. Then, one day, I stumbled upon the beautiful work of @arpeegee, who gave me the keywords I was missing all this time: Kuwahara filter.

I'm done with this one! Learned a lot in the process 😌 Blender / Cycles https://t.co/VCkYOJeBNw

I'm done with this one! Learned a lot in the process 😌
Blender / Cycles https://t.co/VCkYOJeBNw
I'm done with this one! Learned a lot in the process 😌
Blender / Cycles https://t.co/VCkYOJeBNw

Characteristics of the Kuwahara filter

Originally intended to remove noise through smoothing from medical images, the Kuwahara filter is most commonly used today to transform any input into stylistic paintings through post-processing.

Unlike other smoothing filters that reduce noise but also, on the way, blur our edges, the Kuwahara filter is capable of preserving edges and corners. That specificity is what makes it so good at achieving painting-like artistic effects by getting us two crucial visual properties of paintings:

  • ArrowAn icon representing an arrow
    The smoothing erases some of the texture details.
  • ArrowAn icon representing an arrow
    The filter preserves the edges and, in some cases, increases the sharpness compared to the original input.

It achieves this by executing the following steps on each pixel of our input, which in our case would be our underlying 3D scene:

  1. ArrowAn icon representing an arrow
    We center a box around a pixel.
  2. ArrowAn icon representing an arrow
    We divide the box into 4 "sub-boxes" or Sectors.
  3. ArrowAn icon representing an arrow
    We calculate the average color and variance for each Sector.
  4. ArrowAn icon representing an arrow
    We set our pixel to the average color of the Sector with the lowest variance.

I built the widget below to visually represent those sectors surrounding a given pixel at work:

You can see by selecting the edge examples that:

  • ArrowAn icon representing an arrow
    when a pixel sits just outside an edge, the Sector with the lowest variance will always be outside the edge
  • ArrowAn icon representing an arrow
    when a pixel sits just inside an edge, the Sector with the lowest variance will be inside the edge

highlighting the edge-preserving capabilities of this filter.

The widget below showcases how this process of picking the average color of the Sector with the lowest variance impacts an entire image:

Notice how, with a reasonable kernel size, we maintain the overall shapes featured in the underlying image but erase many details like, in this case, transparency. However, we can also observe that with a higher kernel size, the image filter, unfortunately, falls apart and denatures quite drastically our input, indicating that we will have to strike the right balance between kernel size and the strength of our painterly post-processing effect.

Our first implementation of the Kuwahara filter as a post-processing effect

Now that we have a good understanding of the inner workings of the Kuwahara filter, let's implement it as a post-processing effect on top of a React Three Fiber scene. Much like the work done in Moebius-style post-processing and other stylized shaders, we will define our custom post-processing pass using the Pass class from the post-processing package.

This post-processing pass will:

  • ArrowAn icon representing an arrow
    Take our underlying scene as an input.
  • ArrowAn icon representing an arrow
    Implement the Kuwahara filter in its fragment shader.
  • ArrowAn icon representing an arrow
    Output the resulting scene.

The implementation of it in GLSL goes as follows:

Implementation of the Kuwahara filter in GLSL

1
#define SECTOR_COUNT 4
2
3
uniform int kernelSize;
4
uniform sampler2D inputBuffer;
5
uniform vec4 resolution;
6
uniform sampler2D originalTexture;
7
8
varying vec2 vUv;
9
10
void getSectorVarianceAndAverageColor(vec2 offset, int boxSize, out vec3 avgColor, out float variance) {
11
vec3 colorSum = vec3(0.0);
12
vec3 squaredColorSum = vec3(0.0);
13
float sampleCount = 0.0;
14
15
16
for (int y = 0; y < boxSize; y++) {
17
for (int x = 0; x < boxSize; x++) {
18
vec2 sampleOffset = offset + vec2(float(x), float(y));
19
vec3 color = sampleColor(sampleOffset);
20
colorSum += color;
21
squaredColorSum += color * color;
22
sampleCount += 1.0;
23
}
24
}
25
26
// Calculate average color and variance
27
avgColor = colorSum / sampleCount;
28
vec3 varianceRes = (squaredColorSum / sampleCount) - (avgColor * avgColor);
29
variance = dot(varianceRes, vec3(0.299, 0.587, 0.114));
30
}
31
32
void main() {
33
vec3 boxAvgColors[SECTOR_COUNT];
34
float boxVariances[SECTOR_COUNT];
35
36
getSectorVarianceAndAverageColor(vec2(-kernelSize, -kernelSize), kernelSize, boxAvgColors[0], boxVariances[0]);
37
getSectorVarianceAndAverageColor(vec2(0, -kernelSize), kernelSize, boxAvgColors[1], boxVariances[1]);
38
getSectorVarianceAndAverageColor(vec2(-kernelSize, 0), kernelSize, boxAvgColors[2], boxVariances[2]);
39
getSectorVarianceAndAverageColor(vec2(0, 0), kernelSize, boxAvgColors[3], boxVariances[3]);
40
41
float minVariance = boxVariances[0];
42
vec3 finalColor = boxAvgColors[0];
43
44
for (int i = 1; i < SECTOR_COUNT; i++) {
45
if (boxVariances[i] < minVariance) {
46
minVariance = boxVariances[i];
47
finalColor = boxAvgColors[i];
48
}
49
}
50
51
gl_FragColor = vec4(finalColor, 1.0);
52
}

In the code featured above, we can see that:

  1. ArrowAn icon representing an arrow
    We're allocating an array of vec3 and float to store the average color and variance for each of our Sectors, respectively.
  2. ArrowAn icon representing an arrow
    We're computing those stats by sampling the color for every pixel in a given Sector along the x-axis and y-axis and calculating the average and variance using the formulas we saw in the previous part.
  3. ArrowAn icon representing an arrow
    We then set minVariance and finalColor as the variance and average color of the first Sector by default.
  4. ArrowAn icon representing an arrow
    We then loop through the other Sectors to find the Sector with the lowest variance and set its corresponding average color as the output of our fragment shader.

Doing this will result in a version of our scene with fewer texture details while preserving the overall shape of our geometries and also introducing artifacts that look somewhat like brush strokes:

The Papari extension

Through our initial implementation of the Kuwahara filter as a custom post-processing pass, we already achieved a somewhat convincing paint effect. However, it suffers from several drawbacks at high kernel sizes that may make our underlying scene not discernable by the viewer, and the appearance of our "brush strokes" could most likely be improved.

Ideally, we could fix both these issues in one go by allowing larger kernel sizes for a more intense stylized effect for our scene and also have more natural-looking artifacts with a slight tweak of our Kuwahara filter. That is what Giuseppe Papari proposes in his paper Artistic Edge and Corner Enhancing Smoothing. In it, he presents an extension to the Kuwahara filter that:

  1. ArrowAn icon representing an arrow
    Removes the square-shaped kernel in favor of a circular one.
  2. ArrowAn icon representing an arrow
    Improves the weight influence that each pixel has on the filter.

In this part, we will implement some of those improvements with an additional section on potential performance improvements. It may sound absurd but; I like my paintings to run as close to 60fps as possible.

Circular kernel

A circular-shaped kernel allows us to increase the number of Sectors when evaluating the color of a given pixel. Its box-shaped counterpart can only allow for up to four Sectors making more complex edge lines not well-preserved. Through his experiments, Papari found that eight Sectors were the ideal amount to balance output quality and good performance.

Diagram showcasing both box-shaped and circular kernels and highlighting which parts of the image are being smoothed while maintaining the sharpness of the edges.

The widget below showcases the impact that this subtle modification in the implementation of our filter has on a simple image:

Here, we can see that, compared to what we observed in the first part, this version of the Kuwahara filter yields a way better output even at a larger kernel size! The circular shape and the higher number of Sectors allow us to preserve edge lines for more complex angles and patterns.

We can tweak the fragment shader code to implement Papari's circular kernel into our post-processing pass by modifying the parts that compute variance and average color by Sector:

Switching to a circular kernel in our Kuwahara filter implementation

1
//...
2
3
void getSectorVarianceAndAverageColor(float angle, float radius, out vec3 avgColor, out float variance) {
4
vec3 colorSum = vec3(0.0);
5
vec3 squaredColorSum = vec3(0.0);
6
float sampleCount = 0.0;
7
8
for (float r = 1.0; r <= radius; r += 1.0) {
9
for (float a = -0.392699; a <= 0.392699; a += 0.196349) {
10
vec2 sampleOffset = r * vec2(cos(angle + a), sin(angle + a));
11
vec3 color = sampleColor(sampleOffset);
12
colorSum += color;
13
squaredColorSum += color * color;
14
sampleCount += 1.0;
15
}
16
}
17
18
// Calculate average color and variance
19
avgColor = colorSum / sampleCount;
20
vec3 varianceRes = (squaredColorSum / sampleCount) - (avgColor * avgColor);
21
variance = dot(varianceRes, vec3(0.299, 0.587, 0.114));
22
}
23
24
25
void main() {
26
vec3 sectorAvgColors[SECTOR_COUNT];
27
float sectorVariances[SECTOR_COUNT];
28
29
for (int i = 0; i < SECTOR_COUNT; i++) {
30
float angle = float(i) * 6.28318 / float(SECTOR_COUNT);
31
getSectorVarianceAndAverageColor(angle, float(radius), sectorAvgColors[i], sectorVariances[i]);
32
}
33
34
//...
35
}

In the end, the filter's principle is roughly the same. What changes here is which pixel is picked for a given sector and included in the average and variance.

Diagram showcasing how the sampling of a slice of the circular kernel is split.

The demo below uses this extended version of the Kuwahara filter on top of our scene. It lets us use a larger kernel size, yielding a more stylized output while keeping our scene discernable:

A better weighting of colors

The artifacts left by the filter are still quite visible at large kernel sizes. According to Papari's study, this effect is due to the weighting of the colors when calculating the average and the variance of a given Sector, as it does not discriminate any pixel. If a pixel is in a sector, it contributes the same amount to the calculations.

Using Gaussian weights fixes this issue:

  • ArrowAn icon representing an arrow
    It provides a gradual fall-off from the center of the Sector to its edges.
  • ArrowAn icon representing an arrow
    This results in more natural-looking transitions from one Sector to another, thus avoiding abrupt changes that would otherwise make those artifacts very visible.
  • ArrowAn icon representing an arrow
    It emphasizes the central pixels more likely to represent the most accurate sector color, giving them more importance when calculating the average color.
Diagram comparing standard weighting vs Gaussian weighting. The darker the color of a pixel is, the 'heavier' the weight for said pixel is.

To use Gaussian weighting in our filter, we only need to change a few lines from the previous implementation:

Using Gaussian weighting in our Kuwahara filter implementation

1
// ...
2
3
float gaussianWeight(float distance, float sigma) {
4
return exp(-(distance * distance) / (2.0 * sigma * sigma));
5
}
6
7
void getSectorVarianceAndAverageColor(float angle, float radius, out vec3 avgColor, out float variance) {
8
vec3 weightedColorSum = vec3(0.0);
9
vec3 weightedSquaredColorSum = vec3(0.0);
10
float totalWeight = 0.0;
11
12
float sigma = radius / 3.0;
13
14
for (float r = 1.0; r <= radius; r += 1.0) {
15
for (float a = -0.392699; a <= 0.392699; a += 0.196349) {
16
vec2 sampleOffset = r * vec2(cos(angle + a), sin(angle + a));
17
vec3 color = sampleColor(sampleOffset);
18
float weight = gaussianWeight(length(sampleOffset), sigma);
19
20
weightedColorSum += color * weight;
21
weightedSquaredColorSum += color * color * weight;
22
totalWeight += weight;
23
}
24
}
25
26
// Calculate average color and variance
27
avgColor = weightedColorSum / totalWeight;
28
vec3 varianceRes = (weightedSquaredColorSum / totalWeight) - (avgColor * avgColor);
29
variance = dot(varianceRes, vec3(0.299, 0.587, 0.114)); // Convert to luminance
30
}
31
32
// ...

As we can see in the code snippet above:

  1. ArrowAn icon representing an arrow
    We compute the weight using the Gaussian formula.
  2. ArrowAn icon representing an arrow
    We then take the resulting weight into account when calculating our weighted color sum and weighted squared color sum.
  3. ArrowAn icon representing an arrow
    Those weights are then reflected in the final result of the variance and average color for each Sector.

Now, you may wonder whether adding all those nitty-gritty improvements to a filter that already performed quite well is worth it. You can see the result for yourself in the comparison below between a frame of our first scene using the original implementation of the Kuwahara filter and one using the Papari extension:

Before
After
ArrowAn icon representing an arrowArrowAn icon representing an arrow
Comparing the output of our extended Kuwahara filter without and with Gaussian weighting.

The result looks already better than previously, even at a low kernel size! Let's observe now the difference between both implementations for a higher value:

Before
After
ArrowAn icon representing an arrowArrowAn icon representing an arrow
Comparing the output of our extended Kuwahara filter without and with Gaussian weighting at a high kernel size.

Fixing a performance pitfall

The Gaussian function used to calculate the weights of our pixel is, unfortunately, quite expensive to run, even more so in all those nested loops necessary to sample each pixel of a given Sector. This improvement is costing us a lot of frames, and it would be very tempting to disregard it and revert to the original weighting method.

However, a team of smart individuals found a way to squeeze more performance out of the filter while maintaining the advantages given by the Gaussian. In Anisotropic Kuwahara Filtering with PolynomialWeighting Functions, they propose to replace the Gaussian with a polynomial that roughly approximates it.

[(x + ζ) − ηy^2]^2

Using Polynomial weighting in our Kuwahara filter implementation

1
// ...
2
3
float polynomialWeight(float x, float y, float eta, float lambda) {
4
float polyValue = (x + eta) - lambda * (y * y);
5
return max(0.0, polyValue * polyValue);
6
}
7
8
void getSectorVarianceAndAverageColor(float angle, float radius, out vec3 avgColor, out float variance) {
9
vec3 colorSum = vec3(0.0);
10
vec3 squaredColorSum = vec3(0.0);
11
float weightSum = 0.0;
12
13
float eta = 0.1;
14
float lambda = 0.5;
15
16
for (float r = 1.0; r <= radius; r += 1.0) {
17
for (float a = -0.392699; a <= 0.392699; a += 0.196349) {
18
19
vec2 sampleOffset = vec2(r * cos(angle + a), r * sin(angle + a));
20
vec3 color = sampleColor(sampleOffset);
21
float weight = polynomialWeight(sampleOffset.x, sampleOffset.y, eta, lambda);
22
//...
23
}
24
}
25
// ...
26
}
27
28
// ...

When trying out this formula, I noticed the performance improvements highlighted in the paper but also that I could get a satisfying output at a smaller kernel size compared to what we've implemented so far. I'm still not exactly sure why 🤷‍♂️, but I'll take it.

From Structure Tensor to Anisotropic Kuwahara filter

Papari's extension of the Kuwahara filter is already a step towards a lovely painterly shader, however, (I must sound like a broken record at this point) we still can perceive some artifacts left by the filter on the resulting scene. Indeed, the painting effect is applied indiscriminately to the input without any consideration for the overall structure of the current frame. In the case of a true painting, the artist would adapt the flow of the brush based on the direction of the edge. We're not getting this here.

Luckily, Anisotropic Kuwahara Filtering with PolynomialWeighting Functions answers this problem by building upon the concepts of the extended Kuwahara filter we just saw.

In this paper, the authors suggest that we could adapt our circular kernel to the local features of the input during sampling by squeezing and rotating it, thus letting us sample our Sectors more accurately. As a result, our kernel would look more like an ellipsis rather than a circle.

Diagram comparing isotropic kernels vs anisotropic kernels.

Implementing this improvement will require us to follow several steps and make our post-processing effect multi-pass:

  1. ArrowAn icon representing an arrow
    The first pass will take the underlying scene as input and return the structure of the scene.
  2. ArrowAn icon representing an arrow
    Then based on this structure as input we'll apply our Kuwahara filter to the underlying scene.
  3. ArrowAn icon representing an arrow
    Finally, we'll add one extra pass to apply some tone mapping and other details that will make our painterly shader shine.
Diagram showcasing the multi-pass post-processing pipeline details in this section.

Sobel operator and partial derivatives

To obtain the structure of our scene, we can leverage a Sobel operator which I introduced in Moebius-style post-processing and other stylized shaders earlier this year.

Sx = [[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]]

Sy = [[-1, -2, -1], [0, 0, 0], [1, 2, 1]]

However, this time we won't use it for edge detection like we did back then. We will apply our Sobel Matrix for every pixel but instead of a rapid change in intensity, which is characteristic of edges, we want to extract the partial derivatives along the x/y axis representing the rate of change.

Applying our Sobel Matrices when sampling

1
varying vec2 vUv;
2
uniform sampler2D inputBuffer;
3
uniform vec4 resolution;
4
5
// Sobel kernels
6
const mat3 Gx = mat3( -1, -2, -1, 0, 0, 0, 1, 2, 1 ); // x direction kernel
7
const mat3 Gy = mat3( -1, 0, 1, -2, 0, 2, -1, 0, 1 ); // y direction kernel
8
9
vec4 computeStructureTensor(sampler2D inputTexture, vec2 uv) {
10
vec3 tx0y0 = texture2D(inputTexture, uv + vec2(-1, -1) / resolution.xy).rgb;
11
vec3 tx0y1 = texture2D(inputTexture, uv + vec2(-1, 0) / resolution.xy).rgb;
12
vec3 tx0y2 = texture2D(inputTexture, uv + vec2(-1, 1) / resolution.xy).rgb;
13
vec3 tx1y0 = texture2D(inputTexture, uv + vec2( 0, -1) / resolution.xy).rgb;
14
vec3 tx1y1 = texture2D(inputTexture, uv + vec2( 0, 0) / resolution.xy).rgb;
15
vec3 tx1y2 = texture2D(inputTexture, uv + vec2( 0, 1) / resolution.xy).rgb;
16
vec3 tx2y0 = texture2D(inputTexture, uv +