[css-color-4] CSS gamut mapping algorithm clarifications

[ccameron-chromium]

This is with respect to the definition of the CSS gamut mapping algorithm in CSS Gamut Mapping to an RGB Destination.

In particular, the following text:

For the analytical implementation, having found the exact intersection, project outwards (towards higher chroma) along the line of constant lightness until the deltaEOK between the projected point and a clipped version of that point exceeds one JND. Then returned the clipped version of the color as the mapped result.

In effect there are two steps: We start at the original point. We pull chroma in until we hit the border of the gamut, then we push chroma back out until clipping introduces only a certain level of error.

First question: Should we limit the amount we push chroma out in the second phase to never push beyond the original point?

I suspect so, but it might be worth adding text to the spec regarding this. The difference between these two implementations is visible in these images (created using this page):

• Walking out just to the original point
• Walking out possibly beyond the original point

Notice the artifacts that are visible in the blue area on the left.

Second question: Should we reconsider the trade-off of the "push chrome back out" step?

While the "push chroma back out" is motivated in the subsections titled "Excessive Chroma Reduction" and "Chroma Reduction with Local Clipping", the results when visualizing constant luminance suggest that there there is some degradation of results.

The difference can be seen in these two images

• Walking out just to the original point
• No walking out, staying on the gamut border

Notice the artifacts that are visible around the sRGB gamut when the walkout is applied versus when it isn't.

Also note that the "binary search implementation" is "allowed to" return points that are on the border, if the path of the binary search iteration happens to land there. These points (unpredictably located), would suffer the aforementioned chroma reduction that the "chroma walk out" part of the analytic algorithm is added to prevent.

Third question: Should we replace "binary search implementation" with something more precise and efficient?

The spec discusses two implementations, "binary search implementation" and "analytical implementation". I haven't found a code listing for the "analytical implementation", and I generally shy away from analytic solutions for nontrivial problems.

There do exist fast algorithms for computing the intersection of chroma with the border. For instance, this algorithm uses Newton's method to find the gamut boundary (up to floating-point precision), and requires just a few iterations (5 is a lot), and does nothing but adds and multiplies (and one divide) per iteration.

This might be a better algorithm to include in the code listing. I haven't looked into an efficient way to do the "chroma walk out" part of the algorithm, but perhaps we can remove it, depending on the resolution of the second question.

[danburzo]

Hi @ccameron-chromium, as I'm interested in exploring alternatives to the binary search method for gamut mapping, I thought I'd get the ball rolling discussing the issues you raise. Apologies in advance if I'm misconstruing anything.

First question: Should we limit the amount we push chroma out in the second phase to never push beyond the original point?

I agree it's fair to assume the end of the line is at the original chroma, no reason to look beyond.

Second question: Should we reconsider the trade-off of the "push chrome back out" step?

Looking at the provided demo page, it seems the algorithm implemented is a binary search in Oklch, although expressed as Oklab (instead of looking for Chroma, you're looking for an alpha coefficient for the a and b components), with improvements brought over from the analytical solution paragraph ("push chroma back out"), which sounds like combining the two alternative techniques?

Third question: Should we replace "binary search implementation" with something more precise and efficient?

A more precise and efficient solution is hard to argue against. I'll try to extract the salient part of the notebook as a JavaScript function so we can take a better look at the algorithm.

[danburzo]

This is my first pass at the algorithm:

Expand JavaScript code

const OKLAB_TO_LMS = [
	[0.99999999845051981432, 0.39633779217376785678, 0.21580375806075880339],
	[1.0000000088817607767, -0.1055613423236563494, -0.063854174771705903402],
	[1.0000000546724109177, -0.089484182094965759684, -1.2914855378640917399]
];

const LMS_TO_XYZ_D65 = [
	[1.2268798733741557, -0.5578149965554813, 0.28139105017721583],
	[-0.04057576262431372, 1.1122868293970594, -0.07171106666151701],
	[-0.07637294974672142, -0.4214933239627914, 1.5869240244272418]
];

const XYZ_D65_TO_LRGB = [
	[3.2409699419045226, -1.537383177570094, -0.4986107602930034],
	[-0.9692436362808796, 1.8759675015077202, 0.04155505740717559],
	[0.05563007969699366, -0.20397695888897652, 1.0569715142428786]
];

const LMS_TO_LRGB = multiplyMatrices(XYZ_D65_TO_LRGB, LMS_TO_XYZ_D65);

function multiplyMatrices(a, b) {
	let res = [[], [], []];
	for (let i = 0; i < 3; i++) {
		for (let j = 0; j < 3; j++) {
			res[i][j] = a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j];
		}
	}
	return res;
}

function multiplyMatrixWithVector(m, v) {
	return [
		m[0][0] * v[0] + m[0][1] * v[1] + m[0][2] * v[2],
		m[1][0] * v[0] + m[1][1] * v[1] + m[1][2] * v[2],
		m[2][0] * v[0] + m[2][1] * v[1] + m[2][2] * v[2]
	];
}

