1    | /***************************************
2    |   $Header: /cvsroot/petscgraphics/render.c,v 1.15 2004/11/30 14:35:54 hazelsct Exp $
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    | 
10   | /* Build with -DDEBUG for debugging output */
11   | #undef DPRINTF
12   | #ifdef DEBUG
13   | #define DPRINTF(fmt, args...) PetscPrintf (PETSC_COMM_WORLD, "%s: " fmt, __FUNCT__, args)
14   | #else
15   | #define DPRINTF(fmt, args...)
16   | #endif
17   | 
18   | 
19   | #undef __FUNCT__
20   | #define __FUNCT__ "pseudocolor"
21   | 
22   | /*++++++++++++++++++++++++++++++++++++++
23   |   This little function converts a scalar value into an rgb color from red to
24   |   blue.
25   | 
26   |   PetscScalar val Value to convert.
27   | 
28   |   PetscScalar* scale Array with minimum and maximum values in which to scale
29   |   val.
30   | 
31   |   guchar *pixel Address in rgb buffer where this pixel should be painted.
32   |   ++++++++++++++++++++++++++++++++++++++*/
33   | 
34   | static inline void pseudocolor (PetscScalar val, PetscScalar* scale,
35   | 				guchar *pixel)
36   | {
37   |   PetscScalar shade = (val - scale[0]) / (scale[1] - scale[0]);
38   |   /* Old stuff *
39   |   if (shade < 0.2) {      /* red -> yellow *
40   |     *pixel++ = 255;
41   |     *pixel++ = 1275*shade;
42   |     *pixel++ = 0; }
43   |   else if (shade < 0.4) { /* yellow -> green *
44   |     *pixel++ = 510-1275*shade;
45   |     *pixel++ = 255;
46   |     *pixel++ = 0; }
47   |   else if (shade < 0.6) { /* green -> cyan *
48   |     *pixel++ = 0;
49   |     *pixel++ = 255;
50   |     *pixel++ = 1275*shade-510; }
51   |   else if (shade < 0.8) { /* cyan -> blue *
52   |     *pixel++ = 0;
53   |     *pixel++ = 1020-1275*shade;
54   |     *pixel++ = 255; }
55   |   else {                  /* blue -> magenta *
56   |     *pixel++ = 1275*shade-1020;
57   |     *pixel++ = 0;
58   |     *pixel++ = 255; }
59   |   /* New stuff */
60   |   if (shade < 0.25) {      /* red -> yellow */
61   |     *pixel++ = 255;
62   |     *pixel++ = 1020*shade;
63   |     *pixel++ = 0; }
64   |   else if (shade < 0.5) { /* yellow -> green */
65   |     *pixel++ = 510-1020*shade;
66   |     *pixel++ = 255;
67   |     *pixel++ = 0; }
68   |   else if (shade < 0.75) { /* green -> cyan */
69   |     *pixel++ = 0;
70   |     *pixel++ = 255;
71   |     *pixel++ = 1020*shade-510; }
72   |   else { /* cyan -> blue */
73   |     *pixel++ = 0;
74   |     *pixel++ = 1020-1020*shade;
75   |     *pixel++ = 255; }
76   | }
77   | 
78   | 
79   | #undef __FUNCT__
80   | #define __FUNCT__ "pseudoternarycolor"
81   | 
82   | /*++++++++++++++++++++++++++++++++++++++
83   |   This little function converts two ternary fractions into an rgb color, with
84   |   yellow, cyan and magenta indicating the corners.
85   | 
86   |   PetscScalar A First ternary fraction.
87   | 
88   |   PetscScalar B Second ternary fraction.
89   | 
90   |   PetscScalar *scale Array first and second ternary fractions of each of the
91   |   three corner values for scaling.
92   | 
93   |   guchar *pixel  Address in rgb buffer where this pixel should be painted.
94   |   ++++++++++++++++++++++++++++++++++++++*/
95   | 
96   | static inline void pseudoternarycolor
97   | (PetscScalar A, PetscScalar B, PetscScalar *scale, guchar *pixel)
98   | {
99   |   PetscScalar x1, x2, inverse_det;
100  | 
101  |   /* Transform A,B into x1,x2 based on corners */
102  |   inverse_det = 1./((scale[2]-scale[0])*(scale[5]-scale[1]) -
103  | 		    (scale[4]-scale[0])*(scale[3]-scale[1]));
104  |   x1 = ((A-scale[0])*(scale[5]-scale[1]) -
105  | 	(B-scale[1])*(scale[4]-scale[0])) * inverse_det;
106  |   x2 = ((B-scale[1])*(scale[2]-scale[0]) -
107  | 	(A-scale[0])*(scale[3]-scale[1])) * inverse_det;
108  | 
109  |   /* Now colorize */
110  |   *pixel++ = 255*(1.-x1);
111  |   *pixel++ = 255*(x1+x2);
112  |   *pixel++ = 255*(1.-x2);
113  | }
114  | 
115  | 
116  | #undef __FUNCT__
117  | #define __FUNCT__ "pseudohueintcolor"
118  | 
119  | /*++++++++++++++++++++++++++++++++++++++
120  |   This little function converts a vector into an rgb color with hue indicating
121  |   direction (green, yellow, red, blue at 0, 90, 180, 270 degrees) and intensity
122  |   indicating magnitude relative to reference magnitude in scale[1].
123  | 
124  |   PetscScalar vx Vector's
125  |   +latex+$x$-component.
126  |   +html+ <i>x</i>-component.
127  | 
128  |   PetscScalar vy Vector's
129  |   +latex+$y$-component.
130  |   +html+ <i>y</i>-component.
131  | 
132  |   PetscScalar *scale Array whose second entry has the reference magnitude.
133  | 
134  |   guchar *pixel Address in rgb buffer where this pixel should be painted.
135  |   ++++++++++++++++++++++++++++++++++++++*/
136  | 
137  | static inline void pseudohueintcolor
138  | (PetscScalar vx, PetscScalar vy, PetscScalar *scale, guchar *pixel)
139  | {
140  |   PetscScalar mag=sqrt(vx*vx+vy*vy), theta=atan2(vy,vx), red, green, blue;
141  |   if (scale[1] <= 0.)
142  |     {
143  |       *pixel = *(pixel+1) = *(pixel+2) = 0;
144  |       return;
145  |     }
146  |   mag = (mag > scale[1]) ? 1.0 : mag/scale[1];
147  | 
148  |   red = 2./M_PI * ((theta<-M_PI/2.) ? -M_PI/2.-theta :
149  | 		   ((theta<0.) ? 0. : ((theta<M_PI/2.) ? theta : M_PI/2.)));
150  |   green = 2./M_PI * ((theta<-M_PI/2.) ? 0. :
151  | 		     ((theta<0.) ? theta+M_PI/2. :
152  | 		      ((theta<M_PI/2.) ? M_PI/2. : M_PI-theta)));
153  |   blue = 2./M_PI * ((theta<-M_PI/2.) ? theta+M_PI :
154  | 		    ((theta<0.) ? -theta : 0.));
155  | 
156  |   *pixel++ = 255*mag*red;
157  |   *pixel++ = 255*mag*green;
158  |   *pixel++ = 255*mag*blue;
159  | }
160  | 
161  | 
162  | #undef __FUNCT__
163  | #define __FUNCT__ "pseudoshearcolor"
164  | 
165  | /*++++++++++++++++++++++++++++++++++++++
166  |   This little function converts a shear tensor (symmetric, diagonals sum to
167  |   zero) into an rgb color with hue indicating direction of the tensile stress
168  |   (red, yellow, green, cyan, blue, magenta at 0, 30, 60, 90, 120, 150 degrees
169  |   respectively; 180 is equivalent to zero for stress) and intensity
170  |   indicating magnitude relative to reference magnitude in scale[0].
171  | 
172  |   PetscScalar gxx Tensor's
173  |   +latex+$xx$-component.
174  |   +html+ <i>xx</i>-component.
175  | 
176  |   PetscScalar gxy Tensor's
177  |   +latex+$xy$-component.
178  |   +html+ <i>xy</i>-component.
179  | 
180  |   PetscScalar *scale Array whose first entry has the reference magnitude.
181  | 
182  |   guchar *pixel Address in rgb buffer where this pixel should be painted.
183  |   ++++++++++++++++++++++++++++++++++++++*/
184  | 
185  | static inline void pseudoshearcolor
186  | (PetscScalar gxx, PetscScalar gxy, PetscScalar *scale, guchar *pixel)
187  | {
188  |   PetscScalar mag=sqrt(gxx*gxx+gxy*gxy), theta=atan2(gxy,gxx),
189  |     red,green,blue;
190  |   if (scale[0] <= 0.)
191  |     {
192  |       *pixel = *(pixel+1) = *(pixel+2) = 0;
193  |       return;
194  |     }
195  |   mag = (mag > scale[0]) ? 1.0 : mag/scale[0];
196  | 
197  |   red = 3./M_PI * ((theta < -2.*M_PI/3.) ? 0. :
198  | 		   ((theta < -M_PI/3.) ? theta + 2.*M_PI/3. :
199  | 		    ((theta < M_PI/3.) ? M_PI/3. :
200  | 		     ((theta < 2.*M_PI/3.) ? 2.*M_PI/3. - theta: 0.))));
201  |   green = 3./M_PI * ((theta < -M_PI/3.) ? M_PI/3. :
202  | 		     ((theta < 0.) ? -theta :
203  | 		      ((theta < 2.*M_PI/3.) ? 0. : theta - 2.*M_PI/3.)));
204  |   blue = 3./M_PI * ((theta < -2.*M_PI/3) ? -2.*M_PI/3. - theta :
205  | 		    ((theta < 0.) ? 0. :
206  | 		     ((theta < M_PI/3.) ? theta : M_PI/3.)));
207  | 
208  |   *pixel++ = 255*mag*red;
209  |   *pixel++ = 255*mag*green;
210  |   *pixel++ = 255*mag*blue;
211  | }
212  | 
213  | 
214  | #undef __FUNCT__
215  | #define __FUNCT__ "render_scale_2d"
216  | 
217  | /*++++++++++++++++++++++++++++++++++++++
218  |   This draws a little rwidth x rheight image depicting the scale of a field
219  |   variable.
220  | 
221  |   int render_scale_2d It returns zero or an error code.
222  | 
223  |   guchar *rgb Image to draw into.
224  | 
225  |   int rwidth Intended width of the image.
226  | 
227  |   int rheight Intended height of the image.
228  | 
229  |   int bytes_per_pixel Number of bytes per pixel in the image.
230  | 
231  |   field_plot_type fieldtype Type of plot.
232  | 
233  |   int symmetry Symmetry order for vector scale image.
234  |   ++++++++++++++++++++++++++++++++++++++*/
235  | 
236  | int render_scale_2d (guchar *rgb, int rwidth, int rheight, int bytes_per_pixel,
237  | 		     field_plot_type fieldtype, int symmetry)
238  | {
239  |   PetscScalar scale[6];
240  |   int i, j, myheight;
241  | 
242  |   switch (fieldtype)
243  |     {
244  |     case FIELD_SCALAR:
245  |       scale[0]=0.;
246  |       scale[1]=1.;
247  |       for (j=0; j<rheight; j++)
248  | 	for (i=0; i<rwidth; i++)
249  | 	  pseudocolor ((PetscScalar) i/(rwidth-1), scale,
250  | 		       rgb+(j*rwidth + i) * bytes_per_pixel);
251  |       return 0;
252  | 
253  |     case FIELD_TERNARY:
254  |       scale[0] = scale[1] = scale[3] = scale[4] = 0.;
255  |       scale[2] = scale[5] = 1.;
256  |       myheight = ((PetscScalar) rheight > 0.5*sqrt(3.) * rwidth) ?
257  | 	(int) (0.5*sqrt(3.)* rwidth) : rheight;
258  |       for (j=0; j<myheight; j++)
259  | 	for (i=0; i<(int) (2./sqrt(3.)*(myheight-j)); i++)
260  | 	  pseudoternarycolor
261  | 	    ((PetscScalar) i/(myheight*2./sqrt(3.)),
262  | 	     (PetscScalar) j/myheight, scale,
263  | 	     rgb + (((rheight+myheight)/2 -j)*rwidth +i +(int)(sqrt(3.)/3.*j))*
264  | 	     bytes_per_pixel);
265  |       return 0;
266  | 
267  |     case FIELD_VECTOR:
268  |       scale[0] = 0.;
269  |       scale[1] = 1.;
270  |       myheight = ((rheight > rwidth) ? rwidth : rheight)/2;
271  |       for (j=-myheight; j<myheight; j++)
272  | 	for (i=-(int)(myheight*sqrt(1.-((PetscScalar)j*j/myheight/myheight)));
273  | 	     i<(int)(myheight*sqrt(1.-((PetscScalar)j*j/myheight/myheight)));
274  | 	     i++)
275  | 	  {
276  | 	    PetscScalar mag = (PetscScalar) sqrt (i*i+j*j)/myheight,
277  | 	      theta = atan2 (j, i);
278  | 	    pseudohueintcolor
279  | 	      (mag * cos (symmetry*theta), mag * sin (symmetry*theta), scale,
280  | 	       rgb + ((myheight-j)*rwidth + rwidth/2+i)*bytes_per_pixel);
281  | 	  }
282  |       return 0;
283  | 
284  |     case FIELD_TENSOR_SHEAR:
285  |       scale[0] = 1.;
286  |       myheight = ((rheight > rwidth) ? rwidth : rheight)/2;
287  |       for (j=-myheight; j<myheight; j++)
288  | 	for (i=-(int)(myheight*sqrt(1.-((PetscScalar)j*j/myheight/myheight)));
289  | 	     i<(int)(myheight*sqrt(1.-((PetscScalar)j*j/myheight/myheight)));
290  | 	     i++)
291  | 	  {
292  | 	    PetscScalar mag = (PetscScalar) sqrt (i*i+j*j)/myheight,
293  | 	      theta = atan2 (j, i);
294  | 	    pseudoshearcolor
295  | 	      (mag * cos (2*theta), mag * sin (2*theta), scale,
296  | 	       rgb + ((myheight-j)*rwidth + rwidth/2+i)*bytes_per_pixel);
297  | 	  }
298  |       return 0;
299  |     }
300  |   SETERRQ (PETSC_ERR_ARG_OUTOFRANGE, "Field type not yet supported");
301  | }
302  | 
303  | 
304  | #undef __FUNCT__
305  | #define __FUNCT__ "render_rgb_local_2d"
306  | 
307  | /*++++++++++++++++++++++++++++++++++++++
308  |   Render data from global_array into local part of an RGB buffer.  When running
309  |   in parallel, these local buffers should be collected and layered to produce
310  |   the full image.
311  | 
312  |   int render_rgb_local_2d Returns zero or an error code.
313  | 
314  |   guchar *rgb RGB buffer, or location in the RGB buffer, in which to render.
315  | 
316  |   int rwidth Total width of the section of the RGB buffer to draw.
317  | 
318  |   int rheight Total height of the section of the RGB buffer to draw.
319  | 
320  |   int rowskip Total width of the RGB buffer.
321  | 
322  |   int bytes_per_pixel Number of bytes per pixel in this RGB buffer (typically 3
323  |   or 4).
324  | 
325  |   PetscScalar *global_array Local array of global vector data to render.
326  | 
327  |   int num_fields Number of field variables in the array.
328  | 
329  |   int display_field The (first) field we are rendering now.
330  | 
331  |   field_plot_type fieldtype The type of this field.
332  | 
333  |   PetscScalar *scale Array of minimum and maximum values to pass to the
334  |   various pseudocolor functions; if NULL, call auto_scale to determine those
335  |   values.
336  | 
337  |   int nx Width of the array.
338  | 
339  |   int ny Height of the array.
340  | 
341  |   int xs Starting
342  |   +latex+$x$-coordinate
343  |   +html+ <i>x</i>-coordinate
344  |   of the local part of the global vector.
345  | 
346  |   int ys Starting
347  |   +latex+$y$-coordinate
348  |   +html+ <i>y</i>-coordinate
349  |   of the local part of the global vector.
350  | 
351  |   int xm Width of the local part of the global vector.
352  | 
353  |   int ym Height of the local part of the global vector.
354  |   ++++++++++++++++++++++++++++++++++++++*/
355  | 
356  | int render_rgb_local_2d
357  | (guchar *rgb, int rwidth, int rheight, int rowskip, int bytes_per_pixel,
358  |  PetscScalar *global_array, int num_fields, int display_field,
359  |  field_plot_type fieldtype, PetscScalar *scale, int nx,int ny, int xs,int ys,
360  |  int xm,int ym)
361  | {
362  |   int ix,iy, ierr;
363  |   PetscScalar local_scale[6];
364  | 
365  |   /* Determine default min and max if none are provided */
366  |   if (scale == NULL)
367  |     {
368  |       scale = local_scale;
369  |       ierr = auto_scale (global_array, xm*ym, num_fields, display_field,
370  | 			   fieldtype, 2, scale); CHKERRQ (ierr);
371  |     }
372  |   DPRINTF ("Final scale: %g %g %g %g %g %g\n", scale[0],
373  | 	   scale[1], scale[2], scale[3], scale[4], scale[5]);
374  | 
375  |   /* Do the rendering (note switch in inner loop, gotta get rid of that) */
376  |   for (iy=rheight*ys/ny; iy<rheight*(ys+ym)/ny; iy++)
377  |     for (ix=rwidth*xs/nx; ix<rwidth*(xs+xm)/nx; ix++)
378  |       {
379  | 	int vecindex = (((rheight-iy-1)*ny/rheight)*nx +
380  | 			ix*nx/rwidth)*num_fields + display_field;
381  | 	guchar *pixel = rgb + bytes_per_pixel*(iy*rowskip + ix);
382  | 
383  | 	switch (fieldtype)
384  | 	  {
385  | 	  case FIELD_SCALAR:
386  | 	  case FIELD_SCALAR+1:
387  | 	    {
388  | 	      pseudocolor (global_array [vecindex], scale, pixel);
389  | 	      break;
390  | 	    }
391  | 	  case FIELD_TERNARY:
392  | 	    {
393  | 	      pseudoternarycolor
394  | 		(global_array [vecindex], global_array [vecindex+1], scale,
395  | 		 pixel);
396  | 	      break;
397  | 	    }
398  | 	  case FIELD_VECTOR:
399  | 	  case FIELD_VECTOR+1:
400  | 	    {
401  | 	      pseudohueintcolor
402  | 		(global_array [vecindex], global_array [vecindex+1], scale,
403  | 		 pixel);
404  | 	      break;
405  | 	    }
406  | 	  case FIELD_TENSOR_SHEAR:
407  | 	    {
408  | 	      pseudoshearcolor
409  | 		(global_array [vecindex], global_array [vecindex+1], scale,
410  | 		 pixel);
411  | 	      break;
412  | 	    }
413  | 	  default:
414  | 	    SETERRQ (PETSC_ERR_ARG_OUTOFRANGE, "Field type not yet supported");
415  | 	  }
416  |       }
417  |   return 0;
418  | }
419  | 
420  | 
421  | #undef __FUNCT__
422  | #define __FUNCT__ "render_rgb_local_3d"
423  | 
424  | /*++++++++++++++++++++++++++++++++++++++
425  |   Render triangle data into an RGB buffer.  When called in parallel, the
426  |   resulting images should be layered to give the complete picture.  Zooming is
427  |   done by adjusting the ratio of the dir vector to the right vector.
428  |   +latex+\par
429  |   +html+ <p>
430  |   The coordinate transformation is pretty simple.
431  |   +latex+A point $\vec{p}$ in space is transformed to $x,y$ on the screen by
432  |   +latex+representing it as:
433  |   +latex+\begin{equation}
434  |   +latex+  \label{eq:transform}
435  |   +latex+  \vec{p} = \vec{\imath} + a \vec{d} + ax \vec{r} + ay \vec{u},
436  |   +latex+\end{equation}
437  |   +latex+where $\vec{\imath}$ is the observer's point (passed in as {\tt eye}),
438  |   +latex+$\vec{d}$ is the direction of observation (passed in as {\tt dir}),
439  |   +latex+$\vec{r}$ is the rightward direction to the observer (passed in as
440  |   +latex+{\tt right}), and $\vec{u}$ is the upward direction to the observer
441  |   +latex+which is given the direction of $\vec{r}\times\vec{d}$ and the
442  |   +latex+magnitude of $\vec{r}$.
443  |   +html+ A point <u><i>p</i></u> in space is transformed to <i>x, y</i> on the
444  |   +html+ screen by representing it as:
445  |   +html+ <center><u><i>p</i></u> = <u><i>i</i></u> + <i>a <u>d</u></i> +
446  |   +html+ <i>ax <u>r</u></i> + <i>ay <u>u</u></i>,</center>
447  |   +html+ where <u><i>i</i></u> is the observer's point (passed in as
448  |   +html+ <tt>eye</tt>), <u><i>d</i></u> is the direction of observation (passed
449  |   +html+ in as <tt>dir</tt>), <u><i>r</i></u> is the rightward direction to the
450  |   +html+ observer (passed in as <tt>right</tt>), and <u><i>u</i></u> is the
451  |   +html+ upward direction to the observer which is given the direction of
452  |   +html+ the cross product of <u><i>d</i></u> and <u><i>r</i></u> and the
453  |   +html+ magnitude of <u><i>r</i></u>.
454  |   +latex+\par
455  |   +html+ <p>
456  |   +latex+This system in equation \ref{eq:transform} can easily be solved for
457  |   +latex+$x$ and $y$ by first making it into a matrix equation:
458  |   +latex+\begin{equation}
459  |   +latex+  \label{eq:matrix}
460  |   +latex+  \left(\begin{array}{ccc} d_x&r_x&u_x \\ d_y&r_y&u_y \\ d_z&r_z&u_z
461  |   +latex+    \end{array}\right)
462  |   +latex+  \left(\begin{array}{c} a \\ ax \\ ay \end{array}\right) =
463  |   +latex+  \left(\begin{array}{c} p_x-i_x \\ p_y-i_y \\ p_z-i_z
464  |   +latex+    \end{array}\right).
465  |   +latex+\end{equation}
466  |   +latex+Calling the matrix $M$ the inverse of the matrix on the left side of
467  |   +latex+equation \ref{eq:matrix} gives the result:
468  |   +latex+\begin{eqnarray}
469  |   +latex+  a&=& M_{00} (p_x-i_x) + M_{01} (p_y-i_y) + M_{02} (p_z-i_z), \\
470  |   +latex+  x&=& \frac{1}{a}[M_{10}(p_x-i_x)+M_{11}(p_y-i_y)+M_{12}(p_z-i_z)],\\
471  |   +latex+  y&=& \frac{1}{a}[M_{20}(p_x-i_x)+M_{21}(p_y-i_y)+M_{22}(p_z-i_z)].
472  |   +latex+\end{eqnarray}
473  |   +html+ This system can easily be solved for <i>x</i> and <i>y</i> by first
474  |   +html+ making it into a matrix equation:
475  |   +html+ <center>(<i><u>d</u> <u>r</u> <u>u</u></i>) (<i>a, ax, ay</i>) =
476  |   +html+ <u><i>p</i></u> - <u><i>i</i></u>.</center>
477  |   +html+ Calling the matrix <u><i>M</i></u> the inverse of the matrix on the
478  |   +html+ left side gives the result:
479  |   +html+ <center><table>
480  |   +html+ <tr><td><i>a</i></td><td>=</td><td>
481  |   +html+ <i>M</i><sub>00</sub> (<i>p<sub>x</sub></i>-<i>i<sub>x</sub></i>) +
482  |   +html+ <i>M</i><sub>01</sub> (<i>p<sub>y</sub></i>-<i>i<sub>y</sub></i>) +
483  |   +html+ <i>M</i><sub>02</sub> (<i>p<sub>z</sub></i>-<i>i<sub>z</sub></i>),
484  |   +html+ </td></tr><tr><td><i>x</i></td><td>=</td><td>
485  |   +html+ [<i>M</i><sub>10</sub> (<i>p<sub>x</sub></i>-<i>i<sub>x</sub></i>) +
486  |   +html+ <i>M</i><sub>11</sub> (<i>p<sub>y</sub></i>-<i>i<sub>y</sub></i>) +
487  |   +html+ <i>M</i><sub>12</sub> (<i>p<sub>z</sub></i>-<i>i<sub>z</sub></i>)]
488  |   +html+ / <i>a</i>,</td></tr><tr><td><i>y</i></td><td>=</td><td>
489  |   +html+ [<i>M</i><sub>00</sub> (<i>p<sub>x</sub></i>-<i>i<sub>x</sub></i>) +
490  |   +html+ <i>M</i><sub>01</sub> (<i>p<sub>y</sub></i>-<i>i<sub>y</sub></i>) +
491  |   +html+ <i>M</i><sub>02</sub> (<i>p<sub>z</sub></i>-<i>i<sub>z</sub></i>)]
492  |   +html+ / <i>a</i>.</td></tr></table></center>
493  |   +latex+\par
494  |   +html+ <p>
495  |   The triple product of the vector from the light source to the triangle
496  |   centroid and the two vectors making up the triangle edges determines the
497  |   cosine of the angle made by the triangle normal and the incident light, whose
498  |   absolute value is used here to shade the triangle.  At this point, the light
499  |   source is taken as the observer location at
500  |   +latex+``{\tt eye}'',
501  |   +html+ "<tt>eye</tt>",
502  |   but that can easily be modified to use one or more independent light sources.
503  | 
504  |   int render_rgb_local_3d Returns zero or an error code.
505  | 
506  |   guchar *rgb RGB buffer in which to render.
507  | 
508  |   int rwidth Total width of the RGB buffer.
509  | 
510  |   int rheight Total height of the RGB buffer.
511  | 
512  |   int bytes_per_pixel Number of bytes per pixel in this RGB buffer (typically 3
513  |   or 4).
514  | 
515  |   int num_triangles Number of triangles to render.
516  | 
517  |   PetscScalar *vertices Table of coordinates (x1,y1,z1, x2,y2,z2, x3,y3,z3) and
518  |   colors (RGBA 0-1) making thirteen values per triangle.
519  | 
520  |   PetscScalar *eye Point from where we're looking (x,y,z).
521  | 
522  |   PetscScalar *dir Direction we're looking (x,y,z).
523  | 
524  |   PetscScalar *right Rightward direction in physical space (x,y,z).
525  |   ++++++++++++++++++++++++++++++++++++++*/
526  | 
527  | int render_rgb_local_3d
528  | (guchar *rgb, int rwidth, int rheight, int bytes_per_pixel,
529  |  int num_triangles, PetscScalar *vertices,
530  |  PetscScalar *eye, PetscScalar *dir, PetscScalar *right)
531  | {
532  | #ifdef IMLIB2_EXISTS
533  |   int *screen_coords, i;
534  |   PetscScalar *shades, up[3], m00,m01,m02,m10,m11,m12,m20,m21,m22, a;
535  | 
536  |   if (bytes_per_pixel != 4)
537  |     SETERRQ (PETSC_ERR_ARG_SIZ, "Imlib2 rendering requires 4 bytes per pixel");
538  | 
539  |   /* Allocate memory for the shades and integer coordinates */
540  |   i = PetscMalloc (6*num_triangles * sizeof(int), &screen_coords); CHKERRQ (i);
541  |   i = PetscMalloc (num_triangles * sizeof(PetscScalar), &shades); CHKERRQ (i);
542  | 
543  |   /* Calculate the "up" vector as described above, storing the magnitude of dir
544  |      in a */
545  |   a = sqrt (dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]);
546  |   up[0] = (right[1]*dir[2] - right[2]*dir[1]) / a;
547  |   up[1] = (right[2]*dir[0] - right[0]*dir[2]) / a;
548  |   up[2] = (right[0]*dir[1] - right[1]*dir[0]) / a;
549  | 
550  |   /* Calculate the inverse of the dir, right, up columnar matrix, storing the
551  |      determinant of the matrix in a */
552  |   a = (dir[0] * (right[1]*up[2] - right[2]*up[1]) +
553  |        dir[1] * (right[2]*up[0] - right[0]*up[2]) +
554  |        dir[2] * (right[0]*up[1] - right[1]*up[0]));
555  |   m00 = (right[1]*up[2]   - right[2]*up[1])    / a;
556  |   m01 = (right[2]*up[0]   - right[0]*up[2])    / a;
557  |   m02 = (right[0]*up[1]   - right[1]*up[0])    / a;
558  |   m10 = (up[1]   *dir[2]  - up[2]   *dir[1])   / a;
559  |   m11 = (up[2]   *dir[0]  - up[0]   *dir[2])   / a;
560  |   m12 = (up[0]   *dir[1]  - up[1]   *dir[0])   / a;
561  |   m20 = (dir[1] *right[2] - dir[2]  *right[1]) / a;
562  |   m21 = (dir[2] *right[0] - dir[0]  *right[2]) / a;
563  |   m22 = (dir[0] *right[1] - dir[1]  *right[0]) / a;
564  | 
565  |   /* Loop over triangles, transforming their coordinates and setting shading */
566  |   for (i=0; i<num_triangles; i++)
567  |     {
568  |       PetscScalar pmi00,pmi01,pmi02,pmi10,pmi11,pmi12,pmi20,pmi21,pmi22, c[3];
569  | 
570  |       /* Calculate coordinates relative to observer p-i (hence pmi) */
571  |       pmi00 = vertices[13*i]   - eye[0];
572  |       pmi01 = vertices[13*i+1] - eye[1];
573  |       pmi02 = vertices[13*i+2] - eye[2];
574  |       pmi10 = vertices[13*i+3] - eye[0];
575  |       pmi11 = vertices[13*i+4] - eye[1];
576  |       pmi12 = vertices[13*i+5] - eye[2];
577  |       pmi20 = vertices[13*i+6] - eye[0];
578  |       pmi21 = vertices[13*i+7] - eye[1];
579  |       pmi22 = vertices[13*i+8] - eye[2];
580  | 
581  |       /* Calculate x,y for point 0 */
582  |       a = m00*pmi00 + m01*pmi01 + m02*pmi02;
583  |       if (a < 0)
584  | 	screen_coords [6*i+0] = screen_coords [6*i+1] = -1;
585  |       else
586  | 	{
587  | 	  screen_coords [6*i+0] = rwidth/2 *
588  | 	    (1. + (m10*pmi00 + m11*pmi01 + m12*pmi02) / a);
589  | 	  screen_coords [6*i+1] = rheight/2 - rwidth/2 *
590  | 	    (m20*pmi00 + m21*pmi01 + m22*pmi02) / a;
591  | 	}
592  | 
593  |       /* Calculate x,y for point 1 */
594  |       a = m00*pmi10 + m01*pmi11 + m02*pmi12;
595  |       if (a <= 0)
596  | 	screen_coords [6*i+2] = screen_coords [6*i+3] = -1;
597  |       else
598  | 	{
599  | 	  screen_coords [6*i+2] = rwidth/2 *
600  | 	    (1. + (m10*pmi10 + m11*pmi11 + m12*pmi12) / a);
601  | 	  screen_coords [6*i+3] = rheight/2 - rwidth/2 *
602  | 	    (m20*pmi10 + m21*pmi11 + m22*pmi12) / a;
603  | 	}
604  | 
605  |       /* Calculate x,y for point 2 */
606  |       a = m00*pmi20 + m01*pmi21 + m02*pmi22;
607  |       if (a <= 0)
608  | 	screen_coords [6*i+4] = screen_coords [6*i+5] = -1;
609  |       else
610  | 	{
611  | 	  screen_coords [6*i+4] = rwidth/2 *
612  | 	    (1. + (m10*pmi20 + m11*pmi21 + m12*pmi22) / a);
613  | 	  screen_coords [6*i+5] = rheight/2 - rwidth/2 *
614  | 	    (m20*pmi20 + m21*pmi21 + m22*pmi22) / a;
615  | 	}
616  | 
617  | #ifdef DEBUG
618  |       if (i==0)
619  | 	DPRINTF ("First triangle: (%g,%g,%g), (%g,%g,%g), (%g,%g,%g);\n"
620  | 		 "2-D coordinates: (%d,%d), (%d,%d), (%d,%d)\n",
621  | 		 vertices[0],vertices[1],vertices[2],vertices[3],vertices[4],
622  | 		 vertices[5],vertices[6],vertices[7],vertices[8],
623  | 		 screen_coords[0],screen_coords[1],screen_coords[3],
624  | 		 screen_coords[3],screen_coords[4],screen_coords[5]);
625  | #endif
626  | 
627  |       /* Calculate relative triangle centroid, triangle arms, and shade */
628  |       c[0] = (pmi00 + pmi10 + pmi20) / 3.;
629  |       c[1] = (pmi01 + pmi11 + pmi21) / 3.;
630  |       c[2] = (pmi02 + pmi12 + pmi22) / 3.;
631  |       pmi10 -= pmi00;   pmi11 -= pmi01;   pmi12 -= pmi02;
632  |       pmi20 -= pmi00;   pmi21 -= pmi01;   pmi22 -= pmi02;
633  |       shades [i] = PetscAbsScalar (c[0] * (pmi11*pmi22 - pmi12*pmi21) +
634  | 				   c[1] * (pmi12*pmi20 - pmi10*pmi22) +
635  | 				   c[2] * (pmi10*pmi21 - pmi11*pmi20)) /
636  | 	sqrt (c[0]*c[0] + c[1]*c[1] + c[2]*c[2]) /
637  | 	sqrt (pmi10*pmi10 + pmi11*pmi11 + pmi12*pmi12) /
638  | 	sqrt (pmi20*pmi20 + pmi21*pmi21 + pmi22*pmi22);
639  |     }
640  | 
641  |   /* Do the Imlib2 rendering */
642  |   imlib2_render_triangles ((DATA32 *) rgb, rwidth, rheight, num_triangles,
643  | 			   screen_coords, vertices+9, 13, shades, 1);
644  | 
645  |   i = PetscFree (shades); CHKERRQ (i);
646  |   i = PetscFree (screen_coords); CHKERRQ (i);
647  | #else
648  |   SETERRQ (PETSC_ERR_SUP, "3-D rendering requires Imlib2");
649  | #endif
650  | }