00001 /*
00002  * vp_renderA.c
00003  *
00004  * Shear-warp volume rendering algorithm for affine view transformations.
00005  *
00006  * Copyright (c) 1994 The Board of Trustees of The Leland Stanford
00007  * Junior University.  All rights reserved.
00008  *
00009  * Permission to use, copy, modify and distribute this software and its
00010  * documentation for any purpose is hereby granted without fee, provided
00011  * that the above copyright notice and this permission notice appear in
00012  * all copies of this software and that you do not sell the software.
00013  * Commercial licensing is available by contacting the author.
00014  * 
00015  * THE SOFTWARE IS PROVIDED "AS IS" AND WITHOUT WARRANTY OF ANY KIND,
00016  * EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
00017  * WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
00018  *
00019  * Author:
00020  *    Phil Lacroute
00021  *    Computer Systems Laboratory
00022  *    Electrical Engineering Dept.
00023  *    Stanford University
00024  */
00025 
00026 /*
00027  * $Date: 2001/12/17 16:16:23 $
00028  * $Revision: 1.1 $
00029  */
00030 
00031 #include "vp_global.h"
00032 
00033 /*#define DUMP_SHADOW_VOLUME*/
00034 /*#define DUMP_GRAY_VOLUME*/
00035 
00036 extern void VPCompAC00G ANSI_ARGS((vpContext *vpc, int icount, int jcount,
00037     int k, double slice_depth_cueing_dbl, GrayIntPixel *intimage,
00038     double weightTLdbl, double weightBLdbl, double weightTRdbl,
00039     double weightBRdbl, unsigned char *run_lengths, void *voxel_data));
00040 extern void VPCompAR00G ANSI_ARGS((vpContext *vpc, int icount, int jcount,
00041     int k, double slice_depth_cueing_dbl, GrayIntPixel *intimage,
00042     double weightTLdbl, double weightBLdbl, double weightTRdbl,
00043     double weightBRdbl, void *voxel_data, int voxel_istride,
00044     int voxel_jstride));
00045 extern void VPWarpA101N ANSI_ARGS((GrayIntPixel *in_image, int in_width,
00046     int in_height, int in_bytes_per_scan, unsigned char *out_image,
00047     int out_width, int out_height, int out_bytes_per_scan,
00048     vpMatrix3 warp_matrix));
00049 extern void VPWarpA110N ANSI_ARGS((GrayIntPixel *in_image, int in_width,
00050     int in_height, int in_bytes_per_scan, unsigned char *out_image,
00051     int out_width, int out_height, int out_bytes_per_scan,
00052     vpMatrix3 warp_matrix));
00053 extern void VPWarpA111N ANSI_ARGS((GrayIntPixel *in_image, int in_width,
00054     int in_height, int in_bytes_per_scan, unsigned char *out_image,
00055     int out_width, int out_height, int out_bytes_per_scan,
00056     vpMatrix3 warp_matrix));
00057 extern void VPWarpA301N ANSI_ARGS((RGBIntPixel *in_image, int in_width,
00058     int in_height, int in_bytes_per_scan, unsigned char *out_image,
00059     int out_width, int out_height, int out_bytes_per_scan,
00060     vpMatrix3 warp_matrix));
00061 extern void VPWarpA330N ANSI_ARGS((RGBIntPixel *in_image, int in_width,
00062     int in_height, int in_bytes_per_scan, unsigned char *out_image,
00063     int out_width, int out_height, int out_bytes_per_scan,
00064     vpMatrix3 warp_matrix));
00065 extern void VPWarpA331N ANSI_ARGS((RGBIntPixel *in_image, int in_width,
00066     int in_height, int in_bytes_per_scan, unsigned char *out_image,
00067     int out_width, int out_height, int out_bytes_per_scan,
00068     vpMatrix3 warp_matrix));
00069 extern void VPWarpA330R ANSI_ARGS((RGBIntPixel *in_image, int in_width,
00070     int in_height, int in_bytes_per_scan, unsigned char *out_image,
00071     int out_width, int out_height, int out_bytes_per_scan,
00072     vpMatrix3 warp_matrix));
00073 extern void VPWarpA331R ANSI_ARGS((RGBIntPixel *in_image, int in_width,
00074     int in_height, int in_bytes_per_scan, unsigned char *out_image,
00075     int out_width, int out_height, int out_bytes_per_scan,
00076     vpMatrix3 warp_matrix));
00077 
00078 #ifdef STATISTICS
00079 extern int vpResampleCount;
00080 extern int vpCompositeCount;
00081 extern int vpERTSkipCount;
00082 extern int vpERTSkipAgainCount;
00083 extern int vpERTUpdateCount;
00084 extern int vpSpecialZeroSkipCount;
00085 extern int vpRunFragmentCount;
00086 #endif
00087 
00088 /*
00089  * VPRenderAffine
00090  *
00091  * Render a classified volume with an affine viewing transformation.
00092  */
00093 
00094 void
00095 VPRenderAffine(vpc, algorithm, composite_func)
00096 vpContext *vpc;
00097 int algorithm;  /* USE_RLEVOLUME or USE_RAWVOLUME */
00098 void (*composite_func)(); /* function to do the compositing */
00099 {
00100     int icount;                 /* voxels per voxel scanline */
00101     int jcount;                 /* voxel scanlines per voxel slice */
00102     int kcount;                 /* voxel slices in the volume */
00103     int istride;                /* strides for each dimension of raw volume */
00104     int jstride;
00105     int kstride;
00106     int k;                      /* voxel slice index */
00107     int kstart, kstop;          /* values of k for first and last slices */
00108     int kincr;                  /* value to add to k to get to the next slice
00109                                    (either 1 or -1) */
00110     RLEVoxels *rle_voxels;      /* run-length encoded volume */
00111     float slice_u, slice_v;     /* sheared object space coordinates of the
00112                                    top-left corner of the current constant-k
00113                                    slice of the volume data */
00114     int slice_u_int;            /* integer part of slice_u and slice_v */
00115     int slice_v_int;
00116     float slice_u_frac;         /* fractional part of slice_u and slice_v */
00117     float slice_v_frac;
00118     int slice_start_index;      /* index of top-left int. image pixel */
00119     float WgtTL, WgtBL,         /* weights in the range 0..1 which give the */
00120           WgtTR, WgtBR;         /*   fractional contribution of the */
00121                                 /*   neighboring voxels to the current */
00122                                 /*   intermediate image pixel */
00123     int color_channels;         /* number of color channels to compute */
00124     float slice_depth_cueing;   /* depth cueing factor for current slice */
00125     float slice_dc_ratio;       /* multiplier to get depth cueing factor
00126                                    for the next slice */
00127     unsigned char *run_lengths; /* run lengths for slice */
00128     void *voxel_data;           /* voxel data for slice */
00129     void *intimage;             /* first intermediate image pixel for slice */
00130     int scan_offset_index;      /* index into scan_offsets for this slice */
00131     float shadow_slice_u;       /* top-left corner of voxel slice in shadow */
00132     float shadow_slice_v;       /*    buffer coordinates */
00133     int shadow_slice_u_int;     /* integer part of shadow_slice_u/v */
00134     int shadow_slice_v_int;
00135     int shadow_slice_start_index;/* index of top-left shadow buffer pixel */
00136     GrayIntPixel *shadow_image; /* first shadow buffer pixel for slice */
00137     int shadow_k;               /* voxel slice number plus shadow bias */
00138 #ifdef DUMP_SHADOW_VOLUME
00139     unsigned char *shadow_dump;
00140 #endif
00141 #ifdef DUMP_GRAY_VOLUME
00142     unsigned char *gray_dump;
00143 #endif
00144 #ifdef DUMP_SHADOW_VOLUME
00145     int dump_fd;
00146     int dump_value;
00147 #else
00148 #ifdef DUMP_GRAY_VOLUME
00149     int dump_fd;
00150     int dump_value;
00151 #endif
00152 #endif
00153 #ifdef DEBUG
00154     GrayIntPixel *trace_gray_ptr = &vpc->int_image.gray_intim[vpc->trace_u + 
00155                                     vpc->trace_v*vpc->intermediate_width];
00156     RGBIntPixel *trace_rgb_ptr = &vpc->int_image.rgb_intim[vpc->trace_u + 
00157                                     vpc->trace_v*vpc->intermediate_width];
00158     float vox_depth;
00159 #endif
00160     DECLARE_TIME(t0);
00161     DECLARE_TIME(t1);
00162     DECLARE_TIME(tA);
00163     DECLARE_TIME(tB);
00164 
00165 #ifdef STATISTICS
00166     vpResampleCount = 0;
00167     vpCompositeCount = 0;
00168     vpERTSkipCount = 0;
00169     vpERTSkipAgainCount = 0;
00170     vpERTUpdateCount = 0;
00171     vpSpecialZeroSkipCount = 0;
00172     vpRunFragmentCount = 0;
00173 #endif
00174 
00175     GET_TIME(vpc, tA);
00176 
00177     /* initialize for the fast classification algorithm */
00178     if (algorithm == USE_RAWVOLUME && vpc->mm_octree != NULL) {
00179         ASSERT(vpc->raw_voxels != NULL);
00180         GET_TIME(vpc, t0);
00181         VPComputeSummedAreaTable(vpc);
00182         VPClassifyOctree(vpc);
00183         GET_TIME(vpc, t1);
00184         STORE_TIME(vpc, VPTIMER_CLSFY_OCTREE, t0, t1);
00185     }
00186 
00187     /* find size of volume */
00188     if (algorithm == USE_RLEVOLUME) {
00189         switch (vpc->best_view_axis) {
00190         case VP_X_AXIS:
00191             rle_voxels = vpc->rle_x;
00192             break;
00193         case VP_Y_AXIS:
00194             rle_voxels = vpc->rle_y;
00195             break;
00196         case VP_Z_AXIS:
00197             rle_voxels = vpc->rle_z;
00198             break;
00199         default:
00200             VPBug("invalid viewing axis in AffineRender");
00201         }
00202         icount = rle_voxels->ilen;
00203         jcount = rle_voxels->jlen;
00204         kcount = rle_voxels->klen;
00205     } else {
00206         switch (vpc->best_view_axis) {
00207         case VP_X_AXIS:
00208             icount = vpc->ylen;
00209             jcount = vpc->zlen;
00210             kcount = vpc->xlen;
00211             istride = vpc->ystride;
00212             jstride = vpc->zstride;
00213             kstride = vpc->xstride;
00214             break;
00215         case VP_Y_AXIS:
00216             icount = vpc->zlen;
00217             jcount = vpc->xlen;
00218             kcount = vpc->ylen;
00219             istride = vpc->zstride;
00220             jstride = vpc->xstride;
00221             kstride = vpc->ystride;
00222             break;
00223         case VP_Z_AXIS:
00224             icount = vpc->xlen;
00225             jcount = vpc->ylen;
00226             kcount = vpc->zlen;
00227             istride = vpc->xstride;
00228             jstride = vpc->ystride;
00229             kstride = vpc->zstride;
00230             break;
00231         default:
00232             VPBug("invalid viewing axis in AffineRender");
00233         }
00234     }
00235 
00236     GET_TIME(vpc, t0);
00237 
00238     /* initialize intermediate image */
00239     color_channels = vpc->color_channels;
00240     vpc->pad_int_to_maxwidth = 0;
00241     if (color_channels == 1) {
00242         bzero(vpc->int_image.gray_intim, vpc->intermediate_width *
00243               vpc->intermediate_height * sizeof(GrayIntPixel));
00244     } else {
00245         ASSERT(color_channels == 3);
00246         bzero(vpc->int_image.rgb_intim, vpc->intermediate_width *
00247               vpc->intermediate_height * sizeof(RGBIntPixel));
00248     }
00249 
00250     /* initialize shadow buffer */
00251     if (vpc->enable_shadows) {
00252         vpc->pad_shadow_to_maxwidth = 0;
00253         bzero(vpc->shadow_buffer, vpc->shadow_width *
00254               vpc->shadow_height * sizeof(GrayIntPixel));
00255     }
00256 #ifdef DUMP_SHADOW_VOLUME
00257     Alloc(vpc, shadow_dump, char *, vpc->shadow_width * vpc->shadow_height *
00258           kcount, "shadow_dump");
00259 #endif
00260 #ifdef DUMP_GRAY_VOLUME
00261     Alloc(vpc, gray_dump, char *, vpc->intermediate_width * 
00262           vpc->intermediate_height * kcount, "gray_dump");
00263 #endif
00264 
00265     GET_TIME(vpc, t1);
00266     STORE_TIME(vpc, VPTIMER_CLEAR, t0, t1);
00267 
00268     /* initialize depth cueing */
00269     slice_dc_ratio = VPSliceDepthCueRatio(vpc);
00270     slice_depth_cueing = 1.;
00271 #ifdef DEBUG
00272     Debug((vpc, VPDEBUG_DEPTHCUE, "depth cueing at cube corners:\n"));
00273     vox_depth = vpc->depth_000 + 0*vpc->depth_di +
00274         0*vpc->depth_dj + 0*vpc->depth_dk;
00275     if (vox_depth < 0.0)
00276         vox_depth = 0.0;
00277     Debug((vpc, VPDEBUG_DEPTHCUE,
00278            "  %3d %3d %3d: depth = %12.6f, factor = %12.6f\n",
00279            0, 0, 0, vox_depth,
00280            vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - vox_depth))));
00281     vox_depth = vpc->depth_000 + icount*vpc->depth_di +
00282         0*vpc->depth_dj + 0*vpc->depth_dk;
00283     if (vox_depth < 0.0)
00284         vox_depth = 0.0;
00285     Debug((vpc, VPDEBUG_DEPTHCUE,
00286            "  %3d %3d %3d: depth = %12.6f, factor = %12.6f\n",
00287            icount, 0, 0, vox_depth,
00288            vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - vox_depth))));
00289     vox_depth = vpc->depth_000 + icount*vpc->depth_di +
00290         jcount*vpc->depth_dj + 0*vpc->depth_dk;
00291     if (vox_depth < 0.0)
00292         vox_depth = 0.0;
00293     Debug((vpc, VPDEBUG_DEPTHCUE,
00294            "  %3d %3d %3d: depth = %12.6f, factor = %12.6f\n",
00295            icount, jcount, 0, vox_depth,
00296            vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - vox_depth))));
00297     vox_depth = vpc->depth_000 + 0*vpc->depth_di +
00298         jcount*vpc->depth_dj + 0*vpc->depth_dk;
00299     if (vox_depth < 0.0)
00300         vox_depth = 0.0;
00301     Debug((vpc, VPDEBUG_DEPTHCUE,
00302            "  %3d %3d %3d: depth = %12.6f, factor = %12.6f\n",
00303            0, jcount, 0, vox_depth,
00304            vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - vox_depth))));
00305     vox_depth = vpc->depth_000 + 0*vpc->depth_di +
00306         0*vpc->depth_dj + kcount*vpc->depth_dk;
00307     if (vox_depth < 0.0)
00308         vox_depth = 0.0;
00309     Debug((vpc, VPDEBUG_DEPTHCUE,
00310            "  %3d %3d %3d: depth = %12.6f, factor = %12.6f\n",
00311            0, 0, kcount, vox_depth,
00312            vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - vox_depth))));
00313     vox_depth = vpc->depth_000 + icount*vpc->depth_di +
00314         0*vpc->depth_dj + kcount*vpc->depth_dk;
00315     if (vox_depth < 0.0)
00316         vox_depth = 0.0;
00317     Debug((vpc, VPDEBUG_DEPTHCUE,
00318            "  %3d %3d %3d: depth = %12.6f, factor = %12.6f\n",
00319            icount, 0, kcount, vox_depth,
00320            vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - vox_depth))));
00321     vox_depth = vpc->depth_000 + icount*vpc->depth_di +
00322         jcount*vpc->depth_dj + kcount*vpc->depth_dk;
00323     if (vox_depth < 0.0)
00324         vox_depth = 0.0;
00325     Debug((vpc, VPDEBUG_DEPTHCUE,
00326            "  %3d %3d %3d: depth = %12.6f, factor = %12.6f\n",
00327            icount, jcount, kcount, vox_depth,
00328            vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - vox_depth))));
00329     vox_depth = vpc->depth_000 + 0*vpc->depth_di +
00330         jcount*vpc->depth_dj + kcount*vpc->depth_dk;
00331     if (vox_depth < 0.0)
00332         vox_depth = 0.0;
00333     Debug((vpc, VPDEBUG_DEPTHCUE,
00334            "  %3d %3d %3d: depth = %12.6f, factor = %12.6f\n",
00335            0, jcount, kcount, vox_depth,
00336            vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - vox_depth))));
00337 #endif /* DEBUG */
00338 
00339 #ifdef DEBUG
00340     /* initialize pixel tracing */
00341     if (vpc->trace_u != -1) {
00342         if (vpc->trace_u < 0 || vpc->trace_v < 0 ||
00343             vpc->trace_u >= vpc->intermediate_width ||
00344             vpc->trace_v >= vpc->intermediate_height) {
00345             printf("Traced pixel is out of bounds.\n");
00346         } else {
00347             printf("Trace for pixel u=%d, v=%d",
00348                    vpc->trace_u, vpc->trace_v);
00349             if (vpc->enable_shadows)
00350                 printf(", shadow_k=%d", vpc->trace_shadow_k);
00351             printf(" (View %c, slice size %d,%d)\n",
00352                    vpc->best_view_axis + 'X', icount, jcount);
00353             printf("Slice   Slice      TopLft       BotLft       ");
00354             printf("TopRgt       BotRgt    Compos.\n");
00355             printf("       BRX/BRY  Opc/Clr/Wgt  Opc/Clr/Wgt  Opc/Clr/Wgt  ");
00356             printf("Opc/Clr/Wgt  Opc/Clr\n");
00357         }
00358     }
00359 #endif
00360 
00361     /* compute outer loop bounds */
00362     if (vpc->reverse_slice_order) {
00363         kstart = kcount-1;
00364         kstop = -1;
00365         kincr = -1;
00366     } else {
00367         kstart = 0;
00368         kincr = 1;
00369         kstop = kcount;
00370     }
00371     shadow_k = kstart - vpc->shadow_bias * kincr;
00372 
00373     /* loop over slices of the voxel data in front-to-back order */
00374     for (k = kstart; k != kstop; k += kincr) {
00375         ReportStatus(vpc, (double)(k - kstart)/(double)(kstop - kstart));
00376 
00377         /* update shadow buffer */
00378         if (vpc->enable_shadows && shadow_k >= 0 && shadow_k < kcount) {
00379             /* compute coordinates of slice in shadow buffer;
00380                shadow bias determines which slice (usually
00381                a few slices old in order to eliminate self-shadowing) */
00382             shadow_slice_u = vpc->shadow_shear_i * shadow_k + 
00383                 vpc->shadow_trans_i;
00384             shadow_slice_v = vpc->shadow_shear_j * shadow_k +
00385                 vpc->shadow_trans_j;
00386             shadow_slice_u_int = (int)ceil(shadow_slice_u) - 1;
00387             shadow_slice_v_int = (int)ceil(shadow_slice_v) - 1;
00388             shadow_slice_start_index = shadow_slice_u_int +
00389                 shadow_slice_v_int*vpc->shadow_width;
00390             shadow_image = &vpc->shadow_buffer[shadow_slice_start_index];
00391 
00392             /* compute resampling weights for voxel slice in shadow buffer */
00393             slice_u_frac = shadow_slice_u - shadow_slice_u_int;
00394             slice_v_frac = shadow_slice_v - shadow_slice_v_int;
00395             WgtTL = slice_u_frac * slice_v_frac;
00396             WgtBL = slice_u_frac * ((float)1. - slice_v_frac);
00397             WgtTR = ((float)1. - slice_u_frac) * slice_v_frac;
00398             WgtBR = ((float)1. - slice_u_frac) * ((float)1. - slice_v_frac);
00399 
00400             /* composite voxel opacities into shadow buffer */
00401             if (algorithm == USE_RLEVOLUME) {
00402                 scan_offset_index = shadow_k *
00403                     rle_voxels->scan_offsets_per_slice;
00404                 run_lengths = rle_voxels->run_lengths + 
00405                     rle_voxels->scan_offsets[scan_offset_index].first_len;
00406                 voxel_data = (void *)((char *)rle_voxels->data +
00407                     rle_voxels->scan_offsets[scan_offset_index].first_data);
00408                 VPCompAC00G(vpc, icount, jcount, shadow_k, slice_depth_cueing,
00409                             shadow_image, WgtTL, WgtBL, WgtTR, WgtBR,
00410                             run_lengths, voxel_data);
00411             } else {
00412                 voxel_data = (void *)((char *)vpc->raw_voxels + 
00413                                       shadow_k*kstride);
00414                 VPCompAR00G(vpc, icount, jcount, shadow_k, slice_depth_cueing,
00415                             shadow_image, WgtTL, WgtBL, WgtTR, WgtBR,
00416                             voxel_data, istride, jstride);
00417             }
00418         }
00419         shadow_k += kincr;
00420 
00421         /* compute coordinates of top-left corner of voxel slice in
00422            intermediate image */
00423         slice_u = vpc->shear_i * k + vpc->trans_i;
00424         slice_v = vpc->shear_j * k + vpc->trans_j;
00425         slice_u_int = (int)ceil(slice_u) - 1;
00426         slice_v_int = (int)ceil(slice_v) - 1;
00427         slice_start_index = slice_u_int + slice_v_int*vpc->intermediate_width;
00428         if (color_channels == 1)
00429             intimage = &vpc->int_image.gray_intim[slice_start_index];
00430         else
00431             intimage = &vpc->int_image.rgb_intim[slice_start_index];
00432 
00433         /* compute resampling weights for this slice */
00434         slice_u_frac = slice_u - slice_u_int;
00435         slice_v_frac = slice_v - slice_v_int;
00436         WgtTL = slice_u_frac * slice_v_frac;
00437         WgtBL = slice_u_frac * ((float)1. - slice_v_frac);
00438         WgtTR = ((float)1. - slice_u_frac) * slice_v_frac;
00439         WgtBR = ((float)1. - slice_u_frac) * ((float)1. - slice_v_frac);
00440 
00441         /* compute coordinates of voxel slice in shadow buffer */
00442         if (vpc->enable_shadows) {
00443             shadow_slice_u = vpc->shadow_shear_i * k + vpc->shadow_trans_i;
00444             shadow_slice_v = vpc->shadow_shear_j * k + vpc->shadow_trans_j;
00445             shadow_slice_u_int = (int)ceil(shadow_slice_u) - 1;
00446             shadow_slice_v_int = (int)ceil(shadow_slice_v) - 1;
00447             shadow_slice_start_index = shadow_slice_u_int +
00448                 shadow_slice_v_int*vpc->shadow_width;
00449             shadow_image = &vpc->shadow_buffer[shadow_slice_start_index];
00450         }
00451 
00452         /* find voxel data for this slice and composite */
00453         if (algorithm == USE_RLEVOLUME) {
00454             scan_offset_index = k * rle_voxels->scan_offsets_per_slice;
00455             run_lengths = rle_voxels->run_lengths + 
00456                       rle_voxels->scan_offsets[scan_offset_index].first_len;
00457             voxel_data = (void *)((char *)rle_voxels->data +
00458                      rle_voxels->scan_offsets[scan_offset_index].first_data);
00459 #ifdef INDEX_VOLUME
00460             composite_func(vpc, icount, jcount, k, slice_depth_cueing,
00461                            intimage, WgtTL, WgtBL, WgtTR, WgtBR,
00462                            run_lengths, voxel_data,
00463                            rle_voxels->voxel_index + k * icount * jcount,
00464                            shadow_image);
00465 #else
00466             composite_func(vpc, icount, jcount, k, slice_depth_cueing,
00467                            intimage, WgtTL, WgtBL, WgtTR, WgtBR,
00468                            run_lengths, voxel_data, shadow_image);
00469 #endif
00470         } else {
00471             voxel_data = (void *)((char *)vpc->raw_voxels + k*kstride);
00472             composite_func(vpc, icount, jcount, k, slice_depth_cueing,
00473                            intimage, WgtTL, WgtBL, WgtTR, WgtBR,
00474                            voxel_data, istride, jstride, shadow_image);
00475         }
00476 
00477         /* update depth cueing factor */
00478         slice_depth_cueing *= slice_dc_ratio;
00479 
00480 #ifdef DUMP_SHADOW_VOLUME
00481         vpGetImage(vpc, shadow_dump + k * vpc->shadow_width *
00482                    vpc->shadow_height, vpc->shadow_width, vpc->shadow_height,
00483                    vpc->shadow_width, VP_ALPHA, VP_SHADOW_BUFFER);
00484 #endif
00485 #ifdef DUMP_GRAY_VOLUME
00486         vpGetImage(vpc, gray_dump + k * vpc->intermediate_width *
00487                    vpc->intermediate_height, vpc->intermediate_width,
00488                    vpc->intermediate_height, VP_LUMINANCE, VP_IMAGE_BUFFER);
00489 #endif
00490 
00491     }
00492     ReportStatus(vpc, 1.0);
00493 
00494     GET_TIME(vpc, t1);
00495     STORE_TIME(vpc, VPTIMER_COMPOSITE, t0, t1);
00496 
00497 #ifdef DEBUG
00498     /* print traced pixel before depth cueing */
00499     if (vpc->trace_u != -1) {
00500         if (vpc->trace_u >= 0 && vpc->trace_v >= 0 &&
00501             vpc->trace_u < vpc->intermediate_width &&
00502             vpc->trace_v < vpc->intermediate_height) {
00503             if (color_channels == 1) {
00504                 printf("Before depth cueing: opc = %.9f = %d",
00505                        trace_gray_ptr->opcflt*255.,
00506                        (int)(trace_gray_ptr->opcflt*255.));
00507                 printf("   clr = %.9f = %d\n",
00508                        trace_gray_ptr->clrflt,
00509                        (int)trace_gray_ptr->clrflt);
00510             } else {
00511                 printf("Before depth cueing: opc = %14.9f = %3d",
00512                        trace_rgb_ptr->opcflt*255.,
00513                        (int)(trace_rgb_ptr->opcflt*255.));
00514                 printf("   r = %14.9f = %d\n",
00515                        trace_rgb_ptr->rclrflt,
00516                        (int)trace_rgb_ptr->rclrflt);
00517                 printf("                                               ");
00518                 printf("   g = %14.9f = %d\n",
00519                        trace_rgb_ptr->gclrflt,
00520                        (int)trace_rgb_ptr->gclrflt);
00521                 printf("                                               ");
00522                 printf("   b = %14.9f = %d\n",
00523                        trace_rgb_ptr->bclrflt,
00524                        (int)trace_rgb_ptr->bclrflt);
00525             }
00526         }
00527     }
00528 #endif
00529 
00530     /* depth cue the intermediate image */
00531     if (vpc->dc_enable) {
00532         GET_TIME(vpc, t0);
00533         VPDepthCueIntImage(vpc, vpc->reverse_slice_order ? kcount-1 : 0);
00534         GET_TIME(vpc, t1);
00535         STORE_TIME(vpc, VPTIMER_DEPTHCUE, t0, t1);
00536     }
00537 
00538 #ifdef DEBUG
00539     /* print final value of traced pixel */
00540     if (vpc->trace_u != -1) {
00541         if (vpc->trace_u >= 0 && vpc->trace_v >= 0 &&
00542             vpc->trace_u < vpc->intermediate_width &&
00543             vpc->trace_v < vpc->intermediate_height) {
00544             if (color_channels == 1) {
00545                 printf("Final pixel value:   opc = %.9f = %d",
00546                        trace_gray_ptr->opcflt*255.,
00547                        (int)(trace_gray_ptr->opcflt*255.));
00548                 printf("   clr = %.9f = %d\n",
00549                        trace_gray_ptr->clrflt,
00550                        (int)trace_gray_ptr->clrflt);
00551             } else {
00552                 printf("Final pixel value:   opc = %14.9f = %3d",
00553                        trace_rgb_ptr->opcflt*255.,
00554                        (int)(trace_rgb_ptr->opcflt*255.));
00555                 printf("   r = %14.9f = %d\n",
00556                        trace_rgb_ptr->rclrflt,
00557                        (int)trace_rgb_ptr->rclrflt);
00558                 printf("                                               ");
00559                 printf("   g = %14.9f = %d\n",
00560                        trace_rgb_ptr->gclrflt,
00561                        (int)trace_rgb_ptr->gclrflt);
00562                 printf("                                               ");
00563                 printf("   b = %14.9f = %d\n",
00564                        trace_rgb_ptr->bclrflt,
00565                        (int)trace_rgb_ptr->bclrflt);
00566             }
00567         }
00568     }
00569 #endif
00570 
00571     /* warp the intermediate image into the final image */
00572     GET_TIME(vpc, t0);
00573     switch (vpc->pixel_type) {
00574     case VP_ALPHA:
00575         if (color_channels == 1) {
00576             VPWarpA101N(vpc->int_image.gray_intim, vpc->intermediate_width,
00577                         vpc->intermediate_height, sizeof(GrayIntPixel) *
00578                         (vpc->pad_int_to_maxwidth ?
00579                          vpc->max_intermediate_width:vpc->intermediate_width),
00580                         vpc->image, vpc->image_width, vpc->image_height,
00581                         vpc->image_bytes_per_scan, vpc->warp_2d);
00582         } else {
00583             VPWarpA301N(vpc->int_image.rgb_intim, vpc->intermediate_width,
00584                         vpc->intermediate_height, sizeof(RGBIntPixel) *
00585                         (vpc->pad_int_to_maxwidth ?
00586                          vpc->max_intermediate_width:vpc->intermediate_width),
00587                         vpc->image, vpc->image_width, vpc->image_height,
00588                         vpc->image_bytes_per_scan, vpc->warp_2d);
00589         }
00590         break;
00591     case VP_LUMINANCE:
00592         ASSERT(color_channels == 1);
00593         VPWarpA110N(vpc->int_image.gray_intim, vpc->intermediate_width,
00594                     vpc->intermediate_height, sizeof(GrayIntPixel) *
00595                     (vpc->pad_int_to_maxwidth ?
00596                      vpc->max_intermediate_width:vpc->intermediate_width),
00597                     vpc->image, vpc->image_width, vpc->image_height,
00598                     vpc->image_bytes_per_scan, vpc->warp_2d);
00599         break;
00600     case VP_LUMINANCEA:
00601         ASSERT(color_channels == 1);
00602         VPWarpA111N(vpc->int_image.gray_intim, vpc->intermediate_width,
00603                     vpc->intermediate_height, sizeof(GrayIntPixel) *
00604                     (vpc->pad_int_to_maxwidth ?
00605                      vpc->max_intermediate_width:vpc->intermediate_width),
00606                     vpc->image, vpc->image_width, vpc->image_height,
00607                     vpc->image_bytes_per_scan, vpc->warp_2d);
00608         break;
00609     case VP_RGB:
00610         ASSERT(color_channels == 3);
00611         VPWarpA330N(vpc->int_image.rgb_intim, vpc->intermediate_width,
00612                     vpc->intermediate_height, sizeof(RGBIntPixel) *
00613                     (vpc->pad_int_to_maxwidth ?
00614                      vpc->max_intermediate_width:vpc->intermediate_width),
00615                     vpc->image, vpc->image_width, vpc->image_height,
00616                     vpc->image_bytes_per_scan, vpc->warp_2d);
00617         break;
00618     case VP_RGBA:
00619         ASSERT(color_channels == 3);
00620         VPWarpA331N(vpc->int_image.rgb_intim, vpc->intermediate_width,
00621                     vpc->intermediate_height, sizeof(RGBIntPixel) *
00622                     (vpc->pad_int_to_maxwidth ?
00623                      vpc->max_intermediate_width:vpc->intermediate_width),
00624                     vpc->image, vpc->image_width, vpc->image_height,
00625                     vpc->image_bytes_per_scan, vpc->warp_2d);
00626         break;
00627     case VP_BGR:
00628         ASSERT(color_channels == 3);
00629         VPWarpA330R(vpc->int_image.rgb_intim, vpc->intermediate_width,
00630                     vpc->intermediate_height, sizeof(RGBIntPixel) *
00631                     (vpc->pad_int_to_maxwidth ?
00632                      vpc->max_intermediate_width:vpc->intermediate_width),
00633                     vpc->image, vpc->image_width, vpc->image_height,
00634                     vpc->image_bytes_per_scan, vpc->warp_2d);
00635         break;
00636     case VP_ABGR:
00637         ASSERT(color_channels == 3);
00638         VPWarpA331R(vpc->int_image.rgb_intim, vpc->intermediate_width,
00639                     vpc->intermediate_height, sizeof(RGBIntPixel) *
00640                     (vpc->pad_int_to_maxwidth ?
00641                      vpc->max_intermediate_width:vpc->intermediate_width),
00642                     vpc->image, vpc->image_width, vpc->image_height,
00643                     vpc->image_bytes_per_scan, vpc->warp_2d);
00644         break;
00645     default:
00646         VPBug("bad pixel type");
00647     }
00648     GET_TIME(vpc, t1);
00649     STORE_TIME(vpc, VPTIMER_WARP, t0, t1);
00650 
00651     GET_TIME(vpc, tB);
00652     STORE_TIME(vpc, VPTIMER_RENDER, tA, tB);
00653 
00654 #ifdef DUMP_SHADOW_VOLUME
00655     printf("Dumping shadow map images to shadow.dump....");
00656     fflush(stdout);
00657     if ((dump_fd = creat("shadow.dump", 0644)) < 0)
00658         VPBug("open failed");
00659     dump_value = vpc->shadow_width;
00660     write(dump_fd, &dump_value, sizeof(int));
00661     dump_value = vpc->shadow_height;
00662     write(dump_fd, &dump_value, sizeof(int));
00663     dump_value = kcount;
00664     write(dump_fd, &dump_value, sizeof(int));
00665     write(dump_fd, shadow_dump, vpc->shadow_width * vpc->shadow_height *
00666           kcount);
00667     close(dump_fd);
00668     printf("\n");
00669     Dealloc(vpc, shadow_dump);
00670 #endif
00671 
00672 #ifdef DUMP_GRAY_VOLUME
00673     printf("Dumping grayscale intermediate images to gray.dump....");
00674     fflush(stdout);
00675     if ((dump_fd = creat("gray.dump", 0644)) < 0)
00676         VPBug("open failed");
00677     dump_value = vpc->intermediate_width;
00678     write(dump_fd, &dump_value, sizeof(int));
00679     dump_value = vpc->intermediate_height;
00680     write(dump_fd, &dump_value, sizeof(int));
00681     dump_value = kcount;
00682     write(dump_fd, &dump_value, sizeof(int));
00683     write(dump_fd, gray_dump, vpc->intermediate_width *
00684           vpc->intermediate_height * kcount);
00685     close(dump_fd);
00686     printf("\n");
00687     Dealloc(vpc, gray_dump);
00688 #endif
00689 }
00690 
00691 #ifdef DEBUG
00692 /*
00693  * vpPrintRayPath
00694  *
00695  * Print a trace of the voxels that contribute to the pixel specified
00696  * with vpTracePixel.
00697  */
00698 
00699 vpResult
00700 vpPrintRayPath(vpc)
00701 vpContext *vpc;
00702 {
00703     int icount;                 /* voxels per voxel scanline */
00704     int jcount;                 /* voxel scanlines per voxel slice */
00705     int kcount;                 /* voxel slices in the volume */
00706     int k;                      /* voxel slice index */
00707     int kstart, kstop;          /* values of k for first and last slices */
00708     int kincr;                  /* value to add to k to get to the next slice
00709                                    (either 1 or -1) */
00710     float slice_u, slice_v;     /* sheared object space coordinates of the
00711                                    top-left corner of the current constant-k
00712                                    slice of the volume data */
00713     int slice_u_int;            /* integer part of slice_u and slice_v */
00714     int slice_v_int;
00715     float slice_u_frac;         /* fractional part of slice_u and slice_v */
00716     float slice_v_frac;
00717     float WgtTL, WgtBL,         /* weights in the range 0..1 which give the */
00718           WgtTR, WgtBR;         /*   fractional contribution of the */
00719                                 /*   neighboring voxels to the current */
00720                                 /*   intermediate image pixel */
00721     int i, j;                   /* voxel coordinates in current slice of
00722                                    the voxel to the BR of the ray */
00723     int shadow_trace_u;         /* coords. of shadow buffer pixel to trace */
00724     int shadow_trace_v;
00725     int retcode;
00726 
00727     /* check for errors and initialize */
00728     if ((retcode = VPFactorView(vpc)) != VP_OK)
00729         return(retcode);
00730     if (vpc->trace_u < 0 || vpc->trace_v < 0 ||
00731         vpc->trace_u >= vpc->intermediate_width ||
00732         vpc->trace_v >= vpc->intermediate_height) {
00733         printf("Traced pixel is out of bounds.\n");
00734         return(VP_OK);
00735     }
00736 
00737     /* find size of volume */
00738     switch (vpc->best_view_axis) {
00739     case VP_X_AXIS:
00740         icount = vpc->ylen;
00741         jcount = vpc->zlen;
00742         kcount = vpc->xlen;
00743         break;
00744     case VP_Y_AXIS:
00745         icount = vpc->zlen;
00746         jcount = vpc->xlen;
00747         kcount = vpc->ylen;
00748         break;
00749     case VP_Z_AXIS:
00750         icount = vpc->xlen;
00751         jcount = vpc->ylen;
00752         kcount = vpc->zlen;
00753         break;
00754     default:
00755         VPBug("invalid viewing axis in vpPrintRayPath");
00756     }
00757 
00758     /* print column headings */
00759     printf("Ray path for pixel u=%d, v=%d", vpc->trace_u, vpc->trace_v);
00760     if (vpc->enable_shadows)
00761         printf(", shadow_k=%d", vpc->trace_shadow_k);
00762     printf(" (View %c, slice size %d,%d)\n",
00763            vpc->best_view_axis + 'X', icount, jcount);
00764     printf("Slice     TopLft            BotLft            TopRgt");
00765     printf("            BotRgt\n");
00766     printf("      _X_/_Y_/_Z_/Wgt   _X_/_Y_/_Z_/Wgt   _X_/_Y_/_Z_/Wgt");
00767     printf("   _X_/_Y_/_Z_/Wgt\n");
00768 
00769     /* compute outer loop bounds */
00770     if (vpc->reverse_slice_order) {
00771         kstart = kcount-1;
00772         kstop = -1;
00773         kincr = -1;
00774     } else {
00775         kstart = 0;
00776         kincr = 1;
00777         kstop = kcount;
00778     }
00779 
00780     /* loop over slices of the voxel data in front-to-back order */
00781     for (k = kstart; k != kstop; k += kincr) {
00782         /* compute coordinates of top-left corner of voxel slice in
00783            intermediate image */
00784         slice_u = vpc->shear_i * k + vpc->trans_i;
00785         slice_v = vpc->shear_j * k + vpc->trans_j;
00786         slice_u_int = (int)ceil(slice_u) - 1;
00787         slice_v_int = (int)ceil(slice_v) - 1;
00788 
00789         /* compute resampling weights for this slice */
00790         slice_u_frac = slice_u - slice_u_int;
00791         slice_v_frac = slice_v - slice_v_int;
00792         WgtTL = slice_u_frac * slice_v_frac;
00793         WgtBL = slice_u_frac * ((float)1. - slice_v_frac);
00794         WgtTR = ((float)1. - slice_u_frac) * slice_v_frac;
00795         WgtBR = ((float)1. - slice_u_frac) * ((float)1. - slice_v_frac);
00796 
00797         /* compute intersection of the ray with this slice */
00798         i = vpc->trace_u - slice_u_int;
00799         j = vpc->trace_v - slice_v_int;
00800 
00801         /* print ray location at this slice */
00802         printf("[%3d]", k);
00803         switch (vpc->best_view_axis) {
00804         case VP_X_AXIS:
00805             printf("%4d%4d%4d %3d  ", k, i-1, j-1, (int)(WgtTL*100.));
00806             printf("%4d%4d%4d %3d  ", k, i-1,   j, (int)(WgtBL*100.));
00807             printf("%4d%4d%4d %3d  ", k,   i, j-1, (int)(WgtTR*100.));
00808             printf("%4d%4d%4d %3d\n", k,   i,   j, (int)(WgtBR*100.));
00809             break;
00810         case VP_Y_AXIS:
00811             printf("%4d%4d%4d %3d  ", j-1, k, i-1, (int)(WgtTL*100.));
00812             printf("%4d%4d%4d %3d  ",   j, k, i-1, (int)(WgtBL*100.));
00813             printf("%4d%4d%4d %3d  ", j-1, k,   i, (int)(WgtTR*100.));
00814             printf("%4d%4d%4d %3d\n",   j, k,   i, (int)(WgtBR*100.));
00815             break;
00816         case VP_Z_AXIS:
00817             printf("%4d%4d%4d %3d  ", i-1, j-1, k, (int)(WgtTL*100.));
00818             printf("%4d%4d%4d %3d  ", i-1, j,   k, (int)(WgtBL*100.));
00819             printf("%4d%4d%4d %3d  ", i,   j-1, k, (int)(WgtTR*100.));
00820             printf("%4d%4d%4d %3d\n", i,   j,   k, (int)(WgtBR*100.));
00821             break;
00822         }
00823     } /* for k */
00824 
00825     if (!vpc->enable_shadows)
00826         return(VP_OK);
00827 
00828     /* compute coordinates of shadow buffer pixel to trace */
00829     shadow_trace_u = vpc->trace_u + 
00830         (int)ceil(vpc->shadow_shear_i*vpc->trace_shadow_k+vpc->shadow_trans_i)-
00831         (int)ceil(vpc->shear_i * vpc->trace_shadow_k + vpc->trans_i);
00832     shadow_trace_v = vpc->trace_v + 
00833         (int)ceil(vpc->shadow_shear_j*vpc->trace_shadow_k+vpc->shadow_trans_j)-
00834         (int)ceil(vpc->shear_j * vpc->trace_shadow_k + vpc->trans_j);
00835 
00836     /* print column headings for shadow trace */
00837     printf("\nShadow Ray Path (intersecting traced pixel at k=%d):\n",
00838            vpc->trace_shadow_k);
00839     printf("Slice     TopLft            BotLft            TopRgt");
00840     printf("            BotRgt\n");
00841     printf("      _X_/_Y_/_Z_/Wgt   _X_/_Y_/_Z_/Wgt   _X_/_Y_/_Z_/Wgt");
00842     printf("   _X_/_Y_/_Z_/Wgt\n");
00843 
00844     /* loop over slices of the voxel data in front-to-back order */
00845     for (k = kstart; k != kstop; k += kincr) {
00846         /* compute coordinates of top-left corner of voxel slice in
00847            intermediate image */
00848         slice_u = vpc->shadow_shear_i * k + vpc->shadow_trans_i;
00849         slice_v = vpc->shadow_shear_j * k + vpc->shadow_trans_j;
00850         slice_u_int = (int)ceil(slice_u) - 1;
00851         slice_v_int = (int)ceil(slice_v) - 1;
00852 
00853         /* compute resampling weights for this slice */
00854         slice_u_frac = slice_u - slice_u_int;
00855         slice_v_frac = slice_v - slice_v_int;
00856         WgtTL = slice_u_frac * slice_v_frac;
00857         WgtBL = slice_u_frac * ((float)1. - slice_v_frac);
00858         WgtTR = ((float)1. - slice_u_frac) * slice_v_frac;
00859         WgtBR = ((float)1. - slice_u_frac) * ((float)1. - slice_v_frac);
00860 
00861         /* compute intersection of the ray with this slice */
00862         i = shadow_trace_u - slice_u_int;
00863         j = shadow_trace_v - slice_v_int;
00864 
00865         /* print ray location at this slice */
00866         printf("[%3d]", k);
00867         switch (vpc->best_view_axis) {
00868         case VP_X_AXIS:
00869             printf("%4d%4d%4d %3d  ", k, i-1, j-1, (int)(WgtTL*100.));
00870             printf("%4d%4d%4d %3d  ", k, i-1,   j, (int)(WgtBL*100.));
00871             printf("%4d%4d%4d %3d  ", k,   i, j-1, (int)(WgtTR*100.));
00872             printf("%4d%4d%4d %3d\n", k,   i,   j, (int)(WgtBR*100.));
00873             break;
00874         case VP_Y_AXIS:
00875             printf("%4d%4d%4d %3d  ", j-1, k, i-1, (int)(WgtTL*100.));
00876             printf("%4d%4d%4d %3d  ",   j, k, i-1, (int)(WgtBL*100.));
00877             printf("%4d%4d%4d %3d  ", j-1, k,   i, (int)(WgtTR*100.));
00878             printf("%4d%4d%4d %3d\n",   j, k,   i, (int)(WgtBR*100.));
00879             break;
00880         case VP_Z_AXIS:
00881             printf("%4d%4d%4d %3d  ", i-1, j-1, k, (int)(WgtTL*100.));
00882             printf("%4d%4d%4d %3d  ", i-1, j,   k, (int)(WgtBL*100.));
00883             printf("%4d%4d%4d %3d  ", i,   j-1, k, (int)(WgtTR*100.));
00884             printf("%4d%4d%4d %3d\n", i,   j,   k, (int)(WgtBR*100.));
00885             break;
00886         }
00887     } /* for k */
00888     return(VP_OK);
00889 }
00890 #endif /* DEBUG */