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