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:

- The absence of texture/normal detail
- Preserved edges
- 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

### 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:

- The smoothing erases some of the texture details.
- 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:

- We center a box around a pixel.
- We divide the box into 4 "sub-boxes" or
**Sectors**. **We calculate the average color and variance for each Sector**.- 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:

- when a pixel sits just outside an edge, the Sector with the lowest variance will always be outside the edge
- 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:

- Take our underlying scene as an input.
- Implement the Kuwahara filter in its fragment shader.
- Output the resulting scene.

The implementation of it in GLSL goes as follows:

Implementation of the Kuwahara filter in GLSL

1#define SECTOR_COUNT 423uniform int kernelSize;4uniform sampler2D inputBuffer;5uniform vec4 resolution;6uniform sampler2D originalTexture;78varying vec2 vUv;910void getSectorVarianceAndAverageColor(vec2 offset, int boxSize, out vec3 avgColor, out float variance) {11vec3 colorSum = vec3(0.0);12vec3 squaredColorSum = vec3(0.0);13float sampleCount = 0.0;141516for (int y = 0; y < boxSize; y++) {17for (int x = 0; x < boxSize; x++) {18vec2 sampleOffset = offset + vec2(float(x), float(y));19vec3 color = sampleColor(sampleOffset);20colorSum += color;21squaredColorSum += color * color;22sampleCount += 1.0;23}24}2526// Calculate average color and variance27avgColor = colorSum / sampleCount;28vec3 varianceRes = (squaredColorSum / sampleCount) - (avgColor * avgColor);29variance = dot(varianceRes, vec3(0.299, 0.587, 0.114));30}3132void main() {33vec3 boxAvgColors[SECTOR_COUNT];34float boxVariances[SECTOR_COUNT];3536getSectorVarianceAndAverageColor(vec2(-kernelSize, -kernelSize), kernelSize, boxAvgColors[0], boxVariances[0]);37getSectorVarianceAndAverageColor(vec2(0, -kernelSize), kernelSize, boxAvgColors[1], boxVariances[1]);38getSectorVarianceAndAverageColor(vec2(-kernelSize, 0), kernelSize, boxAvgColors[2], boxVariances[2]);39getSectorVarianceAndAverageColor(vec2(0, 0), kernelSize, boxAvgColors[3], boxVariances[3]);4041float minVariance = boxVariances[0];42vec3 finalColor = boxAvgColors[0];4344for (int i = 1; i < SECTOR_COUNT; i++) {45if (boxVariances[i] < minVariance) {46minVariance = boxVariances[i];47finalColor = boxAvgColors[i];48}49}5051gl_FragColor = vec4(finalColor, 1.0);52}

In the code featured above, we can see that:

- We're allocating an array of
`vec3`

and`float`

to store the average color and variance for each of our Sectors, respectively. - 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.
- We then set
`minVariance`

and`finalColor`

as the variance and average color of the first Sector by default. - 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:

- Removes the square-shaped kernel in favor of a circular one.
- 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.

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//...23void getSectorVarianceAndAverageColor(float angle, float radius, out vec3 avgColor, out float variance) {4vec3 colorSum = vec3(0.0);5vec3 squaredColorSum = vec3(0.0);6float sampleCount = 0.0;78for (float r = 1.0; r <= radius; r += 1.0) {9for (float a = -0.392699; a <= 0.392699; a += 0.196349) {10vec2 sampleOffset = r * vec2(cos(angle + a), sin(angle + a));11vec3 color = sampleColor(sampleOffset);12colorSum += color;13squaredColorSum += color * color;14sampleCount += 1.0;15}16}1718// Calculate average color and variance19avgColor = colorSum / sampleCount;20vec3 varianceRes = (squaredColorSum / sampleCount) - (avgColor * avgColor);21variance = dot(varianceRes, vec3(0.299, 0.587, 0.114));22}232425void main() {26vec3 sectorAvgColors[SECTOR_COUNT];27float sectorVariances[SECTOR_COUNT];2829for (int i = 0; i < SECTOR_COUNT; i++) {30float angle = float(i) * 6.28318 / float(SECTOR_COUNT);31getSectorVarianceAndAverageColor(angle, float(radius), sectorAvgColors[i], sectorVariances[i]);32}3334//...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.

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:

