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(0);
21 }
22 
PetscHlsHelper(double m1,double m2,double h)23 PETSC_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 PETSC_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 PETSC_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(0);
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   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(0);
78 }
79 
PetscDrawCmap_Jet(int mapsize,unsigned char R[],unsigned char G[],unsigned char B[])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};
84 
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; 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: r = 0;     g = 0;   b = (m+1)/2; break;
93     case 1: r = 0;     g = m;   b = 1;       break;
94     case 2: r = m;     g = 1;   b = 1-m;     break;
95     case 3: r = 1;     g = 1-m; b = 0;       break;
96     case 4: r = 1-m/2; g = 0;   b = 0;       break;
97     }
98     R[i] = (unsigned char)(255*PetscMin(r,1.0));
99     G[i] = (unsigned char)(255*PetscMin(g,1.0));
100     B[i] = (unsigned char)(255*PetscMin(b,1.0));
101   }
102   PetscFunctionReturn(0);
103 }
104 
PetscDrawCmap_Hot(int mapsize,unsigned char R[],unsigned char G[],unsigned char B[])105 static PetscErrorCode PetscDrawCmap_Hot(int mapsize,unsigned char R[],unsigned char G[],unsigned char B[])
106 {
107   int          i;
108   const double knots[] =  {0, 3/8., 3/4., 1};
109 
110   PetscFunctionBegin;
111   for (i=0; i<mapsize; i++) {
112     double u = (double)i/(mapsize-1);
113     double m, r=0, g=0, b=0; int k = 0;
114     while (k < 2 && u > knots[k+1]) k++;
115     m = (u-knots[k])/(knots[k+1]-knots[k]);
116     switch(k) {
117     case 0: r = m; g = 0; b = 0; break;
118     case 1: r = 1; g = m; b = 0; break;
119     case 2: r = 1; g = 1; b = m; break;
120     }
121     R[i] = (unsigned char)(255*PetscMin(r,1.0));
122     G[i] = (unsigned char)(255*PetscMin(g,1.0));
123     B[i] = (unsigned char)(255*PetscMin(b,1.0));
124   }
125   PetscFunctionReturn(0);
126 }
127 
PetscDrawCmap_Bone(int mapsize,unsigned char R[],unsigned char G[],unsigned char B[])128 static PetscErrorCode PetscDrawCmap_Bone(int mapsize,unsigned char R[],unsigned char G[],unsigned char B[])
129 {
130   int i;
131   PetscFunctionBegin;
132   (void)PetscDrawCmap_Hot(mapsize,R,G,B);
133   for (i=0; i<mapsize; i++) {
134     double u = (double)i/(mapsize-1);
135     double r = (7*u + B[i]/255.0)/8;
136     double g = (7*u + G[i]/255.0)/8;
137     double b = (7*u + R[i]/255.0)/8;
138     R[i] = (unsigned char)(255*PetscMin(r,1.0));
139     G[i] = (unsigned char)(255*PetscMin(g,1.0));
140     B[i] = (unsigned char)(255*PetscMin(b,1.0));
141   }
142   PetscFunctionReturn(0);
143 }
144 
145 #include "../src/sys/classes/draw/utils/cmap/coolwarm.h"
146 #include "../src/sys/classes/draw/utils/cmap/parula.h"
147 #include "../src/sys/classes/draw/utils/cmap/viridis.h"
148 #include "../src/sys/classes/draw/utils/cmap/plasma.h"
149 #include "../src/sys/classes/draw/utils/cmap/inferno.h"
150 #include "../src/sys/classes/draw/utils/cmap/magma.h"
151 
152 static struct {
153   const char           *name;
154   const unsigned char (*data)[3];
155   PetscErrorCode      (*cmap)(int,unsigned char[],unsigned char[],unsigned char[]);
156 } PetscDrawCmapTable[] = {
157   {"hue",      NULL, PetscDrawCmap_Hue },     /* varying hue with constant lightness and saturation */
158   {"gray",     NULL, PetscDrawCmap_Gray},     /* black to white with shades of gray */
159   {"bone",     NULL, PetscDrawCmap_Bone},     /* black to white with gray-blue shades */
160   {"jet",      NULL, PetscDrawCmap_Jet },     /* rainbow-like colormap from NCSA, University of Illinois */
161   {"hot",      NULL, PetscDrawCmap_Hot },     /* black-body radiation */
162   {"coolwarm", PetscDrawCmap_coolwarm, NULL}, /* ParaView default (Cool To Warm with Diverging interpolation) */
163   {"parula",   PetscDrawCmap_parula,   NULL}, /* MATLAB (default since R2014b) */
164   {"viridis",  PetscDrawCmap_viridis,  NULL}, /* matplotlib 1.5 (default since 2.0) */
165   {"plasma",   PetscDrawCmap_plasma,   NULL}, /* matplotlib 1.5 */
166   {"inferno",  PetscDrawCmap_inferno,  NULL}, /* matplotlib 1.5 */
167   {"magma",    PetscDrawCmap_magma,    NULL}, /* matplotlib 1.5 */
168 };
169 
PetscDrawUtilitySetCmap(const char colormap[],int mapsize,unsigned char R[],unsigned char G[],unsigned char B[])170 PetscErrorCode  PetscDrawUtilitySetCmap(const char colormap[],int mapsize,unsigned char R[],unsigned char G[],unsigned char B[])
171 {
172   int             i,j;
173   const char      *cmap_name_list[sizeof(PetscDrawCmapTable)/sizeof(PetscDrawCmapTable[0])];
174   PetscInt        id = 0, count = (PetscInt)(sizeof(cmap_name_list)/sizeof(char*));
175   PetscBool       reverse = PETSC_FALSE, brighten = PETSC_FALSE;
176   PetscReal       beta = 0;
177   PetscErrorCode  ierr;
178 
179   PetscFunctionBegin;
180   for (i=0; i<count; i++) cmap_name_list[i] = PetscDrawCmapTable[i].name;
181   if (colormap && colormap[0]) {
182     PetscBool match = PETSC_FALSE;
183     for (id=0; !match && id<count; id++) {ierr = PetscStrcasecmp(colormap,cmap_name_list[id],&match);CHKERRQ(ierr);}
184     if (!match) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Colormap '%s' not found",colormap);
185   }
186   ierr = PetscOptionsGetEList(NULL,NULL,"-draw_cmap",cmap_name_list,count,&id,NULL);CHKERRQ(ierr);
187   ierr = PetscOptionsGetBool(NULL,NULL,"-draw_cmap_reverse",&reverse,NULL);CHKERRQ(ierr);
188   ierr = PetscOptionsGetReal(NULL,NULL,"-draw_cmap_brighten",&beta,&brighten);CHKERRQ(ierr);
189   if (brighten && (beta <= (PetscReal)-1 || beta >= (PetscReal)+1)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"brighten parameter %g must be in the range (-1,1)",(double)beta);
190 
191   if (PetscDrawCmapTable[id].cmap) {
192     ierr = PetscDrawCmapTable[id].cmap(mapsize,R,G,B);CHKERRQ(ierr);
193   } else {
194     const unsigned char (*rgb)[3] = PetscDrawCmapTable[id].data;
195     if (mapsize != 256-PETSC_DRAW_BASIC_COLORS) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Colormap '%s' with size %d not supported",cmap_name_list[id],mapsize);
196     for (i=0; i<mapsize; i++) {R[i] = rgb[i][0]; G[i] = rgb[i][1]; B[i] = rgb[i][2];}
197   }
198 
199   if (reverse) {
200     i = 0; j = mapsize-1;
201     while (i < j) {
202 #define SWAP(a,i,j) do { unsigned char t = a[i]; a[i] = a[j]; a[j] = t; } while (0)
203       SWAP(R,i,j);
204       SWAP(G,i,j);
205       SWAP(B,i,j);
206 #undef SWAP
207       i++; j--;
208     }
209   }
210 
211   if (brighten) {
212     PetscReal gamma = (beta > 0.0) ? (1 - beta) : (1 / (1 + beta));
213     for (i=0; i<mapsize; i++) {
214       PetscReal r = PetscPowReal((PetscReal)R[i]/255,gamma);
215       PetscReal g = PetscPowReal((PetscReal)G[i]/255,gamma);
216       PetscReal b = PetscPowReal((PetscReal)B[i]/255,gamma);
217       R[i] = (unsigned char)(255*PetscMin(r,(PetscReal)1.0));
218       G[i] = (unsigned char)(255*PetscMin(g,(PetscReal)1.0));
219       B[i] = (unsigned char)(255*PetscMin(b,(PetscReal)1.0));
220     }
221   }
222   PetscFunctionReturn(0);
223 }
224