/*
 * vp_renderB.c
 *
 * Brute-force shear-warp volume rendering algorithm.
 *
 * Copyright (c) 1994 The Board of Trustees of The Leland Stanford
 * Junior University.  All rights reserved.
 *
 * Permission to use, copy, modify and distribute this software and its
 * documentation for any purpose is hereby granted without fee, provided
 * that the above copyright notice and this permission notice appear in
 * all copies of this software and that you do not sell the software.
 * Commercial licensing is available by contacting the author.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS" AND WITHOUT WARRANTY OF ANY KIND,
 * EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
 * WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
 *
 * Author:
 *    Phil Lacroute
 *    Computer Systems Laboratory
 *    Electrical Engineering Dept.
 *    Stanford University
 */

/*
 * $Date: 1994/12/30 23:52:38 $
 * $Revision: 1.28 $
 */

#include "vp_global.h"

static void AffineBruteForceRender ANSI_ARGS((vpContext *vpc));
static void ClassifySlice ANSI_ARGS((vpContext *vpc, int slicenum,
    float *opc_slice));
static void ShadeSlice ANSI_ARGS((vpContext *vpc, int slicenum,
    float *clr_slice));
static void ScaleColors ANSI_ARGS((double scale, float *clr_slice,
    int width, int height, int color_channels));
static void AlphaScaleColors ANSI_ARGS((float *opc_slice, float *clr_slice,
    int width, int height, int color_channels));
static void DepthCueSlice ANSI_ARGS((vpContext *vpc, float *clr_slice,
    int width, int height, int color_channels, double depth_00k,
    double depth_di, double depth_dj));
static void TranslateSlice ANSI_ARGS((float *opc_slice, float *clr_slice,
    int width, int height, double WgtTL_d, double WgtBL_d, double WgtTR_d,
    double WgtBR_d, int color_channels, float *resamp_opc_slice, 
    float *resmp_clr_slice));
static void CompositeSlice ANSI_ARGS((float *resamp_opc, float *resamp_clr,
    int width, int height, int color_channels, void *int_image_ptr,
    int int_image_width, double min_opacity));
static void AffineBruteForceWarp ANSI_ARGS((vpContext *vpc));

/*
 * vpBruteForceRender
 *
 * Render an unclassified volume using the basic shear-warp algorithm
 * without any optimizations (no spatial data structure is used and
 * coherence is ignored).  Use this routine as a standard for
 * correctness checking.
 */

vpResult
vpBruteForceRender(vpc)
vpContext *vpc;
{
    int retcode;

    /* check for errors and initialize */
    if ((retcode = VPCheckRawVolume(vpc)) != VP_OK)
	return(retcode);
    if ((retcode = VPCheckClassifier(vpc)) != VP_OK)
	return(retcode);
    if ((retcode = VPCheckShader(vpc)) != VP_OK)
	return(retcode);
    if ((retcode = VPCheckImage(vpc)) != VP_OK)
	return(retcode);
    if ((retcode = VPFactorView(vpc)) != VP_OK)
	return(retcode);

    /* render */
    if (vpc->affine_view)
	AffineBruteForceRender(vpc);
    else
	return(VPSetError(vpc, VPERROR_BAD_OPTION));

    return(VP_OK);
}

/*
 * AffineBruteForceRender
 *
 * Render an unclassified volume using the brute-force shear-warp
 * algorithm for an affine view transformation.
 */