function multiplyComponentWise(arr1, arr2) {
	return arr1.map((it, i) => it * arr2[i]);
}

function multiplyVectors(arr1, arr2) {
	return arr1.reduce((acc, curr, i) => {
		return acc + curr * arr2[i];
	}, 0);
}

export function GamutMapNewton(original) {
	const zero_a_b = [0, original[1], original[2]];
	let alpha = 1;
	let res_oklab = original.slice();
	let res_rgb = multiplyMatrixWithVector(
		LMS_TO_LRGB,
		multiplyMatrixWithVector(OKLAB_TO_LMS, res_oklab).map(it => it ** 3)
	);
	for (let comp = 0; comp < 3; comp++) {
		if (res_rgb[comp] >= 0 && res_rgb[comp] <= 1) continue;
		let target = res_rgb[comp] > 1 ? 1 : 0;
		for (let iter = 0; iter < 6; iter++) {
			let residual =
				multiplyVectors(
					LMS_TO_LRGB[comp],
					multiplyMatrixWithVector(OKLAB_TO_LMS, res_oklab).map(it => it ** 3)
				) - target;
			if (Math.abs(residual) < 1e-15) break;
			let gradient = multiplyVectors(
				LMS_TO_LRGB[comp],
				multiplyComponentWise(
					multiplyMatrixWithVector(OKLAB_TO_LMS, res_oklab).map(it => 3 * it ** 2),
					multiplyMatrixWithVector(OKLAB_TO_LMS, zero_a_b)
				)
			);
			alpha -= residual / gradient;
			res_oklab[1] = alpha * original[1];
			res_oklab[2] = alpha * original[2];
		}
		res_rgb = multiplyMatrixWithVector(
			LMS_TO_LRGB,
			multiplyMatrixWithVector(OKLAB_TO_LMS, res_oklab).map(it => it ** 3)
		);
	}
	return [res_oklab, res_rgb];
}

I've updated my gamut mapping demo page to include this new algorithm (see the newton (@ccameron) method).

There seem to be some visual glitches outside the sRGB gamut, so I can't exclude translation errors on my part either in the code above, or in the glue code for the demo page, but I hope it's a good starting point.

