Calc more accurate inverse matrices for cat and pre-calc D50 <-> D65

[facelessuser]

When applying the inverse CAT transform, the inverse should be calculated at double precision to match the precision of the values used.

[facelessuser]

Values were calculated with numpy. For reproduction, here is a script:

"""Calculate inverse CAT matrices"""
import numpy as np

np.set_printoptions(precision=16, sign='-', floatmode='fixed')

D50 = [0.3457 / 0.3585, 1.00000, (1.0 - 0.3457 - 0.3585) / 0.3585]
D65 = [0.3127 / 0.3290, 1.00000, (1.0 - 0.3127 - 0.3290) / 0.3290]

bradford = [
    [0.8951000, 0.2664000, -0.1614000],
    [-0.7502000, 1.7135000, 0.0367000],
    [0.0389000, -0.0685000, 1.0296000]
]

von_kries = [
    [0.4002400, 0.7076000, -0.0808100],
    [-0.2263000, 1.1653200, 0.0457000],
    [0.0000000, 0.0000000, 0.9182200]
]

cat02 = [
    [0.7328000, 0.4296000, -0.1624000],
    [-0.7036000, 1.6975000, 0.0061000],
    [0.0030000, 0.0136000, 0.9834000]
]

cat16 = [
    [  0.401288,  0.650173, -0.051461 ],
    [ -0.250268,  1.204414,  0.045854 ],
    [ -0.002079,  0.048952,  0.953127 ]
]

names = [
    'bradford',
    'von_kries',
    'cat02',
    'cat16'
]

if __name__ == '__main__':
    for cat in names:
        m = globals().get(cat)
        print(f'===== {cat} =====')
        print(m)
        print(f'===== {cat} inverse =====')
        print(np.linalg.inv(np.asarray(m, dtype=np.double)))

The recalculated D50 <--> D65 translation was dumped directly from Color.js via:

export function adapt (W1, W2, id = "Bradford") {
	// adapt from a source whitepoint or illuminant W1
	// to a destination whitepoint or illuminant W2,
	// using the given chromatic adaptation transform (CAT)
	// debugger;
	let method = CATs[id];

	let [ρs, γs, βs] = multiplyMatrices(method.toCone_M, W1);
	let [ρd, γd, βd] = multiplyMatrices(method.toCone_M, W2);

	// all practical illuminants have non-zero XYZ so no division by zero can occur below
	let scale = [
		[ρd/ρs,    0,      0      ],
		[0,        γd/γs,  0      ],
		[0,        0,      βd/βs  ]
	];
	// console.log({scale});

	let scaled_cone_M = multiplyMatrices(scale, method.toCone_M);
	let adapt_M	= multiplyMatrices(method.fromCone_M, scaled_cone_M);
	console.log({scaled_cone_M, adapt_M});
	return adapt_M;
};

[LeaVerou]

Hey, thanks for working on this!

I'm wondering if it would make sense to have a build process that calculates these and stores the result in a simple ESM file? Then we could change the precision at will. The downside is it would require someone porting the python code to JS, so it requires more effort. Maybe something for the backlog.

[facelessuser]

@LeaVerou Forgive me, while capable in JS, it is not the language I usually work with. I saw a similar script living in the repo https://github.com/LeaVerou/color.js/blob/main/scripts/RGB_matrix_maker.py, so I figured it was fine reaching for what I know.

I can update with a JS version. A little research turned up this project which seems capable of doing linear algebra: https://mathjs.org/docs/getting_started.html.

I'll provide the script in Scripts. Update dev dependencies to include math.js, then update the matrices with the results from the script. I will likely let someone familiar with the build environment update any automated processes. Is that okay?

So, using that lib, here is a script:

import { create, all } from 'mathjs'
const config = {}
const math = create(all, config)

const bradford = [
    [0.8951000, 0.2664000, -0.1614000],
    [-0.7502000, 1.7135000, 0.0367000],
    [0.0389000, -0.0685000, 1.0296000]
]

const von_kries = [
    [0.4002400, 0.7076000, -0.0808100],
    [-0.2263000, 1.1653200, 0.0457000],
    [0.0000000, 0.0000000, 0.9182200]
]

