1    | /***************************************
2    |   $Header$
3    | 
4    |   This file contains the rendering code for Illuminator, which renders 2-D or
5    |   3-D data into an RGB(A) unsigned char array (using perspective in 3-D).
6    | ***************************************/
7    | 
8    | #include "illuminator.h"
9    | #include "structures.h"
10   | 
11   | /* Build with -DDEBUG for debugging output */
12   | #undef DPRINTF
13   | #ifdef DEBUG
14   | #define DPRINTF(fmt, args...) PetscPrintf (PETSC_COMM_WORLD, "%s: " fmt, __FUNCT__, args)
15   | #else
16   | #define DPRINTF(fmt, args...)
17   | #endif
18   | 
19   | 
20   | #undef __FUNCT__
21   | #define __FUNCT__ "pseudocolor"
22   | 
23   | /*++++++++++++++++++++++++++++++++++++++
24   |   This little function converts a scalar value into an rgb color from red to
25   |   blue.
26   | 
27   |   PetscScalar val Value to convert.
28   | 
29   |   PetscScalar* scale Array with minimum and maximum values in which to scale
30   |   val.
31   | 
32   |   guchar *pixel Address in rgb buffer where this pixel should be painted.
33   |   ++++++++++++++++++++++++++++++++++++++*/
34   | 
35   | static inline void pseudocolor (PetscScalar val, PetscScalar* scale,
36   | 				guchar *pixel)
37   | {
38   |   PetscScalar shade = (val - scale[0]) / (scale[1] - scale[0]);
39   |   /* Old stuff *
40   |   if (shade < 0.2) {      /* red -> yellow *
41   |     *pixel++ = 255;
42   |     *pixel++ = 1275*shade;
43   |     *pixel++ = 0; }
44   |   else if (shade < 0.4) { /* yellow -> green *
45   |     *pixel++ = 510-1275*shade;
46   |     *pixel++ = 255;
47   |     *pixel++ = 0; }
48   |   else if (shade < 0.6) { /* green -> cyan *
49   |     *pixel++ = 0;
50   |     *pixel++ = 255;
51   |     *pixel++ = 1275*shade-510; }
52   |   else if (shade < 0.8) { /* cyan -> blue *
53   |     *pixel++ = 0;
54   |     *pixel++ = 1020-1275*shade;
55   |     *pixel++ = 255; }
56   |   else {                  /* blue -> magenta *
57   |     *pixel++ = 1275*shade-1020;
58   |     *pixel++ = 0;
59   |     *pixel++ = 255; }
60   |   /* New stuff */
61   |   if (shade < 0.25) {      /* red -> yellow */
62   |     *pixel++ = 255;
63   |     *pixel++ = 1020*shade;
64   |     *pixel++ = 0; }
65   |   else if (shade < 0.5) { /* yellow -> green */
66   |     *pixel++ = 510-1020*shade;
67   |     *pixel++ = 255;
68   |     *pixel++ = 0; }
69   |   else if (shade < 0.75) { /* green -> cyan */
70   |     *pixel++ = 0;
71   |     *pixel++ = 255;
72   |     *pixel++ = 1020*shade-510; }
73   |   else { /* cyan -> blue */
74   |     *pixel++ = 0;
75   |     *pixel++ = 1020-1020*shade;
76   |     *pixel++ = 255; }
77   | }
78   | 
79   | 
80   | #undef __FUNCT__
81   | #define __FUNCT__ "pseudoternarycolor"
82   | 
83   | /*++++++++++++++++++++++++++++++++++++++
84   |   This little function converts two ternary fractions into an rgb color, with
85   |   yellow, cyan and magenta indicating the corners.
86   | 
87   |   PetscScalar A First ternary fraction.
88   | 
89   |   PetscScalar B Second ternary fraction.
90   | 
91   |   PetscScalar *scale Array first and second ternary fractions of each of the
92   |   three corner values for scaling.
93   | 
94   |   guchar *pixel  Address in rgb buffer where this pixel should be painted.
95   |   ++++++++++++++++++++++++++++++++++++++*/
96   | 
97   | static inline void pseudoternarycolor
98   | (PetscScalar A, PetscScalar B, PetscScalar *scale, guchar *pixel)
99   | {
100  |   PetscScalar x1, x2, inverse_det;
101  | 
102  |   /* Transform A,B into x1,x2 based on corners */
103  |   inverse_det = 1./((scale[2]-scale[0])*(scale[5]-scale[1]) -
104  | 		    (scale[4]-scale[0])*(scale[3]-scale[1]));
105  |   x1 = ((A-scale[0])*(scale[5]-scale[1]) -
106  | 	(B-scale[1])*(scale[4]-scale[0])) * inverse_det;
107  |   x2 = ((B-scale[1])*(scale[2]-scale[0]) -
108  | 	(A-scale[0])*(scale[3]-scale[1])) * inverse_det;
109  | 
110  |   /* Now colorize */
111  |   /* Simple R-G-B */
112  |   /* *pixel++ = 255*x1;
113  |   *pixel++ = 255*(1.-x1-x2);
114  |   *pixel++ = 255*sqrt(x2); */
115  |   /* *pixel++ = 255*(1.-x1);
116  |   *pixel++ = 255*(x1+x2);
117  |   *pixel++ = 255*(1.-x2); */
118  |   /* Fancy three-square colorization */
119  |   if (2*x1+x2<=1 && 2*x2+x1<=1)
120  |     {
121  |       *pixel++ = 245;
122  |       *pixel++ = 245*(2*x1+3*x1*x2);
123  |       *pixel = 245*(2*x2+3*x1*x2);
124  |     }
125  |   else if (x2>=x1)
126  |     {
127  |       *pixel++ = 245*(2*(1-x1-x2)+3*(1-x1-x2)*x1);
128  |       *pixel++ = 245*(2*x1+3*(1-x1-x2)*x1);
129  |       *pixel = 245;
130  |     }
131  |   else
132  |     {
133  |       *pixel++ = 245*(2.-2*x1-2*x2+3*(1-x1-x2)*x2);
134  |       *pixel++ = 245;
135  |       *pixel = 245*(2*x2+3*(1-x1-x2)*x2);
136  |     }
137  |   /* Even fancier three-square colorization */
138  |   /*if (2*x1+x2<=1 && 4*x2+x1<=1) // 0,0 corner
139  |     {
140  |       *pixel++ = 254;
141  |       *pixel++ = 254*(2*x1+7./6.*x1*x2);
142  |       *pixel = 254*(4*x2+7./4.*x1*x2);
143  |     }
144  |   else if (3*x2>=x1) // 0,1 corner
145  |     {
146  |       *pixel++ = 254*(4./3.*(1-x1-x2)+7./3.*(1-x1-x2)*x1);
147  |       *pixel++ = 254*(4./3.*x1+7./3.*(1-x1-x2)*x1);
148  |       *pixel = 254;
149  |     }
150  |   else // 1,0 corner
151  |     {
152  |       *pixel++ = 254*(2*(1-x1-x2)+7./3.*(1-x1-x2)*x2);
153  |       *pixel++ = 254;
154  |       *pixel = 254*(4*x2+7./3.*(1-x1-x2)*x2);
155  |     }*/
156  | }
157  | 
158  | 
159  | #undef __FUNCT__
160  | #define __FUNCT__ "pseudoternarysquarecolor"
161  | 
162  | /*++++++++++++++++++++++++++++++++++++++
163  |   This little function converts two composition values into an rgb color, with
164  |   blue, red, green and yellow indicating the corners.
165  | 
166  |   PetscScalar A First ternary composition.
167  | 
168  |   PetscScalar B Second ternary composition.
169  | 
170  |   PetscScalar *scale Array: minimum and maximum first and second compositions
171  |   for scaling.
172  | 
173  |   guchar *pixel  Address in rgb buffer where this pixel should be painted.
174  |   ++++++++++++++++++++++++++++++++++++++*/
175  | 
176  | static inline void pseudoternarysquarecolor
177  | (PetscScalar A, PetscScalar B, PetscScalar *scale, guchar *pixel)
178  | {
179  |   PetscScalar x1, x2, inverse_det;
180  | 
181  |   /* Transform A,B into x1,x2 based on corners */
182  |   x1 = (A-scale[0])/(scale[1]-scale[0]);
183  |   x2 = (B-scale[2])/(scale[3]-scale[2]);
184  | 
185  |   /* Now colorize */
186  |   *pixel++ = 255*x2;
187  |   *pixel++ = 255*x1;
188  |   *pixel++ = 255*(1.-x1)*(1.-x2);
189  | }
190  | 
191  | 
192  | #undef __FUNCT__
193  | #define __FUNCT__ "pseudohueintcolor"
194  | 
195  | /*++++++++++++++++++++++++++++++++++++++
196  |   This little function converts a vector into an rgb color with hue indicating
197  |   direction (green, yellow, red, blue at 0, 90, 180, 270 degrees) and intensity
198  |   indicating magnitude relative to reference magnitude in scale[1].
199  | 
200  |   PetscScalar vx Vector's
201  |   +latex+$x$-component.
202  |   +html+ <i>x</i>-component.
203  | 
204  |   PetscScalar vy Vector's
205  |   +latex+$y$-component.
206  |   +html+ <i>y</i>-component.
207  | 
208  |   PetscScalar *scale Array whose second entry has the reference magnitude.
209  | 
210  |   guchar *pixel Address in rgb buffer where this pixel should be painted.
211  |   ++++++++++++++++++++++++++++++++++++++*/
212  | 
213  | static inline void pseudohueintcolor
214  | (PetscScalar vx, PetscScalar vy, PetscScalar *scale, guchar *pixel)
215  | {
216  |   PetscScalar mag=sqrt(vx*vx+vy*vy), theta=atan2(vy,vx), red, green, blue;
217  |   if (scale[1] <= 0.)
218  |     {
219  |       *pixel = *(pixel+1) = *(pixel+2) = 0;
220  |       return;
221  |     }
222  |   mag = (mag > scale[1]) ? 1.0 : mag/scale[1];
223  | 
224  |   red = 2./M_PI * ((theta<-M_PI/2.) ? -M_PI/2.-theta :
225  | 		   ((theta<0.) ? 0. : ((theta<M_PI/2.) ? theta : M_PI/2.)));
226  |   green = 2./M_PI * ((theta<-M_PI/2.) ? 0. :
227  | 		     ((theta<0.) ? theta+M_PI/2. :
228  | 		      ((theta<M_PI/2.) ? M_PI/2. : M_PI-theta)));
229  |   blue = 2./M_PI * ((theta<-M_PI/2.) ? theta+M_PI :
230  | 		    ((theta<0.) ? -theta : 0.));
231  | 
232  |   *pixel++ = 255*mag*red;
233  |   *pixel++ = 255*mag*green;
234  |   *pixel++ = 255*mag*blue;
235  | }
236  | 
237  | 
238  | #undef __FUNCT__
239  | #define __FUNCT__ "pseudoshearcolor"
240  | 
241  | /*++++++++++++++++++++++++++++++++++++++
242  |   This little function converts a shear tensor (symmetric, diagonals sum to
243  |   zero) into an rgb color with hue indicating direction of the tensile stress
244  |   (red, yellow, green, cyan, blue, magenta at 0, 30, 60, 90, 120, 150 degrees
245  |   respectively; 180 is equivalent to zero for stress) and intensity
246  |   indicating magnitude relative to reference magnitude in scale[0].
247  | 
248  |   PetscScalar gxx Tensor's
249  |   +latex+$xx$-component.
250  |   +html+ <i>xx</i>-component.
251  | 
252  |   PetscScalar gxy Tensor's
253  |   +latex+$xy$-component.
254  |   +html+ <i>xy</i>-component.
255  | 
256  |   PetscScalar *scale Array whose first entry has the reference magnitude.
257  | 
258  |   guchar *pixel Address in rgb buffer where this pixel should be painted.
259  |   ++++++++++++++++++++++++++++++++++++++*/
260  | 
261  | static inline void pseudoshearcolor
262  | (PetscScalar gxx, PetscScalar gxy, PetscScalar *scale, guchar *pixel)
263  | {
264  |   PetscScalar mag=sqrt(gxx*gxx+gxy*gxy), theta=atan2(gxy,gxx),
265  |     red,green,blue;
266  |   if (scale[0] <= 0.)
267  |     {
268  |       *pixel = *(pixel+1) = *(pixel+2) = 0;
269  |       return;
270  |     }
271  |   mag = (mag > scale[0]) ? 1.0 : mag/scale[0];
272  | 
273  |   red = 3./M_PI * ((theta < -2.*M_PI/3.) ? 0. :
274  | 		   ((theta < -M_PI/3.) ? theta + 2.*M_PI/3. :
275  | 		    ((theta < M_PI/3.) ? M_PI/3. :
276  | 		     ((theta < 2.*M_PI/3.) ? 2.*M_PI/3. - theta: 0.))));
277  |   green = 3./M_PI * ((theta < -M_PI/3.) ? M_PI/3. :
278  | 		     ((theta < 0.) ? -theta :
279  | 		      ((theta < 2.*M_PI/3.) ? 0. : theta - 2.*M_PI/3.)));
280  |   blue = 3./M_PI * ((theta < -2.*M_PI/3) ? -2.*M_PI/3. - theta :
281  | 		    ((theta < 0.) ? 0. :
282  | 		     ((theta < M_PI/3.) ? theta : M_PI/3.)));
283  | 
284  |   *pixel++ = 255*mag*red;
285  |   *pixel++ = 255*mag*green;
286  |   *pixel++ = 255*mag*blue;
287  | }
288  | 
289  | 
290  | #undef __FUNCT__
291  | #define __FUNCT__ "render_scale_2d"
292  | 
293  | /*++++++++++++++++++++++++++++++++++++++
294  |   This draws a little rwidth x rheight image depicting the scale of a field
295  |   variable.
296  | 
297  |   int render_scale_2d It returns zero or an error code.
298  | 
299  |   IDisplay Disp Display object into which to draw the scale.
300  | 
301  |   field_plot_type fieldtype Type of plot.
302  | 
303  |   int symmetry Symmetry order for vector scale image.
304  |   ++++++++++++++++++++++++++++++++++++++*/
305  | 
306  | int render_scale_2d (IDisplay Disp, field_plot_type fieldtype, int symmetry)
307  | {
308  |   PetscScalar scale[6];
309  |   int i, j, myheight;
310  |   struct idisplay *thedisp = (struct idisplay *) Disp;
311  | 
312  |   switch (fieldtype)
313  |     {
314  |     case FIELD_SCALAR:
315  |       scale[0]=0.;
316  |       scale[1]=1.;
317  |       for (j=0; j<thedisp->rgb_height; j++)
318  | 	for (i=0; i<thedisp->rgb_width; i++)
319  | 	  pseudocolor ((PetscScalar) i/(thedisp->rgb_width-1), scale,
320  | 		       thedisp->rgb+(j*thedisp->rgb_rowskip + i) *
321  | 		       thedisp->rgb_bpp);
322  |       return 0;
323  | 
324  |     case FIELD_VECTOR:
325  |       scale[0] = 0.;
326  |       scale[1] = 1.;
327  |       myheight = ((thedisp->rgb_height > thedisp->rgb_width) ?
328  | 		  thedisp->rgb_width : thedisp->rgb_height)/2;
329  |       for (j=-myheight; j<myheight; j++)
330  | 	for (i=-(int)(myheight*sqrt(1.-((PetscScalar)j*j/myheight/myheight)));
331  | 	     i<(int)(myheight*sqrt(1.-((PetscScalar)j*j/myheight/myheight)));
332  | 	     i++)
333  | 	  {
334  | 	    PetscScalar mag = (PetscScalar) sqrt (i*i+j*j)/myheight,
335  | 	      theta = atan2 (j, i);
336  | 	    pseudohueintcolor
337  | 	      (mag * cos (symmetry*theta), mag * sin (symmetry*theta), scale,
338  | 	       thedisp->rgb + ((myheight-j)*thedisp->rgb_rowskip +
339  | 			       thedisp->rgb_width/2+i)*thedisp->rgb_bpp);
340  | 	  }
341  |       return 0;
342  | 
343  |     case FIELD_TERNARY:
344  |       scale[0] = scale[1] = scale[3] = scale[4] = 0.;
345  |       scale[2] = scale[5] = 1.;
346  |       myheight = ((PetscScalar) thedisp->rgb_height >
347  | 		  0.5*sqrt(3.) * thedisp->rgb_width) ?
348  | 	(int) (0.5*sqrt(3.)* thedisp->rgb_width) : thedisp->rgb_height;
349  |       for (j=0; j<myheight; j++)
350  | 	for (i=0; i<(int) (2./sqrt(3.)*(myheight-j)); i++)
351  | 	  pseudoternarycolor
352  | 	    ((PetscScalar) i/(myheight*2./sqrt(3.)),
353  | 	     (PetscScalar) j/myheight, scale, thedisp->rgb +
354  | 	     (((thedisp->rgb_height+myheight)/2 -j) * thedisp->rgb_rowskip +
355  | 	      i + (int)(sqrt(3.)/3.*j)) * thedisp->rgb_bpp);
356  |       return 0;
357  | 
358  |     case FIELD_TERNARY_SQUARE:
359  |       scale[0] = scale[2] = 0.;
360  |       scale[1] = scale[3] = 1.;
361  |       for (j=0; j<thedisp->rgb_height; j++)
362  | 	for (i=0; i<thedisp->rgb_width; i++)
363  | 	  pseudoternarysquarecolor
364  | 	    ((PetscScalar) i/(thedisp->rgb_width-1),
365  | 	     (PetscScalar) j/(thedisp->rgb_height-1), scale,
366  | 	     thedisp->rgb+(j*thedisp->rgb_rowskip + i) * thedisp->rgb_bpp);
367  |       return 0;
368  | 
369  |     case FIELD_TENSOR_SHEAR:
370  |       scale[0] = 1.;
371  |       myheight = ((thedisp->rgb_height > thedisp->rgb_width) ?
372  | 		  thedisp->rgb_width : thedisp->rgb_height)/2;
373  |       for (j=-myheight; j<myheight; j++)
374  | 	for (i=-(int)(myheight*sqrt(1.-((PetscScalar)j*j/myheight/myheight)));
375  | 	     i<(int)(myheight*sqrt(1.-((PetscScalar)j*j/myheight/myheight)));
376  | 	     i++)
377  | 	  {
378  | 	    PetscScalar mag = (PetscScalar) sqrt (i*i+j*j)/myheight,
379  | 	      theta = atan2 (j, i);
380  | 	    pseudoshearcolor
381  | 	      (mag * cos (2*theta), mag * sin (2*theta), scale,
382  | 	       thedisp->rgb + ((myheight-j)*thedisp->rgb_rowskip +
383  | 			       thedisp->rgb_width/2+i)*thedisp->rgb_bpp);
384  | 	  }
385  |       return 0;
386  |     }
387  |   SETERRQ (PETSC_ERR_ARG_OUTOFRANGE, "Field type not yet supported");
388  | }
389  | 
390  | 
391  | #undef __FUNCT__
392  | #define __FUNCT__ "render_composition_path"
393  | 
394  | 
395  | /*++++++++++++++++++++++++++++++++++++++
396  |   Render a composition path, such as a diffusion path or phase diagram, onto an
397  |   IDisplay image.
398  | 
399  |   int render_composition_path Returns zero or an error code.
400  | 
401  |   IDisplay Disp IDisplay object into which to draw the path.
402  | 
403  |   PetscScalar comp_array Array of compositions for drawing.
404  | 
405  |   int gridpoints Number of gridpoints in the composition array.
406  | 
407  |   int num_fields Total number of fields in the composition array.
408  | 
409  |   field_plot_type fieldtype The type of this field.
410  | 
411  |   PetscScalar *scale Array of minimum and maximum values to pass to the
412  |   various pseudocolor functions; if NULL, call auto_scale to determine those
413  |   values.
414  | 
415  |   PetscScalar red Red color for composition path.
416  | 
417  |   PetscScalar green Green color for composition path.
418  | 
419  |   PetscScalar blue Blue color for composition path.
420  | 
421  |   PetscScalar alpha Alpha channel for composition path.
422  |   ++++++++++++++++++++++++++++++++++++++*/
423  | 
424  | int render_composition_path
425  | (IDisplay Disp, PetscScalar *comp_array, int gridpoints, int num_fields,
426  |  field_plot_type fieldtype, PetscScalar *scale,
427  |  PetscScalar red,PetscScalar green,PetscScalar blue,PetscScalar alpha)
428  | {
429  |   int i, ierr;
430  |   PetscScalar local_scale[6];
431  |   struct idisplay *thedisp = (struct idisplay *) Disp;
432  | 
433  |   /* Sanity checks */
434  |   if (!thedisp)
435  |     SETERRQ (PETSC_ERR_ARG_CORRUPT, "Null display object");
436  |   if (!(thedisp->rgb))
437  |     SETERRQ (PETSC_ERR_ARG_WRONGSTATE,
438  | 	     "Display object has no RGB buffer to draw in");
439  | 
440  |   /* Determine default min and max if none are provided */
441  |   if (scale == NULL)
442  |     {
443  |       scale = local_scale;
444  |       ierr = auto_scale (comp_array, gridpoints, num_fields, 0, fieldtype, 2,
445  | 			 scale); CHKERRQ (ierr);
446  |     }
447  | 
448  |   for (i=0; i<gridpoints; i++)
449  |     {
450  |       PetscScalar A=comp_array [i*num_fields], B=comp_array [i*num_fields+1];
451  |       int sx, sy;
452  |       guchar *spixel;
453  | 
454  |       if (fieldtype == FIELD_TERNARY)
455  | 	{
456  | 	  int myheight = ((PetscScalar) thedisp->rgb_height >
457  | 			  0.5*sqrt(3.)*thedisp->rgb_width) ?
458  | 	    (int) (0.5*sqrt(3.)* thedisp->rgb_width) : thedisp->rgb_height;
459  | 	  PetscScalar inverse_det =
460  | 	    1./((scale[2]-scale[0])*(scale[5]-scale[1]) -
461  | 		(scale[4]-scale[0])*(scale[3]-scale[1]));
462  | 	  sy = myheight*(((B-scale[1])*(scale[2]-scale[0]) -
463  | 			  (A-scale[0])*(scale[3]-scale[1])) *inverse_det);
464  | 	  sx = thedisp->rgb_width* (((A-scale[0])*(scale[5]-scale[1]) -
465  | 				     (B-scale[1])*(scale[4]-scale[0])) *inverse_det)+
466  | 	    thedisp->rgb_width*sy/myheight/2;
467  | 	  sy = (thedisp->rgb_height+myheight)/2 - sy;
468  | 	}
469  |       else if (fieldtype == FIELD_TERNARY_SQUARE)
470  | 	{
471  | 	  sx = thedisp->rgb_width * (A-scale[0])/(scale[1]-scale[0]);
472  | 	  sy = thedisp->rgb_height*(B-scale[2])/(scale[3]-scale[2]);
473  | 	}
474  |       else
475  | 	sx=sy=-1; /* Don't draw diffusion path */
476  | 
477  |       if (sx>=0 && sy>=0 && sx<thedisp->rgb_width &&
478  | 	  sy<thedisp->rgb_height)
479  | 	{
480  | 	  spixel = thedisp->rgb +
481  | 	    thedisp->rgb_bpp*(sy*thedisp->rgb_rowskip + sx);
482  | 	  *spixel = (1.-alpha)* *spixel + red*alpha*255;
483  | 	  *(spixel+1) = (1.-alpha)* *(spixel+1) + green*alpha*255;
484  | 	  *(spixel+2) = (1.-alpha)* *(spixel+2) + blue*alpha*255;
485  | 	}
486  |     }
487  |   return 0;
488  | }
489  | 
490  | 
491  | #undef __FUNCT__
492  | #define __FUNCT__ "render_rgb_local_2d"
493  | 
494  | /*++++++++++++++++++++++++++++++++++++++
495  |   Render data from global_array into local part of an RGB buffer.  When running
496  |   in parallel, these local buffers should be collected and layered to produce
497  |   the full image.
498  | 
499  |   int render_rgb_local_2d Returns zero or an error code.
500  | 
501  |   IDisplay Disp Display object holding the RGB buffer and its characteristics.
502  | 
503  |   PetscScalar *global_array Local array of global vector data to render.
504  | 
505  |   int num_fields Number of field variables in the array.
506  | 
507  |   int display_field The (first) field we are rendering now.
508  | 
509  |   field_plot_type fieldtype The type of this field.
510  | 
511  |   PetscScalar *scale Array of minimum and maximum values to pass to the
512  |   various pseudocolor functions; if NULL, call auto_scale to determine those
513  |   values.
514  | 
515  |   int nx Width of the array.
516  | 
517  |   int ny Height of the array.
518  | 
519  |   int xs Starting
520  |   +latex+$x$-coordinate
521  |   +html+ <i>x</i>-coordinate
522  |   of the local part of the global vector.
523  | 
524  |   int ys Starting
525  |   +latex+$y$-coordinate
526  |   +html+ <i>y</i>-coordinate
527  |   of the local part of the global vector.
528  | 
529  |   int xm Width of the local part of the global vector.
530  | 
531  |   int ym Height of the local part of the global vector.
532  | 
533  |   int transform Transformation flags.
534  | 
535  |   IDisplay SDisp Display object in which to draw the diffusion path (optional,
536  |   ignored if NULL).
537  | 
538  |   PetscScalar dpred Red color for diffusion path points (0-1).
539  | 
540  |   PetscScalar dpgreen Green color for diffusion path points (0-1).
541  | 
542  |   PetscScalar dpblue Blue color for diffusion path points (0-1).
543  | 
544  |   PetscScalar dpalpha Alpha channel for diffusion path points (0-1).
545  |   ++++++++++++++++++++++++++++++++++++++*/
546  | 
547  | int render_rgb_local_2d
548  | (IDisplay Disp, PetscScalar *global_array, int num_fields, int display_field,
549  |  field_plot_type fieldtype, PetscScalar *scale, int nx,int ny, int xs,int ys,
550  |  int xm,int ym, int transform, IDisplay SDisp, PetscScalar dpred,
551  |  PetscScalar dpgreen, PetscScalar dpblue, PetscScalar dpalpha)
552  | {
553  |   int ix,iy, ierr=0;
554  |   PetscScalar local_scale[6];
555  |   struct idisplay *thedisp = (struct idisplay *) Disp;
556  | 
557  |   /* Sanity checks */
558  |   if (!thedisp)
559  |     SETERRQ (PETSC_ERR_ARG_CORRUPT, "Null display object");
560  |   if (!(thedisp->rgb))
561  |     SETERRQ (PETSC_ERR_ARG_WRONGSTATE,
562  | 	     "Display object has no RGB buffer to draw in");
563  | 
564  |   /* Determine default min and max if none are provided */
565  |   if (scale == NULL)
566  |     {
567  |       scale = local_scale;
568  |       ierr = auto_scale (global_array, xm*ym, num_fields, display_field,
569  | 			   fieldtype, 2, scale); CHKERRQ (ierr);
570  |     }
571  |   DPRINTF ("Final scale: %g %g %g %g %g %g\n", scale[0],
572  | 	   scale[1], scale[2], scale[3], scale[4], scale[5]);
573  | 
574  |   /* Do the rendering (note switch in inner loop, gotta get rid of that) */
575  |   for (iy=thedisp->rgb_height*ys/ny; iy<thedisp->rgb_height*(ys+ym)/ny; iy++)
576  |     for (ix=thedisp->rgb_width*xs/nx; ix<thedisp->rgb_width*(xs+xm)/nx; ix++)
577  |       {
578  | 	int localx, localy;
579  | 	if (transform & ROTATE_LEFT)
580  | 	  {
581  | 	    localx = ((transform & FLIP_HORIZONTAL) ?
582  | 		      iy : thedisp->rgb_height-iy-1) * nx/thedisp->rgb_height;
583  | 	    localy = ((transform & FLIP_VERTICAL) ?
584  | 		      ix : thedisp->rgb_width-ix-1) * ny/thedisp->rgb_width;
585  | 	  }
586  | 	else
587  | 	  {
588  | 	    localx = ((transform & FLIP_HORIZONTAL) ?
589  | 		      thedisp->rgb_width-ix-1 : ix) * nx/thedisp->rgb_width;
590  | 	    localy = ((transform & FLIP_VERTICAL) ?
591  | 		      iy : thedisp->rgb_height-iy-1) * ny/thedisp->rgb_height;
592  | 	  }
593  | 
594  | 	int vecindex = (localy*xm + localx)*num_fields + display_field;
595  | 	guchar *pixel = thedisp->rgb +
596  | 	  thedisp->rgb_bpp*(iy*thedisp->rgb_rowskip + ix);
597  | 
598  | 	switch (fieldtype)
599  | 	  {
600  | 	  case FIELD_SCALAR:
601  | 	  case FIELD_SCALAR+1:
602  | 	    {
603  | 	      pseudocolor (global_array [vecindex], scale, pixel);
604  | 	      break;
605  | 	    }
606  | 	  case FIELD_TERNARY:
607  | 	    {
608  | 	      pseudoternarycolor
609  | 		(global_array [vecindex], global_array [vecindex+1], scale,
610  | 		 pixel);
611  | 	      break;
612  | 	    }
613  | 	  case FIELD_TERNARY_SQUARE:
614  | 	    {
615  | 	      pseudoternarysquarecolor
616  | 		(global_array [vecindex], global_array [vecindex+1], scale,
617  | 		 pixel);
618  | 	      break;
619  | 	    }
620  | 	  case FIELD_VECTOR:
621  | 	  case FIELD_VECTOR+1:
622  | 	    {
623  | 	      pseudohueintcolor
624  | 		(global_array [vecindex], global_array [vecindex+1], scale,
625  | 		 pixel);
626  | 	      break;
627  | 	    }
628  | 	  case FIELD_TENSOR_SHEAR:
629  | 	    {
630  | 	      pseudoshearcolor
631  | 		(global_array [vecindex], global_array [vecindex+1], scale,
632  | 		 pixel);
633  | 	      break;
634  | 	    }
635  | 	  default:
636  | 	    SETERRQ (PETSC_ERR_ARG_OUTOFRANGE, "Field type not yet supported");
637  | 	  }
638  |       }
639  | 
640  |   /* (Optional) diffusion path */
641  |   DPRINTF ("Done main rendering, doing scale.\n",0);
642  |   if (SDisp)
643  |     ierr=render_composition_path (SDisp, global_array+display_field, xm*ym,
644  | 				  num_fields, fieldtype, scale,
645  | 				  dpred,dpgreen,dpblue,dpalpha);
646  | 
647  |   return ierr;
648  | }
649  | 
650  | 
651  | #undef __FUNCT__
652  | #define __FUNCT__ "render_rgb_local_3d"
653  | 
654  | /*++++++++++++++++++++++++++++++++++++++
655  |   Render triangle data into an RGB buffer.  When called in parallel, the
656  |   resulting images should be layered to give the complete picture.  Zooming is
657  |   done by adjusting the ratio of the dir vector to the right vector.
658  |   +latex+\par
659  |   +html+ <p>
660  |   The coordinate transformation is pretty simple.
661  |   +latex+A point $\vec{p}$ in space is transformed to $x,y$ on the screen by
662  |   +latex+representing it as:
663  |   +latex+\begin{equation}
664  |   +latex+  \label{eq:transform}
665  |   +latex+  \vec{p} = \vec{\imath} + a \vec{d} + ax \vec{r} + ay \vec{u},
666  |   +latex+\end{equation}
667  |   +latex+where $\vec{\imath}$ is the observer's point (passed in as {\tt eye}),
668  |   +latex+$\vec{d}$ is the direction of observation (passed in as {\tt dir}),
669  |   +latex+$\vec{r}$ is the rightward direction to the observer (passed in as
670  |   +latex+{\tt right}), and $\vec{u}$ is the upward direction to the observer
671  |   +latex+which is given the direction of $\vec{r}\times\vec{d}$ and the
672  |   +latex+magnitude of $\vec{r}$.
673  |   +html+ A point <u><i>p</i></u> in space is transformed to <i>x, y</i> on the
674  |   +html+ screen by representing it as:
675  |   +html+ <center><u><i>p</i></u> = <u><i>i</i></u> + <i>a <u>d</u></i> +
676  |   +html+ <i>ax <u>r</u></i> + <i>ay <u>u</u></i>,</center>
677  |   +html+ where <u><i>i</i></u> is the observer's point (passed in as
678  |   +html+ <tt>eye</tt>), <u><i>d</i></u> is the direction of observation (passed
679  |   +html+ in as <tt>dir</tt>), <u><i>r</i></u> is the rightward direction to the
680  |   +html+ observer (passed in as <tt>right</tt>), and <u><i>u</i></u> is the
681  |   +html+ upward direction to the observer which is given the direction of
682  |   +html+ the cross product of <u><i>d</i></u> and <u><i>r</i></u> and the
683  |   +html+ magnitude of <u><i>r</i></u>.
684  |   +latex+\par
685  |   +html+ <p>
686  |   +latex+This system in equation \ref{eq:transform} can easily be solved for
687  |   +latex+$x$ and $y$ by first making it into a matrix equation:
688  |   +latex+\begin{equation}
689  |   +latex+  \label{eq:matrix}
690  |   +latex+  \left(\begin{array}{ccc} d_x&r_x&u_x \\ d_y&r_y&u_y \\ d_z&r_z&u_z
691  |   +latex+    \end{array}\right)
692  |   +latex+  \left(\begin{array}{c} a \\ ax \\ ay \end{array}\right) =
693  |   +latex+  \left(\begin{array}{c} p_x-i_x \\ p_y-i_y \\ p_z-i_z
694  |   +latex+    \end{array}\right).
695  |   +latex+\end{equation}
696  |   +latex+Calling the matrix $M$ the inverse of the matrix on the left side of
697  |   +latex+equation \ref{eq:matrix} gives the result:
698  |   +latex+\begin{eqnarray}
699  |   +latex+  a&=& M_{00} (p_x-i_x) + M_{01} (p_y-i_y) + M_{02} (p_z-i_z),
700  |   +latex+  \\ x&=& \frac{1}{a}[M_{10}(p_x-i_x)+M_{11}(p_y-i_y)+M_{12}(p_z-i_z)],
701  |   +latex+  \\ y&=& \frac{1}{a}[M_{20}(p_x-i_x)+M_{21}(p_y-i_y)+M_{22}(p_z-i_z)].
702  |   +latex+\end{eqnarray}
703  |   +html+ This system can easily be solved for <i>x</i> and <i>y</i> by first
704  |   +html+ making it into a matrix equation:
705  |   +html+ <center>(<i><u>d</u> <u>r</u> <u>u</u></i>) (<i>a, ax, ay</i>) =
706  |   +html+ <u><i>p</i></u> - <u><i>i</i></u>.</center>
707  |   +html+ Calling the matrix <u><i>M</i></u> the inverse of the matrix on the
708  |   +html+ left side gives the result:
709  |   +html+ <center><table>
710  |   +html+ <tr><td><i>a</i></td><td>=</td><td>
711  |   +html+ <i>M</i><sub>00</sub> (<i>p<sub>x</sub></i>-<i>i<sub>x</sub></i>) +
712  |   +html+ <i>M</i><sub>01</sub> (<i>p<sub>y</sub></i>-<i>i<sub>y</sub></i>) +
713  |   +html+ <i>M</i><sub>02</sub> (<i>p<sub>z</sub></i>-<i>i<sub>z</sub></i>),
714  |   +html+ </td></tr><tr><td><i>x</i></td><td>=</td><td>
715  |   +html+ [<i>M</i><sub>10</sub> (<i>p<sub>x</sub></i>-<i>i<sub>x</sub></i>) +
716  |   +html+ <i>M</i><sub>11</sub> (<i>p<sub>y</sub></i>-<i>i<sub>y</sub></i>) +
717  |   +html+ <i>M</i><sub>12</sub> (<i>p<sub>z</sub></i>-<i>i<sub>z</sub></i>)]
718  |   +html+ / <i>a</i>,</td></tr><tr><td><i>y</i></td><td>=</td><td>
719  |   +html+ [<i>M</i><sub>00</sub> (<i>p<sub>x</sub></i>-<i>i<sub>x</sub></i>) +
720  |   +html+ <i>M</i><sub>01</sub> (<i>p<sub>y</sub></i>-<i>i<sub>y</sub></i>) +
721  |   +html+ <i>M</i><sub>02</sub> (<i>p<sub>z</sub></i>-<i>i<sub>z</sub></i>)]
722  |   +html+ / <i>a</i>.</td></tr></table></center>
723  |   +latex+\par
724  |   +html+ <p>
725  |   The triple product of the vector from the light source to the triangle
726  |   centroid and the two vectors making up the triangle edges determines the
727  |   cosine of the angle made by the triangle normal and the incident light, whose
728  |   absolute value is used here to shade the triangle.  At this point, the light
729  |   source is taken as the observer location at
730  |   +latex+``{\tt eye}'',
731  |   +html+ "<tt>eye</tt>",
732  |   but that can easily be modified to use one or more independent light sources.
733  | 
734  |   int render_rgb_local_3d Returns zero or an error code.
735  | 
736  |   IDisplay Disp Display object holding the RGB buffer and its characteristics.
737  | 
738  |   ISurface Surf ISurface object containing triangles to render using imlib2
739  |   backend.
740  | 
741  |   PetscScalar *eye Point from where we're looking (x,y,z).
742  | 
743  |   PetscScalar *dir Direction we're looking (x,y,z).
744  | 
745  |   PetscScalar *right Rightward direction in physical space (x,y,z).
746  |   ++++++++++++++++++++++++++++++++++++++*/
747  | 
748  | int render_rgb_local_3d
749  | (IDisplay Disp, ISurface Surf, PetscScalar *eye, PetscScalar *dir,
750  |  PetscScalar *right)
751  | {
752  |   int *screen_coords, i;
753  |   PetscScalar *shades, up[3], m00,m01,m02,m10,m11,m12,m20,m21,m22, a;
754  |   struct isurface *thesurf = (struct isurface *) Surf;
755  |   struct idisplay *thedisp = (struct idisplay *) Disp;
756  | 
757  |   /* Sanity checks */
758  |   if (!thedisp)
759  |     SETERRQ (PETSC_ERR_ARG_CORRUPT, "Null display object");
760  |   if (!(thedisp->rgb))
761  |     SETERRQ (PETSC_ERR_ARG_WRONGSTATE,
762  | 	     "Display object has no RGB buffer to draw in");
763  | 
764  |   if (thedisp->rgb_bpp != 4)
765  |     SETERRQ (PETSC_ERR_ARG_SIZ, "Imlib2 rendering requires 4 bytes per pixel");
766  | 
767  |   /* Allocate memory for the shades and integer coordinates */
768  |   i = PetscMalloc (6*thesurf->num_triangles * sizeof(int), &screen_coords);
769  |   CHKERRQ (i);
770  |   i = PetscMalloc (thesurf->num_triangles * sizeof(PetscScalar), &shades);
771  |   CHKERRQ (i);
772  | 
773  |   /* Calculate the "up" vector as described above, storing the magnitude of dir
774  |      in a */
775  |   a = sqrt (dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]);
776  |   up[0] = (right[1]*dir[2] - right[2]*dir[1]) / a;
777  |   up[1] = (right[2]*dir[0] - right[0]*dir[2]) / a;
778  |   up[2] = (right[0]*dir[1] - right[1]*dir[0]) / a;
779  | 
780  |   /* Calculate the inverse of the dir, right, up columnar matrix, storing the
781  |      determinant of the matrix in a */
782  |   a = (dir[0] * (right[1]*up[2] - right[2]*up[1]) +
783  |        dir[1] * (right[2]*up[0] - right[0]*up[2]) +
784  |        dir[2] * (right[0]*up[1] - right[1]*up[0]));
785  |   m00 = (right[1]*up[2]   - right[2]*up[1])    / a;
786  |   m01 = (right[2]*up[0]   - right[0]*up[2])    / a;
787  |   m02 = (right[0]*up[1]   - right[1]*up[0])    / a;
788  |   m10 = (up[1]   *dir[2]  - up[2]   *dir[1])   / a;
789  |   m11 = (up[2]   *dir[0]  - up[0]   *dir[2])   / a;
790  |   m12 = (up[0]   *dir[1]  - up[1]   *dir[0])   / a;
791  |   m20 = (dir[1] *right[2] - dir[2]  *right[1]) / a;
792  |   m21 = (dir[2] *right[0] - dir[0]  *right[2]) / a;
793  |   m22 = (dir[0] *right[1] - dir[1]  *right[0]) / a;
794  | 
795  |   /* Loop over triangles, transforming their coordinates and setting shading */
796  |   for (i=0; i<thesurf->num_triangles; i++)
797  |     {
798  |       PetscScalar pmi00,pmi01,pmi02,pmi10,pmi11,pmi12,pmi20,pmi21,pmi22, c[3];
799  | 
800  |       /* Calculate coordinates relative to observer p-i (hence pmi) */
801  |       pmi00 = thesurf->vertices[13*i]   - eye[0];
802  |       pmi01 = thesurf->vertices[13*i+1] - eye[1];
803  |       pmi02 = thesurf->vertices[13*i+2] - eye[2];
804  |       pmi10 = thesurf->vertices[13*i+3] - eye[0];
805  |       pmi11 = thesurf->vertices[13*i+4] - eye[1];
806  |       pmi12 = thesurf->vertices[13*i+5] - eye[2];
807  |       pmi20 = thesurf->vertices[13*i+6] - eye[0];
808  |       pmi21 = thesurf->vertices[13*i+7] - eye[1];
809  |       pmi22 = thesurf->vertices[13*i+8] - eye[2];
810  | 
811  |       /* Calculate x,y for point 0 */
812  |       a = m00*pmi00 + m01*pmi01 + m02*pmi02;
813  |       if (a < 0)
814  | 	screen_coords [6*i+0] = screen_coords [6*i+1] = -1;
815  |       else
816  | 	{
817  | 	  screen_coords [6*i+0] = thedisp->rgb_width/2 *
818  | 	    (1. + (m10*pmi00 + m11*pmi01 + m12*pmi02) / a);
819  | 	  screen_coords [6*i+1] = thedisp->rgb_height/2 - thedisp->rgb_width/2 *
820  | 	    (m20*pmi00 + m21*pmi01 + m22*pmi02) / a;
821  | 	}
822  | 
823  |       /* Calculate x,y for point 1 */
824  |       a = m00*pmi10 + m01*pmi11 + m02*pmi12;
825  |       if (a <= 0)
826  | 	screen_coords [6*i+2] = screen_coords [6*i+3] = -1;
827  |       else
828  | 	{
829  | 	  screen_coords [6*i+2] = thedisp->rgb_width/2 *
830  | 	    (1. + (m10*pmi10 + m11*pmi11 + m12*pmi12) / a);
831  | 	  screen_coords [6*i+3] = thedisp->rgb_height/2 - thedisp->rgb_width/2 *
832  | 	    (m20*pmi10 + m21*pmi11 + m22*pmi12) / a;
833  | 	}
834  | 
835  |       /* Calculate x,y for point 2 */
836  |       a = m00*pmi20 + m01*pmi21 + m02*pmi22;
837  |       if (a <= 0)
838  | 	screen_coords [6*i+4] = screen_coords [6*i+5] = -1;
839  |       else
840  | 	{
841  | 	  screen_coords [6*i+4] = thedisp->rgb_width/2 *
842  | 	    (1. + (m10*pmi20 + m11*pmi21 + m12*pmi22) / a);
843  | 	  screen_coords [6*i+5] = thedisp->rgb_height/2 - thedisp->rgb_width/2 *
844  | 	    (m20*pmi20 + m21*pmi21 + m22*pmi22) / a;
845  | 	}
846  | 
847  | #ifdef DEBUG
848  |       if (i==0)
849  | 	DPRINTF("First triangle: (%g,%g,%g), (%g,%g,%g), (%g,%g,%g);\n"
850  | 		"2-D coordinates: (%d,%d), (%d,%d), (%d,%d)\n",
851  | 		thesurf->vertices[0],thesurf->vertices[1],thesurf->vertices[2],
852  | 		thesurf->vertices[3],thesurf->vertices[4],thesurf->vertices[5],
853  | 		thesurf->vertices[6],thesurf->vertices[7],thesurf->vertices[8],
854  | 		screen_coords[0],screen_coords[1],screen_coords[3],
855  | 		screen_coords[3],screen_coords[4],screen_coords[5]);
856  | #endif
857  | 
858  |       /* Calculate relative triangle centroid, triangle arms, and shade */
859  |       c[0] = (pmi00 + pmi10 + pmi20) / 3.;
860  |       c[1] = (pmi01 + pmi11 + pmi21) / 3.;
861  |       c[2] = (pmi02 + pmi12 + pmi22) / 3.;
862  |       pmi10 -= pmi00;   pmi11 -= pmi01;   pmi12 -= pmi02;
863  |       pmi20 -= pmi00;   pmi21 -= pmi01;   pmi22 -= pmi02;
864  |       shades [i] = PetscAbsScalar (c[0] * (pmi11*pmi22 - pmi12*pmi21) +
865  | 				   c[1] * (pmi12*pmi20 - pmi10*pmi22) +
866  | 				   c[2] * (pmi10*pmi21 - pmi11*pmi20)) /
867  | 	sqrt (c[0]*c[0] + c[1]*c[1] + c[2]*c[2]) /
868  | 	sqrt (pmi10*pmi10 + pmi11*pmi11 + pmi12*pmi12) /
869  | 	sqrt (pmi20*pmi20 + pmi21*pmi21 + pmi22*pmi22);
870  |     }
871  | 
872  | #ifdef IMLIB2_EXISTS
873  |   /* Do the Imlib2 rendering */
874  |   imlib2_render_triangles
875  |     ((DATA32 *) thedisp->rgb, thedisp->rgb_width, thedisp->rgb_height, thesurf->num_triangles, screen_coords,
876  |      thesurf->vertices+9, 13, shades, 1);
877  | 
878  |   /* Imlib2 switches red and blue for some reason */
879  |   for (i=0; i<thedisp->rgb_height; i++)
880  |     {
881  |       int j;
882  |       for (j = i*thedisp->rgb_rowskip;
883  | 	   j < i*thedisp->rgb_rowskip + thedisp->rgb_width; j++)
884  | 	{
885  | 	  guchar tmp = thedisp->rgb [4*j];
886  | 	  thedisp->rgb[4*j] = thedisp->rgb[4*j+2];
887  | 	  thedisp->rgb[4*j+2] = tmp;
888  | 	}
889  |     }
890  | #else
891  |   SETERRQ (PETSC_ERR_SUP, "3-D rendering requires Imlib2");
892  | #endif
893  | 
894  |   i = PetscFree (shades); CHKERRQ (i);
895  |   i = PetscFree (screen_coords); CHKERRQ (i);
896  | }