static void
AffineBruteForceRender(vpc)
vpContext *vpc;
{
    int icount;			/* voxels per voxel scanline */
    int jcount;			/* voxel scanlines per voxel slice */
    int kcount;			/* voxel slices in the volume */
    int k;			/* voxel slice index */
    int kstart, kstop;		/* values of k for first and last slices */
    int kincr;			/* value to add to k to get to the next slice
				   (either 1 or -1) */
    float slice_u, slice_v;	/* sheared object space coordinates of the
				   top-left corner of the current constant-k
				   slice of the volume data */
    int slice_u_int;		/* integer part of slice_u and slice_v */
    int slice_v_int;
    float slice_u_frac;		/* fractional part of slice_u and slice_v */
    float slice_v_frac;
    int slice_start_index;	/* index of top-left int. image pixel */
    float WgtTL, WgtBL,		/* weights in the range 0..1 which give the */
	  WgtTR, WgtBR;		/*   fractional contribution of the */
    				/*   neighboring voxels to the current */
    			        /*   intermediate image pixel */
    int color_channels;		/* number of color channels to compute */
    float *opc_slice;		/* opacities after correction for viewpoint */
    float *resamp_opc_slice;	/* opacities after resampling */
    float *clr_slice;		/* colors for current voxel slice */
    float *resamp_clr_slice;	/* colors after resampling */
#ifdef FAST_DEPTH_CUEING
    float slice_depth_cueing;	/* depth cueing factor for current slice */
    float slice_dc_ratio;	/* multiplier to get depth cueing factor
				   for the next slice */
#endif
    void *intim;		/* intermediate image pointer */

    Debug((vpc, VPDEBUG_RENDER, "Algorithm: affine brute force\n"));

    /* find size of volume */
    switch (vpc->best_view_axis) {
    case VP_X_AXIS:
	icount = vpc->ylen;
	jcount = vpc->zlen;
	kcount = vpc->xlen;
	break;
    case VP_Y_AXIS:
	icount = vpc->zlen;
	jcount = vpc->xlen;
	kcount = vpc->ylen;
	break;
    case VP_Z_AXIS:
	icount = vpc->xlen;
	jcount = vpc->ylen;
	kcount = vpc->zlen;
	break;
    default:
	VPBug("invalid viewing axis in AffineBruteForceRender");
    }

    /* initialize intermediate image */
    color_channels = vpc->color_channels;
    vpc->pad_int_to_maxwidth = 0;
    if (color_channels == 1) {
	bzero(vpc->int_image.gray_intim, vpc->intermediate_width *
	      vpc->intermediate_height * sizeof(GrayIntPixel));
    } else {
	ASSERT(color_channels == 3);
	bzero(vpc->int_image.rgb_intim, vpc->intermediate_width *
	      vpc->intermediate_height * sizeof(RGBIntPixel));
    }

    /* allocate memory for shaded and resampled voxel slices */
    Alloc(vpc, opc_slice, float *, icount*jcount*sizeof(float), "opc_slice");
    Alloc(vpc, resamp_opc_slice, float *, (icount+1)*(jcount+1)*sizeof(float),
	  "resamp_opc_slice");
    Alloc(vpc, clr_slice, float *, color_channels*icount*jcount*sizeof(float),
	  "clr_slice");
    Alloc(vpc, resamp_clr_slice, float *,
	  color_channels*(icount+1)*(jcount+1)*sizeof(float),
	  "resamp_clr_slice");

#ifdef FAST_DEPTH_CUEING
    /* initialize depth cueing */
    if (vpc->dc_enable) {
	slice_dc_ratio = VPSliceDepthCueRatio(vpc);
	slice_depth_cueing = 1.;
    }
#endif

    /* compute outer loop bounds */
    if (vpc->reverse_slice_order) {
	kstart = kcount-1;
	kstop = -1;
	kincr = -1;
    } else {
	kstart = 0;
	kincr = 1;
	kstop = kcount;
    }

    /* loop over slices of the voxel data in front-to-back order */
    for (k = kstart; k != kstop; k += kincr) {
	ReportStatus(vpc, (double)(k - kstart) / (double)(kstop - kstart));

	/* compute coordinates of top-left corner of slice in
	   sheared object space */
	slice_u = vpc->shear_i * k + vpc->trans_i;
	slice_v = vpc->shear_j * k + vpc->trans_j;
	slice_u_int = (int)ceil(slice_u) - 1;
	slice_v_int = (int)ceil(slice_v) - 1;

	/* compute resampling weights for this slice */
	slice_u_frac = slice_u - slice_u_int;
	slice_v_frac = slice_v - slice_v_int;
	WgtTL = slice_u_frac * slice_v_frac;
	WgtBL = slice_u_frac * ((float)1. - slice_v_frac);
	WgtTR = ((float)1. - slice_u_frac) * slice_v_frac;
	WgtBR = ((float)1. - slice_u_frac) * ((float)1. - slice_v_frac);

	/* classify the slice of voxels */
	ClassifySlice(vpc, k, opc_slice);

	/* shade the slice of voxels */
	ShadeSlice(vpc, k, clr_slice);

	/* perform depth cueing on the slice */
	if (vpc->dc_enable) {
#ifdef FAST_DEPTH_CUEING
	    ScaleColors(slice_depth_cueing, clr_slice, icount, jcount,
			color_channels);
	    slice_depth_cueing *= slice_dc_ratio;
#else
	    DepthCueSlice(vpc, clr_slice, icount, jcount, color_channels,
			  vpc->depth_000 + k*vpc->depth_dk,
			  vpc->depth_di, vpc->depth_dj);
#endif
	}

	/* weight the voxels colors by the voxel opacities */
	AlphaScaleColors(opc_slice, clr_slice, icount, jcount, color_channels);

	/* resample the slice of voxels */
	TranslateSlice(opc_slice, clr_slice, icount, jcount,
		       WgtTL, WgtBL, WgtTR, WgtBR, color_channels,
		       resamp_opc_slice, resamp_clr_slice);

	/* composite the slice of resampled voxels */
	slice_start_index = slice_u_int + slice_v_int*vpc->intermediate_width;
	if (color_channels == 1)
	    intim = &vpc->int_image.gray_intim[slice_start_index];
	else
	    intim = &vpc->int_image.rgb_intim[slice_start_index];
	CompositeSlice(resamp_opc_slice, resamp_clr_slice, icount+1, jcount+1,
		       color_channels, intim, vpc->intermediate_width,
		       vpc->min_opacity);
    }
    ReportStatus(vpc, 1.0);

#ifdef FAST_DEPTH_CUEING
    /* depth cue the intermediate image */
    if (vpc->dc_enable)
	VPDepthCueIntImage(vpc, vpc->reverse_slice_order ? kcount-1 : 0);
#endif

    /* warp the intermediate image into the final image */
    AffineBruteForceWarp(vpc);

    /* clean up */
    Dealloc(vpc, opc_slice);
    Dealloc(vpc, resamp_opc_slice);
    Dealloc(vpc, clr_slice);
    Dealloc(vpc, resamp_clr_slice);
}

