[BUG]: sCAM assumes a white point that contradicts with the IPT LMS matrix that is being used

[facelessuser]

Description

I was looking over the implementation paper of sCAM, and was verifying my results against Colour Science's, and I am a little confused about one point.

sCAM uses the IPT LMS transformation matrix. The sCAM paper states that the color should be in D65 when converting to Iab, but does not state what D65 is. What is interesting is that it states the HPE matrix is used, but what is listed is the modified one that IPT uses, which was specifically adapted for a white point of [0.9504, 1.0, 1.0889] (which is specifically stated in the IPT paper).

TO_LMS = [
    [0.4002, 0.7075, -0.0807],
    [-0.2280, 1.1500, 0.06124],
    [0.0, 0.0, 0.9184]
]

FROM_LMS = LMS_TO_XYZ = [
    [1.8502, -1.1383, 0.2385],
    [0.3668, 0.6439, -0.0107],
    [0.0, 0.0, 1.0889]
]

I will note that the IPT paper specifies 1.8501 in the inverse matrix, but this appears to be a mistake, as the inverse clearly yields a 1.8502.

>>> import numpy as np
>>> m = [
...     [0.4002, 0.7075, -0.0807],
...     [-0.2280, 1.1500, 0.06124],
...     [0.0, 0.0, 0.9184]
... ]
>>> np.linalg.inv(m)
array([[ 1.85024294, -1.13830164,  0.23848454],
       [ 0.36683078,  0.64388454, -0.01070149],
       [ 0.        ,  0.        ,  1.08885017]])

Now the forward transform creates noise in anything greater than 16-bit due to rounding.

>>> w = [0.9504, 1.0, 1.0889]
>>> np.matmul(m, w)
array([0.99997585, 0.99999304, 1.00004576])

But if we used the inverse transform, we can see that it yields a good LMS transform of [1, 1, 1] to the expected white point past 16-bit. It is clear that the transformation matrix was tailored to [0.9504, 1.0, 1.0889].

>>> mi = [
...     [1.8502, -1.1383, 0.2385],
...     [0.3668, 0.6439, -0.0107],
...     [0.0, 0.0, 1.0889]
... ]
>>> np.matmul(mi, [1, 1, 1])
array([0.9504, 1.    , 1.0889])
>>> np.matmul(np.linalg.inv(mi), [0.9504, 1.0, 1.0889])
array([1., 1., 1.])

Colour Science seems to assume a white point of:

TVS_D65_sCAM = np.array([0.95047, 1.00000, 1.08883])

And while I am aware this is a commonly used white point, and am aware of how it is derived, it seems to be an odd choice for the LMS matrix being used.

Is [0.95047, 1.00000, 1.08883] truly the white point that should be used? Or should it be [0.9504, 1.0, 1.0889]?

Currently, even when inputting [0.95047, 1.00000, 1.08883] into the sCAM conversion, you will not get 100% lightness in sCAM 1) due to the reverse LMS transform being generated from the forward transform which generates noise in anything greater than 16-bit, and 2) the LMS transform is not meant for the white point being used (unless it was adapted to the D65 variant being used).

If the appropriate white point (in the context of the IPT LMS matrix) was used, and the forward transform was generated from the reverse transform, you would then yield a clean 100% lightness in sCAM for the same white if it was the input.

Anyway, I thought it was worth bringing up the topic.

Code for Reproduction

Exception Message

Environment Information

[Jackchou00]

Hi @facelessuser ,

Thanks for the detailed analysis. My understanding is as follows:

Since you can't pre-define a D65 white point and a forward matrix at four or five effective digits to yield a precise [1, 1, 1] LMS response. That's maybe a limitation for all color model in paper (because you always need to give a clear matrix in paper), and the model development may not dive that deeply into such tiny precision errors.

I think it will be interesting to use a pre-defined inverse matrix and use [1, 1, 1] to yield the optimal input D65 white point (like [0.9504, 1, 1.0889] in IPT), just like what you have done above. Then, the forward matrix can be created using numpy’s inv function from the inverse matrix to avoid precision issues in the pre-defined matrix. This idea is worth learning for me, and the affect for model performance should be very minor.

That's a great perspective from engineering and industry (I guess), I’d really appreciate hearing thoughts from experts and thank you again.

Jack.

[facelessuser]

Since you can't pre-define a D65 white point and a forward matrix at four or five effective digits to yield a precise [1, 1, 1] LMS response. That's maybe a limitation for all color model in paper (because you always need to give a clear matrix in paper), and the model development may not dive that deeply into such tiny precision errors.

Yep, I totally understand this. I've seen this in other cases with color as well. In some papers, both the forward and reverse transform are provided, rounded to some digits, which is perfectly reasonable. People usually take the forward transform and invert it without a thought, and it can sometimes yield noise in higher precision values.

CSS had a similar issue with Oklab. Once they generated proper matrices for their version of D65 (per the author's recommendation), they still had noise in their 64-bit floats, most noticeable with achromatic colors. The matrices to get from LMS ** 1/3 to Oklab and back were provided as 32-bit forward and reverse matrices in the official Oklab write-up. People just assume that taking the forward LMS to Oklab matrix and inverting it gave the best results, and in 32-bit precision, it does just fine, but in 64-bit, you just get noise past 8 decimals. The reverse Oklab to LMS ** 1/3 matrix gave a perfect transform in 64-bit and 32-bit from Oklab [1, 0, 0] to LMS ** 1/3 of [1, 1, 1]. In this case, the 32-bit reverse transfrom gave better results in higher precision floating point. This was the case in CSS, and they were able to get cleaner results in 64-bit once they started creating the forward transfrom by inverting the reverse transform.

I think it will be interesting to use a pre-defined inverse matrix and use [1, 1, 1] to yield the optimal input D65 white point (like [0.9504, 1, 1.0889] in IPT), just like what you have done above. Then, the forward matrix can be created using numpy’s inv function from the inverse matrix to avoid precision issues in the pre-defined matrix. This idea is worth learning for me, and the affect for model performance should be very minor.

I think you can kind of adapt the current LMS matrix to work with a new white point. If [0.95047, 1.00000, 1.08883] is the preferred assumption for the Iab transform, you can take the inverse IPT matrix and probably adapt it as follows. This should give you clean LMS transforms with similar results.

>>> import numpy as np
>>> m = [
...     [1.8502, -1.1383, 0.2385],
...     [0.3668, 0.6439, -0.0107],
...     [0.0, 0.0, 1.0889]
... ]
>>> w = [0.95047, 1.00000, 1.08883]
>>> lms = np.linalg.solve(m, w)
>>> lms2xyz = np.multiply(m, lms)
>>> xyz2lms = np.linalg.inv(lms2xyz)
>>> lms2xyz.tolist()
[[1.8502622865264113, -1.138276954540003, 0.23848466801359167], [0.36681234823148184, 0.6438869639183941, -0.01069931214987602], [0.0, 0.0, 1.08883]]
>>> xyz2lms.tolist()
[[0.4002047620631536, 0.7074904187862833, -0.08070409428877665], [-0.22799040324166622, 1.150021385162843, 0.061236973086949954], [0.0, 0.0, 0.9184170164304805]]