Fix README.md
[gencolormap.git] / colormap.cpp
1 /*
2  * Copyright (C) 2015, 2016 Computer Graphics Group, University of Siegen
3  * Written by Martin Lambers <martin.lambers@uni-siegen.de>
4  *
5  * Permission is hereby granted, free of charge, to any person obtaining a copy
6  * of this software and associated documentation files (the "Software"), to deal
7  * in the Software without restriction, including without limitation the rights
8  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
9  * copies of the Software, and to permit persons to whom the Software is
10  * furnished to do so, subject to the following conditions:
11  *
12  * The above copyright notice and this permission notice shall be included in
13  * all copies or substantial portions of the Software.
14  *
15  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL THE
18  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
20  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
21  * SOFTWARE.
22  */
23
24 #include <algorithm>
25 #include <vector>
26 #include <limits>
27 #include <cmath>
28 #include <cstring>
29
30 #include "colormap.hpp"
31
32 /* Notes about the color spaces used internally:
33  *
34  * - We use D65 white everywhere
35  * - RGB means linear RGB; we also have sRGB
36  * - RGB and sRGB values are in [0,1]
37  * - XYZ, LUV, and similar values are in the original range (not normalized);
38  *   often this is [0,100]
39  * - All angles (for hue) are measured in radians
40  */
41
42 namespace ColorMap {
43
44 /* Generic helpers */
45
46 static const float pi = M_PI;
47 static const float twopi = 2.0 * M_PI;
48
49 static float clamp(float x, float lo, float hi)
50 {
51     return std::min(std::max(x, lo), hi);
52 }
53
54 static float hue_diff(float h0, float h1)
55 {
56     float t = std::fabs(h1 - h0);
57     return (t < pi ? t : twopi - t);
58 }
59
60 static float uchar_to_float(unsigned char x)
61 {
62     return x / 255.0f;
63 }
64
65 static unsigned char float_to_uchar(float x, bool* clipped = NULL)
66 {
67     int v = std::round(x * 255.0f);
68     if (v < 0) {
69         v = 0;
70         if (clipped)
71             *clipped = true;
72     } else if (v > 255) {
73         v = 255;
74         if (clipped)
75             *clipped = true;
76     } else {
77         if (clipped)
78             *clipped = false;
79     }
80     return v;
81 }
82
83 /* A color triplet class without assumptions about the color space */
84
85 class triplet {
86 public:
87     struct {
88         union { float x, l,       m, r; };
89         union { float y, a, u, c, s, g; };
90         union { float z, b, v, h      ; };
91     };
92
93     triplet() {}
94     triplet(float t0, float t1, float t2) : x(t0), y(t1), z(t2) {}
95     triplet(const triplet& t) : x(t.x), y(t.y), z(t.z) {}
96 };
97
98 triplet operator+(triplet t0, triplet t1)
99 {
100     return triplet(t0.x + t1.x, t0.y + t1.y, t0.z + t1.z);
101 }
102
103 triplet operator*(float s, triplet t)
104 {
105     return triplet(s * t.x, s * t.y, s * t.z);
106 }
107
108 /* XYZ and related color spaces helper functions and values */
109
110 static float u_prime(triplet xyz)
111 {
112     return 4.0f * xyz.x / (xyz.x + 15.0f * xyz.y + 3.0f * xyz.z);
113 }
114
115 static float v_prime(triplet xyz)
116 {
117     return 9.0f * xyz.y / (xyz.x + 15.0f * xyz.y + 3.0f * xyz.z);
118 }
119
120 static const triplet d65_xyz = triplet(95.047f, 100.000f, 108.883f);
121 static const float d65_u_prime = u_prime(d65_xyz);
122 static const float d65_v_prime = v_prime(d65_xyz);
123
124 static triplet adjust_y(triplet xyz, float new_y)
125 {
126     float sum = xyz.x + xyz.y + xyz.z;
127     // keep old chromaticity in terms of x, y
128     float x = xyz.x / sum;
129     float y = xyz.y / sum;
130     // apply new luminance
131     float r = new_y / y;
132     return triplet(r * x, new_y, r * (1.0f - x - y));
133 }
134
135 /* Color space conversion: LCH <-> LUV */
136
137 static triplet lch_to_luv(triplet lch)
138 {
139     return triplet(lch.l, lch.c * std::cos(lch.h), lch.c * std::sin(lch.h));
140 }
141
142 static triplet luv_to_lch(triplet luv)
143 {
144     triplet lch(luv.l, std::hypot(luv.u, luv.v), std::atan2(luv.v, luv.u));
145     if (lch.h < 0.0f)
146         lch.h += twopi;
147     return lch;
148 }
149
150 static float lch_saturation(float l, float c)
151 {
152     return c / std::max(l, 1e-8f);
153 }
154
155 static float lch_chroma(float l, float s)
156 {
157     return s * l;
158 }
159
160 /* Color space conversion: LUV <-> XYZ */
161
162 static triplet luv_to_xyz(triplet luv)
163 {
164     triplet xyz;
165     float u_prime = luv.u / (13.0f * luv.l) + d65_u_prime;
166     float v_prime = luv.v / (13.0f * luv.l) + d65_v_prime;
167     if (luv.l <= 8.0f) {
168         xyz.y = d65_xyz.y * luv.l * (3.0f * 3.0f * 3.0f / (29.0f * 29.0f * 29.0f));
169     } else {
170         float tmp = (luv.l + 16.0f) / 116.0f;
171         xyz.y = d65_xyz.y * tmp * tmp * tmp;
172     }
173     xyz.x = xyz.y * (9.0f * u_prime) / (4.0f * v_prime);
174     xyz.z = xyz.y * (12.0f - 3.0f * u_prime - 20.0f * v_prime) / (4.0f * v_prime);
175     return xyz;
176 }
177
178 static triplet xyz_to_luv(triplet xyz)
179 {
180     triplet luv;
181     float y_ratio = xyz.y / d65_xyz.y;
182     if (y_ratio <= (6.0f * 6.0f * 6.0f) / (29.0f * 29.0f * 29.0f)) {
183         luv.l = (29.0f * 29.0f * 29.0f) / (3.0f * 3.0f * 3.0f) * y_ratio;
184     } else {
185         luv.l = 116.0f * std::cbrt(y_ratio) - 16.0f;
186     }
187     luv.u = 13.0f * luv.l * (u_prime(xyz) - d65_u_prime);
188     luv.v = 13.0f * luv.l * (v_prime(xyz) - d65_v_prime);
189     return luv;
190 }
191
192 static float luv_saturation(triplet luv)
193 {
194     return lch_saturation(luv.l, std::hypot(luv.u, luv.v));
195 }
196
197 /* Color space conversion: LAB <-> XYZ */
198
199 static float lab_invf(float t)
200 {
201     if (t > 6.0f / 29.0f)
202         return t * t * t;
203     else
204         return (3.0f * 6.0f * 6.0f) / (29.0f * 29.0f) * (t - 4.0f / 29.0f);
205 }
206
207 static triplet lab_to_xyz(triplet lab)
208 {
209     float t = (lab.l + 16.0f) / 116.0f;
210     return triplet(d65_xyz.x * lab_invf(t + lab.a / 500.0f),
211                    d65_xyz.y * lab_invf(t),
212                    d65_xyz.z * lab_invf(t - lab.b / 200.0f));
213 }
214
215 static float lab_f(float t)
216 {
217     if (t > (6.0f * 6.0f * 6.0f) / (29.0f * 29.0f * 29.0f))
218         return std::cbrt(t);
219     else
220         return (29.0f * 29.0f) / (3.0f * 6.0f * 6.0f) * t + 4.0f / 29.0f;
221 }
222
223 static triplet xyz_to_lab(triplet xyz)
224 {
225     triplet fxyz = triplet(lab_f(xyz.x / d65_xyz.x), lab_f(xyz.y / d65_xyz.y), lab_f(xyz.z / d65_xyz.z));
226     return triplet(116.0f * fxyz.y - 16.0f, 500.0f * (fxyz.x - fxyz.y), 200.0f * (fxyz.y - fxyz.z));
227 }
228
229 /* Color space conversion: RGB <-> XYZ */
230
231 static triplet rgb_to_xyz(triplet rgb)
232 {
233     return 100.0f * triplet(
234             (0.4124f * rgb.r + 0.3576f * rgb.g + 0.1805f * rgb.b),
235             (0.2126f * rgb.r + 0.7152f * rgb.g + 0.0722f * rgb.b),
236             (0.0193f * rgb.r + 0.1192f * rgb.g + 0.9505f * rgb.b));
237 }
238
239 static triplet xyz_to_rgb(triplet xyz)
240 {
241     return 0.01f * triplet(
242             (+3.2406255f * xyz.x - 1.5372080f * xyz.y - 0.4986286f * xyz.z),
243             (-0.9689307f * xyz.x + 1.8757561f * xyz.y + 0.0415175f * xyz.z),
244             (+0.0557101f * xyz.x - 0.2040211f * xyz.y + 1.0569959f * xyz.z));
245 }
246
247 /* Color space conversion: RGB <-> sRGB */
248
249 static float rgb_to_srgb_helper(float x)
250 {
251     return (x <= 0.0031308f ? (x * 12.92f) : (1.055f * std::pow(x, 1.0f / 2.4f) - 0.055f));
252 }
253
254 static triplet rgb_to_srgb(triplet rgb)
255 {
256     return triplet(rgb_to_srgb_helper(rgb.r), rgb_to_srgb_helper(rgb.g), rgb_to_srgb_helper(rgb.b));
257 }
258
259 static float srgb_to_rgb_helper(float x)
260 {
261     return (x <= 0.04045f ? (x / 12.92f) : std::pow((x + 0.055f) / 1.055f, 2.4f));
262 }
263
264 static triplet srgb_to_rgb(triplet srgb)
265 {
266     return triplet(srgb_to_rgb_helper(srgb.r), srgb_to_rgb_helper(srgb.g), srgb_to_rgb_helper(srgb.b));
267 }
268
269 /* Helpers for the conversion to colormap entries */
270
271 static bool xyz_to_colormap(triplet xyz, unsigned char* colormap)
272 {
273     bool clipped[3];
274     triplet srgb = rgb_to_srgb(xyz_to_rgb(xyz));
275     colormap[0] = float_to_uchar(srgb.r, clipped + 0);
276     colormap[1] = float_to_uchar(srgb.g, clipped + 1);
277     colormap[2] = float_to_uchar(srgb.b, clipped + 2);
278     return clipped[0] || clipped[1] || clipped[2];
279 }
280
281 static bool luv_to_colormap(triplet luv, unsigned char* colormap)
282 {
283     return xyz_to_colormap(luv_to_xyz(luv), colormap);
284 }
285
286 static bool lch_to_colormap(triplet lch, unsigned char* colormap)
287 {
288     return luv_to_colormap(lch_to_luv(lch), colormap);
289 }
290
291 static bool lab_to_colormap(triplet lab, unsigned char* colormap)
292 {
293     return xyz_to_colormap(lab_to_xyz(lab), colormap);
294 }
295
296 /* Various helpers */
297
298 static float srgb_to_lch_hue(triplet srgb)
299 {
300     return luv_to_lch(xyz_to_luv(rgb_to_xyz(srgb_to_rgb(srgb)))).h;
301 }
302
303 // Compute most saturated color that fits into the sRGB
304 // cube for the given LCH hue value. This is the core
305 // of the Wijffelaars paper.
306 static triplet most_saturated_in_srgb(float lch_hue)
307 {
308     /* Static values, only computed once */
309     static float h[] = {
310         srgb_to_lch_hue(triplet(1, 0, 0)),
311         srgb_to_lch_hue(triplet(1, 1, 0)),
312         srgb_to_lch_hue(triplet(0, 1, 0)),
313         srgb_to_lch_hue(triplet(0, 1, 1)),
314         srgb_to_lch_hue(triplet(0, 0, 1)),
315         srgb_to_lch_hue(triplet(1, 0, 1))
316     };
317
318     /* RGB values and variable pointers to them */
319     int i, j, k;
320     if (lch_hue < h[0]) {
321         i = 2;
322         j = 1;
323         k = 0;
324     } else if (lch_hue < h[1]) {
325         i = 1;
326         j = 2;
327         k = 0;
328     } else if (lch_hue < h[2]) {
329         i = 0;
330         j = 2;
331         k = 1;
332     } else if (lch_hue < h[3]) {
333         i = 2;
334         j = 0;
335         k = 1;
336     } else if (lch_hue < h[4]) {
337         i = 1;
338         j = 0;
339         k = 2;
340     } else if (lch_hue < h[5]) {
341         i = 0;
342         j = 1;
343         k = 2;
344     } else {
345         i = 2;
346         j = 1;
347         k = 0;
348     }
349
350     /* Compute the missing component */
351     float srgb[3];
352     float M[3][3] = {
353         { 0.4124f, 0.3576f, 0.1805f },
354         { 0.2126f, 0.7152f, 0.0722f },
355         { 0.0193f, 0.1192f, 0.9505f }
356     };
357     float alpha = -std::sin(lch_hue);
358     float beta = std::cos(lch_hue);
359     float T = alpha * d65_u_prime + beta * d65_v_prime;
360     srgb[j] = 0.0f;
361     srgb[k] = 1.0f;
362     float q0 = T * (M[0][k] + 15.0f * M[1][k] + 3.0f * M[2][k]) - (4.0f * alpha * M[0][k] + 9.0f * beta * M[1][k]);
363     float q1 = T * (M[0][i] + 15.0f * M[1][i] + 3.0f * M[2][i]) - (4.0f * alpha * M[0][i] + 9.0f * beta * M[1][i]);
364     srgb[i] = rgb_to_srgb_helper(clamp(- q0 / q1, 0.0f, 1.0f));
365
366     /* Convert back to LUV */
367     return xyz_to_luv(rgb_to_xyz(srgb_to_rgb(triplet(srgb[0], srgb[1], srgb[2]))));
368 }
369
370 static float s_max(float l, float h)
371 {
372     triplet pmid = most_saturated_in_srgb(h);
373     triplet pend = triplet(0.0f, 0.0f, 0.0f);
374     if (l > pmid.l)
375         pend.l = 100.0f;
376     float alpha = (pend.l - l) / (pend.l - pmid.l);
377     float pmids = luv_saturation(pmid);
378     float pends = luv_saturation(pend);
379     return alpha * (pmids - pends) + pends;
380 }
381
382 static triplet get_bright_point()
383 {
384     static triplet pb(-1.0f, -1.0f, -1.0f);
385     if (pb.l < 0.0f) {
386         pb = xyz_to_luv(rgb_to_xyz(triplet(1.0f, 1.0f, 0.0f)));
387     }
388     return pb;
389 }
390
391 static float mix_hue(float alpha, float h0, float h1)
392 {
393     float M = std::fmod(pi + h1 - h0, twopi) - pi;
394     return std::fmod(h0 + alpha * M, twopi);
395 }
396
397 static void get_color_points(float hue, float saturation, float warmth,
398         triplet pb, float pb_hue, float pb_saturation,
399         triplet* p0, triplet* p1, triplet* p2,
400         triplet* q0, triplet* q1, triplet* q2)
401 {
402     *p0 = lch_to_luv(triplet(0.0f, 0.0f, hue));
403     *p1 = most_saturated_in_srgb(hue);
404     triplet p2_lch;
405     p2_lch.l = (1.0f - warmth) * 100.0f + warmth * pb.l;
406     p2_lch.h = mix_hue(warmth, hue, pb_hue);
407     p2_lch.c = lch_chroma(p2_lch.l, std::min(s_max(p2_lch.l, p2_lch.h), warmth * saturation * pb_saturation));
408     *p2 = lch_to_luv(p2_lch);
409     *q0 = (1.0f - saturation) * (*p0) + saturation * (*p1);
410     *q2 = (1.0f - saturation) * (*p2) + saturation * (*p1);
411     *q1 = 0.5f * ((*q0) + (*q2));
412 }
413
414 static triplet b(triplet b0, triplet b1, triplet b2, float t)
415 {
416     float a = (1.0f - t) * (1.0f - t);
417     float b = 2.0f * (1.0f - t) * t;
418     float c = t * t;
419     return a * b0 + b * b1 + c * b2;
420 }
421
422 static float inv_b(float b0, float b1, float b2, float v)
423 {
424     return (b0 - b1 + std::sqrt(std::max(b1 * b1 - b0 * b2 + (b0 - 2.0f * b1 + b2) * v, 0.0f)))
425         / (b0 - 2.0f * b1 + b2);
426 }
427
428 static triplet get_colormap_entry(float t,
429         triplet p0, triplet p2,
430         triplet q0, triplet q1, triplet q2,
431         float contrast, float brightness)
432 {
433     float l = 125 - 125 * std::pow(0.2f, (1.0f - contrast) * brightness + t * contrast);
434     float T = (l <= q1.l ? 0.5f * inv_b(p0.l, q0.l, q1.l, l) : 0.5f * inv_b(q1.l, q2.l, p2.l, l) + 0.5f);
435     return (T <= 0.5f ? b(p0, q0, q1, 2.0f * T) : b(q1, q2, p2, 2.0f * (T - 0.5f)));
436 }
437
438 /* Brewer-like color maps */
439
440 float BrewerSequentialDefaultContrastForSmallN(int n)
441 {
442     return std::min(0.88f, 0.34f + 0.06f * n);
443 }
444
445 int BrewerSequential(int n, unsigned char* colormap, float hue,
446         float contrast, float saturation, float brightness, float warmth)
447 {
448     triplet pb, p0, p1, p2, q0, q1, q2;
449     pb = get_bright_point();
450     triplet pb_lch = luv_to_lch(pb);
451     float pbs = lch_saturation(pb_lch.l, pb_lch.c);
452     get_color_points(hue, saturation, warmth, pb, pb_lch.h, pbs, &p0, &p1, &p2, &q0, &q1, &q2);
453
454     int clipped = 0;
455     for (int i = 0; i < n; i++) {
456         float t = i / (n - 1.0f);
457         triplet c = get_colormap_entry(t, p0, p2, q0, q1, q2, contrast, brightness);
458         if (luv_to_colormap(c, colormap + 3 * i))
459             clipped++;
460     }
461     return clipped;
462 }
463
464 float BrewerDivergingDefaultContrastForSmallN(int n)
465 {
466     return std::min(0.88f, 0.34f + 0.06f * n);
467 }
468
469 int BrewerDiverging(int n, unsigned char* colormap, float hue, float divergence,
470         float contrast, float saturation, float brightness, float warmth)
471 {
472     float hue1 = hue + divergence;
473     if (hue1 >= twopi)
474         hue1 -= twopi;
475
476     triplet pb;
477     triplet p00, p01, p02, q00, q01, q02;
478     triplet p10, p11, p12, q10, q11, q12;
479     pb = get_bright_point();
480     triplet pb_lch = luv_to_lch(pb);
481     float pbs = lch_saturation(pb_lch.l, pb_lch.c);
482     get_color_points(hue,  saturation, warmth, pb, pb_lch.h, pbs, &p00, &p01, &p02, &q00, &q01, &q02);
483     get_color_points(hue1, saturation, warmth, pb, pb_lch.h, pbs, &p10, &p11, &p12, &q10, &q11, &q12);
484
485     int clipped = 0;
486     for (int i = 0; i < n; i++) {
487         triplet c;
488         if (n % 2 == 1 && i == n / 2) {
489             // compute neutral color in the middle of the map
490             triplet c0 = get_colormap_entry(1.0f, p00, p02, q00, q01, q02, contrast, brightness);
491             triplet c1 = get_colormap_entry(1.0f, p10, p12, q10, q11, q12, contrast, brightness);
492             if (n <= 9) {
493                 // for discrete color maps, use an extra neutral color
494                 float c0s = luv_saturation(c0);
495                 float c1s = luv_saturation(c1);
496                 float sn = 0.5f * (c0s + c1s) * warmth;
497                 c.l = 0.5f * (c0.l + c1.l);
498                 float cc = lch_chroma(c.l, std::min(s_max(c.l, pb_lch.h), sn));
499                 c = lch_to_luv(triplet(c.l, cc, pb_lch.h));
500             } else {
501                 // for continuous color maps, use an average, since the extra neutral color looks bad
502                 c = 0.5f * (c0 + c1);
503             }
504         } else {
505             float t = i / (n - 1.0f);
506             if (i < n / 2) {
507                 float tt = 2.0f * t;
508                 c = get_colormap_entry(tt, p00, p02, q00, q01, q02, contrast, brightness);
509             } else {
510                 float tt = 2.0f * (1.0f - t);
511                 c = get_colormap_entry(tt, p10, p12, q10, q11, q12, contrast, brightness);
512             }
513         }
514         if (luv_to_colormap(c, colormap + 3 * i))
515             clipped++;
516     }
517     return clipped;
518 }
519
520 int BrewerQualitative(int n, unsigned char* colormap, float hue, float divergence,
521         float contrast, float saturation, float brightness)
522 {
523     // Get all information about yellow
524     static triplet ylch(-1.0f, -1.0f, -1.0f);
525     if (ylch.l < 0.0f) {
526         ylch = luv_to_lch(xyz_to_luv(rgb_to_xyz(triplet(1.0f, 1.0f, 0.0f))));
527     }
528
529     // Get saturation of red (maximum possible saturation)
530     static float rs = -1.0f;
531     if (rs < 0.0f) {
532         triplet rluv = xyz_to_luv(rgb_to_xyz(triplet(1.0f, 0.0f, 0.0f)));
533         rs = luv_saturation(rluv);
534     }
535
536     // Derive parameters of the method
537     float eps = hue / twopi;
538     float r = divergence / twopi;
539     float l0 = brightness * ylch.l;
540     float l1 = (1.0f - contrast) * l0;
541
542     // Generate colors
543     int clipped = 0;
544     for (int i = 0; i < n; i++) {
545         float t = i / (n - 1.0f);
546         float ch = std::fmod(twopi * (eps + t * r), twopi);
547         float alpha = hue_diff(ch, ylch.h) / pi;
548         float cl = (1.0f - alpha) * l0 + alpha * l1;
549         float cs = std::min(s_max(cl, ch), saturation * rs);
550         triplet c = lch_to_luv(triplet(cl, lch_chroma(cl, cs), ch));
551         if (luv_to_colormap(c, colormap + 3 * i))
552             clipped++;
553     }
554     return clipped;
555 }
556
557 /* Perceptually linear (PL) */
558
559 int PLSequentialLightness(int n, unsigned char* colormap, float saturation, float hue)
560 {
561     triplet lch;
562     lch.h = hue;
563     int clipped = 0;
564     for (int i = 0; i < n; i++) {
565         float t = i / (n - 1.0f);
566         lch.l = std::max(0.01f, t * 100.0f);
567         lch.c = lch_chroma(lch.l, saturation * 5.0f * (1.0f - t));
568         if (lch_to_colormap(lch, colormap + 3 * i))
569             clipped++;
570     }
571     return clipped;
572 }
573
574 int PLSequentialSaturation(int n, unsigned char* colormap,
575         float lightness, float saturation, float hue)
576 {
577     triplet lch;
578     lch.l = std::max(0.01f, lightness * 100.0f);
579     lch.h = hue;
580     int clipped = 0;
581     for (int i = 0; i < n; i++) {
582         float t = i / (n - 1.0f);
583         lch.c = lch_chroma(lch.l, saturation * 5.0f * (1.0f - t));
584         if (lch_to_colormap(lch, colormap + 3 * i))
585             clipped++;
586     }
587     return clipped;
588 }
589
590 int PLSequentialRainbow(int n, unsigned char* colormap, float hue, float rotations, float saturation)
591 {
592     triplet lch;
593     int clipped = 0;
594     for (int i = 0; i < n; i++) {
595         float t = i / (n - 1.0f);
596         lch.l = std::max(0.01f, t * 100.0f);
597         lch.c = lch_chroma(lch.l, (1.0f - t) * saturation);
598         lch.h = hue + t * rotations * twopi;
599         if (lch_to_colormap(lch, colormap + 3 * i))
600             clipped++;
601     }
602     return clipped;
603 }
604
605 static float plancks_law(float temperature, float lambda)
606 {
607     const float c = 299792458.0f;     // speed of light in vacuum
608     //const float c = 299700000.0f;     // speed of visible light in air
609     const float h = 6.626070041e-34f; // Planck's constant
610     const float k = 1.38064853e-23f;  // Boltzmann constant
611     return 2.0f * h * c * c * std::pow(lambda, -5.0f)
612         / (std::exp(h * c / (lambda * k * temperature)) - 1.0f);
613
614 }
615
616 #if 0
617 static float sqr(float x) { return x * x; }
618 static triplet color_matching_function_approx(float lambda)
619 {
620     // Analytic approximation of the CIE 1964 10 deg standard observer CMF.
621     // Taken from the paper "Simple Analytic Approximations to the CIE XYZ
622     // Color Matching Functions" by Wyman, Sloan, Shirley, JCGT 2(2), 2013.
623     lambda *= 1e9f; // convert to nanometers
624     triplet xyz;
625     xyz.x = 0.398f * std::exp(-1250.0f * sqr(std::log((lambda + 570.1f) / 1014.0f)))
626           + 1.132f * std::exp(-234.0f  * sqr(std::log((1338.0f - lambda) / 743.5f)));
627     xyz.y = 1.011f * std::exp(-0.5f    * sqr((lambda - 556.1f) / 46.14f));
628     xyz.z = 2.060f * std::exp(-32.0f   * sqr(std::log((lambda - 265.8f) / 180.4f)));
629     return xyz;
630 }
631 #endif
632
633 static triplet color_matching_function(int lambda /* in nanometers */)
634 {
635     // Tables in 5nm resolution for the CIE 1931 2 deg Standard Observer CMF,
636     // from lambda=380 to lambda=780. Obtained from
637     // http://www.cie.co.at/publ/abst/datatables15_2004/CIE_sel_colorimetric_tables.xls
638     // on 2016-02-10.
639     static const triplet cmf_xyz[] = {
640         triplet(0.001368, 0.000039, 0.006450),
641         triplet(0.002236, 0.000064, 0.010550),
642         triplet(0.004243, 0.000120, 0.020050),
643         triplet(0.007650, 0.000217, 0.036210),
644         triplet(0.014310, 0.000396, 0.067850),
645         triplet(0.023190, 0.000640, 0.110200),
646         triplet(0.043510, 0.001210, 0.207400),
647         triplet(0.077630, 0.002180, 0.371300),
648         triplet(0.134380, 0.004000, 0.645600),
649         triplet(0.214770, 0.007300, 1.039050),
650         triplet(0.283900, 0.011600, 1.385600),
651         triplet(0.328500, 0.016840, 1.622960),
652         triplet(0.348280, 0.023000, 1.747060),
653         triplet(0.348060, 0.029800, 1.782600),
654         triplet(0.336200, 0.038000, 1.772110),
655         triplet(0.318700, 0.048000, 1.744100),
656         triplet(0.290800, 0.060000, 1.669200),
657         triplet(0.251100, 0.073900, 1.528100),
658         triplet(0.195360, 0.090980, 1.287640),
659         triplet(0.142100, 0.112600, 1.041900),
660         triplet(0.095640, 0.139020, 0.812950),
661         triplet(0.057950, 0.169300, 0.616200),
662         triplet(0.032010, 0.208020, 0.465180),
663         triplet(0.014700, 0.258600, 0.353300),
664         triplet(0.004900, 0.323000, 0.272000),
665         triplet(0.002400, 0.407300, 0.212300),
666         triplet(0.009300, 0.503000, 0.158200),
667         triplet(0.029100, 0.608200, 0.111700),
668         triplet(0.063270, 0.710000, 0.078250),
669         triplet(0.109600, 0.793200, 0.057250),
670         triplet(0.165500, 0.862000, 0.042160),
671         triplet(0.225750, 0.914850, 0.029840),
672         triplet(0.290400, 0.954000, 0.020300),
673         triplet(0.359700, 0.980300, 0.013400),
674         triplet(0.433450, 0.994950, 0.008750),
675         triplet(0.512050, 1.000000, 0.005750),
676         triplet(0.594500, 0.995000, 0.003900),
677         triplet(0.678400, 0.978600, 0.002750),
678         triplet(0.762100, 0.952000, 0.002100),
679         triplet(0.842500, 0.915400, 0.001800),
680         triplet(0.916300, 0.870000, 0.001650),
681         triplet(0.978600, 0.816300, 0.001400),
682         triplet(1.026300, 0.757000, 0.001100),
683         triplet(1.056700, 0.694900, 0.001000),
684         triplet(1.062200, 0.631000, 0.000800),
685         triplet(1.045600, 0.566800, 0.000600),
686         triplet(1.002600, 0.503000, 0.000340),
687         triplet(0.938400, 0.441200, 0.000240),
688         triplet(0.854450, 0.381000, 0.000190),
689         triplet(0.751400, 0.321000, 0.000100),
690         triplet(0.642400, 0.265000, 0.000050),
691         triplet(0.541900, 0.217000, 0.000030),
692         triplet(0.447900, 0.175000, 0.000020),
693         triplet(0.360800, 0.138200, 0.000010),
694         triplet(0.283500, 0.107000, 0.000000),
695         triplet(0.218700, 0.081600, 0.000000),
696         triplet(0.164900, 0.061000, 0.000000),
697         triplet(0.121200, 0.044580, 0.000000),
698         triplet(0.087400, 0.032000, 0.000000),
699         triplet(0.063600, 0.023200, 0.000000),
700         triplet(0.046770, 0.017000, 0.000000),
701         triplet(0.032900, 0.011920, 0.000000),
702         triplet(0.022700, 0.008210, 0.000000),
703         triplet(0.015840, 0.005723, 0.000000),
704         triplet(0.011359, 0.004102, 0.000000),
705         triplet(0.008111, 0.002929, 0.000000),
706         triplet(0.005790, 0.002091, 0.000000),
707         triplet(0.004109, 0.001484, 0.000000),
708         triplet(0.002899, 0.001047, 0.000000),
709         triplet(0.002049, 0.000740, 0.000000),
710         triplet(0.001440, 0.000520, 0.000000),
711         triplet(0.001000, 0.000361, 0.000000),
712         triplet(0.000690, 0.000249, 0.000000),
713         triplet(0.000476, 0.000172, 0.000000),
714         triplet(0.000332, 0.000120, 0.000000),
715         triplet(0.000235, 0.000085, 0.000000),
716         triplet(0.000166, 0.000060, 0.000000),
717         triplet(0.000117, 0.000042, 0.000000),
718         triplet(0.000083, 0.000030, 0.000000),
719         triplet(0.000059, 0.000021, 0.000000),
720         triplet(0.000042, 0.000015, 0.000000),
721     };
722     triplet xyz(0, 0, 0);
723     if (lambda >= 380 && lambda <= 780) {
724         int i = (lambda - 380) / 5;
725         xyz = cmf_xyz[i];
726         if (lambda % 5 != 0) {
727             int i1 = i + 1;
728             triplet xyz1 = cmf_xyz[i1];
729             float alpha = (lambda % 5) / 5.0f;
730             xyz = (1.0f - alpha) * xyz + alpha * xyz1;
731         }
732     }
733     return xyz;
734 }
735
736 int PLSequentialBlackBody(int n, unsigned char* colormap, float temperature, float range, float saturation)
737 {
738     int clipped = 0;
739     for (int i = 0; i < n; i++) {
740         float fract = i / (n - 1.0f);
741         // Black body temperature for this color map entry
742         float t = temperature + fract * range;
743         // Integrate radiance over the visible spectrum; according
744         // to literature, sampling at 10nm intervals is enough.
745         triplet xyz(0, 0, 0);
746         int stepsize = 5;
747         float s = stepsize * 1e-9f; // stepsize in meters
748         for (int lambda = 360; lambda <= 830; lambda += stepsize) {
749             float l = lambda * 1e-9f; // lambda in in meters
750             float radiosity = pi * plancks_law(t, l);
751             //xyz = xyz + s * radiosity * color_matching_function_approx(l);
752             xyz = xyz + s * radiosity * color_matching_function(lambda);
753         }
754         triplet lch = luv_to_lch(xyz_to_luv(adjust_y(xyz, 10.0f)));
755         lch.l = std::max(0.01f, fract * 100.0f);
756         lch.c = lch_chroma(lch.l, (1.0f - fract) * saturation);
757         if (lch_to_colormap(lch, colormap + 3 * i))
758             clipped++;
759     }
760     return clipped;
761 }
762
763 int PLDivergingLightness(int n, unsigned char* colormap, float lightness, float saturation, float hue, float divergence)
764 {
765     triplet lch;
766     int clipped = 0;
767     for (int i = 0; i < n; i++) {
768         float t = i / (n - 1.0f);
769         lch.l = std::max(0.01f, 100.0f * lightness + 100.0f * (1.0f - lightness) * (t <= 0.5f ? 2.0f * t : 2.0f * (1.0f - t)));
770         float s = saturation * 5.0f * (t <= 0.5f ? 2.0f * (0.5f - t) : 2.0f * (t - 0.5f));
771         lch.c = lch_chroma(lch.l, s);
772         lch.h = (t <= 0.5f ? hue : hue + divergence);
773         if (lch_to_colormap(lch, colormap + 3 * i))
774             clipped++;
775     }
776     return clipped;
777 }
778
779 int PLDivergingSaturation(int n, unsigned char* colormap,
780         float lightness, float saturation, float hue, float divergence)
781 {
782     triplet lch;
783     lch.l = std::max(0.01f, lightness * 100.0f);
784     int clipped = 0;
785     for (int i = 0; i < n; i++) {
786         float t = i / (n - 1.0f);
787         float s = saturation * 5.0f * (t <= 0.5f ? 2.0f * (0.5f - t) : 2.0f * (t - 0.5f));
788         lch.c = lch_chroma(lch.l, s);
789         lch.h = (t <= 0.5f ? hue : hue + divergence);
790         if (lch_to_colormap(lch, colormap + 3 * i))
791             clipped++;
792     }
793     return clipped;
794 }
795
796 int PLQualitativeHue(int n, unsigned char* colormap,
797         float lightness, float saturation, float hue)
798 {
799     float divergence = twopi * (n - 1.0f) / n;
800     triplet lch;
801     lch.l = std::max(0.01f, lightness * 100.0f);
802     lch.c = lch_chroma(lch.l, saturation * 5.0f);
803     int clipped = 0;
804     for (int i = 0; i < n; i++) {
805         float t = i / (n - 1.0f);
806         lch.h = hue + t * divergence;
807         if (lch_to_colormap(lch, colormap + 3 * i))
808             clipped++;
809     }
810     return clipped;
811 }
812
813 /* CubeHelix */
814
815 int CubeHelix(int n, unsigned char* colormap, float hue,
816         float rot, float saturation, float gamma)
817 {
818     int clipped = 0;
819     for (int i = 0; i < n; i++) {
820         float fract = i / (n - 1.0f);
821         float angle = twopi * (hue / 3.0f + 1.0f + rot * fract);
822         fract = std::pow(fract, gamma);
823         float amp = saturation * fract * (1.0f - fract) / 2.0f;
824         float s = std::sin(angle);
825         float c = std::cos(angle);
826         triplet srgb(
827                 fract + amp * (-0.14861f * c + 1.78277f * s),
828                 fract + amp * (-0.29227f * c - 0.90649f * s),
829                 fract + amp * (1.97294f * c));
830         bool clipped_[3];
831         colormap[3 * i + 0] = float_to_uchar(srgb.r, clipped_ + 0);
832         colormap[3 * i + 1] = float_to_uchar(srgb.g, clipped_ + 1);
833         colormap[3 * i + 2] = float_to_uchar(srgb.b, clipped_ + 2);
834         if (clipped_[0] || clipped_[1] || clipped_[2])
835             clipped++;
836     }
837     return clipped;
838 }
839
840 /* Moreland */
841
842 static triplet lab_to_msh(triplet lab)
843 {
844     triplet msh;
845     msh.m = std::sqrt(lab.l * lab.l + lab.a * lab.a + lab.b * lab.b);
846     msh.s = (msh.m > 0.001f) ? std::acos(lab.l / msh.m) : 0.0f;
847     msh.h = (msh.s > 0.001f) ? std::atan2(lab.b, lab.a) : 0.0f;
848     return msh;
849 }
850
851 static triplet msh_to_lab(triplet msh)
852 {
853     return triplet(
854             msh.m * std::cos(msh.s),
855             msh.m * std::sin(msh.s) * std::cos(msh.h),
856             msh.m * std::sin(msh.s) * std::sin(msh.h));
857 }
858
859 static float adjust_hue(triplet msh, float unsaturated_m)
860 {
861     if (msh.m >= unsaturated_m - 0.1f) {
862         return msh.h;
863     } else {
864         float hue_spin = msh.s * std::sqrt(unsaturated_m * unsaturated_m - msh.m * msh.m)
865             / (msh.m * std::sin(msh.s));
866         if (msh.h > -pi / 3.0f)
867             return msh.h + hue_spin;
868         else
869             return msh.h - hue_spin;
870     }
871 }
872
873 int Moreland(int n, unsigned char* colormap,
874         unsigned char sr0, unsigned char sg0, unsigned char sb0,
875         unsigned char sr1, unsigned char sg1, unsigned char sb1)
876 {
877     triplet omsh0 = lab_to_msh(xyz_to_lab(rgb_to_xyz(srgb_to_rgb(triplet(
878                             uchar_to_float(sr0), uchar_to_float(sg0), uchar_to_float(sb0))))));
879     triplet omsh1 = lab_to_msh(xyz_to_lab(rgb_to_xyz(srgb_to_rgb(triplet(
880                             uchar_to_float(sr1), uchar_to_float(sg1), uchar_to_float(sb1))))));
881     bool place_white = (omsh0.s >= 0.05f && omsh1.s >= 0.05f && hue_diff(omsh0.h, omsh1.h) > pi / 3.0f);
882     float mmid = std::max(std::max(omsh0.m, omsh1.m), 88.0f);
883
884     int clipped = 0;
885     for (int i = 0; i < n; i++) {
886         triplet msh0 = omsh0;
887         triplet msh1 = omsh1;
888         float t = i / (n - 1.0f);
889         if (place_white) {
890             if (t < 0.5f) {
891                 msh1.m = mmid;
892                 msh1.s = 0.0f;
893                 msh1.h = 0.0f;
894                 t *= 2.0f;
895             } else {
896                 msh0.m = mmid;
897                 msh0.s = 0.0f;
898                 msh0.h = 0.0f;
899                 t = 2.0f * t - 1.0f;
900             }
901         }
902         if (msh0.s < 0.05f && msh1.s >= 0.05f) {
903             msh0.h = adjust_hue(msh1, msh0.m);
904         } else if (msh0.s >= 0.05f && msh1.s < 0.05f) {
905             msh1.h = adjust_hue(msh0, msh1.m);
906         }
907         triplet msh = (1.0f - t) * msh0 + t * msh1;
908         if (lab_to_colormap(msh_to_lab(msh), colormap + 3 * i))
909             clipped++;
910     }
911     return clipped;
912 }
913
914 /* McNames */
915
916 static void cart2pol(float x, float y, float* theta, float* rho)
917 {
918     *theta = std::atan2(y, x);
919     *rho = std::hypot(x, y);
920 }
921
922 static void pol2cart(float theta, float rho, float* x, float* y)
923 {
924     *x = rho * std::cos(theta);
925     *y = rho * std::sin(theta);
926 }
927
928 static float windowfunc(float t)
929 {
930     static const float ww = std::sqrt(3.0f / 8.0f);
931     /* triangular window function: */
932 #if 0
933     if (t <= 0.5f) {
934         return ww * 2.0f * t;
935     } else {
936         return ww * 2.0f * (1.0f - t);
937     }
938 #endif
939     /* window function based on cosh: */
940 #if 1
941     static const float acosh2 = std::acosh(2.0f);
942     return 0.95f * ww * (2.0f - std::cosh(acosh2 * (2.0f * t - 1.0f)));
943 #endif
944 }
945
946 int McNames(int n, unsigned char* colormap, float periods)
947 {
948     static const float sqrt3 = std::sqrt(3.0f);
949     static const float a12 = std::asin(1.0f / sqrt3);
950     static const float a23 = pi / 4.0f;
951
952     int clipped = 0;
953     for (int i = 0; i < n; i++) {
954         float t = 1.0f - i / (n - 1.0f);
955         float w = windowfunc(t);
956         float tt = (1.0f - t) * sqrt3;
957         float ttt = (tt - sqrt3 / 2.0f) * periods * twopi / sqrt3;
958
959         float r0, g0, b0, r1, g1, b1, r2, g2, b2;
960         float ag, rd;
961         r0 = tt;
962         g0 = w * std::cos(ttt);
963         b0 = w * std::sin(ttt);
964         cart2pol(r0, g0, &ag, &rd);
965         pol2cart(ag + a12, rd, &r1, &g1);
966         b1 = b0;
967         cart2pol(r1, b1, &ag, &rd);
968         pol2cart(ag + a23, rd, &r2, &b2);
969         g2 = g1;
970
971         bool clipped_[3];
972         colormap[3 * i + 0] = float_to_uchar(r2, clipped_ + 0);
973         colormap[3 * i + 1] = float_to_uchar(g2, clipped_ + 1);
974         colormap[3 * i + 2] = float_to_uchar(b2, clipped_ + 2);
975         if (clipped_[0] || clipped_[1] || clipped_[2])
976             clipped++;
977     }
978     return clipped;
979 }
980
981 }