00001 /*
00002  * vp_shade.c
00003  *
00004  * Routines to implement the Phong shading equation using a lookup-table.
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 /*
00034  * Lookup Table Shading and Normal Vector Encoding
00035  *
00036  * The shader defined in this file implements the Phong shading equation
00037  * (I = Ia + Id*(N.L) + Is*(N.H)^n) using lookup tables.  To use this
00038  * shader you must include a "normal" field in each voxel.  This field
00039  * is computed by preprocessing the volume to estimate a surface normal
00040  * vector for each voxel (using the central-difference gradient operator)
00041  * and then encoding the normal in an index (13 bits in this program).
00042  * An index is stored in the normal field of each voxel.   At rendering
00043  * time the index is used to look up a color for the voxel in a table.
00044  *
00045  * There are many ways a normal vector can be encoded using the 13 bits of
00046  * the index.  A straight-forward method is to use 6 bits for the x component
00047  * of the vector, 6 bits for the y component, and 1 bit to indicate the sign
00048  * of the z component.  Assuming that the vector is normalized the z
00049  * component can be reconstructed from x and y.  Unfortunately this method
00050  * results in an uneven distribution of code points: the distance between
00051  * exactly-representable vectors is much smaller near N = +z or -z than
00052  * N = +x, -x, +y or -z (where x, y and z are the unit vectors).  This
00053  * can result in significant quantization error near the "equator", or more
00054  * table storage than necessary near the "poles".
00055  *
00056  * The normal vector encoding scheme used here is derived from a recursive
00057  * tesselation of a regular octahedron (an eight-sided solid with
00058  * equilateral triangular faces).  Consider subdividing each triangular
00059  * face into four smaller equilateral triangles, and then subdividing those
00060  * triangles recursively until a sufficiently large number of triangles
00061  * have been generated.  The representable normal vectors are the vectors
00062  * connecting the center of the solid to the vertices of the triangles.
00063  * The distribution of these vectors is not perfectly uniform (the density
00064  * is lower at the "mid-latitudes"), but the variation is relatively small.
00065  *
00066  * Each normal vector is now assigned a unique index as follows.  Let the
00067  * origin be at the center of the solid, and the x, y and z axes each
00068  * intersect a vertex of the original octahedron.  Use one bit of the index
00069  * to indicate the sign of the z component of the normal.  Now project the
00070  * subdivided triangle vertices onto the x-y plane.  This forms a square
00071  * grid of dots rotated 45 degrees relative to the x-y axes.  Starting at
00072  * (x=0, y=max), assign sequential integers to each normal vector projection,
00073  * proceeding left-to-right and top-to-bottom.  Finally, append the sign bit
00074  * for the z component to the least-significant end of the integer to get
00075  * the normal vector encoding.
00076  *
00077  * This scheme is useful because it is easy to compute the vector components
00078  * from an index and conversely to find the closest index for a given vector,
00079  * yet the distribution of representable vectors is pretty uniform.
00080  *
00081  * XXX better method is to rotate 45 degrees (M = [1 -1 ; 1 1]) and then
00082  * assign points in scanline order; no lookup tables needed; implement this!
00083  *
00084  * The layout of the shading lookup table is as follows:
00085  *    float shade_table[MAX_NORMAL+1][materials][color_channels]
00086  * where materials is the number of materials and color_channels is 1
00087  * for grayscale intensities or 3 for RGB intensities (stored in R,G,B
00088  * order).
00089  */
00090 
00091 /* define normal index parameters; if you change these then you
00092    must change VP_NORM_MAX in volpack.h */
00093 #define NORM_13                 /* use 13 bit normals */
00094 #undef NORM_15                  /* don't use 15 bit normals */
00095 
00096 #ifdef NORM_13                  /* parameters for 13 bit normals */
00097 #define NORM_N          44      /*   N parameter for normal computation */
00098 #define NORM_BITS       13      /*   number of bits to encode a normal:
00099                                      (1+ceil(log2(2*(N*N+N+1)))) */
00100 #define MAX_NORMAL      7923    /*   maximum normal index */
00101 #endif
00102 
00103 #ifdef NORM_15                  /* parameters for 15 bit normals */
00104 #define NORM_N          90
00105 #define NORM_BITS       15
00106 #define MAX_NORMAL      16131
00107 #endif
00108 
00109 
00110 /* static lookup tables (computed only once, used for all vpContexts) */
00111 static int NormalTablesInitialized = 0; /* set to 1 after initialization */
00112 static short *NormalPy; /* NormalPy[py] = normal index for the normal
00113                                    whose projection in the x-y plane is (0,py)
00114                                    and whose z component is +
00115                                    (py = -NORM_N to +NORM_N) */
00116 static short NormalPyStorage[1+2*NORM_N];       /* storage for NormalPy */
00117 static short *NormalXLimit;     /* max abs(x) allowed for a given y */
00118 static short NormalXLimitStorage[1+2*NORM_N]; /* storage for NormalXLimit */
00119 static void InitNormalTables ANSI_ARGS((void));
00120 
00121 /*
00122  * InitNormalTables
00123  *
00124  * Initialize lookup tables for computing normal indices.
00125  */
00126 
00127 static void
00128 InitNormalTables()
00129 {
00130     int y, xcount, codepoint;
00131     int sum;
00132     double value;
00133 
00134     /* initialize NormalPy */
00135     xcount = 1;
00136     codepoint = 2;
00137     NormalPy = &NormalPyStorage[NORM_N];
00138     NormalXLimit = &NormalXLimitStorage[NORM_N];
00139     for (y = -NORM_N; y <= NORM_N; y++) {
00140         NormalPy[y] = codepoint + (((xcount-1)/2) << 1);
00141         codepoint += (xcount << 1);
00142         NormalXLimit[y] = (xcount-1)/2;
00143         if (y < 0)
00144             xcount += 2;
00145         else
00146             xcount -= 2;
00147     }
00148 
00149     NormalTablesInitialized = 1;
00150 }
00151 
00152 /*
00153  * vpNormalIndex
00154  *
00155  * Return the best normal index for the given normal vector.
00156  */
00157 
00158 int
00159 vpNormalIndex(nx, ny, nz)
00160 double nx, ny, nz;
00161 {
00162     int n, x, y, xlimit;
00163     double denom, denominv;
00164 
00165     if (!NormalTablesInitialized)
00166         InitNormalTables();
00167     denom = (nx < 0) ? -nx : nx;
00168     denom += (ny < 0) ? -ny : ny;
00169     denom += (nz < 0) ? -nz : nz;
00170     denominv = (double)NORM_N / denom;
00171     x = (int)rint(nx * denominv);
00172     y = (int)rint(ny * denominv);
00173 
00174     /* clamp x */
00175     xlimit = NormalXLimit[y];
00176     if (x < 0) {
00177         if (-x > xlimit)
00178             x = -xlimit;
00179     } else {
00180         if (x > xlimit)
00181             x = xlimit;
00182     }
00183 
00184     n = NormalPy[y] + (x << 1);
00185     if (nz < 0)
00186         n |= 1;
00187     ASSERT(n >= 0 && n <= VP_NORM_MAX);
00188     return(n);
00189 }
00190 
00191 /*
00192  * vpNormal
00193  *
00194  * Compute normal vector components given a normal vector index.
00195  */
00196 
00197 vpResult
00198 vpNormal(n, nx, ny, nz)
00199 int n;                  /* normal index */
00200 double *nx, *ny, *nz;   /* storage for result */
00201 {
00202     int py, px, pz, pxlimit2;
00203     double pxd, pyd, pzd, plength;
00204 
00205     if (!NormalTablesInitialized)
00206         InitNormalTables();
00207     for (py = -NORM_N; py <= NORM_N; py++) {
00208         pxlimit2 = 2 * ((py<0) ? (NORM_N + py) : (NORM_N - py));
00209         if (NormalPy[py] - pxlimit2 <= n &&
00210             NormalPy[py] + pxlimit2 + 1 >= n) {
00211             break;
00212         }
00213     }
00214     if (py > NORM_N) {
00215         return(VPERROR_BAD_VALUE);
00216     } else {
00217         px = (n - NormalPy[py]) >> 1;
00218         pz = NORM_N;
00219         if (px < 0)
00220             pz += px;
00221         else
00222             pz -= px;
00223         if (py < 0)
00224             pz += py;
00225         else
00226             pz -= py;
00227         if (n & 1)
00228             pz = -pz;
00229         pxd = (double)px;
00230         pyd = (double)py;
00231         pzd = (double)pz;
00232         plength = 1. / sqrt(pxd*pxd+pyd*pyd+pzd*pzd);
00233         *nx = pxd * plength;
00234         *ny = pyd * plength;
00235         *nz = pzd * plength;
00236     }
00237     return(VP_OK);
00238 }
00239 
00240 /*
00241  * vpScanlineNormals
00242  *
00243  * Compute normals and/or gradients for a scanline of voxels.
00244  */
00245 
00246 vpResult
00247 vpScanlineNormals(vpc, length, scalar_data, scalar_minus_y,
00248                   scalar_plus_y, scalar_minus_z, scalar_plus_z,
00249                   voxel_data, scalar_field, grad_field, norm_field)
00250 vpContext *vpc;         /* context */
00251 int length;             /* number of scalars in scanline */
00252 unsigned char *scalar_data;     /* scanline of scalar data */
00253 unsigned char *scalar_minus_y;  /* adjacent scanline of scalar data (-y) */
00254 unsigned char *scalar_plus_y;   /* adjacent scanline of scalar data (+y) */
00255 unsigned char *scalar_minus_z;  /* adjacent scanline of scalar data (-z) */
00256 unsigned char *scalar_plus_z;   /* adjacent scanline of scalar data (+z) */
00257 void *voxel_data;       /* location to store first voxel */
00258 int scalar_field;       /* voxel field for scalar, or VP_SKIP_FIELD */
00259 int grad_field;         /* voxel field for gradient, or VP_SKIP_FIELD */
00260 int norm_field;         /* voxel field for normal, or VP_SKIP_FIELD */
00261 {
00262     int x;                      /* voxel index */
00263     double grad_x;              /* components of the gradient vector */
00264     double grad_y;
00265     double grad_z;
00266     double twice_grad_mag;      /* twice the magnitude of the gradient */
00267     int grad;                   /* gradient magnitude */
00268     int norm;                   /* normal index */
00269     int edge;                   /* true if this scanline is on the edge
00270                                    of the volume */
00271     int voxel_size;             /* size of a voxel in bytes */
00272     int scalar_offset;          /* byte offset for scalar in voxel */
00273     int grad_offset;            /* byte offset for gradient in voxel */
00274     int norm_offset;            /* byte offset for normal in voxel */
00275     char *voxel;                /* pointer to current voxel */
00276     int retcode;
00277 
00278     /* check for errors */
00279     if ((retcode = VPCheckVoxelFields(vpc)) != VP_OK)
00280         return(retcode);
00281     if (scalar_field != VP_SKIP_FIELD) {
00282         if (scalar_field < 0 || scalar_field >= vpc->num_voxel_fields)
00283             return(VPSetError(vpc, VPERROR_BAD_VALUE));
00284         if (vpc->field_size[scalar_field] != VP_SCALAR_SIZE)
00285             return(VPSetError(vpc, VPERROR_BAD_VALUE));
00286         scalar_offset = vpc->field_offset[scalar_field];
00287     }
00288     if (grad_field != VP_SKIP_FIELD) {
00289         if (grad_field < 0 || grad_field >= vpc->num_voxel_fields)
00290             return(VPSetError(vpc, VPERROR_BAD_VALUE));
00291         if (vpc->field_size[grad_field] != VP_GRAD_SIZE)
00292             return(VPSetError(vpc, VPERROR_BAD_VALUE));
00293         grad_offset = vpc->field_offset[grad_field];
00294     }
00295     if (norm_field != VP_SKIP_FIELD) {
00296         if (norm_field < 0 || norm_field >= vpc->num_voxel_fields)
00297             return(VPSetError(vpc, VPERROR_BAD_VALUE));
00298         if (vpc->field_size[norm_field] != VP_NORM_SIZE)
00299             return(VPSetError(vpc, VPERROR_BAD_VALUE));
00300         norm_offset = vpc->field_offset[norm_field];
00301     }
00302     voxel_size = vpc->raw_bytes_per_voxel;
00303 
00304     /* compute the scanline */
00305     voxel = voxel_data;
00306     if (scalar_minus_y == NULL || scalar_plus_y == NULL ||
00307         scalar_minus_z == NULL || scalar_plus_z == NULL) {
00308         edge = 1;
00309     } else {
00310         edge = 0;
00311     }
00312     for (x = 0; x < length; x++) {
00313         /* compute gradient and normal for voxel x and store */
00314         if (edge || x == 0 || x == length-1) {
00315             if (scalar_field != VP_SKIP_FIELD)
00316                 ByteField(voxel, scalar_offset) = scalar_data[x];
00317             if (grad_field != VP_SKIP_FIELD)
00318                 ByteField(voxel, grad_offset) = 0;
00319             if (norm_field != VP_SKIP_FIELD)
00320                 ShortField(voxel, norm_offset) = 0;
00321         } else {
00322             grad_x = (int)scalar_data[x+1] - (int)scalar_data[x-1];
00323             grad_y = (int)scalar_plus_y[x] - (int)scalar_minus_y[x];
00324             grad_z = (int)scalar_plus_z[x] - (int)scalar_minus_z[x];
00325             twice_grad_mag = sqrt(grad_x*grad_x+grad_y*grad_y+grad_z*grad_z);
00326             if (scalar_field != VP_SKIP_FIELD)
00327                 ByteField(voxel, scalar_offset) = scalar_data[x];
00328             if (grad_field != VP_SKIP_FIELD) {
00329                 grad = (int)rint(0.5 * twice_grad_mag);
00330                 ASSERT(grad >= 0 && grad <= VP_GRAD_MAX);
00331                 ByteField(voxel, grad_offset) = grad;
00332             }
00333             if (norm_field != VP_SKIP_FIELD) {
00334                 if (twice_grad_mag < VP_EPS)
00335                     norm = 0;
00336                 else
00337                     norm = vpNormalIndex(grad_x / twice_grad_mag,
00338                                          grad_y / twice_grad_mag,
00339                                          grad_z / twice_grad_mag);
00340                 ShortField(voxel, norm_offset) = norm;
00341             }
00342         }
00343 
00344         /* go on to next voxel */
00345         voxel += voxel_size;
00346     }
00347     return(VP_OK);
00348 }
00349 
00350 /*
00351  * vpVolumeNormals
00352  *
00353  * Compute normals and/or gradients for a volume.  Result is stored
00354  * in raw_voxels in the current context.
00355  */
00356 
00357 vpResult
00358 vpVolumeNormals(vpc, scalar_data, length, scalar_field, grad_field, norm_field)
00359 vpContext *vpc;         /* context */
00360 unsigned char *scalar_data;     /* 3D array of scalar data */
00361 int length;             /* number of scalars in scalar_data */
00362 int scalar_field;       /* voxel field for scalar, or VP_SKIP_FIELD */
00363 int grad_field;         /* voxel field for gradient, or VP_SKIP_FIELD */
00364 int norm_field;         /* voxel field for normal, or VP_SKIP_FIELD */
00365 {
00366     int xlen, ylen, zlen;       /* volume dimensions */
00367     int y, z;                   /* loop indices */
00368     unsigned char *scalar;      /* pointer to current scalar */
00369     int scalar_ystride;         /* stride to next scalar scanline */
00370     int scalar_zstride;         /* stride to next scalar slice */
00371     char *voxel;                /* pointer to current voxel */
00372     int voxel_ystride;          /* stride to next voxel scanline */
00373     int voxel_zstride;          /* stride to next voxel slice */
00374     unsigned char *s_py, *s_my, *s_pz, *s_mz; /* ptrs to adjacent scans */
00375     int retcode;
00376 
00377     /* check for errors */
00378     if ((retcode = VPCheckRawVolume(vpc)) != VP_OK)
00379         return(retcode);
00380     xlen = vpc->xlen;
00381     ylen = vpc->ylen;
00382     zlen = vpc->zlen;
00383     if (xlen * ylen * zlen != length)
00384         return(VPSetError(vpc, VPERROR_BAD_SIZE));
00385 
00386     /* initialize */
00387     scalar = scalar_data;
00388     scalar_ystride = xlen;
00389     scalar_zstride = xlen*ylen;
00390     voxel = vpc->raw_voxels;
00391     voxel_ystride = vpc->ystride;
00392     voxel_zstride = vpc->zstride;
00393 
00394     /* compute volume data */
00395     for (z = 0; z < zlen; z++) {
00396         ReportStatus(vpc, (double)z / (double)zlen);
00397         for (y = 0; y < ylen; y++) {
00398             s_my = (y == 0)      ? NULL : scalar - scalar_ystride;
00399             s_py = (y == ylen-1) ? NULL : scalar + scalar_ystride;
00400             s_mz = (z == 0)      ? NULL : scalar - scalar_zstride;
00401             s_pz = (z == zlen-1) ? NULL : scalar + scalar_zstride;
00402             retcode = vpScanlineNormals(vpc, xlen, scalar, s_my, s_py,
00403                                         s_mz, s_pz, voxel, scalar_field,
00404                                         grad_field, norm_field);
00405             if (retcode != VP_OK)
00406                 return(retcode);
00407             scalar += scalar_ystride;
00408             voxel += voxel_ystride;
00409         }
00410         scalar += scalar_zstride - ylen*scalar_ystride;
00411         voxel += voxel_zstride - ylen*voxel_ystride;
00412     }
00413     ReportStatus(vpc, 1.0);
00414     return(VP_OK);
00415 }
00416 
00417 /*
00418  * vpShadeTable
00419  *
00420  * Compute a shading lookup table for the current lighting and
00421  * model matrix.
00422  */
00423 
00424 vpResult
00425 vpShadeTable(vpc)
00426 vpContext *vpc;
00427 {
00428     int num_lights;             /* number of enabled lights */
00429     vpVector3 light_color[VP_MAX_LIGHTS]; /* light colors */
00430     vpVector4 obj_light[VP_MAX_LIGHTS]; /* light_vector in object space */
00431     vpVector4 obj_highlight[VP_MAX_LIGHTS]; /* halfway-vector */
00432     vpVector4 obj_viewpoint;    /* viewpoint in object coordinates */
00433     vpMatrix4 a;                /* linear system matrix */
00434     double *rhs[VP_MAX_LIGHTS+1];/* right-hand-side/solution vectors */
00435     int px, py, pz;             /* code point coordinates (integers) */
00436     double pxd, pyd, pzd;       /* code point coordinates (doubles) */
00437     double plength;
00438     int pxlimit;                /* maximum absolute value of px for this py */
00439     double nx, ny, nz;          /* normal vector components */
00440     double n_dot_v_xy;          /* normal_vector dot obj_viewpoint (x&y) */
00441     double n_dot_v_z;           /* normal_vector dot obj_viewpoint (z) */
00442     double n_dot_v;             /* normal_vector dot obj_viewpoint */
00443     double n_dot_l_xy;          /* normal_vector dot light_vector (x&y) */
00444     double n_dot_l_z;           /* normal_vector dot light_vector (z) */
00445     double n_dot_l;             /* normal_vector dot light_vector */
00446     double n_dot_h_xy;          /* normal_vector dot halfway_vector (x&y) */
00447     double n_dot_h_z;           /* normal_vector dot halfway_vector (z) */
00448     double n_dot_h;             /* normal_vector dot halfway_vector */
00449     float r, g, b;              /* intensities to store in shade table */
00450     float *table_row;           /* pointer to table row for current normal */
00451     float *table;               /* pointer to table entry */
00452     float *shadow_table_row;    /* pointer to table row for current normal */
00453     float *shadow_table;        /* pointer to shadow table entry */
00454     int surface_side;           /* EXT_SURFACE or INT_SURFACE */
00455     int znegative;              /* true iff nz is negative */
00456     int light_both_sides;       /* use two-sided lighting */
00457     int reverse_surface_sides;  /* reverse interior and exterior */
00458     int color_channels;         /* number of color channels */
00459     int num_materials;          /* number of materials */
00460     int retcode;
00461     double *matl_props;         /* material properties */
00462     int enable_shadows;         /* true if shadows are enabled */
00463     int shadow_light;           /* light index for light casting shadows */
00464     int clamp;                  /* true if table entries should be clamped */
00465     int c, l, m;
00466 #ifdef DEBUG
00467     vpVector4 tmpv;
00468 #endif
00469     DECLARE_TIME(t0);
00470     DECLARE_TIME(t1);
00471 
00472     /* error checking */
00473     if (vpc->shading_mode != LOOKUP_SHADER)
00474         return(VP_OK);
00475     if ((retcode = VPCheckShader(vpc)) != VP_OK)
00476         return(retcode);
00477     if ((retcode = VPCheckShadows(vpc)) != VP_OK)
00478         return(retcode);
00479     ASSERT(vpc->color_channels == 1 || vpc->color_channels == 3);
00480 
00481     /* start timer */
00482     GET_TIME(vpc, t0);
00483 
00484     /* transform viewpoint vector and light vectors to object space */
00485     vpSetVector4(obj_viewpoint, 0., 0., 1., 1.);
00486     rhs[0] = obj_viewpoint;
00487     num_lights = 0;
00488     for (c = 0; c < VP_MAX_LIGHTS; c++) {
00489         if (vpc->light_enable[c]) {
00490             bcopy(vpc->light_color[c], light_color[num_lights],
00491                   sizeof(vpVector3));
00492             bcopy(vpc->light_vector[c], obj_light[num_lights],
00493                   sizeof(vpVector4));
00494             rhs[num_lights+1] = obj_light[num_lights];
00495             num_lights++;
00496         }
00497     }
00498     bcopy(vpc->transforms[VP_MODEL], a, sizeof(vpMatrix4));
00499     retcode = vpSolveSystem4(a, rhs, num_lights+1);
00500     if (retcode != VP_OK)
00501         return(retcode);
00502 
00503 #ifdef DEBUG
00504     /* make sure the solver gave the right answer */
00505     vpMatrixVectorMult4(tmpv, vpc->transforms[VP_MODEL], obj_viewpoint);
00506     if (fabs(tmpv[0]) > VP_EPS || fabs(tmpv[1]) > VP_EPS ||
00507         fabs(tmpv[2] - 1.) > VP_EPS || fabs(tmpv[3] - 1.) > VP_EPS) {
00508         printf("\n");
00509         printf("Modelview:\n");
00510         printf("    %12.8f %12.8f %12.8f %12.8f\n",
00511            vpc->transforms[VP_MODEL][0][0], vpc->transforms[VP_MODEL][0][1],
00512            vpc->transforms[VP_MODEL][0][2], vpc->transforms[VP_MODEL][0][3]);
00513         printf("    %12.8f %12.8f %12.8f %12.8f\n",
00514            vpc->transforms[VP_MODEL][1][0], vpc->transforms[VP_MODEL][1][1],
00515            vpc->transforms[VP_MODEL][1][2], vpc->transforms[VP_MODEL][1][3]);
00516         printf("    %12.8f %12.8f %12.8f %12.8f\n",
00517            vpc->transforms[VP_MODEL][2][0], vpc->transforms[VP_MODEL][2][1],
00518            vpc->transforms[VP_MODEL][2][2], vpc->transforms[VP_MODEL][2][3]);
00519         printf("    %12.8f %12.8f %12.8f %12.8f\n",
00520            vpc->transforms[VP_MODEL][3][0], vpc->transforms[VP_MODEL][3][1],
00521            vpc->transforms[VP_MODEL][3][2], vpc->transforms[VP_MODEL][3][3]);
00522         VPBug("SolveSystem failed on viewpoint");
00523     }
00524     l = 0;
00525     for (c = 0; c < VP_MAX_LIGHTS; c++) {
00526         if (vpc->light_enable[c]) {
00527             vpMatrixVectorMult4(tmpv, vpc->transforms[VP_MODEL], obj_light[l]);
00528             if (fabs(tmpv[0] - vpc->light_vector[c][0]) > VP_EPS ||
00529                 fabs(tmpv[1] - vpc->light_vector[c][1]) > VP_EPS ||
00530                 fabs(tmpv[2] - vpc->light_vector[c][2]) > VP_EPS ||
00531                 fabs(tmpv[3] - vpc->light_vector[c][3]) > VP_EPS)
00532                 VPBug("SolveSystem failed on light %d\n", c);
00533             l++;
00534         }
00535     }
00536 #endif
00537 
00538     /* compute highlight vectors */
00539     for (l = 0; l < num_lights; l++) {
00540         obj_highlight[l][0] = obj_light[l][0] + obj_viewpoint[0];
00541         obj_highlight[l][1] = obj_light[l][1] + obj_viewpoint[1];
00542         obj_highlight[l][2] = obj_light[l][2] + obj_viewpoint[2];
00543         vpNormalize3(obj_highlight[l]);
00544     }
00545 
00546     /* initialize options */
00547     light_both_sides = vpc->light_both_sides;
00548     reverse_surface_sides = vpc->reverse_surface_sides;
00549     color_channels = vpc->color_channels;
00550     num_materials = vpc->num_materials;
00551     table = vpc->shade_color_table;
00552     enable_shadows = vpc->enable_shadows;
00553     if (enable_shadows) {
00554         shadow_table = vpc->shadow_color_table;
00555         shadow_light = vpc->shadow_light_num - VP_LIGHT0;
00556         bzero(shadow_table, vpc->shadow_color_table_size);
00557     } else {
00558         shadow_table = NULL;
00559         shadow_light = -1;
00560     }
00561     clamp = vpc->clamp_shade_table;
00562 
00563     /* store shade table entries for the zero-vector */
00564     for (znegative = 0; znegative <= 1; znegative++) {
00565         if (znegative) {
00566             if (reverse_surface_sides)
00567                 surface_side = EXT_SURFACE;
00568             else
00569                 surface_side = INT_SURFACE;
00570         } else {
00571             if (reverse_surface_sides)
00572                 surface_side = INT_SURFACE;
00573             else
00574                 surface_side = EXT_SURFACE;
00575         }
00576         for (m = 0; m < num_materials; m++) {
00577             matl_props = vpc->matl_props[m][surface_side];
00578             *table++ = matl_props[MATL_AMB_R];
00579             if (color_channels == 3) {
00580                 *table++ = matl_props[MATL_AMB_G];
00581                 *table++ = matl_props[MATL_AMB_B];
00582             }
00583         }
00584     }
00585     table_row = table;
00586     if (enable_shadows) {
00587         for (znegative = 0; znegative <= 1; znegative++) {
00588             for (m = 0; m < num_materials; m++) {
00589                 *shadow_table++ = 0;
00590                 if (color_channels == 3) {
00591                     *shadow_table++ = 0;
00592                     *shadow_table++ = 0;
00593                 }
00594             }
00595         }
00596     }
00597     shadow_table_row = shadow_table;
00598 
00599     /* compute the shade table entries for nonzero normals */
00600     for (py = -NORM_N; py <= NORM_N; py++) {
00601         pxlimit = (py < 0) ? (NORM_N + py) : (NORM_N - py);
00602         pz = -1;
00603         pxd = (double)(-pxlimit-1);
00604         pyd = (double)py;
00605         pzd = (double)(-1);
00606         for (px = -pxlimit; px <= pxlimit; px++) {
00607             /* compute normal vector components for this code point */
00608             pxd += 1.0;
00609             if (px <= 0) {
00610                 pz++;
00611                 pzd += 1.0;
00612             } else {
00613                 pz--;
00614                 pzd -= 1.0;
00615             }
00616             plength = 1. / sqrt(pxd*pxd + pyd*pyd + pzd*pzd);
00617             nx = pxd * plength;
00618             ny = pyd * plength;
00619             nz = pzd * plength;
00620 
00621             /* compute n dot v (for determining surface side) */
00622             n_dot_v_xy = nx*obj_viewpoint[0] + ny*obj_viewpoint[1];
00623             n_dot_v_z = nz*obj_viewpoint[2];
00624 
00625             /* store ambient light terms */
00626             table = table_row;
00627             for (znegative = 0; znegative <= 1; znegative++) {
00628                 if (znegative)
00629                     n_dot_v = n_dot_v_xy - n_dot_v_z;
00630                 else
00631                     n_dot_v = n_dot_v_xy + n_dot_v_z;
00632                 if (reverse_surface_sides)
00633                     n_dot_v = -n_dot_v;
00634                 if (n_dot_v >= 0)
00635                     surface_side = EXT_SURFACE;
00636                 else
00637                     surface_side = INT_SURFACE;
00638                 for (m = 0; m < num_materials; m++) {
00639                     matl_props = vpc->matl_props[m][surface_side];
00640                     *table++ = matl_props[MATL_AMB_R];
00641                     if (color_channels == 3) {
00642                         *table++ = matl_props[MATL_AMB_G];
00643                         *table++ = matl_props[MATL_AMB_B];
00644                     }
00645                 }
00646             }
00647 
00648             /* loop over lights */
00649             for (l = 0; l < num_lights; l++) {
00650                 if (l == shadow_light)
00651                     table = shadow_table_row;
00652                 else
00653                     table = table_row;
00654 
00655                 /* compute n dot l and n dot h */
00656                 n_dot_l_xy = nx*obj_light[l][0] + ny*obj_light[l][1];
00657                 n_dot_l_z = nz*obj_light[l][2];
00658                 n_dot_h_xy = nx*obj_highlight[l][0] + ny*obj_highlight[l][1];
00659                 n_dot_h_z = nz*obj_highlight[l][2];
00660 
00661                 /* loop over the two signs for z */
00662                 for (znegative = 0; znegative <= 1; znegative++) {
00663                     if (znegative) {
00664                         n_dot_v = n_dot_v_xy - n_dot_v_z;
00665                         n_dot_l = n_dot_l_xy - n_dot_l_z;
00666                         n_dot_h = n_dot_h_xy - n_dot_h_z;
00667                     } else {
00668                         n_dot_v = n_dot_v_xy + n_dot_v_z;
00669                         n_dot_l = n_dot_l_xy + n_dot_l_z;
00670                         n_dot_h = n_dot_h_xy + n_dot_h_z;
00671                     }
00672                     if (reverse_surface_sides) {
00673                         n_dot_v = -n_dot_v;
00674                         n_dot_l = -n_dot_l;
00675                         n_dot_h = -n_dot_h;
00676                     }
00677                     if (n_dot_v >= 0)
00678                         surface_side = EXT_SURFACE;
00679                     else
00680                         surface_side = INT_SURFACE;
00681                     if (light_both_sides) {
00682                         n_dot_l = fabs(n_dot_l);
00683                         n_dot_h = fabs(n_dot_h);
00684                     } else if (surface_side == EXT_SURFACE) {
00685                         n_dot_l = MAX(n_dot_l, 0.0);
00686                         n_dot_h = MAX(n_dot_h, 0.0);
00687                     } else {
00688                         n_dot_l = MAX(-n_dot_l, 0.0);
00689                         n_dot_h = MAX(-n_dot_h, 0.0);
00690                     }
00691 
00692                     /* loop over material types */
00693                     for (m = 0; m < num_materials; m++) {
00694                         matl_props = vpc->matl_props[m][surface_side];
00695                         *table++ += light_color[l][0]*(matl_props[MATL_DIFF_R]*
00696                                     n_dot_l + matl_props[MATL_SPEC_R]*
00697                                     pow(n_dot_h, matl_props[MATL_SHINY]));
00698                         if (color_channels == 3) {
00699                             *table++ += light_color[l][1]*
00700                                     (matl_props[MATL_DIFF_G]*
00701                                     n_dot_l + matl_props[MATL_SPEC_G]*
00702                                     pow(n_dot_h, matl_props[MATL_SHINY]));
00703                             *table++ += light_color[l][2]*
00704                                     (matl_props[MATL_DIFF_B]*
00705                                     n_dot_l + matl_props[MATL_SPEC_B]*
00706                                     pow(n_dot_h, matl_props[MATL_SHINY]));
00707                         }
00708                     } /* for m */
00709                 } /* for znegative */
00710             } /* for l */
00711 
00712             /* clamp */
00713             if (clamp) {
00714                 if (enable_shadows) {
00715                     table = table_row;
00716                     shadow_table = shadow_table_row;
00717                     for (znegative = 0; znegative <= 1; znegative++) {
00718                         for (m = 0; m < num_materials; m++) {
00719                             for (c = 0; c < color_channels; c++) {
00720                                 if (*table > 255.)
00721                                     *table = 255.;
00722                                 if (*table + *shadow_table > 255.)
00723                                     *shadow_table = 255. - *table;
00724                                 shadow_table++;
00725                                 table++;
00726                             }
00727                         }
00728                     }
00729                 } else {
00730                     table = table_row;
00731                     for (znegative = 0; znegative <= 1; znegative++) {
00732                         for (m = 0; m < num_materials; m++) {
00733                             for (c = 0; c < color_channels; c++) {
00734                                 if (*table > 255.)
00735                                     *table = 255.;
00736                                 table++;
00737                             }
00738                         }
00739                     }
00740                 }
00741             }
00742 
00743             if (num_materials == 1) {
00744                 table_row += 2*color_channels;
00745             } else {
00746                 if (color_channels == 1)
00747                     table_row += 2*num_materials;
00748                 else
00749                     table_row += 6*num_materials;
00750             }
00751 
00752             if (enable_shadows) {
00753                 if (num_materials == 1) {
00754                     shadow_table_row += 2*color_channels;
00755                 } else {
00756                     if (color_channels == 1)
00757                         shadow_table_row += 2*num_materials;
00758                     else
00759                         shadow_table_row += 6*num_materials;
00760                 }
00761             }
00762         } /* for px */
00763     } /* for py */
00764 
00765     /* stop timer */
00766     GET_TIME(vpc, t1);
00767     STORE_TIME(vpc, VPTIMER_SHADE, t0, t1);
00768 
00769     return(VP_OK);
00770 }