xref: /petsc/src/sys/classes/draw/utils/cmap.c (revision 83c9f7f220f3448cfcd2a76a64462ac242b37309)
1 #include <petscsys.h> /*I "petscsys.h" I*/
2 #include <petscdraw.h>
3 
4 /*
5     Set up a color map, using uniform separation in hue space.
6     Map entries are Red, Green, Blue.
7     Values are "gamma" corrected.
8  */
9 
10 /*
11    Gamma is a monitor dependent value.  The value here is an
12    approximate that gives somewhat better results than Gamma = 1.
13  */
14 static PetscReal Gamma = 2.0;
15 
16 PetscErrorCode PetscDrawUtilitySetGamma(PetscReal g) {
17   PetscFunctionBegin;
18   Gamma = g;
19   PetscFunctionReturn(0);
20 }
21 
22 static inline double PetscHlsHelper(double m1, double m2, double h) {
23   while (h > 1.0) h -= 1.0;
24   while (h < 0.0) h += 1.0;
25   if (h < 1 / 6.0) return m1 + (m2 - m1) * h * 6;
26   if (h < 1 / 2.0) return m2;
27   if (h < 2 / 3.0) return m1 + (m2 - m1) * (2 / 3.0 - h) * 6;
28   return m1;
29 }
30 
31 static inline void PetscHlsToRgb(double h, double l, double s, double *r, double *g, double *b) {
32   if (s > 0.0) {
33     double m2 = l <= 0.5 ? l * (1.0 + s) : l + s - (l * s);
34     double m1 = 2 * l - m2;
35     *r        = PetscHlsHelper(m1, m2, h + 1 / 3.);
36     *g        = PetscHlsHelper(m1, m2, h);
37     *b        = PetscHlsHelper(m1, m2, h - 1 / 3.);
38   } else {
39     /* ignore hue */
40     *r = *g = *b = l;
41   }
42 }
43 
44 static inline void PetscGammaCorrect(double *r, double *g, double *b) {
45   PetscReal igamma = 1 / Gamma;
46   *r               = (double)PetscPowReal((PetscReal)*r, igamma);
47   *g               = (double)PetscPowReal((PetscReal)*g, igamma);
48   *b               = (double)PetscPowReal((PetscReal)*b, igamma);
49 }
50 
51 static PetscErrorCode PetscDrawCmap_Hue(int mapsize, unsigned char R[], unsigned char G[], unsigned char B[]) {
52   int    i;
53   double maxhue = 212.0 / 360, lightness = 0.5, saturation = 1.0;
54 
55   PetscFunctionBegin;
56   for (i = 0; i < mapsize; i++) {
57     double hue = maxhue * (double)i / (mapsize - 1), r, g, b;
58     PetscHlsToRgb(hue, lightness, saturation, &r, &g, &b);
59     PetscGammaCorrect(&r, &g, &b);
60     R[i] = (unsigned char)(255 * PetscMin(r, 1.0));
61     G[i] = (unsigned char)(255 * PetscMin(g, 1.0));
62     B[i] = (unsigned char)(255 * PetscMin(b, 1.0));
63   }
64   PetscFunctionReturn(0);
65 }
66 
67 static PetscErrorCode PetscDrawCmap_Gray(int mapsize, unsigned char R[], unsigned char G[], unsigned char B[]) {
68   int i;
69   PetscFunctionBegin;
70   for (i = 0; i < mapsize; i++) R[i] = G[i] = B[i] = (unsigned char)((255.0 * i) / (mapsize - 1));
71   PetscFunctionReturn(0);
72 }
73 
74 static PetscErrorCode PetscDrawCmap_Jet(int mapsize, unsigned char R[], unsigned char G[], unsigned char B[]) {
75   int          i;
76   const double knots[] = {0, 1 / 8., 3 / 8., 5 / 8., 7 / 8., 1};
77 
78   PetscFunctionBegin;
79   for (i = 0; i < mapsize; i++) {
80     double u = (double)i / (mapsize - 1);
81     double m, r = 0, g = 0, b = 0;
82     int    k = 0;
83     while (k < 4 && u > knots[k + 1]) k++;
84     m = (u - knots[k]) / (knots[k + 1] - knots[k]);
85     switch (k) {
86     case 0:
87       r = 0;
88       g = 0;
89       b = (m + 1) / 2;
90       break;
91     case 1:
92       r = 0;
93       g = m;
94       b = 1;
95       break;
96     case 2:
97       r = m;
98       g = 1;
99       b = 1 - m;
100       break;
101     case 3:
102       r = 1;
103       g = 1 - m;
104       b = 0;
105       break;
106     case 4:
107       r = 1 - m / 2;
108       g = 0;
109       b = 0;
110       break;
111     }
112     R[i] = (unsigned char)(255 * PetscMin(r, 1.0));
113     G[i] = (unsigned char)(255 * PetscMin(g, 1.0));
114     B[i] = (unsigned char)(255 * PetscMin(b, 1.0));
115   }
116   PetscFunctionReturn(0);
117 }
118 
119 static PetscErrorCode PetscDrawCmap_Hot(int mapsize, unsigned char R[], unsigned char G[], unsigned char B[]) {
120   int          i;
121   const double knots[] = {0, 3 / 8., 3 / 4., 1};
122 
123   PetscFunctionBegin;
124   for (i = 0; i < mapsize; i++) {
125     double u = (double)i / (mapsize - 1);
126     double m, r = 0, g = 0, b = 0;
127     int    k = 0;
128     while (k < 2 && u > knots[k + 1]) k++;
129     m = (u - knots[k]) / (knots[k + 1] - knots[k]);
130     switch (k) {
131     case 0:
132       r = m;
133       g = 0;
134       b = 0;
135       break;
136     case 1:
137       r = 1;
138       g = m;
139       b = 0;
140       break;
141     case 2:
142       r = 1;
143       g = 1;
144       b = m;
145       break;
146     }
147     R[i] = (unsigned char)(255 * PetscMin(r, 1.0));
148     G[i] = (unsigned char)(255 * PetscMin(g, 1.0));
149     B[i] = (unsigned char)(255 * PetscMin(b, 1.0));
150   }
151   PetscFunctionReturn(0);
152 }
153 
154 static PetscErrorCode PetscDrawCmap_Bone(int mapsize, unsigned char R[], unsigned char G[], unsigned char B[]) {
155   int i;
156   PetscFunctionBegin;
157   (void)PetscDrawCmap_Hot(mapsize, R, G, B);
158   for (i = 0; i < mapsize; i++) {
159     double u = (double)i / (mapsize - 1);
160     double r = (7 * u + B[i] / 255.0) / 8;
161     double g = (7 * u + G[i] / 255.0) / 8;
162     double b = (7 * u + R[i] / 255.0) / 8;
163     R[i]     = (unsigned char)(255 * PetscMin(r, 1.0));
164     G[i]     = (unsigned char)(255 * PetscMin(g, 1.0));
165     B[i]     = (unsigned char)(255 * PetscMin(b, 1.0));
166   }
167   PetscFunctionReturn(0);
168 }
169 
170 #include "../src/sys/classes/draw/utils/cmap/coolwarm.h"
171 #include "../src/sys/classes/draw/utils/cmap/parula.h"
172 #include "../src/sys/classes/draw/utils/cmap/viridis.h"
173 #include "../src/sys/classes/draw/utils/cmap/plasma.h"
174 #include "../src/sys/classes/draw/utils/cmap/inferno.h"
175 #include "../src/sys/classes/draw/utils/cmap/magma.h"
176 
177 static struct {
178   const char *name;
179   const unsigned char (*data)[3];
180   PetscErrorCode (*cmap)(int, unsigned char[], unsigned char[], unsigned char[]);
181 } PetscDrawCmapTable[] = {
182   {"hue",      NULL,                   PetscDrawCmap_Hue }, /* varying hue with constant lightness and saturation */
183   {"gray",     NULL,                   PetscDrawCmap_Gray}, /* black to white with shades of gray */
184   {"bone",     NULL,                   PetscDrawCmap_Bone}, /* black to white with gray-blue shades */
185   {"jet",      NULL,                   PetscDrawCmap_Jet }, /* rainbow-like colormap from NCSA, University of Illinois */
186   {"hot",      NULL,                   PetscDrawCmap_Hot }, /* black-body radiation */
187   {"coolwarm", PetscDrawCmap_coolwarm, NULL              }, /* ParaView default (Cool To Warm with Diverging interpolation) */
188   {"parula",   PetscDrawCmap_parula,   NULL              }, /* MATLAB (default since R2014b) */
189   {"viridis",  PetscDrawCmap_viridis,  NULL              }, /* matplotlib 1.5 (default since 2.0) */
190   {"plasma",   PetscDrawCmap_plasma,   NULL              }, /* matplotlib 1.5 */
191   {"inferno",  PetscDrawCmap_inferno,  NULL              }, /* matplotlib 1.5 */
192   {"magma",    PetscDrawCmap_magma,    NULL              }, /* matplotlib 1.5 */
193 };
194 
195 PetscErrorCode PetscDrawUtilitySetCmap(const char colormap[], int mapsize, unsigned char R[], unsigned char G[], unsigned char B[]) {
196   int         i, j;
197   const char *cmap_name_list[PETSC_STATIC_ARRAY_LENGTH(PetscDrawCmapTable)];
198   PetscInt    id = 0, count = (PetscInt)(sizeof(cmap_name_list) / sizeof(char *));
199   PetscBool   reverse = PETSC_FALSE, brighten = PETSC_FALSE;
200   PetscReal   beta = 0;
201 
202   PetscFunctionBegin;
203   for (i = 0; i < count; i++) cmap_name_list[i] = PetscDrawCmapTable[i].name;
204   if (colormap && colormap[0]) {
205     PetscBool match = PETSC_FALSE;
206     for (id = 0; !match && id < count; id++) PetscCall(PetscStrcasecmp(colormap, cmap_name_list[id], &match));
207     PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Colormap '%s' not found", colormap);
208   }
209   PetscCall(PetscOptionsGetEList(NULL, NULL, "-draw_cmap", cmap_name_list, count, &id, NULL));
210   PetscCall(PetscOptionsGetBool(NULL, NULL, "-draw_cmap_reverse", &reverse, NULL));
211   PetscCall(PetscOptionsGetReal(NULL, NULL, "-draw_cmap_brighten", &beta, &brighten));
212   PetscCheck(!brighten || (beta > (PetscReal)-1 && beta < (PetscReal) + 1), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "brighten parameter %g must be in the range (-1,1)", (double)beta);
213 
214   if (PetscDrawCmapTable[id].cmap) {
215     PetscCall(PetscDrawCmapTable[id].cmap(mapsize, R, G, B));
216   } else {
217     const unsigned char(*rgb)[3] = PetscDrawCmapTable[id].data;
218     PetscCheck(mapsize == 256 - PETSC_DRAW_BASIC_COLORS, PETSC_COMM_SELF, PETSC_ERR_SUP, "Colormap '%s' with size %d not supported", cmap_name_list[id], mapsize);
219     for (i = 0; i < mapsize; i++) {
220       R[i] = rgb[i][0];
221       G[i] = rgb[i][1];
222       B[i] = rgb[i][2];
223     }
224   }
225 
226   if (reverse) {
227     i = 0;
228     j = mapsize - 1;
229     while (i < j) {
230 #define SWAP(a, i, j) \
231   do { \
232     unsigned char t = a[i]; \
233     a[i]            = a[j]; \
234     a[j]            = t; \
235   } while (0)
236       SWAP(R, i, j);
237       SWAP(G, i, j);
238       SWAP(B, i, j);
239 #undef SWAP
240       i++;
241       j--;
242     }
243   }
244 
245   if (brighten) {
246     PetscReal gamma = (beta > 0.0) ? (1 - beta) : (1 / (1 + beta));
247     for (i = 0; i < mapsize; i++) {
248       PetscReal r = PetscPowReal((PetscReal)R[i] / 255, gamma);
249       PetscReal g = PetscPowReal((PetscReal)G[i] / 255, gamma);
250       PetscReal b = PetscPowReal((PetscReal)B[i] / 255, gamma);
251       R[i]        = (unsigned char)(255 * PetscMin(r, (PetscReal)1.0));
252       G[i]        = (unsigned char)(255 * PetscMin(g, (PetscReal)1.0));
253       B[i]        = (unsigned char)(255 * PetscMin(b, (PetscReal)1.0));
254     }
255   }
256   PetscFunctionReturn(0);
257 }
258