/*
 * ClassifySlice
 *
 * Classify a slice of voxels.
 */

static void
ClassifySlice(vpc, slicenum, opc_slice)
vpContext *vpc;
int slicenum;
float *opc_slice;
{
    switch (vpc->best_view_axis) {
    case VP_X_AXIS:
	VPClassifyBlock(vpc, 1, slicenum, 0, 0, slicenum, vpc->ylen-1,
			vpc->zlen-1, opc_slice, 0, sizeof(float),
			vpc->ylen*sizeof(float));
	break;
    case VP_Y_AXIS:
	VPClassifyBlock(vpc, 1, 0, slicenum, 0, vpc->xlen-1, slicenum,
			vpc->zlen-1, opc_slice, vpc->zlen*sizeof(float),
			0, sizeof(float));
	break;
    case VP_Z_AXIS:
	VPClassifyBlock(vpc, 1, 0, 0, slicenum, vpc->xlen-1, vpc->ylen-1,
			slicenum, opc_slice, sizeof(float),
			vpc->xlen*sizeof(float), 0);
	break;
    }
}

/*
 * ShadeSlice
 *
 * Shade a slice of voxels.
 */

static void
ShadeSlice(vpc, slicenum, clr_slice)
vpContext *vpc;
int slicenum;
float *clr_slice;
{
    int color_bytes;

    color_bytes = sizeof(float) * vpc->color_channels;

    switch (vpc->best_view_axis) {
    case VP_X_AXIS:
	VPShadeBlock(vpc, slicenum, 0, 0, slicenum, vpc->ylen-1, vpc->zlen-1,
		     clr_slice, 0, color_bytes, vpc->ylen*color_bytes);
	break;
    case VP_Y_AXIS:
	VPShadeBlock(vpc, 0, slicenum, 0, vpc->xlen-1, slicenum, vpc->zlen-1,
		     clr_slice, vpc->zlen*color_bytes, 0, color_bytes);
	break;
    case VP_Z_AXIS:
	VPShadeBlock(vpc, 0, 0, slicenum, vpc->xlen-1, vpc->ylen-1, slicenum,
		     clr_slice, color_bytes, vpc->xlen*color_bytes, 0);
	break;
    }
}

