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