const cat02 = [
    [0.7328000, 0.4296000, -0.1624000],
    [-0.7036000, 1.6975000, 0.0061000],
    [0.0030000, 0.0136000, 0.9834000]
]

const cat16 = [
    [  0.401288,  0.650173, -0.051461 ],
    [ -0.250268,  1.204414,  0.045854 ],
    [ -0.002079,  0.048952,  0.953127 ]
]

console.log('===== bradford =====')
console.log(bradford)
console.log('===== bradford inverse =====')
console.log(math.inv(bradford))

console.log('===== von kries =====')
console.log(von_kries)
console.log('===== von kries inverse =====')
console.log(math.inv(von_kries))

console.log('===== cat02 =====')
console.log(cat02)
console.log('===== cat02 inverse =====')
console.log(math.inv(cat02))

console.log('===== cat16 =====')
console.log(cat16)
console.log('===== cat16 inverse =====')
console.log(math.inv(cat16))

Here are the results:

===== bradford =====
[
  [ 0.8951, 0.2664, -0.1614 ],
  [ -0.7502, 1.7135, 0.0367 ],
  [ 0.0389, -0.0685, 1.0296 ]
]
===== bradford inverse =====
[
  [ 0.9869929054667121, -0.14705425642099013, 0.15996265166373122 ],
  [ 0.4323052697233945, 0.5183602715367774, 0.049291228212855594 ],
  [ -0.00852866457517732, 0.04004282165408486, 0.96848669578755 ]
]
===== von kries =====
[
  [ 0.40024, 0.7076, -0.08081 ],
  [ -0.2263, 1.16532, 0.0457 ],
  [ 0, 0, 0.91822 ]
]
===== von kries inverse =====
[
  [ 1.8599363874558397, -1.1293816185800916, 0.21989740959619328 ],
  [ 0.3611914362417676, 0.6388124632850422, -0.000006370596838649899 ],
  [ 0, 0, 1.0890636230968613 ]
]
===== cat02 =====
[
  [ 0.7328, 0.4296, -0.1624 ],
  [ -0.7036, 1.6975, 0.0061 ],
  [ 0.003, 0.0136, 0.9834 ]
]
===== cat02 inverse =====
[
  [ 1.0961238208355142, -0.27886900021828726, 0.18274517938277307 ],
  [ 0.4543690419753592, 0.4735331543074117, 0.07209780371722911 ],
  [ -0.009627608738429355, -0.00569803121611342, 1.0153256399545427 ]
]
===== cat16 =====
[
  [ 0.401288, 0.650173, -0.051461 ],
  [ -0.250268, 1.204414, 0.045854 ],
  [ -0.002079, 0.048952, 0.953127 ]
]
===== cat16 inverse =====
[
  [ 1.862067855087233, -1.0112546305316845, 0.14918677544445172 ],
  [ 0.3875265432361372, 0.6214474419314753, -0.008973985167612521 ],
  [ -0.01584149884933386, -0.03412293802851557, 1.0499644368778496 ]
]

[facelessuser]

I updated with a JS solution and included the script along with the necessary math.js dependency. My aim is mainly to fix the matrices. Any kind of refactoring, if desired, should be possible using math.js. I've provided the reproduction script along with the dependency.

[facelessuser]

I will state that there may be better JS matrix libraries out there, if a better one is discovered, I completely understand. If this project would like to roll its own inverse method, it is completely doable, though a bit more involved. I've done so with my own project as I do not currently rely on numpy, but implement matrix inverse doing LU decomposition. So a simplied version could be ported here, but that would most likely be something for another PR.

[svgeesus]

Then we could change the precision at will.

I see no reason to do this. The value we are currently using was calculated at lower precision; this PR does it at higher precision and we will use that.

[svgeesus]

Thanks for doing this!

[romainmenke]

These are very different values.

The other changes seem to be the same value at higher precision, but this one is completely different.

Is this correct?

When we use these values in postcss-preset-env it dramatically changes the color results.

[facelessuser]

@romainmenke great catch! Accidentally copied in the adapted M cone response as the inverse. I have a fix here: #360. This one copies the correct matrix, directly calculated by color.js when a D50 <-> D65 conversion occurs.