Understanding LMS-based Color Blindness Simulations¶
Diving deeper in the theory of Viénot, Brettel & Mollon's method to understand what's behind the scene.
- toc: true
- badges: true
- comments: true
- categories: [simulation]
- permalink: /understanding-cvd-simulation/
- hide: false
- sticky_rank: 2
# hide
%load_ext autoreload
%autoreload 2
# hide
%pip install -q pandas Geometry3D colour-science
%pip install -q --upgrade plotly daltonlens
Note: you may need to restart the kernel to use updated packages. Note: you may need to restart the kernel to use updated packages.
# hide_input
import plotly.io as pio
if pio.renderers.default != 'colab':
pio.renderers.default = 'notebook_connected'
from IPython.display import display, HTML
js = '<script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" integrity="sha512-c3Nl8+7g4LMSTdrm621y7kf9v3SDPnhxLNhcjFJbKECVnmZHTdo+IRO05sNLTH/D3vA6u1X32ehoLC7WFVdheg==" crossorigin="anonymous"></script>'
display(HTML(js))
Introduction¶
This post is a followup to our Review of Open Source Color Blindness Simulations and assumes that you are familiar with color vision deficiencies (CVD) simulation. The goal of this post is to dive into how these methods actually work. To visualize and play with the models interactively this page is actually generated from a Jupyter notebook (using fastpages) and can be run in e.g. Google colab to play with it.
The interactive 3D plots make it hard to read on mobile, a large screen is recommended.
{% cite brettel_computerized_1997 %}, {% cite vienot_digital_1999 %} and {% cite machado_physiologically_based_2009 %} all rely on the same low-level principles to simulate CVD. For all of them the final code will basically just apply a sRGB/gamma function, one or two 3x3 matrices, undo the gamma, and call it a day. But we're going to analyze that pipeline and try to understand how to compute those matrices and get a sense of what they do.
The Machado 2009 method is actually a bit more involved as it uses a second stage of modeling (opponent-color theory), so we will not cover it here and only focus on Viénot 1999 and Brettel 1997. Both focus on full dichromacy, so we'll start with that. It's worth noting that Machado 2009 has similar results for dichromacy anyway, and the extension to anomalous trichromacy will be discussed in section What about anomalous trichromacy?.
Let's first get an overview of the process. Protanopes, deuteranopes and tritanopes respectively lack the L, M, and S cone cells. Thus a natural approach to simulate dichromacy is to go to the LMS color-space and somehow remove the component coming from the missing cone cells. The good news is that going from an RGB color space to LMS can be done with a linear transformation (matrix multiplication). But as usual there are many ways to do it and it'll deserve its own section.
One extra complexity is that images are typically encoded in sRGB to get properly displayed on digital monitors. So our first step will be to go from sRGB to linear RGB, and then from there to LMS.
Once in the LMS color space, colors that only differ by the coordinate of the axis corresponding to the missing cone cannot be differentiated by a dichromat. This creates so-called confusion lines, which are, by construction, parallel to the missing axis in the LMS color space. So according to this model a protanope cannot distinguish colors that only differ by their L component, a deuteranope those that only differ by their M component, and tritanopes those that only differ by their S component.
To illustrate this let's build a plot with a few confusion lines in the LMS space, and their counterpart on the linear RGB space. The code to generate the confusion lines is in the Appendix. We're going to just draw 5 precomputed lines, starting with protanopia (missing L cones). To have a reference we'll also add the gamut (range) of colors that can be produced by a digital monitor. It's a cube in the RGB color space, and a parallelepiped in the LMS space. The 8 vertices are black (K = RGB[0,0,0]), white (W = RGB[1, 1, 1]), red (R = RGB[1, 0, 0]), green (G = RGB[0,1,0]), blue (B = RGB[0,0,1]), yellow (Y = RGB[1,1,0]), magenta (M = RGB[1,0,1]) and cyan (C = RGB[0,1,1]).
The plots are in 3D and interactive, so feel free to play with the point of view!
# collapse-hide
# This cell defines all our imports and a few utilities
import itertools
import numpy as np
import plotly.express as px
import plotly.graph_objs as go
from plotly.subplots import make_subplots
import pandas as pd
import Geometry3D as geo3d
import daltonlens.convert as convert
import daltonlens.generate as generate
import daltonlens.simulate as simulate
import daltonlens.geometry as geometry
np.set_printoptions(suppress=True)
def printMatrix(name, m):
mat_str = '\n '.join(np.array_repr(m).split('\n'))
print (f'{name} = \n {mat_str}')
# Utility for plotly imshow
hide_image_axes = dict(yaxis_visible=False, yaxis_showticklabels=False, xaxis_visible=False, xaxis_showticklabels=False, margin=dict(l=0, r=0, b=0, t=0))
# collapse-hide
def show_confusion_lines(title, lines):
rgb_cube = pd.DataFrame.from_records([
('K', 'black', 0, 0, 0, 0. , 0. , 0. , '#000000'),
('B', 'blue', 0, 0, 1, 0.03596577, 0.03620635, 0.01528089, '#0000ff'),
('M', 'magenta', 1, 0, 1, 0.21482533, 0.07001029, 0.01559176, '#ff00ff'),
('R', 'red', 1, 0, 0, 0.17885956, 0.03380394, 0.00031087, '#ff0000'),
('G', 'green', 0, 1, 0, 0.43997117, 0.27515242, 0.00191661, '#00ff00'),
('C', 'cyan', 0, 1, 1, 0.47593694, 0.31135877, 0.0171975 , '#00ffff'),
('W', 'white', 1, 1, 1, 0.6547965 , 0.34516271, 0.01750837, '#ffffff'),
('Y', 'yellow', 1, 1, 0, 0.61883073, 0.30895636, 0.00222748, '#ffff00')],
columns = ['short_name', 'name', 'R', 'G', 'B', 'L', 'M', 'S', 'sRGB_hex'])
rgb_cube.index = rgb_cube['short_name']
rgb_cube.index.name = 'key'
# display(HTML("<h2>RGB Cube</h2>"), rgb_cube)
fig = make_subplots(1,2,
subplot_titles=('RGB Color Space', 'LMS Color Space'),
specs=[[{"type": "scene"}, {"type": "scene"}]])
fig.add_scatter3d(x=rgb_cube['R'], y=rgb_cube['G'], z=rgb_cube['B'],
text=rgb_cube.index, mode="markers+text",
showlegend=False,
marker=dict(color=rgb_cube['sRGB_hex']),
row=1, col=1)
fig.add_scatter3d(x=rgb_cube['L'], y=rgb_cube['M'], z=rgb_cube['S'],
text=rgb_cube.index, mode="markers+text",
showlegend=False,
marker=dict(color=rgb_cube['sRGB_hex']),
row=1, col=2)
fig.update_layout(title=title)
fig.update_layout(dragmode='orbit',
scene = dict(
xaxis_title='R',
yaxis_title='G',
zaxis_title='B'))
fig.update_layout(scene2 = dict(
xaxis_title='L',
yaxis_title='M',
zaxis_title='S'))
# Add the parallelepiped lines
df = rgb_cube.loc[['K','B','M','R', 'K', 'G', 'Y', 'W', 'C', 'G', 'Y', 'R', 'M', 'W', 'C', 'B']]
fig.add_scatter3d(x=df['R'], y=df['G'], z=df['B'], mode='lines', showlegend=False, row=1, col=1)
fig.add_scatter3d(x=df['L'], y=df['M'], z=df['S'], mode='lines', showlegend=False, row=1, col=2)
for line in lines:
fig.add_scatter3d(x=line['R'], y=line['G'], z=line['B'], showlegend=False,
marker=dict(color=line['sRGB_hex']),
row=1, col=1)
fig.add_scatter3d(x=line['L'], y=line['M'], z=line['S'], showlegend=False,
marker=dict(color=line['sRGB_hex']),
row=1, col=2)
fig.show()
# Confusion lines for Deficiency.PROTAN
protan_confusion_lines = [
pd.DataFrame.from_records([
(0.2299, 0.0690, 0.0034, 1.0000, 0.0987, 0.1964, '#fe587a'),
(0.2093, 0.0690, 0.0034, 0.8333, 0.1198, 0.1971, '#eb617a'),
(0.1887, 0.0690, 0.0034, 0.6667, 0.1409, 0.1979, '#d5687a'),
(0.1681, 0.0690, 0.0034, 0.5000, 0.1620, 0.1986, '#bb6f7b'),
(0.1475, 0.0690, 0.0034, 0.3334, 0.1831, 0.1994, '#9c767b'),
(0.1269, 0.0690, 0.0034, 0.1667, 0.2042, 0.2002, '#717c7b'),
(0.1063, 0.0690, 0.0034, 0.0001, 0.2253, 0.2009, '#00827b')],
columns = ['L', 'M', 'S', 'R', 'G', 'B', 'sRGB_hex']),
pd.DataFrame.from_records([
(0.2658, 0.1724, 0.0084, 0.0000, 0.5633, 0.5023, '#00c5bb'),
(0.2864, 0.1724, 0.0084, 0.1667, 0.5422, 0.5015, '#71c2bb'),
(0.3070, 0.1724, 0.0084, 0.3333, 0.5211, 0.5008, '#9cbfbb'),
(0.3276, 0.1724, 0.0084, 0.5000, 0.5000, 0.5000, '#bbbbbb'),
(0.3482, 0.1724, 0.0084, 0.6666, 0.4789, 0.4992, '#d5b7bb'),
(0.3688, 0.1724, 0.0084, 0.8333, 0.4578, 0.4985, '#ebb4bb'),
(0.3894, 0.1724, 0.0084, 0.9999, 0.4367, 0.4977, '#feb0bb')],
columns = ['L', 'M', 'S', 'R', 'G', 'B', 'sRGB_hex']),
pd.DataFrame.from_records([
(0.2770, 0.1055, 0.0115, 1.0000, 0.1550, 0.7466, '#fe6de0'),
(0.2565, 0.1055, 0.0115, 0.8333, 0.1761, 0.7474, '#eb74e0'),
(0.2359, 0.1055, 0.0115, 0.6667, 0.1972, 0.7481, '#d57ae0'),
(0.2153, 0.1055, 0.0115, 0.5000, 0.2183, 0.7489, '#bb80e0'),
(0.1947, 0.1055, 0.0115, 0.3334, 0.2394, 0.7496, '#9c86e0'),
(0.1741, 0.1055, 0.0115, 0.1667, 0.2605, 0.7504, '#718be0'),
(0.1535, 0.1055, 0.0115, 0.0001, 0.2816, 0.7511, '#0090e0')],
columns = ['L', 'M', 'S', 'R', 'G', 'B', 'sRGB_hex']),
pd.DataFrame.from_records([
(0.5017, 0.2393, 0.0053, 1.0000, 0.7183, 0.2489, '#fedc88'),
(0.4811, 0.2393, 0.0053, 0.8333, 0.7394, 0.2496, '#ebdf88'),
(0.4605, 0.2393, 0.0053, 0.6667, 0.7605, 0.2504, '#d5e289'),
(0.4399, 0.2393, 0.0053, 0.5000, 0.7816, 0.2511, '#bbe489'),
(0.4193, 0.2393, 0.0053, 0.3334, 0.8027, 0.2519, '#9ce789'),
(0.3987, 0.2393, 0.0053, 0.1667, 0.8239, 0.2526, '#71ea89'),
(0.3781, 0.2393, 0.0053, 0.0001, 0.8450, 0.2534, '#00ec89')],
columns = ['L', 'M', 'S', 'R', 'G', 'B', 'sRGB_hex']),
pd.DataFrame.from_records([
(0.4253, 0.2758, 0.0135, 0.0000, 0.9013, 0.8036, '#00f3e7'),
(0.4459, 0.2758, 0.0135, 0.1667, 0.8802, 0.8029, '#71f1e7'),
(0.4665, 0.2758, 0.0135, 0.3333, 0.8591, 0.8021, '#9ceee7'),
(0.4871, 0.2758, 0.0135, 0.5000, 0.8380, 0.8014, '#bbebe7'),
(0.5077, 0.2758, 0.0135, 0.6666, 0.8169, 0.8006, '#d5e9e7'),
(0.5283, 0.2758, 0.0135, 0.8333, 0.7958, 0.7998, '#ebe6e7'),
(0.5488, 0.2758, 0.0135, 0.9999, 0.7747, 0.7991, '#fee3e6')],
columns = ['L', 'M', 'S', 'R', 'G', 'B', 'sRGB_hex']),
]
show_confusion_lines('Examples of protan confusion lines', protan_confusion_lines)
The colors of some of the 5 lines might be easier to discriminate for some protan, but overall with this point size they should be pretty difficult to differentiate on a monitor screen.
See how the confusion lines are parallel to the L axis in the LMS color space, but not aligned with any axis in the RGB space. As expected for protanopes, the colors along the confusion lines vary mostly by their red and green components and are pretty much orthogonal to the blue axis.
Let's also visualize some lines for deuteranopes (parallel to the M axis) and tritanopes (parallel to the S axis).
# collapse-hide
# Confusion lines for Deficiency.DEUTAN
deutan_confusion_lines = [
pd.DataFrame.from_records([
(0.1310, 0.0843, 0.0034, 0.0000, 0.2828, 0.1937, '#009079'),
(0.1310, 0.0756, 0.0034, 0.1139, 0.2357, 0.1973, '#5e857a'),
(0.1310, 0.0668, 0.0034, 0.2277, 0.1885, 0.2009, '#83787b'),
(0.1310, 0.0581, 0.0034, 0.3416, 0.1414, 0.2045, '#9d697c'),
(0.1310, 0.0494, 0.0034, 0.4554, 0.0943, 0.2081, '#b3567d'),
(0.1310, 0.0407, 0.0034, 0.5693, 0.0471, 0.2117, '#c63d7e'),
(0.1310, 0.0319, 0.0034, 0.6831, 0.0000, 0.2153, '#d7007f')],
columns = ['L', 'M', 'S', 'R', 'G', 'B', 'sRGB_hex']),
pd.DataFrame.from_records([
(0.3276, 0.1341, 0.0084, 1.0000, 0.2930, 0.5158, '#fe93be'),
(0.3276, 0.1468, 0.0084, 0.8333, 0.3620, 0.5105, '#eba2bd'),
(0.3276, 0.1596, 0.0084, 0.6667, 0.4310, 0.5053, '#d5afbc'),
(0.3276, 0.1724, 0.0084, 0.5000, 0.5000, 0.5000, '#bbbbbb'),
(0.3276, 0.1852, 0.0084, 0.3334, 0.5690, 0.4947, '#9cc6ba'),
(0.3276, 0.1979, 0.0084, 0.1667, 0.6380, 0.4895, '#71d1b9'),
(0.3276, 0.2107, 0.0084, 0.0001, 0.7069, 0.4842, '#00dab8')],
columns = ['L', 'M', 'S', 'R', 'G', 'B', 'sRGB_hex']),
pd.DataFrame.from_records([
(0.1844, 0.1247, 0.0115, 0.0000, 0.3535, 0.7421, '#00a0df'),
(0.1844, 0.1138, 0.0115, 0.1423, 0.2946, 0.7466, '#6993e0'),
(0.1844, 0.1029, 0.0115, 0.2846, 0.2357, 0.7511, '#9185e0'),
(0.1844, 0.0920, 0.0115, 0.4270, 0.1768, 0.7556, '#ae74e1'),
(0.1844, 0.0811, 0.0115, 0.5693, 0.1178, 0.7601, '#c660e1'),
(0.1844, 0.0702, 0.0115, 0.7116, 0.0589, 0.7646, '#db44e2'),
(0.1844, 0.0593, 0.0115, 0.8539, 0.0000, 0.7691, '#ed00e3')],
columns = ['L', 'M', 'S', 'R', 'G', 'B', 'sRGB_hex']),
pd.DataFrame.from_records([
(0.4708, 0.2201, 0.0053, 1.0000, 0.6465, 0.2579, '#fed28a'),
(0.4708, 0.2310, 0.0053, 0.8577, 0.7054, 0.2534, '#eeda89'),
(0.4708, 0.2419, 0.0053, 0.7154, 0.7643, 0.2489, '#dbe288'),
(0.4708, 0.2528, 0.0053, 0.5730, 0.8232, 0.2444, '#c7ea87'),
(0.4708, 0.2637, 0.0053, 0.4307, 0.8822, 0.2399, '#aff186'),
(0.4708, 0.2746, 0.0053, 0.2884, 0.9411, 0.2354, '#92f885'),
(0.4708, 0.2855, 0.0053, 0.1461, 1.0000, 0.2309, '#6afe84')],
columns = ['L', 'M', 'S', 'R', 'G', 'B', 'sRGB_hex']),
pd.DataFrame.from_records([
(0.5241, 0.3128, 0.0135, 0.3168, 1.0000, 0.7847, '#98fee5'),
(0.5241, 0.3041, 0.0135, 0.4307, 0.9529, 0.7883, '#aff9e5'),
(0.5241, 0.2954, 0.0135, 0.5445, 0.9057, 0.7919, '#c2f4e6'),
(0.5241, 0.2867, 0.0135, 0.6584, 0.8586, 0.7955, '#d4eee6'),
(0.5241, 0.2780, 0.0135, 0.7722, 0.8115, 0.7991, '#e3e8e7'),
(0.5241, 0.2692, 0.0135, 0.8861, 0.7644, 0.8027, '#f1e2e7'),
(0.5241, 0.2605, 0.0135, 1.0000, 0.7172, 0.8063, '#fedce7')],
columns = ['L', 'M', 'S', 'R', 'G', 'B', 'sRGB_hex']),
]
show_confusion_lines('Examples of deutan confusion lines', deutan_confusion_lines)
# collapse-hide
# Confusion lines for Deficiency.TRITAN
tritan_confusion_lines = [
pd.DataFrame.from_records([
(0.1310, 0.0690, 0.0005, 0.1663, 0.2328, 0.0000, '#718400'),
(0.1310, 0.0690, 0.0029, 0.1944, 0.2055, 0.1667, '#797d71'),
(0.1310, 0.0690, 0.0053, 0.2224, 0.1782, 0.3333, '#81759c'),
(0.1310, 0.0690, 0.0077, 0.2505, 0.1509, 0.5000, '#896cbb'),
(0.1310, 0.0690, 0.0101, 0.2785, 0.1236, 0.6666, '#8f62d5'),
(0.1310, 0.0690, 0.0125, 0.3066, 0.0963, 0.8333, '#9657eb'),
(0.1310, 0.0690, 0.0149, 0.3346, 0.0689, 0.9999, '#9c4afe')],
columns = ['L', 'M', 'S', 'R', 'G', 'B', 'sRGB_hex']),
pd.DataFrame.from_records([
(0.3276, 0.1724, 0.0156, 0.5842, 0.4181, 1.0000, '#c9adfe'),
(0.3276, 0.1724, 0.0132, 0.5561, 0.4454, 0.8333, '#c4b2eb'),
(0.3276, 0.1724, 0.0108, 0.5281, 0.4727, 0.6667, '#c0b6d5'),
(0.3276, 0.1724, 0.0084, 0.5000, 0.5000, 0.5000, '#bbbbbb'),
(0.3276, 0.1724, 0.0060, 0.4720, 0.5273, 0.3334, '#b6c09c'),
(0.3276, 0.1724, 0.0036, 0.4439, 0.5546, 0.1667, '#b1c471'),
(0.3276, 0.1724, 0.0012, 0.4159, 0.5819, 0.0001, '#acc800')],
columns = ['L', 'M', 'S', 'R', 'G', 'B', 'sRGB_hex']),
pd.DataFrame.from_records([
(0.1844, 0.1055, 0.0007, 0.1238, 0.3729, 0.0000, '#62a400'),
(0.1844, 0.1055, 0.0031, 0.1518, 0.3456, 0.1667, '#6c9e71'),
(0.1844, 0.1055, 0.0055, 0.1799, 0.3183, 0.3333, '#75989c'),
(0.1844, 0.1055, 0.0079, 0.2079, 0.2910, 0.5000, '#7d92bb'),
(0.1844, 0.1055, 0.0103, 0.2360, 0.2637, 0.6666, '#858cd5'),
(0.1844, 0.1055, 0.0127, 0.2640, 0.2364, 0.8333, '#8c85eb'),
(0.1844, 0.1055, 0.0151, 0.2921, 0.2091, 0.9999, '#937efe')],
columns = ['L', 'M', 'S', 'R', 'G', 'B', 'sRGB_hex']),
pd.DataFrame.from_records([
(0.4708, 0.2393, 0.0161, 0.8762, 0.6271, 1.0000, '#f0cffe'),
(0.4708, 0.2393, 0.0137, 0.8482, 0.6544, 0.8333, '#edd3eb'),
(0.4708, 0.2393, 0.0113, 0.8201, 0.6817, 0.6667, '#e9d7d5'),
(0.4708, 0.2393, 0.0089, 0.7921, 0.7090, 0.5000, '#e6dbbb'),
(0.4708, 0.2393, 0.0065, 0.7640, 0.7363, 0.3334, '#e2de9c'),
(0.4708, 0.2393, 0.0041, 0.7360, 0.7636, 0.1667, '#dee271'),
(0.4708, 0.2393, 0.0017, 0.7079, 0.7909, 0.0001, '#dae500')],
columns = ['L', 'M', 'S', 'R', 'G', 'B', 'sRGB_hex']),
pd.DataFrame.from_records([
(0.5241, 0.2758, 0.0019, 0.6654, 0.9311, 0.0000, '#d4f700'),
(0.5241, 0.2758, 0.0043, 0.6934, 0.9038, 0.1667, '#d8f371'),
(0.5241, 0.2758, 0.0067, 0.7215, 0.8765, 0.3333, '#dcf09c'),
(0.5241, 0.2758, 0.0091, 0.7495, 0.8492, 0.5000, '#e0edbb'),
(0.5241, 0.2758, 0.0115, 0.7776, 0.8218, 0.6666, '#e4e9d5'),
(0.5241, 0.2758, 0.0139, 0.8056, 0.7945, 0.8333, '#e7e6eb'),
(0.5241, 0.2758, 0.0163, 0.8337, 0.7672, 0.9999, '#ebe2fe')],
columns = ['L', 'M', 'S', 'R', 'G', 'B', 'sRGB_hex']),
]
show_confusion_lines('Examples of tritan confusion lines', tritan_confusion_lines)
The deutan confusion lines are somewhat similar to the protan confusion lines (this is why both are called red-green deficiency), but the tritan ones are mostly along the blue axis.
At this stage if you are affected by a color vision deficiency you might already have a more precise idea of what kind of dichromacy better models your color vision. For me the colors on the protan confusion lines look pretty similar at such a small point size, but I can better distinguish the deutan lines and very well the tritan ones, confirming that I am a protan.
Dichromats confusion lines are sometimes shown in the CIE XYZ color space (e.g. on color-blindness.com), usually drawing the lines on a so-called chromaticity diagram. These diagram are a 2d projection of the CIE XYZ color space for a constant maximal luminance. I always got confused about those because I could clearly see differences along these lines, especially when they transition from green to red. The reason is that the confusion lines are actually 3d in CIE XYZ and the luminance needs to get adjusted as the hue/chroma changes along the line to actually get confused by a dichromat.
So now that we can compute these confusion lines for each kind of dichromacy, can we simulate how a person with a CVD sees? Since it's related to perception it's impossible to know how another person "feels" the colors, but to understand the loss of color sensitivity we can transform images so that all the colors on each confusion line projects to a single one, essentially making the color space 2D. If the models are correct then the transformed image should still appear identical to a dichromat, and a trichromat will get a sense of the color differences that get lost for the dichromat.
The remaining question is what color should we pick to represent each confusion line? That will be the subject of the Simulating dichromacy in LMS section.
Finally once the LMS colors are projected along the confusion lines, we can transform them back to the linearRGB and then sRGB color spaces to get the final simulated image. The overall process is summarized in the diagram below, and we're going to dive into each part in the next sections.
From sRGB to LMS¶
In this section we're going to review the steps to go to that nice LMS color space. Several models have been proposed in the literature, so there is not just one way to do it.
From sRGB to linearRGB¶
Most digital images are encoded with 3 values per pixel: red, green and blue. But what is often neglected is that this "RGB" usually corresponds to the sRGB standard that precisely defines how these values should be interpreted. sRGB is not a linear space because it tries to compensate the non-linear gamma factor applied by monitors, so that the image captured by a digital camera and then rendered on a monitor will look like the original scene.
We'll call the RGB color space after removing that non-linearity the linearRGB color space. The transformation roughly corresponds to applying a gamma exponent of 2.4 on each channel, with a special linear slope for small values. So unlike all the other color transformations of that notebook it is non-linear and cannot be expressed as a matrix.
The good news is that sRGB is a well-defined standard now, so we just need to implement the function given by the standard. For more details see the wikipedia page.
Since both Brettel 1997 and Viénot 1999 were published before the sRGB standard became mainstream, they use a different model in their paper. Brettel 1997 calibrated a CRT monitor with a spectroradiometer to estime the spectral distribution of the primaries, while Viénot used a constant gamma of 2.2 to convert to linearRGB. Fortunately once in that space it uses the same BT.709 chromaticities as the sRGB standard, so it does not need to be adapted. However it's worth noting that even if the Viénot 1999 paper explicitly mentions the required non-linear gamma transform for their CRT monitor, the most widely replicated matlab code from {% cite fidaner_analysis_2005 %} does not take it into account and directly treat the input image as linear RGB. As a consequence this is also the case of most of the source code listed on websites like daltonize.org.
Does it matter? Yes!¶
As you can see in the image below forgetting the sRGB -> linearRGB transform makes a big difference and it can invalidate perceptual experiments. The first row has 4 sRGB colors, and the second row shows the same colors wrongly interpreted as linear RGB.
# collapse-hide
im_sRGB = np.array([[[127,0,0], [64,128,0], [0, 0, 192], [64, 64, 192]]], dtype=float) / 255.0
px.imshow(np.vstack([im_sRGB, # First row is the sRGB content
convert.linearRGB_from_sRGB(im_sRGB)]), height=200).update_layout(hide_image_axes)
And the impact is also significant on the final simulation, here is a comparison of Viénot 1999 implemented with and without the proper sRGB decoding / encoding. Note how a bunch of colors appear too dark in the simulation without sRGB.
| Original | Viénot 1999 | Viénot 1999 without sRGB |
|---|---|---|
 |
 |
 |
