1 /****************************************************************************
2 * radiosit.cpp
3 *
4 * This module contains all radiosity calculation functions.
5 *
6 * This file was written by Jim McElhiney.
7 *
8 * from Persistence of Vision(tm) Ray Tracer version 3.6.
9 * Copyright 1991-2003 Persistence of Vision Team
10 * Copyright 2003-2004 Persistence of Vision Raytracer Pty. Ltd.
11 *---------------------------------------------------------------------------
12 * NOTICE: This source code file is provided so that users may experiment
13 * with enhancements to POV-Ray and to port the software to platforms other
14 * than those supported by the POV-Ray developers. There are strict rules
15 * regarding how you are permitted to use this file. These rules are contained
16 * in the distribution and derivative versions licenses which should have been
17 * provided with this file.
18 *
19 * These licences may be found online, linked from the end-user license
20 * agreement that is located at http://www.povray.org/povlegal.html
21 *---------------------------------------------------------------------------
22 * This program is based on the popular DKB raytracer version 2.12.
23 * DKBTrace was originally written by David K. Buck.
24 * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
25 *---------------------------------------------------------------------------
26 * $File: //depot/povray/3.6-release/source/radiosit.cpp $
27 * $Revision: #3 $
28 * $Change: 3032 $
29 * $DateTime: 2004/08/02 18:43:41 $
30 * $Author: chrisc $
31 * $Log$
32 *****************************************************************************/
33
34 /************************************************************************
35 * Radiosity calculation routies.
36 *
37 * (This does not work the way that most radiosity programs do, but it accomplishes
38 * the diffuse interreflection integral the hard way and produces similar results. It
39 * is called radiosity here to avoid confusion with ambient and diffuse, which
40 * already have well established meanings within POV).
41 * Inspired by the paper "A Ray Tracing Solution for Diffuse Interreflection"
42 * by Ward, Rubinstein, and Clear, in Siggraph '88 proceedings.
43 *
44 * Basic Idea: Never use a constant ambient term. Instead,
45 * - For first pixel, cast a whole bunch of rays in different directions
46 * from the object intersection point to see what the diffuse illumination
47 * really is. Save this value, after estimating its
48 * degree of reusability. (Method 1)
49 * - For second and subsequent pixels,
50 * - If there are one or more nearby values already computed,
51 * average them and use the result (Method 2), else
52 * - Use method 1.
53 *
54 * Implemented by and (c) 1994-6 Jim McElhiney, mcelhiney@acm.org or 71201,1326
55 * All standard POV distribution rights granted. All other rights reserved.
56 *************************************************************************/
57
58 #include <string.h>
59 #include <algorithm>
60
61 #include "frame.h"
62 #include "lighting.h"
63 #include "vector.h"
64 #include "povray.h"
65 #include "render.h"
66 #include "texture.h"
67 #include "octree.h"
68 #include "radiosit.h"
69 #include "ray.h"
70 #include "colour.h"
71 #include "pov_util.h"
72
73 BEGIN_POV_NAMESPACE
74
75 USING_POV_BASE_NAMESPACE
76
77 int firstRadiosityPass;
78
79 /*****************************************************************************
80 * Local preprocessor defines
81 ******************************************************************************/
82
83 #define RAD_GRADIENT 1
84 #define SAW_METHOD 1
85 /* #define SIGMOID_METHOD 1 */
86 /* #define SHOW_SAMPLE_SPOTS 1 */ /* try this! bright spots at sample pts */
87 /* #define LOW_COUNT_BRIGHT 1 */ /* this will highlight areas of low density if no extra samples are taken in the final pass */
88
89 /*****************************************************************************
90 * Local typedefs
91 ******************************************************************************/
92
93
94
95 /*****************************************************************************
96 * Local variables
97 ******************************************************************************/
98
99
100 long ra_reuse_count = 0; // GLOBAL VARIABLE
101 long ra_gather_count = 0; // GLOBAL VARIABLE
102
103 int Radiosity_Trace_Level = 1; // GLOBAL VARIABLE
104
105 COLOUR Radiosity_Gather_Total; // GLOBAL VARIABLE
106 long Radiosity_Gather_Total_Count; // GLOBAL VARIABLE
107 COLOUR Radiosity_Setting_Total; // GLOBAL VARIABLE
108 long Radiosity_Setting_Total_Count; // GLOBAL VARIABLE
109
110 #ifdef RADSTATS
111 extern long ot_blockcount; // GLOBAL VARIABLE
112 long ot_seenodecount = 0; // GLOBAL VARIABLE
113 long ot_seeblockcount = 0; // GLOBAL VARIABLE
114 long ot_doblockcount = 0; // GLOBAL VARIABLE
115 long ot_dotokcount = 0; // GLOBAL VARIABLE
116 long ot_lastcount = 0; // GLOBAL VARIABLE
117 long ot_lowerrorcount = 0; // GLOBAL VARIABLE
118 #endif
119
120 VECTOR *fast_rad_samples = NULL; // GLOBAL VARIABLE
121
122 OT_NODE *ot_root = NULL; // GLOBAL VARIABLE
123
124 /* This (and all other changing globals) should really be in an execution
125 * context structure passed down the execution call tree as a parameter to
126 * each function. This would allow for a multiprocessor/multithreaded version.
127 */
128 OStream *ot_fd = NULL; // GLOBAL VARIABLE
129
130
131 /*****************************************************************************
132 * Static functions
133 ******************************************************************************/
134
135 static int ra_reuse (VECTOR IPoint, VECTOR S_Normal, COLOUR Illuminance);
136 static int ra_average_near (OT_BLOCK *block, void *void_info);
137 static void ra_gather (VECTOR IPoint, VECTOR Raw_Normal, VECTOR LayNormal2, COLOUR Illuminance, DBL Weight);
138 static void VUnpack (VECTOR dest_vec, const BYTE_XYZ *pack);
139
140
141
142 /*****************************************************************************
143 *
144 * FUNCTION
145 *
146 * Compute_Ambient
147 *
148 * INPUT
149 *
150 * OUTPUT
151 *
152 * RETURNS
153 *
154 * AUTHOUR
155 *
156 * Jim McElhiney
157 *
158 * DESCRIPTION
159 *
160 * Main entry point for calculated diffuse illumination
161 *
162 * CHANGES
163 *
164 * --- 1994 : Creation.
165 *
166 ******************************************************************************/
167 /* the colour to be calculated */
168 /* maximum possible contribution to pixel colour */
169 /* NK rad 22 Nov 1999 - added LayNormal */
Compute_Ambient(VECTOR IPoint,VECTOR Raw_Normal,VECTOR LayNormal,COLOUR Ambient_Colour,DBL Weight)170 int Compute_Ambient(VECTOR IPoint, VECTOR Raw_Normal, VECTOR LayNormal, COLOUR Ambient_Colour, DBL Weight)
171 {
172 int retval, reuse;
173 DBL grey, save_bound;
174
175 save_bound = opts.Radiosity_Error_Bound;
176 if ( Weight < .25 )
177 {
178 opts.Radiosity_Error_Bound += (.25 - Weight);
179 }
180 /* NK rad 22 Nov 1999 - switched to LayNormal */
181 reuse = ra_reuse(IPoint, LayNormal, Ambient_Colour);
182 opts.Radiosity_Error_Bound = save_bound;
183
184
185 /* allow more samples on final pass - unless user says not to */
186 if((reuse >= opts.Radiosity_Nearest_Count ) ||
187 ((firstRadiosityPass == false) && (reuse > 0) && !opts.Radiosity_Add_On_Final_Trace))
188 {
189 ra_reuse_count++;
190 retval = 0;
191 #ifdef LOW_COUNT_BRIGHT
192 /* use this for testing - it will tell you where too few are found */
193 if (reuse<opts.Radiosity_Nearest_Count)
194 {
195 Ambient_Colour[0] = 4.;
196 Ambient_Colour[1] = 4.;
197 Ambient_Colour[2] = 4.;
198 }
199 #endif
200 }
201 else
202 {
203 ra_gather(IPoint, Raw_Normal, LayNormal, Ambient_Colour, Weight);
204
205 /* NK rad - always use reuse - avoids bright/dark dots */
206 /* NK rad 22 Nov 1999 - switched to LayNormal */
207 reuse=ra_reuse(IPoint, LayNormal, Ambient_Colour);
208
209 ra_gather_count++; /* keep a running count */
210
211 retval = 1;
212 }
213
214 grey = GREY_SCALE(Ambient_Colour);
215
216 /* note grey spelling: american options structure with worldbeat calculations! */
217 Ambient_Colour[pRED] = opts.Radiosity_Gray * grey + Ambient_Colour[pRED] * (1.-opts.Radiosity_Gray);
218 Ambient_Colour[pGREEN] = opts.Radiosity_Gray * grey + Ambient_Colour[pGREEN] * (1.-opts.Radiosity_Gray);
219 Ambient_Colour[pBLUE] = opts.Radiosity_Gray * grey + Ambient_Colour[pBLUE] * (1.-opts.Radiosity_Gray);
220
221 /* Scale up by current brightness factor prior to return */
222 VScale(Ambient_Colour, Ambient_Colour, opts.Radiosity_Brightness);
223
224 return(retval);
225 }
226
227
228
229 /*****************************************************************************
230 *
231 * FUNCTION
232 *
233 * ra_reuse
234 *
235 * INPUT
236 *
237 * OUTPUT
238 *
239 * RETURNS
240 *
241 * AUTHOUR
242 *
243 * Jim McElhiney
244 *
245 * DESCRIPTION
246 *
247 * Returns whether or not there were some prestored values close enough to
248 * reuse.
249 *
250 * CHANGES
251 *
252 * --- 1994 : Creation.
253 *
254 ******************************************************************************/
255
ra_reuse(VECTOR IPoint,VECTOR S_Normal,COLOUR Illuminance)256 static int ra_reuse(VECTOR IPoint, VECTOR S_Normal, COLOUR Illuminance)
257 {
258 int i;
259 WT_AVG gather;
260
261 if (ot_root != NULL)
262 {
263 Make_Colour(gather.Weights_Times_Illuminances, 0.0, 0.0, 0.0);
264
265 gather.Weights = 0.0;
266
267 Assign_Vector(gather.P, IPoint);
268 Assign_Vector(gather.N, S_Normal);
269
270 gather.Weights_Count = 0;
271 gather.Good_Count = 0;
272 gather.Close_Count = 0;
273 gather.Current_Error_Bound = opts.Radiosity_Error_Bound;
274
275 for (i = 1; i < Radiosity_Trace_Level; i++)
276 {
277 gather.Current_Error_Bound *= 2;
278 }
279
280 /*
281 * Go through the tree calculating a weighted average of all of the
282 * usable points near this one
283 */
284
285 ot_dist_traverse(ot_root, IPoint, Radiosity_Trace_Level,
286 ra_average_near, (void *)&gather);
287
288 /* Did we get any nearby points we could reuse? */
289
290 if (gather.Good_Count > 0)
291 {
292 /* NK rad - Average together all of the samples (sums were returned by
293 ot_dist_traverse). We are using nearest_count as a lower bound,
294 not an upper bound.
295 */
296
297 VInverseScale(Illuminance, gather.Weights_Times_Illuminances, gather.Weights);
298 }
299 }
300 else
301 {
302 gather.Good_Count = 0; /* No tree, so no reused values */
303 }
304
305 return(gather.Good_Count);
306 }
307
308
309
310 /*****************************************************************************
311 *
312 * FUNCTION
313 *
314 * ra_average_near
315 *
316 * INPUT
317 *
318 * OUTPUT
319 *
320 * RETURNS
321 *
322 * AUTHOUR
323 *
324 * Jim McElhiney
325 *
326 * DESCRIPTION
327 *
328 * Tree traversal function used by ra_reuse()
329 * Calculate the weight of this cached value, taking into account how far
330 * it is from our test point, and the difference in surface normal angles.
331 *
332 * Given a node with an old cached value, check to see if it is reusable, and
333 * aggregate its info into the weighted average being built during the tree
334 * traversal. block contains Point, Normal, Illuminance,
335 * Harmonic_Mean_Distance
336 *
337 * CHANGES
338 *
339 * --- 1994 : Creation.
340 *
341 ******************************************************************************/
342
ra_average_near(OT_BLOCK * block,void * void_info)343 static int ra_average_near(OT_BLOCK *block, void *void_info)
344 {
345 /* int ind, i;*/
346 WT_AVG *info = (WT_AVG *) void_info;
347 VECTOR half, delta, delta_unit;
348 COLOUR tc, prediction;
349 DBL ri, error_reuse, dir_diff, in_front, dist, weight, square_dist, dr, dg, db;
350 DBL error_reuse_rotate, error_reuse_translate, inverse_dist, cos_diff_from_nearest;
351 DBL quickcheck_rad;
352
353
354 #ifdef RADSTATS
355 ot_doblockcount++;
356 #endif
357
358 VSub(delta, info->P, block->Point); /* a = b - c, which is test p minus old pt */
359
360 square_dist = VSumSqr(delta);
361
362 quickcheck_rad = (DBL)block->Harmonic_Mean_Distance * info->Current_Error_Bound;
363
364 /* first we do a tuning test--this func gets called a LOT */
365 if (square_dist < quickcheck_rad * quickcheck_rad)
366 {
367
368 dist = sqrt(square_dist);
369 ri = (DBL)block->Harmonic_Mean_Distance;
370
371
372 if ( dist > .000001 )
373 {
374 inverse_dist = 1./dist;
375 VScale(delta_unit, delta, inverse_dist); /* this is a normalization */
376
377 /* This block reduces the radius of influence when it points near the nearest
378 surface found during sampling. */
379 VDot( cos_diff_from_nearest, block->To_Nearest_Surface, delta_unit);
380 if ( cos_diff_from_nearest > 0. )
381 {
382 ri = cos_diff_from_nearest * (DBL)block->Nearest_Distance +
383 (1.-cos_diff_from_nearest) * ri;
384 }
385 }
386
387 if (dist < ri * info->Current_Error_Bound)
388 {
389 VDot(dir_diff, info->N, block->S_Normal);
390
391 /* NB error_reuse varies from 0 to 3.82 (1+ 2 root 2) */
392
393 error_reuse_translate = dist / ri;
394
395 error_reuse_rotate = 2.0 * sqrt(fabs(1.0 - dir_diff));
396
397 error_reuse = error_reuse_translate + error_reuse_rotate;
398
399 /* is this old point within a reasonable error distance? */
400
401 if (error_reuse < info->Current_Error_Bound)
402 {
403 #ifdef RADSTATS
404 ot_lowerrorcount++;
405 #endif
406
407 if (dist > 0.000001)
408 {
409 /*
410 * Make sure that the old point is not in front of this point, the
411 * old surface might shadow this point and make the result
412 * meaningless
413 */
414 VHalf(half, info->N, block->S_Normal);
415 VNormalizeEq(half); /* needed so check can be to constant */
416 VDot(in_front, delta_unit, half);
417 }
418 else
419 {
420 in_front = 1.0;
421 }
422
423 /*
424 * Theory: eliminate the use of old points well in front of our
425 * new point we are calculating, but not ones which are just a little
426 * tiny bit in front. This (usually) avoids eliminating points on the
427 * same surface by accident.
428 */
429
430 if (in_front > (-0.05))
431 {
432 #ifdef RADSTATS
433 ot_dotokcount++;
434 #endif
435
436 #ifdef SIGMOID_METHOD
437 weight = error_reuse / info->Current_Error_Bound; /* 0 < t < 1 */
438 weight = (cos(weight * M_PI) + 1.0) * 0.5; /* 0 < w < 1 */
439 #endif
440 #ifdef SAW_METHOD
441 weight = 1.0 - (error_reuse / info->Current_Error_Bound); /* 0 < t < 1 */
442 weight = sqrt(sqrt(weight)); /* less splotchy */
443 /*weight = sqrt(sqrt(sqrt(weight))); maybe even less splotchy */
444 /*weight = weight*weight*weight*weight*weight; more splotchy */
445 #endif
446
447 if ( weight > 0.001 )
448 { /* avoid floating point oddities near zero */
449
450 /* This is the block where we use the gradient to improve the prediction */
451 dr = delta[X] * block->drdx + delta[Y] * block->drdy + delta[Z] * block->drdz;
452 dg = delta[X] * block->dgdx + delta[Y] * block->dgdy + delta[Z] * block->dgdz;
453 db = delta[X] * block->dbdx + delta[Y] * block->dbdy + delta[Z] * block->dbdz;
454 #ifndef RAD_GRADIENT
455 dr = dg = db = 0.;
456 #endif
457 #if 0
458 /* Ensure that the total change in colour is a reasonable magnitude */
459 if ( dr > .1 ) dr = .1; else if ( dr < -.1 ) dr = -.1;
460 if ( dg > .1 ) dg = .1; else if ( dg < -.1 ) dg = -.1;
461 if ( db > .1 ) db = .1; else if ( db < -.1 ) db = -.1;
462 #endif
463 // NK 6-May-2003 removed clipping - not sure why it was here in the
464 // first place, but it sure causes problems for HDR scenes, and removing
465 // it doesn't seem to cause problems for non-HRD scenes.
466 // But we want to make sure that our deltas don't cause a positive illumination
467 // to go below zero, while allowing negative illuminations to stay negative.
468 if ((dr + block->Illuminance[pRED] < 0.0) && block->Illuminance[pRED]>0.0) dr = -block->Illuminance[pRED];
469 if ((dg + block->Illuminance[pGREEN] < 0.0) && block->Illuminance[pGREEN]>0.0) dg = -block->Illuminance[pGREEN];
470 if ((db + block->Illuminance[pBLUE] < 0.0) && block->Illuminance[pBLUE]>0.0) db = -block->Illuminance[pBLUE];
471
472 prediction[pRED] = block->Illuminance[pRED] + dr;
473 //prediction[pRED] = max(prediction[pRED], (COLC)0.0), (COLC)1.0);
474 prediction[pGREEN] = block->Illuminance[pGREEN] + dg;
475 //prediction[pGREEN] = max(prediction[pGREEN], (COLC)0.0), (COLC)1.0);
476 prediction[pBLUE] = block->Illuminance[pBLUE] + db;
477 //prediction[pBLUE] = max(prediction[pBLUE], (COLC)0.0), (COLC)1.0);
478
479 #ifdef SHOW_SAMPLE_SPOTS
480 if ( dist < opts.Radiosity_Dist_Max * .015 ) {
481 prediction[pRED] = prediction[pGREEN] = prediction[pBLUE] = 3.;
482 }
483 #endif
484
485 /* The predicted colour is an extrapolation based on the old value */
486 VScale(tc, prediction, weight);
487
488 VAddEq(info->Weights_Times_Illuminances, tc);
489
490 info->Weights += weight;
491
492 info->Weights_Count++;
493
494 info->Good_Count++;
495
496 /* NK rad - it fit in the error bound, so keep it. We use all
497 that fit the error bounding criteria. There is no need to put
498 a maximum on the number of samples that are averaged.
499 */
500 }
501 }
502 }
503 }
504 }
505
506 return(1);
507 }
508
509
510
511 /*****************************************************************************
512 *
513 * FUNCTION
514 *
515 * ra_gather
516 *
517 * INPUT
518 * IPoint - a point at which the illumination is needed
519 * Raw_Normal - the surface normal (not perturbed by the current layer) at that point
520 * Illuminance - a place to put the return result
521 * Weight - the weight of this point in final output, to drive ADC_Bailout
522 *
523 * OUTPUT
524 * The average colour of light of objects visible from the specified point.
525 * The colour is returned in the Illuminance parameter.
526 *
527 *
528 * RETURNS
529 *
530 * AUTHOUR
531 *
532 * Jim McElhiney
533 *
534 * DESCRIPTION
535 * Gather up the incident light and average it.
536 * Return the results in Illuminance, and also cache them for later.
537 * Note that last parameter is similar to weight parameter used
538 * to control ADC_Bailout as a parameter to Trace(), but it also
539 * takes into account that this subsystem calculates only ambient
540 * values. Therefore, coming in at the top level, the value might
541 * be 0.3 if the first object hit had an ambient of 0.3, whereas
542 * Trace() would have been passed a parameter of 1.0 (since it
543 * calculates the whole pixel value).
544 *
545 * CHANGES
546 *
547 * --- 1994 : Creation.
548 *
549 ******************************************************************************/
ra_gather(VECTOR IPoint,VECTOR Raw_Normal,VECTOR LayNormal2,COLOUR Illuminance,DBL Weight)550 static void ra_gather(VECTOR IPoint, VECTOR Raw_Normal, VECTOR LayNormal2, COLOUR Illuminance, DBL Weight)
551 {
552 extern FRAME Frame;
553 int i, hit, Current_Radiosity_Count;
554 unsigned int Save_Quality_Flags, Save_Options;
555 VECTOR random_vec, direction, up, min_dist_vec;
556 int save_Max_Trace_Level;
557 DBL Inverse_Distance_Sum, depth, mean_dist, weight, save_min_reuse,
558 drdxs, dgdxs, dbdxs, drdys, dgdys, dbdys, drdzs, dgdzs, dbdzs,
559 depth_weight_for_this_gradient, dxsquared, dysquared, dzsquared,
560 constant_term, deemed_depth, min_dist, reuse_dist_min, to_eye,
561 sum_of_inverse_dist, sum_of_dist, average_dist, gradient_count;
562 COLOUR Colour_Sums, Temp_Colour;
563 RAY New_Ray;
564 OT_BLOCK *block;
565 OT_ID id;
566 int sampleNum;
567 VECTOR n2,n3;
568
569
570 int save_nearest_count;
571
572 DBL max_ill;
573 DBL save_dist_max;
574 DBL save_adc_bailout;
575
576 VECTOR LayNormal;
577
578 /* we might change laynormal, so make sure we're using our local copy */
579 Assign_Vector(LayNormal,LayNormal2);
580
581 /*
582 * A bit of theory: The goal is to create a set of "random" direction rays
583 * so that the probability of close-to-normal versus close-to-tangent rolls
584 * off in a cos-theta curve, where theta is the deviation from normal.
585 * That is, lots of rays close to normal, and very few close to tangent.
586 * You also want to have all of the rays be evenly spread, no matter how
587 * many you want to use. The lookup array has an array of points carefully
588 * chosen to meet all of these criteria.
589 */
590
591
592 /* The number of rays to trace varies with our recursion depth */
593
594 Current_Radiosity_Count = opts.Radiosity_Count;
595 save_min_reuse = opts.Radiosity_Min_Reuse;
596 save_nearest_count = opts.Radiosity_Nearest_Count;
597 save_dist_max = opts.Radiosity_Dist_Max;
598
599 /* NK rad - use different adc_bailout for radiosity calculations */
600 save_adc_bailout = ADC_Bailout;
601 if (Radiosity_Trace_Level==1)
602 {
603 ADC_Bailout *= opts.Radiosity_ADC_Bailout;
604 }
605
606 /* NK rad - compute dist_max on the fly */
607 /* we really should use the ray's footprint here... but that means we need
608 to keep track of the ray footprints, first */
609 VDist(opts.Radiosity_Dist_Max, Frame.Camera->Location, IPoint);
610 opts.Radiosity_Dist_Max *= 0.2;
611
612 /*for ( i=1; i<Radiosity_Trace_Level; i++ )*/
613 if (Radiosity_Trace_Level>1)
614 {
615 Current_Radiosity_Count /= 3;
616 opts.Radiosity_Min_Reuse *= 2.;
617 opts.Radiosity_Nearest_Count = 2;
618 opts.Radiosity_Dist_Max *= 2;
619 }
620 if (Radiosity_Trace_Level>2)
621 {
622 Current_Radiosity_Count /= 2;
623 opts.Radiosity_Nearest_Count = 1;
624 opts.Radiosity_Dist_Max *= 2;
625 }
626 if (Current_Radiosity_Count<5) Current_Radiosity_Count=5;
627
628 /* Save some global stuff which we have to change for now */
629 save_Max_Trace_Level = Max_Trace_Level;
630
631 // adjust the max_trace_level
632 Max_Trace_Level = Trace_Level + opts.Radiosity_Recursion_Limit + 1;
633 // but make sure it doesn't exceed the original max trace level
634 if (Max_Trace_Level>save_Max_Trace_Level) Max_Trace_Level = save_Max_Trace_Level;
635 // if (Max_Trace_Level>MAX_TRACE_LEVEL_LIMIT) Max_Trace_Level = MAX_TRACE_LEVEL_LIMIT;
636
637 /* Since we'll be calculating averages, zero the accumulators */
638
639 Make_Colour(Colour_Sums, 0., 0., 0.);
640 Inverse_Distance_Sum = 0.;
641
642
643 min_dist = BOUND_HUGE;
644
645 if ( fabs(fabs(LayNormal[Z])- 1.) < .1 ) {
646 /* too close to vertical for comfort, so use cross product with horizon */
647 up[X] = 0.; up[Y] = 1.; up[Z] = 0.;
648 }
649 else
650 {
651 up[X] = 0.; up[Y] = 0.; up[Z] = 1.;
652 }
653
654 VCross(n2, LayNormal, up); VNormalizeEq(n2);
655 VCross(n3, LayNormal, n2); VNormalizeEq(n3);
656
657
658 /* Note that this max() forces at least one ray to be shot.
659 Otherwise, the loop does nothing, since every call to
660 Trace() just bails out immediately! */
661 weight = max(ADC_Bailout, Weight/(DBL)Current_Radiosity_Count);
662
663 /* Initialized the accumulators for the integrals which will be come the rad gradient */
664 drdxs = dgdxs = dbdxs = drdys = dgdys = dbdys = drdzs = dgdzs = dbdzs = 0.;
665 sum_of_inverse_dist = sum_of_dist = gradient_count = 0.;
666
667 for (i = sampleNum = hit = 0; i < Current_Radiosity_Count; i++)
668 {
669 //DBL scale;
670 DBL rayOk=-1;
671 int lockupTest = 0;
672
673 /* loop through here choosing rays until we get one that is not behind
674 the surface */
675
676 while(rayOk<=0 && lockupTest<1600)
677 {
678 lockupTest++;
679 /*/Increase_Counter(stats[Gather_Performed_Count]);*/
680 Assign_Vector(random_vec, fast_rad_samples[sampleNum++]);
681 //scale = 1.0;
682
683 /* don't go beyond range of samples */
684 if (sampleNum>=1600) sampleNum = 0;
685
686 /* if we've taken too many incorrect samples, use raw_normal instead of laynormal */
687 if (sampleNum>Current_Radiosity_Count*5)
688 {
689 Assign_Vector(LayNormal,Raw_Normal);
690 }
691
692 if ( fabs(LayNormal[Z] - 1.) < .001 ) /* pretty well straight Z, folks */
693 {
694 /* we are within 1/20 degree of pointing in the Z axis. */
695 /* use all vectors as is--they're precomputed this way */
696 Assign_Vector(direction, random_vec);
697 }
698 else
699 {
700 direction[X] = n2[X]*random_vec[X] + n3[X]*random_vec[Y] + LayNormal[X]*random_vec[Z];
701 direction[Y] = n2[Y]*random_vec[X] + n3[Y]*random_vec[Y] + LayNormal[Y]*random_vec[Z];
702 direction[Z] = n2[Z]*random_vec[X] + n3[Z]*random_vec[Y] + LayNormal[Z]*random_vec[Z];
703 }
704
705 /* make sure we don't go behind raw_normal */
706 VDot(rayOk,direction,Raw_Normal);
707 }
708
709 /* Build a ray pointing in the chosen direction */
710 Make_Colour(Temp_Colour, 0.0, 0.0, 0.0);
711 Initialize_Ray_Containers(&New_Ray);
712 Assign_Vector(New_Ray.Initial, IPoint);
713 Assign_Vector(New_Ray.Direction, direction);
714
715 /* save some flags that must be set to a different value during the trace() */
716 Save_Quality_Flags = opts.Quality_Flags;
717 Save_Options = opts.Options;
718
719 opts.Radiosity_Quality = 6;
720
721 #ifdef SAFE_BUT_SLOW
722 opts.Quality_Flags = Quality_Values[opts.Radiosity_Quality];
723 #else
724 /* Set up a custom quality level with no area lights or light buffer */
725 opts.Options &= ~USE_LIGHT_BUFFER;
726 opts.Quality_Flags &= ~Q_AREA_LIGHT;
727 if(!opts.Radiosity_Use_Media)
728 {
729 opts.Quality_Flags &= ~Q_VOLUME;
730 }
731 #endif
732
733 /* Go down in recursion, trace the result, and come back up */
734
735 Trace_Level++;
736 Radiosity_Trace_Level++;
737 depth = Trace(&New_Ray, Temp_Colour, weight);
738 Radiosity_Trace_Level--;
739 Trace_Level--;
740
741
742 /* NK rad - each sample is limited to a user-specified brightness */
743 /* this is necessary to fix problems splotchiness caused by very
744 bright objects */
745 /* changed lighting.c to ignore phong/specular if tracing radiosity beam */
746 max_ill = max3(Temp_Colour[0],Temp_Colour[1],Temp_Colour[2]);
747
748 if(max_ill>opts.Maximum_Sample_Brightness
749 && opts.Maximum_Sample_Brightness > 0.0)
750 {
751 max_ill = opts.Maximum_Sample_Brightness/max_ill;
752 VScaleEq(Temp_Colour, max_ill);
753 }
754
755 /* let's wait and scale it after we do the gradients
756 VScaleEq(Temp_Colour, scale);*/
757
758 /* Add into illumination gradient integrals */
759
760 deemed_depth = depth;
761 if (deemed_depth < opts.Radiosity_Dist_Max * 10.)
762 {
763 depth_weight_for_this_gradient = 1. / deemed_depth;
764 sum_of_inverse_dist += 1. / deemed_depth;
765 sum_of_dist += deemed_depth;
766 gradient_count++;
767
768 dxsquared = direction[X] * direction[X]; if (direction[X] < 0.) dxsquared = -dxsquared;
769 dysquared = direction[Y] * direction[Y]; if (direction[Y] < 0.) dysquared = -dysquared;
770 dzsquared = direction[Z] * direction[Z]; if (direction[Z] < 0.) dzsquared = -dzsquared;
771 drdxs += dxsquared * Temp_Colour[pRED] * depth_weight_for_this_gradient;
772 dgdxs += dxsquared * Temp_Colour[pGREEN] * depth_weight_for_this_gradient;
773 dbdxs += dxsquared * Temp_Colour[pBLUE] * depth_weight_for_this_gradient;
774 drdys += dysquared * Temp_Colour[pRED] * depth_weight_for_this_gradient;
775 dgdys += dysquared * Temp_Colour[pGREEN] * depth_weight_for_this_gradient;
776 dbdys += dysquared * Temp_Colour[pBLUE] * depth_weight_for_this_gradient;
777 drdzs += dzsquared * Temp_Colour[pRED] * depth_weight_for_this_gradient;
778 dgdzs += dzsquared * Temp_Colour[pGREEN] * depth_weight_for_this_gradient;
779 dbdzs += dzsquared * Temp_Colour[pBLUE] * depth_weight_for_this_gradient;
780 }
781
782
783 if (depth > opts.Radiosity_Dist_Max)
784 {
785 depth = opts.Radiosity_Dist_Max;
786 }
787 else
788 {
789 #ifdef RADSTATS
790 hit++;
791 #endif
792 }
793
794 if (depth < min_dist)
795 {
796 min_dist = depth;
797 Assign_Vector(min_dist_vec, direction);
798 }
799
800 opts.Quality_Flags = Save_Quality_Flags;
801 opts.Options = Save_Options;
802
803 /* Add into total illumination integral */
804
805 /* Ok, now we will scale the color (previous scale is commented out) */
806 VAddEq(Colour_Sums, Temp_Colour);
807 /* scale is always 1.0 here, so we just addeq instead of scaling!
808 VAddScaledEq(Colour_Sums, scale, Temp_Colour); */
809
810 Inverse_Distance_Sum += 1.0 / depth;
811 } /* end ray sampling loop */
812
813
814 /*
815 * Use the accumulated values to calculate the averages needed. The sphere
816 * of influence of this primary-method sample point is based on the
817 * harmonic mean distance to the points encountered. (An harmonic mean is
818 * the inverse of the mean of the inverses).
819 */
820
821 mean_dist = 1.0 / (Inverse_Distance_Sum / (DBL) Current_Radiosity_Count);
822
823 VInverseScale(Illuminance, Colour_Sums, (DBL) Current_Radiosity_Count);
824
825 /* Keep a running total of the final Illuminances we calculated */
826 if ( Radiosity_Trace_Level == 1) {
827 VAddEq(Radiosity_Gather_Total, Illuminance);
828 Radiosity_Gather_Total_Count++;
829 }
830
831
832 /* We want to cached this block for later reuse. But,
833 * if ground units not big enough, meaning that the value has very
834 * limited reuse potential, forget it.
835 */
836
837 if (Radiosity_Trace_Level == 1 || Current_Radiosity_Count>=5)
838 if (mean_dist > (opts.Radiosity_Dist_Max * 0.0001))
839 {
840
841 /*
842 * Theory: we don't want to calculate a primary method ray loop at every
843 * point along the inside edges, so a minimum effectivity is practical.
844 * It is expressed as a fraction of the distance to the eyepoint. 1/2%
845 * is a good number. This enhancement was Greg Ward's idea, but the use
846 * of % units is my idea. [JDM]
847 */
848
849 VDist(to_eye, Frame.Camera->Location, IPoint);
850 reuse_dist_min = to_eye * opts.Radiosity_Min_Reuse;
851 if (mean_dist < reuse_dist_min)
852 {
853 mean_dist = reuse_dist_min;
854 }
855
856
857 /* figure out the block id */
858 ot_index_sphere(IPoint, mean_dist * opts.Real_Radiosity_Error_Bound, &id);
859
860
861 #ifdef RADSTATS
862 ot_blockcount++;
863 #endif
864
865 /* After end of ray loop, we've decided that this point is worth storing */
866 /* Allocate a block, and fill it with values for reuse in cacheing later */
867 block = (OT_BLOCK *)POV_MALLOC(sizeof(OT_BLOCK), "octree block");
868 memset(block, 0, sizeof(OT_BLOCK));
869
870 /* beta */
871 if ( gradient_count > 10)
872 {
873 average_dist = sum_of_dist / gradient_count;
874 constant_term = 1.00 / (sum_of_inverse_dist * average_dist );
875
876 block->drdx = (float)(drdxs * constant_term);
877 block->dgdx = (float)(dgdxs * constant_term);
878 block->dbdx = (float)(dbdxs * constant_term);
879 block->drdy = (float)(drdys * constant_term);
880 block->dgdy = (float)(dgdys * constant_term);
881 block->dbdy = (float)(dbdys * constant_term);
882 block->drdz = (float)(drdzs * constant_term);
883 block->dgdz = (float)(dgdzs * constant_term);
884 block->dbdz = (float)(dbdzs * constant_term);
885 }
886
887
888 /* Fill up the values in the octree (ot_) cache block */
889
890 Assign_RGB(block->Illuminance, Illuminance);
891 Assign_Vector(block->To_Nearest_Surface, min_dist_vec);
892 block->Harmonic_Mean_Distance = (float)mean_dist;
893 block->Nearest_Distance = (float)min_dist;
894 block->Bounce_Depth = (short)Radiosity_Trace_Level;
895 Assign_Vector(block->Point, IPoint);
896 Assign_Vector(block->S_Normal, LayNormal);
897 block->next = NULL;
898
899 /* store the info block in the oct tree */
900 ot_ins(&ot_root, block, &id);
901
902 /* In case the rendering is suspended, save the cache tree values to a file */
903 if ( opts.Radiosity_File_SaveWhileRendering && (ot_fd != NULL) ) {
904 ot_write_block(block, ot_fd);
905 }
906 }
907
908 /* Put things back where they were in recursion depth */
909 Max_Trace_Level = save_Max_Trace_Level;
910 opts.Radiosity_Min_Reuse = save_min_reuse;
911 opts.Radiosity_Nearest_Count = save_nearest_count;
912 opts.Radiosity_Dist_Max = save_dist_max;
913 /* NK rad - put back adc_bailout */
914 ADC_Bailout = save_adc_bailout;
915 /* NK ---- */
916
917 }
918
919
920 /*****************************************************************************
921 *
922 * FUNCTION VUnpack() - Unpacks "pack_vec" into "dest_vec" and normalizes it.
923 *
924 * INPUT
925 *
926 * OUTPUT
927 *
928 * RETURNS Nothing
929 *
930 * AUTHOUR Jim McElhiney
931 *
932 * DESCRIPTION
933 *
934 * The precomputed radiosity rays are packed into a lookup array with one byte
935 * for each of dx, dy, and dz. dx and dy are scaled from the range (-1. to 1.),
936 * and dz is scaled from the range (0. to 1.), and both are stored in the range
937 * 0 to 255.
938 *
939 * The reason for this function is that it saves a bit of memory. There are 2000
940 * entries in the table, and packing them saves 21 bytes each, or 42KB.
941 *
942 * CHANGES
943 *
944 * --- Jan 1996 : Creation.
945 *
946 ******************************************************************************/
VUnpack(VECTOR dest_vec,const BYTE_XYZ * pack_vec)947 static void VUnpack(VECTOR dest_vec, const BYTE_XYZ * pack_vec)
948 {
949 dest_vec[X] = ((double)pack_vec->x * (1./ 255.))*2.-1.;
950 dest_vec[Y] = ((double)pack_vec->y * (1./ 255.))*2.-1.;
951 dest_vec[Z] = ((double)pack_vec->z * (1./ 255.));
952
953 VNormalizeEq(dest_vec); /* already good to about 1%, but we can do better */
954 }
955
956
957 /*****************************************************************************
958 *
959 * FUNCTION Initialize_Radiosity_Code
960 *
961 * INPUT Nothing.
962 *
963 * OUTPUT Sets various global states used by radiosity. Notably,
964 * ot_fd - the file identifier of the file used to save radiosity values
965 *
966 * RETURNS 1 for Success, 0 for failure (e.g., could not open cache file)
967 *
968 * AUTHOUR Jim McElhiney
969 *
970 * DESCRIPTION
971 *
972 * CHANGES
973 *
974 * --- Jan 1996 : Creation.
975 *
976 ******************************************************************************/
Initialize_Radiosity_Code()977 bool Initialize_Radiosity_Code()
978 {
979 bool retval;
980 bool used_existing_file;
981 IStream *fd;
982 char rad_cache_filename[256];
983 int i;
984
985 retval = true; /* assume the best */
986
987 fast_rad_samples = (VECTOR *)POV_MALLOC(sizeof(VECTOR) * 1600, "Radiosity sample data");
988
989 for(i = 0; i < 1600; i++)
990 {
991 VUnpack(fast_rad_samples[i], &rad_samples[i]);
992 }
993
994 // always clear these even if radiosity isn't enabled, as otherwise
995 // we get misleading statistics on subsequent non-radiosity renders
996 opts.Radiosity_Preview_Done = 0;
997 ra_gather_count = 0;
998 ra_reuse_count = 0;
999
1000 #ifdef RADSTATS
1001 ot_seenodecount = 0;
1002 ot_seeblockcount = 0;
1003 ot_doblockcount = 0;
1004 ot_dotokcount = 0;
1005 ot_lowerrorcount = 0;
1006 ot_lastcount = 0;
1007 #endif
1008
1009 if ( opts.Radiosity_Enabled)
1010 {
1011 if ( opts.Radiosity_Dist_Max == 0. )
1012 {
1013 /* User hasn't picked a radiosity dist max, so pick one automatically. */
1014 VDist(opts.Radiosity_Dist_Max, Frame.Camera->Location,
1015 Frame.Camera->Look_At);
1016 opts.Radiosity_Dist_Max *= 0.2;
1017 }
1018
1019 if ( ot_fd != NULL ) /* if already open for some unknown reason, close it */
1020 {
1021 delete ot_fd;
1022 ot_fd = 0;
1023 }
1024
1025 /* build the file name for the radiosity cache file */
1026 strcpy(rad_cache_filename, opts.Scene_Name);
1027 strcat(rad_cache_filename, RADIOSITY_CACHE_EXTENSION);
1028
1029 opts.Real_Radiosity_Error_Bound = opts.Radiosity_Error_Bound;
1030
1031 /* NK rad */
1032 if (opts.Radiosity_Load_File_Name)
1033 {
1034 fd = New_Checked_IStream(opts.Radiosity_Load_File_Name, POV_File_Data_RCA);
1035 if ( fd != NULL) {
1036 ot_read_file(fd);
1037 delete fd;
1038 }
1039 POV_FREE(opts.Radiosity_Load_File_Name);
1040 opts.Radiosity_Load_File_Name = NULL;
1041 }
1042 /* NK ---- */
1043
1044 used_existing_file = false;
1045 if ( ((opts.Options & CONTINUE_TRACE) && opts.Radiosity_File_ReadOnContinue) ||
1046 opts.Radiosity_File_AlwaysReadAtStart )
1047 {
1048 fd = New_Checked_IStream(rad_cache_filename, POV_File_Data_RCA); /* "myname.rca" */
1049 if ( fd != NULL) {
1050 used_existing_file = ot_read_file(fd);
1051 retval &= used_existing_file;
1052 delete fd;
1053 }
1054 }
1055 else
1056 {
1057 DELETE_FILE(rad_cache_filename); /* default case, force a clean start */
1058 }
1059
1060 if ( opts.Radiosity_File_SaveWhileRendering )
1061 {
1062 /* If we are writing a file, but not using what's there, we truncate,
1063 since we conclude that what is there is bad.
1064 But, if we are also using what's there, then it must be good, so
1065 we just append to it.
1066 */
1067 ot_fd = New_Checked_OStream(rad_cache_filename, POV_File_Data_RCA, used_existing_file);
1068 retval &= (ot_fd != NULL);
1069 }
1070
1071 }
1072
1073 return retval;
1074 }
1075
1076
1077 /*****************************************************************************
1078 *
1079 * FUNCTION Deinitialize_Radiosity_Code()
1080 *
1081 * INPUT Nothing.
1082 *
1083 * OUTPUT Sets various global states used by radiosity. Notably,
1084 * ot_fd - the file identifier of the file used to save radiosity values
1085 *
1086 * RETURNS 1 for total success, 0 otherwise (e.g., could not save cache tree)
1087 *
1088 * AUTHOUR Jim McElhiney
1089 *
1090 * DESCRIPTION
1091 * Wrap up and free any radiosity-specific features.
1092 * Note that this function is safe to call even if radiosity was not on.
1093 *
1094 * CHANGES
1095 *
1096 * --- Jan 1996 : Creation.
1097 *
1098 ******************************************************************************/
Deinitialize_Radiosity_Code()1099 bool Deinitialize_Radiosity_Code()
1100 {
1101 bool retval;
1102 char rad_cache_filename[256];
1103 OStream *fd;
1104
1105 retval = true; /* assume the best */
1106
1107 if ( opts.Radiosity_Enabled)
1108 {
1109 /* if the global file identifier is set, close it */
1110 if ( ot_fd != NULL ) {
1111 delete ot_fd;
1112 ot_fd = NULL;
1113 }
1114
1115
1116 /* build the file name for the radiosity cache file */
1117 strcpy(rad_cache_filename, opts.Scene_Name);
1118 strcat(rad_cache_filename, RADIOSITY_CACHE_EXTENSION);
1119
1120
1121 /* If user has not asked us to save the radiosity cache file, delete it */
1122 if ( opts.Radiosity_File_SaveWhileRendering &&
1123 !(opts.Radiosity_File_KeepAlways || (Stop_Flag && opts.Radiosity_File_KeepOnAbort) ) )
1124 {
1125 DELETE_FILE(rad_cache_filename);
1126 }
1127
1128 /* after-the-fact version. This is an alternative to putting a call to
1129 ot_write_node after the call to ot_ins in ra_gather().
1130 The on-the-fly version (all of the code which uses ot_fd) is superior
1131 in that you will get partial results if you restart your rendering
1132 with a different resolution or camera angle. This version is superior
1133 in that your rendering goes a lot quicker.
1134 */
1135
1136 /* NK rad */
1137 if (opts.Radiosity_Save_File_Name)
1138 {
1139 fd = New_Checked_OStream(opts.Radiosity_Save_File_Name, POV_File_Data_RCA, false);
1140 if ( fd != NULL ) {
1141 ot_save_tree(ot_root, fd);
1142 delete fd;
1143 }
1144 POV_FREE(opts.Radiosity_Save_File_Name);
1145 opts.Radiosity_Save_File_Name = NULL;
1146 }
1147 /* NK ---- */
1148
1149 if (!(opts.Radiosity_File_KeepAlways || (Stop_Flag && opts.Radiosity_File_KeepOnAbort)) &&
1150 !opts.Radiosity_File_SaveWhileRendering && ot_root != NULL )
1151 {
1152 fd = New_Checked_OStream(rad_cache_filename, POV_File_Data_RCA, false);
1153
1154 if ( fd != NULL ) {
1155 retval &= ot_save_tree(ot_root, fd);
1156
1157 delete fd;
1158 }
1159 else
1160 {
1161 retval = false;
1162 }
1163 }
1164
1165
1166 /* Note that multiframe animations should call this free function if they have
1167 moving objects and want correct results.
1168 They should NOT call this function if they have no moving objects (like
1169 fly-throughs) and want speed
1170 */
1171 if ( ot_root != NULL ) {
1172 retval &= ot_free_tree(&ot_root); /* this zeroes the root pointer */
1173 }
1174
1175 }
1176
1177 if(fast_rad_samples != NULL)
1178 {
1179 POV_FREE(fast_rad_samples);
1180 }
1181
1182 return retval;
1183 }
1184
1185 END_POV_NAMESPACE
1186