Assuming the algorithm is sound and the glitches can be ironed out, what I'm observing is that it's very close to (indistinguishable from) basic binary search (the chroma-reduce method) when rendering lch() / oklch() / etc. canvases, but seems to converge faster to the solution. (I'll update the demo page to include lab() / oklab() canvases)

[ccameron-chromium]

Second question: Should we reconsider the trade-off of the "push chrome back out" step?

Looking at the provided demo page, it seems the algorithm implemented is a binary search in Oklch, although expressed as Oklab (instead of looking for Chroma, you're looking for an alpha coefficient for the a and b components)

Yes, this is just to avoid switching between polar and rectangular coordinates.

with improvements brought over from the analytical solution paragraph ("push chroma back out"), which sounds like combining the two alternative techniques?

This paragraph complicates things a lot.

For the analytical implementation, having found the exact intersection, project outwards (towards higher chroma) along the line of constant lightness until the deltaEOK between the projected point and a clipped version of that point exceeds one JND. Then returned the clipped version of the color as the mapped result.

The text in the "Excessive Chroma Reduction" and "Chroma Reduction with Local Clipping" sections gives a motivation for this, but I would like to see concrete examples of the form "gamut map color X to gamut G with and without this step and see the difference".

Separately, the term "analytical" should be avoided. It generally refers to closed-form solutions (as opposed to iterative solutions), which are very rarely desirable. The Newton-based solution I mentioned is not analytic -- it's in the same class of iterative solutions as binary search. Both solutions will converge to the exact solution (up to floating point precision), the only difference is the number of steps that are needed. The idea of adding an extra step to solutions that arrive at the exact solution "too quickly" seems inappropriate and indicative of deeper problems.

I would like to write up the code that is also capable of finding the point on the surface of the gamut that is closest to the input point in a weighted L2 sense. So, if the input point is (L,a,b), find the point (L*,a*,b*) that is within the gamut (or on its boundary) such that we minimize their distance (L-L*)^2+(a-a*)^2+(b-b*)^2.

Again, various text gives motivation for the "only change chroma" scheme, but I would want to see concrete examples where it demonstrates superiority. In the example cases that I've tried (e.g, "red-redder" from here), I am not convinced that the result is desirable.

[svgeesus]

Yes, this is just to avoid switching between polar and rectangular coordinates.

Its a nice optimization and gives the same result.

Separately, the term "analytical" should be avoided. I

Agreed, it was not a good choice. I guess I meant "geometrically" ie the intersection of a line and a polyhedron.

The idea of adding an extra step to solutions that arrive at the exact solution "too quickly" seems inappropriate and indicative of deeper problems.

No, its an awareness of a problem, which is fairly well known in the gamut mapping literature, caused by shallow and concave gamut boundaries. Instead of intersecting a (zero-width) line with the boundary, we actually want to intersect a (0.5 JND radius) cylinder.

[danburzo]

I have deployed another demo page to illustrate the Oklab color space specifically, when gamut mapping is performed in Oklch with a JND (just-noticeable difference) of 0.02: Gamut-mapping Oklab to sRGB

You can compare:

• clip: simple clipping
• chroma-reduce: maximum chroma along the line of constant lightness & hue, obtained through binary search
• culori.toGamut(), fuzzy and css-color-4-smooth are three algorithms with minute, inconsequential variations on the theme of changing the binary search test to match "roughly in gamut" (explained below)
• css-color-4 is the current binary search algorithm defined in the spec (no longer relevant, spec was changed)

The "roughly in gamut" test changes the "in gamut" test from each step of the binary search from if (inGamut(color)) to if (inGamut(color) || delta(color, clipped(color)) < JND), which I believe to be the spirit of the adjustments suggested by the spec: when the line of constant lightness & hue hovers just above the gamut boundary, you can grab a nearby in-gamut color if it's not too dissimilar and gives you better chroma.

The idea of adding an extra step to solutions that arrive at the exact solution "too quickly" seems inappropriate and indicative of deeper problems.

The adjustment proposed can be readily baked into the binary search, but is not at all a factor for other methods, e.g. intersecting the line with a polygonal / polyhedral gamut boundary. (Is it applicable to the Newton method? I don't have an intuition for this).

I would like to write up the code that is also capable of finding the point on the surface of the gamut that is closest to the input point in a weighted L2 sense. […] Again, various text gives motivation for the "only change chroma" scheme, but I would want to see concrete examples where it demonstrates superiority.

Experimentally, people tended to prefer gamut mapping with these dividing weights (lightness/chroma/hue): [1, 2.6, 1.3] (as per my reading of Ján Morovič's Color Gamut Mapping), that is: chroma adjustment is most tolerable, followed by hue; people judged changes in lightness least favorably. So true MINDE with a weighted Euclidean distance in a perceptually-uniformed color space would probably give the best results, but is probably more involved algorithmically. (Morovič gives an algorithm that makes use of gamut boundary polygons of constant hue.)

Chroma reduction with constant hue and lightness is straightforward algorithmically and acceptable perceptually, but hardly 'superior' (unless you compare it to naive clipping). As @svgeesus notes, its shortcomings are well-documented. Björn Ottosson has experimented with keeping the hue constant & projecting along different lines here, but evaluation was done on photographic images, which have different trade-offs.

So we're actually talking several questions:

1. Is chroma reduction the best overall gamut mapping technique for the spec to provide?
2. Is binary search the best algorithm to express the technique?
3. Is the 'roughly in gamut' enhancement a net positive addition to the algorithm? Can it be expressed in terms of the algorithm provided?

[danburzo]

In terms of clarifying the prose, I think it would be beneficial for section 13.1.5. Chroma Reduction with Local Clipping to explain the algorithm as finding the point of maximum chroma along the segment between the color and its achromatic version that's either in gamut or has a clipped version that is not noticeably different (ΔE < JND).

Section 13.2. CSS Gamut Mapping to an RGB Destination could then explain how the adjustment can be baked into the binary search method and how other (iterative/analytical) methods may need to compensate by "pushing chroma out" — which, by lack of mathematical imagination, I reckon is still some sort of binary search, unless it is similarly amenable to faster methods.

[danburzo]

A more concerning effect of chroma reduction in Oklch is the series of unexpected artifacts produced when gamut-mapping lch() colors, the cause of which would be interesting to trace.

[ccameron-chromium]

This is my first pass at the algorithm:

There seem to be some visual glitches outside the sRGB gamut, so I can't exclude translation errors on my part either in the code above, or in the glue code for the demo page, but I hope it's a good starting point.

I encountered this as well, and found that the problem was that I wasn't starting Newton iteration from near enough to the surface to converge. This is particularly bad in the non-convex areas of the the image the various gamuts in oklab space.

The solution that I started towards was to first intersect with a polyhedral approximation of the gamut and then refine using Newton iteration.

(I also found that the polyhedral approximation is extremely close to the gamut and might be better).

[Artoria2e5]

Umm @danburzo, do you mind adding the geometric (formerly "analytical") solutions of BjO to the comparison too? I have ham-fistedly translated his code to TS at https://gist.github.com/Artoria2e5/9c7ba0bcda480b5bc2ae0b0ffe0bfb91 for a use-and-throw-away color-picking job; you can probably do a lot better with the multiplyMatrices function already found in the Newton thing.

It's not cool to talk about a geometric solution existing without ever showing an example of it, IMO. Sure this one's got very strong dependence on the shape of sRGB in particular, but... at least we can figure out whether the result is more desirable.

[danburzo]

@Artoria2e5 There are more gamut mapping methods implemented in the color.js Gamut Mapping Playground (repo), but not all implemented in the library itself.