[BUG]: Numerical precision issues with CIECAM02 and CIECAM16 for L_A = 0 and Y_b = 0

[mesvam]

Description

When L_A = 0 and/or Y_b = 0, CIECAM02 and CIECAM16 models give incorrect results

Code for Reproduction

import colour
import pandas as pd
sRGB = [
    (0, 0, 0), # black
    (0, 1, 0), # green
    (1, 1, 1)  # white
]
cam = colour.convert(sRGB, 'sRGB', 'CIECAM16',
    XYZ_w=colour.TVS_ILLUMINANTS["CIE 1931 2 Degree Standard Observer"]["D65"]/100,
    L_A=0, Y_b=20,
)
pd.DataFrame(cam, columns=['J', 'C', 'h', 's', 'Q', 'M', 'H', 'HC'])

C:\Users\yy51\AppData\Local\Programs\Miniconda3\Lib\site-packages\colour\algebra\common.py:488: RuntimeWarning: divide by zero encountered in power
  a_p = np.sign(a) * np.abs(a) ** p
C:\Users\yy51\AppData\Local\Programs\Miniconda3\Lib\site-packages\colour\algebra\common.py:488: RuntimeWarning: invalid value encountered in scalar multiply
  a_p = np.sign(a) * np.abs(a) ** p
     J             C    h    s    Q    M         H  HC
0  1.0  3.446651e-14  0.5  0.0  0.0  0.0  0.561821 NaN
1  1.0  3.446651e-14  0.5  0.0  0.0  0.0  0.561821 NaN
2  1.0  3.446651e-14  0.5  0.0  0.0  0.0  0.561821 NaN

All outputs are identical no matter the input. Looking a the output at L_A = 1e-30 suggests J, h, H should reach a limit, only C, s, Q, M should go to zero.

cam = colour.convert(sRGB, 'sRGB', 'CIECAM16',
    XYZ_w=colour.TVS_ILLUMINANTS["CIE 1931 2 Degree Standard Observer"]["D65"]/100,
    L_A=1e-30, Y_b=20,
)
pd.DataFrame(cam, columns=['J', 'C', 'h', 's', 'Q', 'M', 'H', 'HC'])

              J             C         h             s             Q             M         H  HC
0  8.734361e-08  1.018622e-17  0.500000  3.855345e-07  2.167138e-12  3.221165e-25  0.561821 NaN
1  7.875197e-01  3.343898e-10  0.395617  4.031120e-05  6.507311e-09  1.057433e-17  0.443577 NaN
2  1.000007e+00  1.354441e-11  0.583176  7.642647e-06  7.332842e-09  4.283117e-19  0.666321 NaN

Similarly for Y_b=0

cam = colour.convert(sRGB, 'sRGB', 'CIECAM16',
    XYZ_w=colour.TVS_ILLUMINANTS["CIE 1931 2 Degree Standard Observer"]["D65"]/100,
    L_A=64*0.2/np.pi, Y_b=0,
)
pd.DataFrame(cam, columns=['J', 'C', 'h', 's', 'Q', 'M', 'H', 'HC'])

     J    C         h    s    Q    M         H  HC
0  0.0  0.0  0.500000  0.0  0.0  0.0  0.561821 NaN
1  0.0  0.0  0.395067  0.0  0.0  0.0  0.443016 NaN
2  0.0  0.0  0.583308  0.0  0.0  0.0  0.666476 NaN

Testing Y_b=1e-30 suggests, J, h, H should reach a limit, s should go to zero, and C, Q, M should go to infinity.

cam = colour.convert(sRGB, 'sRGB', 'CIECAM16',
    XYZ_w=colour.TVS_ILLUMINANTS["CIE 1931 2 Degree Standard Observer"]["D65"]/100,
    L_A=64*0.2/np.pi, Y_b=1e-30,
)
pd.DataFrame(cam, columns=['J', 'C', 'h', 's', 'Q', 'M', 'H', 'HC'])

              J             C         h             s             Q             M         H  HC
