1 /*
2 * vp_view.c
3 *
4 * Routines to compute quantities derived from the view transformation.
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.29 $
29 */
30
31 #include "vp_global.h"
32
33 static int FactorAffineView ANSI_ARGS((vpContext *vpc, vpMatrix4 vm));
34 static int FactorPerspectiveView ANSI_ARGS((vpContext *vpc, vpMatrix4 vm));
35 static void ComputeAffineOpacityCorrection ANSI_ARGS((vpContext *vpc,
36 double shear_i, double shear_j, float table[VP_OPACITY_MAX+1]));
37 static void CheckRenderBuffers ANSI_ARGS((vpContext *vpc));
38 static void ComputeLightViewTransform ANSI_ARGS((vpContext *vpc,vpMatrix4 vm));
39 static int FactorLightView ANSI_ARGS((vpContext *vpc, vpMatrix4 vm));
40 static void CheckShadowBuffer ANSI_ARGS((vpContext *vpc));
41
42 /*
43 * VPFactorView
44 *
45 * Factor the viewing matrix.
46 */
47
48 vpResult
VPFactorView(vpc)49 VPFactorView(vpc)
50 vpContext *vpc;
51 {
52 vpMatrix4 vm;
53 int retcode;
54
55 if (vpc->factored_view_ready) {
56 CheckRenderBuffers(vpc);
57 return(VP_OK);
58 }
59
60 /* compute the overall view transformation */
61 VPComputeViewTransform(vpc, vm);
62
63 /* check if transformation is affine and factor it */
64 if (fabs(vm[3][0]) < VP_EPS && fabs(vm[3][1]) < VP_EPS &&
65 fabs(vm[3][2]) < VP_EPS && fabs(vm[3][3]-1.) < VP_EPS) {
66 if ((retcode = FactorAffineView(vpc, vm)) != VP_OK)
67 return(retcode);
68 ComputeAffineOpacityCorrection(vpc, vpc->shear_i, vpc->shear_j,
69 vpc->affine_opac_correct);
70 } else {
71 FactorPerspectiveView(vpc, vm);
72 }
73 CheckRenderBuffers(vpc);
74
75 /* compute viewing transformation from the point of view of the light
76 source (for calculating shadows) */
77 if ((retcode = VPCheckShadows(vpc)) != VP_OK)
78 return(retcode);
79 if (vpc->enable_shadows) {
80 ComputeLightViewTransform(vpc, vm);
81 if ((retcode = FactorLightView(vpc, vm)) != VP_OK)
82 return(retcode);
83 ComputeAffineOpacityCorrection(vpc, vpc->shadow_shear_i,
84 vpc->shadow_shear_j,
85 vpc->shadow_opac_correct);
86 CheckShadowBuffer(vpc);
87 }
88 return(VP_OK);
89 }
90
91 /*
92 * VPComputeViewTransform
93 *
94 * Compute the overall view transformation.
95 */
96
97 void
VPComputeViewTransform(vpc,vm)98 VPComputeViewTransform(vpc, vm)
99 vpContext *vpc;
100 vpMatrix4 vm; /* storage for result */
101 {
102 vpMatrix4 prematrix; /* transform volume to unit cube */
103 vpMatrix4 viewportm; /* viewport matrix */
104 vpMatrix4 tmp1m, tmp2m, tmp3m; /* temporary matrices */
105 int maxdim;
106 double scale;
107
108 /* compute prematrix */
109 vpIdentity4(prematrix);
110 maxdim = vpc->xlen;
111 if (vpc->ylen > maxdim)
112 maxdim = vpc->ylen;
113 if (vpc->zlen > maxdim)
114 maxdim = vpc->zlen;
115 scale = 1. / (double)maxdim;
116 prematrix[0][0] = scale;
117 prematrix[1][1] = scale;
118 prematrix[2][2] = scale;
119 prematrix[0][3] = -vpc->xlen * scale * 0.5;
120 prematrix[1][3] = -vpc->ylen * scale * 0.5;
121 prematrix[2][3] = -vpc->zlen * scale * 0.5;
122
123 /* compute the viewport matrix */
124 vpIdentity4(viewportm);
125 viewportm[0][0] = 0.5 * vpc->image_width;
126 viewportm[0][3] = 0.5 * vpc->image_width;
127 viewportm[1][1] = 0.5 * vpc->image_height;
128 viewportm[1][3] = 0.5 * vpc->image_height;
129 viewportm[2][2] = -0.5; /* minus sign: switch to left-handed coords. */
130 viewportm[2][3] = 0.5;
131
132 /* compute the view matrix */
133 vpMatrixMult4(tmp1m, vpc->transforms[VP_MODEL], prematrix);
134 vpMatrixMult4(tmp2m, vpc->transforms[VP_VIEW], tmp1m);
135 vpMatrixMult4(tmp3m, vpc->transforms[VP_PROJECT], tmp2m);
136 vpMatrixMult4(vm, viewportm, tmp3m);
137 }
138
139 /*
140 * FactorAffineView
141 *
142 * Factor an affine viewing matrix into two parts:
143 * 1) A shear and translation which map object coordinates into
144 * intermediate image coordinates.
145 * 2) An affine warp which transforms intermediate image coordinates to
146 * image coordinates.
147 * Return value is a result code.
148 */
149
150 static int
FactorAffineView(vpc,vm)151 FactorAffineView(vpc, vm)
152 vpContext *vpc;
153 vpMatrix4 vm;
154 {
155 vpMatrix4 p; /* permutation matrix */
156 vpMatrix4 pvm; /* permutation of the viewing matrix */
157 int icount, jcount, kcount; /* dimensions of volume in rotated space */
158 vpVector4 xobj, yobj, zobj;
159 vpVector4 xim, yim, zim;
160 double x, y, z, denom;
161
162 Debug((vpc, VPDEBUG_VIEW, "FactorAffineView\n"));
163
164 vpc->affine_view = 1;
165
166 /*
167 * Transform the unit x, y and z object-coordinate vectors into image
168 * space and see which one is most aligned with the view direction
169 * (which is the z-axis in image coordinates).
170 */
171 vpSetVector4(xobj, 1., 0., 0., 0.);
172 vpSetVector4(yobj, 0., 1., 0., 0.);
173 vpSetVector4(zobj, 0., 0., 1., 0.);
174 vpMatrixVectorMult4(xim, vm, xobj);
175 vpMatrixVectorMult4(yim, vm, yobj);
176 vpMatrixVectorMult4(zim, vm, zobj);
177 x = fabs((vpNormalize3(xim) == VPERROR_SINGULAR) ? 0. : xim[2]);
178 y = fabs((vpNormalize3(yim) == VPERROR_SINGULAR) ? 0. : yim[2]);
179 z = fabs((vpNormalize3(zim) == VPERROR_SINGULAR) ? 0. : zim[2]);
180 if (x >= y) {
181 if (x >= z) {
182 vpc->best_view_axis = VP_X_AXIS;
183 } else {
184 vpc->best_view_axis = VP_Z_AXIS;
185 }
186 } else {
187 if (y >= z) {
188 vpc->best_view_axis = VP_Y_AXIS;
189 } else {
190 vpc->best_view_axis = VP_Z_AXIS;
191 }
192 }
193 switch (vpc->axis_override) {
194 case VP_X_AXIS:
195 if (x < VP_EPS)
196 return(VPSetError(vpc, VPERROR_SINGULAR));
197 vpc->best_view_axis = VP_X_AXIS;
198 break;
199 case VP_Y_AXIS:
200 if (y < VP_EPS)
201 return(VPSetError(vpc, VPERROR_SINGULAR));
202 vpc->best_view_axis = VP_Y_AXIS;
203 break;
204 case VP_Z_AXIS:
205 if (z < VP_EPS)
206 return(VPSetError(vpc, VPERROR_SINGULAR));
207 vpc->best_view_axis = VP_Z_AXIS;
208 break;
209 default:
210 break;
211 }
212
213 /* permute the rows of the viewing matrix so that the third axis is
214 most parallel to the viewing direction */
215 bzero(p, sizeof(vpMatrix4));
216 switch (vpc->best_view_axis) {
217 case VP_X_AXIS:
218 p[0][2] = 1.;
219 p[1][0] = 1.;
220 p[2][1] = 1.;
221 p[3][3] = 1.;
222 icount = vpc->ylen;
223 jcount = vpc->zlen;
224 kcount = vpc->xlen;
225 break;
226 case VP_Y_AXIS:
227 p[0][1] = 1.;
228 p[1][2] = 1.;
229 p[2][0] = 1.;
230 p[3][3] = 1.;
231 icount = vpc->zlen;
232 jcount = vpc->xlen;
233 kcount = vpc->ylen;
234 break;
235 case VP_Z_AXIS:
236 p[0][0] = 1.;
237 p[1][1] = 1.;
238 p[2][2] = 1.;
239 p[3][3] = 1.;
240 icount = vpc->xlen;
241 jcount = vpc->ylen;
242 kcount = vpc->zlen;
243 break;
244 }
245 vpMatrixMult4(pvm, vm, p);
246
247 /* compute the shear coefficients */
248 denom = pvm[0][0]*pvm[1][1] - pvm[0][1]*pvm[1][0];
249 if (fabs(denom) < VP_EPS)
250 return(VPSetError(vpc, VPERROR_SINGULAR));
251 vpc->shear_i = (pvm[0][2]*pvm[1][1] - pvm[0][1]*pvm[1][2]) / denom;
252 vpc->shear_j = (pvm[0][0]*pvm[1][2] - pvm[1][0]*pvm[0][2]) / denom;
253 if (pvm[2][0]*vpc->shear_i + pvm[2][1]*vpc->shear_j - pvm[2][2] > 0)
254 vpc->reverse_slice_order = 0;
255 else
256 vpc->reverse_slice_order = 1;
257
258 /* compute the intermediate image size */
259 vpc->intermediate_width = icount + 1 + (int)ceil((kcount-1)*
260 fabs(vpc->shear_i));
261 vpc->intermediate_height = jcount + 1 + (int)ceil((kcount-1)*
262 fabs(vpc->shear_j));
263
264 /* compute the translation coefficients */
265 if (vpc->shear_i >= 0.)
266 vpc->trans_i = 1.;
267 else
268 vpc->trans_i = 1. - vpc->shear_i * (kcount - 1);
269 if (vpc->shear_j >= 0.)
270 vpc->trans_j = 1.;
271 else
272 vpc->trans_j = 1. - vpc->shear_j * (kcount - 1);
273
274 /* compute the depth coefficients */
275 vpc->depth_di = pvm[2][0];
276 vpc->depth_dj = pvm[2][1];
277 vpc->depth_dk = pvm[2][2];
278 vpc->depth_000 = pvm[2][3];
279
280 /* compute the mapping from compositing space to image space */
281 vpc->warp_2d[0][0] = pvm[0][0];
282 vpc->warp_2d[0][1] = pvm[0][1];
283 vpc->warp_2d[0][2] = pvm[0][3] - pvm[0][0]*vpc->trans_i -
284 pvm[0][1]*vpc->trans_j;
285 vpc->warp_2d[1][0] = pvm[1][0];
286 vpc->warp_2d[1][1] = pvm[1][1];
287 vpc->warp_2d[1][2] = pvm[1][3] - pvm[1][0]*vpc->trans_i -
288 pvm[1][1]*vpc->trans_j;
289 vpc->warp_2d[2][0] = 0.;
290 vpc->warp_2d[2][1] = 0.;
291 vpc->warp_2d[2][2] = 1.;
292
293 vpc->factored_view_ready = 1;
294
295 Debug((vpc, VPDEBUG_VIEW, " best_view_axis: %c%c\n",
296 vpc->reverse_slice_order ? '-' : '+',
297 vpc->best_view_axis == VP_X_AXIS ? 'x' :
298 (vpc->best_view_axis == VP_Y_AXIS ? 'y' : 'z')));
299 Debug((vpc, VPDEBUG_VIEW, " shear factors: %g %g\n",
300 vpc->shear_i, vpc->shear_j));
301 Debug((vpc, VPDEBUG_VIEW, " translation: %g %g\n",
302 vpc->trans_i, vpc->trans_j));
303 Debug((vpc, VPDEBUG_VIEW, " depth: d000: %g\n", vpc->depth_000));
304 Debug((vpc, VPDEBUG_VIEW, " di: %g\n", vpc->depth_di));
305 Debug((vpc, VPDEBUG_VIEW, " dj: %g\n", vpc->depth_dj));
306 Debug((vpc, VPDEBUG_VIEW, " dk: %g\n", vpc->depth_dk));
307 Debug((vpc, VPDEBUG_VIEW, " intermediate image size: %d %d\n",
308 vpc->intermediate_width, vpc->intermediate_height));
309 return(VP_OK);
310 }
311
312 /*
313 * FactorPerspectiveView
314 *
315 * Factor a perspective view matrix into two parts:
316 * 1) A shear, translation and scale which map object coordinates into
317 * intermediate image coordinates.
318 * 2) A perspective warp which transforms intermediate image coordinates to
319 * image coordinates.
320 */
321
322 static int
FactorPerspectiveView(vpc,vm)323 FactorPerspectiveView(vpc, vm)
324 vpContext *vpc;
325 vpMatrix4 vm;
326 {
327 vpc->affine_view = 0;
328 return(VP_OK);
329
330 #ifdef notdef
331 Matrix4 p; /* permutation matrix */
332 Matrix4 pvm; /* permutation of the viewing matrix */
333 Matrix4 m2d; /* final warp */
334 Matrix4 t;
335 double alpha1, alpha2, alpha3, alpha4;
336 int icount, jcount, kcount; /* dimensions of volume in rotated space */
337 Vector4 xobj, yobj, zobj;
338 double x, y, z, denom;
339 double i0, j0, i1, j1;
340 double imin, imax, jmin, jmax;
341
342 Debug((DEBUG_VIEW, "FactorPerspectiveView\n"));
343
344 rbuf->perspective_proj = 1;
345
346 /*
347 * Transform the unit x, y and z object-coordinate vectors into image
348 * space and see which one is most aligned with the view direction
349 * (which is the z-axis in image coordinates).
350 */
351 xobj[0] = 1.; xobj[1] = 0.; xobj[2] = 0.; xobj[3] = 0.;
352 yobj[0] = 0.; yobj[1] = 1.; yobj[2] = 0.; yobj[3] = 0.;
353 zobj[0] = 0.; zobj[1] = 0.; zobj[2] = 1.; zobj[3] = 0.;
354 TransformVector4(xobj, view->view_matrix);
355 TransformVector4(yobj, view->view_matrix);
356 TransformVector4(zobj, view->view_matrix);
357 /* normalize each vector to unit length and compare the absolute value
358 of the z component; note that the w component drops out */
359 xobj[2] = fabs(xobj[2]);
360 yobj[2] = fabs(yobj[2]);
361 zobj[2] = fabs(zobj[2]);
362 x = (xobj[2] < EPS) ? 0. :
363 (xobj[2] / sqrt(xobj[0]*xobj[0] + xobj[1]*xobj[1] + xobj[2]*xobj[2]));
364 y = (yobj[2] < EPS) ? 0. :
365 (yobj[2] / sqrt(yobj[0]*yobj[0] + yobj[1]*yobj[1] + yobj[2]*yobj[2]));
366 z = (zobj[2] < EPS) ? 0. :
367 (zobj[2] / sqrt(zobj[0]*zobj[0] + zobj[1]*zobj[1] + zobj[2]*zobj[2]));
368 if (x >= y) {
369 if (x >= z) {
370 rbuf->best_view_axis = VP_XAXIS;
371 } else {
372 rbuf->best_view_axis = VP_ZAXIS;
373 }
374 } else {
375 if (y >= z) {
376 rbuf->best_view_axis = VP_YAXIS;
377 } else {
378 rbuf->best_view_axis = VP_ZAXIS;
379 }
380 }
381
382 /* permute the rows of the viewing matrix so that the third axis is
383 most parallel to the viewing direction */
384 bzero(p, sizeof(Matrix4));
385 switch (rbuf->best_view_axis) {
386 case VP_XAXIS:
387 p[0][2] = 1.;
388 p[1][0] = 1.;
389 p[2][1] = 1.;
390 p[3][3] = 1.;
391 icount = ylen;
392 jcount = zlen;
393 kcount = xlen;
394 break;
395 case VP_YAXIS:
396 p[0][1] = 1.;
397 p[1][2] = 1.;
398 p[2][0] = 1.;
399 p[3][3] = 1.;
400 icount = zlen;
401 jcount = xlen;
402 kcount = ylen;
403 break;
404 case VP_ZAXIS:
405 p[0][0] = 1.;
406 p[1][1] = 1.;
407 p[2][2] = 1.;
408 p[3][3] = 1.;
409 icount = xlen;
410 jcount = ylen;
411 kcount = zlen;
412 break;
413 default:
414 VPBug("wierd value for best_view_axis in FactorPerspectiveView\n");
415 }
416 MatrixMult4(pvm, view->view_matrix, p);
417
418 /* compute the magic alpha coefficients */
419 alpha1 = pvm[3][1] * (pvm[0][3]*pvm[1][2] - pvm[0][2]*pvm[1][3]) +
420 pvm[3][2] * (pvm[0][1]*pvm[1][3] - pvm[0][3]*pvm[1][1]) +
421 pvm[3][3] * (pvm[0][2]*pvm[1][1] - pvm[0][1]*pvm[1][2]);
422 alpha2 = pvm[3][0] * (pvm[0][2]*pvm[1][3] - pvm[0][3]*pvm[1][2]) +
423 pvm[3][2] * (pvm[0][3]*pvm[1][0] - pvm[0][0]*pvm[1][3]) +
424 pvm[3][3] * (pvm[0][0]*pvm[1][2] - pvm[0][2]*pvm[1][0]);
425 alpha3 = pvm[3][0] * (pvm[0][1]*pvm[1][2] - pvm[0][2]*pvm[1][1]) +
426 pvm[3][1] * (pvm[0][2]*pvm[1][0] - pvm[0][0]*pvm[1][2]) +
427 pvm[3][2] * (pvm[0][0]*pvm[1][1] - pvm[0][1]*pvm[1][0]);
428 alpha4 = pvm[3][0] * (pvm[0][1]*pvm[1][3] - pvm[0][3]*pvm[1][1]) +
429 pvm[3][1] * (pvm[0][3]*pvm[1][0] - pvm[0][0]*pvm[1][3]) +
430 pvm[3][3] * (pvm[0][0]*pvm[1][1] - pvm[0][1]*pvm[1][0]);
431
432 /* determine the order of the slices */
433 if (pvm[2][2] - (pvm[2][0]*alpha1 + pvm[2][1]*alpha2)/(alpha3+alpha4) > 0)
434 rbuf->reverse_k_order = 1;
435 else
436 rbuf->reverse_k_order = 0;
437
438 /* compute the scale coefficients */
439 rbuf->w_factor = alpha3 / alpha4;
440 if (rbuf->reverse_k_order)
441 rbuf->normalize_scale = 1. + rbuf->w_factor*(kcount-1);
442 else
443 rbuf->normalize_scale = 1.;
444
445 /* compute the bounding box of the image in compositing space */
446 denom = 1. / (alpha4 + alpha3*(kcount-1));
447 i0 = rbuf->normalize_scale*alpha1*(kcount-1) * denom;
448 j0 = rbuf->normalize_scale*alpha2*(kcount-1) * denom;
449 i1 = rbuf->normalize_scale*(alpha4*icount + alpha1*(kcount-1)) * denom;
450 j1 = rbuf->normalize_scale*(alpha4*jcount + alpha2*(kcount-1)) * denom;
451 imin = MIN(0, i0);
452 imax = MAX(rbuf->normalize_scale*icount, i1);
453 jmin = MIN(0, j0);
454 jmax = MAX(rbuf->normalize_scale*jcount, j1);
455
456 /* compute the size of the intermediate image */
457 rbuf->intermediate_width = (int)ceil(imax - imin);
458 rbuf->intermediate_height = (int)ceil(jmax - jmin);
459
460 /* compute the translation and shear coefficients */
461 rbuf->shear_i = (rbuf->normalize_scale*alpha1 - alpha3*imin) / alpha4;
462 rbuf->shear_j = (rbuf->normalize_scale*alpha2 - alpha3*jmin) / alpha4;
463 rbuf->trans_i = -imin;
464 rbuf->trans_j = -jmin;
465
466 /* compute the depth coefficients */
467 rbuf->depth_di = pvm[2][0];
468 rbuf->depth_dj = pvm[2][1];
469 rbuf->depth_dk = pvm[2][2];
470 rbuf->depth_000 = pvm[2][3];
471 rbuf->w_di = pvm[3][0];
472 rbuf->w_dj = pvm[3][1];
473 rbuf->w_dk = pvm[3][2];
474 rbuf->w_000 = pvm[3][3];
475
476 /* compute the mapping from compositing space to image space */
477 Identity4(t);
478 t[0][0] = 1. / rbuf->normalize_scale;
479 t[1][1] = 1. / rbuf->normalize_scale;
480 t[0][2] = -alpha1 / alpha4;
481 t[1][2] = -alpha2 / alpha4;
482 t[3][2] = -alpha3 / alpha4;
483 t[0][3] = imin / rbuf->normalize_scale;
484 t[1][3] = jmin / rbuf->normalize_scale;
485 MatrixMult4(m2d, pvm, t);
486
487 rbuf->warp_2d[0][0] = m2d[0][0];
488 rbuf->warp_2d[1][0] = m2d[1][0];
489 rbuf->warp_2d[2][0] = m2d[3][0];
490 rbuf->warp_2d[0][1] = m2d[0][1];
491 rbuf->warp_2d[1][1] = m2d[1][1];
492 rbuf->warp_2d[2][1] = m2d[3][1];
493 rbuf->warp_2d[0][2] = m2d[0][3];
494 rbuf->warp_2d[1][2] = m2d[1][3];
495 rbuf->warp_2d[2][2] = m2d[3][3];
496 #endif /* notdef */
497 }
498
499 /*
500 * ComputeAffineOpacityCorrection
501 *
502 * Precompute a lookup table which corrects opacity for an affine viewing
503 * transformation. (Opacity correction accounts for variations in the
504 * apparent thickness of a voxel depending on viewpoint.)
505 */
506
507 static void
ComputeAffineOpacityCorrection(vpc,shear_i,shear_j,table)508 ComputeAffineOpacityCorrection(vpc, shear_i, shear_j, table)
509 vpContext *vpc;
510 double shear_i;
511 double shear_j;
512 float table[VP_OPACITY_MAX+1];
513 {
514 float voxel_size;
515 int i;
516
517 Debug((vpc, VPDEBUG_OPCCORRECT,
518 "Computing affine opacity correction table.\n"));
519 voxel_size = sqrt(1 + shear_i*shear_i + shear_j*shear_j);
520 for (i = 0; i <= VP_OPACITY_MAX; i++) {
521 #ifdef NO_OPAC_CORRECT
522 table[i] = (double)i / (double)VP_OPACITY_MAX;
523 #else
524 table[i] = 1.-pow(1.-(double)i/(double)VP_OPACITY_MAX,voxel_size);
525 #endif
526 }
527 }
528
529 /*
530 * CheckRenderBuffers
531 *
532 * Resize the buffers used during rendering, if necessary.
533 */
534
535 static void
CheckRenderBuffers(vpc)536 CheckRenderBuffers(vpc)
537 vpContext *vpc;
538 {
539 int new_max_width, new_max_height, new_max_scan;
540 int resize = 0;
541
542 /* determine if resizing is necessary */
543 if (vpc->intermediate_width > vpc->max_intermediate_width) {
544 new_max_width = MAX(vpc->intermediate_width,
545 vpc->int_image_width_hint);
546 resize = 1;
547 } else {
548 new_max_width = MAX(vpc->max_intermediate_width,
549 vpc->int_image_width_hint);
550 }
551 if (vpc->intermediate_height > vpc->max_intermediate_height) {
552 new_max_height = MAX(vpc->intermediate_height,
553 vpc->int_image_height_hint);
554 resize = 1;
555 } else {
556 new_max_height = MAX(vpc->max_intermediate_height,
557 vpc->int_image_height_hint);
558 }
559 new_max_scan = vpc->xlen;
560 if (vpc->ylen > new_max_scan)
561 new_max_scan = vpc->ylen;
562 if (vpc->zlen > new_max_scan)
563 new_max_scan = vpc->zlen;
564 if (new_max_scan > vpc->max_scan_length)
565 resize = 1;
566 if (vpc->color_channels != vpc->intermediate_color_channels)
567 resize = 1;
568
569 /* resize */
570 if (resize)
571 VPResizeRenderBuffers(vpc, new_max_width, new_max_height,new_max_scan);
572 }
573
574 /*
575 * VPResizeRenderBuffers
576 *
577 * Resize the rendering buffers.
578 */
579
580 void
VPResizeRenderBuffers(vpc,max_width,max_height,max_scan)581 VPResizeRenderBuffers(vpc, max_width, max_height, max_scan)
582 vpContext *vpc;
583 int max_width; /* new width of the intermediate image */
584 int max_height; /* new height of the intermediate image */
585 int max_scan; /* new max. scanline length */
586 {
587 /* free old buffers */
588 if (vpc->int_image.gray_intim != NULL) {
589 Debug((vpc, VPDEBUG_RBUF, "Freeing old RenderBuffer(%d,%d,%d,%d)\n",
590 vpc->max_intermediate_width, vpc->max_intermediate_height,
591 vpc->max_scan_length, vpc->intermediate_color_channels));
592 Dealloc(vpc, vpc->int_image.gray_intim);
593 }
594
595 /* allocate new buffers */
596 Debug((vpc, VPDEBUG_RBUF, "Allocating RenderBuffer(%d,%d,%d,%d)\n",
597 max_width, max_height, max_scan, vpc->color_channels));
598 vpc->max_intermediate_width = max_width;
599 vpc->max_intermediate_height = max_height;
600 vpc->max_scan_length = max_scan;
601 vpc->intermediate_color_channels = vpc->color_channels;
602 if (max_width > 0) {
603 if (vpc->color_channels == 1) {
604 Alloc(vpc, vpc->int_image.gray_intim, GrayIntPixel *,
605 max_width * max_height * sizeof(GrayIntPixel), "int_image");
606 } else {
607 Alloc(vpc, vpc->int_image.rgb_intim, RGBIntPixel *,
608 max_width * max_height * sizeof(RGBIntPixel), "int_image");
609 }
610 } else {
611 vpc->int_image.gray_intim = NULL;
612 }
613 }
614
615 /*
616 * VPResizeDepthCueTable
617 *
618 * Resize the depth cueing table.
619 */
620
621 void
VPResizeDepthCueTable(vpc,entries,copy)622 VPResizeDepthCueTable(vpc, entries, copy)
623 vpContext *vpc;
624 int entries; /* new number of table entries */
625 int copy; /* if true, copy old entries */
626 {
627 float *new_dc_table;
628
629 Debug((vpc, VPDEBUG_DEPTHCUE, "resizing dctable to %d entries (%s)\n",
630 entries, copy ? "copy" : "nocopy"));
631 if (entries == 0) {
632 if (vpc->dc_table != NULL) {
633 Dealloc(vpc, vpc->dc_table);
634 vpc->dc_table = NULL;
635 }
636 vpc->dc_table_len = 0;
637 } else {
638 Alloc(vpc, new_dc_table, float *, entries * sizeof(float), "dc_table");
639 if (vpc->dc_table != NULL) {
640 if (copy && vpc->dc_table_len > 0) {
641 bcopy(vpc->dc_table, new_dc_table,
642 MIN(vpc->dc_table_len, entries) * sizeof(float));
643 }
644 Dealloc(vpc, vpc->dc_table);
645 }
646 vpc->dc_table = new_dc_table;
647 vpc->dc_table_len = entries;
648 }
649 }
650
651 /*
652 * VPComputeDepthCueTable
653 *
654 * Compute entries in the depth cueing lookup table.
655 */
656
657 void
VPComputeDepthCueTable(vpc,first,last)658 VPComputeDepthCueTable(vpc, first, last)
659 vpContext *vpc;
660 int first; /* first entry to compute */
661 int last; /* last entry to compute */
662 {
663 int c;
664 double delta_depth, front_factor, density;
665
666 Debug((vpc, VPDEBUG_DEPTHCUE, "computing dctable entries %d to %d\n",
667 first, last));
668 delta_depth = vpc->dc_quantization;
669 front_factor = vpc->dc_front_factor;
670 density = vpc->dc_density;
671 for (c = first; c <= last; c++)
672 vpc->dc_table[c] = front_factor * exp(-density*(1.0 - c*delta_depth));
673 }
674
675 /*
676 * VPSliceDepthCueRatio
677 *
678 * Return the ratio of the depth cueing factor for two adjacent slices
679 * for an affine view. A constant factor is applied to all voxels in a
680 * slice, and then a fixup is applied to the pixels of the intermediate
681 * image. This produces the correct answer without having to compute
682 * the depth of each voxel.
683 */
684
685 float
VPSliceDepthCueRatio(vpc)686 VPSliceDepthCueRatio(vpc)
687 vpContext *vpc;
688 {
689 float delta_depth; /* change in depth between adjacent slices */
690 float slice_dc_ratio; /* return value */
691
692 if (!vpc->dc_enable)
693 return(1.);
694 delta_depth = vpc->depth_dk - vpc->depth_di*vpc->shear_i -
695 vpc->depth_dj*vpc->shear_j;
696 if (vpc->reverse_slice_order)
697 delta_depth = -delta_depth;
698 /* slice_dc_ratio = exp(-vpc->dc_density * (-delta_depth)) */
699 slice_dc_ratio = exp(vpc->dc_density * delta_depth);
700 Debug((vpc, VPDEBUG_DEPTHCUE, "slice_dc_ratio = %f\n", slice_dc_ratio));
701 return(slice_dc_ratio);
702 }
703
704 /*
705 * VPDepthCueIntImage
706 *
707 * Perform depth cueing on the intermediate image.
708 */
709
710 void
VPDepthCueIntImage(vpc,slicenum)711 VPDepthCueIntImage(vpc, slicenum)
712 vpContext *vpc;
713 int slicenum; /* slice number corresponding to location of int. image */
714 {
715 float pixel_depth_quant; /* depth of current pixel in image
716 (multiplied by depth_quant) */
717 int pixel_depth_int; /* pixel_depth truncated to an integer */
718 float left_depth; /* depth of pixel on left edge of current
719 scanline in image */
720 float left_depth_quant; /* left_depth * depth_quant */
721 float *dc_table; /* depth cueing table */
722 float depth_di, depth_dj; /* change in depth for a unit change in each
723 rotated object space coordinate */
724 float depth_quant; /* number of quantization levels for depth */
725 float depth_di_quant; /* depth_di * depth_quant */
726 float depth_dj_quant; /* depth_dj * depth_quant */
727 float max_depth; /* maximum (closest) depth in image */
728 int max_depth_int; /* maximum quantized depth */
729 int i, j; /* intermediate image coordinates */
730 float slice_u, slice_v; /* sheared object space coordinates */
731 GrayIntPixel *gray_intim; /* image data (grayscale) */
732 RGBIntPixel *rgb_intim; /* image data (RGB) */
733 int width, height; /* size of intermediate image */
734 int c;
735 #ifdef DEBUG
736 float pix_depth;
737 #endif
738
739 Debug((vpc, VPDEBUG_DEPTHCUE, "depth cueing intermediate image\n"));
740
741 /* check the size of the depth cueing table and enlarge if necessary */
742 width = vpc->intermediate_width;
743 height = vpc->intermediate_height;
744 depth_quant = 1.0 / vpc->dc_quantization;
745 depth_di = vpc->depth_di;
746 depth_dj = vpc->depth_dj;
747 slice_u = vpc->shear_i * slicenum + vpc->trans_i;
748 slice_v = vpc->shear_j * slicenum + vpc->trans_j;
749 left_depth = vpc->depth_000 + vpc->depth_dk*slicenum -
750 slice_u*depth_di - slice_v*depth_dj;
751 if (depth_di > 0) {
752 if (depth_dj > 0) {
753 max_depth = left_depth + depth_di*width + depth_dj*height;
754 } else {
755 max_depth = left_depth + depth_di*width;
756 }
757 } else {
758 if (depth_dj > 0) {
759 max_depth = left_depth + depth_dj*height;
760 } else {
761 max_depth = left_depth;
762 }
763 }
764 max_depth_int = max_depth * depth_quant;
765 if (max_depth_int >= vpc->dc_table_len) {
766 c = vpc->dc_table_len;
767 VPResizeDepthCueTable(vpc, max_depth_int+1, 1);
768 VPComputeDepthCueTable(vpc, c, vpc->dc_table_len-1);
769 }
770 dc_table = vpc->dc_table;
771 depth_di_quant = depth_di * depth_quant;
772 depth_dj_quant = depth_dj * depth_quant;
773 left_depth_quant = left_depth * depth_quant;
774
775 #ifdef DEBUG
776 Debug((vpc, VPDEBUG_DEPTHCUE, "depth cueing at image corners:\n"));
777 pix_depth = left_depth + 0*depth_di + 0*depth_dj;
778 pixel_depth_int = (int)(pix_depth * depth_quant);
779 if (pixel_depth_int < 0) pixel_depth_int = 0;
780 Debug((vpc, VPDEBUG_DEPTHCUE,
781 " %3d %3d: depth = %10.6f, factor = %10.6f, table[%d] = %10.6f\n",
782 0, 0, pix_depth,
783 vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - pix_depth)),
784 pixel_depth_int, dc_table[pixel_depth_int]));
785 pix_depth = left_depth + width*depth_di + 0*depth_dj;
786 pixel_depth_int = (int)(pix_depth * depth_quant);
787 if (pixel_depth_int < 0) pixel_depth_int = 0;
788 Debug((vpc, VPDEBUG_DEPTHCUE,
789 " %3d %3d: depth = %10.6f, factor = %10.6f, table[%d] = %10.6f\n",
790 width, 0, pix_depth,
791 vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - pix_depth)),
792 pixel_depth_int, dc_table[pixel_depth_int]));
793 pix_depth = left_depth + width*depth_di + height*depth_dj;
794 pixel_depth_int = (int)(pix_depth * depth_quant);
795 if (pixel_depth_int < 0) pixel_depth_int = 0;
796 Debug((vpc, VPDEBUG_DEPTHCUE,
797 " %3d %3d: depth = %10.6f, factor = %10.6f, table[%d] = %10.6f\n",
798 width, height, pix_depth,
799 vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - pix_depth)),
800 pixel_depth_int, dc_table[pixel_depth_int]));
801 pix_depth = left_depth + 0*depth_di + height*depth_dj;
802 pixel_depth_int = (int)(pix_depth * depth_quant);
803 if (pixel_depth_int < 0) pixel_depth_int = 0;
804 Debug((vpc, VPDEBUG_DEPTHCUE,
805 " %3d %3d: depth = %10.6f, factor = %10.6f, table[%d] = %10.6f\n",
806 0, height, pix_depth,
807 vpc->dc_front_factor * exp(-vpc->dc_density * (1.0 - pix_depth)),
808 pixel_depth_int, dc_table[pixel_depth_int]));
809 #endif /* DEBUG */
810
811 /* foreach pixel, compute depth and scale color by dc factor */
812 if (vpc->color_channels == 1) {
813 gray_intim = vpc->int_image.gray_intim;
814 for (j = height; j > 0; j--) {
815 pixel_depth_quant = left_depth_quant;
816 left_depth_quant += depth_dj_quant;
817 for (i = width; i > 0; i--) {
818 pixel_depth_int = pixel_depth_quant;
819 pixel_depth_quant += depth_di_quant;
820 if (pixel_depth_int < 0)
821 pixel_depth_int = 0;
822 if (pixel_depth_int >= vpc->dc_table_len) {
823 VPBug("VPDepthCueIntImage: depth too large (%d >= %d)",
824 pixel_depth_int, vpc->dc_table_len);
825 }
826 gray_intim->clrflt *= dc_table[pixel_depth_int];
827 gray_intim++;
828 } /* for i */
829 } /* for j */
830 } else {
831 rgb_intim = vpc->int_image.rgb_intim;
832 for (j = height; j > 0; j--) {
833 pixel_depth_quant = left_depth_quant;
834 left_depth_quant += depth_dj_quant;
835 for (i = width; i > 0; i--) {
836 pixel_depth_int = pixel_depth_quant;
837 pixel_depth_quant += depth_di_quant;
838 if (pixel_depth_int < 0)
839 pixel_depth_int = 0;
840 if (pixel_depth_int >= vpc->dc_table_len) {
841 VPBug("VPDepthCueIntImage: depth too large (%d >= %d)",
842 pixel_depth_int, vpc->dc_table_len);
843 }
844 rgb_intim->rclrflt *= dc_table[pixel_depth_int];
845 rgb_intim->gclrflt *= dc_table[pixel_depth_int];
846 rgb_intim->bclrflt *= dc_table[pixel_depth_int];
847 rgb_intim++;
848 } /* for i */
849 } /* for j */
850 }
851 }
852
853 /*
854 * ComputeLightViewTransform
855 *
856 * Compute the view transformation from the point of view of the one
857 * light source that produces shadows.
858 */
859
860 static void
ComputeLightViewTransform(vpc,vm)861 ComputeLightViewTransform(vpc, vm)
862 vpContext *vpc;
863 vpMatrix4 vm; /* storage for result */
864 {
865 vpMatrix4 prematrix; /* transform volume to unit cube */
866 vpMatrix4 viewportm; /* viewport matrix */
867 vpVector3 vpn; /* view plane normal */
868 vpVector3 vup; /* view up vector */
869 vpVector3 tmp1v, tmp2v; /* temporary vectors */
870 vpMatrix4 view; /* transform world coordinates to
871 eye coordinates, with view
872 direction equal to light vector */
873 vpMatrix4 tmp1m, tmp2m; /* temporary matrices */
874 double lx, ly, lz; /* components of light vector */
875 int light_index;
876 int maxdim;
877 double scale;
878
879 /* check for errors */
880 ASSERT(vpc->shadow_light_num >= VP_LIGHT0 &&
881 vpc->shadow_light_num <= VP_LIGHT5);
882 ASSERT(vpc->light_enable[vpc->shadow_light_num - VP_LIGHT0]);
883
884 /* compute prematrix */
885 vpIdentity4(prematrix);
886 maxdim = vpc->xlen;
887 if (vpc->ylen > maxdim)
888 maxdim = vpc->ylen;
889 if (vpc->zlen > maxdim)
890 maxdim = vpc->zlen;
891 scale = 1. / (double)maxdim;
892 prematrix[0][0] = scale;
893 prematrix[1][1] = scale;
894 prematrix[2][2] = scale;
895 prematrix[0][3] = -vpc->xlen * scale * 0.5;
896 prematrix[1][3] = -vpc->ylen * scale * 0.5;
897 prematrix[2][3] = -vpc->zlen * scale * 0.5;
898
899 /* compute the world-to-eye coordinate transformation */
900 light_index = vpc->shadow_light_num - VP_LIGHT0;
901 lx = vpc->light_vector[light_index][0];
902 ly = vpc->light_vector[light_index][1];
903 lz = vpc->light_vector[light_index][2];
904 vpSetVector3(vpn, lx, ly, lz);
905 if (fabs(lx) < fabs(ly)) {
906 if (fabs(lx) < fabs(lz)) {
907 vpSetVector3(vup, 1.0, 0.0, 0.0);
908 } else {
909 vpSetVector3(vup, 0.0, 0.0, 1.0);
910 }
911 } else {
912 if (fabs(ly) < fabs(lz)) {
913 vpSetVector3(vup, 0.0, 1.0, 0.0);
914 } else {
915 vpSetVector3(vup, 0.0, 0.0, 1.0);
916 }
917 }
918 vpCrossProduct(tmp1v, vup, vpn);
919 vpCrossProduct(tmp2v, vpn, tmp1v);
920
921 vpIdentity4(view);
922 view[0][0] = tmp1v[0];
923 view[0][1] = tmp1v[1];
924 view[0][2] = tmp1v[2];
925 view[1][0] = tmp2v[0];
926 view[1][1] = tmp2v[1];
927 view[1][2] = tmp2v[2];
928 view[2][0] = vpn[0];
929 view[2][1] = vpn[1];
930 view[2][2] = vpn[2];
931
932 /* initialize matrix to switch to left-handed coords. */
933 vpIdentity4(viewportm);
934 viewportm[2][2] = -1.0;
935
936 /* compute the view matrix */
937 vpMatrixMult4(tmp1m, vpc->transforms[VP_MODEL], prematrix);
938 vpMatrixMult4(tmp2m, view, tmp1m);
939 vpMatrixMult4(vm, viewportm, tmp2m);
940 }
941
942 /*
943 * FactorLightView
944 *
945 * Factor an affine viewing matrix that specifies the view seen by
946 * a light source. Most of the parameters of the factorization are
947 * taken from the factorization of the normal viewing matrix.
948 */
949
950 static int
FactorLightView(vpc,vm)951 FactorLightView(vpc, vm)
952 vpContext *vpc;
953 vpMatrix4 vm;
954 {
955 vpMatrix4 p; /* permutation matrix */
956 vpMatrix4 pvm; /* permutation of the viewing matrix */
957 int icount, jcount, kcount; /* dimensions of volume in rotated space */
958 double denom;
959
960 Debug((vpc, VPDEBUG_SHADOW, "FactorLightView\n"));
961
962 /* permute the rows of the viewing matrix according to the best viewing
963 axis for the viewing direction */
964 bzero(p, sizeof(vpMatrix4));
965 switch (vpc->best_view_axis) {
966 case VP_X_AXIS:
967 p[0][2] = 1.;
968 p[1][0] = 1.;
969 p[2][1] = 1.;
970 p[3][3] = 1.;
971 icount = vpc->ylen;
972 jcount = vpc->zlen;
973 kcount = vpc->xlen;
974 break;
975 case VP_Y_AXIS:
976 p[0][1] = 1.;
977 p[1][2] = 1.;
978 p[2][0] = 1.;
979 p[3][3] = 1.;
980 icount = vpc->zlen;
981 jcount = vpc->xlen;
982 kcount = vpc->ylen;
983 break;
984 case VP_Z_AXIS:
985 p[0][0] = 1.;
986 p[1][1] = 1.;
987 p[2][2] = 1.;
988 p[3][3] = 1.;
989 icount = vpc->xlen;
990 jcount = vpc->ylen;
991 kcount = vpc->zlen;
992 break;
993 }
994 vpMatrixMult4(pvm, vm, p);
995
996 /* compute the shear coefficients */
997 denom = pvm[0][0]*pvm[1][1] - pvm[0][1]*pvm[1][0];
998 if (fabs(denom) < VP_EPS)
999 return(VPSetError(vpc, VPERROR_SINGULAR));
1000 vpc->shadow_shear_i = (pvm[0][2]*pvm[1][1] - pvm[0][1]*pvm[1][2]) / denom;
1001 vpc->shadow_shear_j = (pvm[0][0]*pvm[1][2] - pvm[1][0]*pvm[0][2]) / denom;
1002
1003 /* check that light direction is compatible with compositing direction */
1004 if (pvm[2][0]*vpc->shadow_shear_i + pvm[2][1]*vpc->shadow_shear_j -
1005 pvm[2][2] > 0) {
1006 if (vpc->reverse_slice_order != 0)
1007 return(VPERROR_BAD_SHADOW);
1008 } else {
1009 if (vpc->reverse_slice_order != 1)
1010 return(VPERROR_BAD_SHADOW);
1011 }
1012
1013 /* compute the shadow buffer image size */
1014 vpc->shadow_width = (int)ceil((kcount-1)*fabs(vpc->shadow_shear_i)) +
1015 icount + 1;
1016 vpc->shadow_height = (int)ceil((kcount-1)*fabs(vpc->shadow_shear_j)) +
1017 jcount + 1;
1018
1019 /* compute the translation coefficients */
1020 if (vpc->shadow_shear_i >= 0.)
1021 vpc->shadow_trans_i = 1.;
1022 else
1023 vpc->shadow_trans_i = 1. - vpc->shadow_shear_i * (kcount - 1);
1024 if (vpc->shadow_shear_j >= 0.)
1025 vpc->shadow_trans_j = 1.;
1026 else
1027 vpc->shadow_trans_j = 1. - vpc->shadow_shear_j * (kcount - 1);
1028
1029 Debug((vpc, VPDEBUG_SHADOW, " shadow shear factors: %g %g\n",
1030 vpc->shadow_shear_i, vpc->shadow_shear_j));
1031 Debug((vpc, VPDEBUG_SHADOW, " shadow translation: %g %g\n",
1032 vpc->shadow_trans_i, vpc->shadow_trans_j));
1033
1034 return(VP_OK);
1035 }
1036
1037 /*
1038 * CheckShadowBuffer
1039 *
1040 * Resize the shadow buffer, if necessary.
1041 */
1042
1043 static void
CheckShadowBuffer(vpc)1044 CheckShadowBuffer(vpc)
1045 vpContext *vpc;
1046 {
1047 int new_max_width, new_max_height;
1048 int resize = 0;
1049
1050 /* determine if resizing is necessary */
1051 if (vpc->shadow_width > vpc->max_shadow_width) {
1052 new_max_width = MAX(vpc->shadow_width, vpc->shadow_width_hint);
1053 resize = 1;
1054 } else {
1055 new_max_width = MAX(vpc->max_shadow_width, vpc->shadow_width_hint);
1056 }
1057 if (vpc->shadow_height > vpc->max_shadow_height) {
1058 new_max_height = MAX(vpc->shadow_height, vpc->shadow_height_hint);
1059 resize = 1;
1060 } else {
1061 new_max_height = MAX(vpc->max_shadow_height, vpc->shadow_height_hint);
1062 }
1063
1064 /* resize */
1065 if (resize)
1066 VPResizeShadowBuffer(vpc, new_max_width, new_max_height);
1067 }
1068
1069 /*
1070 * VPResizeShadowBuffer
1071 *
1072 * Resize the shadow buffer.
1073 */
1074
1075 void
VPResizeShadowBuffer(vpc,max_width,max_height)1076 VPResizeShadowBuffer(vpc, max_width, max_height)
1077 vpContext *vpc;
1078 int max_width; /* new width of the intermediate image */
1079 int max_height; /* new height of the intermediate image */
1080 {
1081 /* free old buffers */
1082 if (vpc->shadow_buffer != NULL) {
1083 Dealloc(vpc, vpc->shadow_buffer);
1084 }
1085
1086 /* allocate new buffers */
1087 vpc->max_shadow_width = max_width;
1088 vpc->max_shadow_height = max_height;
1089 if (max_width > 0) {
1090 Alloc(vpc, vpc->shadow_buffer, GrayIntPixel *,
1091 max_width * max_height * sizeof(GrayIntPixel), "shadow_buffer");
1092 } else {
1093 vpc->shadow_buffer = NULL;
1094 }
1095 }
1096