- It provides a gradual fall-off from the center of the Sector to its edges.
- This results in more natural-looking transitions from one Sector to another, thus avoiding abrupt changes that would otherwise make those artifacts very visible.
- It
*emphasizes*the central pixels more likely to represent the**most accurate**sector color, giving them more importance when calculating the average color.

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// ...23float gaussianWeight(float distance, float sigma) {4return exp(-(distance * distance) / (2.0 * sigma * sigma));5}67void getSectorVarianceAndAverageColor(float angle, float radius, out vec3 avgColor, out float variance) {8vec3 weightedColorSum = vec3(0.0);9vec3 weightedSquaredColorSum = vec3(0.0);10float totalWeight = 0.0;1112float sigma = radius / 3.0;1314for (float r = 1.0; r <= radius; r += 1.0) {15for (float a = -0.392699; a <= 0.392699; a += 0.196349) {16vec2 sampleOffset = r * vec2(cos(angle + a), sin(angle + a));17vec3 color = sampleColor(sampleOffset);18float weight = gaussianWeight(length(sampleOffset), sigma);1920weightedColorSum += color * weight;21weightedSquaredColorSum += color * color * weight;22totalWeight += weight;23}24}2526// Calculate average color and variance27avgColor = weightedColorSum / totalWeight;28vec3 varianceRes = (weightedSquaredColorSum / totalWeight) - (avgColor * avgColor);29variance = dot(varianceRes, vec3(0.299, 0.587, 0.114)); // Convert to luminance30}3132// ...

As we can see in the code snippet above:

- We compute the weight using the Gaussian formula.
- We then take the resulting weight into account when calculating our
**weighted color sum**and**weighted squared color sum**. - 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:

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:

### 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// ...23float polynomialWeight(float x, float y, float eta, float lambda) {4float polyValue = (x + eta) - lambda * (y * y);5return max(0.0, polyValue * polyValue);6}78void getSectorVarianceAndAverageColor(float angle, float radius, out vec3 avgColor, out float variance) {9vec3 colorSum = vec3(0.0);10vec3 squaredColorSum = vec3(0.0);11float weightSum = 0.0;1213float eta = 0.1;14float lambda = 0.5;1516for (float r = 1.0; r <= radius; r += 1.0) {17for (float a = -0.392699; a <= 0.392699; a += 0.196349) {1819vec2 sampleOffset = vec2(r * cos(angle + a), r * sin(angle + a));20vec3 color = sampleColor(sampleOffset);21float weight = polynomialWeight(sampleOffset.x, sampleOffset.y, eta, lambda);22//...23}24}25// ...26}2728// ...

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.

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

- The first pass will take the underlying scene as input and return the
**structure of the scene**. - Then based on this structure as input we'll apply our Kuwahara filter to the underlying scene.
- Finally, we'll add one extra pass to apply some tone mapping and other details that will make our painterly shader shine.

### 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

1varying vec2 vUv;2uniform sampler2D inputBuffer;3uniform vec4 resolution;45// Sobel kernels6const mat3 Gx = mat3( -1, -2, -1, 0, 0, 0, 1, 2, 1 ); // x direction kernel7const mat3 Gy = mat3( -1, 0, 1, -2, 0, 2, -1, 0, 1 ); // y direction kernel89vec4 computeStructureTensor(sampler2D inputTexture, vec2 uv) {10vec3 tx0y0 = texture2D(inputTexture, uv + vec2(-1, -1) / resolution.xy).rgb;11vec3 tx0y1 = texture2D(inputTexture, uv + vec2(-1, 0) / resolution.xy).rgb;12vec3 tx0y2 = texture2D(inputTexture, uv + vec2(-1, 1) / resolution.xy).rgb;13vec3 tx1y0 = texture2D(inputTexture, uv + vec2( 0, -1) / resolution.xy).rgb;14vec3 tx1y1 = texture2D(inputTexture, uv + vec2( 0, 0) / resolution.xy).rgb;15vec3 tx1y2 = texture2D(inputTexture, uv + vec2( 0, 1) / resolution.xy).rgb;16vec3 tx2y0 = texture2D(inputTexture, uv +