0  9.187890e-19  1.147554e-17  0.500000  6.667251e-08  1.866260e-03  8.295944e-18  0.561821 NaN
1  8.352566e-01  3.861502e+05  0.395067  3.960840e-01  1.779403e+06  2.791572e+05  0.443016 NaN
2  1.000007e+00  1.070548e+04  0.583308  6.304728e-02  1.947000e+06  7.739247e+03  0.666476 NaN

Environment Information

===============================================================================
*                                                                             *
*   Interpreter :                                                             *
*       python : 3.11.5 | packaged by Anaconda, Inc. | (main, Sep 11 2023,    *
*   13:26:23) [MSC v.1916 64 bit (AMD64)]                                     *
*                                                                             *
*   colour-science.org :                                                      *
*       colour : 0.4.3                                                        *
*                                                                             *
*   Runtime :                                                                 *
*       imageio : 2.31.4                                                      *
*       matplotlib : 3.8.0                                                    *
*       networkx : 3.1                                                        *
*       numpy : 1.25.2                                                        *
*       pandas : 2.1.1                                                        *
*       scipy : 1.11.3                                                        *
*       tqdm : 4.65.0                                                         *
*                                                                             *
===============================================================================

defaultdict(dict,
            {'Interpreter': {'python': '3.11.5 | packaged by Anaconda, Inc. | (main, Sep 11 2023, 13:26:23) [MSC v.1916 64 bit (AMD64)]'},
             'colour-science.org': {'colour': '0.4.3'},
             'Runtime': {'imageio': '2.31.4',
              'matplotlib': '3.8.0',
              'networkx': '3.1',
              'numpy': '1.25.2',
              'pandas': '2.1.1',
              'scipy': '1.11.3',
              'tqdm': '4.65.0'}})

[KelSolaar]

Hello,

Those are warnings not errors and I do not really see any incorrect results here. CIECAM02 and CIECAM16 are models fitted on experimental data so it is not really a good idea to use an adapting field luminance or background luminous factor of 0. Those models are known to have some numerical precision issues.

Cheers,

Thomas

[mesvam]

I understand that the model may not be valid at these extremes, but even so, there's no mathematical reason for a discontinuity from 1e-30 to 0, especially since the output approaches a clear limit. It may still be of interest to provide a more sane output in the limit of total darkness.

[mesvam]

Also, @KelSolaar, would you happen to have any additional insight on my other question here: #1226? Even without this discontinuity issue, I'm a bit confused about the behavior of the model, which basically says that color goes away at low L_A, but subjectively I can still see color when surrounded by a dark environment. Am I misinterpreting the meaning of adapting field luminance L_A and background luminance Y_b?

[tjdcs]

Hi @mesvam

I’m on vacation right now but I will happily type up an explanation when I have access to my work stuff again later this week.

In short… the implementation of the model is correct. The model has some particularly nuanced and confusing parameters.

Cheers,

[tjdcs]

Hi @mesvam

In order to understand the model's instability for these inputs / outputs it's important to understand what they are modeling. If you aren't familiar with it I highly recommend Mark Fairchild's "Color Appearance Models" 3rd ed.

Let's get into the rabbit hole a bit...

First L_A and Y_b are two ways of controlling the same underlying parameter, what is the luminance of J = 100? Or what is often called "diffuse white". "J" is the correlate of "lightness", the same as L in CIELAB. A white piece of paper with no fluorescent whiteners has a "Lightness" of about 100. A grey patch that reflects about 20% of all light has a J value of around 50. This is due to the non-linear response of the eye.

So to transform from radiometric measurements, we need an anchor point. In the 70s and 80s when R.W.G. Hunt was first doing color appearance modeling at Kodak he used the grey surrounding card as the reference point. Y_b denotes the reflectance % of the surrounding card, typically 20%. L_A is the measured luminance of the surround. My desk, the luminance of 20% grey is is 6.8 nits this morning. If Y_b is 20 and L_A is 6.8 this implies a diffuse white point or adapting white point of 34 nits.

It's more common today to directly measure diffuse white with a spectrometer and use that measurement to derive the L_A value by dividing by 5. If that isn't an option, other experimenters will measure the illuminance using an illuminance meter, then divide by 5*pi.

If you are working on already normalized XYZ data already, then leaving L_A = 20 and Y_b = 20 is the correct input.

As you decrease L_a, that is the same as saying there is no light at all. In which case, there is no color.