/*
 * ScaleColors
 *
 * Weight voxel colors by a constant factor for the whole slice.
 */

static void
ScaleColors(scale, clr_slice, width, height, color_channels)
double scale;
float *clr_slice;
int width;
int height;
int color_channels;
{
    int i, j;
    float s;

    s = scale;
    for (j = 0; j < height; j++) {
	for (i = 0; i < width; i++) {
	    if (color_channels == 1) {
		clr_slice[0] *= s;
	    } else {
		clr_slice[0] *= s;
		clr_slice[1] *= s;
		clr_slice[2] *= s;
	    }
	    clr_slice += color_channels;
	}
    }
}

/*
 * AlphaScaleColors
 *
 * Weight voxel colors by voxels opacities.
 */

static void
AlphaScaleColors(opc_slice, clr_slice, width, height, color_channels)
float *opc_slice;	/* 2D array of opacities (width by height) */
float *clr_slice;	/* 2D array of colors (width by height) */
int width;		/* size of voxel slice */
int height;
int color_channels;	/* number of color channels in clr_slice */
{
    int i, j;

    for (j = 0; j < height; j++) {
	for (i = 0; i < width; i++) {
	    if (color_channels == 1) {
		clr_slice[0] *= opc_slice[0];
	    } else {
		clr_slice[0] *= opc_slice[0];
		clr_slice[1] *= opc_slice[0];
		clr_slice[2] *= opc_slice[0];
	    }
	    clr_slice += color_channels;
	    opc_slice++;
	}
    }
}

/*
 * DepthCueSlice
 *
 * Apply depth cueing factor to each voxel in a slice.
 */

static void
DepthCueSlice(vpc, clr_slice, width, height, color_channels,
	      depth_00k, depth_di, depth_dj)
vpContext *vpc;
float *clr_slice;
int width;
int height;
int color_channels;
double depth_00k;		/* depth of top-left voxel in slice */
double depth_di, depth_dj;	/* change in depth for unit change in
				   i/j directions */
{
    int i, j;
    double depth, depth_0jk, factor;
    double dc_front_factor, dc_density;

    dc_front_factor = vpc->dc_front_factor;
    dc_density = vpc->dc_density;
    depth_0jk = depth_00k;
    for (j = 0; j < height; j++) {
	depth = depth_0jk;
	for (i = 0; i < width; i++) {
	    if (depth < 0.0)
		factor = dc_front_factor * exp(-dc_density);
	    else
		factor = dc_front_factor * exp(-dc_density * (1.0 - depth));
	    if (color_channels == 1) {
		clr_slice[0] *= factor;
	    } else {
		clr_slice[0] *= factor;
		clr_slice[1] *= factor;
		clr_slice[2] *= factor;
	    }
	    clr_slice += color_channels;
	    depth += depth_di;
	}
	depth_0jk += depth_dj;
    }
}

/*
 * TranslateSlice
 *
 * Translate and resample a slice of voxels.
 */

