Studying Gamut Clipping

Introduction

Recently, I was studying a technique called gamut clipping from this blog post. This technique is used for handling out of gamut color and bring them back to a valid range, which helps reducing hue shift and color clipping. From that blog post, it explains the concept clearly, but I was struggling to understand the sample code the author provided. So this blog post will describe what I have learnt when studying the gamut clipping source code. Also, I have written a Shadertoy sample to help me to understand the problem.

Drawing the valid sRGB color with Oklab chroma as horizontal axis
Oklab lightness as vertical axis with hue value displayed at upper right corner

Overview of the gamut clipping

This section briefly describe the steps to perform gamut clipping, feel free to skip it if you have read the original gamut clipping blog post already. The technique first start to convert the out of gamut color (e.g. those pixels with values >1.0 or <0.0 in sRGB) into the Oklab color space (which can be expressed with lightness, hue and chroma). Then we try to project this out of gamut color along a straight line to the"triangular" gamut boundary (like below picture). To calculate this intersection, we need to find the cusp coordinates at this particular hue slice. The author use curve fitting to approximate cusp coordinates with polynomial equation. With the cusp coordinates calculated, the intersection point can be found with numerical approximation using Halley's method. Finally that intersection point can be converted back to sRGB color space.

Showing an out of gamut color point being projected back to a valid value along a straight line.	Displaying all color including those sRGB clipped color, out of gamut color will result in hue shift (e.g. blue hue displayed as purple).

Finding maximum saturation

To approximate the cusp, instead of using the hue value and returning the cusp (chroma, lightness) coordinates directly. The author tries to use the maximum saturation value to find the cusp coordinates. This is the first thing I don't understand when reading the source code. Why saturation is related to the cusp? To understand this, we can start from the definition of saturation first. Both terms chroma and saturation describe colorfulness, but saturation is "somehow normalized" and "not affected" by the brightness. A chroma/lightness/saturation relationship can be defined like this:

Saturation definition for CIELAB space from Wiki,
we use the same definition for Oklab space too.

Reading the author's Colab source code, he tried to approximate the max saturation with a polynomial function using a=cos(hue) and b=sin(hue):

Then these (saturation, hue) values will be converted to linear sRGB space and minimize it for R=0, G=0 and B=0, which obtained 3 polynomial equations to approximate the maximum saturation.

But why we can obtain maximum saturation when either R / G / B=0? To have a rough understanding of it, I decided to open a color picker with HSL slider to "simulate" how the RGB values changes:

A color picker with maximum saturation and changing hue slowly.
There will always be a 0 in the RGB value.

First I choose the most saturated red color (255, 0, 0) in sRGB, which yield a HSL value (0, 240, 120). Then I change the hue value slowly to observe how the RGB value changes. From the above animated gif, we can see that there will always be a 0 in either R or G or B channel. So I guess the author is using this property and the Oklab space has similar property.

To further understand how the saturation looks like for all hue slices in Oklab, I plotted the saturation value, seems like it always have the largest value at the lower "triangle" edge:

We can also draw the line for those pixels with R/G/B value = 0:

The lower "triangle" edge is actually the "clipping" lines when R=0 or G=0 or B=0 (and switching between these 3 lines)! That's why the author tried to use 3 different polynomials to approximate the max saturation. Then the next problem is how to pick one of the three approximated polynomials. From the above animated gif, we know that the "clipping line" changes when sRGB values = (1, 0, 0), (0, 1,0) or (0, 0, 1). We can see this from the Colab sample code which find the r_dir/g_dir/b_dir, and use this to pick one of the three approximated polynomials.

The r_dir/g_dir/b_dir is the "half vector" between the "opposite hue" like the below figure. And this direction vector can be used with dot product to check whether a hue value (in form of cos(hue), sin(hue)) is lying within which r_h/g_h/b_h range.

Visualizing the r/g/b_dir' vector (before multiplying the scalar constant for fast dot product check).

The formula for checking which hue range the a, b hue value is in, which can be derived from dot product.

Another problem I don't understand is: when the author tried to minimize the RGB value to 0, he raised the error value to the power 10 in the e_R(x)/e_G(x)/e_B(x) function.

When I tried to throw random values to the initial guess value in scipy.optimize.minimize() (e.g. use all 0 or 1 as initial guess), the resulting approximated curve is not that good...

