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