static void
TranslateSlice(opc_slice, clr_slice, width, height,
	       WgtTL_d, WgtBL_d, WgtTR_d, WgtBR_d,
	       color_channels, resamp_opc_slice, resamp_clr_slice)
float *opc_slice;	/* 2D array of opacities (width by height) */
float *clr_slice;	/* 2D array of colors (width by height) */
int width;		/* size of voxel slice */
int height;
double WgtTL_d;		/* resampling weights */
double WgtBL_d;
double WgtTR_d;
double WgtBR_d;
int color_channels;	/* number of color channels in clr_slice */
float *resamp_opc_slice;/* 2D array for storing resampled opacities
			   (width+1 by height+1) */
float *resamp_clr_slice;/* 2D array for storing resampled colors
			   (width+1 by height+1) */
{
    int i, j;
    float WgtTL, WgtBL, WgtTR, WgtBR;
    float OpcAcc, RClrAcc, GClrAcc, BClrAcc;

    WgtTL = WgtTL_d;
    WgtBL = WgtBL_d;
    WgtTR = WgtTR_d;
    WgtBR = WgtBR_d;
    for (j = 0; j <= height; j++) {
	for (i = 0; i <= width; i++) {
	    OpcAcc = 0.;
	    RClrAcc = 0.;
	    GClrAcc = 0.;
	    BClrAcc = 0.;
	    if (i > 0 && j > 0) {
		OpcAcc += WgtTL * opc_slice[-1-width];
		if (color_channels == 1) {
		    RClrAcc += WgtTL * clr_slice[-1-width];
		} else {
		    RClrAcc += WgtTL * clr_slice[3*(-1-width)];
		    GClrAcc += WgtTL * clr_slice[3*(-1-width)+1];
		    BClrAcc += WgtTL * clr_slice[3*(-1-width)+2];
		}
	    }
	    if (i > 0 && j < height) {
		OpcAcc += WgtBL * opc_slice[-1];
		if (color_channels == 1) {
		    RClrAcc += WgtBL * clr_slice[-1];
		} else {
		    RClrAcc += WgtBL * clr_slice[3*(-1)];
		    GClrAcc += WgtBL * clr_slice[3*(-1)+1];
		    BClrAcc += WgtBL * clr_slice[3*(-1)+2];
		}
	    }
	    if (i < width && j > 0) {
		OpcAcc += WgtTR * opc_slice[-width];
		if (color_channels == 1) {
		    RClrAcc += WgtTR * clr_slice[-width];
		} else {
		    RClrAcc += WgtTR * clr_slice[3*(-width)];
		    GClrAcc += WgtTR * clr_slice[3*(-width)+1];
		    BClrAcc += WgtTR * clr_slice[3*(-width)+2];
		}
	    }
	    if (i < width && j < height) {
		OpcAcc += WgtBR * opc_slice[0];
		if (color_channels == 1) {
		    RClrAcc += WgtBR * clr_slice[0];
		} else {
		    RClrAcc += WgtBR * clr_slice[3*(0)];
		    GClrAcc += WgtBR * clr_slice[3*(0)+1];
		    BClrAcc += WgtBR * clr_slice[3*(0)+2];
		}
	    }
	    *resamp_opc_slice = OpcAcc;
	    if (color_channels == 1) {
		*resamp_clr_slice = RClrAcc;
	    } else {
		resamp_clr_slice[0] = RClrAcc;
		resamp_clr_slice[1] = GClrAcc;
		resamp_clr_slice[2] = BClrAcc;
	    }
	    resamp_opc_slice++;
	    resamp_clr_slice += color_channels;
	    if (i != width) {
		opc_slice++;
		clr_slice += color_channels;;
	    }
	}
    }
}

/*
 * CompositeSlice
 *
 * Composite a resampled slice of voxels into the intermediate image.
 */

static void
CompositeSlice(resamp_opc, resamp_clr, width, height, color_channels,
	       int_image_ptr, int_image_width, min_opacity)