However when I changed the pow(x, 10) to abs(x), using random initial guess values can have a more predictable result (though the error is not as small as the author's approximation). May be I will use abs(x) instead to optimize for different color spaces in the future.

Finding cusp from maximum saturation

With the above polynomial approximation, we can find the "lower clipping line" of the gamut "triangle" (i.e. with either R/G/B value = 0). The cusp must be on this "clipping line". We only need to identify the lightness value at the cusp. Looking at the author's find_cusp() function, it seems like the cusp will have R or G or B value equals to 1:

To have a better understanding, we can repeat the above color picker "experiment", when changing the hue value at maximum saturation, beside there always having a 0 value in R/G/B, there will always having a 255 value in R/G/B too!

The same color picker gif as above,
pay attention to the RGB values,
in addition to always having a 0 value,
there is a 255 value too.

Plotting the line when R/G/B value = 1 for all hue slices in Oklab space:

When R/G/B value = 1, it is the upper clipping line of the "gamut triangle"! So the cusp is the intersection between lower clipping line with R/G/B=0 and the upper clipping line with R/G/B=1.

So, to obtain the maximum lightness and maximum saturation, we can scale the lightness so that when converted back to sRGB, one of the RGB value will = 1 (taking cubic root is needed because of the Oklab color space LMS definition).

Finding gamut intersection

With the gamut cusp coordinates, target project coordinates and out of gamut coordinates, we can find the intersection coordinates at the gamut boundary during projection. First we need to check whether the intersection is in the upper half or the lower half of the valid range gamut "triangle". This can be determined by below formula which can be derived from checking whether the projection line is on the left/right side to the cusp line using cross product.

If in the lower half, it is just a simple line-line intersection. If in the upper half, we first approximate the intersection with line-line intersection, and then refine the answer with Halley's method. Since we know the upper clipping line is R/G/B=1, we can see the author is using this property when using Halley's method in the following code snippet:

and the author use all 3 clipping curves and take the minimum value to select the closest curves that clip the gamut.

If small error is acceptable, we can use 1 clipping curve instead of 3. The upper clipping curve changes happens roughly at RGB value=(0, 1, 1), (1, 0, 1) and (1, 1, 0). We can use the same method as picking the R/G/B=0 curve during maximum saturation calculation which relies on the r_dir, g_dir, b_dir. We can compute the c_dir, m_dir and y_dir (which correspond to the cyan, magenta and yellow direction). Those coefficients can be found in my Shadertoy source code. Because the upper clipping curve is not straight line, we may need 2 clipping lines to compute the correct answer for some hue values:

Need 2 upper clipping lines for this hue slice

Conclusion

In this post, I have described the process of learning the gamut clipping technique. With the help of the Shadertoy sample, we can see that the gamut boundary is the line with R/G/B=0 and R/G/B=1. And the gamut cusp is the intersection between R/G/B=0 line and R/G/B=1 line. But I still have some questions left in my mind like: Is it meaningful to use negative chroma values? Is the gamut clipping works well in pratice, or need a gamut compression method? May be I will implement gamut clipping in my toy path tracer to see how it looks like.

Zooming out the graph, do the negative values (e.g. -ve chroma, -ve lightness) have meaning?
The G=1 clipping line "wrap around" to negative lightness, does it have meaning too?

Reference

[1] https://bottosson.github.io/posts/gamutclipping/

[2] https://bottosson.github.io/posts/oklab/

[3] https://en.wikipedia.org/wiki/Colorfulness#CIELUV_and_CIELAB

Importance sampling visible wavelength

Introduction

It has been half a year since my last post. Due to the pandemic and political environment in Hong Kong, I don't have much time/mood to work on my hobby path tracer... And until recently, I tried to get back to this hobby, may be it is better to work on some small tasks first. One thing I am not satisfied in previous spectral path tracing post is using 3 different cosh curves(with peak at the XYZ color matching functions) to importance sample the visible wavelength instead of 1. So I decided to revise it and find another PDF to take random wavelength samples. A demo with the updated importance sampled wavelength can be downloaded here and the python code used for generating the PDF can be viewed here (inspired by Bart Wronski's tweet to use Colab).

First Failed Try

The basic idea is to create a function with peak values at the same location as the color matching functions. I decided to use the sum of the XYZ color matching curves as PDF.

