1 /*
2 * vp_renderA.c
3 *
4 * Shear-warp volume rendering algorithm for affine view transformations.
5 *
6 * Copyright (c) 1994 The Board of Trustees of The Leland Stanford
7 * Junior University. All rights reserved.
8 *
9 * Permission to use, copy, modify and distribute this software and its
10 * documentation for any purpose is hereby granted without fee, provided
11 * that the above copyright notice and this permission notice appear in
12 * all copies of this software and that you do not sell the software.
13 * Commercial licensing is available by contacting the author.
14 *
15 * THE SOFTWARE IS PROVIDED "AS IS" AND WITHOUT WARRANTY OF ANY KIND,
16 * EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
17 * WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
18 *
19 * Author:
20 * Phil Lacroute
21 * Computer Systems Laboratory
22 * Electrical Engineering Dept.
23 * Stanford University
24 */
25
26 /*
27 * $Date: 1994/12/30 23:52:38 $
28 * $Revision: 1.23 $
29 */
30
31 #include "vp_global.h"
32
33 /*#define DUMP_SHADOW_VOLUME*/
34 /*#define DUMP_GRAY_VOLUME*/
35
36 extern void VPCompAC00G ANSI_ARGS((vpContext *vpc, int icount, int jcount,
37 int k, double slice_depth_cueing_dbl, GrayIntPixel *intimage,
38 double weightTLdbl, double weightBLdbl, double weightTRdbl,
39 double weightBRdbl, unsigned char *run_lengths, void *voxel_data));
40 extern void VPCompAR00G ANSI_ARGS((vpContext *vpc, int icount, int jcount,
41 int k, double slice_depth_cueing_dbl, GrayIntPixel *intimage,
42 double weightTLdbl, double weightBLdbl, double weightTRdbl,
43 double weightBRdbl, void *voxel_data, int voxel_istride,
44 int voxel_jstride));
45 extern void VPWarpA101N ANSI_ARGS((GrayIntPixel *in_image, int in_width,
46 int in_height, int in_bytes_per_scan, unsigned char *out_image,
47 int out_width, int out_height, int out_bytes_per_scan,
48 vpMatrix3 warp_matrix));
49 extern void VPWarpA110N ANSI_ARGS((GrayIntPixel *in_image, int in_width,
50 int in_height, int in_bytes_per_scan, unsigned char *out_image,
51 int out_width, int out_height, int out_bytes_per_scan,
52 vpMatrix3 warp_matrix));
53 extern void VPWarpA111N ANSI_ARGS((GrayIntPixel *in_image, int in_width,
54 int in_height, int in_bytes_per_scan, unsigned char *out_image,
55 int out_width, int out_height, int out_bytes_per_scan,
56 vpMatrix3 warp_matrix));
57 extern void VPWarpA301N ANSI_ARGS((RGBIntPixel *in_image, int in_width,
58 int in_height, int in_bytes_per_scan, unsigned char *out_image,
59 int out_width, int out_height, int out_bytes_per_scan,
60 vpMatrix3 warp_matrix));
61 extern void VPWarpA330N ANSI_ARGS((RGBIntPixel *in_image, int in_width,
62 int in_height, int in_bytes_per_scan, unsigned char *out_image,
63 int out_width, int out_height, int out_bytes_per_scan,
64 vpMatrix3 warp_matrix));
65 extern void VPWarpA331N ANSI_ARGS((RGBIntPixel *in_image, int in_width,
66 int in_height, int in_bytes_per_scan, unsigned char *out_image,
67 int out_width, int out_height, int out_bytes_per_scan,
68 vpMatrix3 warp_matrix));
69 extern void VPWarpA330R ANSI_ARGS((RGBIntPixel *in_image, int in_width,
70 int in_height, int in_bytes_per_scan, unsigned char *out_image,
71 int out_width, int out_height, int out_bytes_per_scan,
72 vpMatrix3 warp_matrix));
73 extern void VPWarpA331R ANSI_ARGS((RGBIntPixel *in_image, int in_width,
74 int in_height, int in_bytes_per_scan, unsigned char *out_image,
75 int out_width, int out_height, int out_bytes_per_scan,
76 vpMatrix3 warp_matrix));
77
78 #ifdef STATISTICS
79 extern int vpResampleCount;
80 extern int vpCompositeCount;
81 extern int vpERTSkipCount;
82 extern int vpERTSkipAgainCount;
83 extern int vpERTUpdateCount;
84 extern int vpSpecialZeroSkipCount;
85 extern int vpRunFragmentCount;
86 #endif
87
88 /*
89 * VPRenderAffine
90 *
91 * Render a classified volume with an affine viewing transformation.
92 */
93
94 void
VPRenderAffine(vpc,algorithm,composite_func)95 VPRenderAffine(vpc, algorithm, composite_func)
96 vpContext *vpc;
97 int algorithm; /* USE_RLEVOLUME or USE_RAWVOLUME */
98 void (*composite_func)(); /* function to do the compositing */
99 {
100 int icount; /* voxels per voxel scanline */
101 int jcount; /* voxel scanlines per voxel slice */
102 int kcount; /* voxel slices in the volume */
103 int istride; /* strides for each dimension of raw volume */
104 int jstride;
105 int kstride;
106 int k; /* voxel slice index */
107 int kstart, kstop; /* values of k for first and last slices */
108 int kincr; /* value to add to k to get to the next slice
109 (either 1 or -1) */
110 RLEVoxels *rle_voxels; /* run-length encoded volume */
111 float slice_u, slice_v; /* sheared object space coordinates of the
112 top-left corner of the current constant-k
113 slice of the volume data */
114 int slice_u_int; /* integer part of slice_u and slice_v */
115 int slice_v_int;
116 float slice_u_frac; /* fractional part of slice_u and slice_v */
117 float slice_v_frac;
118 int slice_start_index; /* index of top-left int. image pixel */
119 float WgtTL, WgtBL, /* weights in the range 0..1 which give the */
120 WgtTR, WgtBR; /* fractional contribution of the */
121 /* neighboring voxels to the current */
122 /* intermediate image pixel */
123 int color_channels; /* number of color channels to compute */
124 float slice_depth_cueing; /* depth cueing factor for current slice */
125 float slice_dc_ratio; /* multiplier to get depth cueing factor
126 for the next slice */
127 unsigned char *run_lengths; /* run lengths for slice */
128 void *voxel_data; /* voxel data for slice */
129 void *intimage; /* first intermediate image pixel for slice */
130 int scan_offset_index; /* index into scan_offsets for this slice */
131 float shadow_slice_u; /* top-left corner of voxel slice in shadow */
132 float shadow_slice_v; /* buffer coordinates */
133 int shadow_slice_u_int; /* integer part of shadow_slice_u/v */
134 int shadow_slice_v_int;
135 int shadow_slice_start_index;/* index of top-left shadow buffer pixel */
136 GrayIntPixel *shadow_image; /* first shadow buffer pixel for slice */
137 int shadow_k; /* voxel slice number plus shadow bias */
138 #ifdef DUMP_SHADOW_VOLUME
139 unsigned char *shadow_dump;
140 #endif
141 #ifdef DUMP_GRAY_VOLUME
142 unsigned char *gray_dump;
143 #endif
144 #ifdef DUMP_SHADOW_VOLUME
145 int dump_fd;
146 int dump_value;
147 #else
148 #ifdef DUMP_GRAY_VOLUME
149 int dump_fd;
150 int dump_value;
151 #endif
152 #endif
153 #ifdef DEBUG
154 GrayIntPixel *trace_gray_ptr = &vpc->int_image.gray_intim[vpc->trace_u +
155 vpc->trace_v*vpc->intermediate_width];
156 RGBIntPixel *trace_rgb_ptr = &vpc->int_image.rgb_intim[vpc->trace_u +
157 vpc->trace_v*vpc->intermediate_width];
158 float vox_depth;
159 #endif
160 DECLARE_TIME(t0);
161 DECLARE_TIME(t1);
162 DECLARE_TIME(tA);
163 DECLARE_TIME(tB);
164
165 #ifdef STATISTICS
166 vpResampleCount = 0;
167 vpCompositeCount = 0;
168 vpERTSkipCount = 0;
169 vpERTSkipAgainCount = 0;
170 vpERTUpdateCount = 0;
171 vpSpecialZeroSkipCount = 0;
172 vpRunFragmentCount = 0;
173 #endif
174
175 GET_TIME(vpc, tA);
176
177 /* initialize for the fast classification algorithm */
178 if (algorithm == USE_RAWVOLUME && vpc->mm_octree != NULL) {
179 ASSERT(vpc->raw_voxels != NULL);
180 GET_TIME(vpc, t0);
181 VPComputeSummedAreaTable(vpc);
182 VPClassifyOctree(vpc);
183 GET_TIME(vpc, t1);
184 STORE_TIME(vpc, VPTIMER_CLSFY_OCTREE, t0, t1);
185 }
186
187 /* find size of volume */
188 if (algorithm == USE_RLEVOLUME) {
189 switch (vpc->best_view_axis) {
190 case VP_X_AXIS:
191 rle_voxels = vpc->rle_x;
192 break;
193 case VP_Y_AXIS:
194 rle_voxels = vpc->rle_y;
195 break;
196 case VP_Z_AXIS:
197 rle_voxels = vpc->rle_z;
198 break;
199 default:
200 VPBug("invalid viewing axis in AffineRender");
201 }
202 icount = rle_voxels->ilen;
203 jcount = rle_voxels->jlen;
204 kcount = rle_voxels->klen;
205 } else {
206 switch (vpc->best_view_axis) {
207 case VP_X_AXIS:
208 icount = vpc->ylen;
209 jcount = vpc->zlen;
210 kcount = vpc->xlen;
211 istride = vpc->ystride;
212 jstride = vpc->zstride;
213 kstride = vpc->xstride;
214 break;
215 case VP_Y_AXIS:
216 icount = vpc->zlen;
217 jcount = vpc->xlen;
218 kcount = vpc->ylen;
219 istride = vpc->zstride;
220 jstride = vpc->xstride;
221 kstride = vpc->ystride;
222 break;
223 case VP_Z_AXIS:
224 icount = vpc->xlen;
225 jcount = vpc->ylen;
226 kcount = vpc->zlen;
227 istride = vpc->xstride;
228 jstride = vpc->ystride;
229 kstride = vpc->zstride;
230 break;
231 default:
232 VPBug("invalid viewing axis in AffineRender");
233 }
234 }
235
236 GET_TIME(vpc, t0);
237
238 /* initialize intermediate image */
239 color_channels = vpc->color_channels;
240 vpc->pad_int_to_maxwidth = 0;
241 if (color_channels == 1) {
242 bzero(vpc->int_image.gray_intim, vpc->intermediate_width *
243 vpc->intermediate_height * sizeof(GrayIntPixel));
244 } else {
245 ASSERT(color_channels == 3);
246 bzero(vpc->int_image.rgb_intim, vpc->intermediate_width *
247 vpc->intermediate_height * sizeof(RGBIntPixel));
248 }
249
250 /* initialize shadow buffer */
251 if (vpc->enable_shadows) {
252 vpc->pad_shadow_to_maxwidth = 0;
253 bzero(vpc->shadow_buffer, vpc->shadow_width *
254 vpc->shadow_height * sizeof(GrayIntPixel));
255 }
256 #ifdef DUMP_SHADOW_VOLUME
257 Alloc(vpc, shadow_dump, char *, vpc->shadow_width * vpc->shadow_height *
258 kcount, "shadow_dump");
259 #endif
260 #ifdef DUMP_GRAY_VOLUME
261 Alloc(vpc, gray_dump, char *, vpc->intermediate_width *
262 vpc->intermediate_height * kcount, "gray_dump");
263 #endif
264
265 GET_TIME(vpc, t1);
266 STORE_TIME(vpc, VPTIMER_CLEAR, t0, t1);
267
268 /* initialize depth cueing */
269 slice_dc_ratio = VPSliceDepthCueRatio(vpc);
270 slice_depth_cueing = 1.;
271 #ifdef DEBUG
272 Debug((vpc, VPDEBUG_DEPTHCUE, "depth cueing at cube corners:\n"));
273 vox_depth = vpc->depth_000 + 0*vpc->depth_di +
274 0*vpc->depth_dj + 0*vpc->depth_dk;
275 if (vox_depth < 0.0)
276 vox_depth = 0.0;
277 Debug((vpc, VPDEBUG_DEPTHCUE,
278 " %3d %3d %3d: depth = %12.6f, factor = %12.6f\n",
279 0, 0, 0, vox_depth,
280 vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - vox_depth))));
281 vox_depth = vpc->depth_000 + icount*vpc->depth_di +
282 0*vpc->depth_dj + 0*vpc->depth_dk;
283 if (vox_depth < 0.0)
284 vox_depth = 0.0;
285 Debug((vpc, VPDEBUG_DEPTHCUE,
286 " %3d %3d %3d: depth = %12.6f, factor = %12.6f\n",
287 icount, 0, 0, vox_depth,
288 vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - vox_depth))));
289 vox_depth = vpc->depth_000 + icount*vpc->depth_di +
290 jcount*vpc->depth_dj + 0*vpc->depth_dk;
291 if (vox_depth < 0.0)
292 vox_depth = 0.0;
293 Debug((vpc, VPDEBUG_DEPTHCUE,
294 " %3d %3d %3d: depth = %12.6f, factor = %12.6f\n",
295 icount, jcount, 0, vox_depth,
296 vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - vox_depth))));
297 vox_depth = vpc->depth_000 + 0*vpc->depth_di +
298 jcount*vpc->depth_dj + 0*vpc->depth_dk;
299 if (vox_depth < 0.0)
300 vox_depth = 0.0;
301 Debug((vpc, VPDEBUG_DEPTHCUE,
302 " %3d %3d %3d: depth = %12.6f, factor = %12.6f\n",
303 0, jcount, 0, vox_depth,
304 vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - vox_depth))));
305 vox_depth = vpc->depth_000 + 0*vpc->depth_di +
306 0*vpc->depth_dj + kcount*vpc->depth_dk;
307 if (vox_depth < 0.0)
308 vox_depth = 0.0;
309 Debug((vpc, VPDEBUG_DEPTHCUE,
310 " %3d %3d %3d: depth = %12.6f, factor = %12.6f\n",
311 0, 0, kcount, vox_depth,
312 vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - vox_depth))));
313 vox_depth = vpc->depth_000 + icount*vpc->depth_di +
314 0*vpc->depth_dj + kcount*vpc->depth_dk;
315 if (vox_depth < 0.0)
316 vox_depth = 0.0;
317 Debug((vpc, VPDEBUG_DEPTHCUE,
318 " %3d %3d %3d: depth = %12.6f, factor = %12.6f\n",
319 icount, 0, kcount, vox_depth,
320 vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - vox_depth))));
321 vox_depth = vpc->depth_000 + icount*vpc->depth_di +
322 jcount*vpc->depth_dj + kcount*vpc->depth_dk;
323 if (vox_depth < 0.0)
324 vox_depth = 0.0;
325 Debug((vpc, VPDEBUG_DEPTHCUE,
326 " %3d %3d %3d: depth = %12.6f, factor = %12.6f\n",
327 icount, jcount, kcount, vox_depth,
328 vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - vox_depth))));
329 vox_depth = vpc->depth_000 + 0*vpc->depth_di +
330 jcount*vpc->depth_dj + kcount*vpc->depth_dk;
331 if (vox_depth < 0.0)
332 vox_depth = 0.0;
333 Debug((vpc, VPDEBUG_DEPTHCUE,
334 " %3d %3d %3d: depth = %12.6f, factor = %12.6f\n",
335 0, jcount, kcount, vox_depth,
336 vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - vox_depth))));
337 #endif /* DEBUG */
338
339 #ifdef DEBUG
340 /* initialize pixel tracing */
341 if (vpc->trace_u != -1) {
342 if (vpc->trace_u < 0 || vpc->trace_v < 0 ||
343 vpc->trace_u >= vpc->intermediate_width ||
344 vpc->trace_v >= vpc->intermediate_height) {
345 printf("Traced pixel is out of bounds.\n");
346 } else {
347 printf("Trace for pixel u=%d, v=%d",
348 vpc->trace_u, vpc->trace_v);
349 if (vpc->enable_shadows)
350 printf(", shadow_k=%d", vpc->trace_shadow_k);
351 printf(" (View %c, slice size %d,%d)\n",
352 vpc->best_view_axis + 'X', icount, jcount);
353 printf("Slice Slice TopLft BotLft ");
354 printf("TopRgt BotRgt Compos.\n");
355 printf(" BRX/BRY Opc/Clr/Wgt Opc/Clr/Wgt Opc/Clr/Wgt ");
356 printf("Opc/Clr/Wgt Opc/Clr\n");
357 }
358 }
359 #endif
360
361 /* compute outer loop bounds */
362 if (vpc->reverse_slice_order) {
363 kstart = kcount-1;
364 kstop = -1;
365 kincr = -1;
366 } else {
367 kstart = 0;
368 kincr = 1;
369 kstop = kcount;
370 }
371 shadow_k = kstart - vpc->shadow_bias * kincr;
372
373 /* loop over slices of the voxel data in front-to-back order */
374 for (k = kstart; k != kstop; k += kincr) {
375 ReportStatus(vpc, (double)(k - kstart)/(double)(kstop - kstart));
376
377 /* update shadow buffer */
378 if (vpc->enable_shadows && shadow_k >= 0 && shadow_k < kcount) {
379 /* compute coordinates of slice in shadow buffer;
380 shadow bias determines which slice (usually
381 a few slices old in order to eliminate self-shadowing) */
382 shadow_slice_u = vpc->shadow_shear_i * shadow_k +
383 vpc->shadow_trans_i;
384 shadow_slice_v = vpc->shadow_shear_j * shadow_k +
385 vpc->shadow_trans_j;
386 shadow_slice_u_int = (int)ceil(shadow_slice_u) - 1;
387 shadow_slice_v_int = (int)ceil(shadow_slice_v) - 1;
388 shadow_slice_start_index = shadow_slice_u_int +
389 shadow_slice_v_int*vpc->shadow_width;
390 shadow_image = &vpc->shadow_buffer[shadow_slice_start_index];
391
392 /* compute resampling weights for voxel slice in shadow buffer */
393 slice_u_frac = shadow_slice_u - shadow_slice_u_int;
394 slice_v_frac = shadow_slice_v - shadow_slice_v_int;
395 WgtTL = slice_u_frac * slice_v_frac;
396 WgtBL = slice_u_frac * ((float)1. - slice_v_frac);
397 WgtTR = ((float)1. - slice_u_frac) * slice_v_frac;
398 WgtBR = ((float)1. - slice_u_frac) * ((float)1. - slice_v_frac);
399
400 /* composite voxel opacities into shadow buffer */
401 if (algorithm == USE_RLEVOLUME) {
402 scan_offset_index = shadow_k *
403 rle_voxels->scan_offsets_per_slice;
404 run_lengths = rle_voxels->run_lengths +
405 rle_voxels->scan_offsets[scan_offset_index].first_len;
406 voxel_data = (void *)((char *)rle_voxels->data +
407 rle_voxels->scan_offsets[scan_offset_index].first_data);
408 VPCompAC00G(vpc, icount, jcount, shadow_k, slice_depth_cueing,
409 shadow_image, WgtTL, WgtBL, WgtTR, WgtBR,
410 run_lengths, voxel_data);
411 } else {
412 voxel_data = (void *)((char *)vpc->raw_voxels +
413 shadow_k*kstride);
414 VPCompAR00G(vpc, icount, jcount, shadow_k, slice_depth_cueing,
415 shadow_image, WgtTL, WgtBL, WgtTR, WgtBR,
416 voxel_data, istride, jstride);
417 }
418 }
419 shadow_k += kincr;
420
421 /* compute coordinates of top-left corner of voxel slice in
422 intermediate image */
423 slice_u = vpc->shear_i * k + vpc->trans_i;
424 slice_v = vpc->shear_j * k + vpc->trans_j;
425 slice_u_int = (int)ceil(slice_u) - 1;
426 slice_v_int = (int)ceil(slice_v) - 1;
427 slice_start_index = slice_u_int + slice_v_int*vpc->intermediate_width;
428 if (color_channels == 1)
429 intimage = &vpc->int_image.gray_intim[slice_start_index];
430 else
431 intimage = &vpc->int_image.rgb_intim[slice_start_index];
432
433 /* compute resampling weights for this slice */
434 slice_u_frac = slice_u - slice_u_int;
435 slice_v_frac = slice_v - slice_v_int;
436 WgtTL = slice_u_frac * slice_v_frac;
437 WgtBL = slice_u_frac * ((float)1. - slice_v_frac);
438 WgtTR = ((float)1. - slice_u_frac) * slice_v_frac;
439 WgtBR = ((float)1. - slice_u_frac) * ((float)1. - slice_v_frac);
440
441 /* compute coordinates of voxel slice in shadow buffer */
442 if (vpc->enable_shadows) {
443 shadow_slice_u = vpc->shadow_shear_i * k + vpc->shadow_trans_i;
444 shadow_slice_v = vpc->shadow_shear_j * k + vpc->shadow_trans_j;
445 shadow_slice_u_int = (int)ceil(shadow_slice_u) - 1;
446 shadow_slice_v_int = (int)ceil(shadow_slice_v) - 1;
447 shadow_slice_start_index = shadow_slice_u_int +
448 shadow_slice_v_int*vpc->shadow_width;
449 shadow_image = &vpc->shadow_buffer[shadow_slice_start_index];
450 }
451
452 /* find voxel data for this slice and composite */
453 if (algorithm == USE_RLEVOLUME) {
454 scan_offset_index = k * rle_voxels->scan_offsets_per_slice;
455 run_lengths = rle_voxels->run_lengths +
456 rle_voxels->scan_offsets[scan_offset_index].first_len;
457 voxel_data = (void *)((char *)rle_voxels->data +
458 rle_voxels->scan_offsets[scan_offset_index].first_data);
459 #ifdef INDEX_VOLUME
460 composite_func(vpc, icount, jcount, k, slice_depth_cueing,
461 intimage, WgtTL, WgtBL, WgtTR, WgtBR,
462 run_lengths, voxel_data,
463 rle_voxels->voxel_index + k * icount * jcount,
464 shadow_image);
465 #else
466 composite_func(vpc, icount, jcount, k, slice_depth_cueing,
467 intimage, WgtTL, WgtBL, WgtTR, WgtBR,
468 run_lengths, voxel_data, shadow_image);
469 #endif
470 } else {
471 voxel_data = (void *)((char *)vpc->raw_voxels + k*kstride);
472 composite_func(vpc, icount, jcount, k, slice_depth_cueing,
473 intimage, WgtTL, WgtBL, WgtTR, WgtBR,
474 voxel_data, istride, jstride, shadow_image);
475 }
476
477 /* update depth cueing factor */
478 slice_depth_cueing *= slice_dc_ratio;
479
480 #ifdef DUMP_SHADOW_VOLUME
481 vpGetImage(vpc, shadow_dump + k * vpc->shadow_width *
482 vpc->shadow_height, vpc->shadow_width, vpc->shadow_height,
483 vpc->shadow_width, VP_ALPHA, VP_SHADOW_BUFFER);
484 #endif
485 #ifdef DUMP_GRAY_VOLUME
486 vpGetImage(vpc, gray_dump + k * vpc->intermediate_width *
487 vpc->intermediate_height, vpc->intermediate_width,
488 vpc->intermediate_height, VP_LUMINANCE, VP_IMAGE_BUFFER);
489 #endif
490
491 }
492 ReportStatus(vpc, 1.0);
493
494 GET_TIME(vpc, t1);
495 STORE_TIME(vpc, VPTIMER_COMPOSITE, t0, t1);
496
497 #ifdef DEBUG
498 /* print traced pixel before depth cueing */
499 if (vpc->trace_u != -1) {
500 if (vpc->trace_u >= 0 && vpc->trace_v >= 0 &&
501 vpc->trace_u < vpc->intermediate_width &&
502 vpc->trace_v < vpc->intermediate_height) {
503 if (color_channels == 1) {
504 printf("Before depth cueing: opc = %.9f = %d",
505 trace_gray_ptr->opcflt*255.,
506 (int)(trace_gray_ptr->opcflt*255.));
507 printf(" clr = %.9f = %d\n",
508 trace_gray_ptr->clrflt,
509 (int)trace_gray_ptr->clrflt);
510 } else {
511 printf("Before depth cueing: opc = %14.9f = %3d",
512 trace_rgb_ptr->opcflt*255.,
513 (int)(trace_rgb_ptr->opcflt*255.));
514 printf(" r = %14.9f = %d\n",
515 trace_rgb_ptr->rclrflt,
516 (int)trace_rgb_ptr->rclrflt);
517 printf(" ");
518 printf(" g = %14.9f = %d\n",
519 trace_rgb_ptr->gclrflt,
520 (int)trace_rgb_ptr->gclrflt);
521 printf(" ");
522 printf(" b = %14.9f = %d\n",
523 trace_rgb_ptr->bclrflt,
524 (int)trace_rgb_ptr->bclrflt);
525 }
526 }
527 }
528 #endif
529
530 /* depth cue the intermediate image */
531 if (vpc->dc_enable) {
532 GET_TIME(vpc, t0);
533 VPDepthCueIntImage(vpc, vpc->reverse_slice_order ? kcount-1 : 0);
534 GET_TIME(vpc, t1);
535 STORE_TIME(vpc, VPTIMER_DEPTHCUE, t0, t1);
536 }
537
538 #ifdef DEBUG
539 /* print final value of traced pixel */
540 if (vpc->trace_u != -1) {
541 if (vpc->trace_u >= 0 && vpc->trace_v >= 0 &&
542 vpc->trace_u < vpc->intermediate_width &&
543 vpc->trace_v < vpc->intermediate_height) {
544 if (color_channels == 1) {
545 printf("Final pixel value: opc = %.9f = %d",
546 trace_gray_ptr->opcflt*255.,
547 (int)(trace_gray_ptr->opcflt*255.));
548 printf(" clr = %.9f = %d\n",
549 trace_gray_ptr->clrflt,
550 (int)trace_gray_ptr->clrflt);
551 } else {
552 printf("Final pixel value: opc = %14.9f = %3d",
553 trace_rgb_ptr->opcflt*255.,
554 (int)(trace_rgb_ptr->opcflt*255.));
555 printf(" r = %14.9f = %d\n",
556 trace_rgb_ptr->rclrflt,
557 (int)trace_rgb_ptr->rclrflt);
558 printf(" ");
559 printf(" g = %14.9f = %d\n",
560 trace_rgb_ptr->gclrflt,
561 (int)trace_rgb_ptr->gclrflt);
562 printf(" ");
563 printf(" b = %14.9f = %d\n",
564 trace_rgb_ptr->bclrflt,
565 (int)trace_rgb_ptr->bclrflt);
566 }
567 }
568 }
569 #endif
570
571 /* warp the intermediate image into the final image */
572 GET_TIME(vpc, t0);
573 switch (vpc->pixel_type) {
574 case VP_ALPHA:
575 if (color_channels == 1) {
576 VPWarpA101N(vpc->int_image.gray_intim, vpc->intermediate_width,
577 vpc->intermediate_height, sizeof(GrayIntPixel) *
578 (vpc->pad_int_to_maxwidth ?
579 vpc->max_intermediate_width:vpc->intermediate_width),
580 vpc->image, vpc->image_width, vpc->image_height,
581 vpc->image_bytes_per_scan, vpc->warp_2d);
582 } else {
583 VPWarpA301N(vpc->int_image.rgb_intim, vpc->intermediate_width,
584 vpc->intermediate_height, sizeof(RGBIntPixel) *
585 (vpc->pad_int_to_maxwidth ?
586 vpc->max_intermediate_width:vpc->intermediate_width),
587 vpc->image, vpc->image_width, vpc->image_height,
588 vpc->image_bytes_per_scan, vpc->warp_2d);
589 }
590 break;
591 case VP_LUMINANCE:
592 ASSERT(color_channels == 1);
593 VPWarpA110N(vpc->int_image.gray_intim, vpc->intermediate_width,
594 vpc->intermediate_height, sizeof(GrayIntPixel) *
595 (vpc->pad_int_to_maxwidth ?
596 vpc->max_intermediate_width:vpc->intermediate_width),
597 vpc->image, vpc->image_width, vpc->image_height,
598 vpc->image_bytes_per_scan, vpc->warp_2d);
599 break;
600 case VP_LUMINANCEA:
601 ASSERT(color_channels == 1);
602 VPWarpA111N(vpc->int_image.gray_intim, vpc->intermediate_width,
603 vpc->intermediate_height, sizeof(GrayIntPixel) *
604 (vpc->pad_int_to_maxwidth ?
605 vpc->max_intermediate_width:vpc->intermediate_width),
606 vpc->image, vpc->image_width, vpc->image_height,
607 vpc->image_bytes_per_scan, vpc->warp_2d);
608 break;
609 case VP_RGB:
610 ASSERT(color_channels == 3);
611 VPWarpA330N(vpc->int_image.rgb_intim, vpc->intermediate_width,
612 vpc->intermediate_height, sizeof(RGBIntPixel) *
613 (vpc->pad_int_to_maxwidth ?
614 vpc->max_intermediate_width:vpc->intermediate_width),
615 vpc->image, vpc->image_width, vpc->image_height,
616 vpc->image_bytes_per_scan, vpc->warp_2d);
617 break;
618 case VP_RGBA:
619 ASSERT(color_channels == 3);
620 VPWarpA331N(vpc->int_image.rgb_intim, vpc->intermediate_width,
621 vpc->intermediate_height, sizeof(RGBIntPixel) *
622 (vpc->pad_int_to_maxwidth ?
623 vpc->max_intermediate_width:vpc->intermediate_width),
624 vpc->image, vpc->image_width, vpc->image_height,
625 vpc->image_bytes_per_scan, vpc->warp_2d);
626 break;
627 case VP_BGR:
628 ASSERT(color_channels == 3);
629 VPWarpA330R(vpc->int_image.rgb_intim, vpc->intermediate_width,
630 vpc->intermediate_height, sizeof(RGBIntPixel) *
631 (vpc->pad_int_to_maxwidth ?
632 vpc->max_intermediate_width:vpc->intermediate_width),
633 vpc->image, vpc->image_width, vpc->image_height,
634 vpc->image_bytes_per_scan, vpc->warp_2d);
635 break;
636 case VP_ABGR:
637 ASSERT(color_channels == 3);
638 VPWarpA331R(vpc->int_image.rgb_intim, vpc->intermediate_width,
639 vpc->intermediate_height, sizeof(RGBIntPixel) *
640 (vpc->pad_int_to_maxwidth ?
641 vpc->max_intermediate_width:vpc->intermediate_width),
642 vpc->image, vpc->image_width, vpc->image_height,
643 vpc->image_bytes_per_scan, vpc->warp_2d);
644 break;
645 default:
646 VPBug("bad pixel type");
647 }
648 GET_TIME(vpc, t1);
649 STORE_TIME(vpc, VPTIMER_WARP, t0, t1);
650
651 GET_TIME(vpc, tB);
652 STORE_TIME(vpc, VPTIMER_RENDER, tA, tB);
653
654 #ifdef DUMP_SHADOW_VOLUME
655 printf("Dumping shadow map images to shadow.dump....");
656 fflush(stdout);
657 if ((dump_fd = creat("shadow.dump", 0644)) < 0)
658 VPBug("open failed");
659 dump_value = vpc->shadow_width;
660 write(dump_fd, &dump_value, sizeof(int));
661 dump_value = vpc->shadow_height;
662 write(dump_fd, &dump_value, sizeof(int));
663 dump_value = kcount;
664 write(dump_fd, &dump_value, sizeof(int));
665 write(dump_fd, shadow_dump, vpc->shadow_width * vpc->shadow_height *
666 kcount);
667 close(dump_fd);
668 printf("\n");
669 Dealloc(vpc, shadow_dump);
670 #endif
671
672 #ifdef DUMP_GRAY_VOLUME
673 printf("Dumping grayscale intermediate images to gray.dump....");
674 fflush(stdout);
675 if ((dump_fd = creat("gray.dump", 0644)) < 0)
676 VPBug("open failed");
677 dump_value = vpc->intermediate_width;
678 write(dump_fd, &dump_value, sizeof(int));
679 dump_value = vpc->intermediate_height;
680 write(dump_fd, &dump_value, sizeof(int));
681 dump_value = kcount;
682 write(dump_fd, &dump_value, sizeof(int));
683 write(dump_fd, gray_dump, vpc->intermediate_width *
684 vpc->intermediate_height * kcount);
685 close(dump_fd);
686 printf("\n");
687 Dealloc(vpc, gray_dump);
688 #endif
689 }
690
691 #ifdef DEBUG
692 /*
693 * vpPrintRayPath
694 *
695 * Print a trace of the voxels that contribute to the pixel specified
696 * with vpTracePixel.
697 */
698
699 vpResult
vpPrintRayPath(vpc)700 vpPrintRayPath(vpc)
701 vpContext *vpc;
702 {
703 int icount; /* voxels per voxel scanline */
704 int jcount; /* voxel scanlines per voxel slice */
705 int kcount; /* voxel slices in the volume */
706 int k; /* voxel slice index */
707 int kstart, kstop; /* values of k for first and last slices */
708 int kincr; /* value to add to k to get to the next slice
709 (either 1 or -1) */
710 float slice_u, slice_v; /* sheared object space coordinates of the
711 top-left corner of the current constant-k
712 slice of the volume data */
713 int slice_u_int; /* integer part of slice_u and slice_v */
714 int slice_v_int;
715 float slice_u_frac; /* fractional part of slice_u and slice_v */
716 float slice_v_frac;
717 float WgtTL, WgtBL, /* weights in the range 0..1 which give the */
718 WgtTR, WgtBR; /* fractional contribution of the */
719 /* neighboring voxels to the current */
720 /* intermediate image pixel */
721 int i, j; /* voxel coordinates in current slice of
722 the voxel to the BR of the ray */
723 int shadow_trace_u; /* coords. of shadow buffer pixel to trace */
724 int shadow_trace_v;
725 int retcode;
726
727 /* check for errors and initialize */
728 if ((retcode = VPFactorView(vpc)) != VP_OK)
729 return(retcode);
730 if (vpc->trace_u < 0 || vpc->trace_v < 0 ||
731 vpc->trace_u >= vpc->intermediate_width ||
732 vpc->trace_v >= vpc->intermediate_height) {
733 printf("Traced pixel is out of bounds.\n");
734 return(VP_OK);
735 }
736
737 /* find size of volume */
738 switch (vpc->best_view_axis) {
739 case VP_X_AXIS:
740 icount = vpc->ylen;
741 jcount = vpc->zlen;
742 kcount = vpc->xlen;
743 break;
744 case VP_Y_AXIS:
745 icount = vpc->zlen;
746 jcount = vpc->xlen;
747 kcount = vpc->ylen;
748 break;
749 case VP_Z_AXIS:
750 icount = vpc->xlen;
751 jcount = vpc->ylen;
752 kcount = vpc->zlen;
753 break;
754 default:
755 VPBug("invalid viewing axis in vpPrintRayPath");
756 }
757
758 /* print column headings */
759 printf("Ray path for pixel u=%d, v=%d", vpc->trace_u, vpc->trace_v);
760 if (vpc->enable_shadows)
761 printf(", shadow_k=%d", vpc->trace_shadow_k);
762 printf(" (View %c, slice size %d,%d)\n",
763 vpc->best_view_axis + 'X', icount, jcount);
764 printf("Slice TopLft BotLft TopRgt");
765 printf(" BotRgt\n");
766 printf(" _X_/_Y_/_Z_/Wgt _X_/_Y_/_Z_/Wgt _X_/_Y_/_Z_/Wgt");
767 printf(" _X_/_Y_/_Z_/Wgt\n");
768
769 /* compute outer loop bounds */
770 if (vpc->reverse_slice_order) {
771 kstart = kcount-1;
772 kstop = -1;
773 kincr = -1;
774 } else {
775 kstart = 0;
776 kincr = 1;
777 kstop = kcount;
778 }
779
780 /* loop over slices of the voxel data in front-to-back order */
781 for (k = kstart; k != kstop; k += kincr) {
782 /* compute coordinates of top-left corner of voxel slice in
783 intermediate image */
784 slice_u = vpc->shear_i * k + vpc->trans_i;
785 slice_v = vpc->shear_j * k + vpc->trans_j;
786 slice_u_int = (int)ceil(slice_u) - 1;
787 slice_v_int = (int)ceil(slice_v) - 1;
788
789 /* compute resampling weights for this slice */
790 slice_u_frac = slice_u - slice_u_int;
791 slice_v_frac = slice_v - slice_v_int;
792 WgtTL = slice_u_frac * slice_v_frac;
793 WgtBL = slice_u_frac * ((float)1. - slice_v_frac);
794 WgtTR = ((float)1. - slice_u_frac) * slice_v_frac;
795 WgtBR = ((float)1. - slice_u_frac) * ((float)1. - slice_v_frac);
796
797 /* compute intersection of the ray with this slice */
798 i = vpc->trace_u - slice_u_int;
799 j = vpc->trace_v - slice_v_int;
800
801 /* print ray location at this slice */
802 printf("[%3d]", k);
803 switch (vpc->best_view_axis) {
804 case VP_X_AXIS:
805 printf("%4d%4d%4d %3d ", k, i-1, j-1, (int)(WgtTL*100.));
806 printf("%4d%4d%4d %3d ", k, i-1, j, (int)(WgtBL*100.));
807 printf("%4d%4d%4d %3d ", k, i, j-1, (int)(WgtTR*100.));
808 printf("%4d%4d%4d %3d\n", k, i, j, (int)(WgtBR*100.));
809 break;
810 case VP_Y_AXIS:
811 printf("%4d%4d%4d %3d ", j-1, k, i-1, (int)(WgtTL*100.));
812 printf("%4d%4d%4d %3d ", j, k, i-1, (int)(WgtBL*100.));
813 printf("%4d%4d%4d %3d ", j-1, k, i, (int)(WgtTR*100.));
814 printf("%4d%4d%4d %3d\n", j, k, i, (int)(WgtBR*100.));
815 break;
816 case VP_Z_AXIS:
817 printf("%4d%4d%4d %3d ", i-1, j-1, k, (int)(WgtTL*100.));
818 printf("%4d%4d%4d %3d ", i-1, j, k, (int)(WgtBL*100.));
819 printf("%4d%4d%4d %3d ", i, j-1, k, (int)(WgtTR*100.));
820 printf("%4d%4d%4d %3d\n", i, j, k, (int)(WgtBR*100.));
821 break;
822 }
823 } /* for k */
824
825 if (!vpc->enable_shadows)
826 return(VP_OK);
827
828 /* compute coordinates of shadow buffer pixel to trace */
829 shadow_trace_u = vpc->trace_u +
830 (int)ceil(vpc->shadow_shear_i*vpc->trace_shadow_k+vpc->shadow_trans_i)-
831 (int)ceil(vpc->shear_i * vpc->trace_shadow_k + vpc->trans_i);
832 shadow_trace_v = vpc->trace_v +
833 (int)ceil(vpc->shadow_shear_j*vpc->trace_shadow_k+vpc->shadow_trans_j)-
834 (int)ceil(vpc->shear_j * vpc->trace_shadow_k + vpc->trans_j);
835
836 /* print column headings for shadow trace */
837 printf("\nShadow Ray Path (intersecting traced pixel at k=%d):\n",
838 vpc->trace_shadow_k);
839 printf("Slice TopLft BotLft TopRgt");
840 printf(" BotRgt\n");
841 printf(" _X_/_Y_/_Z_/Wgt _X_/_Y_/_Z_/Wgt _X_/_Y_/_Z_/Wgt");
842 printf(" _X_/_Y_/_Z_/Wgt\n");
843
844 /* loop over slices of the voxel data in front-to-back order */
845 for (k = kstart; k != kstop; k += kincr) {
846 /* compute coordinates of top-left corner of voxel slice in
847 intermediate image */
848 slice_u = vpc->shadow_shear_i * k + vpc->shadow_trans_i;
849 slice_v = vpc->shadow_shear_j * k + vpc->shadow_trans_j;
850 slice_u_int = (int)ceil(slice_u) - 1;
851 slice_v_int = (int)ceil(slice_v) - 1;
852
853 /* compute resampling weights for this slice */
854 slice_u_frac = slice_u - slice_u_int;
855 slice_v_frac = slice_v - slice_v_int;
856 WgtTL = slice_u_frac * slice_v_frac;
857 WgtBL = slice_u_frac * ((float)1. - slice_v_frac);
858 WgtTR = ((float)1. - slice_u_frac) * slice_v_frac;
859 WgtBR = ((float)1. - slice_u_frac) * ((float)1. - slice_v_frac);
860
861 /* compute intersection of the ray with this slice */
862 i = shadow_trace_u - slice_u_int;
863 j = shadow_trace_v - slice_v_int;
864
865 /* print ray location at this slice */
866 printf("[%3d]", k);
867 switch (vpc->best_view_axis) {
868 case VP_X_AXIS:
869 printf("%4d%4d%4d %3d ", k, i-1, j-1, (int)(WgtTL*100.));
870 printf("%4d%4d%4d %3d ", k, i-1, j, (int)(WgtBL*100.));
871 printf("%4d%4d%4d %3d ", k, i, j-1, (int)(WgtTR*100.));
872 printf("%4d%4d%4d %3d\n", k, i, j, (int)(WgtBR*100.));
873 break;
874 case VP_Y_AXIS:
875 printf("%4d%4d%4d %3d ", j-1, k, i-1, (int)(WgtTL*100.));
876 printf("%4d%4d%4d %3d ", j, k, i-1, (int)(WgtBL*100.));
877 printf("%4d%4d%4d %3d ", j-1, k, i, (int)(WgtTR*100.));
878 printf("%4d%4d%4d %3d\n", j, k, i, (int)(WgtBR*100.));
879 break;
880 case VP_Z_AXIS:
881 printf("%4d%4d%4d %3d ", i-1, j-1, k, (int)(WgtTL*100.));
882 printf("%4d%4d%4d %3d ", i-1, j, k, (int)(WgtBL*100.));
883 printf("%4d%4d%4d %3d ", i, j-1, k, (int)(WgtTR*100.));
884 printf("%4d%4d%4d %3d\n", i, j, k, (int)(WgtBR*100.));
885 break;
886 }
887 } /* for k */
888 return(VP_OK);
889 }
890 #endif /* DEBUG */
891