/* * vp_shade.c * * Routines to implement the Phong shading equation using a lookup-table. * * 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$ * $Revision$ */ #include "vp_global.h" /* * Lookup Table Shading and Normal Vector Encoding * * The shader defined in this file implements the Phong shading equation * (I = Ia + Id*(N.L) + Is*(N.H)^n) using lookup tables. To use this * shader you must include a "normal" field in each voxel. This field * is computed by preprocessing the volume to estimate a surface normal * vector for each voxel (using the central-difference gradient operator) * and then encoding the normal in an index (13 bits in this program). * An index is stored in the normal field of each voxel. At rendering * time the index is used to look up a color for the voxel in a table. * * There are many ways a normal vector can be encoded using the 13 bits of * the index. A straight-forward method is to use 6 bits for the x component * of the vector, 6 bits for the y component, and 1 bit to indicate the sign * of the z component. Assuming that the vector is normalized the z * component can be reconstructed from x and y. Unfortunately this method * results in an uneven distribution of code points: the distance between * exactly-representable vectors is much smaller near N = +z or -z than * N = +x, -x, +y or -z (where x, y and z are the unit vectors). This * can result in significant quantization error near the "equator", or more * table storage than necessary near the "poles". * * The normal vector encoding scheme used here is derived from a recursive * tesselation of a regular octahedron (an eight-sided solid with * equilateral triangular faces). Consider subdividing each triangular * face into four smaller equilateral triangles, and then subdividing those * triangles recursively until a sufficiently large number of triangles * have been generated. The representable normal vectors are the vectors * connecting the center of the solid to the vertices of the triangles. * The distribution of these vectors is not perfectly uniform (the density * is lower at the "mid-latitudes"), but the variation is relatively small. * * Each normal vector is now assigned a unique index as follows. Let the * origin be at the center of the solid, and the x, y and z axes each * intersect a vertex of the original octahedron. Use one bit of the index * to indicate the sign of the z component of the normal. Now project the * subdivided triangle vertices onto the x-y plane. This forms a square * grid of dots rotated 45 degrees relative to the x-y axes. Starting at * (x=0, y=max), assign sequential integers to each normal vector projection, * proceeding left-to-right and top-to-bottom. Finally, append the sign bit * for the z component to the least-significant end of the integer to get * the normal vector encoding. * * This scheme is useful because it is easy to compute the vector components * from an index and conversely to find the closest index for a given vector, * yet the distribution of representable vectors is pretty uniform. * * XXX better method is to rotate 45 degrees (M = [1 -1 ; 1 1]) and then * assign points in scanline order; no lookup tables needed; implement this! * * The layout of the shading lookup table is as follows: * float shade_table[MAX_NORMAL+1][materials][color_channels] * where materials is the number of materials and color_channels is 1 * for grayscale intensities or 3 for RGB intensities (stored in R,G,B * order). */ /* define normal index parameters; if you change these then you must change VP_NORM_MAX in volpack.h */ #define NORM_13 /* use 13 bit normals */ #undef NORM_15 /* don't use 15 bit normals */ #ifdef NORM_13 /* parameters for 13 bit normals */ #define NORM_N 44 /* N parameter for normal computation */ #define NORM_BITS 13 /* number of bits to encode a normal: (1+ceil(log2(2*(N*N+N+1)))) */ #define MAX_NORMAL 7923 /* maximum normal index */ #endif #ifdef NORM_15 /* parameters for 15 bit normals */ #define NORM_N 90 #define NORM_BITS 15 #define MAX_NORMAL 16131 #endif /* static lookup tables (computed only once, used for all vpContexts) */ static int NormalTablesInitialized = 0; /* set to 1 after initialization */ static short *NormalPy; /* NormalPy[py] = normal index for the normal whose projection in the x-y plane is (0,py) and whose z component is + (py = -NORM_N to +NORM_N) */ static short NormalPyStorage[1+2*NORM_N]; /* storage for NormalPy */ static short *NormalXLimit; /* max abs(x) allowed for a given y */ static short NormalXLimitStorage[1+2*NORM_N]; /* storage for NormalXLimit */ static void InitNormalTables ANSI_ARGS((void)); /* * InitNormalTables * * Initialize lookup tables for computing normal indices. */ static void InitNormalTables() { int y, xcount, codepoint; int sum; double value; /* initialize NormalPy */ xcount = 1; codepoint = 2; NormalPy = &NormalPyStorage[NORM_N]; NormalXLimit = &NormalXLimitStorage[NORM_N]; for (y = -NORM_N; y <= NORM_N; y++) { NormalPy[y] = codepoint + (((xcount-1)/2) << 1); codepoint += (xcount << 1); NormalXLimit[y] = (xcount-1)/2; if (y < 0) xcount += 2; else xcount -= 2; } NormalTablesInitialized = 1; } /* * vpNormalIndex * * Return the best normal index for the given normal vector. */ int vpNormalIndex(nx, ny, nz) double nx, ny, nz; { int n, x, y, xlimit; double denom, denominv; if (!NormalTablesInitialized) InitNormalTables(); denom = (nx < 0) ? -nx : nx; denom += (ny < 0) ? -ny : ny; denom += (nz < 0) ? -nz : nz; denominv = (double)NORM_N / denom; x = (int)rint(nx * denominv); y = (int)rint(ny * denominv); /* clamp x */ xlimit = NormalXLimit[y]; if (x < 0) { if (-x > xlimit) x = -xlimit; } else { if (x > xlimit) x = xlimit; } n = NormalPy[y] + (x << 1); if (nz < 0) n |= 1; ASSERT(n >= 0 && n <= VP_NORM_MAX); return(n); } /* * vpNormal * * Compute normal vector components given a normal vector index. */ vpResult vpNormal(n, nx, ny, nz) int n; /* normal index */ double *nx, *ny, *nz; /* storage for result */ { int py, px, pz, pxlimit2; double pxd, pyd, pzd, plength; if (!NormalTablesInitialized) InitNormalTables(); for (py = -NORM_N; py <= NORM_N; py++) { pxlimit2 = 2 * ((py<0) ? (NORM_N + py) : (NORM_N - py)); if (NormalPy[py] - pxlimit2 <= n && NormalPy[py] + pxlimit2 + 1 >= n) { break; } } if (py > NORM_N) { return(VPERROR_BAD_VALUE); } else { px = (n - NormalPy[py]) >> 1; pz = NORM_N; if (px < 0) pz += px; else pz -= px; if (py < 0) pz += py; else pz -= py; if (n & 1) pz = -pz; pxd = (double)px; pyd = (double)py; pzd = (double)pz; plength = 1. / sqrt(pxd*pxd+pyd*pyd+pzd*pzd); *nx = pxd * plength; *ny = pyd * plength; *nz = pzd * plength; } return(VP_OK); } /* * vpScanlineNormals * * Compute normals and/or gradients for a scanline of voxels. */ vpResult vpScanlineNormals(vpc, length, scalar_data, scalar_minus_y, scalar_plus_y, scalar_minus_z, scalar_plus_z, voxel_data, scalar_field, grad_field, norm_field) vpContext *vpc; /* context */ int length; /* number of scalars in scanline */ unsigned char *scalar_data; /* scanline of scalar data */ unsigned char *scalar_minus_y; /* adjacent scanline of scalar data (-y) */ unsigned char *scalar_plus_y; /* adjacent scanline of scalar data (+y) */ unsigned char *scalar_minus_z; /* adjacent scanline of scalar data (-z) */ unsigned char *scalar_plus_z; /* adjacent scanline of scalar data (+z) */ void *voxel_data; /* location to store first voxel */ int scalar_field; /* voxel field for scalar, or VP_SKIP_FIELD */ int grad_field; /* voxel field for gradient, or VP_SKIP_FIELD */ int norm_field; /* voxel field for normal, or VP_SKIP_FIELD */ { int x; /* voxel index */ double grad_x; /* components of the gradient vector */ double grad_y; double grad_z; double twice_grad_mag; /* twice the magnitude of the gradient */ int grad; /* gradient magnitude */ int norm; /* normal index */ int edge; /* true if this scanline is on the edge of the volume */ int voxel_size; /* size of a voxel in bytes */ int scalar_offset=0; /* byte offset for scalar in voxel */ int grad_offset=0; /* byte offset for gradient in voxel */ int norm_offset=0; /* byte offset for normal in voxel */ char *voxel; /* pointer to current voxel */ int retcode; /* check for errors */ if ((retcode = VPCheckVoxelFields(vpc)) != VP_OK) return(retcode); if (scalar_field != VP_SKIP_FIELD) { if (scalar_field < 0 || scalar_field >= vpc->num_voxel_fields) return(VPSetError(vpc, VPERROR_BAD_VALUE)); if (vpc->field_size[scalar_field] != VP_SCALAR_SIZE) return(VPSetError(vpc, VPERROR_BAD_VALUE)); scalar_offset = vpc->field_offset[scalar_field]; } if (grad_field != VP_SKIP_FIELD) { if (grad_field < 0 || grad_field >= vpc->num_voxel_fields) return(VPSetError(vpc, VPERROR_BAD_VALUE)); if (vpc->field_size[grad_field] != VP_GRAD_SIZE) return(VPSetError(vpc, VPERROR_BAD_VALUE)); grad_offset = vpc->field_offset[grad_field]; } if (norm_field != VP_SKIP_FIELD) { if (norm_field < 0 || norm_field >= vpc->num_voxel_fields) return(VPSetError(vpc, VPERROR_BAD_VALUE)); if (vpc->field_size[norm_field] != VP_NORM_SIZE) return(VPSetError(vpc, VPERROR_BAD_VALUE)); norm_offset = vpc->field_offset[norm_field]; } voxel_size = vpc->raw_bytes_per_voxel; /* compute the scanline */ voxel = voxel_data; if (scalar_minus_y == NULL || scalar_plus_y == NULL || scalar_minus_z == NULL || scalar_plus_z == NULL) { edge = 1; } else { edge = 0; } for (x = 0; x < length; x++) { /* compute gradient and normal for voxel x and store */ if (edge || x == 0 || x == length-1) { if (scalar_field != VP_SKIP_FIELD) ByteField(voxel, scalar_offset) = scalar_data[x]; if (grad_field != VP_SKIP_FIELD) ByteField(voxel, grad_offset) = 0; if (norm_field != VP_SKIP_FIELD) ShortField(voxel, norm_offset) = 0; } else { grad_x = (int)scalar_data[x+1] - (int)scalar_data[x-1]; grad_y = (int)scalar_plus_y[x] - (int)scalar_minus_y[x]; grad_z = (int)scalar_plus_z[x] - (int)scalar_minus_z[x]; twice_grad_mag = sqrt(grad_x*grad_x+grad_y*grad_y+grad_z*grad_z); if (scalar_field != VP_SKIP_FIELD) ByteField(voxel, scalar_offset) = scalar_data[x]; if (grad_field != VP_SKIP_FIELD) { grad = (int)rint(0.5 * twice_grad_mag); ASSERT(grad >= 0 && grad <= VP_GRAD_MAX); ByteField(voxel, grad_offset) = grad; } if (norm_field != VP_SKIP_FIELD) { if (twice_grad_mag < VP_EPS) norm = 0; else norm = vpNormalIndex(grad_x / twice_grad_mag, grad_y / twice_grad_mag, grad_z / twice_grad_mag); ShortField(voxel, norm_offset) = norm; } } /* go on to next voxel */ voxel += voxel_size; } return(VP_OK); } /* * vpVolumeNormals * * Compute normals and/or gradients for a volume. Result is stored * in raw_voxels in the current context. */ vpResult vpVolumeNormals(vpc, scalar_data, length, scalar_field, grad_field, norm_field) vpContext *vpc; /* context */ unsigned char *scalar_data; /* 3D array of scalar data */ int length; /* number of scalars in scalar_data */ int scalar_field; /* voxel field for scalar, or VP_SKIP_FIELD */ int grad_field; /* voxel field for gradient, or VP_SKIP_FIELD */ int norm_field; /* voxel field for normal, or VP_SKIP_FIELD */ { int xlen, ylen, zlen; /* volume dimensions */ int y, z; /* loop indices */ unsigned char *scalar; /* pointer to current scalar */ int scalar_ystride; /* stride to next scalar scanline */ int scalar_zstride; /* stride to next scalar slice */ char *voxel; /* pointer to current voxel */ int voxel_ystride; /* stride to next voxel scanline */ int voxel_zstride; /* stride to next voxel slice */ unsigned char *s_py, *s_my, *s_pz, *s_mz; /* ptrs to adjacent scans */ int retcode; /* check for errors */ if ((retcode = VPCheckRawVolume(vpc)) != VP_OK) return(retcode); xlen = vpc->xlen; ylen = vpc->ylen; zlen = vpc->zlen; if (xlen * ylen * zlen != length) return(VPSetError(vpc, VPERROR_BAD_SIZE)); /* initialize */ scalar = scalar_data; scalar_ystride = xlen; scalar_zstride = xlen*ylen; voxel = vpc->raw_voxels; voxel_ystride = vpc->ystride; voxel_zstride = vpc->zstride; /* compute volume data */ for (z = 0; z < zlen; z++) { ReportStatus(vpc, (double)z / (double)zlen); for (y = 0; y < ylen; y++) { s_my = (y == 0) ? NULL : scalar - scalar_ystride; s_py = (y == ylen-1) ? NULL : scalar + scalar_ystride; s_mz = (z == 0) ? NULL : scalar - scalar_zstride; s_pz = (z == zlen-1) ? NULL : scalar + scalar_zstride; retcode = vpScanlineNormals(vpc, xlen, scalar, s_my, s_py, s_mz, s_pz, voxel, scalar_field, grad_field, norm_field); if (retcode != VP_OK) return(retcode); scalar += scalar_ystride; voxel += voxel_ystride; } scalar += scalar_zstride - ylen*scalar_ystride; voxel += voxel_zstride - ylen*voxel_ystride; } ReportStatus(vpc, 1.0); return(VP_OK); } /* * vpShadeTable * * Compute a shading lookup table for the current lighting and * model matrix. */ vpResult vpShadeTable(vpc) vpContext *vpc; { int num_lights; /* number of enabled lights */ vpVector3 light_color[VP_MAX_LIGHTS]; /* light colors */ vpVector4 obj_light[VP_MAX_LIGHTS]; /* light_vector in object space */ vpVector4 obj_highlight[VP_MAX_LIGHTS]; /* halfway-vector */ vpVector4 obj_viewpoint; /* viewpoint in object coordinates */ vpMatrix4 a; /* linear system matrix */ double *rhs[VP_MAX_LIGHTS+1];/* right-hand-side/solution vectors */ int px, py, pz; /* code point coordinates (integers) */ double pxd, pyd, pzd; /* code point coordinates (doubles) */ double plength; int pxlimit; /* maximum absolute value of px for this py */ double nx, ny, nz; /* normal vector components */ double n_dot_v_xy; /* normal_vector dot obj_viewpoint (x&y) */ double n_dot_v_z; /* normal_vector dot obj_viewpoint (z) */ double n_dot_v; /* normal_vector dot obj_viewpoint */ double n_dot_l_xy; /* normal_vector dot light_vector (x&y) */ double n_dot_l_z; /* normal_vector dot light_vector (z) */ double n_dot_l; /* normal_vector dot light_vector */ double n_dot_h_xy; /* normal_vector dot halfway_vector (x&y) */ double n_dot_h_z; /* normal_vector dot halfway_vector (z) */ double n_dot_h; /* normal_vector dot halfway_vector */ float r, g, b; /* intensities to store in shade table */ float *table_row; /* pointer to table row for current normal */ float *table; /* pointer to table entry */ float *shadow_table_row; /* pointer to table row for current normal */ float *shadow_table; /* pointer to shadow table entry */ int surface_side; /* EXT_SURFACE or INT_SURFACE */ int znegative; /* true iff nz is negative */ int light_both_sides; /* use two-sided lighting */ int reverse_surface_sides; /* reverse interior and exterior */ int color_channels; /* number of color channels */ int num_materials; /* number of materials */ int retcode; double *matl_props; /* material properties */ int enable_shadows; /* true if shadows are enabled */ int shadow_light; /* light index for light casting shadows */ int clamp; /* true if table entries should be clamped */ int c, l, m; #ifdef DEBUG vpVector4 tmpv; #endif DECLARE_TIME(t0); DECLARE_TIME(t1); /* error checking */ if (vpc->shading_mode != LOOKUP_SHADER) return(VP_OK); if ((retcode = VPCheckShader(vpc)) != VP_OK) return(retcode); if ((retcode = VPCheckShadows(vpc)) != VP_OK) return(retcode); ASSERT(vpc->color_channels == 1 || vpc->color_channels == 3); /* start timer */ GET_TIME(vpc, t0); /* transform viewpoint vector and light vectors to object space */ vpSetVector4(obj_viewpoint, 0., 0., 1., 1.); rhs[0] = obj_viewpoint; num_lights = 0; for (c = 0; c < VP_MAX_LIGHTS; c++) { if (vpc->light_enable[c]) { bcopy(vpc->light_color[c], light_color[num_lights], sizeof(vpVector3)); bcopy(vpc->light_vector[c], obj_light[num_lights], sizeof(vpVector4)); rhs[num_lights+1] = obj_light[num_lights]; num_lights++; } } bcopy(vpc->transforms[VP_MODEL], a, sizeof(vpMatrix4)); retcode = vpSolveSystem4(a, rhs, num_lights+1); if (retcode != VP_OK) return(retcode); #ifdef DEBUG /* make sure the solver gave the right answer */ vpMatrixVectorMult4(tmpv, vpc->transforms[VP_MODEL], obj_viewpoint); if (fabs(tmpv[0]) > VP_EPS || fabs(tmpv[1]) > VP_EPS || fabs(tmpv[2] - 1.) > VP_EPS || fabs(tmpv[3] - 1.) > VP_EPS) { printf("\n"); printf("Modelview:\n"); printf(" %12.8f %12.8f %12.8f %12.8f\n", vpc->transforms[VP_MODEL][0][0], vpc->transforms[VP_MODEL][0][1], vpc->transforms[VP_MODEL][0][2], vpc->transforms[VP_MODEL][0][3]); printf(" %12.8f %12.8f %12.8f %12.8f\n", vpc->transforms[VP_MODEL][1][0], vpc->transforms[VP_MODEL][1][1], vpc->transforms[VP_MODEL][1][2], vpc->transforms[VP_MODEL][1][3]); printf(" %12.8f %12.8f %12.8f %12.8f\n", vpc->transforms[VP_MODEL][2][0], vpc->transforms[VP_MODEL][2][1], vpc->transforms[VP_MODEL][2][2], vpc->transforms[VP_MODEL][2][3]); printf(" %12.8f %12.8f %12.8f %12.8f\n", vpc->transforms[VP_MODEL][3][0], vpc->transforms[VP_MODEL][3][1], vpc->transforms[VP_MODEL][3][2], vpc->transforms[VP_MODEL][3][3]); VPBug("SolveSystem failed on viewpoint"); } l = 0; for (c = 0; c < VP_MAX_LIGHTS; c++) { if (vpc->light_enable[c]) { vpMatrixVectorMult4(tmpv, vpc->transforms[VP_MODEL], obj_light[l]); if (fabs(tmpv[0] - vpc->light_vector[c][0]) > VP_EPS || fabs(tmpv[1] - vpc->light_vector[c][1]) > VP_EPS || fabs(tmpv[2] - vpc->light_vector[c][2]) > VP_EPS || fabs(tmpv[3] - vpc->light_vector[c][3]) > VP_EPS) VPBug("SolveSystem failed on light %d\n", c); l++; } } #endif /* compute highlight vectors */ for (l = 0; l < num_lights; l++) { obj_highlight[l][0] = obj_light[l][0] + obj_viewpoint[0]; obj_highlight[l][1] = obj_light[l][1] + obj_viewpoint[1]; obj_highlight[l][2] = obj_light[l][2] + obj_viewpoint[2]; vpNormalize3(obj_highlight[l]); } /* initialize options */ light_both_sides = vpc->light_both_sides; reverse_surface_sides = vpc->reverse_surface_sides; color_channels = vpc->color_channels; num_materials = vpc->num_materials; table = vpc->shade_color_table; enable_shadows = vpc->enable_shadows; if (enable_shadows) { shadow_table = vpc->shadow_color_table; shadow_light = vpc->shadow_light_num - VP_LIGHT0; bzero(shadow_table, vpc->shadow_color_table_size); } else { shadow_table = NULL; shadow_light = -1; } clamp = vpc->clamp_shade_table; /* store shade table entries for the zero-vector */ for (znegative = 0; znegative <= 1; znegative++) { if (znegative) { if (reverse_surface_sides) surface_side = EXT_SURFACE; else surface_side = INT_SURFACE; } else { if (reverse_surface_sides) surface_side = INT_SURFACE; else surface_side = EXT_SURFACE; } for (m = 0; m < num_materials; m++) { matl_props = vpc->matl_props[m][surface_side]; *table++ = matl_props[MATL_AMB_R]; if (color_channels == 3) { *table++ = matl_props[MATL_AMB_G]; *table++ = matl_props[MATL_AMB_B]; } } } table_row = table; if (enable_shadows) { for (znegative = 0; znegative <= 1; znegative++) { for (m = 0; m < num_materials; m++) { *shadow_table++ = 0; if (color_channels == 3) { *shadow_table++ = 0; *shadow_table++ = 0; } } } } shadow_table_row = shadow_table; /* compute the shade table entries for nonzero normals */ for (py = -NORM_N; py <= NORM_N; py++) { pxlimit = (py < 0) ? (NORM_N + py) : (NORM_N - py); pz = -1; pxd = (double)(-pxlimit-1); pyd = (double)py; pzd = (double)(-1); for (px = -pxlimit; px <= pxlimit; px++) { /* compute normal vector components for this code point */ pxd += 1.0; if (px <= 0) { pz++; pzd += 1.0; } else { pz--; pzd -= 1.0; } plength = 1. / sqrt(pxd*pxd + pyd*pyd + pzd*pzd); nx = pxd * plength; ny = pyd * plength; nz = pzd * plength; /* compute n dot v (for determining surface side) */ n_dot_v_xy = nx*obj_viewpoint[0] + ny*obj_viewpoint[1]; n_dot_v_z = nz*obj_viewpoint[2]; /* store ambient light terms */ table = table_row; for (znegative = 0; znegative <= 1; znegative++) { if (znegative) n_dot_v = n_dot_v_xy - n_dot_v_z; else n_dot_v = n_dot_v_xy + n_dot_v_z; if (reverse_surface_sides) n_dot_v = -n_dot_v; if (n_dot_v >= 0) surface_side = EXT_SURFACE; else surface_side = INT_SURFACE; for (m = 0; m < num_materials; m++) { matl_props = vpc->matl_props[m][surface_side]; *table++ = matl_props[MATL_AMB_R]; if (color_channels == 3) { *table++ = matl_props[MATL_AMB_G]; *table++ = matl_props[MATL_AMB_B]; } } } /* loop over lights */ for (l = 0; l < num_lights; l++) { if (l == shadow_light) table = shadow_table_row; else table = table_row; /* compute n dot l and n dot h */ n_dot_l_xy = nx*obj_light[l][0] + ny*obj_light[l][1]; n_dot_l_z = nz*obj_light[l][2]; n_dot_h_xy = nx*obj_highlight[l][0] + ny*obj_highlight[l][1]; n_dot_h_z = nz*obj_highlight[l][2]; /* loop over the two signs for z */ for (znegative = 0; znegative <= 1; znegative++) { if (znegative) { n_dot_v = n_dot_v_xy - n_dot_v_z; n_dot_l = n_dot_l_xy - n_dot_l_z; n_dot_h = n_dot_h_xy - n_dot_h_z; } else { n_dot_v = n_dot_v_xy + n_dot_v_z; n_dot_l = n_dot_l_xy + n_dot_l_z; n_dot_h = n_dot_h_xy + n_dot_h_z; } if (reverse_surface_sides) { n_dot_v = -n_dot_v; n_dot_l = -n_dot_l; n_dot_h = -n_dot_h; } if (n_dot_v >= 0) surface_side = EXT_SURFACE; else surface_side = INT_SURFACE; if (light_both_sides) { n_dot_l = fabs(n_dot_l); n_dot_h = fabs(n_dot_h); } else if (surface_side == EXT_SURFACE) { n_dot_l = MAX(n_dot_l, 0.0); n_dot_h = MAX(n_dot_h, 0.0); } else { n_dot_l = MAX(-n_dot_l, 0.0); n_dot_h = MAX(-n_dot_h, 0.0); } /* loop over material types */ for (m = 0; m < num_materials; m++) { matl_props = vpc->matl_props[m][surface_side]; *table++ += light_color[l][0]*(matl_props[MATL_DIFF_R]* n_dot_l + matl_props[MATL_SPEC_R]* pow(n_dot_h, matl_props[MATL_SHINY])); if (color_channels == 3) { *table++ += light_color[l][1]* (matl_props[MATL_DIFF_G]* n_dot_l + matl_props[MATL_SPEC_G]* pow(n_dot_h, matl_props[MATL_SHINY])); *table++ += light_color[l][2]* (matl_props[MATL_DIFF_B]* n_dot_l + matl_props[MATL_SPEC_B]* pow(n_dot_h, matl_props[MATL_SHINY])); } } /* for m */ } /* for znegative */ } /* for l */ /* clamp */ if (clamp) { if (enable_shadows) { table = table_row; shadow_table = shadow_table_row; for (znegative = 0; znegative <= 1; znegative++) { for (m = 0; m < num_materials; m++) { for (c = 0; c < color_channels; c++) { if (*table > 255.) *table = 255.; if (*table + *shadow_table > 255.) *shadow_table = 255. - *table; shadow_table++; table++; } } } } else { table = table_row; for (znegative = 0; znegative <= 1; znegative++) { for (m = 0; m < num_materials; m++) { for (c = 0; c < color_channels; c++) { if (*table > 255.) *table = 255.; table++; } } } } } if (num_materials == 1) { table_row += 2*color_channels; } else { if (color_channels == 1) table_row += 2*num_materials; else table_row += 6*num_materials; } if (enable_shadows) { if (num_materials == 1) { shadow_table_row += 2*color_channels; } else { if (color_channels == 1) shadow_table_row += 2*num_materials; else shadow_table_row += 6*num_materials; } } } /* for px */ } /* for py */ /* stop timer */ GET_TIME(vpc, t1); STORE_TIME(vpc, VPTIMER_SHADE, t0, t1); return(VP_OK); }