Black line is the Sum of XYZ curves

To simplify calculation, analytical approximation of the XYZ curves can be used. A common approximation can be found here, which seems hard to be integrated(due to the λ square) to create the CDF. So a rational function is used instead:

The X curve is approximated with 1 peak and dropped the small peak at around 450nm, because we want to compute the "sum of XYZ" curve, that missing peak can be compensated by scaling the Z curves. The approximated curves looks like:

Light colored curves are the approximated function

Grey color curve is the approximated PDF

The approximated PDF is roughly similar to the sum of XYZ color matching curves. But I have made a mistake: Although the rational function can be integrated to create the CDF, I don't know how to compute the inverse of CDF(which is needed by the inverse method to draw random samples using uniform random numbers). So I need to find another way...

Second Try

From previous section, although I don't know how to find the inverse of the approximated CDF, out of curiosity, I still want to know how the approximated CDF looks like, so I plot the graph:

Black line is the original CDF
Grey line is the approximated CDF

It looks like some smoothstep functions added together, with 1 base smoothstep curve in range [380, 780] with 2 smaller smoothstep curves (in range around [400, 500] and [500, 650]) added on top of the base curve. May be I can approximate this CDF with some kind of polynomial function. After some trail and error, an approximated CDF is found (details of CDF and PDF can be found in the python code):

The above function divides the visible wavelength spectrum into 4 intervals to form a piecewise function. Since smoothstep is a cubic function, adding another smoothstep function is still another cubic function, which can be inverted and differentiated. And the approximated "smoothstep CDF/PDF" curve looks like:

Blue line is the approximated "smoothstep CDF" Black line is the original CDF

Blue line is the approximated "smoothstep PDF" Black line is the original PDF

Although the "smoothstep CDF" looks smooth, its PDF is not (not C1 continuous at around 500nm). But it has 2 peaks value at around 450nm  and 600nm, may be let's try to render some images to see how it behaves.

Result

The same Sponza scene is rendered with 3 different wavelength sampling functions: uniform, cosh and smoothstep (with stratification of random number disabled, color noise are more noticeable when zoomed in):

uniform weighted sampling

cosh weighted sampling

smoothstep weighted sampling

Both cosh and smoothstep wavelength sampling method show less color noise than uniform sampling method, with the smoothstep PDF slightly better than cosh function. Seems the C1 discontinuity of the PDF does not affect rendering very much. A demo can be downloaded here to see how it looks in real-time.

Conclusion

This post described an approximated function to importance sample the visible wavelength by using the sum of the color matching functions, which reduce the color noise slightly. The approximated CDF is composed of cubic piecewise functions. The python code used for generating the polynomial coefficients can be found here (with some unused testing code too. e.g. I have tried to use 1 linear base function with 2 smaller smoothstep functions added on top, but the result is not much better...). Although the approximated PDF is not C1 continuous, it does not affect the rendering very much. If someone knows more about how the C1 continuity of the PDF affect rendering, please leave a comment below. Thank you.

References

[1] https://en.wikipedia.org/wiki/CIE_1931_color_space#Color_matching_functions

[2] Simple Analytic Approximations to the CIE XYZColor Matching Functions

[3] https://stackoverflow.com/questions/13328676/c-solving-cubic-equations

sRGB/ACEScg Luminance Comparison

Introduction

When I was searching information about rendering in different color spaces, I came across that using wider color primaries  (e.g. ACEScg instead of sRGB/Rec709) to perform lighting calculation will have a result closer to spectral rendering. But will this affect the brightness of the bounced light? So I decided to find it out. (The math is a bit lengthy, please feel free to skip to the result.)

 

Comparison method

To predict the brightness of the rendered image, we can consider the reflected light color after n bounces. To simplify the problem, we assume all the surfaces are diffuse material. We can derive a formula for the RGB color vector c after n light bounces.

To calculate lighting in different color spaces, we need to convert the albedo and initial light color to our desired color gamut by multiplying with a matrix M (For rendering in sRGB/Rec709, M is identity matrix).

Finally, we can calculate the luminance of the bounced light by computing the dot product between the color vector c and luminance vector Y of the color space (i.e. Y is the second row vector of the conversion matrix from a color space to XYZ space, with chromatic adaptation).

Now, we have an equation to compute the brightness of a reflected ray after n bounces in arbitrary color space.

