00001 /*
00002  * vp_renderB.c
00003  *
00004  * Brute-force shear-warp volume rendering algorithm.
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: 1999/02/18 00:26:22 $
00028  * $Revision: 1.1 $
00029  */
00030 
00031 #include "vp_global.h"
00032 
00033 static void AffineBruteForceRender ANSI_ARGS((vpContext *vpc));
00034 static void ClassifySlice ANSI_ARGS((vpContext *vpc, int slicenum,
00035     float *opc_slice));
00036 static void ShadeSlice ANSI_ARGS((vpContext *vpc, int slicenum,
00037     float *clr_slice));
00038 static void ScaleColors ANSI_ARGS((double scale, float *clr_slice,
00039     int width, int height, int color_channels));
00040 static void AlphaScaleColors ANSI_ARGS((float *opc_slice, float *clr_slice,
00041     int width, int height, int color_channels));
00042 static void DepthCueSlice ANSI_ARGS((vpContext *vpc, float *clr_slice,
00043     int width, int height, int color_channels, double depth_00k,
00044     double depth_di, double depth_dj));
00045 static void TranslateSlice ANSI_ARGS((float *opc_slice, float *clr_slice,
00046     int width, int height, double WgtTL_d, double WgtBL_d, double WgtTR_d,
00047     double WgtBR_d, int color_channels, float *resamp_opc_slice, 
00048     float *resmp_clr_slice));
00049 static void CompositeSlice ANSI_ARGS((float *resamp_opc, float *resamp_clr,
00050     int width, int height, int color_channels, void *int_image_ptr,
00051     int int_image_width, double min_opacity));
00052 static void AffineBruteForceWarp ANSI_ARGS((vpContext *vpc));
00053 
00054 #define RWCOX                  /* attempts to speed up a little */
00055 #ifdef RWCOX
00056 static float max_opacity ;
00057 #endif
00058 
00059 /*
00060  * vpBruteForceRender
00061  *
00062  * Render an unclassified volume using the basic shear-warp algorithm
00063  * without any optimizations (no spatial data structure is used and
00064  * coherence is ignored).  Use this routine as a standard for
00065  * correctness checking.
00066  */
00067 
00068 vpResult
00069 vpBruteForceRender(vpc)
00070 vpContext *vpc;
00071 {
00072     int retcode;
00073 
00074     /* check for errors and initialize */
00075     if ((retcode = VPCheckRawVolume(vpc)) != VP_OK)
00076         return(retcode);
00077     if ((retcode = VPCheckClassifier(vpc)) != VP_OK)
00078         return(retcode);
00079     if ((retcode = VPCheckShader(vpc)) != VP_OK)
00080         return(retcode);
00081     if ((retcode = VPCheckImage(vpc)) != VP_OK)
00082         return(retcode);
00083     if ((retcode = VPFactorView(vpc)) != VP_OK)
00084         return(retcode);
00085 
00086 #ifdef RWCOX
00087     max_opacity = vpc->max_opacity ;
00088 #endif
00089 
00090     /* render */
00091     if (vpc->affine_view)
00092         AffineBruteForceRender(vpc);
00093     else
00094         return(VPSetError(vpc, VPERROR_BAD_OPTION));
00095 
00096     return(VP_OK);
00097 }
00098 
00099 /*
00100  * AffineBruteForceRender
00101  *
00102  * Render an unclassified volume using the brute-force shear-warp
00103  * algorithm for an affine view transformation.
00104  */
00105 
00106 static void
00107 AffineBruteForceRender(vpc)
00108 vpContext *vpc;
00109 {
00110     int icount;                 /* voxels per voxel scanline */
00111     int jcount;                 /* voxel scanlines per voxel slice */
00112     int kcount;                 /* voxel slices in the volume */
00113     int k;                      /* voxel slice index */
00114     int kstart, kstop;          /* values of k for first and last slices */
00115     int kincr;                  /* value to add to k to get to the next slice
00116                                    (either 1 or -1) */
00117     float slice_u, slice_v;     /* sheared object space coordinates of the
00118                                    top-left corner of the current constant-k
00119                                    slice of the volume data */
00120     int slice_u_int;            /* integer part of slice_u and slice_v */
00121     int slice_v_int;
00122     float slice_u_frac;         /* fractional part of slice_u and slice_v */
00123     float slice_v_frac;
00124     int slice_start_index;      /* index of top-left int. image pixel */
00125     float WgtTL, WgtBL,         /* weights in the range 0..1 which give the */
00126           WgtTR, WgtBR;         /*   fractional contribution of the */
00127                                 /*   neighboring voxels to the current */
00128                                 /*   intermediate image pixel */
00129     int color_channels;         /* number of color channels to compute */
00130     float *opc_slice;           /* opacities after correction for viewpoint */
00131     float *resamp_opc_slice;    /* opacities after resampling */
00132     float *clr_slice;           /* colors for current voxel slice */
00133     float *resamp_clr_slice;    /* colors after resampling */
00134 #ifdef FAST_DEPTH_CUEING
00135     float slice_depth_cueing;   /* depth cueing factor for current slice */
00136     float slice_dc_ratio;       /* multiplier to get depth cueing factor
00137                                    for the next slice */
00138 #endif
00139     void *intim;                /* intermediate image pointer */
00140 
00141     Debug((vpc, VPDEBUG_RENDER, "Algorithm: affine brute force\n"));
00142 
00143     /* find size of volume */
00144     switch (vpc->best_view_axis) {
00145     case VP_X_AXIS:
00146         icount = vpc->ylen;
00147         jcount = vpc->zlen;
00148         kcount = vpc->xlen;
00149         break;
00150     case VP_Y_AXIS:
00151         icount = vpc->zlen;
00152         jcount = vpc->xlen;
00153         kcount = vpc->ylen;
00154         break;
00155     case VP_Z_AXIS:
00156         icount = vpc->xlen;
00157         jcount = vpc->ylen;
00158         kcount = vpc->zlen;
00159         break;
00160     default:
00161         VPBug("invalid viewing axis in AffineBruteForceRender");
00162     }
00163 
00164     /* initialize intermediate image */
00165     color_channels = vpc->color_channels;
00166     vpc->pad_int_to_maxwidth = 0;
00167     if (color_channels == 1) {
00168         bzero(vpc->int_image.gray_intim, vpc->intermediate_width *
00169               vpc->intermediate_height * sizeof(GrayIntPixel));
00170     } else {
00171         ASSERT(color_channels == 3);
00172         bzero(vpc->int_image.rgb_intim, vpc->intermediate_width *
00173               vpc->intermediate_height * sizeof(RGBIntPixel));
00174     }
00175 
00176     /* allocate memory for shaded and resampled voxel slices */
00177     Alloc(vpc, opc_slice, float *, icount*jcount*sizeof(float), "opc_slice");
00178     Alloc(vpc, resamp_opc_slice, float *, (icount+1)*(jcount+1)*sizeof(float),
00179           "resamp_opc_slice");
00180     Alloc(vpc, clr_slice, float *, color_channels*icount*jcount*sizeof(float),
00181           "clr_slice");
00182     Alloc(vpc, resamp_clr_slice, float *,
00183           color_channels*(icount+1)*(jcount+1)*sizeof(float),
00184           "resamp_clr_slice");
00185 
00186 #ifdef FAST_DEPTH_CUEING
00187     /* initialize depth cueing */
00188     if (vpc->dc_enable) {
00189         slice_dc_ratio = VPSliceDepthCueRatio(vpc);
00190         slice_depth_cueing = 1.;
00191     }
00192 #endif
00193 
00194     /* compute outer loop bounds */
00195     if (vpc->reverse_slice_order) {
00196         kstart = kcount-1;
00197         kstop = -1;
00198         kincr = -1;
00199     } else {
00200         kstart = 0;
00201         kincr = 1;
00202         kstop = kcount;
00203     }
00204 
00205     /* loop over slices of the voxel data in front-to-back order */
00206     for (k = kstart; k != kstop; k += kincr) {
00207         ReportStatus(vpc, (double)(k - kstart) / (double)(kstop - kstart));
00208 
00209         /* compute coordinates of top-left corner of slice in
00210            sheared object space */
00211         slice_u = vpc->shear_i * k + vpc->trans_i;
00212         slice_v = vpc->shear_j * k + vpc->trans_j;
00213         slice_u_int = (int)ceil(slice_u) - 1;
00214         slice_v_int = (int)ceil(slice_v) - 1;
00215 
00216         /* compute resampling weights for this slice */
00217         slice_u_frac = slice_u - slice_u_int;
00218         slice_v_frac = slice_v - slice_v_int;
00219         WgtTL = slice_u_frac * slice_v_frac;
00220         WgtBL = slice_u_frac * ((float)1. - slice_v_frac);
00221         WgtTR = ((float)1. - slice_u_frac) * slice_v_frac;
00222         WgtBR = ((float)1. - slice_u_frac) * ((float)1. - slice_v_frac);
00223 
00224         /* classify the slice of voxels */
00225         ClassifySlice(vpc, k, opc_slice);
00226 
00227         /* shade the slice of voxels */
00228         ShadeSlice(vpc, k, clr_slice);
00229 
00230         /* perform depth cueing on the slice */
00231         if (vpc->dc_enable) {
00232 #ifdef FAST_DEPTH_CUEING
00233             ScaleColors(slice_depth_cueing, clr_slice, icount, jcount,
00234                         color_channels);
00235             slice_depth_cueing *= slice_dc_ratio;
00236 #else
00237             DepthCueSlice(vpc, clr_slice, icount, jcount, color_channels,
00238                           vpc->depth_000 + k*vpc->depth_dk,
00239                           vpc->depth_di, vpc->depth_dj);
00240 #endif
00241         }
00242 
00243         /* weight the voxels colors by the voxel opacities */
00244         AlphaScaleColors(opc_slice, clr_slice, icount, jcount, color_channels);
00245 
00246         /* resample the slice of voxels */
00247         TranslateSlice(opc_slice, clr_slice, icount, jcount,
00248                        WgtTL, WgtBL, WgtTR, WgtBR, color_channels,
00249                        resamp_opc_slice, resamp_clr_slice);
00250 
00251         /* composite the slice of resampled voxels */
00252         slice_start_index = slice_u_int + slice_v_int*vpc->intermediate_width;
00253         if (color_channels == 1)
00254             intim = &vpc->int_image.gray_intim[slice_start_index];
00255         else
00256             intim = &vpc->int_image.rgb_intim[slice_start_index];
00257         CompositeSlice(resamp_opc_slice, resamp_clr_slice, icount+1, jcount+1,
00258                        color_channels, intim, vpc->intermediate_width,
00259                        vpc->min_opacity);
00260     }
00261     ReportStatus(vpc, 1.0);
00262 
00263 #ifdef FAST_DEPTH_CUEING
00264     /* depth cue the intermediate image */
00265     if (vpc->dc_enable)
00266         VPDepthCueIntImage(vpc, vpc->reverse_slice_order ? kcount-1 : 0);
00267 #endif
00268 
00269     /* warp the intermediate image into the final image */
00270     AffineBruteForceWarp(vpc);
00271 
00272     /* clean up */
00273     Dealloc(vpc, opc_slice);
00274     Dealloc(vpc, resamp_opc_slice);
00275     Dealloc(vpc, clr_slice);
00276     Dealloc(vpc, resamp_clr_slice);
00277 }
00278 
00279 /*
00280  * ClassifySlice
00281  *
00282  * Classify a slice of voxels.
00283  */
00284 
00285 static void
00286 ClassifySlice(vpc, slicenum, opc_slice)
00287 vpContext *vpc;
00288 int slicenum;
00289 float *opc_slice;
00290 {
00291     switch (vpc->best_view_axis) {
00292     case VP_X_AXIS:
00293         VPClassifyBlock(vpc, 1, slicenum, 0, 0, slicenum, vpc->ylen-1,
00294                         vpc->zlen-1, opc_slice, 0, sizeof(float),
00295                         vpc->ylen*sizeof(float));
00296         break;
00297     case VP_Y_AXIS:
00298         VPClassifyBlock(vpc, 1, 0, slicenum, 0, vpc->xlen-1, slicenum,
00299                         vpc->zlen-1, opc_slice, vpc->zlen*sizeof(float),
00300                         0, sizeof(float));
00301         break;
00302     case VP_Z_AXIS:
00303         VPClassifyBlock(vpc, 1, 0, 0, slicenum, vpc->xlen-1, vpc->ylen-1,
00304                         slicenum, opc_slice, sizeof(float),
00305                         vpc->xlen*sizeof(float), 0);
00306         break;
00307     }
00308 }
00309 
00310 /*
00311  * ShadeSlice
00312  *
00313  * Shade a slice of voxels.
00314  */
00315 
00316 static void
00317 ShadeSlice(vpc, slicenum, clr_slice)
00318 vpContext *vpc;
00319 int slicenum;
00320 float *clr_slice;
00321 {
00322     int color_bytes;
00323 
00324     color_bytes = sizeof(float) * vpc->color_channels;
00325 
00326     switch (vpc->best_view_axis) {
00327     case VP_X_AXIS:
00328         VPShadeBlock(vpc, slicenum, 0, 0, slicenum, vpc->ylen-1, vpc->zlen-1,
00329                      clr_slice, 0, color_bytes, vpc->ylen*color_bytes);
00330         break;
00331     case VP_Y_AXIS:
00332         VPShadeBlock(vpc, 0, slicenum, 0, vpc->xlen-1, slicenum, vpc->zlen-1,
00333                      clr_slice, vpc->zlen*color_bytes, 0, color_bytes);
00334         break;
00335     case VP_Z_AXIS:
00336         VPShadeBlock(vpc, 0, 0, slicenum, vpc->xlen-1, vpc->ylen-1, slicenum,
00337                      clr_slice, color_bytes, vpc->xlen*color_bytes, 0);
00338         break;
00339     }
00340 }
00341 
00342 /*
00343  * ScaleColors
00344  *
00345  * Weight voxel colors by a constant factor for the whole slice.
00346  */
00347 
00348 static void
00349 ScaleColors(scale, clr_slice, width, height, color_channels)
00350 double scale;
00351 float *clr_slice;
00352 int width;
00353 int height;
00354 int color_channels;
00355 {
00356     int i, j;
00357     float s;
00358 
00359     s = scale;
00360 
00361 #ifndef RWCOX
00362     for (j = 0; j < height; j++) {
00363         for (i = 0; i < width; i++) {
00364             if (color_channels == 1) {
00365                 clr_slice[0] *= s;
00366             } else {
00367                 clr_slice[0] *= s;
00368                 clr_slice[1] *= s;
00369                 clr_slice[2] *= s;
00370             }
00371             clr_slice += color_channels;
00372         }
00373     }
00374 #else
00375     if( color_channels == 1 ){
00376      for (j = 0; j < height; j++) {
00377         for (i = 0; i < width; i++) {
00378                 clr_slice[0] *= s;
00379             clr_slice ++;
00380         }
00381      }
00382     } else {
00383      for (j = 0; j < height; j++) {
00384         for (i = 0; i < width; i++) {
00385                 clr_slice[0] *= s;
00386                 clr_slice[1] *= s;
00387                 clr_slice[2] *= s;
00388             clr_slice += 3;
00389         }
00390      }
00391    }
00392 #endif
00393 }
00394 
00395 /*
00396  * AlphaScaleColors
00397  *
00398  * Weight voxel colors by voxels opacities.
00399  */
00400 
00401 static void
00402 AlphaScaleColors(opc_slice, clr_slice, width, height, color_channels)
00403 float *opc_slice;       /* 2D array of opacities (width by height) */
00404 float *clr_slice;       /* 2D array of colors (width by height) */
00405 int width;              /* size of voxel slice */
00406 int height;
00407 int color_channels;     /* number of color channels in clr_slice */
00408 {
00409     int i, j;
00410 
00411 #ifndef RWCOX
00412     for (j = 0; j < height; j++) {
00413         for (i = 0; i < width; i++) {
00414             if (color_channels == 1) {
00415                 clr_slice[0] *= opc_slice[0];
00416             } else {
00417                 clr_slice[0] *= opc_slice[0];
00418                 clr_slice[1] *= opc_slice[0];
00419                 clr_slice[2] *= opc_slice[0];
00420             }
00421             clr_slice += color_channels;
00422             opc_slice++;
00423         }
00424     }
00425 #else
00426     if( color_channels == 1 ){
00427      for (j = 0; j < height; j++) {
00428         for (i = 0; i < width; i++) {
00429                 clr_slice[0] *= opc_slice[0];
00430             clr_slice++;
00431             opc_slice++;
00432         }
00433      }
00434     } else {
00435      for (j = 0; j < height; j++) {
00436         for (i = 0; i < width; i++) {
00437                 clr_slice[0] *= opc_slice[0];
00438                 clr_slice[1] *= opc_slice[0];
00439                 clr_slice[2] *= opc_slice[0];
00440             clr_slice += 3;
00441             opc_slice++;
00442         }
00443      }
00444    }
00445 #endif
00446 }
00447 
00448 /*
00449  * DepthCueSlice
00450  *
00451  * Apply depth cueing factor to each voxel in a slice.
00452  */
00453 
00454 static void
00455 DepthCueSlice(vpc, clr_slice, width, height, color_channels,
00456               depth_00k, depth_di, depth_dj)
00457 vpContext *vpc;
00458 float *clr_slice;
00459 int width;
00460 int height;
00461 int color_channels;
00462 double depth_00k;               /* depth of top-left voxel in slice */
00463 double depth_di, depth_dj;      /* change in depth for unit change in
00464                                    i/j directions */
00465 {
00466     int i, j;
00467     double depth, depth_0jk, factor;
00468     double dc_front_factor, dc_density;
00469 
00470     dc_front_factor = vpc->dc_front_factor;
00471     dc_density = vpc->dc_density;
00472     depth_0jk = depth_00k;
00473     for (j = 0; j < height; j++) {
00474         depth = depth_0jk;
00475         for (i = 0; i < width; i++) {
00476             if (depth < 0.0)
00477                 factor = dc_front_factor * exp(-dc_density);
00478             else
00479                 factor = dc_front_factor * exp(-dc_density * (1.0 - depth));
00480             if (color_channels == 1) {
00481                 clr_slice[0] *= factor;
00482             } else {
00483                 clr_slice[0] *= factor;
00484                 clr_slice[1] *= factor;
00485                 clr_slice[2] *= factor;
00486             }
00487             clr_slice += color_channels;
00488             depth += depth_di;
00489         }
00490         depth_0jk += depth_dj;
00491     }
00492 }
00493 
00494 /*
00495  * TranslateSlice
00496  *
00497  * Translate and resample a slice of voxels.
00498  */
00499 
00500 static void
00501 TranslateSlice(opc_slice, clr_slice, width, height,
00502                WgtTL_d, WgtBL_d, WgtTR_d, WgtBR_d,
00503                color_channels, resamp_opc_slice, resamp_clr_slice)
00504 float *opc_slice;       /* 2D array of opacities (width by height) */
00505 float *clr_slice;       /* 2D array of colors (width by height) */
00506 int width;              /* size of voxel slice */
00507 int height;
00508 double WgtTL_d;         /* resampling weights */
00509 double WgtBL_d;
00510 double WgtTR_d;
00511 double WgtBR_d;
00512 int color_channels;     /* number of color channels in clr_slice */
00513 float *resamp_opc_slice;/* 2D array for storing resampled opacities
00514                            (width+1 by height+1) */
00515 float *resamp_clr_slice;/* 2D array for storing resampled colors
00516                            (width+1 by height+1) */
00517 {
00518     int i, j;
00519     float WgtTL, WgtBL, WgtTR, WgtBR;
00520     float OpcAcc, RClrAcc, GClrAcc, BClrAcc;
00521 
00522     WgtTL = WgtTL_d;
00523     WgtBL = WgtBL_d;
00524     WgtTR = WgtTR_d;
00525     WgtBR = WgtBR_d;
00526     for (j = 0; j <= height; j++) {
00527         for (i = 0; i <= width; i++) {
00528             OpcAcc = 0.;
00529             RClrAcc = 0.;
00530             GClrAcc = 0.;
00531             BClrAcc = 0.;
00532             if (i > 0 && j > 0) {
00533                 OpcAcc += WgtTL * opc_slice[-1-width];
00534                 if (color_channels == 1) {
00535                     RClrAcc += WgtTL * clr_slice[-1-width];
00536                 } else {
00537                     RClrAcc += WgtTL * clr_slice[3*(-1-width)];
00538                     GClrAcc += WgtTL * clr_slice[3*(-1-width)+1];
00539                     BClrAcc += WgtTL * clr_slice[3*(-1-width)+2];
00540                 }
00541             }
00542             if (i > 0 && j < height) {
00543                 OpcAcc += WgtBL * opc_slice[-1];
00544                 if (color_channels == 1) {
00545                     RClrAcc += WgtBL * clr_slice[-1];
00546                 } else {
00547                     RClrAcc += WgtBL * clr_slice[3*(-1)];
00548                     GClrAcc += WgtBL * clr_slice[3*(-1)+1];
00549                     BClrAcc += WgtBL * clr_slice[3*(-1)+2];
00550                 }
00551             }
00552             if (i < width && j > 0) {
00553                 OpcAcc += WgtTR * opc_slice[-width];
00554                 if (color_channels == 1) {
00555                     RClrAcc += WgtTR * clr_slice[-width];
00556                 } else {
00557                     RClrAcc += WgtTR * clr_slice[3*(-width)];
00558                     GClrAcc += WgtTR * clr_slice[3*(-width)+1];
00559                     BClrAcc += WgtTR * clr_slice[3*(-width)+2];
00560                 }
00561             }
00562             if (i < width && j < height) {
00563                 OpcAcc += WgtBR * opc_slice[0];
00564                 if (color_channels == 1) {
00565                     RClrAcc += WgtBR * clr_slice[0];
00566                 } else {
00567                     RClrAcc += WgtBR * clr_slice[3*(0)];
00568                     GClrAcc += WgtBR * clr_slice[3*(0)+1];
00569                     BClrAcc += WgtBR * clr_slice[3*(0)+2];
00570                 }
00571             }
00572             *resamp_opc_slice = OpcAcc;
00573             if (color_channels == 1) {
00574                 *resamp_clr_slice = RClrAcc;
00575             } else {
00576                 resamp_clr_slice[0] = RClrAcc;
00577                 resamp_clr_slice[1] = GClrAcc;
00578                 resamp_clr_slice[2] = BClrAcc;
00579             }
00580             resamp_opc_slice++;
00581             resamp_clr_slice += color_channels;
00582             if (i != width) {
00583                 opc_slice++;
00584                 clr_slice += color_channels;;
00585             }
00586         }
00587     }
00588 }
00589 
00590 /*
00591  * CompositeSlice
00592  *
00593  * Composite a resampled slice of voxels into the intermediate image.
00594  */
00595 
00596 static void
00597 CompositeSlice(resamp_opc, resamp_clr, width, height, color_channels,
00598                int_image_ptr, int_image_width, min_opacity)
00599 float *resamp_opc;      /* array of resampled opacities (width by height) */
00600 float *resamp_clr;      /* array of resampled colors (width by height) */
00601 int width;              /* size of resampled voxel arrays */
00602 int height;
00603 int color_channels;     /* number of color channels */
00604 void *int_image_ptr;    /* pointer to intermediate image pixel corresponding
00605                            to top-left resampled voxel */
00606 int int_image_width;    /* number of pixels in intermediate image scanline */
00607 double min_opacity;     /* low opacity threshold */
00608 {
00609     int i, j;
00610     float old_opc, old_r, old_g, old_b;
00611     float new_opc, new_r, new_g, new_b;
00612     GrayIntPixel *gray_intim;
00613     RGBIntPixel *rgb_intim;
00614 
00615     if (color_channels == 1)
00616         gray_intim = int_image_ptr;
00617     else
00618         rgb_intim = int_image_ptr;
00619     for (j = 0; j < height; j++) {
00620         for (i = 0; i < width; i++) {
00621             if (*resamp_opc > min_opacity) {
00622                 if (color_channels == 1) {
00623                     old_opc = gray_intim->opcflt;
00624 #ifdef RWCOX
00625                     if( old_opc < max_opacity ){
00626 #endif
00627                     old_r = gray_intim->clrflt;
00628                     new_opc = old_opc + *resamp_opc * ((float)1. - old_opc);
00629                     new_r = old_r + *resamp_clr * ((float)1. - old_opc);
00630                     gray_intim->opcflt = new_opc;
00631                     gray_intim->clrflt = new_r;
00632 #ifdef RWCOX
00633                     }
00634 #endif
00635                 } else {
00636                     old_opc = rgb_intim->opcflt;
00637 #ifdef RWCOX
00638                     if( old_opc < max_opacity ){
00639 #endif
00640                     old_r = rgb_intim->rclrflt;
00641                     old_g = rgb_intim->gclrflt;
00642                     old_b = rgb_intim->bclrflt;
00643                     new_opc = old_opc + *resamp_opc * ((float)1. - old_opc);
00644                     new_r = old_r + resamp_clr[0] * ((float)1. - old_opc);
00645                     new_g = old_g + resamp_clr[1] * ((float)1. - old_opc);
00646                     new_b = old_b + resamp_clr[2] * ((float)1. - old_opc);
00647                     rgb_intim->opcflt = new_opc;
00648                     rgb_intim->rclrflt = new_r;
00649                     rgb_intim->gclrflt = new_g;
00650                     rgb_intim->bclrflt = new_b;
00651 #ifdef RWCOX
00652                     }
00653 #endif
00654                 }
00655                 
00656             }
00657             resamp_opc++;
00658             if (color_channels == 1) {
00659                 resamp_clr++;
00660                 gray_intim++;
00661             } else {
00662                 resamp_clr += 3;
00663                 rgb_intim++;
00664             }
00665         } /* for i */
00666         if (color_channels == 1)
00667             gray_intim += int_image_width - width;
00668         else
00669             rgb_intim += int_image_width - width;
00670     } /* for j */
00671 }
00672 
00673 /*
00674  * AffineBruteForceWarp
00675  *
00676  * Warp the intermediate image into the final image (brute-force version,
00677  * affine transformations only).
00678  */
00679 
00680 static void
00681 AffineBruteForceWarp(vpc)
00682 vpContext *vpc;
00683 {
00684     unsigned char *int_image;   /* pointer to start of intermediate image
00685                                    (GrayIntPixel or RGBIntPixel) */
00686     int int_width;              /* size of intermediate image */
00687     int int_height;
00688     int int_scanbytes;          /* bytes per scanline in intermediate image */
00689     unsigned char *image;       /* final image pixel */
00690     int i, j;                   /* coordinates of final image pixel */
00691     float int_i_flt, int_j_flt; /* position of final image pixel in
00692                                    intermediate image coordinates */
00693     float int_i_int, int_j_int; /* truncated int_i_flt, int_j_flt */
00694     int int_i, int_j;           /* integer int_i_int, int_j_int */
00695     double alpha_i, alpha_j;    /* separable interpolation weights */
00696     double wgt;                 /* interpolation weight */
00697     GrayIntPixel *gray_pix;     /* intermediate image pixel (grayscale) */
00698     RGBIntPixel *rgb_pix;       /* intermediate image pixel (RGB) */
00699     double denom;
00700     double ma, mb, mc, md, me, mf;
00701     float r, g, b, alpha;
00702     int r_int, g_int, b_int, alpha_int;
00703     int color_channels;         /* number of color channels in int. image */
00704     int pixel_type;             /* type of output image pixel */
00705 
00706     /* initialize */
00707     color_channels = vpc->color_channels;
00708     pixel_type = vpc->pixel_type;
00709     int_width = vpc->intermediate_width;
00710     int_height = vpc->intermediate_height;
00711     if (vpc->color_channels == 1) {
00712         int_image = (unsigned char *)vpc->int_image.gray_intim;
00713         if (vpc->pad_int_to_maxwidth)
00714             int_scanbytes = vpc->max_intermediate_width*sizeof(GrayIntPixel);
00715         else
00716             int_scanbytes = vpc->intermediate_width*sizeof(GrayIntPixel);
00717     } else {
00718         int_image = (unsigned char *)vpc->int_image.rgb_intim;
00719         if (vpc->pad_int_to_maxwidth)
00720             int_scanbytes = vpc->max_intermediate_width*sizeof(RGBIntPixel);
00721         else
00722             int_scanbytes = vpc->intermediate_width*sizeof(RGBIntPixel);
00723     }
00724 
00725     /* compute transformation from final image pixel to intermediate
00726        image pixel */
00727     denom = 1. / (vpc->warp_2d[0][0]*vpc->warp_2d[1][1] -
00728                   vpc->warp_2d[0][1]*vpc->warp_2d[1][0]);
00729     ma = vpc->warp_2d[1][1] * denom;
00730     mb = -vpc->warp_2d[0][1] * denom;
00731     mc = (vpc->warp_2d[0][1]*vpc->warp_2d[1][2] - 
00732           vpc->warp_2d[1][1]*vpc->warp_2d[0][2]) * denom;
00733     md = -vpc->warp_2d[1][0] * denom;
00734     me = vpc->warp_2d[0][0] * denom;
00735     mf = (vpc->warp_2d[1][0]*vpc->warp_2d[0][2] - 
00736           vpc->warp_2d[0][0]*vpc->warp_2d[1][2]) * denom;
00737 
00738     /* loop over the pixels of the final image */
00739     for (j = 0; j < vpc->image_height; j++) {
00740         image = (unsigned char *)vpc->image + j*vpc->image_bytes_per_scan;
00741         for (i = 0; i < vpc->image_width; i++) {
00742             /* reverse-map final image pixel into intermediate image */
00743             int_i_flt = ma*i + mb*j + mc;
00744             int_j_flt = md*i + me*j + mf;
00745 
00746             /* compute interpolation weights */
00747             int_i_int = floor(int_i_flt);
00748             int_j_int = floor(int_j_flt);
00749             alpha_i = int_i_flt - int_i_int;
00750             alpha_j = int_j_flt - int_j_int;
00751             int_i = (int)int_i_int;
00752             int_j = (int)int_j_int;
00753 
00754             /* interpolate */
00755             r = 0;
00756             g = 0;
00757             b = 0;
00758             alpha = 0;
00759             if (int_i >= 0 && int_i < int_width &&
00760                 int_j >= 0 && int_j < int_height) {
00761                 wgt = (1. - alpha_i) * (1. - alpha_j);
00762                 if (color_channels == 1) {
00763                     gray_pix = (GrayIntPixel *)(int_image + int_j*
00764                                                 int_scanbytes) + int_i;
00765                     r += gray_pix->clrflt*wgt;
00766                     alpha += gray_pix->opcflt*wgt;
00767                 } else {
00768                     rgb_pix = (RGBIntPixel *)(int_image + int_j*
00769                                               int_scanbytes) + int_i;
00770                     r += rgb_pix->rclrflt*wgt;
00771                     g += rgb_pix->gclrflt*wgt;
00772                     b += rgb_pix->bclrflt*wgt;
00773                     alpha += rgb_pix->opcflt*wgt;
00774                 }
00775             }
00776             if (int_i >= 0 && int_i < int_width &&
00777                 int_j >= -1 && int_j < int_height-1) {
00778                 wgt = (1. - alpha_i) * alpha_j;
00779                 if (color_channels == 1) {
00780                     gray_pix = (GrayIntPixel *)(int_image + (int_j+1)*
00781                                                 int_scanbytes) + int_i;
00782                     r += gray_pix->clrflt*wgt;
00783                     alpha += gray_pix->opcflt*wgt;
00784                 } else {
00785                     rgb_pix = (RGBIntPixel *)(int_image + (int_j+1)*
00786                                               int_scanbytes) + int_i;
00787                     r += rgb_pix->rclrflt*wgt;
00788                     g += rgb_pix->gclrflt*wgt;
00789                     b += rgb_pix->bclrflt*wgt;
00790                     alpha += rgb_pix->opcflt*wgt;
00791                 }
00792             }
00793             if (int_i >= -1 && int_i < int_width-1 &&
00794                 int_j >= 0 && int_j < int_height) {
00795                 wgt = alpha_i * (1. - alpha_j);
00796                 if (color_channels == 1) {
00797                     gray_pix = (GrayIntPixel *)(int_image + int_j*
00798                                                 int_scanbytes) + int_i+1;
00799                     r += gray_pix->clrflt*wgt;
00800                     alpha += gray_pix->opcflt*wgt;
00801                 } else {
00802                     rgb_pix = (RGBIntPixel *)(int_image + int_j*
00803                                               int_scanbytes) + int_i+1;
00804                     r += rgb_pix->rclrflt*wgt;
00805                     g += rgb_pix->gclrflt*wgt;
00806                     b += rgb_pix->bclrflt*wgt;
00807                     alpha += rgb_pix->opcflt*wgt;
00808                 }
00809             }
00810             if (int_i >= -1 && int_i < int_width-1 &&
00811                 int_j >= -1 && int_j < int_height-1) {
00812                 wgt = alpha_i * alpha_j;
00813                 if (color_channels == 1) {
00814                     gray_pix = (GrayIntPixel *)(int_image + (int_j+1)*
00815                                                 int_scanbytes) + int_i+1;
00816                     r += gray_pix->clrflt*wgt;
00817                     alpha += gray_pix->opcflt*wgt;
00818                 } else {
00819                     rgb_pix = (RGBIntPixel *)(int_image + (int_j+1)*
00820                                               int_scanbytes) + int_i+1;
00821                     r += rgb_pix->rclrflt*wgt;
00822                     g += rgb_pix->gclrflt*wgt;
00823                     b += rgb_pix->bclrflt*wgt;
00824                     alpha += rgb_pix->opcflt*wgt;
00825                 }
00826             }
00827 
00828             /* clamp the pixel */
00829             if (alpha > 255.)
00830                 alpha_int = 255;
00831             else if (alpha < 0.)
00832                 alpha_int = 0;
00833             else
00834                 alpha_int = alpha;
00835             if (r > 255.)
00836                 r_int = 255;
00837             else if (r < 0)
00838                 r_int = 0;
00839             else
00840                 r_int = r;
00841             
00842             if (color_channels == 3) {
00843                 if (g > 255.)
00844                     g_int = 255;
00845                 else if (g < 0.)
00846                     g_int = 0;
00847                 else
00848                     g_int = g;
00849                 if (b > 255.)
00850                     b_int = 255;
00851                 else if (b < 0.)
00852                     b_int = 0;
00853                 else
00854                     b_int = b;
00855             }
00856 
00857             /* store the pixel */
00858             switch (pixel_type) {
00859             case VP_ALPHA:
00860                 *image++ = alpha_int;
00861                 break;
00862             case VP_LUMINANCE:
00863                 *image++ = r_int;
00864                 break;
00865             case VP_LUMINANCEA:
00866                 *image++ = r_int;
00867                 *image++ = alpha_int;
00868                 break;
00869             case VP_RGB:
00870                 *image++ = r_int;
00871                 *image++ = g_int;
00872                 *image++ = b_int;
00873                 break;
00874             case VP_RGBA:
00875                 *image++ = r_int;
00876                 *image++ = g_int;
00877                 *image++ = b_int;
00878                 *image++ = alpha_int;
00879                 break;
00880             case VP_BGR:
00881                 *image++ = b_int;
00882                 *image++ = g_int;
00883                 *image++ = r_int;
00884                 break;
00885             case VP_ABGR:
00886                 *image++ = alpha_int;
00887                 *image++ = b_int;
00888                 *image++ = g_int;
00889                 *image++ = r_int;
00890                 break;
00891             default:
00892                 VPBug("bad pixel type");
00893             }
00894         } /* for i */
00895     } /* for j */
00896 }
00897 
00898 #ifdef DEBUG
00899 StoreFloatImage(data, width, height, scale, filename)
00900 float *data;            /* array of input data */
00901 int width, height;      /* size of array */
00902 double scale;           /* factor for scaling pixel values */
00903 char *filename;         /* name of file to store result */
00904 {
00905     unsigned char *image, *imptr;
00906     int i, j;
00907 
00908     image = (unsigned char *)malloc(width*height);
00909     imptr = image;
00910     for (j = 0; j < height; j++) {
00911         for (i = 0; i < width; i++) {
00912             *imptr++ = (int)rint(scale * *data++);
00913         }
00914     }
00915     VprWriteGrayscaleTIFF(filename, width, height, width, image);
00916     free(image);
00917 }
00918 #endif