Implement CubeHelix color maps, and reorganize command line and GUI accordingly.
[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 /* XYZ and related color spaces helper functions and values */
55
56 static float u_prime(float x, float y, float z)
57 {
58     return 4.0f * x / (x + 15.0f * y + 3.0f * z);
59 }
60
61 static float v_prime(float x, float y, float z)
62 {
63     return 9.0f * y / (x + 15.0f * y + 3.0f * z);
64 }
65
66 static const float d65_x =  95.047f;
67 static const float d65_y = 100.000f;
68 static const float d65_z = 108.883f;
69 static const float d65_u_prime = u_prime(d65_x, d65_y, d65_z);
70 static const float d65_v_prime = v_prime(d65_x, d65_y, d65_z);
71
72 /* Color space conversion: LCH <-> LUV */
73
74 static float lch_saturation(float l, float c)
75 {
76     return c / std::max(l, 1e-8f);
77 }
78
79 static float lch_chroma(float l, float s)
80 {
81     return s * l;
82 }
83
84 static void lch_to_luv(float c, float h, float* u, float* v)
85 {
86     *u = c * std::cos(h);
87     *v = c * std::sin(h);
88 }
89
90 static void luv_to_lch(float u, float v, float* c, float* h)
91 {
92     *c = std::hypot(u, v);
93     *h = std::atan2(v, u);
94     if (*h < 0.0f)
95         *h += twopi;
96 }
97
98 static float luv_saturation(float l, float u, float v)
99 {
100     return lch_saturation(l, std::hypot(u, v));
101 }
102
103 /* Color space conversion: LUV <-> XYZ */
104
105 static void luv_to_xyz(float l, float u, float v, float* x, float* y, float* z)
106 {
107     float u_prime = u / (13.0f * l) + d65_u_prime;
108     float v_prime = v / (13.0f * l) + d65_v_prime;
109     if (l <= 8.0f) {
110         *y = d65_y * l * (3.0f * 3.0f * 3.0f / (29.0f * 29.0f * 29.0f));
111     } else {
112         float tmp = (l + 16.0f) / 116.0f;
113         *y = d65_y * tmp * tmp * tmp;
114     }
115     *x = (*y) * (9.0f * u_prime) / (4.0f * v_prime);
116     *z = (*y) * (12.0f - 3.0f * u_prime - 20.0f * v_prime) / (4.0f * v_prime);
117 }
118
119 static void xyz_to_luv(float x, float y, float z, float* l, float* u, float* v)
120 {
121     float y_ratio = y / d65_y;
122     if (y_ratio <= (6.0f * 6.0f * 6.0f) / (29.0f * 29.0f * 29.0f)) {
123         *l = (29.0f * 29.0f * 29.0f) / (3.0f * 3.0f * 3.0f) * y_ratio;
124     } else {
125         *l = 116.0f * std::cbrt(y_ratio) - 16.0f;
126     }
127     *u = 13.0f * (*l) * (u_prime(x, y, z) - d65_u_prime);
128     *v = 13.0f * (*l) * (v_prime(x, y, z) - d65_v_prime);
129 }
130
131 /* Color space conversion: RGB <-> XYZ */
132
133 static void rgb_to_xyz(float r, float g, float b, float* x, float* y, float *z)
134 {
135     *x = (0.4124f * r + 0.3576f * g + 0.1805f * b) * 100.0f;
136     *y = (0.2126f * r + 0.7152f * g + 0.0722f * b) * 100.0f;
137     *z = (0.0193f * r + 0.1192f * g + 0.9505f * b) * 100.0f;
138 }
139
140 static void xyz_to_rgb(float x, float y, float z, float* r, float* g, float* b)
141 {
142     *r = clamp((+3.2406255f * x - 1.5372080f * y - 0.4986286f * z) / 100.0f, 0.0f, 1.0f);
143     *g = clamp((-0.9689307f * x + 1.8757561f * y + 0.0415175f * z) / 100.0f, 0.0f, 1.0f);
144     *b = clamp((+0.0557101f * x - 0.2040211f * y + 1.0569959f * z) / 100.0f, 0.0f, 1.0f);
145 }
146
147 /* Color space conversion: RGB <-> sRGB */
148
149 static float rgb_to_srgb_helper(float x)
150 {
151     return (x <= 0.0031308f ? (x * 12.92f) : (1.055f * std::pow(x, 1.0f / 2.4f) - 0.055f));
152 }
153
154 static void rgb_to_srgb(float r, float g, float b, float* sr, float* sg, float* sb)
155 {
156     *sr = rgb_to_srgb_helper(r);
157     *sg = rgb_to_srgb_helper(g);
158     *sb = rgb_to_srgb_helper(b);
159 }
160
161 static float srgb_to_rgb_helper(float x)
162 {
163     return (x <= 0.04045f ? (x / 12.92f) : std::pow((x + 0.055f) / 1.055f, 2.4f));
164 }
165
166 static void srgb_to_rgb(float sr, float sg, float sb, float* r, float* g, float* b)
167 {
168     *r = srgb_to_rgb_helper(sr);
169     *g = srgb_to_rgb_helper(sg);
170     *b = srgb_to_rgb_helper(sb);
171 }
172
173 /* Helpers for LUV colors */
174
175 typedef struct {
176     float l;
177     float u;
178     float v;
179 } LUVColor;
180
181 LUVColor operator+(LUVColor a, LUVColor b)
182 {
183     LUVColor c = { .l = a.l + b.l, .u = a.u + b.u, .v = a.v + b.v };
184     return c;
185 }
186
187 LUVColor operator*(float a, LUVColor b)
188 {
189     LUVColor c = { .l = a * b.l, .u = a * b.u, .v = a * b.v };
190     return c;
191 }
192
193 static float srgb_to_lch_hue(float sr, float sg, float sb)
194 {
195     float r, g, b;
196     srgb_to_rgb(sr, sg, sb, &r, &g, &b);
197     float x, y, z;
198     rgb_to_xyz(r, g, b, &x, &y, &z);
199     float l, u, v;
200     xyz_to_luv(x, y, z, &l, &u, &v);
201     float c, h;
202     luv_to_lch(u, v, &c, &h);
203     return h;
204 }
205
206 // Compute most saturated color that fits into the sRGB
207 // cube for the given LCH hue value. This is the core
208 // of the paper.
209 static LUVColor most_saturated_in_srgb(float hue)
210 {
211     /* Static values, only computed once */
212     static float h[] = {
213         srgb_to_lch_hue(1, 0, 0),
214         srgb_to_lch_hue(1, 1, 0),
215         srgb_to_lch_hue(0, 1, 0),
216         srgb_to_lch_hue(0, 1, 1),
217         srgb_to_lch_hue(0, 0, 1),
218         srgb_to_lch_hue(1, 0, 1)
219     };
220
221     /* RGB values and variable pointers to them */
222     int i, j, k;
223     if (hue < h[0]) {
224         i = 2;
225         j = 1;
226         k = 0;
227     } else if (hue < h[1]) {
228         i = 1;
229         j = 2;
230         k = 0;
231     } else if (hue < h[2]) {
232         i = 0;
233         j = 2;
234         k = 1;
235     } else if (hue < h[3]) {
236         i = 2;
237         j = 0;
238         k = 1;
239     } else if (hue < h[4]) {
240         i = 1;
241         j = 0;
242         k = 2;
243     } else if (hue < h[5]) {
244         i = 0;
245         j = 1;
246         k = 2;
247     } else {
248         i = 2;
249         j = 1;
250         k = 0;
251     }
252
253     /* Compute the missing component */
254     float srgb[3];
255     float M[3][3] = {
256         { 0.4124f, 0.3576f, 0.1805f },
257         { 0.2126f, 0.7152f, 0.0722f },
258         { 0.0193f, 0.1192f, 0.9505f }
259     };
260     float alpha = -std::sin(hue);
261     float beta = std::cos(hue);
262     float T = alpha * d65_u_prime + beta * d65_v_prime;
263     srgb[j] = 0.0f;
264     srgb[k] = 1.0f;
265     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]);
266     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]);
267     srgb[i] = rgb_to_srgb_helper(clamp(- q0 / q1, 0.0f, 1.0f));
268
269     /* Convert back to LUV */
270     float r, g, b;
271     srgb_to_rgb(srgb[0], srgb[1], srgb[2], &r, &g, &b);
272     float x, y, z;
273     rgb_to_xyz(r, g, b, &x, &y, &z);
274     float l, u, v;
275     xyz_to_luv(x, y, z, &l, &u, &v);
276     LUVColor color = { .l = l, .u = u, .v = v };
277     return color;
278 }
279
280 static float Smax(float l, float h)
281 {
282     LUVColor pmid = most_saturated_in_srgb(h);
283     LUVColor pend = { .l = 0.0f, .u = 0.0f, .v = 0.0f };
284     if (l > pmid.l)
285         pend.l = 100.0f;
286     float alpha = (pend.l - l) / (pend.l - pmid.l);
287     float pmids = luv_saturation(pmid.l, pmid.u, pmid.v);
288     float pends = luv_saturation(pend.l, pend.u, pend.v);
289     return alpha * (pmids - pends) + pends;
290 }
291
292 static LUVColor get_bright_point()
293 {
294     static LUVColor pb = { .l = -1.0f, .u = -1.0f, .v = -1.0f };
295     if (pb.l < 0.0f) {
296         float x, y, z;
297         rgb_to_xyz(1, 1, 0, &x, &y, &z);
298         float l, u, v;
299         xyz_to_luv(x, y, z, &l, &u, &v);
300         pb.l = l;
301         pb.u = u;
302         pb.v = v;
303     }
304     return pb;
305 }
306
307 static float mix_hue(float alpha, float h0, float h1)
308 {
309     float M = std::fmod(pi + h1 - h0, twopi) - pi;
310     return std::fmod(h0 + alpha * M, twopi);
311 }
312
313 static void get_color_points(float hue, float saturation, float warmth,
314         LUVColor pb, float pb_hue, float pb_saturation,
315         LUVColor* p0, LUVColor* p1, LUVColor* p2,
316         LUVColor* q0, LUVColor* q1, LUVColor* q2)
317 {
318     p0->l = 0.0f;
319     lch_to_luv(0.0f, hue, &(p0->u), &(p0->v));
320     *p1 = most_saturated_in_srgb(hue);
321     float p2l = (1.0f - warmth) * 100.0f + warmth * pb.l;
322     float p2h = mix_hue(warmth, hue, pb_hue);
323     float p2c = lch_chroma(p2l, std::min(Smax(p2l, p2h), warmth * saturation * pb_saturation));
324     p2->l = p2l;
325     lch_to_luv(p2c, p2h, &(p2->u), &(p2->v));
326     *q0 = (1.0f - saturation) * (*p0) + saturation * (*p1);
327     *q2 = (1.0f - saturation) * (*p2) + saturation * (*p1);
328     *q1 = 0.5f * ((*q0) + (*q2));
329 }
330
331 static LUVColor B(LUVColor b0, LUVColor b1, LUVColor b2, float t)
332 {
333     float a = (1.0f - t) * (1.0f - t);
334     float b = 2.0f * (1.0f - t) * t;
335     float c = t * t;
336     LUVColor color = a * b0 + b * b1 + c * b2;
337     return color;
338 }
339
340 static float invB(float b0, float b1, float b2, float v)
341 {
342     return (b0 - b1 + std::sqrt(std::max(b1 * b1 - b0 * b2 + (b0 - 2.0f * b1 + b2) * v, 0.0f)))
343         / (b0 - 2.0f * b1 + b2);
344 }
345
346 static LUVColor get_colormap_entry(float t,
347         LUVColor p0, LUVColor p2,
348         LUVColor q0, LUVColor q1, LUVColor q2,
349         float contrast, float brightness)
350 {
351     float l = 125 - 125 * std::pow(0.2f, (1.0f - contrast) * brightness + t * contrast);
352     float T = (l <= q1.l ? 0.5f * invB(p0.l, q0.l, q1.l, l) : 0.5f * invB(q1.l, q2.l, p2.l, l) + 0.5f);
353     return (T <= 0.5f ? B(p0, q0, q1, 2.0f * T) : B(q1, q2, p2, 2.0f * (T - 0.5f)));
354 }
355
356 static void convert_colormap_entry(LUVColor color, unsigned char* srgb)
357 {
358     float x, y, z;
359     luv_to_xyz(color.l, color.u, color.v, &x, &y, &z);
360     float r, g, b;
361     xyz_to_rgb(x, y, z, &r, &g, &b);
362     float sr, sg, sb;
363     rgb_to_srgb(r, g, b, &sr, &sg, &sb);
364     srgb[0] = std::round(sr * 255.0f);
365     srgb[1] = std::round(sg * 255.0f);
366     srgb[2] = std::round(sb * 255.0f);
367 }
368
369 /* Public functions: Brewer-like color maps */
370
371 float BrewerSequentialDefaultContrastForSmallN(int n)
372 {
373     return std::min(0.88f, 0.34f + 0.06f * n);
374 }
375
376 void BrewerSequential(int n, unsigned char* colormap, float hue,
377         float contrast, float saturation, float brightness, float warmth)
378 {
379     LUVColor pb, p0, p1, p2, q0, q1, q2;
380     pb = get_bright_point();
381     float pbc, pbh, pbs;
382     luv_to_lch(pb.u, pb.v, &pbc, &pbh);
383     pbs = lch_saturation(pb.l, pbc);
384     get_color_points(hue, saturation, warmth, pb, pbh, pbs, &p0, &p1, &p2, &q0, &q1, &q2);
385
386     for (int i = 0; i < n; i++) {
387         float t = (n - 1 - i) / (n - 1.0f);
388         LUVColor c = get_colormap_entry(t, p0, p2, q0, q1, q2, contrast, brightness);
389         convert_colormap_entry(c, colormap + 3 * i);
390     }
391 }
392
393 float BrewerDivergingDefaultContrastForSmallN(int n)
394 {
395     return std::min(0.88f, 0.34f + 0.06f * n);
396 }
397
398 void BrewerDiverging(int n, unsigned char* colormap, float hue, float divergence,
399         float contrast, float saturation, float brightness, float warmth)
400 {
401     float hue1 = hue + divergence;
402     if (hue1 >= twopi)
403         hue1 -= twopi;
404
405     LUVColor pb;
406     LUVColor p00, p01, p02, q00, q01, q02;
407     LUVColor p10, p11, p12, q10, q11, q12;
408     pb = get_bright_point();
409     float pbc, pbh, pbs;
410     luv_to_lch(pb.u, pb.v, &pbc, &pbh);
411     pbs = lch_saturation(pb.l, pbc);
412     get_color_points(hue,  saturation, warmth, pb, pbh, pbs, &p00, &p01, &p02, &q00, &q01, &q02);
413     get_color_points(hue1, saturation, warmth, pb, pbh, pbs, &p10, &p11, &p12, &q10, &q11, &q12);
414
415     for (int i = 0; i < n; i++) {
416         LUVColor c;
417         if (n % 2 == 1 && i == n / 2) {
418             // compute neutral color in the middle of the map
419             LUVColor c0 = get_colormap_entry(1.0f, p00, p02, q00, q01, q02, contrast, brightness);
420             LUVColor c1 = get_colormap_entry(1.0f, p10, p12, q10, q11, q12, contrast, brightness);
421             if (n <= 9) {
422                 // for discrete color maps, use an extra neutral color
423                 float c0s = luv_saturation(c0.l, c0.u, c0.v);
424                 float c1s = luv_saturation(c1.l, c1.u, c1.v);
425                 float sn = 0.5f * (c0s + c1s) * warmth;
426                 c.l = 0.5f * (c0.l + c1.l);
427                 float cc = lch_chroma(c.l, std::min(Smax(c.l, pbh), sn));
428                 lch_to_luv(cc, pbh, &(c.u), &(c.v));
429             } else {
430                 // for continuous color maps, use an average, since the extra neutral color looks bad
431                 c = 0.5f * (c0 + c1);
432             }
433         } else {
434             float t = i / (n - 1.0f);
435             if (i < n / 2) {
436                 float tt = 2.0f * t;
437                 c = get_colormap_entry(tt, p00, p02, q00, q01, q02, contrast, brightness);
438             } else {
439                 float tt = 2.0f * (1.0f - t);
440                 c = get_colormap_entry(tt, p10, p12, q10, q11, q12, contrast, brightness);
441             }
442         }
443         convert_colormap_entry(c, colormap + 3 * i);
444     }
445 }
446
447 static float HueDiff(float h0, float h1)
448 {
449     float t = std::fabs(h1 - h0);
450     return (t < pi ? t : twopi - t);
451 }
452
453 void BrewerQualitative(int n, unsigned char* colormap, float hue, float divergence,
454         float contrast, float saturation, float brightness)
455 {
456     // Get all information about yellow
457     static float yl = -1.0f, yh = -1.0f;
458     if (yl < 0.0f) {
459         float yx, yy, yz;
460         rgb_to_xyz(1, 1, 0, &yx, &yy, &yz);
461         float yu, yv;
462         xyz_to_luv(yx, yy, yz, &yl, &yu, &yv);
463         float yc;
464         luv_to_lch(yu, yv, &yc, &yh);
465     }
466
467     // Get saturation of red (maximum possible saturation)
468     static float rs = -1.0f;
469     if (rs < 0.0f) {
470         float rx, ry, rz;
471         rgb_to_xyz(1, 0, 0, &rx, &ry, &rz);
472         float rl, ru, rv;
473         xyz_to_luv(rx, ry, rz, &rl, &ru, &rv);
474         rs = luv_saturation(rl, ru, rv);
475     }
476
477     // Derive parameters of the method
478     float eps = hue / twopi;
479     float r = divergence / twopi;
480     float l0 = brightness * yl;
481     float l1 = (1.0f - contrast) * l0;
482
483     // Generate colors
484     for (int i = 0; i < n; i++) {
485         float t = i / (n - 1.0f);
486         float ch = std::fmod(twopi * (eps + t * r), twopi);
487         float alpha = HueDiff(ch, yh) / pi;
488         float cl = (1.0f - alpha) * l0 + alpha * l1;
489         float cs = std::min(Smax(cl, ch), saturation * rs);
490         LUVColor c;
491         c.l = cl;
492         lch_to_luv(lch_chroma(cl, cs), ch, &(c.u), &(c.v));
493         convert_colormap_entry(c, colormap + 3 * i);
494     }
495 }
496
497 /* Public functions: CubeHelix */
498
499 int CubeHelix(int n, unsigned char* colormap, float hue,
500         float rot, float saturation, float gamma)
501 {
502     int clippings = 0;
503     for (int i = 0; i < n; i++) {
504         float fract = i / (n - 1.0f);
505         float angle = twopi * (hue / 3.0f + 1.0f + rot * fract);
506         fract = std::pow(fract, gamma);
507         float amp = saturation * fract * (1.0f - fract) / 2.0f;
508         float s = std::sin(angle);
509         float c = std::cos(angle);
510         float r = fract + amp * (-0.14861f * c + 1.78277f * s);
511         float g = fract + amp * (-0.29227f * c - 0.90649f * s);
512         float b = fract + amp * (1.97294f * c);
513         bool clipped = false;
514         if (r < 0.0f) {
515             r = 0.0f;
516             clipped = true;
517         } else if (r > 1.0f) {
518             r = 1.0f;
519             clipped = true;
520         }
521         if (g < 0.0f) {
522             g = 0.0f;
523             clipped = true;
524         } else if (g > 1.0f) {
525             g = 1.0f;
526             clipped = true;
527         }
528         if (b < 0.0f) {
529             b = 0.0f;
530             clipped = true;
531         } else if (b > 1.0f) {
532             b = 1.0f;
533             clipped = true;
534         }
535         if (clipped)
536             clippings++;
537         colormap[3 * i + 0] = r * 255.0f;
538         colormap[3 * i + 1] = g * 255.0f;
539         colormap[3 * i + 2] = b * 255.0f;
540     }
541     return clippings;
542 }
543
544 }