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