Now we should also look at the "Colorfulness" and "Chroma" correlates, M and C respectively. First, I'll talk about the definitions of these.

Chroma is the colorfulness of a stimulus relative to the brightness of the scene / adaptation. Colorfulness is absolute in that it is not scaled by adaptation to the brightness of the scene. So again, as L_A goes towards zero colorfulness also goes towards 0. There is no light, or the amount of light is decreasing and so the stimulus must appear less colorful also. Chroma, being relative to adaptation, is more relevant in our visual system.

If you sit in a dark room, M might decrease to very low values because there is not a lot of chromatic appearance on an absolute scale, compared to a brighter room. But relative to sitting there in that room with your adapted visual system, there might still be enough chromatic stuff for you to see colors still. Chroma accounts for this. colorfulness (M) does not.

For appearance manipulation it's much more likely that you will have some success using J and C. Hopefully this also explains the usage of L_a better.

[mesvam]

@tjdcs Thanks for the detailed explanation!

Oh dear, I should not have gone down this rabbit hole, haha. I went and read up on a bunch of papers, and finally realized color theory is not at all an exact science and can get pretty messy.

As far as this issue is concerned though, from reading https://cielab.xyz/pdf/CIECAM02_and_Its_Recent_Developments.pdf, it seems like in the forward mode, L_A = 0 causes F_L = 0, which causes numerical issues when you take the power of 0^0.42 later. If you take 0^0.42 = 0, the model (or at least the forward mode) is otherwise perfectly stable.

[tjdcs]

Nah, the world needs more people thinking about this stuff. It affects things like CSS specification... w3c/csswg-drafts#9449

Just this morning a newly minted color scientist held their dissertation defense on some VAST improvements to CIECAM. (Dr. Luke Hellwig)