Grey Material Test

To further simplify the problem, we assume all the surfaces are using the same material:

Then assuming all the surfaces are grey in color:

Now, the luminance equation is simpler to understand.

Substituting the matrix M and luminance vector Y for sRGB color gamut, the equation is very simple:

Then we do the same thing for ACEScg. Surprisingly, there are some constants roughly equal to one, so we can approximate them with constant one and the result will roughly equal to the luminance equation of sRGB.

As both equations are roughly equal, the rendered images in sRGB and ACEScg should be similar. Let's try to render images in sRGB and ACEScg to see the result (images are path traced with sRGB and ACEScg primaries, and then displayed in sRGB).

Path traced in sRGB

Path traced in ACEScg

Both images looks very similar! So rendering in different color spaces with grey material will not change the brightness of the image. At least, the difference is very small after tone mapped to a displayable range.

Red Material Test

Now, let's try to use red material instead of grey material to see how the luminance changes (where k is a variable to control how 'red' the material is):

But the equation is still a bit complex,  so we further assume the initial light color is white in color.

Then we perform the same steps in last section, substituting M and Y into luminance equation.

sRGB luminance equation

ACEScg luminance equation

Unfortunately, both equations are a bit too complex to compare, having 2 variables k and n... May be let's try to plot some graphs to see how those variables affect the luminance, with number of bounced light = 3 and 5 (i.e. n=3 and n=5, skipping the N dot L part because both equations have such term). From the graphs below: when k increases (i.e. the red color is getting more saturated, with RGB value closer to (1, 0, 0) ), luminance difference will increase, hence sRGB will have a larger luminance value than ACEScg.

Then comparing the images rendered in sRGB and ACEScg:

Path traced in sRGB

Path traced in ACEScg

The indirectly lit area looks much brighter when rendered in sRGB. This makes sense because for any red color, its red channel value will be closer to one (while green/blue values will be closer to 0) when represented in sRGB compared to be represented in ACEScg. After several multiplication, the reflected light value should be larger when computation is done in sRGB.

 

RGB Material Test

How about using different colored material this time? Assuming 1/3 of the light bounced on surface that is red material, 1/3 is green material, and 1/3 is blue material. 

Like previous 2 sections, substituting M and Y, the luminance equation becomes:

sRGB luminance equation	ACEScg luminance equation

 

And then plotting graphs to see how k and n varies:

The result is different this time. The sRGB luminance is smaller than ACEScg luminance, and the difference increases when both k and n increases. So the bounced light will be darker when rendered in sRGB.

Let's try rendering some images to see whether this is true. Although we cannot force the ray to bounce with 1/3 red/green/blue material exactly, I roughly assigned 1/3 of the material in red/green/blue color in the scene.

Path traced in sRGB

Path traced in ACEScg

From the screen shots above, the indirectly lit red material looks darker when rendered in sRGB (especially the curtains on the ground floor), while the differences for the blue and green material are small. We can think of the result like previous section, for a given red color, when it is represented in sRGB, its red channel value is closer to one, however, its blue and green channels values are closer to 0 compared to be represented in ACEScg (same for both blue and green material). So after several multiplication with different color material, the RGB values in sRGB may becomes closer to 0 because different color material cancel out each other (e.g. when light bounced on red and green albedo surface (1, 0, 0) , (0, 1, 0) in sRGB, the reflected light will be zero, while the same color represented in ACEScg, the color will not be "zeroed out"), resulting in darker image.

Conclusion

After testing with different assumptions, the brightness of images when rendered in sRGB can be darker, roughly equal or brighter than rendered in ACEScg. The brightness difference depends on material used in the scene. If the scene uses grey material only, brightness will be equal. If material has similar color (e.g. all red material), sRGB image will be brighter. If the scene has more color variation, the sRGB image may becomes darker. And turns out this conclusion can be arrived without doing such lengthy math. We can think of the same color value represented in sRGB and ACEScg space: Is the RGB value closer to 0 or 1 when represented in the color space? Will the RGB values 'cancel' each other when performing lighting calculation in the color space? So I was too stupid to figure out this simple answer early and instead worked on such lengthy math... >.<

Reference

[1] https://chrisbrejon.com/cg-cinematography/chapter-1-5-academy-color-encoding-system-aces/