1 #include <omp.h> 2 3 /*************************************************************************/ 4 /* */ 5 /* Copyright (c) 1994 Stanford University */ 6 /* */ 7 /* All rights reserved. */ 8 /* */ 9 /* Permission is given to use, copy, and modify this software for any */ 10 /* non-commercial purpose as long as this copyright notice is not */ 11 /* removed. All other uses, including redistribution in whole or in */ 12 /* part, are forbidden without prior written permission. */ 13 /* */ 14 /* This software is provided with absolutely no warranty and no */ 15 /* support. */ 16 /* */ 17 /*************************************************************************/ 18 19 /************************************************************************* 20 * * 21 * adaptive.c: Render dataset via raytracing. * 22 * * 23 *************************************************************************/ 24 25 #include "incl.h" 26 27 float invjacobian[NM][NM]; /* Jacobian matrix showing object space */ 28 /* d{x,y,z} per image space d{x,y,z} */ 29 /* [0][0] is dx(object)/dx(image) */ 30 /* [0][2] is dz(object)/dx(image) */ 31 /* [2][0] is dx(object)/dz(image) */ 32 float invinvjacobian[NM][NM]; /* [i][j] = 1.0 / invjacobian[i][j] */ 33 /* For gathering statistics: */ 34 int num_rays_traced; /* number of calls to Trace_Ray */ 35 int num_traced_rays_hit_volume; /* number of traced rays that hit volume */ 36 int num_samples_trilirped; /* number of samples trilirped */ 37 int itest; 38 39 #define RAY_TRACED ((MAX_PIXEL+1)/2) /* Ray traced at this pixel */ 40 #define START_RAY 1 41 #define INTERPOLATED ((MAX_PIXEL+1)/32) /* This pixel interpolated */ 42 43 44 45 #include "anl.h" 46 47 void Ray_Trace(int my_node) { 48 int outx, outy, outz; 49 int i, j; 50 unsigned int starttime, stoptime, exectime, exectime1; 51 int pid; 52 char cmd[FILENAME_STRING_SIZE]; 53 54 /* Assumptions made by ray tracer: */ 55 /* o Frustrum clipping is performed. */ 56 /* All viewing frustums will be handled correctly. */ 57 /* o Contributions are obtained only from nearest 8 neighbors. */ 58 /* If downsizing was specified, some input voxels will be */ 59 /* unsampled, but upsizing may be specified and will be */ 60 /* handled correctly. */ 61 62 /* Compute inverse Jacobian matrix from */ 63 /* coordinates of output map unit voxel in object space, */ 64 /* then make a copy of object space d{x,y,z} per image space dz, */ 65 /* which controls number of ray samples per object space voxel */ 66 /* (i.e. Z-sampling rate), to allow recomputation per region. */ 67 for (i = 0; i < NM; i++) { 68 invjacobian[X][i] = uout_invvertex[0][0][1][i] 69 - uout_invvertex[0][0][0][i]; 70 invjacobian[Y][i] = uout_invvertex[0][1][0][i] 71 - uout_invvertex[0][0][0][i]; 72 invjacobian[Z][i] = uout_invvertex[1][0][0][i] 73 - uout_invvertex[0][0][0][i]; 74 } 75 76 /* Compute multiplicative inverse of inverse Jacobian matrix. */ 77 /* If any Jacobian is zero, compute no inverse for that element. */ 78 /* This test must be repeated before access to any inverse element. */ 79 for (i = 0; i < NM; i++) { 80 for (j = 0; j < NM; j++) { 81 if (ABS(invjacobian[i][j]) > SMALL) 82 invinvjacobian[i][j] = 1.0 / invjacobian[i][j]; 83 } 84 } 85 86 num_rays_traced = 0; 87 num_traced_rays_hit_volume = 0; 88 num_samples_trilirped = 0; 89 90 /* Invoke adaptive or non-adaptive ray tracer */ 91 if (adaptive) { 92 93 {;}; 94 95 {long time(); (starttime) = time(0);}; 96 Pre_Shade(my_node); 97 98 {;}; 99 Global->Counter--; 100 {;}; 101 while (Global->Counter) 102 ; 103 104 Ray_Trace_Adaptively(my_node); 105 106 {long time(); (stoptime) = time(0);}; 107 108 {;}; 109 110 mclock(stoptime, starttime, &exectime); 111 112 /* If adaptively ray tracing and highest sampling size is greater */ 113 /* than lowest size for volume data if polygon list exists or */ 114 /* display pixel size if it does not, recursively interpolate to */ 115 /* fill in any missing samples down to lowest size for volume data. */ 116 117 if (highest_sampling_boxlen > 1) { 118 119 {;}; 120 121 {long time(); (starttime) = time(0);}; 122 Interpolate_Recursively(my_node); 123 124 {long time(); (stoptime) = time(0);}; 125 126 {;}; 127 128 mclock(stoptime, starttime, &exectime1); 129 } 130 } else { 131 132 {;}; 133 134 {long time(); (starttime) = time(0);}; 135 136 Pre_Shade(my_node); 137 138 {;}; 139 Global->Counter--; 140 {;}; 141 while (Global->Counter) 142 ; 143 144 Ray_Trace_Non_Adaptively(my_node); 145 146 {long time(); (stoptime) = time(0);}; 147 148 {;}; 149 150 mclock(stoptime, starttime, &exectime); 151 exectime1 = 0; 152 } 153 154 {;}; 155 printf("%3d\t%3d\t%6u\t%6u\t%6d\t%6d\t%8d\n", my_node, frame, exectime, 156 exectime1, num_rays_traced, num_traced_rays_hit_volume, 157 num_samples_trilirped); 158 159 {;}; 160 161 {;}; 162 } 163 164 void Ray_Trace_Adaptively(int my_node) 165 { 166 int i,outx,outy,yindex,xindex; 167 168 int num_xqueue,num_yqueue,num_queue,lnum_xblocks,lnum_yblocks,lnum_blocks; 169 int xstart,xstop,ystart,ystop,local_node,work; 170 171 itest = 0; 172 173 num_xqueue = ROUNDUP((float)image_len[X]/(float)image_section[X]); 174 num_yqueue = ROUNDUP((float)image_len[Y]/(float)image_section[Y]); 175 num_queue = num_xqueue * num_yqueue; 176 lnum_xblocks = ROUNDUP((float)num_xqueue/(float)block_xlen); 177 lnum_yblocks = ROUNDUP((float)num_yqueue/(float)block_ylen); 178 lnum_blocks = lnum_xblocks * lnum_yblocks; 179 local_node = my_node; 180 Global->Queue[local_node][0] = 0; 181 while (Global->Queue[num_nodes][0] > 0) { 182 xstart = (local_node % image_section[X]) * num_xqueue; 183 xstart = ROUNDUP((float)xstart/(float)highest_sampling_boxlen); 184 xstart = xstart * highest_sampling_boxlen; 185 xstop = MIN(xstart+num_xqueue,image_len[X]); 186 ystart = (local_node / image_section[X]) * num_yqueue; 187 ystart = ROUNDUP((float)ystart/(float)highest_sampling_boxlen); 188 ystart = ystart * highest_sampling_boxlen; 189 ystop = MIN(ystart+num_yqueue,image_len[Y]); 190 {;}; 191 work = Global->Queue[local_node][0]; 192 Global->Queue[local_node][0] += 1; 193 {;}; 194 while (work < lnum_blocks) { 195 xindex = xstart + (work%lnum_xblocks)*block_xlen; 196 yindex = ystart + (work/lnum_xblocks)*block_ylen; 197 for (outy=yindex; outy<yindex+block_ylen && outy<ystop; 198 outy+=highest_sampling_boxlen) { 199 for (outx=xindex; outx<xindex+block_xlen && outx<xstop; 200 outx+=highest_sampling_boxlen) { 201 202 /* Trace rays within square box of highest sampling size */ 203 /* whose lower-left corner is current image space location. */ 204 Ray_Trace_Adaptive_Box(outx,outy,highest_sampling_boxlen); 205 } 206 } 207 {;}; 208 work = Global->Queue[local_node][0]; 209 Global->Queue[local_node][0] += 1; 210 {;}; 211 } 212 if (my_node == local_node) { 213 {;}; 214 Global->Queue[num_nodes][0]--; 215 {;}; 216 } 217 local_node = (local_node+1)%num_nodes; 218 while (Global->Queue[local_node][0] >= lnum_blocks && 219 Global->Queue[num_nodes][0] > 0) 220 local_node = (local_node+1)%num_nodes; 221 } 222 223 } 224 225 void Ray_Trace_Adaptive_Box(outx, outy, boxlen) 226 int outx, outy, boxlen; 227 { 228 int i,j; 229 int half_boxlen; 230 int min_volume_color,max_volume_color; 231 float foutx,fouty; 232 volatile int imask; 233 234 PIXEL *pixel_address; 235 236 /* Trace rays from all four corners of the box into the map, */ 237 /* being careful not to exceed the boundaries of the output image, */ 238 /* and using a flag array to avoid retracing any rays. */ 239 /* For diagnostic display, flag is set to a light gray. */ 240 /* If mipmapping, flag is light gray minus current mipmap level. */ 241 /* If mipmapping and ray has already been traced, */ 242 /* retrace it if current mipmap level is lower than */ 243 /* mipmap level in effect when ray was last traced, */ 244 /* thus replacing crude approximation with better one. */ 245 /* Meanwhile, compute minimum and maximum geometry/volume colors. */ 246 /* If polygon list exists, compute geometry-only colors */ 247 /* and volume-attenuated geometry-only colors as well. */ 248 /* If stochastic sampling and box is smaller than a display pixel, */ 249 /* distribute the rays uniformly across a square centered on the */ 250 /* nominal ray location and of size equal to the image array spacing.*/ 251 /* This scheme interpolates the jitter size / sample spacing ratio */ 252 /* from zero at one sample per display pixel, avoiding jitter noise, */ 253 /* to one at the maximum sampling rate, insuring complete coverage, */ 254 /* all the while building on previously traced rays where possible. */ 255 /* The constant radius also prevents overlap of jitter squares from */ 256 /* successive subdivision levels, preventing sample clumping noise. */ 257 258 min_volume_color = MAX_PIXEL; 259 max_volume_color = MIN_PIXEL; 260 261 for (i=0; i<=boxlen && outy+i<image_len[Y]; i+=boxlen) { 262 for (j=0; j<=boxlen && outx+j<image_len[X]; j+=boxlen) { 263 264 /*reschedule processes here if rescheduling only at synch points on simulator*/ 265 266 if (MASK_IMAGE(outy+i,outx+j) == 0) { 267 268 /*reschedule processes here if rescheduling only at synch points on simulator*/ 269 270 MASK_IMAGE(outy+i,outx+j) = START_RAY; 271 272 /*reschedule processes here if rescheduling only at synch points on simulator*/ 273 274 foutx = (float)(outx+j); 275 fouty = (float)(outy+i); 276 pixel_address = IMAGE_ADDRESS(outy+i,outx+j); 277 278 /*reschedule processes here if rescheduling only at synch points on simulator*/ 279 280 Trace_Ray(outx+j,outy+i,foutx,fouty,pixel_address); 281 282 /*reschedule processes here if rescheduling only at synch points on simulator*/ 283 284 MASK_IMAGE(outy+i,outx+j) = RAY_TRACED; 285 } 286 } 287 } 288 for (i=0; i<=boxlen && outy+i<image_len[Y]; i+=boxlen) { 289 for (j=0; j<=boxlen && outx+j<image_len[X]; j+=boxlen) { 290 imask = MASK_IMAGE(outy+i,outx+j); 291 292 /*reschedule processes here if rescheduling only at synch points on simulator*/ 293 294 while (imask == START_RAY) { 295 296 /*reschedule processes here if rescheduling only at synch points on simulator*/ 297 298 imask = MASK_IMAGE(outy+i,outx+j); 299 300 /*reschedule processes here if rescheduling only at synch points on simulator*/ 301 302 } 303 min_volume_color = MIN(IMAGE(outy+i,outx+j),min_volume_color); 304 max_volume_color = MAX(IMAGE(outy+i,outx+j),max_volume_color); 305 } 306 } 307 308 /* If size of current box is above lowest size for volume data and */ 309 /* magnitude of geometry/volume color difference is significant, or */ 310 /* size of current box is above lowest size for geometric data and */ 311 /* magnitudes of geometry-only and volume-attenuated geometry-only */ 312 /* are both significant, thus detecting only visible geometry events,*/ 313 /* invoke this function recursively to trace rays within the */ 314 /* four equal-sized square sub-boxes enclosed by the current box, */ 315 /* being careful not to exceed the boundaries of the output image. */ 316 /* Use of geometry-only color difference suppressed in accordance */ 317 /* with hybrid.trf as published in IEEE CG&A, March, 1990. */ 318 if (boxlen > lowest_volume_boxlen && 319 max_volume_color - min_volume_color >= 320 volume_color_difference) { 321 half_boxlen = boxlen >> 1; 322 for (i=0; i<boxlen && outy+i<image_len[Y]; i+=half_boxlen) { 323 for (j=0; j<boxlen && outx+j<image_len[X]; j+=half_boxlen) { 324 Ray_Trace_Adaptive_Box(outx+j,outy+i,half_boxlen); 325 } 326 } 327 } 328 } 329 330 void Ray_Trace_Non_Adaptively(int my_node) 331 { 332 int i,outx,outy,xindex,yindex; 333 float foutx,fouty; 334 PIXEL *pixel_address; 335 336 int num_xqueue,num_yqueue,num_queue,lnum_xblocks,lnum_yblocks,lnum_blocks; 337 int xstart,xstop,ystart,ystop,local_node,work; 338 339 num_xqueue = ROUNDUP((float)image_len[X]/(float)image_section[X]); 340 num_yqueue = ROUNDUP((float)image_len[Y]/(float)image_section[Y]); 341 num_queue = num_xqueue * num_yqueue; 342 lnum_xblocks = ROUNDUP((float)num_xqueue/(float)block_xlen); 343 lnum_yblocks = ROUNDUP((float)num_yqueue/(float)block_ylen); 344 lnum_blocks = lnum_xblocks * lnum_yblocks; 345 local_node = my_node; 346 Global->Queue[local_node][0] = 0; 347 while (Global->Queue[num_nodes][0] > 0) { 348 xstart = (local_node % image_section[X]) * num_xqueue; 349 xstop = MIN(xstart+num_xqueue,image_len[X]); 350 ystart = (local_node / image_section[X]) * num_yqueue; 351 ystop = MIN(ystart+num_yqueue,image_len[Y]); 352 {;}; 353 work = Global->Queue[local_node][0]++; 354 {;}; 355 while (work < lnum_blocks) { 356 int my; 357 xindex = xstart + (work%lnum_xblocks)*block_xlen; 358 yindex = ystart + (work/lnum_xblocks)*block_ylen; 359 my = ystop < yindex + block_ylen? ystop :yindex + block_ylen; 360 #pragma omp parallel for private (outx, fouty, foutx, pixel_address) schedule(dynamic) 361 for (outy=yindex; outy<my; outy++) { 362 int mx = xstop < xindex+block_xlen? xstop : xindex+block_xlen; 363 for (outx=xindex; outx<mx; outx++) { 364 365 /* Trace ray from specified image space location into map. */ 366 /* Stochastic sampling is as described in adaptive code. */ 367 foutx = (float)(outx); 368 fouty = (float)(outy); 369 pixel_address = IMAGE_ADDRESS(outy,outx); 370 Trace_Ray(outx,outy,foutx,fouty,pixel_address); 371 } 372 } 373 {;}; 374 work = Global->Queue[local_node][0]++; 375 {;}; 376 } 377 if (my_node == local_node) { 378 {;}; 379 Global->Queue[num_nodes][0]--; 380 {;}; 381 } 382 local_node = (local_node+1)%num_nodes; 383 while (Global->Queue[local_node][0] >= lnum_blocks && 384 Global->Queue[num_nodes][0] > 0) 385 local_node = (local_node+1)%num_nodes; 386 } 387 } 388 389 void Ray_Trace_Fast_Non_Adaptively(int my_node) 390 { 391 int i,outx,outy,xindex,yindex; 392 float foutx,fouty; 393 PIXEL *pixel_address; 394 395 for (i=0; i<num_blocks; i+=num_nodes) { 396 yindex = ((my_node+i)/num_xblocks)*block_ylen; 397 xindex = ((my_node+i)%num_xblocks)*block_xlen; 398 399 for (outy=yindex; outy<yindex+block_ylen && 400 outy<image_len[Y]; outy+=lowest_volume_boxlen) { 401 for (outx=xindex; outx<xindex+block_xlen && 402 outx<image_len[X]; outx+=lowest_volume_boxlen) { 403 404 /* Trace ray from specified image space location into map. */ 405 /* Stochastic sampling is as described in adaptive code. */ 406 MASK_IMAGE(outy,outx) += RAY_TRACED; 407 foutx = (float)(outx); 408 fouty = (float)(outy); 409 pixel_address = IMAGE_ADDRESS(outy,outx); 410 Trace_Ray(outx,outy,foutx,fouty,pixel_address); 411 num_rays_traced++; 412 } 413 } 414 } 415 } 416 417 void Interpolate_Recursively(int my_node) 418 { 419 int i,outx,outy,xindex,yindex; 420 421 for (i=0; i<num_blocks; i+=num_nodes) { 422 yindex = ((my_node+i)/num_xblocks)*block_ylen; 423 xindex = ((my_node+i)%num_xblocks)*block_xlen; 424 425 for (outy=yindex; outy<yindex+block_ylen && 426 outy<image_len[Y]; outy+=highest_sampling_boxlen) { 427 for (outx=xindex; outx<xindex+block_xlen && 428 outx<image_len[X]; outx+=highest_sampling_boxlen) { 429 430 /* Fill in image within square box of highest sampling size */ 431 /* whose lower-left corner is current image space location. */ 432 Interpolate_Recursive_Box(outx,outy,highest_sampling_boxlen); 433 } 434 } 435 } 436 } 437 438 void Interpolate_Recursive_Box(outx, outy, boxlen) 439 int outx, outy, boxlen; 440 { 441 int i,j; 442 int half_boxlen; 443 int corner_color[2][2],color; 444 int outx_plus_boxlen,outy_plus_boxlen; 445 446 float one_over_boxlen; 447 float xalpha,yalpha; 448 float one_minus_xalpha,one_minus_yalpha; 449 450 /* Fill in the four pixels at the midpoints of the sides and at */ 451 /* the center of the box by bilirping between the four corners, */ 452 /* being careful not to exceed the boundaries of the output image, */ 453 /* and using a flag array to avoid recomputing any pixels. */ 454 /* For diagnostic display, flag is set to a dark gray. */ 455 /* By making interpolation follow ray tracing, no pixels along the */ 456 /* perimeter are filled in using interpolation, just to be replaced */ 457 /* by a more accurate color when adjacent boxes are ray traced. */ 458 /* By making the interpolation recursive, it is able to fill in */ 459 /* large boxes using all pixels available along their perimeters, */ 460 /* rather than just at their four corners. This prevents "creases", */ 461 /* or, stated another way, insures zero-order continuity everywhere. */ 462 half_boxlen = boxlen >> 1; 463 one_over_boxlen = 1.0 / (float)boxlen; 464 465 outx_plus_boxlen = outx+boxlen < image_len[X] ? outx+boxlen : outx; 466 outy_plus_boxlen = outy+boxlen < image_len[Y] ? outy+boxlen : outy; 467 468 corner_color[0][0] = IMAGE(outy,outx); 469 corner_color[0][1] = IMAGE(outy,outx_plus_boxlen); 470 corner_color[1][0] = IMAGE(outy_plus_boxlen,outx); 471 corner_color[1][1] = IMAGE(outy_plus_boxlen,outx_plus_boxlen); 472 for (i=0; i<=boxlen && outy+i<image_len[Y]; i+=half_boxlen) { 473 yalpha = (float)i * one_over_boxlen; 474 one_minus_yalpha = 1.0 - yalpha; 475 for (j=0; j<=boxlen && outx+j<image_len[X]; j+=half_boxlen) { 476 xalpha = (float)j * one_over_boxlen; 477 one_minus_xalpha = 1.0 - xalpha; 478 if (MASK_IMAGE(outy+i,outx+j) == 0) { 479 color = corner_color[0][0]*one_minus_xalpha*one_minus_yalpha+ 480 corner_color[0][1]*xalpha*one_minus_yalpha+ 481 corner_color[1][0]*one_minus_xalpha*yalpha+ 482 corner_color[1][1]*xalpha*yalpha; 483 IMAGE(outy+i,outx+j) = color; 484 MASK_IMAGE(outy+i,outx+j) += INTERPOLATED; 485 } 486 } 487 } 488 489 /* If size of sub-boxes is above lowest size for volume data */ 490 /* if polygon list exists or display pixel size if it does not, */ 491 /* invoke this function recursively to fill in image within the */ 492 /* four equal-sized square sub-boxes enclosed by the current box, */ 493 /* being careful not to exceed the boundaries of the output image. */ 494 if (half_boxlen > 1) { 495 for (i=0; i<boxlen && outy+i<image_len[Y]; i+=half_boxlen) { 496 for (j=0; j<boxlen && outx+j<image_len[X]; j+=half_boxlen) { 497 Interpolate_Recursive_Box(outx+j,outy+i,half_boxlen); 498 } 499 } 500 } 501 } 502