From linearRGB to CIE XYZ¶
So far we only mentioned the LMS color space, saying that we'd just go from linearRGB to LMS. But to understand the literature around that transform it is necessary to add an intermediate step which is the CIE XYZ color space. The reason is historical, this color space was standardized in 1931 and it has been used as the main space for color science research since. The main reason is that it can (almost) represent the full range of colors that can be perceived by an average human observer with just 3 values, is device-independent, and has a well-defined relationship w.r.t the input light spectrum.
The [sRGB wikipedia page](https://en.wikipedia.org/wiki/SRGB#The_reverse_transformation_(sRGB_to_CIE_XYZ) gives the XYZ_from_linearRGB matrix defined by the standard itself. But we need to be careful, as this assumes the CIE XYZ 1931 color space. Some LMS_from_XYZ transforms might assume a different XYZ color space. For example, as noted in Viénot 1999, the LMS from XYZ transform given by {% cite smith_spectral_1975 %} is defined on top of the Judd-Vos corrected XYZ space. As another example the {% cite stockman_spectral_2000 %} one assumes the XYZ space from the Stiles and Burch 10-deg CMFs adjusted to 2-deg.
This section was initially wrong and got updated in December 2021. Before I was assuming that Viénot 1999 was using an outdated model for CRT monitors and I ignored the Judd-Vos correction. But only the non-linear gamma was outdated, the XYZ from linearRGB was using the same BT.709 standard, so that conversion is still valid. Many thanks to Paul Maximov who pointed this out and helped me understand that part.
It's easy to get this wrong, but eventually Viénot proposed a pretty good approximation with a single XYZJuddVos_from_linearRGB matrix, which can then get combined with the LMS_from_RGB by Smith & Pokorny. They evaluated their approximation properly, so refer to the paper for the full details.
# collapse-hide
printMatrix('sRGB XYZJuddVos_from_linearRGB', convert.LMSModel_sRGB_SmithPokorny75().XYZ_from_linearRGB*100)
sRGB XYZJuddVos_from_linearRGB =
array([[40.9568 , 35.5041 , 17.9167 ],
[21.3389 , 70.6743 , 7.9868 ],
[ 1.86297, 11.462 , 91.2367 ]])
From CIE XYZ to LMS¶
That conversion is less clear as many more transforms have been proposed, with different motivations, and even different source XYZ color spaces as we've just seen.
The reference conversion for CVD simulation is probably still {% cite smith_spectral_1975 %}, which is the one used by {% cite vienot_digital_1999 %}. The earlier seminal work of {% cite meyer_color_defective_1988 %} used another transform from {% cite estevez1981fundamental %} but it already mentioned {% cite smith_spectral_1975 %} as being a good candidate. Another popular reference model in the literature is the Hunt-Pointer-Esterez one {% cite hunt_reproduction_2004 %}. It is used a lot for chromatic adaptation (predicting how a perceived color will change with the illumination conditions), but less for CVD simulation experiments.
{% cite stockman_spectral_2000 %} made a thorough comparison of the models, and they also propose a new matrix (for a 10 degree standard observer). But overall they are in good agreement with {% cite smith_spectral_1975 %} once adjusted for a 2 degree observer, with larger differences compared to Hunt-Pointer-Esterez. All 3 are probably reasonable options, but since the CVD simulation evaluations were most often made and evaluated with Smith & Pokorny it's probably still the most solid choice, despite its age. It is also interesting to note that both Pokorny and Viénot continued to promote that transform in more recent papers {% cite vienot_colorimetry_2013 %} {% cite smith_color_2003 %}, and other recent papers in the color science research literature do the same {% cite machado_physiologically_based_2009 %} {% cite simon_liedtke_using_2016 %}.
Most open source code also uses the Smith and Pokorny since they did not modify the original code, but some did. Ixora.io uses the Hunt-Pointer-Estevez matrix and daltonize.py uses the sharpened transformation matrix of CIECAM02. These two matrices are on the LMS Wikipedia page and look attractive as they are more recent. But at least the "sharpened" matrices variants were not designed to match the physiological cone responses but to perform hue-preserving chromatic adaptation (see the discussion p181 of {% cite fairchild2013color %}). The wikipedia page also says A "spectually sharpened" matrix is believed to improve chromatic adaptation especially for blue colors, but does not work as a real cone-describing LMS space for later human vision processing. So I would argue that we should keep using the Smith & Pokorny transform until a proper validation of newer matrices gets done in the context of CVD simulation.
The XYZ<->LMS matrices of the different models can look very different, but the main reason is the normalization factor. There is no unique "scale" for each axis, so each model picked a convention. For example Smith & Pokorny made it so L+S=Y and S/Y=1 at the 400nm wavelengths. For Hunt-Pointer-Esterez they normalize it so that L=M=S for the equal-energy illuminant E (X=Y=Z=1/3), which means that the sum of each rows equals to 1. One reason for this scale variance is that the human vision system can actually change it independently for each axis depending on the illumination conditions. This process is called chromatic adaptation and explains why we still perceive a white sheet of paper as white under an incandescent indoor lighting that makes it actually yellow-ish. It's a bit similar to the white balance correction of digital cameras. These normalizations do not really matter for CVD simulation if we only do linear operations in it, the scale factor will cancel out when we go back to RGB.
Below are a couple examples of those LMS_from_XYZ matrices.
# collapse-hide
printMatrix('Smith & Pokorny (used by Brettel & Viénot, recommended) LMS_from_XYZJuddVos', convert.LMSModel_sRGB_SmithPokorny75().LMS_from_XYZ); print()
printMatrix('Hunt-Pointer-Estevez LMS_from_XYZ', convert.LMSModel_sRGB_HuntPointerEstevez().LMS_from_XYZ); print()
Smith & Pokorny (used by Brettel & Viénot, recommended) LMS_from_XYZJuddVos =
array([[ 0.15514, 0.54312, -0.03286],
[-0.15514, 0.45684, 0.03286],
[ 0. , 0. , 0.01608]])
Hunt-Pointer-Estevez LMS_from_XYZ =
array([[ 0.4002, 0.7076, -0.0808],
[-0.2263, 1.1653, 0.0457],
[ 0. , 0. , 0.9182]])
From linearRGB to LMS¶
Since both XYZ_from_linearRGB and LMS_from_XYZ are linear transforms, we can just combine the two matrices to get a final LMS_from_linearRGB = LMS_from_XYZ . XYZ_from_linearRGB transform. Here are some of the models implemented in DaltonLens-Python, scaled by 100 for an easier comparison with existing open source code:
# collapse-hide
printMatrix('sRGB Viénot (recommended)\nLMS_from_linearRGB', convert.LMSModel_sRGB_SmithPokorny75().LMS_from_linearRGB*100.0); print()
printMatrix('Hunt-Pointer-Esterez\nLMS_from_linearRGB', convert.LMSModel_sRGB_HuntPointerEstevez().LMS_from_linearRGB*100.0); print()
sRGB Viénot (recommended)
LMS_from_linearRGB =
array([[17.88240413, 43.51609057, 4.11934969],
[ 3.45564232, 27.15538246, 3.86713084],
[ 0.02995656, 0.18430896, 1.46708614]])
Hunt-Pointer-Esterez
LMS_from_linearRGB =
array([[31.3989492 , 63.95129383, 4.64975462],
[15.53714069, 75.78944616, 8.67014186],
[ 1.77515606, 10.94420944, 87.25692246]])
Simulating dichromacy in LMS (Brettel 1997 and Viénot 1999)¶
Let's dive into the seminal paper of {% cite brettel_computerized_1997 %}. {% cite vienot_digital_1999 %} then simplified it for the protanopia and deuteranopia cases, reducing the whole pipeline to a single matrix multiplication. We'll see it right after.
Projection planes in LMS¶
Now that we are in a color space that better matches the human physiology, the question becomes how to transform the LMS coordinates to simulate what a color blind person perceives. As discussed in the Introduction, we want all the colors along the confusion lines to project to the same one, so that a dichromat won't see a difference in the transformed image, but a trichromat will see how much of the color information is lost by the dichromat.
As Ixora.io points out, a first idea is to just set the deficient axis to 0. But that won't work for several reasons:
The new point with L=0, S=0 or M=0 is likely to be outside of the sRGB gamut (the parallelepiped we visualized before in LMS), so most of these transformed color would not be valid once transformed back to sRGB.
Experiments on unilateral dichromats (people who have one normal trichromat eye and one dichromat eye) show that they perceive some colors similarly with both eyes. Unsurprisingly among these invariant stimulus we find shades of gray (black to white). Then the similar colors for the two eyes are observed along the 475nm (blue-greenish) and 575nm (yellow) axes for protans and deutans, and along the 485nm (cyan) and 660nm (red) axes for tritans. {% cite meyer_color_defective_1988 %} has solid references for some of these experiments, made in the 1940s {% cite judd_standard_1944 %}. Intuitively it makes sense that the dichromat eye of a protan or deutan does a good job on colors that only differ by the amount of blue (S cone). Indeed in LMS going from black to blue or white to yellow is just a change in the S coordinate.
Based on that that data it is possible to design a projection that respect these constraints. {% cite brettel_computerized_1997 %} proposed the first modern algorithm to compute that. It consists in creating two half-planes around the Black-White diagonal. For protanopes and deuteranopes the wings go through the LMS coordinates of the 475nm and 575nm wavelengths, and for tritanopes towards the coordinates of the 485nm and 660nm wavelengths. By construction colors that lie on these invariant axes are thus preserved.
While black just has a (0,0,0) coordinate in LMS, "white" needs a more careful definition. In {% cite brettel_computerized_1997 %} they take the white point of the equal-energy illuminant as the neutral axis. That illuminant (called E in CIE XYZ terminology) has the same amount of energy on all the visible light wavelengths, so it presents a perfectly even color. It differs a bit from the "white" of a sRGB monitor which is the white point of the D65 illuminant, which roughly corresponds to the stimulus we'd perceive as white under average outdoor daylight. To stay close to the best theoretical value we'll consider Black to E as the neutral axis that dichromats perceive similarly as trichromats.
In practice most implementations actually use D65 white instead of
Eas the neutral axis, this way it's just sRGB (255,255,255). This is what Viénot 1999 did, but this simplification can apply to the original Brettel 1997 approach too and I recommend it.
Let's first visualize how these planes look like in 3D for protans/deutans (they share the same plane) and for tritans. Note that the colors corresponding to pure 475nm, 485nm, 575nm and 660nm cannot be represented on a sRGB monitor, so I use their closest approximation in the plot by de-saturating them (adding "white" to them until they fit in the gamut).
# collapse-hide
lms_model = convert.LMSModel_sRGB_SmithPokorny75()
from daltonlens import cmfs
def wavelength_to_XYZ_JuddVos(l):
XYZ = colour.wavelength_to_XYZ(475, cmfs)
xyY = convert.xyY_from_XYZ(XYZ)
convert.xy_vos1978_from_xy_CIE1931(xyz.x, xyz.y)
xyz_475 = cmfs.wavelength_to_XYZ_JuddVos(475)
xyz_575 = cmfs.wavelength_to_XYZ_JuddVos(575)
xyz_485 = cmfs.wavelength_to_XYZ_JuddVos(485)
xyz_660 = cmfs.wavelength_to_XYZ_JuddVos(660)
# From https://www.johndcook.com/wavelength_to_RGB.html
# Approximate since the XYZ coordinates are not actually in the RGB gamut
# The correspondence has to be done projecting in the direction of the
# whitepoint until reaching the gamut.
# http://www.fourmilab.ch/documents/specrend/
lms_475 = lms_model.LMS_from_XYZ @ xyz_475
lms_575 = lms_model.LMS_from_XYZ @ xyz_575
lms_485 = lms_model.LMS_from_XYZ @ xyz_485
lms_660 = lms_model.LMS_from_XYZ @ xyz_660
# The equal-energy white point. By construction of CIE XYZ
# it has X=Y=Z. The normalization does not matter to define
# the diagonal direction, picking 0.8 to make it close to
# sRGB white.
xyz_E = [0.8, 0.8, 0.8]
rgb_E = lms_model.linearRGB_from_XYZ @ xyz_E
srgb_E = convert.sRGB_from_linearRGB([rgb_E])
lms_E = lms_model.LMS_from_XYZ @ xyz_E
rgb_wavelengths = convert.apply_color_matrix(np.array([xyz_475, xyz_575, xyz_485, xyz_660]), lms_model.linearRGB_from_XYZ)
rgb_wavelengths = convert.desaturate_linearRGB_to_fit_in_gamut(rgb_wavelengths)
srgb_wavelengths = convert.sRGB_from_linearRGB(rgb_wavelengths)
rgb_475, rgb_575, rgb_485, rgb_660 = srgb_wavelengths
rgb_gamut = [
(np.array([0.0, 0.0, 0.0]), 'K', 'black'),
(np.array([0.0, 0.0, 1.0]), 'B', 'blue'),
(np.array([1.0, 0.0, 1.0]), 'M', 'magenta'),
(np.array([1.0, 0.0, 0.0]), 'R', 'red'),
(np.array([0.0, 1.0, 0.0]), 'G', 'green'),
(np.array([0.0, 1.0, 1.0]), 'C', 'cyan'),
(np.array([1.0, 1.0, 1.0]), 'W', 'white'),
(np.array([1.0, 1.0, 0.0]), 'Y', 'yellow'),
]
def dataframe_row(rgb, lms, short_name:str, name: str):
int_rgb255 = (rgb*255.0).astype(np.uint8)
return [name, short_name, *rgb, *lms, f'#{int_rgb255[0]:02x}{int_rgb255[1]:02x}{int_rgb255[2]:02x}']
def lms_of_key(df, key):
return np.array(df.loc[key, ['L','M','S']], dtype=float)
# RGB parallelepiped gamut in LMS
df_rows = [dataframe_row(row[0], lms_model.LMS_from_linearRGB @ row[0], row[1], row[2]) for row in rgb_gamut]
df_rows.append(dataframe_row(rgb_E, lms_E, 'E', 'Equal-Energy'))
df_rows.append(dataframe_row(rgb_475, lms_475, '475nm', '475nm'))
df_rows.append(dataframe_row(rgb_575, lms_575, '575nm', '575nm'))
df_rows.append(dataframe_row(rgb_485, lms_485, '485nm', '485nm'))
df_rows.append(dataframe_row(rgb_660, lms_660, '660nm', '660nm'))
KBMRGCWY = pd.DataFrame(df_rows, columns=['name', 'short_name', 'R','G','B','L','M','S','colorhex'])
KBMRGCWY.set_index('short_name', inplace=True)
KBMRGCWY['short_name'] = KBMRGCWY.index
KBMRGCWY_red_green = KBMRGCWY[0:-2] # drop 485 and 660nm
KBMRGCWY_tritan = KBMRGCWY.drop(KBMRGCWY.index[-4:-2]) # drop 475 and 575nm
figs = [
px.scatter_3d(KBMRGCWY_red_green, x='L', y='M', z='S',
color='colorhex', color_discrete_map='identity', text='short_name',
width=800, height=800,
title='Protan and Deutan Projection Half-Planes according to (Brettel & Mollon, 1997)'),
px.scatter_3d(KBMRGCWY_tritan, x='L', y='M', z='S',
color='colorhex', color_discrete_map='identity', text='short_name',
width=800, height=800,
title='Tritan Projection Half-Planes according to (Brettel & Mollon, 1997)')
]
# Add the parallelepiped lines
df = KBMRGCWY.loc[['K','B','M','R', 'K', 'G', 'Y', 'W', 'C', 'G', 'Y', 'R', 'M', 'W', 'C', 'B']]
for fig in figs:
fig.add_scatter3d(x=df['L'], y=df['M'], z=df['S'], mode='lines', showlegend=False)
# Add the black-E diagonal, dashed
df = KBMRGCWY.loc[['K','E']]
for fig in figs:
fig.add_scatter3d(x=df['L'], y=df['M'], z=df['S'], mode='lines', showlegend=False, line=dict(width=8, dash='dash', color='gray'), row=1, col=1)
def add_half_plane(fig, df, wavelength_key, col, U_norm, V_norm, delaunayaxis='x'):
# The plane is defined by two vectors starting from the origin
# U is the black-white diagonal, and V is
U = lms_of_key(df, 'E') - lms_of_key(df, 'K')
U = U_norm * (U / np.linalg.norm(U))
V = lms_of_key(df, wavelength_key) - lms_of_key(df, 'K')
N = np.cross(U,V)
V = np.cross(N,U)
V = V_norm * (V / np.linalg.norm(V))
u,v = np.mgrid[0:1:50j, 0:1:50j]
xyz = u[...,np.newaxis]*U + v[...,np.newaxis]*V
x = xyz[...,0].flatten()
y = xyz[...,1].flatten()
z = xyz[...,2].flatten()
surfacecolor = convert.apply_color_matrix(xyz, lms_model.linearRGB_from_LMS)
surfacecolor = convert.desaturate_linearRGB_to_fit_in_gamut(surfacecolor)
surfacecolor = surfacecolor.reshape(-1,3)
#display(surfacecolor)
surfacecolor = convert.sRGB_from_linearRGB(surfacecolor)
fig.add_mesh3d(z=z, x=x, y=y, opacity=1.0, delaunayaxis=delaunayaxis, alphahull=-1, vertexcolor=surfacecolor,
lighting=dict(ambient=1.0,specular=0.0,diffuse=0.0),
showlegend=False)
df = KBMRGCWY_red_green
add_half_plane(figs[0], df, '475nm', 1, np.linalg.norm(lms_of_key(df, 'E'))*1.5, lms_of_key(df, 'B')[2]*3)
add_half_plane(figs[0], df, '575nm', 1, np.linalg.norm(lms_of_key(df, 'E'))*1.5, lms_of_key(df, 'B')[2]*3)
df = KBMRGCWY_tritan
add_half_plane(figs[1], df, '485nm', 2, np.linalg.norm(lms_of_key(df, 'E'))*1.5, np.linalg.norm(lms_of_key(df, 'B'))*2, 'y')
add_half_plane(figs[1], df, '660nm', 2, np.linalg.norm(lms_of_key(df, 'E'))*1.5, np.linalg.norm(lms_of_key(df, 'R'))*0.5, 'y')
camera = dict(
up=dict(x=0, y=0, z=1),
center=dict(x=0, y=0, z=0),
eye=dict(x=-1.75, y=-1.75, z=0.75),
)
for fig in figs:
fig.update_layout(scene_camera=camera, dragmode='orbit')
fig.show()