1 /*
2 * vp_shade.c
3 *
4 * Routines to implement the Phong shading equation using a lookup-table.
5 *
6 * Copyright (c) 1994 The Board of Trustees of The Leland Stanford
7 * Junior University. All rights reserved.
8 *
9 * Permission to use, copy, modify and distribute this software and its
10 * documentation for any purpose is hereby granted without fee, provided
11 * that the above copyright notice and this permission notice appear in
12 * all copies of this software and that you do not sell the software.
13 * Commercial licensing is available by contacting the author.
14 *
15 * THE SOFTWARE IS PROVIDED "AS IS" AND WITHOUT WARRANTY OF ANY KIND,
16 * EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
17 * WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
18 *
19 * Author:
20 * Phil Lacroute
21 * Computer Systems Laboratory
22 * Electrical Engineering Dept.
23 * Stanford University
24 */
25
26 /*
27 * $Date: 1994/12/30 23:52:38 $
28 * $Revision: 1.27 $
29 */
30
31 #include "vp_global.h"
32
33 /*
34 * Lookup Table Shading and Normal Vector Encoding
35 *
36 * The shader defined in this file implements the Phong shading equation
37 * (I = Ia + Id*(N.L) + Is*(N.H)^n) using lookup tables. To use this
38 * shader you must include a "normal" field in each voxel. This field
39 * is computed by preprocessing the volume to estimate a surface normal
40 * vector for each voxel (using the central-difference gradient operator)
41 * and then encoding the normal in an index (13 bits in this program).
42 * An index is stored in the normal field of each voxel. At rendering
43 * time the index is used to look up a color for the voxel in a table.
44 *
45 * There are many ways a normal vector can be encoded using the 13 bits of
46 * the index. A straight-forward method is to use 6 bits for the x component
47 * of the vector, 6 bits for the y component, and 1 bit to indicate the sign
48 * of the z component. Assuming that the vector is normalized the z
49 * component can be reconstructed from x and y. Unfortunately this method
50 * results in an uneven distribution of code points: the distance between
51 * exactly-representable vectors is much smaller near N = +z or -z than
52 * N = +x, -x, +y or -z (where x, y and z are the unit vectors). This
53 * can result in significant quantization error near the "equator", or more
54 * table storage than necessary near the "poles".
55 *
56 * The normal vector encoding scheme used here is derived from a recursive
57 * tesselation of a regular octahedron (an eight-sided solid with
58 * equilateral triangular faces). Consider subdividing each triangular
59 * face into four smaller equilateral triangles, and then subdividing those
60 * triangles recursively until a sufficiently large number of triangles
61 * have been generated. The representable normal vectors are the vectors
62 * connecting the center of the solid to the vertices of the triangles.
63 * The distribution of these vectors is not perfectly uniform (the density
64 * is lower at the "mid-latitudes"), but the variation is relatively small.
65 *
66 * Each normal vector is now assigned a unique index as follows. Let the
67 * origin be at the center of the solid, and the x, y and z axes each
68 * intersect a vertex of the original octahedron. Use one bit of the index
69 * to indicate the sign of the z component of the normal. Now project the
70 * subdivided triangle vertices onto the x-y plane. This forms a square
71 * grid of dots rotated 45 degrees relative to the x-y axes. Starting at
72 * (x=0, y=max), assign sequential integers to each normal vector projection,
73 * proceeding left-to-right and top-to-bottom. Finally, append the sign bit
74 * for the z component to the least-significant end of the integer to get
75 * the normal vector encoding.
76 *
77 * This scheme is useful because it is easy to compute the vector components
78 * from an index and conversely to find the closest index for a given vector,
79 * yet the distribution of representable vectors is pretty uniform.
80 *
81 * XXX better method is to rotate 45 degrees (M = [1 -1 ; 1 1]) and then
82 * assign points in scanline order; no lookup tables needed; implement this!
83 *
84 * The layout of the shading lookup table is as follows:
85 * float shade_table[MAX_NORMAL+1][materials][color_channels]
86 * where materials is the number of materials and color_channels is 1
87 * for grayscale intensities or 3 for RGB intensities (stored in R,G,B
88 * order).
89 */
90
91 /* define normal index parameters; if you change these then you
92 must change VP_NORM_MAX in volpack.h */
93 #define NORM_13 /* use 13 bit normals */
94 #undef NORM_15 /* don't use 15 bit normals */
95
96 #ifdef NORM_13 /* parameters for 13 bit normals */
97 #define NORM_N 44 /* N parameter for normal computation */
98 #define NORM_BITS 13 /* number of bits to encode a normal:
99 (1+ceil(log2(2*(N*N+N+1)))) */
100 #define MAX_NORMAL 7923 /* maximum normal index */
101 #endif
102
103 #ifdef NORM_15 /* parameters for 15 bit normals */
104 #define NORM_N 90
105 #define NORM_BITS 15
106 #define MAX_NORMAL 16131
107 #endif
108
109
110 /* static lookup tables (computed only once, used for all vpContexts) */
111 static int NormalTablesInitialized = 0; /* set to 1 after initialization */
112 static short *NormalPy; /* NormalPy[py] = normal index for the normal
113 whose projection in the x-y plane is (0,py)
114 and whose z component is +
115 (py = -NORM_N to +NORM_N) */
116 static short NormalPyStorage[1+2*NORM_N]; /* storage for NormalPy */
117 static short *NormalXLimit; /* max abs(x) allowed for a given y */
118 static short NormalXLimitStorage[1+2*NORM_N]; /* storage for NormalXLimit */
119 static void InitNormalTables ANSI_ARGS((void));
120
121 /*
122 * InitNormalTables
123 *
124 * Initialize lookup tables for computing normal indices.
125 */
126
127 static void
InitNormalTables()128 InitNormalTables()
129 {
130 int y, xcount, codepoint;
131 int sum;
132 double value;
133
134 /* initialize NormalPy */
135 xcount = 1;
136 codepoint = 2;
137 NormalPy = &NormalPyStorage[NORM_N];
138 NormalXLimit = &NormalXLimitStorage[NORM_N];
139 for (y = -NORM_N; y <= NORM_N; y++) {
140 NormalPy[y] = codepoint + (((xcount-1)/2) << 1);
141 codepoint += (xcount << 1);
142 NormalXLimit[y] = (xcount-1)/2;
143 if (y < 0)
144 xcount += 2;
145 else
146 xcount -= 2;
147 }
148
149 NormalTablesInitialized = 1;
150 }
151
152 /*
153 * vpNormalIndex
154 *
155 * Return the best normal index for the given normal vector.
156 */
157
158 int
vpNormalIndex(nx,ny,nz)159 vpNormalIndex(nx, ny, nz)
160 double nx, ny, nz;
161 {
162 int n, x, y, xlimit;
163 double denom, denominv;
164
165 if (!NormalTablesInitialized)
166 InitNormalTables();
167 denom = (nx < 0) ? -nx : nx;
168 denom += (ny < 0) ? -ny : ny;
169 denom += (nz < 0) ? -nz : nz;
170 denominv = (double)NORM_N / denom;
171 x = (int)rint(nx * denominv);
172 y = (int)rint(ny * denominv);
173
174 /* clamp x */
175 xlimit = NormalXLimit[y];
176 if (x < 0) {
177 if (-x > xlimit)
178 x = -xlimit;
179 } else {
180 if (x > xlimit)
181 x = xlimit;
182 }
183
184 n = NormalPy[y] + (x << 1);
185 if (nz < 0)
186 n |= 1;
187 ASSERT(n >= 0 && n <= VP_NORM_MAX);
188 return(n);
189 }
190
191 /*
192 * vpNormal
193 *
194 * Compute normal vector components given a normal vector index.
195 */
196
197 vpResult
vpNormal(n,nx,ny,nz)198 vpNormal(n, nx, ny, nz)
199 int n; /* normal index */
200 double *nx, *ny, *nz; /* storage for result */
201 {
202 int py, px, pz, pxlimit2;
203 double pxd, pyd, pzd, plength;
204
205 if (!NormalTablesInitialized)
206 InitNormalTables();
207 for (py = -NORM_N; py <= NORM_N; py++) {
208 pxlimit2 = 2 * ((py<0) ? (NORM_N + py) : (NORM_N - py));
209 if (NormalPy[py] - pxlimit2 <= n &&
210 NormalPy[py] + pxlimit2 + 1 >= n) {
211 break;
212 }
213 }
214 if (py > NORM_N) {
215 return(VPERROR_BAD_VALUE);
216 } else {
217 px = (n - NormalPy[py]) >> 1;
218 pz = NORM_N;
219 if (px < 0)
220 pz += px;
221 else
222 pz -= px;
223 if (py < 0)
224 pz += py;
225 else
226 pz -= py;
227 if (n & 1)
228 pz = -pz;
229 pxd = (double)px;
230 pyd = (double)py;
231 pzd = (double)pz;
232 plength = 1. / sqrt(pxd*pxd+pyd*pyd+pzd*pzd);
233 *nx = pxd * plength;
234 *ny = pyd * plength;
235 *nz = pzd * plength;
236 }
237 return(VP_OK);
238 }
239
240 /*
241 * vpScanlineNormals
242 *
243 * Compute normals and/or gradients for a scanline of voxels.
244 */
245
246 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)247 vpScanlineNormals(vpc, length, scalar_data, scalar_minus_y,
248 scalar_plus_y, scalar_minus_z, scalar_plus_z,
249 voxel_data, scalar_field, grad_field, norm_field)
250 vpContext *vpc; /* context */
251 int length; /* number of scalars in scanline */
252 unsigned char *scalar_data; /* scanline of scalar data */
253 unsigned char *scalar_minus_y; /* adjacent scanline of scalar data (-y) */
254 unsigned char *scalar_plus_y; /* adjacent scanline of scalar data (+y) */
255 unsigned char *scalar_minus_z; /* adjacent scanline of scalar data (-z) */
256 unsigned char *scalar_plus_z; /* adjacent scanline of scalar data (+z) */
257 void *voxel_data; /* location to store first voxel */
258 int scalar_field; /* voxel field for scalar, or VP_SKIP_FIELD */
259 int grad_field; /* voxel field for gradient, or VP_SKIP_FIELD */
260 int norm_field; /* voxel field for normal, or VP_SKIP_FIELD */
261 {
262 int x; /* voxel index */
263 double grad_x; /* components of the gradient vector */
264 double grad_y;
265 double grad_z;
266 double twice_grad_mag; /* twice the magnitude of the gradient */
267 int grad; /* gradient magnitude */
268 int norm; /* normal index */
269 int edge; /* true if this scanline is on the edge
270 of the volume */
271 int voxel_size; /* size of a voxel in bytes */
272 int scalar_offset; /* byte offset for scalar in voxel */
273 int grad_offset; /* byte offset for gradient in voxel */
274 int norm_offset; /* byte offset for normal in voxel */
275 char *voxel; /* pointer to current voxel */
276 int retcode;
277
278 /* check for errors */
279 if ((retcode = VPCheckVoxelFields(vpc)) != VP_OK)
280 return(retcode);
281 if (scalar_field != VP_SKIP_FIELD) {
282 if (scalar_field < 0 || scalar_field >= vpc->num_voxel_fields)
283 return(VPSetError(vpc, VPERROR_BAD_VALUE));
284 if (vpc->field_size[scalar_field] != VP_SCALAR_SIZE)
285 return(VPSetError(vpc, VPERROR_BAD_VALUE));
286 scalar_offset = vpc->field_offset[scalar_field];
287 }
288 if (grad_field != VP_SKIP_FIELD) {
289 if (grad_field < 0 || grad_field >= vpc->num_voxel_fields)
290 return(VPSetError(vpc, VPERROR_BAD_VALUE));
291 if (vpc->field_size[grad_field] != VP_GRAD_SIZE)
292 return(VPSetError(vpc, VPERROR_BAD_VALUE));
293 grad_offset = vpc->field_offset[grad_field];
294 }
295 if (norm_field != VP_SKIP_FIELD) {
296 if (norm_field < 0 || norm_field >= vpc->num_voxel_fields)
297 return(VPSetError(vpc, VPERROR_BAD_VALUE));
298 if (vpc->field_size[norm_field] != VP_NORM_SIZE)
299 return(VPSetError(vpc, VPERROR_BAD_VALUE));
300 norm_offset = vpc->field_offset[norm_field];
301 }
302 voxel_size = vpc->raw_bytes_per_voxel;
303
304 /* compute the scanline */
305 voxel = voxel_data;
306 if (scalar_minus_y == NULL || scalar_plus_y == NULL ||
307 scalar_minus_z == NULL || scalar_plus_z == NULL) {
308 edge = 1;
309 } else {
310 edge = 0;
311 }
312 for (x = 0; x < length; x++) {
313 /* compute gradient and normal for voxel x and store */
314 if (edge || x == 0 || x == length-1) {
315 if (scalar_field != VP_SKIP_FIELD)
316 ByteField(voxel, scalar_offset) = scalar_data[x];
317 if (grad_field != VP_SKIP_FIELD)
318 ByteField(voxel, grad_offset) = 0;
319 if (norm_field != VP_SKIP_FIELD)
320 ShortField(voxel, norm_offset) = 0;
321 } else {
322 grad_x = (int)scalar_data[x+1] - (int)scalar_data[x-1];
323 grad_y = (int)scalar_plus_y[x] - (int)scalar_minus_y[x];
324 grad_z = (int)scalar_plus_z[x] - (int)scalar_minus_z[x];
325 twice_grad_mag = sqrt(grad_x*grad_x+grad_y*grad_y+grad_z*grad_z);
326 if (scalar_field != VP_SKIP_FIELD)
327 ByteField(voxel, scalar_offset) = scalar_data[x];
328 if (grad_field != VP_SKIP_FIELD) {
329 grad = (int)rint(0.5 * twice_grad_mag);
330 ASSERT(grad >= 0 && grad <= VP_GRAD_MAX);
331 ByteField(voxel, grad_offset) = grad;
332 }
333 if (norm_field != VP_SKIP_FIELD) {
334 if (twice_grad_mag < VP_EPS)
335 norm = 0;
336 else
337 norm = vpNormalIndex(grad_x / twice_grad_mag,
338 grad_y / twice_grad_mag,
339 grad_z / twice_grad_mag);
340 ShortField(voxel, norm_offset) = norm;
341 }
342 }
343
344 /* go on to next voxel */
345 voxel += voxel_size;
346 }
347 return(VP_OK);
348 }
349
350 /*
351 * vpVolumeNormals
352 *
353 * Compute normals and/or gradients for a volume. Result is stored
354 * in raw_voxels in the current context.
355 */
356
357 vpResult
vpVolumeNormals(vpc,scalar_data,length,scalar_field,grad_field,norm_field)358 vpVolumeNormals(vpc, scalar_data, length, scalar_field, grad_field, norm_field)
359 vpContext *vpc; /* context */
360 unsigned char *scalar_data; /* 3D array of scalar data */
361 int length; /* number of scalars in scalar_data */
362 int scalar_field; /* voxel field for scalar, or VP_SKIP_FIELD */
363 int grad_field; /* voxel field for gradient, or VP_SKIP_FIELD */
364 int norm_field; /* voxel field for normal, or VP_SKIP_FIELD */
365 {
366 int xlen, ylen, zlen; /* volume dimensions */
367 int y, z; /* loop indices */
368 unsigned char *scalar; /* pointer to current scalar */
369 int scalar_ystride; /* stride to next scalar scanline */
370 int scalar_zstride; /* stride to next scalar slice */
371 char *voxel; /* pointer to current voxel */
372 int voxel_ystride; /* stride to next voxel scanline */
373 int voxel_zstride; /* stride to next voxel slice */
374 unsigned char *s_py, *s_my, *s_pz, *s_mz; /* ptrs to adjacent scans */
375 int retcode;
376
377 /* check for errors */
378 if ((retcode = VPCheckRawVolume(vpc)) != VP_OK)
379 return(retcode);
380 xlen = vpc->xlen;
381 ylen = vpc->ylen;
382 zlen = vpc->zlen;
383 if (xlen * ylen * zlen != length)
384 return(VPSetError(vpc, VPERROR_BAD_SIZE));
385
386 /* initialize */
387 scalar = scalar_data;
388 scalar_ystride = xlen;
389 scalar_zstride = xlen*ylen;
390 voxel = vpc->raw_voxels;
391 voxel_ystride = vpc->ystride;
392 voxel_zstride = vpc->zstride;
393
394 /* compute volume data */
395 for (z = 0; z < zlen; z++) {
396 ReportStatus(vpc, (double)z / (double)zlen);
397 for (y = 0; y < ylen; y++) {
398 s_my = (y == 0) ? NULL : scalar - scalar_ystride;
399 s_py = (y == ylen-1) ? NULL : scalar + scalar_ystride;
400 s_mz = (z == 0) ? NULL : scalar - scalar_zstride;
401 s_pz = (z == zlen-1) ? NULL : scalar + scalar_zstride;
402 retcode = vpScanlineNormals(vpc, xlen, scalar, s_my, s_py,
403 s_mz, s_pz, voxel, scalar_field,
404 grad_field, norm_field);
405 if (retcode != VP_OK)
406 return(retcode);
407 scalar += scalar_ystride;
408 voxel += voxel_ystride;
409 }
410 scalar += scalar_zstride - ylen*scalar_ystride;
411 voxel += voxel_zstride - ylen*voxel_ystride;
412 }
413 ReportStatus(vpc, 1.0);
414 return(VP_OK);
415 }
416
417 /*
418 * vpShadeTable
419 *
420 * Compute a shading lookup table for the current lighting and
421 * model matrix.
422 */
423
424 vpResult
vpShadeTable(vpc)425 vpShadeTable(vpc)
426 vpContext *vpc;
427 {
428 int num_lights; /* number of enabled lights */
429 vpVector3 light_color[VP_MAX_LIGHTS]; /* light colors */
430 vpVector4 obj_light[VP_MAX_LIGHTS]; /* light_vector in object space */
431 vpVector4 obj_highlight[VP_MAX_LIGHTS]; /* halfway-vector */
432 vpVector4 obj_viewpoint; /* viewpoint in object coordinates */
433 vpMatrix4 a; /* linear system matrix */
434 double *rhs[VP_MAX_LIGHTS+1];/* right-hand-side/solution vectors */
435 int px, py, pz; /* code point coordinates (integers) */
436 double pxd, pyd, pzd; /* code point coordinates (doubles) */
437 double plength;
438 int pxlimit; /* maximum absolute value of px for this py */
439 double nx, ny, nz; /* normal vector components */
440 double n_dot_v_xy; /* normal_vector dot obj_viewpoint (x&y) */
441 double n_dot_v_z; /* normal_vector dot obj_viewpoint (z) */
442 double n_dot_v; /* normal_vector dot obj_viewpoint */
443 double n_dot_l_xy; /* normal_vector dot light_vector (x&y) */
444 double n_dot_l_z; /* normal_vector dot light_vector (z) */
445 double n_dot_l; /* normal_vector dot light_vector */
446 double n_dot_h_xy; /* normal_vector dot halfway_vector (x&y) */
447 double n_dot_h_z; /* normal_vector dot halfway_vector (z) */
448 double n_dot_h; /* normal_vector dot halfway_vector */
449 float r, g, b; /* intensities to store in shade table */
450 float *table_row; /* pointer to table row for current normal */
451 float *table; /* pointer to table entry */
452 float *shadow_table_row; /* pointer to table row for current normal */
453 float *shadow_table; /* pointer to shadow table entry */
454 int surface_side; /* EXT_SURFACE or INT_SURFACE */
455 int znegative; /* true iff nz is negative */
456 int light_both_sides; /* use two-sided lighting */
457 int reverse_surface_sides; /* reverse interior and exterior */
458 int color_channels; /* number of color channels */
459 int num_materials; /* number of materials */
460 int retcode;
461 double *matl_props; /* material properties */
462 int enable_shadows; /* true if shadows are enabled */
463 int shadow_light; /* light index for light casting shadows */
464 int clamp; /* true if table entries should be clamped */
465 int c, l, m;
466 #ifdef DEBUG
467 vpVector4 tmpv;
468 #endif
469 DECLARE_TIME(t0);
470 DECLARE_TIME(t1);
471
472 /* error checking */
473 if (vpc->shading_mode != LOOKUP_SHADER)
474 return(VP_OK);
475 if ((retcode = VPCheckShader(vpc)) != VP_OK)
476 return(retcode);
477 if ((retcode = VPCheckShadows(vpc)) != VP_OK)
478 return(retcode);
479 ASSERT(vpc->color_channels == 1 || vpc->color_channels == 3);
480
481 /* start timer */
482 GET_TIME(vpc, t0);
483
484 /* transform viewpoint vector and light vectors to object space */
485 vpSetVector4(obj_viewpoint, 0., 0., 1., 1.);
486 rhs[0] = obj_viewpoint;
487 num_lights = 0;
488 for (c = 0; c < VP_MAX_LIGHTS; c++) {
489 if (vpc->light_enable[c]) {
490 bcopy(vpc->light_color[c], light_color[num_lights],
491 sizeof(vpVector3));
492 bcopy(vpc->light_vector[c], obj_light[num_lights],
493 sizeof(vpVector4));
494 rhs[num_lights+1] = obj_light[num_lights];
495 num_lights++;
496 }
497 }
498 bcopy(vpc->transforms[VP_MODEL], a, sizeof(vpMatrix4));
499 retcode = vpSolveSystem4(a, rhs, num_lights+1);
500 if (retcode != VP_OK)
501 return(retcode);
502
503 #ifdef DEBUG
504 /* make sure the solver gave the right answer */
505 vpMatrixVectorMult4(tmpv, vpc->transforms[VP_MODEL], obj_viewpoint);
506 if (fabs(tmpv[0]) > VP_EPS || fabs(tmpv[1]) > VP_EPS ||
507 fabs(tmpv[2] - 1.) > VP_EPS || fabs(tmpv[3] - 1.) > VP_EPS) {
508 printf("\n");
509 printf("Modelview:\n");
510 printf(" %12.8f %12.8f %12.8f %12.8f\n",
511 vpc->transforms[VP_MODEL][0][0], vpc->transforms[VP_MODEL][0][1],
512 vpc->transforms[VP_MODEL][0][2], vpc->transforms[VP_MODEL][0][3]);
513 printf(" %12.8f %12.8f %12.8f %12.8f\n",
514 vpc->transforms[VP_MODEL][1][0], vpc->transforms[VP_MODEL][1][1],
515 vpc->transforms[VP_MODEL][1][2], vpc->transforms[VP_MODEL][1][3]);
516 printf(" %12.8f %12.8f %12.8f %12.8f\n",
517 vpc->transforms[VP_MODEL][2][0], vpc->transforms[VP_MODEL][2][1],
518 vpc->transforms[VP_MODEL][2][2], vpc->transforms[VP_MODEL][2][3]);
519 printf(" %12.8f %12.8f %12.8f %12.8f\n",
520 vpc->transforms[VP_MODEL][3][0], vpc->transforms[VP_MODEL][3][1],
521 vpc->transforms[VP_MODEL][3][2], vpc->transforms[VP_MODEL][3][3]);
522 VPBug("SolveSystem failed on viewpoint");
523 }
524 l = 0;
525 for (c = 0; c < VP_MAX_LIGHTS; c++) {
526 if (vpc->light_enable[c]) {
527 vpMatrixVectorMult4(tmpv, vpc->transforms[VP_MODEL], obj_light[l]);
528 if (fabs(tmpv[0] - vpc->light_vector[c][0]) > VP_EPS ||
529 fabs(tmpv[1] - vpc->light_vector[c][1]) > VP_EPS ||
530 fabs(tmpv[2] - vpc->light_vector[c][2]) > VP_EPS ||
531 fabs(tmpv[3] - vpc->light_vector[c][3]) > VP_EPS)
532 VPBug("SolveSystem failed on light %d\n", c);
533 l++;
534 }
535 }
536 #endif
537
538 /* compute highlight vectors */
539 for (l = 0; l < num_lights; l++) {
540 obj_highlight[l][0] = obj_light[l][0] + obj_viewpoint[0];
541 obj_highlight[l][1] = obj_light[l][1] + obj_viewpoint[1];
542 obj_highlight[l][2] = obj_light[l][2] + obj_viewpoint[2];
543 vpNormalize3(obj_highlight[l]);
544 }
545
546 /* initialize options */
547 light_both_sides = vpc->light_both_sides;
548 reverse_surface_sides = vpc->reverse_surface_sides;
549 color_channels = vpc->color_channels;
550 num_materials = vpc->num_materials;
551 table = vpc->shade_color_table;
552 enable_shadows = vpc->enable_shadows;
553 if (enable_shadows) {
554 shadow_table = vpc->shadow_color_table;
555 shadow_light = vpc->shadow_light_num - VP_LIGHT0;
556 bzero(shadow_table, vpc->shadow_color_table_size);
557 } else {
558 shadow_table = NULL;
559 shadow_light = -1;
560 }
561 clamp = vpc->clamp_shade_table;
562
563 /* store shade table entries for the zero-vector */
564 for (znegative = 0; znegative <= 1; znegative++) {
565 if (znegative) {
566 if (reverse_surface_sides)
567 surface_side = EXT_SURFACE;
568 else
569 surface_side = INT_SURFACE;
570 } else {
571 if (reverse_surface_sides)
572 surface_side = INT_SURFACE;
573 else
574 surface_side = EXT_SURFACE;
575 }
576 for (m = 0; m < num_materials; m++) {
577 matl_props = vpc->matl_props[m][surface_side];
578 *table++ = matl_props[MATL_AMB_R];
579 if (color_channels == 3) {
580 *table++ = matl_props[MATL_AMB_G];
581 *table++ = matl_props[MATL_AMB_B];
582 }
583 }
584 }
585 table_row = table;
586 if (enable_shadows) {
587 for (znegative = 0; znegative <= 1; znegative++) {
588 for (m = 0; m < num_materials; m++) {
589 *shadow_table++ = 0;
590 if (color_channels == 3) {
591 *shadow_table++ = 0;
592 *shadow_table++ = 0;
593 }
594 }
595 }
596 }
597 shadow_table_row = shadow_table;
598
599 /* compute the shade table entries for nonzero normals */
600 for (py = -NORM_N; py <= NORM_N; py++) {
601 pxlimit = (py < 0) ? (NORM_N + py) : (NORM_N - py);
602 pz = -1;
603 pxd = (double)(-pxlimit-1);
604 pyd = (double)py;
605 pzd = (double)(-1);
606 for (px = -pxlimit; px <= pxlimit; px++) {
607 /* compute normal vector components for this code point */
608 pxd += 1.0;
609 if (px <= 0) {
610 pz++;
611 pzd += 1.0;
612 } else {
613 pz--;
614 pzd -= 1.0;
615 }
616 plength = 1. / sqrt(pxd*pxd + pyd*pyd + pzd*pzd);
617 nx = pxd * plength;
618 ny = pyd * plength;
619 nz = pzd * plength;
620
621 /* compute n dot v (for determining surface side) */
622 n_dot_v_xy = nx*obj_viewpoint[0] + ny*obj_viewpoint[1];
623 n_dot_v_z = nz*obj_viewpoint[2];
624
625 /* store ambient light terms */
626 table = table_row;
627 for (znegative = 0; znegative <= 1; znegative++) {
628 if (znegative)
629 n_dot_v = n_dot_v_xy - n_dot_v_z;
630 else
631 n_dot_v = n_dot_v_xy + n_dot_v_z;
632 if (reverse_surface_sides)
633 n_dot_v = -n_dot_v;
634 if (n_dot_v >= 0)
635 surface_side = EXT_SURFACE;
636 else
637 surface_side = INT_SURFACE;
638 for (m = 0; m < num_materials; m++) {
639 matl_props = vpc->matl_props[m][surface_side];
640 *table++ = matl_props[MATL_AMB_R];
641 if (color_channels == 3) {
642 *table++ = matl_props[MATL_AMB_G];
643 *table++ = matl_props[MATL_AMB_B];
644 }
645 }
646 }
647
648 /* loop over lights */
649 for (l = 0; l < num_lights; l++) {
650 if (l == shadow_light)
651 table = shadow_table_row;
652 else
653 table = table_row;
654
655 /* compute n dot l and n dot h */
656 n_dot_l_xy = nx*obj_light[l][0] + ny*obj_light[l][1];
657 n_dot_l_z = nz*obj_light[l][2];
658 n_dot_h_xy = nx*obj_highlight[l][0] + ny*obj_highlight[l][1];
659 n_dot_h_z = nz*obj_highlight[l][2];
660
661 /* loop over the two signs for z */
662 for (znegative = 0; znegative <= 1; znegative++) {
663 if (znegative) {
664 n_dot_v = n_dot_v_xy - n_dot_v_z;
665 n_dot_l = n_dot_l_xy - n_dot_l_z;
666 n_dot_h = n_dot_h_xy - n_dot_h_z;
667 } else {
668 n_dot_v = n_dot_v_xy + n_dot_v_z;
669 n_dot_l = n_dot_l_xy + n_dot_l_z;
670 n_dot_h = n_dot_h_xy + n_dot_h_z;
671 }
672 if (reverse_surface_sides) {
673 n_dot_v = -n_dot_v;
674 n_dot_l = -n_dot_l;
675 n_dot_h = -n_dot_h;
676 }
677 if (n_dot_v >= 0)
678 surface_side = EXT_SURFACE;
679 else
680 surface_side = INT_SURFACE;
681 if (light_both_sides) {
682 n_dot_l = fabs(n_dot_l);
683 n_dot_h = fabs(n_dot_h);
684 } else if (surface_side == EXT_SURFACE) {
685 n_dot_l = MAX(n_dot_l, 0.0);
686 n_dot_h = MAX(n_dot_h, 0.0);
687 } else {
688 n_dot_l = MAX(-n_dot_l, 0.0);
689 n_dot_h = MAX(-n_dot_h, 0.0);
690 }
691
692 /* loop over material types */
693 for (m = 0; m < num_materials; m++) {
694 matl_props = vpc->matl_props[m][surface_side];
695 *table++ += light_color[l][0]*(matl_props[MATL_DIFF_R]*
696 n_dot_l + matl_props[MATL_SPEC_R]*
697 pow(n_dot_h, matl_props[MATL_SHINY]));
698 if (color_channels == 3) {
699 *table++ += light_color[l][1]*
700 (matl_props[MATL_DIFF_G]*
701 n_dot_l + matl_props[MATL_SPEC_G]*
702 pow(n_dot_h, matl_props[MATL_SHINY]));
703 *table++ += light_color[l][2]*
704 (matl_props[MATL_DIFF_B]*
705 n_dot_l + matl_props[MATL_SPEC_B]*
706 pow(n_dot_h, matl_props[MATL_SHINY]));
707 }
708 } /* for m */
709 } /* for znegative */
710 } /* for l */
711
712 /* clamp */
713 if (clamp) {
714 if (enable_shadows) {
715 table = table_row;
716 shadow_table = shadow_table_row;
717 for (znegative = 0; znegative <= 1; znegative++) {
718 for (m = 0; m < num_materials; m++) {
719 for (c = 0; c < color_channels; c++) {
720 if (*table > 255.)
721 *table = 255.;
722 if (*table + *shadow_table > 255.)
723 *shadow_table = 255. - *table;
724 shadow_table++;
725 table++;
726 }
727 }
728 }
729 } else {
730 table = table_row;
731 for (znegative = 0; znegative <= 1; znegative++) {
732 for (m = 0; m < num_materials; m++) {
733 for (c = 0; c < color_channels; c++) {
734 if (*table > 255.)
735 *table = 255.;
736 table++;
737 }
738 }
739 }
740 }
741 }
742
743 if (num_materials == 1) {
744 table_row += 2*color_channels;
745 } else {
746 if (color_channels == 1)
747 table_row += 2*num_materials;
748 else
749 table_row += 6*num_materials;
750 }
751
752 if (enable_shadows) {
753 if (num_materials == 1) {
754 shadow_table_row += 2*color_channels;
755 } else {
756 if (color_channels == 1)
757 shadow_table_row += 2*num_materials;
758 else
759 shadow_table_row += 6*num_materials;
760 }
761 }
762 } /* for px */
763 } /* for py */
764
765 /* stop timer */
766 GET_TIME(vpc, t1);
767 STORE_TIME(vpc, VPTIMER_SHADE, t0, t1);
768
769 return(VP_OK);
770 }
771