float *resamp_opc;	/* array of resampled opacities (width by height) */
float *resamp_clr;	/* array of resampled colors (width by height) */
int width;		/* size of resampled voxel arrays */
int height;
int color_channels;	/* number of color channels */
void *int_image_ptr;	/* pointer to intermediate image pixel corresponding
			   to top-left resampled voxel */
int int_image_width;	/* number of pixels in intermediate image scanline */
double min_opacity;	/* low opacity threshold */
{
    int i, j;
    float old_opc, old_r, old_g, old_b;
    float new_opc, new_r, new_g, new_b;
    GrayIntPixel *gray_intim;
    RGBIntPixel *rgb_intim;

    if (color_channels == 1)
	gray_intim = int_image_ptr;
    else
	rgb_intim = int_image_ptr;
    for (j = 0; j < height; j++) {
	for (i = 0; i < width; i++) {
	    if (*resamp_opc > min_opacity) {
		if (color_channels == 1) {
		    old_opc = gray_intim->opcflt;
		    old_r = gray_intim->clrflt;
		    new_opc = old_opc + *resamp_opc * ((float)1. - old_opc);
		    new_r = old_r + *resamp_clr * ((float)1. - old_opc);
		    gray_intim->opcflt = new_opc;
		    gray_intim->clrflt = new_r;
		} else {
		    old_opc = rgb_intim->opcflt;
		    old_r = rgb_intim->rclrflt;
		    old_g = rgb_intim->gclrflt;
		    old_b = rgb_intim->bclrflt;
		    new_opc = old_opc + *resamp_opc * ((float)1. - old_opc);
		    new_r = old_r + resamp_clr[0] * ((float)1. - old_opc);
		    new_g = old_g + resamp_clr[1] * ((float)1. - old_opc);
		    new_b = old_b + resamp_clr[2] * ((float)1. - old_opc);
		    rgb_intim->opcflt = new_opc;
		    rgb_intim->rclrflt = new_r;
		    rgb_intim->gclrflt = new_g;
		    rgb_intim->bclrflt = new_b;
		}
		
	    }
	    resamp_opc++;
	    if (color_channels == 1) {
		resamp_clr++;
		gray_intim++;
	    } else {
		resamp_clr += 3;
		rgb_intim++;
	    }
	} /* for i */
	if (color_channels == 1)
	    gray_intim += int_image_width - width;
	else
	    rgb_intim += int_image_width - width;
    } /* for j */
}

/*
 * AffineBruteForceWarp
 *
 * Warp the intermediate image into the final image (brute-force version,
 * affine transformations only).
 */

