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