1 /* dichromat.i
2 * Color transforms to help visualize what colorblind people see.
3 *
4 * Based on
5 * Brettel, Hans, Vienot, Francoise, and Mollon, John D.,
6 * "Computerized simulation of color appearance for dichromats,"
7 * J. Opt. Soc. Am. A Vol. 14, No.10(October), p. 2647-2655 (1997)
8 * http://perso.telecom-paristech.fr/~brettel
9 *
10 * Meyer, Gary W., and Greenberg, Donald P.,
11 * "Color-Defective Vision and Computer Graphics Displays,"
12 * IEEE Computer Graphics & Applications, September 1988, p. 28-40
13 *
14 * See vischeck.com to download tools for java.
15 * See daltonize.org for more tools, discussion.
16 * Good background reading at http://www.cvrl.org/
17 */
18 /* Copyright (c) 2013, David H. Munro.
19 * All rights reserved.
20 * This file is part of yorick (http://yorick.sourceforge.net).
21 * Read the accompanying LICENSE file for details.
22 */
23 require, "lab2rgb.i";
24
25 func dichromat(type, rgb, g, b, cmax=)
26 /* DOCUMENT drgb = dichromat(type, rgb)
27 * or drgb = dichromat(type, r, g, b)
28 *
29 * Return colors DRGB perceived by a dichromat of TYPE corresponding
30 * to screen colors RGB = [r,g,b], or R, G, B. The return value has
31 * the same dimensions as RGB or [R,G,B]. There are three classes
32 * of dichromats (colorblind people): Protrans, Deuterans, and
33 * Tritans, according to whether they lack long, medium, or short
34 * wavelength cone cells (or which type of cell is damaged). For
35 * good background material, see
36 * http://en.wikipedia.org/wiki/Color_blindness or
37 * http://cvrl.ioo.ucl.ac.uk/ and for details of the algorithm
38 * implemented here see Brettel et al., JOSA 14, 2647 (1997),
39 * http://perso.telecom-paristech.fr/~brettel. By default, DRGB is
40 * a char array, but with cmax=1, you can get a double array. The
41 * input RGB can be a char array, or any integer datatype with a
42 * maximum of 255, or a floating point array with maximum 1.0. Both
43 * input and output are sRGB values, that is, the RGB of your
44 * monitor; see the rgb_s2l and rgb_l2s functions for conversions to
45 * physically linear ("gamma corrected") RGB values.
46 *
47 * TYPE is 1 for protran, 2 for deuteran, or 3 for tritan dichromatism.
48 * You may also pass any string beginning with "p", "d", or "t",
49 * respectively. The dichromat function projects the three dimensional
50 * RGB onto a two dimensional subspace, so that there are only two primary
51 * colors instead of three. Note that tritan dichromatism is quite rare,
52 * while some male in a room of 20 probably suffers from either
53 * deuteran or (less likely) protan dichromacy.
54 *
55 * Use the dichromap function to quickly check a pseudocolor picture
56 * made with pli, plf, or plfc to see how it appears to a dichromat.
57 *
58 * SEE ALSO: dichromap, dichromat_setup, rgb_s2l, rgb_l2s
59 */
60 {
61 rgb = rgb_s2l(rgb, g, b);
62 if (structof(type) == string)
63 type = where(strcase(0,strpart(type,1:1)) == ["p", "d", "t"]);
64 type = (numberof(type)==1)? type(1) : 0;
65 if (type<1 || type>3) error, "type must be one of 1,2,3 or p,d,t";
66 proj = _dichromatp(,,,type); /* select projection operators */
67 test = _dichromatt(,type); /* select projection test */
68 test = double(rgb(..,+) * test(+) >= 0.);
69 rgb = rgb(..,+) * proj(,+,); /* perform both projections */
70 rgb = test*rgb(..,1) + (1.-test)*rgb(..,2); /* pick proper projection */
71 return rgb_l2s(rgb, cmax=cmax);
72 }
73
dichromap(type)74 func dichromap(type)
75 /* DOCUMENT dichromap, type
76 * or dichromap
77 * Modify the current colormap to see how the plot (pli, plf, or
78 * plfc) would appear to a dichromat of the given TYPE, which is 1,
79 * 2, 3 (or "p", "d", "t") to simulate protan, deuteran, or tritan
80 * vision. See help, dichromat for more. Omitting TYPE toggles the
81 * original color map with the modified one provided the current
82 * colormap is still the one installed by dichromap.
83 * SEE ALSO: dichromat, cmap
84 */
85 {
86 require, "cmap.i";
87 map = cmap();
88 if (is_void(map)) error, "no current colormap";
89 if (numberof(map)==numberof(_dichro_map))
90 match = allof(map==_dichro_map) | (allof(map==_dichro_orig)<<1);
91 if (is_void(type)) {
92 if (!match) error, "current map not installed by dichromap";
93 } else {
94 if (match == 1) map = _dichro_orig;
95 match = 2;
96 _dichro_orig = map;
97 _dichro_map = transpose(dichromat(type, transpose(map)));
98 }
99 if (match == 1) cmap, _dichro_orig;
100 else if (match == 2) cmap, _dichro_map;
101 }
102
103 func dichromat_setup(xylms, xypri, &test, recurse=)
104 /* DOCUMENT dichromat_setup, xylms, xypri
105 * Set up the dichromatic projection matrices and selection tests for
106 * the dichromat function. XYLMS is either a 2x3 or 3x3 array of the
107 * CIE xy or XYZ color coordinates of the L, M, and S spectral responses
108 * of the three retinal cone cells responsible for color vision. There
109 * is some dispute about exactly what colors these are, since the
110 * experiments to find them are difficult. XYPRI is either a 2x2x3 or
111 * 3x2x3 of the xy or XYZ color coordinates of the pairs of primary
112 * trichromatic colors which will be used to represent what the dichromat
113 * sees. The third index is [protan,deuteran,tritan]. These should be
114 * chosen so that the colors returned by the dichromat function are
115 * indistinguishable from the input colors by that particular type of
116 * dichromat.
117 * SEE ALSO: dichromat
118 */
119 {
120 if (!recurse) {
121 d = dimsof(xylms);
122 if (d(2) == 2) {
123 d = array(0., d+[0,1,0]);
124 d(1:2,) = xylms; d(3,) = 1. - xylms(sum,);
125 xylms = d;
126 }
127 d = dimsof(xypri);
128 if (d(2) == 2) {
129 d = array(0., d+[0,1,0,0]);
130 d(1:2,,) = xypri; d(3,,) = 1. - xypri(sum,,);
131 xypri = d;
132 }
133 local t1, t2, t3;
134 p1 = dichromat_setup(xylms, xypri(,,1), t1, recurse=1);
135 p2 = dichromat_setup(xylms(,[2,3,1]), xypri(,,2), t2, recurse=1);
136 p3 = dichromat_setup(xylms(,[3,1,2]), xypri(,,3), t3, recurse=1);
137 extern _dichromatp, _dichromatt;
138 _dichromatp = [p1, p2, p3];
139 _dichromatt = [t1, t2, t3];
140
141 } else {
142 /* xylms permuted so that missing cone is xylms(,1),
143 * and xypri is the pair of primaries to use for the result
144 * both have full XYZ color coordinates
145 * also need xyz_white from lab2rgb.i as neutral color
146 * (sRGB monitor definition requires this to be D65 white)
147 */
148 /* xyz = xylms(,+) * lms(+) */
149 lmsxyz = LUsolve(xylms); /* lms = lmsxyz(,+) * xyz(+) */
150 p3 = xylms(,+:2:3) * lmsxyz(+:2:3,); /* project out L axis */
151 /* adding anything times xylms(,1) imperceptible to dichromat */
152 /* first, we project white into the MS plane - this line in
153 * MS is the test for which of the two dichromat primaries
154 * we will use; colors that project to one side of the line
155 * get one primary, while those projecting to the other side
156 * get the other color
157 */
158 lmsw = lmsxyz(,+) * xyz_white(+);
159 test = lmsw(3)*lmsxyz(2,) - lmsw(2)*lmsxyz(3,);
160 t1 = test(+) * xypri(+,1);
161 t2 = test(+) * xypri(+,2);
162 if (t1*t2 >= 0.) error, "bad xypri, two primaries on same side of white";
163 /* now we have to figure out how much xylms(,1) to add to p3:
164 * in LMS space, want to add back L = some linear combination of M,S
165 * L = a*M + b*S, where a and b are chosen to make white and one
166 * primary map to themselves:
167 * lmsw(1) = lmsw(+:2:3) * [a,b](+)
168 * lmspri(1,i) = lmspri(+:2:3,i) * [a,b](+)
169 */
170 lmspri = lmsxyz(,+) * xypri(+,);
171 ab1 = LUsolve(transpose([lmsw(2:3),lmspri(2:3,1)]), [lmsw(1),lmspri(1,1)]);
172 ab2 = LUsolve(transpose([lmsw(2:3),lmspri(2:3,2)]), [lmsw(1),lmspri(1,2)]);
173 p1 = p3 + xylms(,1) * (ab1(-,+) * lmsxyz(+:2:3,));
174 p2 = p3 + xylms(,1) * (ab2(-,+) * lmsxyz(+:2:3,));
175 /* finally, transform everything to lRGB coordinates from XYZ
176 * rgb = _rgb_xyz(,+) * xyz(+)
177 * xyz = _xyz_rgb(,+) * rgb(+)
178 */
179 test = test(+) * _xyz_rgb(+,);
180 p1 = (_rgb_xyz(,+) * p1(+,))(,+) * _xyz_rgb(+,);
181 p2 = (_rgb_xyz(,+) * p2(+,))(,+) * _xyz_rgb(+,);
182 return (t1 > 0.)? [p1, p2] : [p2, p1];
183 }
184 }
185
186 /* Brettel has suggestions for xypri: For both protan and deuteran
187 * dichromats, he chooses monochromatic 475 (blue) and 575 (yellow) as
188 * the two colors, while for tritan dichromats, he chooses 485 (green)
189 * and 660 (red).
190 *
191 * Meyer says protan axis is 473 nm to D65 to 574 nm,
192 * deuteran axis is 477 to D65 to 578, tritan axis is 490 to D65 to 610.
193 *
194 * VC vischeck.com uses somewhat different choices
195 */
196 xypri_brettel = [[[0.109594,0.0868425], [0.478775,0.520202]],
197 [[0.109594,0.0868425], [0.478775,0.520202]],
198 [[0.0687059,0.200723], [0.729969,0.270031]]];
199 xypri_meyer = [[[0.116091,0.073853], [0.47197,0.526967]],
200 [[0.103188,0.102897], [0.49913,0.499907]],
201 [[0.0453907,0.294976],[0.665764,0.334011]]];
202 xypri_vc = [[[0.123165,0.088141], [0.476373,0.542367]],
203 [[0.123165,0.088141], [0.476373,0.542367]],
204 [[0.111112,0.207644], [0.840774,0.202248]]];
205
206 /* Meyer "Color-Defective Vision and Computer Graphics Displays" (1988)
207 * Hunt-Pointer-Estevez (HPE) used in RLAB color appearance model
208 * Bradford used in CIE CAT02, CAT97s, CNCCAT97
209 * spectral sharpening does not change xylms
210 * normalizing to D65 white does not change xylms
211 * VC = VisCheck, vischeck.com
212 * Daltonize = daltonize.org, very close to Meyer
213 */
214 xylms_meyer = [[0.735,0.265], [ 1.14, -0.14], [0.171,-0.003]];
215 xylms_hpe = [[0.837,0.163], [ 2.30, -1.30], [0.168, 0]];
216 xylms_bradford = [[0.700,0.306], [-0.286,1.196], [0.137, 0.043]];
217 xylms_vc = [[0.883854,0.161895], [-13.6498,13.7753], [0.153684,-0.0288546]];
218 xylms_daltonize = [[0.7465, 0.2535], [1.4, -0.4], [0.1748, 0.]];
219
220 dichromat_setup, xylms_vc, xypri_vc;
221
222 /* after daltonize.py from daltonize.org and http://moinmo.in/AccessibleMoin */
223 func daltonize(type, rgb, g, b, cmax=)
224 {
225 rgb = rgb_s2l(rgb, g, b);
226 lrgb_skip = 1;
227 err = rgb - dichromat(type, rgb);
228 lrgb_skip = [];
229 rgb += err(..,+) * daltonizer(+,);
230 return rgb_l2s(min(max(rgb, 0.), 1.), cmax=cmax);
231 }
232 daltonizer = [[0,0,0], [0.7,1,0], [0.7,0,1]];
233