static void
AffineBruteForceWarp(vpc)
vpContext *vpc;
{
    unsigned char *int_image;	/* pointer to start of intermediate image
				   (GrayIntPixel or RGBIntPixel) */
    int int_width;		/* size of intermediate image */
    int int_height;
    int int_scanbytes;		/* bytes per scanline in intermediate image */
    unsigned char *image;	/* final image pixel */
    int i, j;			/* coordinates of final image pixel */
    float int_i_flt, int_j_flt;	/* position of final image pixel in
				   intermediate image coordinates */
    float int_i_int, int_j_int;	/* truncated int_i_flt, int_j_flt */
    int int_i, int_j;		/* integer int_i_int, int_j_int */
    double alpha_i, alpha_j;	/* separable interpolation weights */
    double wgt;			/* interpolation weight */
    GrayIntPixel *gray_pix;	/* intermediate image pixel (grayscale) */
    RGBIntPixel *rgb_pix;	/* intermediate image pixel (RGB) */
    double denom;
    double ma, mb, mc, md, me, mf;
    float r, g, b, alpha;
    int r_int, g_int, b_int, alpha_int;
    int color_channels;		/* number of color channels in int. image */
    int pixel_type;		/* type of output image pixel */

    /* initialize */
    color_channels = vpc->color_channels;
    pixel_type = vpc->pixel_type;
    int_width = vpc->intermediate_width;
    int_height = vpc->intermediate_height;
    if (vpc->color_channels == 1) {
	int_image = (unsigned char *)vpc->int_image.gray_intim;
	if (vpc->pad_int_to_maxwidth)
	    int_scanbytes = vpc->max_intermediate_width*sizeof(GrayIntPixel);
	else
	    int_scanbytes = vpc->intermediate_width*sizeof(GrayIntPixel);
    } else {
	int_image = (unsigned char *)vpc->int_image.rgb_intim;
	if (vpc->pad_int_to_maxwidth)
	    int_scanbytes = vpc->max_intermediate_width*sizeof(RGBIntPixel);
	else
	    int_scanbytes = vpc->intermediate_width*sizeof(RGBIntPixel);
    }

    /* compute transformation from final image pixel to intermediate
       image pixel */
    denom = 1. / (vpc->warp_2d[0][0]*vpc->warp_2d[1][1] -
		  vpc->warp_2d[0][1]*vpc->warp_2d[1][0]);
    ma = vpc->warp_2d[1][1] * denom;
    mb = -vpc->warp_2d[0][1] * denom;
    mc = (vpc->warp_2d[0][1]*vpc->warp_2d[1][2] - 
	  vpc->warp_2d[1][1]*vpc->warp_2d[0][2]) * denom;
    md = -vpc->warp_2d[1][0] * denom;
    me = vpc->warp_2d[0][0] * denom;
    mf = (vpc->warp_2d[1][0]*vpc->warp_2d[0][2] - 
	  vpc->warp_2d[0][0]*vpc->warp_2d[1][2]) * denom;

    /* loop over the pixels of the final image */
    for (j = 0; j < vpc->image_height; j++) {
	image = (unsigned char *)vpc->image + j*vpc->image_bytes_per_scan;
	for (i = 0; i < vpc->image_width; i++) {
	    /* reverse-map final image pixel into intermediate image */
	    int_i_flt = ma*i + mb*j + mc;
	    int_j_flt = md*i + me*j + mf;

	    /* compute interpolation weights */
	    int_i_int = floor(int_i_flt);
	    int_j_int = floor(int_j_flt);
	    alpha_i = int_i_flt - int_i_int;
	    alpha_j = int_j_flt - int_j_int;
	    int_i = (int)int_i_int;
	    int_j = (int)int_j_int;

	    /* interpolate */
	    r = 0;
	    g = 0;
	    b = 0;
	    alpha = 0;
	    if (int_i >= 0 && int_i < int_width &&
		int_j >= 0 && int_j < int_height) {
		wgt = (1. - alpha_i) * (1. - alpha_j);
		if (color_channels == 1) {
		    gray_pix = (GrayIntPixel *)(int_image + int_j*
						int_scanbytes) + int_i;
		    r += gray_pix->clrflt*wgt;
		    alpha += gray_pix->opcflt*wgt;
		} else {
		    rgb_pix = (RGBIntPixel *)(int_image + int_j*
					      int_scanbytes) + int_i;
		    r += rgb_pix->rclrflt*wgt;
		    g += rgb_pix->gclrflt*wgt;
		    b += rgb_pix->bclrflt*wgt;
		    alpha += rgb_pix->opcflt*wgt;
		}
	    }
	    if (int_i >= 0 && int_i < int_width &&
		int_j >= -1 && int_j < int_height-1) {
		wgt = (1. - alpha_i) * alpha_j;
		if (color_channels == 1) {
		    gray_pix = (GrayIntPixel *)(int_image + (int_j+1)*
						int_scanbytes) + int_i;
		    r += gray_pix->clrflt*wgt;
		    alpha += gray_pix->opcflt*wgt;
		} else {
		    rgb_pix = (RGBIntPixel *)(int_image + (int_j+1)*
					      int_scanbytes) + int_i;
		    r += rgb_pix->rclrflt*wgt;
		    g += rgb_pix->gclrflt*wgt;
		    b += rgb_pix->bclrflt*wgt;
		    alpha += rgb_pix->opcflt*wgt;
		}
	    }
	    if (int_i >= -1 && int_i < int_width-1 &&
		int_j >= 0 && int_j < int_height) {
		wgt = alpha_i * (1. - alpha_j);
		if (color_channels == 1) {
		    gray_pix = (GrayIntPixel *)(int_image + int_j*
						int_scanbytes) + int_i+1;
		    r += gray_pix->clrflt*wgt;
		    alpha += gray_pix->opcflt*wgt;
		} else {
		    rgb_pix = (RGBIntPixel *)(int_image + int_j*
					      int_scanbytes) + int_i+1;
		    r += rgb_pix->rclrflt*wgt;
		    g += rgb_pix->gclrflt*wgt;
		    b += rgb_pix->bclrflt*wgt;
		    alpha += rgb_pix->opcflt*wgt;
		}
	    }
	    if (int_i >= -1 && int_i < int_width-1 &&
		int_j >= -1 && int_j < int_height-1) {
		wgt = alpha_i * alpha_j;
		if (color_channels == 1) {
		    gray_pix = (GrayIntPixel *)(int_image + (int_j+1)*
						int_scanbytes) + int_i+1;
		    r += gray_pix->clrflt*wgt;
		    alpha += gray_pix->opcflt*wgt;
		} else {
		    rgb_pix = (RGBIntPixel *)(int_image + (int_j+1)*
					      int_scanbytes) + int_i+1;
		    r += rgb_pix->rclrflt*wgt;
		    g += rgb_pix->gclrflt*wgt;
		    b += rgb_pix->bclrflt*wgt;
		    alpha += rgb_pix->opcflt*wgt;
		}
	    }

	    /* clamp the pixel */
	    if (alpha > 255.)
		alpha_int = 255;
	    else if (alpha < 0.)
		alpha_int = 0;
	    else
		alpha_int = alpha;
	    if (r > 255.)
		r_int = 255;
	    else if (r < 0)
		r_int = 0;
	    else
		r_int = r;
	    
	    if (color_channels == 3) {
		if (g > 255.)
		    g_int = 255;
		else if (g < 0.)
		    g_int = 0;
		else
		    g_int = g;
		if (b > 255.)
		    b_int = 255;
		else if (b < 0.)
		    b_int = 0;
		else
		    b_int = b;
	    }

	    /* store the pixel */
	    switch (pixel_type) {
	    case VP_ALPHA:
		*image++ = alpha_int;
		break;
	    case VP_LUMINANCE:
		*image++ = r_int;
		break;
	    case VP_LUMINANCEA:
		*image++ = r_int;
		*image++ = alpha_int;
		break;
	    case VP_RGB:
		*image++ = r_int;
		*image++ = g_int;
		*image++ = b_int;
		break;
	    case VP_RGBA:
		*image++ = r_int;
		*image++ = g_int;
		*image++ = b_int;
		*image++ = alpha_int;
		break;
	    case VP_BGR:
		*image++ = b_int;
		*image++ = g_int;
		*image++ = r_int;
		break;
	    case VP_ABGR:
		*image++ = alpha_int;
		*image++ = b_int;
		*image++ = g_int;
		*image++ = r_int;
		break;
	    default:
		VPBug("bad pixel type");
	    }
	} /* for i */
    } /* for j */
}

#ifdef DEBUG
StoreFloatImage(data, width, height, scale, filename)
float *data;		/* array of input data */
int width, height;	/* size of array */
double scale;		/* factor for scaling pixel values */
char *filename;		/* name of file to store result */
{
    unsigned char *image, *imptr;
    int i, j;

    image = (unsigned char *)malloc(width*height);
    imptr = image;
    for (j = 0; j < height; j++) {
	for (i = 0; i < width; i++) {
	    *imptr++ = (int)rint(scale * *data++);
	}
    }
    VprWriteGrayscaleTIFF(filename, width, height, width, image);
    free(image);
}
#endif