Actual source code: cmap.c

  1: #include <petscsys.h>
  2: #include <petscdraw.h>

  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:  */

 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;

 16: PetscErrorCode PetscDrawUtilitySetGamma(PetscReal g)
 17: {
 18:   PetscFunctionBegin;
 19:   Gamma = g;
 20:   PetscFunctionReturn(PETSC_SUCCESS);
 21: }

 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: }

 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: }

 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: }

 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;

 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: }

 72: static PetscErrorCode PetscDrawCmap_Gray(int mapsize, unsigned char R[], unsigned char G[], unsigned char B[])
 73: {
 74:   int i;
 75:   PetscFunctionBegin;
 76:   for (i = 0; i < mapsize; i++) R[i] = G[i] = B[i] = (unsigned char)((255.0 * i) / (mapsize - 1));
 77:   PetscFunctionReturn(PETSC_SUCCESS);
 78: }

 80: static PetscErrorCode PetscDrawCmap_Jet(int mapsize, unsigned char R[], unsigned char G[], unsigned char B[])
 81: {
 82:   int          i;
 83:   const double knots[] = {0, 1 / 8., 3 / 8., 5 / 8., 7 / 8., 1};

 85:   PetscFunctionBegin;
 86:   for (i = 0; i < mapsize; i++) {
 87:     double u = (double)i / (mapsize - 1);
 88:     double m, r = 0, g = 0, b = 0;
 89:     int    k = 0;
 90:     while (k < 4 && u > knots[k + 1]) k++;
 91:     m = (u - knots[k]) / (knots[k + 1] - knots[k]);
 92:     switch (k) {
 93:     case 0:
 94:       r = 0;
 95:       g = 0;
 96:       b = (m + 1) / 2;
 97:       break;
 98:     case 1:
 99:       r = 0;
100:       g = m;
101:       b = 1;
102:       break;
103:     case 2:
104:       r = m;
105:       g = 1;
106:       b = 1 - m;
107:       break;
108:     case 3:
109:       r = 1;
110:       g = 1 - m;
111:       b = 0;
112:       break;
113:     case 4:
114:       r = 1 - m / 2;
115:       g = 0;
116:       b = 0;
117:       break;
118:     }
119:     R[i] = (unsigned char)(255 * PetscMin(r, 1.0));
120:     G[i] = (unsigned char)(255 * PetscMin(g, 1.0));
121:     B[i] = (unsigned char)(255 * PetscMin(b, 1.0));
122:   }
123:   PetscFunctionReturn(PETSC_SUCCESS);
124: }

126: static PetscErrorCode PetscDrawCmap_Hot(int mapsize, unsigned char R[], unsigned char G[], unsigned char B[])
127: {
128:   int          i;
129:   const double knots[] = {0, 3 / 8., 3 / 4., 1};

131:   PetscFunctionBegin;
132:   for (i = 0; i < mapsize; i++) {
133:     double u = (double)i / (mapsize - 1);
134:     double m, r = 0, g = 0, b = 0;
135:     int    k = 0;
136:     while (k < 2 && u > knots[k + 1]) k++;
137:     m = (u - knots[k]) / (knots[k + 1] - knots[k]);
138:     switch (k) {
139:     case 0:
140:       r = m;
141:       g = 0;
142:       b = 0;
143:       break;
144:     case 1:
145:       r = 1;
146:       g = m;
147:       b = 0;
148:       break;
149:     case 2:
150:       r = 1;
151:       g = 1;
152:       b = m;
153:       break;
154:     }
155:     R[i] = (unsigned char)(255 * PetscMin(r, 1.0));
156:     G[i] = (unsigned char)(255 * PetscMin(g, 1.0));
157:     B[i] = (unsigned char)(255 * PetscMin(b, 1.0));
158:   }
159:   PetscFunctionReturn(PETSC_SUCCESS);
160: }

162: static PetscErrorCode PetscDrawCmap_Bone(int mapsize, unsigned char R[], unsigned char G[], unsigned char B[])
163: {
164:   int i;
165:   PetscFunctionBegin;
166:   (void)PetscDrawCmap_Hot(mapsize, R, G, B);
167:   for (i = 0; i < mapsize; i++) {
168:     double u = (double)i / (mapsize - 1);
169:     double r = (7 * u + B[i] / 255.0) / 8;
170:     double g = (7 * u + G[i] / 255.0) / 8;
171:     double b = (7 * u + R[i] / 255.0) / 8;
172:     R[i]     = (unsigned char)(255 * PetscMin(r, 1.0));
173:     G[i]     = (unsigned char)(255 * PetscMin(g, 1.0));
174:     B[i]     = (unsigned char)(255 * PetscMin(b, 1.0));
175:   }
176:   PetscFunctionReturn(PETSC_SUCCESS);
177: }