In this repository we are focused on accurate implementations of the various color science papers that we implement. Not necessarily the best or most useful implementations, or the "best" models (although we'd like to have a few of those around)

For various reasons... your other issue about J100, M0, h0 not returning the white point of the RGB color space... is probably actually not as "inaccurate" as you think. In color science we often talk about "Achromatic appearance" and actually, the particular tristimulus ratio(s) that appear achromatic are not always (rarely) the same as the white point.

It will greatly help you to start conceptualizing color science and color algorithms as "Color appearance modeling" and "color containers" or "color algorithm containers" or "color volumes." Often times the engineering or simplicity benefits of a particular color volume like RGBHSV have distinct advantages for particular applications. But they are not perceptual and they don't reflect color appearance. A mix of transforms and doing the right things in the right volume is the role of a color algorithms architect.

If you want to read more about achromatic appearance, Che Shen, Mark Fairchild, and Kevin Smet are good authors to get to know.

[tjdcs]

Sorry one last comment here. "Not an exact science" is one of my least favorite phrases to here about color science. Actually, it is an exact science modeling phenomenon that have a lot of noise or a lot of variance and a lot of independent variables. The methods of psychophysics are objective and based on centuries of scientific work. Work that also founded the fields of psychology, neuroscience, and many aspects of neurobiology. All imaging systems were first conceived relative to color science. It's just a complicated and small (number of researchers) field and filled with 1000x as many authors writing nonsense about color than there are researchers. Almost like the amount of noise in the field is similar to the underlying topic. :)

Small pet peeve. 😁

[mesvam]

Sorry one last comment here. "Not an exact science" is one of my least favorite phrases to here about color science.

Sorry, I didn't mean to suggest it's not scientific. I trust the method is scientific and the results are sound. And I understand human perception is variable and messy so the models will never be perfect.

But what really caught me off guard was realizing that the specifications aren't even mathematically consistent with itself, nevermind with each other!

Even something as basic as the D65 white point is in contention!

colour.xyY_to_XYZ(
    np.append(colour.CCS_ILLUMINANTS['CIE 1931 2 Degree Standard Observer']['D65'], 1)
)*100
array([  95.04559271,  100.        ,  108.90577508])

colour.TVS_ILLUMINANTS['CIE 1931 2 Degree Standard Observer']['D65']
array([  95.04,  100.  ,  108.88])

and if you try to compute it from the spectrum, it doesn't match either

# D65 from CCT
c2 = 6.62607015e-34 * 2.99792458e8 / 1.380649e-23 # second radiation constant
D65_CCT = 6500 * c2 / 0.014380
D65_xy = colour.temperature.CCT_to_xy_CIE_D(D65_CCT)
D65_XYZ = colour.xyY_to_XYZ(np.append(D65_xy, 1))
print(D65_XYZ*100)
[  95.01560357  100.          108.81851717]

# D65 from SPD
D65_sd = colour.sd_CIE_illuminant_D_series(D65_xy, M1_M2_rounding=True)
D65_XYZ = colour.sd_to_XYZ(D65_sd)
D65_XYZ = D65_XYZ / D65_XYZ[1] * 100
print(D65_XYZ)
[  95.04642202  100.          108.89689566]

I haven't deviated from spec, but I get 4 different values!

I mean yea it's not off by a lot, and I'm not asking for a precise definition of something as variable as daylight, but even if not accurate, at least make it consistent please!

When even basic linear transforms are not invertible, I gotta say it's a bit of a shock.

[mesvam]

particular tristimulus ratio(s) that appear achromatic are not always (rarely) the same as the white point

Just blew my mind again! How is the white point defined then, if it does not have to be achromatic? Any specific papers you recommend? I do not work in this field, so I'm not familiar with this, but just have a particular deep interest in this

[mesvam]

Correct me if I'm wrong, but even taking these issues aside, CIECAM model (or at least the official specification of it) is not exactly invertible right?

[mesvam]

In this repository we are focused on accurate implementations of the various color science papers that we implement.

Precisely why this package caught my attention. By far, the most complete and accurate implementation of the standard models. I can tell there is attention to detail.

[tjdcs]

I've never run into a use case where I can't invert CIECAM16, but I haven't checked these strange edge cases and also don't use them to model un-related colors. If these discontinuities are the only issue, patching them with a simple check would probably be acceptable. @KelSolaar thoughts? (PR's welcome 🤗)

[tjdcs]

I mean yea it's not off by a lot, and I'm not asking for a precise definition of something as variable as daylight, but even if not accurate, at least make it consistent please!

Numerical consistency is a big question through out the 90s and 2000s in CIE (and somewhat today too)

but technically the CIE only defined D65 based on a particular data table at 5nm. It is the authoritative D65... But daylight is insanely variable. And if you use that data-table value, use recommended methods to interpolate to 1nm, and compute CCT... D65 actually has a CCT of like 6504K (actually, depending on the exact choices of numerical methods you can get 6503 or 6502) but real daylight is insanely variable, and sensitivity is quite low to these changes so the small numerical differences don't practically matter.

But they do always drive chemists and physicist crazy. We have a noisy environment, and we can't even predict basic phenomenon like simultaneous contrast or neon color spreading yet. Or the differences between people, which can be as much as 12dE in extreme cases. So small issues like this are not a bother.

It just means you can't use ==

[tjdcs]

There is a developers / advanced users group if you want to chat more easily. Send me a message to tucker at tjdcs dot dev and I'll request an invite for you.

[mesvam]

But they do always drive chemists and physicist crazy

You got me :P, my pet peeve

[mesvam]

Nope, I was wrong. J = 100 (A/A_w)^(c*z) seems to be the problem here. A and A_w should go to 0, but doesn't due to numerical precision, so instead you get A/A_w = 1. But I think there is a solution for the limit

[mesvam]

I've analytically determined that $$\lim_{L_A \to 0} \frac{A}{A_w} = \frac{2{R^\prime}^{0.42}+ {G^\prime}^{0.42} + {B^\prime}^{0.42}/10}{2{R^\prime_{w}}^{0.42}+ {G^\prime_{w}}^{0.42} + {B^\prime_{w}}^{0.42}/10}$$

So if we were to fix this edge case, we'd manually set A = 0, A_w = 0, and $$J = 100 \left( \frac{2{R^\prime}^{0.42}+ {G^\prime}^{0.42} + {B^\prime}^{0.42}/10}{2{R^\prime_{w}}^{0.42}+ {G^\prime_{w}}^{0.42} + {B^\prime_{w}}^{0.42}/10} \right)^{cz}$$