179: #include "../src/sys/classes/draw/utils/cmap/coolwarm.h"
180: #include "../src/sys/classes/draw/utils/cmap/parula.h"
181: #include "../src/sys/classes/draw/utils/cmap/viridis.h"
182: #include "../src/sys/classes/draw/utils/cmap/plasma.h"
183: #include "../src/sys/classes/draw/utils/cmap/inferno.h"
184: #include "../src/sys/classes/draw/utils/cmap/magma.h"

186: static struct {
187:   const char *name;
188:   const unsigned char (*data)[3];
189:   PetscErrorCode (*cmap)(int, unsigned char[], unsigned char[], unsigned char[]);
190: } PetscDrawCmapTable[] = {
191:   {"hue",      NULL,                   PetscDrawCmap_Hue }, /* varying hue with constant lightness and saturation */
192:   {"gray",     NULL,                   PetscDrawCmap_Gray}, /* black to white with shades of gray */
193:   {"bone",     NULL,                   PetscDrawCmap_Bone}, /* black to white with gray-blue shades */
194:   {"jet",      NULL,                   PetscDrawCmap_Jet }, /* rainbow-like colormap from NCSA, University of Illinois */
195:   {"hot",      NULL,                   PetscDrawCmap_Hot }, /* black-body radiation */
196:   {"coolwarm", PetscDrawCmap_coolwarm, NULL              }, /* ParaView default (Cool To Warm with Diverging interpolation) */
197:   {"parula",   PetscDrawCmap_parula,   NULL              }, /* MATLAB (default since R2014b) */
198:   {"viridis",  PetscDrawCmap_viridis,  NULL              }, /* matplotlib 1.5 (default since 2.0) */
199:   {"plasma",   PetscDrawCmap_plasma,   NULL              }, /* matplotlib 1.5 */
200:   {"inferno",  PetscDrawCmap_inferno,  NULL              }, /* matplotlib 1.5 */
201:   {"magma",    PetscDrawCmap_magma,    NULL              }, /* matplotlib 1.5 */
202: };

204: PetscErrorCode PetscDrawUtilitySetCmap(const char colormap[], int mapsize, unsigned char R[], unsigned char G[], unsigned char B[])
205: {
206:   int         i, j;
207:   const char *cmap_name_list[PETSC_STATIC_ARRAY_LENGTH(PetscDrawCmapTable)];
208:   PetscInt    id = 0, count = (PetscInt)(sizeof(cmap_name_list) / sizeof(char *));
209:   PetscBool   reverse = PETSC_FALSE, brighten = PETSC_FALSE;
210:   PetscReal   beta = 0;

212:   PetscFunctionBegin;
213:   for (i = 0; i < count; i++) cmap_name_list[i] = PetscDrawCmapTable[i].name;
214:   if (colormap && colormap[0]) {
215:     PetscBool match = PETSC_FALSE;
216:     for (id = 0; !match && id < count; id++) PetscCall(PetscStrcasecmp(colormap, cmap_name_list[id], &match));
217:     PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Colormap '%s' not found", colormap);
218:   }
219:   PetscCall(PetscOptionsGetEList(NULL, NULL, "-draw_cmap", cmap_name_list, count, &id, NULL));
220:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-draw_cmap_reverse", &reverse, NULL));
221:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-draw_cmap_brighten", &beta, &brighten));
222:   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);

224:   if (PetscDrawCmapTable[id].cmap) {
225:     PetscCall(PetscDrawCmapTable[id].cmap(mapsize, R, G, B));
226:   } else {
227:     const unsigned char(*rgb)[3] = PetscDrawCmapTable[id].data;
228:     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);
229:     for (i = 0; i < mapsize; i++) {
230:       R[i] = rgb[i][0];
231:       G[i] = rgb[i][1];
232:       B[i] = rgb[i][2];
233:     }
234:   }

236:   if (reverse) {
237:     i = 0;
238:     j = mapsize - 1;
239:     while (i < j) {
240: #define SWAP(a, i, j) \
241:   do { \
242:     unsigned char t = a[i]; \
243:     a[i]            = a[j]; \
244:     a[j]            = t; \
245:   } while (0)
246:       SWAP(R, i, j);
247:       SWAP(G, i, j);
248:       SWAP(B, i, j);
249: #undef SWAP
250:       i++;
251:       j--;
252:     }
253:   }

255:   if (brighten) {
256:     PetscReal gamma = (beta > 0.0) ? (1 - beta) : (1 / (1 + beta));
257:     for (i = 0; i < mapsize; i++) {
258:       PetscReal r = PetscPowReal((PetscReal)R[i] / 255, gamma);
259:       PetscReal g = PetscPowReal((PetscReal)G[i] / 255, gamma);
260:       PetscReal b = PetscPowReal((PetscReal)B[i] / 255, gamma);
261:       R[i]        = (unsigned char)(255 * PetscMin(r, (PetscReal)1.0));
262:       G[i]        = (unsigned char)(255 * PetscMin(g, (PetscReal)1.0));
263:       B[i]        = (unsigned char)(255 * PetscMin(b, (PetscReal)1.0));
264:     }
265:   }
266:   PetscFunctionReturn(PETSC_SUCCESS);
267: }