1 /****************************************************************************
2 * photons.cpp
3 *
4 * Author: Nathan Kopp
5 *
6 * This module implements Photon Mapping.
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: photons.cpp,v $
34 * $Revision: 1.19 $
35 * $Author: chris $
36 *
37 *****************************************************************************/
38
39
40 #include "frame.h"
41 #include "userio.h"
42 #include "povray.h"
43 #include "vector.h"
44 #include "texture.h" /* for FRAND() */
45 #include "matrices.h"
46 #include "objects.h"
47 #include "csg.h"
48 #include "octree.h"
49 #include "radiosit.h"
50 #include "photons.h"
51 #include "povms.h"
52 #include "povmsend.h"
53 #include "ray.h"
54 #include "pov_util.h"
55
56 #include <algorithm>
57
58 BEGIN_POV_NAMESPACE
59
60 /* ------------------------------------------------------ */
61 /* global variables */
62 /* ------------------------------------------------------ */
63 int backtraceFlag; // GLOBAL VARIABLE
64
65 PHOTON_OPTIONS photonOptions; // GLOBAL VARIABLE
66
67 int InitBacktraceWasCalled; // GLOBAL VARIABLE
68
69 // statistics helpers
70 int gPhotonStat_i = 0; // GLOBAL VARIABLE
71 int gPhotonStat_x_samples = 0; // GLOBAL VARIABLE
72 int gPhotonStat_y_samples = 0; // GLOBAL VARIABLE
73 int gPhotonStat_end = 0; // GLOBAL VARIABLE
74
75 /* ------------------------------------------------------ */
76 /* external variables */
77 /* ------------------------------------------------------ */
78 extern int Trace_Level; // GLOBAL VARIABLE
79 extern int disp_elem; // GLOBAL VARIABLE
80 extern int disp_nelems; // GLOBAL VARIABLE
81
82 #ifdef MOTION_BLUR_PATCH
83 extern int TimeStamp; // GLOBAL VARIABLE
84 #endif
85
86 /* ------------------------------------------------------ */
87 /* static functions */
88 /* ------------------------------------------------------ */
89 static PHOTON* AllocatePhoton(PHOTON_MAP *map);
90 static void FreePhotonMemory();
91 static void InitPhotonMemory();
92 static void sortAndSubdivide(int start, int end, int sorted);
93 static void buildTree(PHOTON_MAP *map);
94 static void ShootPhotonsAtObject(OBJECT *Object, LIGHT_SOURCE *Light, int count);
95 static void SearchThroughObjects(OBJECT *Object, LIGHT_SOURCE *Light, bool count);
96 static int savePhotonMap(void);
97 static int loadPhotonMap(void);
98 static void swapPhotons(int a, int b);
99 static void PQInsert(PHOTON *photon, DBL d);
100 #ifndef REMOVE_NOT_USED_CODE_WARNINGS_PATCH
101 static void PQDelMax(void);
102 #endif
103 static void gatherPhotonsRec(int start, int end);
104 static void setGatherOptions(PHOTON_MAP *map, int mediaMap);
105
106 /* ------------------------------------------------------ */
107 /* static variables */
108 /* ------------------------------------------------------ */
109 /*
110 These static variables are used to conserve stack space during
111 extensive recursion when gathering photons. All of these static
112 variables end in "_s".
113 */
114 static PHOTON **map_s; /* photon map */ // GLOBAL VARIABLE
115 static DBL size_sq_s; /* search radius squared */ // GLOBAL VARIABLE
116 static DBL Size_s; /* search radius (static) */ // GLOBAL VARIABLE
117 static DBL sqrt_dmax_s, dmax_s; /* dynamic search radius... current maximum */ // GLOBAL VARIABLE
118 static int TargetNum_s; /* how many to gather */ // GLOBAL VARIABLE
119 static DBL *pt_s; /* point around which we are gathering */ // GLOBAL VARIABLE
120 static int numfound_s; /* number of photons found */ // GLOBAL VARIABLE
121 static DBL *norm_s; /* surface normal */ // GLOBAL VARIABLE
122 static DBL flattenFactor; /* amount to flatten the spher to make it */ // GLOBAL VARIABLE
123 /* an ellipsoid when gathering photons */
124 /* zero = no flatten, one = regular */
125 static DBL photonCountEstimate; // GLOBAL VARIABLE
126
127 /*****************************************************************************
128
129 FUNCTION
130
131 CheckPassThru()
132
133 Checks to make sure that pass-through, high-density, and refraction
134 are not simultaneously selected. If all three are turned on, we need
135 to choose an appropriate one to turn off.
136
137 Preconditions:
138 'o' is an initialized object
139 'flag' is PH_PASSTHRU_FLAG, PH_TARGET_FLAG, or PH_RFR_ON_FLAG
140 (this is which flag was set most recently)
141
142 Postconditions:
143 One of these flags in 'o' is turned off, since they cannot all be turned on.
144
145 ******************************************************************************/
146
CheckPassThru(OBJECT * o,int flag)147 void CheckPassThru(OBJECT *o, int flag)
148 {
149 if( Test_Flag(o, PH_PASSTHRU_FLAG) &&
150 Test_Flag(o, PH_TARGET_FLAG) &&
151 !Test_Flag(o, PH_RFR_OFF_FLAG) )
152 {
153 switch (flag)
154 {
155 case PH_PASSTHRU_FLAG:
156 Warning(0, "Cannot use pass_through with refraction & target.\nTurning off refraction.");
157 Set_Flag(o, PH_RFR_OFF_FLAG);
158 Clear_Flag(o, PH_RFR_ON_FLAG);
159 break;
160
161 case PH_TARGET_FLAG:
162 if(Test_Flag(o, PH_RFR_ON_FLAG))
163 {
164 Warning(0, "Cannot use pass_through with refraction & target.\nTurning off pass_through.");
165 Clear_Flag(o,PH_PASSTHRU_FLAG);
166 }
167 else
168 {
169 Warning(0, "Cannot use pass_through with refraction & target.\nTurning off refraction.");
170 Set_Flag(o, PH_RFR_OFF_FLAG);
171 Clear_Flag(o, PH_RFR_ON_FLAG);
172 }
173 break;
174
175 case PH_RFR_ON_FLAG:
176 Warning(0, "Cannot use pass_through with refraction & target.\nTurning off pass_through.");
177 Clear_Flag(o, PH_PASSTHRU_FLAG);
178 break;
179 }
180 }
181 }
182
183 /*****************************************************************************
184
185 FUNCTION
186
187 InitBacktraceEverything()
188
189 Allocates memory.
190 Initializes all photon mapping stuff.
191 Does not create the photon map.
192
193 Preconditions: InitBacktraceEverything() not yet called
194 or
195 both InitBacktraceEverything() and FreeBacktraceEverything() called
196
197 Postconditions:
198 If photonOptions.photonsEnabled is true, then
199 memory for photon mapping is allocated.
200 else
201 nothing is done
202
203 ******************************************************************************/
204
InitBacktraceEverything()205 void InitBacktraceEverything()
206 {
207 int i;
208 double theta;
209
210 if (photonOptions.photonsEnabled)
211 {
212 InitBacktraceWasCalled = true;
213
214 photonOptions.photonMap.head = NULL;
215 photonOptions.photonMap.numPhotons = 0;
216 photonOptions.photonMap.numBlocks = 0;
217 #ifdef GLOBAL_PHOTONS
218 photonOptions.globalPhotonMap.head = NULL;
219 photonOptions.globalPhotonMap.numPhotons = 0;
220 photonOptions.globalPhotonMap.numBlocks = 0;
221 #endif
222 photonOptions.mediaPhotonMap.head = NULL;
223 photonOptions.mediaPhotonMap.numPhotons = 0;
224 photonOptions.mediaPhotonMap.numBlocks = 0;
225
226 photonOptions.photonGatherList = (PHOTON**)POV_MALLOC(sizeof(PHOTON *)*photonOptions.maxGatherCount, "Photon Map Info");
227 photonOptions.photonDistances = (DBL *)POV_MALLOC(sizeof(DBL)*photonOptions.maxGatherCount, "Photon Map Info");
228
229 InitPhotonMemory();
230
231 /* create the sin/cos arrays for speed */
232 /* range is -127..+127 => 0..254 */
233 photonOptions.sinTheta = (DBL *)POV_MALLOC(sizeof(DBL)*255, "Photon Map Info");
234 photonOptions.cosTheta = (DBL *)POV_MALLOC(sizeof(DBL)*255, "Photon Map Info");
235 for(i=0; i<255; i++)
236 {
237 theta = (double)(i-127)*M_PI/127.0;
238 photonOptions.sinTheta[i] = sin(theta);
239 photonOptions.cosTheta[i] = cos(theta);
240 }
241 }
242 }
243
244 /* savePhotonMap()
245
246 Saves the caustic photon map to a file.
247
248 Preconditions:
249 InitBacktraceEverything was called
250 the photon map has been built and balanced
251 photonOptions.fileName contains the filename to save
252
253 Postconditions:
254 Returns 1 if success, 0 if failure.
255 If success, the photon map has been written to the file.
256 */
savePhotonMap()257 static int savePhotonMap()
258 {
259 PHOTON *ph;
260 FILE *f;
261 int i, err;
262
263 f = fopen(photonOptions.fileName, "wb");
264 if (!f) return 0;
265
266 /* caustic photons */
267 fwrite(&photonOptions.photonMap.numPhotons, sizeof(photonOptions.photonMap.numPhotons),1,f);
268 if (photonOptions.photonMap.numPhotons>0 && photonOptions.photonMap.head)
269 {
270 for(i=0; i<photonOptions.photonMap.numPhotons; i++)
271 {
272 ph = &(PHOTON_AMF(photonOptions.photonMap.head, i));
273 err = fwrite(ph, sizeof(PHOTON), 1, f);
274
275 if (err<=0)
276 {
277 /* fwrite returned an error! */
278 fclose(f);
279 return 0;
280 }
281 }
282 }
283
284 #ifdef GLOBAL_PHOTONS
285 /* global photons */
286 fwrite(&photonOptions.globalPhotonMap.numPhotons, sizeof(photonOptions.globalPhotonMap.numPhotons),1,f);
287 if (photonOptions.globalPhotonMap.numPhotons>0 && photonOptions.globalPhotonMap.head)
288 {
289 for(i=0; i<photonOptions.globalPhotonMap.numPhotons; i++)
290 {
291 ph = &(PHOTON_AMF(photonOptions.globalPhotonMap.head, i));
292 err = fwrite(ph, sizeof(PHOTON), 1, f);
293
294 if (err<=0)
295 {
296 /* fwrite returned an error! */
297 fclose(f);
298 return 0;
299 }
300 }
301 }
302 #endif
303
304 /* media photons */
305 fwrite(&photonOptions.mediaPhotonMap.numPhotons, sizeof(photonOptions.mediaPhotonMap.numPhotons),1,f);
306 if (photonOptions.mediaPhotonMap.numPhotons>0 && photonOptions.mediaPhotonMap.head)
307 {
308 for(i=0; i<photonOptions.mediaPhotonMap.numPhotons; i++)
309 {
310 ph = &(PHOTON_AMF(photonOptions.mediaPhotonMap.head, i));
311 err = fwrite(ph, sizeof(PHOTON), 1, f);
312
313 if (err<=0)
314 {
315 /* fwrite returned an error! */
316 fclose(f);
317 return 0;
318 }
319 }
320 }
321
322 fclose(f);
323 return true;
324 }
325
326 /* loadPhotonMap()
327
328 Loads the caustic photon map from a file.
329
330 Preconditions:
331 InitBacktraceEverything was called
332 the photon map is empty
333 photonOptions.fileName contains the filename to load
334
335 Postconditions:
336 Returns 1 if success, 0 if failure.
337 If success, the photon map has been loaded from the file.
338 If failure then the render should stop with an error
339 */
loadPhotonMap()340 static int loadPhotonMap()
341 {
342 int i;
343 int err;
344 PHOTON *ph;
345 FILE *f;
346 int numph;
347
348 if (!photonOptions.photonsEnabled) return 0;
349
350 f = fopen(photonOptions.fileName, "rb");
351 if (!f) return 0;
352
353 fread(&numph, sizeof(numph),1,f);
354
355 for(i=0; i<numph; i++)
356 {
357 ph = AllocatePhoton(&photonOptions.photonMap);
358 err = fread(ph, sizeof(PHOTON), 1, f);
359
360 if (err<=0)
361 {
362 /* fread returned an error! */
363 fclose(f);
364 return 0;
365 }
366 }
367
368 if (!feof(f)) /* for backwards file format compatibility */
369 {
370
371 #ifdef GLOBAL_PHOTONS
372 /* global photons */
373 fread(&numph, sizeof(numph),1,f);
374 for(i=0; i<numph; i++)
375 {
376 ph = AllocatePhoton(&photonOptions.globalPhotonMap);
377 err = fread(ph, sizeof(PHOTON), 1, f);
378
379 if (err<=0)
380 {
381 /* fread returned an error! */
382 fclose(f);
383 return 0;
384 }
385 }
386 #endif
387
388 /* media photons */
389 fread(&numph, sizeof(numph),1,f);
390 for(i=0; i<numph; i++)
391 {
392 ph = AllocatePhoton(&photonOptions.mediaPhotonMap);
393 err = fread(ph, sizeof(PHOTON), 1, f);
394
395 if (err<=0)
396 {
397 /* fread returned an error! */
398 fclose(f);
399 return 0;
400 }
401 }
402
403 }
404
405 fclose(f);
406 return true;
407 }
408
409 /*****************************************************************************
410
411 FUNCTION
412
413 FreeBacktraceEverything()
414
415 Preconditions:
416 if photonOptions.photonsEnabled is true, then InitBacktraceEverything()
417 must have been called
418
419 PostConditions:
420 if photonOptions.photonsEnabled is true, then
421 photon memory is freed
422 sets photonOptions.photonsEnabled to false
423 else
424 does nothing
425
426 ******************************************************************************/
427
FreeBacktraceEverything()428 void FreeBacktraceEverything()
429 {
430 if (!InitBacktraceWasCalled) return;
431
432 if (photonOptions.photonsEnabled)
433 {
434 /* free everything that we allocated */
435
436 if(photonOptions.photonGatherList)
437 POV_FREE(photonOptions.photonGatherList);
438 photonOptions.photonGatherList = NULL;
439
440 if(photonOptions.photonDistances)
441 POV_FREE(photonOptions.photonDistances);
442 photonOptions.photonDistances = NULL;
443
444 if (photonOptions.sinTheta)
445 POV_FREE(photonOptions.sinTheta);
446 photonOptions.sinTheta = NULL;
447
448 if (photonOptions.cosTheta)
449 POV_FREE(photonOptions.cosTheta);
450 photonOptions.cosTheta = NULL;
451
452 FreePhotonMemory();
453 photonOptions.photonsEnabled = false;
454 }
455 }
456
457 /*****************************************************************************
458
459 FUNCTION
460
461 AllocatePhoton(PHOTON_MAP *map)
462 allocates a photon
463
464 Photons are allocated in blocks. map->head is a
465 dynamically-created array of these blocks. The blocks themselves
466 are allocated as they are needed.
467
468 Preconditions:
469 InitBacktraceEverything was called
470
471 Postconditions:
472 Marks another photon as allocated (and allocates another block of
473 photons if necessary).
474 Returns a pointer to the new photon.
475 This will be the next available photon in array.
476
477 ******************************************************************************/
478
AllocatePhoton(PHOTON_MAP * map)479 PHOTON* AllocatePhoton(PHOTON_MAP *map)
480 {
481 int i,j,k;
482
483 /* array mapping funciton */
484 /* !!!!!!!!!!! warning
485 This code does the same function as the macro PHOTON_AMF
486 It is done here separatly instead of using the macro for
487 speed reasons (to avoid duplicate operations). If the
488 macro is changed, this MUST ALSO BE CHANGED!
489 */
490 i=(map->numPhotons & PHOTON_BLOCK_MASK);
491 j=(map->numPhotons >> (PHOTON_BLOCK_POWER));
492
493 /* new photon */
494 map->numPhotons++;
495
496 if(j == map->numBlocks)
497 {
498 /* the base array is too small, we need to reallocate it */
499 PHOTON **newMap;
500 newMap = (PHOTON **)POV_MALLOC(sizeof(PHOTON *)*map->numBlocks*2, "photons");
501 map->numBlocks*=2;
502
503 /* copy entries */
504 for(k=0; k<j; k++)
505 newMap[k] = map->head[k];
506
507 /* set new entries to zero */
508 for(k=j; k<map->numBlocks; k++)
509 newMap[k] = NULL;
510
511 /* free old map and put the new map in place */
512 POV_FREE(map->head);
513 map->head = newMap;
514 }
515
516 if(map->head[j] == NULL)
517 /* allocate a new block of photons */
518 map->head[j] = (PHOTON *)POV_MALLOC(sizeof(PHOTON)*PHOTON_BLOCK_SIZE, "photons");
519
520 return &(map->head[j][i]);
521 }
522
523 /*****************************************************************************
524
525 FUNCTION
526
527 InitPhotonMemory()
528
529 Initializes photon memory.
530 Must only be called by InitBacktraceEverything().
531
532 ******************************************************************************/
533
InitPhotonMemory()534 static void InitPhotonMemory()
535 {
536 int k;
537
538 /* allocate the base array */
539 photonOptions.photonMap.numPhotons = 0;
540 photonOptions.photonMap.numBlocks = INITIAL_BASE_ARRAY_SIZE;
541 photonOptions.photonMap.head = (PHOTON_BLOCK *)POV_MALLOC(sizeof(PHOTON_BLOCK *)*INITIAL_BASE_ARRAY_SIZE, "photons");
542
543 /* zero the array */
544 for(k=0; k<photonOptions.photonMap.numBlocks; k++)
545 photonOptions.photonMap.head[k] = NULL;
546
547 #ifdef GLOBAL_PHOTONS
548 /* ------------ global photons ----------------*/
549 /* allocate the base array */
550 photonOptions.globalPhotonMap.numPhotons = 0;
551 photonOptions.globalPhotonMap.numBlocks = INITIAL_BASE_ARRAY_SIZE;
552 photonOptions.globalPhotonMap.head = (PHOTON_BLOCK *)POV_MALLOC(sizeof(PHOTON_BLOCK *)*INITIAL_BASE_ARRAY_SIZE, "photons");
553
554 /* zero the array */
555 for(k=0; k<photonOptions.globalPhotonMap.numBlocks; k++)
556 photonOptions.globalPhotonMap.head[k] = NULL;
557 #endif
558
559 /* ------------ media photons ----------------*/
560 /* allocate the base array */
561 photonOptions.mediaPhotonMap.numPhotons = 0;
562 photonOptions.mediaPhotonMap.numBlocks = INITIAL_BASE_ARRAY_SIZE;
563 photonOptions.mediaPhotonMap.head = (PHOTON_BLOCK *)POV_MALLOC(sizeof(PHOTON_BLOCK *)*INITIAL_BASE_ARRAY_SIZE, "photons");
564
565 /* zero the array */
566 for(k=0; k<photonOptions.mediaPhotonMap.numBlocks; k++)
567 photonOptions.mediaPhotonMap.head[k] = NULL;
568 }
569
570 /*****************************************************************************
571
572 FUNCTION
573
574 FreePhotonMemory()
575
576 Frees all allocated blocks and the base array.
577 Must be called only by FreeBacktraceEverything()
578
579 ******************************************************************************/
580
FreePhotonMemory()581 static void FreePhotonMemory()
582 {
583 int j;
584
585 /* if already freed then stop now */
586 if (photonOptions.photonMap.head==NULL)
587 return;
588
589 /* file name to load or save caustic photon map */
590 if ( photonOptions.fileName )
591 {
592 POV_FREE(photonOptions.fileName);
593 photonOptions.fileName=NULL;
594 }
595
596 /* free all non-NULL arrays */
597 for(j=0; j<photonOptions.photonMap.numBlocks; j++)
598 {
599 if(photonOptions.photonMap.head[j] != NULL)
600 {
601 POV_FREE(photonOptions.photonMap.head[j]);
602 }
603 }
604
605 /* free the base array */
606 POV_FREE(photonOptions.photonMap.head);
607 photonOptions.photonMap.head = NULL;
608
609 #ifdef GLOBAL_PHOTONS
610 /* ---------------- global photons -------------- */
611 /* if already freed then stop now */
612 if (photonOptions.globalPhotonMap.head==NULL)
613 return;
614
615 /* free all non-NULL arrays */
616 for(j=0; j<photonOptions.globalPhotonMap.numBlocks; j++)
617 {
618 if(photonOptions.globalPhotonMap.head[j] != NULL)
619 {
620 POV_FREE(photonOptions.globalPhotonMap.head[j]);
621 }
622 }
623
624 /* free the base array */
625 POV_FREE(photonOptions.globalPhotonMap.head);
626 photonOptions.globalPhotonMap.head = NULL;
627 #endif
628
629 /* ---------------- media photons -------------- */
630 /* if already freed then stop now */
631 if (photonOptions.mediaPhotonMap.head==NULL)
632 return;
633
634 /* free all non-NULL arrays */
635 for(j=0; j<photonOptions.mediaPhotonMap.numBlocks; j++)
636 {
637 if(photonOptions.mediaPhotonMap.head[j] != NULL)
638 {
639 POV_FREE(photonOptions.mediaPhotonMap.head[j]);
640 }
641 }
642
643 /* free the base array */
644 POV_FREE(photonOptions.mediaPhotonMap.head);
645 photonOptions.mediaPhotonMap.head = NULL;
646
647 }
648
649
650 /*****************************************************************************
651 *
652 * FUNCTION
653 *
654 * cubic_spline (copied from point.c for use with light attenuation )
655 *
656 * INPUT
657 *
658 * OUTPUT
659 *
660 * RETURNS
661 *
662 * AUTHOR
663 *
664 * POV-Ray Team
665 *
666 * DESCRIPTION
667 *
668 * Cubic spline that has tangents of slope 0 at x == low and at x == high.
669 * For a given value "pos" between low and high the spline value is returned.
670 *
671 * CHANGES
672 *
673 * -
674 *
675 ******************************************************************************/
676
cubic_spline(DBL low,DBL high,DBL pos)677 static DBL cubic_spline(DBL low, DBL high, DBL pos)
678 {
679 /* Check to see if the position is within the proper boundaries. */
680 if (pos < low)
681 return(0.0);
682 else
683 {
684 if (pos >= high)
685 return(1.0);
686 }
687
688 /* Normalize to the interval [0...1]. */
689 pos = (pos - low) / (high - low);
690
691 /* See where it is on the cubic curve. */
692 return(3 - 2 * pos) * pos * pos;
693 }
694
695 /*****************************************************************************
696
697 FUNCTION
698
699 SearchThroughObjects()
700
701 Searches through 'object' and all siblings and children of 'object' to
702 locate objects with PH_TARGET_FLAG set. This flag means that the object
703 receives photons.
704
705 Preconditions:
706 Photon mapping initialized (InitBacktraceEverything() called)
707 'Object' is a object (with or without siblings)
708 'Light' is a light source in the scene
709
710 Postconditions:
711 Photons may have been shot at the object, its siblings, or its children.
712
713 ******************************************************************************/
714
SearchThroughObjects(OBJECT * Object,LIGHT_SOURCE * Light,bool count)715 static void SearchThroughObjects(OBJECT *Object, LIGHT_SOURCE *Light, bool count)
716 {
717 OBJECT *Sib;
718
719 /* check this object and all siblings */
720 for(Sib = Object; Sib != NULL; Sib = Sib -> Sibling)
721 {
722 if(Test_Flag(Sib, PH_TARGET_FLAG) &&
723 !(Sib->Type & LIGHT_SOURCE_OBJECT))
724 {
725 /* do not shoot photons if global lights are turned off for object */
726 if(!Test_Flag(Sib, NO_GLOBAL_LIGHTS_FLAG))
727 ShootPhotonsAtObject(Sib, Light, count);
728
729 Do_Cooperate(1);
730
731 Check_User_Abort(false);
732 }
733 /* if it has children, check them too */
734 else
735 {
736 #ifdef MOTION_BLUR_PATCH
737 if ((Sib->Type & IS_COMPOUND_OBJECT) && !(Sib->Type & MOTION_BLUR_OBJECT))
738 #else
739 if((Sib->Type & IS_COMPOUND_OBJECT))
740 #endif
741 {
742 SearchThroughObjects(((CSG *)Sib)->Children, Light, count);
743 }
744 }
745 }
746 }
747
748 /*****************************************************************************
749
750 FUNCTION
751
752 ShootPhotonsAtObject()
753
754 Shoots photons from 'Light' to 'Object'
755
756 Preconditions:
757 Photon mapping initialized (InitBacktraceEverything() called)
758 'Object' is a object with PH_TARGET_FLAG set
759 'Light' is a light source in the scene
760
761 Possible future expansion:
762 Object==NULL means create a global photon map.
763
764 Postconditions:
765 Photons may have been shot at the object, depending on a variety
766 of flags
767
768 ******************************************************************************/
769
ShootPhotonsAtObject(OBJECT * Object,LIGHT_SOURCE * Light,int count)770 static void ShootPhotonsAtObject(OBJECT *Object, LIGHT_SOURCE *Light, int count)
771 {
772 RAY Ray; /* ray that we shoot */
773 COLOUR Colour, PhotonColour; /* light color and photon color */
774 int i; /* counter */
775 DBL theta, phi; /* rotation angles */
776 DBL dtheta, dphi; /* deltas for theta and phi */
777 DBL jittheta, jitphi; /* jittered versions of theta and phi */
778 DBL mintheta,maxtheta,minphi,maxphi;
779 /* these are minimum and maximum for theta and
780 phi for the spiral shooting */
781 #ifdef AVOID_MIGHT_BE_USED_UNINITIALIZED_WARNINGS_PATCH
782 DBL dist=0.0; /* distance from light to bounding sphere */
783 DBL rad=0.0; /* radius of bounding sphere */
784 #else
785 DBL dist; /* distance from light to bounding sphere */
786 DBL rad; /* radius of bounding sphere */
787 #endif
788 DBL st,ct; /* cos(theta) & sin(theta) for rotation */
789 DBL Attenuation; /* light attenuation for spotlight */
790 DBL costheta_spot;
791 VECTOR up, left, ctr, toctr, v; /* vectors to determine direction of shot */
792 TRANSFORM Trans; /* transformation for rotation */
793 int mergedFlags=0; /* merged flags to see if we should shoot photons */
794 int notComputed=true; /* have the ray containers been computed for this point yet?*/
795 int hitAtLeastOnce = false; /* have we hit the object at least once - for autostop stuff */
796
797 /* get the light source colour */
798 Assign_Colour(Colour, Light->Colour);
799
800 /* set global variable stuff */
801 photonOptions.Light = Light;
802 photonOptions.photonObject = Object;
803
804 /* first, check on various flags... make sure all is a go for this object */
805 if(Object)
806 {
807 mergedFlags = Light->Flags | photonOptions.photonObject->Flags;
808 photonOptions.lightFlags = Light->Flags;
809 }
810 else
811 {
812 mergedFlags = photonOptions.lightFlags = PH_RFR_ON_FLAG | PH_RFL_ON_FLAG; /*Light->Flags;*/
813 }
814
815 if (!(((mergedFlags & PH_RFR_ON_FLAG) && !(mergedFlags & PH_RFR_OFF_FLAG)) ||
816 ((mergedFlags & PH_RFL_ON_FLAG) && !(mergedFlags & PH_RFL_OFF_FLAG))))
817 /* it is a no-go for this object... bail out now */
818 return;
819
820
821 if(Object)
822 {
823 /* find bounding sphere based on bounding box */
824 ctr[X] = photonOptions.photonObject->BBox.Lower_Left[X] + photonOptions.photonObject->BBox.Lengths[X] / 2.0;
825 ctr[Y] = photonOptions.photonObject->BBox.Lower_Left[Y] + photonOptions.photonObject->BBox.Lengths[Y] / 2.0;
826 ctr[Z] = photonOptions.photonObject->BBox.Lower_Left[Z] + photonOptions.photonObject->BBox.Lengths[Z] / 2.0;
827 VSub(v, ctr,photonOptions.photonObject->BBox.Lower_Left);
828 VLength(rad, v);
829
830 /* find direction from object to bounding sphere */
831 VSub(toctr, ctr, Light->Center);
832 VLength(dist, toctr);
833
834 VNormalizeEq(toctr);
835 if ( fabs(fabs(toctr[Z])- 1.) < .1 ) {
836 /* too close to vertical for comfort, so use cross product with horizon */
837 up[X] = 0.; up[Y] = 1.; up[Z] = 0.;
838 }
839 else
840 {
841 up[X] = 0.; up[Y] = 0.; up[Z] = 1.;
842 }
843
844 /* find "left", which is vector perpendicular to toctr */
845 if(Light->Parallel)
846 {
847 /* for parallel lights, left is actually perpendicular to the direction of the
848 light source */
849 VCross(left, Light->Direction, up); VNormalizeEq(left);
850 }
851 else
852 {
853 VCross(left, toctr, up); VNormalizeEq(left);
854 }
855
856
857 /*
858 light dist ctr
859 * ------------------ +
860 ---___ |
861 ---___ | rad
862 ---_|
863
864 */
865
866 /* calculate the spacial separation (spread) */
867 photonOptions.photonSpread = photonOptions.photonObject->Ph_Density*photonOptions.surfaceSeparation;
868
869 /* if rays aren't parallel, divide by dist so we get separation at a distance of 1 unit */
870 if (!Light->Parallel)
871 {
872 photonOptions.photonSpread /= dist;
873 }
874
875 if (count)
876 {
877 /* try to guess the number of photons */
878 DBL x=rad / (photonOptions.photonObject->Ph_Density*photonOptions.surfaceSeparation);
879 x=x*x*M_PI;
880
881 #if(1)
882 if ( ((mergedFlags & PH_RFR_ON_FLAG) && !(mergedFlags & PH_RFR_OFF_FLAG)) &&
883 ((mergedFlags & PH_RFL_ON_FLAG) && !(mergedFlags & PH_RFL_OFF_FLAG)) )
884 {
885 x *= 1.5; /* assume 2 times as many photons with both reflection & refraction */
886 }
887
888 if ( !Test_Flag(photonOptions.photonObject, PH_IGNORE_PHOTONS_FLAG) )
889 {
890 if ( ((mergedFlags & PH_RFR_ON_FLAG) && !(mergedFlags & PH_RFR_OFF_FLAG)) )
891 {
892 if ( ((mergedFlags & PH_RFL_ON_FLAG) && !(mergedFlags & PH_RFL_OFF_FLAG)) )
893 x *= 3; /* assume 3 times as many photons if ignore_photons not used */
894 else
895 x *= 2; /* assume less for only refraction */
896 }
897 }
898
899 x *= 0.5; /* assume 1/2 of photons hit target object */
900 #endif
901
902 photonCountEstimate += x;
903 return;
904 }
905 }
906 else
907 {
908 #ifdef GLOBAL_PHOTONS
909 /* set up for global photon map */
910 mintheta = 0;
911 maxtheta = M_PI;
912 /* determine photonSpread from a number? */
913 photonOptions.photonSpread = Light->Ph_Density*photonOptions.globalSeparation;
914 Make_Vector(up, 1,0,0);
915 Make_Vector(left, 0,1,0);
916 Make_Vector(toctr, 0,0,1);
917 dist = 1.0;
918
919 if (count)
920 {
921 /* try to guess the number of photons */
922 photonCountEstimate += M_PI*4.0/(photonOptions.photonSpread*photonOptions.photonSpread);
923 return;
924 }
925 #else
926 Error("Internal Error - global photons have been disabled."); /* we should never get here if gobal photons are not enabled */
927 #endif
928 }
929
930 /* adjust spread if we are using an area light */
931 if(Light->Area_Light && Light->Photon_Area_Light)
932 {
933 photonOptions.photonSpread *= sqrt((DBL)(Light->Area_Size1*Light->Area_Size2));
934 }
935
936 /* set the photon density - calculate angular density from spacial */
937 if(Light->Parallel)
938 {
939 /* OK, here we fake things a bit for parallel lights. Theta is not really theta.
940 It really represents the radius... but why re-code an entire loop. For POV 4.0
941 this should be re-written as an abstract class with polymorphism. */
942 dtheta = photonOptions.photonSpread;
943 }
944 else
945 {
946
947 /* calculate delta theta */
948 dtheta = atan(photonOptions.photonSpread);
949 }
950
951 /* if photonSpread <= 0.0, we must return or we'll get stuck in an infinite loop! */
952 if (photonOptions.photonSpread <= 0.0)
953 return;
954
955 mintheta = 0;
956 if (Light->Parallel)
957 {
958 maxtheta = rad;
959 }
960 else if (dist>=rad)
961 {
962 maxtheta = atan(rad/dist);
963 }
964 else
965 {
966 maxtheta = M_PI;
967 if (fabs(dist)<EPSILON)
968 {
969 Make_Vector(up, 1,0,0);
970 Make_Vector(left, 0,1,0);
971 Make_Vector(toctr, 0,0,1);
972 }
973 dist = rad;
974 }
975
976
977 /* ---------------------------------------------
978 main ray-shooting loop
979 --------------------------------------------- */
980 i = 0;
981 notComputed = true;
982 for(theta=mintheta; theta<maxtheta; theta+=dtheta)
983 {
984 photonOptions.hitObject = false;
985
986 if (theta<EPSILON)
987 {
988 dphi=2*M_PI;
989 }
990 else
991 {
992 /* remember that for area lights, "theta" really means "radius" */
993 if (Light->Parallel)
994 {
995 dphi = dtheta / theta;
996 }
997 else
998 {
999 dphi=dtheta/sin(theta);
1000 }
1001 }
1002
1003 minphi = -M_PI + dphi*FRAND()*0.5;
1004 maxphi = M_PI - dphi/2 + (minphi+M_PI);
1005 for(phi=minphi; phi<maxphi; phi+=dphi)
1006 {
1007 int x_samples,y_samples;
1008 int area_x, area_y;
1009 /* ------------------- shoot one photon ------------------ */
1010
1011 /* jitter theta & phi */
1012 jitphi = phi + (dphi)*(FRAND() - 0.5)*1.0*photonOptions.jitter;
1013 jittheta = theta + (dtheta)*(FRAND() - 0.5)*1.0*photonOptions.jitter;
1014
1015 /* actually, shoot multiple samples for area light */
1016 if(Light->Area_Light && Light->Photon_Area_Light && !Light->Parallel)
1017 {
1018 x_samples = Light->Area_Size1;
1019 y_samples = Light->Area_Size2;
1020 }
1021 else
1022 {
1023 x_samples = 1;
1024 y_samples = 1;
1025 }
1026
1027 for(area_x=0; area_x<x_samples; area_x++)
1028 for(area_y=0; area_y<y_samples; area_y++)
1029 {
1030
1031 Assign_Vector(Ray.Initial,Light->Center);
1032
1033 if (Light->Area_Light && !Light->Parallel)
1034 {
1035 /* ------------- area light ----------- */
1036 /* we need to make new up, left, and toctr vectors so we can
1037 do proper rotations of theta and phi about toctr. The
1038 ray's initial point and ending points are both jittered to
1039 produce the area-light effect. */
1040 DBL Jitter_u, Jitter_v, ScaleFactor;
1041 VECTOR NewAxis1, NewAxis2;
1042
1043 /* we must recompute the media containers (new start point) */
1044 notComputed = true;
1045
1046 /*
1047 Jitter_u = (int)(FRAND()*Light->Area_Size1);
1048 Jitter_v = (int)(FRAND()*Light->Area_Size2);
1049 */
1050 Jitter_u = area_x; /*+(0.5*FRAND() - 0.25);*/
1051 Jitter_v = area_y; /*+(0.5*FRAND() - 0.25);*/
1052
1053 if (Light->Area_Size1 > 1 && x_samples>1)
1054 {
1055 ScaleFactor = Jitter_u/(DBL)(Light->Area_Size1 - 1) - 0.5;
1056 VScale (NewAxis1, Light->Axis1, ScaleFactor);
1057 }
1058 else
1059 {
1060 Make_Vector(NewAxis1, 0.0, 0.0, 0.0);
1061 }
1062
1063 if (Light->Area_Size2 > 1 && y_samples>1)
1064 {
1065 ScaleFactor = Jitter_v/(DBL)(Light->Area_Size2 - 1) - 0.5;
1066 VScale (NewAxis2, Light->Axis2, ScaleFactor);
1067 }
1068 else
1069 {
1070 Make_Vector(NewAxis2, 0.0, 0.0, 0.0);
1071 }
1072
1073 /* need a new toctr & left */
1074 VAddEq(Ray.Initial, NewAxis1);
1075 VAddEq(Ray.Initial, NewAxis2);
1076
1077 VSub(toctr, ctr, Ray.Initial);
1078 VLength(dist, toctr);
1079
1080 VNormalizeEq(toctr);
1081 if ( fabs(fabs(toctr[Z])- 1.) < .1 ) {
1082 /* too close to vertical for comfort, so use cross product with horizon */
1083 up[X] = 0.; up[Y] = 1.; up[Z] = 0.;
1084 }
1085 else
1086 {
1087 up[X] = 0.; up[Y] = 0.; up[Z] = 1.;
1088 }
1089 VCross(left, toctr, up); VNormalizeEq(left);
1090
1091 if (fabs(dist)<EPSILON)
1092 {
1093 Make_Vector(up, 1,0,0);
1094 Make_Vector(left, 0,1,0);
1095 Make_Vector(toctr, 0,0,1);
1096 }
1097 }
1098
1099 #ifdef AVOID_MIGHT_BE_USED_UNINITIALIZED_WARNINGS_PATCH
1100 DBL dist_of_initial_from_center=0;
1101 #else
1102 DBL dist_of_initial_from_center;
1103 #endif
1104
1105 if (Light->Parallel)
1106 {
1107 DBL a;
1108 VECTOR v;
1109 /* assign the direction */
1110 Assign_Vector(Ray.Direction,Light->Direction);
1111
1112 /* project ctr onto plane defined by Direction & light location */
1113
1114 VDot(a,Ray.Direction, toctr);
1115 VScale(v,Ray.Direction, -a*dist); /* MAYBE NEEDS TO BE NEGATIVE! */
1116
1117 VAdd(Ray.Initial, ctr, v);
1118
1119 /* move point along "left" distance theta (remember theta means rad) */
1120 VScale(v,left,jittheta);
1121
1122 /* rotate pt around Ray.Direction by phi */
1123 /* use POV funcitons... slower but easy */
1124 Compute_Axis_Rotation_Transform(&Trans,Light->Direction,jitphi);
1125 MTransPoint(v, v, &Trans);
1126
1127 VAddEq(Ray.Initial, v);
1128
1129 // compute the length of "v" if we're going to use it
1130 if (Light->Light_Type == CYLINDER_SOURCE)
1131 {
1132 VECTOR initial_from_center;
1133 VSub(initial_from_center, Ray.Initial, Light->Center);
1134 VLength(dist_of_initial_from_center, initial_from_center);
1135 }
1136 }
1137 else
1138 {
1139 /* rotate toctr by theta around up */
1140 st = sin(jittheta);
1141 ct = cos(jittheta);
1142 /* use fast rotation */
1143 v[X] = -st*left[X] + ct*toctr[X];
1144 v[Y] = -st*left[Y] + ct*toctr[Y];
1145 v[Z] = -st*left[Z] + ct*toctr[Z];
1146
1147 /* then rotate by phi around toctr */
1148 /* use POV funcitons... slower but easy */
1149 Compute_Axis_Rotation_Transform(&Trans,toctr,jitphi);
1150 MTransPoint(Ray.Direction, v, &Trans);
1151 }
1152
1153 /* ------ attenuation for spot/cylinder (copied from point.c) ---- */
1154 Attenuation = 1.0;
1155
1156 /* ---------- spot light --------- */
1157 if (Light->Light_Type == SPOT_SOURCE)
1158 {
1159 VDot(costheta_spot, Ray.Direction, Light->Direction);
1160
1161 if (costheta_spot > 0.0)
1162 {
1163 Attenuation = pow(costheta_spot, Light->Coeff);
1164
1165 if (Light->Radius > 0.0)
1166 Attenuation *= cubic_spline(Light->Falloff, Light->Radius, costheta_spot);
1167
1168 }
1169 else
1170 Attenuation = 0.0;
1171 }
1172 /* ---------- cylinder light ----------- */
1173 else if (Light->Light_Type == CYLINDER_SOURCE)
1174 {
1175 DBL k, len;
1176
1177 VDot(k, Ray.Direction, Light->Direction);
1178
1179 if (k > 0.0)
1180 {
1181 len = dist_of_initial_from_center;
1182
1183 if (len < Light->Falloff)
1184 {
1185 DBL dist = 1.0 - len / Light->Falloff;
1186 Attenuation = pow(dist, Light->Coeff);
1187
1188 if (Light->Radius > 0.0 && len > Light->Radius)
1189 Attenuation *= cubic_spline(0.0, 1.0 - Light->Radius / Light->Falloff, dist);
1190
1191 }
1192 else
1193 Attenuation = 0.0;
1194 }
1195 else
1196 Attenuation = 0.0;
1197 }
1198
1199 /* set up defaults for reflection, refraction */
1200 photonOptions.passThruPrev = true;
1201 photonOptions.passThruThis = false;
1202
1203 photonOptions.photonDepth = 0.0;
1204 Trace_Level = 1;
1205 Total_Depth = 0.0;
1206 Increase_Counter(stats[Number_Of_Photons_Shot]);
1207
1208 /* attenuate for area light extra samples */
1209 Attenuation/=(x_samples*y_samples);
1210
1211 /* compute photon color from light source & attenuation */
1212
1213 PhotonColour[0] = Colour[0]*Attenuation;
1214 PhotonColour[1] = Colour[1]*Attenuation;
1215 PhotonColour[2] = Colour[2]*Attenuation;
1216 PhotonColour[3] = 0.0;
1217 PhotonColour[4] = 0.0;
1218
1219 if (Attenuation<0.00001) continue;
1220
1221 /* handle the projected_through object if it exists */
1222 if (Light->Projected_Through_Object != NULL)
1223 {
1224 /* try to intersect ray with projected-through object */
1225 INTERSECTION Intersect;
1226
1227 Intersect.Object = NULL;
1228 if ( Intersection( &Intersect, Light->Projected_Through_Object, &Ray ) )
1229 {
1230 /* we must recompute the media containers (new start point) */
1231 notComputed = true;
1232
1233 /* we did hit it, so find the 'real' starting point of the ray */
1234 /* find the farthest intersection */
1235 VAddScaledEq(Ray.Initial,Intersect.Depth+EPSILON, Ray.Direction);
1236 photonOptions.photonDepth += Intersect.Depth+EPSILON;
1237 while(Intersection( &Intersect, Light->Projected_Through_Object, &Ray ) )
1238 {
1239 VAddScaledEq(Ray.Initial, Intersect.Depth+EPSILON, Ray.Direction);
1240 photonOptions.photonDepth += Intersect.Depth+EPSILON;
1241 }
1242 }
1243 else
1244 {
1245 /* we didn't hit it, so stop now */
1246 continue;
1247 }
1248
1249 }
1250
1251 /* As mike said, "fire photon torpedo!" */
1252 Initialize_Ray_Containers(&Ray);
1253 initialize_ray_container_state(&Ray, notComputed);
1254 notComputed = false;
1255 disp_elem = 0; /* for dispersion */
1256 disp_nelems = 0; /* for dispersion */
1257
1258 #ifdef MOTION_BLUR_PATCH
1259 TimeStamp = (opts.motionBlurCount+1)/2; /* must be non-zero */
1260 #endif
1261 Trace(&Ray, PhotonColour, 1.0);
1262
1263 /* display here */
1264 i++;
1265 if ((i%100) == 0)
1266 {
1267 gPhotonStat_i = i;
1268 gPhotonStat_x_samples = x_samples;
1269 gPhotonStat_y_samples = y_samples;
1270 Send_ProgressUpdate(PROGRESS_BUILDING_PHOTON_MAPS);
1271 #ifdef FASTER_PROGRESS_DISPLAY_PATCH
1272 Do_Cooperate(0);
1273 #else
1274 Check_User_Abort(false);
1275 #endif
1276 }
1277
1278 } /* end of multiple samples */
1279 }
1280
1281 /* if we didn't hit anything and we're past the autostop angle, then
1282 we should stop
1283
1284 as per suggestion from Smellenberg, changed autostop to a percentage
1285 of the object's bounding sphere. */
1286
1287 /* suggested by Pabs, we only use autostop if we have it it once */
1288 if (photonOptions.hitObject) hitAtLeastOnce=true;
1289
1290 if (hitAtLeastOnce && !photonOptions.hitObject && photonOptions.photonObject)
1291 if (theta>photonOptions.autoStopPercent*maxtheta) break;
1292 } /* end of rays loop */
1293 }
1294
1295 /*****************************************************************************
1296
1297 FUNCTION
1298
1299 BuildPhotonMaps()
1300
1301 This is the primary function for building photon maps.
1302
1303 Preconditions:
1304 Photon memory is allocated (InitBacktraceEverything has been called).
1305 The entire scene has been parsed.
1306 The ray-tracer has been completely initialized.
1307
1308 Postconditions:
1309 The photon map is built based on options specified in the scene file.
1310
1311 ******************************************************************************/
1312
BuildPhotonMaps(void)1313 void BuildPhotonMaps(void)
1314 {
1315 LIGHT_SOURCE *Light; /* light source to use */
1316 LIGHT_GROUP_LIGHT *Light_Group_Light;
1317 int old_mtl; /* saved max_trace_level */
1318 DBL old_adc; /* saved adc_bailout */
1319
1320 /* if not enabled, then return */
1321 if (!photonOptions.photonsEnabled)
1322 return;
1323
1324 /* should we load the photon map instead of building? */
1325 if (photonOptions.fileName && photonOptions.loadFile)
1326 {
1327 /* status bar for user */
1328 Send_Progress("Loading Photon Maps", PROGRESS_LOADING_PHOTON_MAPS);
1329 if (!loadPhotonMap())
1330 {
1331 Error("Could not load photon map (%s)",photonOptions.fileName);
1332 }
1333
1334 Do_Cooperate(0);
1335
1336 /* don't build */
1337 /* but do set photon options automatically */
1338
1339 /* ----------- surface photons ------------- */
1340 if (photonOptions.photonMap.numPhotons>0)
1341 {
1342 setGatherOptions(&photonOptions.photonMap, false);
1343 }
1344
1345 #ifdef GLOBAL_PHOTONS
1346 /* ----------- global photons ------------- */
1347 if (photonOptions.globalPhotonMap.numPhotons>0)
1348 {
1349 setGatherOptions(&photonOptions.globalPhotonMap, false);
1350 }
1351 #endif
1352
1353 /* ----------- media photons ------------- */
1354 if (photonOptions.mediaPhotonMap.numPhotons>0)
1355 {
1356 setGatherOptions(&photonOptions.mediaPhotonMap, true);
1357 }
1358
1359 return;
1360 }
1361
1362 /* set flag so POV knows we're in backtrace step */
1363 backtraceFlag = 1;
1364
1365 /* status bar for user */
1366 Send_Progress("Building Photon Maps", PROGRESS_BUILDING_PHOTON_MAPS);
1367
1368 /* save adc_bailout and max_trace_level */
1369 old_adc = ADC_Bailout;
1370 old_mtl = Max_Trace_Level;
1371
1372 /* use the photon-specific ones if they were specified */
1373 if (photonOptions.Max_Trace_Level>=0)
1374 Max_Trace_Level = photonOptions.Max_Trace_Level;
1375
1376 if (photonOptions.ADC_Bailout>=0)
1377 ADC_Bailout = photonOptions.ADC_Bailout;
1378
1379 /* COUNT THE PHOTONS */
1380 if(photonOptions.surfaceCount>0)
1381 {
1382 DBL factor;
1383 photonCountEstimate = 0.0;
1384
1385 // global lights
1386 photonOptions.Light_Is_Global = true;
1387 for (Light = Frame.Light_Sources;
1388 Light != NULL;
1389 Light = Light->Next_Light_Source)
1390 if (Light->Light_Type != FILL_LIGHT_SOURCE)
1391 {
1392 SearchThroughObjects(Frame.Objects, Light, true);
1393 }
1394
1395 // light_group lights
1396 photonOptions.Light_Is_Global = false;
1397 for (Light_Group_Light = Frame.Light_Group_Lights;
1398 Light_Group_Light != NULL;
1399 Light_Group_Light = Light_Group_Light->Next)
1400 {
1401 Light = Light_Group_Light->Light;
1402 if (Light->Light_Type != FILL_LIGHT_SOURCE)
1403 {
1404 SearchThroughObjects(Frame.Objects, Light, true);
1405 }
1406 }
1407
1408 factor = (DBL)photonCountEstimate/photonOptions.surfaceCount;
1409 factor = sqrt(factor);
1410 photonOptions.surfaceSeparation *= factor;
1411 }
1412
1413 /* COUNT THE GLOBAL PHOTONS */
1414 if(photonOptions.globalCount>0)
1415 {
1416 DBL factor;
1417 photonCountEstimate = 0.0;
1418
1419 photonOptions.Light_Is_Global = true;
1420 for (Light = Frame.Light_Sources;
1421 Light != NULL;
1422 Light = Light->Next_Light_Source)
1423 if (Light->Light_Type != FILL_LIGHT_SOURCE)
1424 {
1425 ShootPhotonsAtObject(NULL, Light, true);
1426 }
1427
1428 // light_group lights
1429 photonOptions.Light_Is_Global = false;
1430 for (Light_Group_Light = Frame.Light_Group_Lights;
1431 Light_Group_Light != NULL;
1432 Light_Group_Light = Light_Group_Light->Next)
1433 {
1434 Light = Light_Group_Light->Light;
1435 if (Light->Light_Type != FILL_LIGHT_SOURCE)
1436 {
1437 ShootPhotonsAtObject(NULL, Light, true);
1438 }
1439 }
1440
1441 factor = (DBL)photonCountEstimate/photonOptions.globalCount;
1442 factor = sqrt(factor);
1443 photonOptions.globalSeparation *= factor;
1444
1445 Do_Cooperate(1);
1446 }
1447
1448 // there is a world out there that wants some attention [trf]
1449 Do_Cooperate(0);
1450
1451 /* loop through global light sources */
1452 photonOptions.Light_Is_Global = true;
1453 for (Light = Frame.Light_Sources;
1454 Light != NULL;
1455 Light = Light->Next_Light_Source)
1456 if (Light->Light_Type != FILL_LIGHT_SOURCE)
1457 {
1458 if (Light->Light_Type == CYLINDER_SOURCE && !Light->Parallel)
1459 {
1460 Warning(0,"Cylinder lights should be parallel when used with photons.");
1461 }
1462
1463 /* do global lighting here if it is ever implemented */
1464 if (Test_Flag(Light, PH_TARGET_FLAG) && (photonOptions.globalCount>0))
1465 ShootPhotonsAtObject(NULL, Light, false);
1466
1467 /* do object-specific lighting */
1468 SearchThroughObjects(Frame.Objects, Light, false);
1469 }
1470
1471 // loop through light_group light sources
1472 photonOptions.Light_Is_Global = false;
1473 for (Light_Group_Light = Frame.Light_Group_Lights;
1474 Light_Group_Light != NULL;
1475 Light_Group_Light = Light_Group_Light->Next)
1476 {
1477 Light = Light_Group_Light->Light;
1478
1479 if (Light->Light_Type == CYLINDER_SOURCE && !Light->Parallel)
1480 {
1481 Warning(0,"Cylinder lights should be parallel when used with photons.");
1482 }
1483
1484 /* do global lighting here if it is ever implemented */
1485 if (Test_Flag(Light, PH_TARGET_FLAG) && (photonOptions.globalCount>0))
1486 ShootPhotonsAtObject(NULL, Light, false);
1487
1488 /* do object-specific lighting */
1489 SearchThroughObjects(Frame.Objects, Light, false);
1490 }
1491
1492 /* clear this flag */
1493 backtraceFlag = 0;
1494
1495 /* restore saved variables */
1496 ADC_Bailout = old_adc;
1497 Max_Trace_Level = old_mtl;
1498
1499 /* now actually build the kd-tree by sorting the array of photons */
1500 if (photonOptions.photonMap.numPhotons>0)
1501 {
1502 buildTree(&photonOptions.photonMap);
1503 setGatherOptions(&photonOptions.photonMap, false);
1504 }
1505
1506 #ifdef GLOBAL_PHOTONS
1507 /* ----------- global photons ------------- */
1508 if (photonOptions.globalPhotonMap.numPhotons>0)
1509 {
1510 buildTree(&photonOptions.globalPhotonMap);
1511 setGatherOptions(&photonOptions.globalPhotonMap, false);
1512 }
1513 #endif
1514
1515 /* ----------- media photons ------------- */
1516 if (photonOptions.mediaPhotonMap.numPhotons>0)
1517 {
1518 buildTree(&photonOptions.mediaPhotonMap);
1519 setGatherOptions(&photonOptions.mediaPhotonMap, true);
1520 }
1521
1522 if (photonOptions.photonMap.numPhotons+
1523 #ifdef GLOBAL_PHOTONS
1524 photonOptions.globalPhotonMap.numPhotons+
1525 #endif
1526 photonOptions.mediaPhotonMap.numPhotons > 0)
1527 {
1528 /* should we load the photon map now that it is built? */
1529 if (photonOptions.fileName && !photonOptions.loadFile)
1530 {
1531 /* status bar for user */
1532 Send_Progress("Saving Photon Maps", PROGRESS_SAVING_PHOTON_MAPS);
1533 if (!savePhotonMap())
1534 {
1535 Warning(0,"Could not save photon map.");
1536 }
1537 }
1538 }
1539 else
1540 {
1541 if (photonOptions.fileName && !photonOptions.loadFile)
1542 {
1543 Warning(0,"Could not save photon map - no photons!");
1544 }
1545 }
1546
1547 // good idea to make sure all warnings and errors arrive frontend now [trf]
1548 Do_Cooperate(0);
1549 }
1550
1551 /*****************************************************************************
1552
1553 FUNCTION
1554
1555 addSurfacePhoton()
1556
1557 Adds a photon to the array of photons.
1558
1559 Preconditions:
1560 InitBacktraceEverything() was called
1561 'Point' is the intersection point to store the photon
1562 'Origin' is the origin of the light ray
1563 'LightCol' is the color of the light propogated through the scene
1564 'RawNorm' is the raw normal of the surface intersected
1565
1566 Postconditions:
1567 Another photon is allocated (by AllocatePhoton())
1568 The information passed in (as well as photonOptions.photonDepth)
1569 is stored in the photon data structure.
1570
1571 ******************************************************************************/
1572
addSurfacePhoton(VECTOR Point,VECTOR Origin,COLOUR LightCol,VECTOR)1573 void addSurfacePhoton(VECTOR Point, VECTOR Origin, COLOUR LightCol, VECTOR /*RawNorm*/)
1574 {
1575 PHOTON *Photon;
1576 COLOUR LightCol2;
1577 DBL Attenuation;
1578 VECTOR d;
1579 DBL d_len, phi, theta;
1580 PHOTON_MAP *map;
1581
1582 /* first, compensate for POV's weird light attenuation */
1583 if ((photonOptions.Light->Fade_Power > 0.0) && (fabs(photonOptions.Light->Fade_Distance) > EPSILON))
1584 {
1585 Attenuation = 2.0 / (1.0 + pow(photonOptions.photonDepth / photonOptions.Light->Fade_Distance, photonOptions.Light->Fade_Power));
1586 }
1587 else
1588 Attenuation = 1;
1589
1590 VScale(LightCol2, LightCol, Attenuation);
1591
1592 if(!photonOptions.Light->Parallel)
1593 {
1594 VScaleEq(LightCol2, photonOptions.photonDepth*photonOptions.photonDepth);
1595 }
1596
1597 VScaleEq(LightCol2, photonOptions.photonSpread*photonOptions.photonSpread);
1598
1599 /* if too dark, maybe we should stop here */
1600
1601 #ifdef GLOBAL_PHOTONS
1602 if(photonOptions.photonObject==NULL)
1603 {
1604 map = &photonOptions.globalPhotonMap;
1605 Increase_Counter(stats[Number_Of_Global_Photons_Stored]);
1606 }
1607 else
1608 #endif
1609 {
1610 map = &photonOptions.photonMap;
1611 Increase_Counter(stats[Number_Of_Photons_Stored]);
1612 }
1613
1614
1615 /* allocate the photon */
1616 Photon = AllocatePhoton(map);
1617
1618 /* convert photon from three floats to 4 bytes */
1619 colour2photonRgbe(Photon->Colour, LightCol2);
1620
1621 /* store the location */
1622 Assign_Vector(Photon->Loc, Point);
1623
1624 /* now determine rotation angles */
1625 VSub(d,Origin, Point);
1626 VNormalizeEq(d);
1627 d_len = sqrt(d[X]*d[X]+d[Z]*d[Z]);
1628
1629 phi = acos(d[X]/d_len);
1630 if (d[Z]<0) phi = -phi;
1631
1632 theta = acos(d_len);
1633 if (d[Y]<0) theta = -theta;
1634
1635 /* cram these rotation angles into two signed bytes */
1636 Photon->theta=(signed char)(theta*127.0/M_PI);
1637 Photon->phi=(signed char)(phi*127.0/M_PI);
1638
1639 }
1640
1641 /*****************************************************************************
1642
1643 FUNCTION
1644
1645 addMediaPhoton()
1646
1647 Adds a photon to the array of photons.
1648
1649 Preconditions:
1650 InitBacktraceEverything() was called
1651 'Point' is the intersection point to store the photon
1652 'Origin' is the origin of the light ray
1653 'LightCol' is the color of the light propogated through the scene
1654
1655 Postconditions:
1656 Another photon is allocated (by AllocatePhoton())
1657 The information passed in (as well as photonOptions.photonDepth)
1658 is stored in the photon data structure.
1659
1660 ******************************************************************************/
1661
addMediaPhoton(VECTOR Point,VECTOR Origin,COLOUR LightCol,DBL depthDiff)1662 void addMediaPhoton(VECTOR Point, VECTOR Origin, COLOUR LightCol, DBL depthDiff)
1663 {
1664 PHOTON *Photon;
1665 COLOUR LightCol2;
1666 DBL Attenuation;
1667 VECTOR d;
1668 DBL d_len, phi, theta;
1669
1670 /* first, compensate for POV's weird light attenuation */
1671 if ((photonOptions.Light->Fade_Power > 0.0) && (fabs(photonOptions.Light->Fade_Distance) > EPSILON))
1672 {
1673 Attenuation = 2.0 / (1.0 + pow((photonOptions.photonDepth+depthDiff) / photonOptions.Light->Fade_Distance, photonOptions.Light->Fade_Power));
1674 }
1675 else
1676 Attenuation = 1;
1677
1678 #if 0
1679 VScale(LightCol2, LightCol, photonOptions.photonSpread*photonOptions.photonSpread);
1680 if(!photonOptions.Light->Parallel)
1681 {
1682 VScaleEq(LightCol2, (photonOptions.photonDepth+depthDiff)*(photonOptions.photonDepth+depthDiff)*Attenuation);
1683 }
1684 #else
1685 VScale(LightCol2, LightCol, Attenuation);
1686 if(!photonOptions.Light->Parallel)
1687 {
1688 VScaleEq(LightCol2, (photonOptions.photonDepth+depthDiff) *
1689 (photonOptions.photonDepth+depthDiff));
1690 }
1691 VScaleEq(LightCol2, photonOptions.photonSpread*photonOptions.photonSpread);
1692 #endif
1693
1694 /* if too dark, maybe we should stop here */
1695
1696 /* allocate the photon */
1697 if(photonOptions.photonObject==NULL) return;
1698
1699 Increase_Counter(stats[Number_Of_Media_Photons_Stored]);
1700
1701 Photon = AllocatePhoton(&photonOptions.mediaPhotonMap);
1702
1703 /* convert photon from three floats to 4 bytes */
1704 colour2photonRgbe(Photon->Colour, LightCol2);
1705
1706 /* store the location */
1707 Assign_Vector(Photon->Loc, Point);
1708
1709 /* now determine rotation angles */
1710 VSub(d,Origin, Point);
1711 VNormalizeEq(d);
1712 d_len = sqrt(d[X]*d[X]+d[Z]*d[Z]);
1713
1714 phi = acos(d[X]/d_len);
1715 if (d[Z]<0) phi = -phi;
1716
1717 theta = acos(d_len);
1718 if (d[Y]<0) theta = -theta;
1719
1720 /* cram these rotation angles into two signed bytes */
1721 Photon->theta=(signed char)(theta*127.0/M_PI);
1722 Photon->phi=(signed char)(phi*127.0/M_PI);
1723
1724 }
1725
1726 /* ====================================================================== */
1727 /* ====================================================================== */
1728 /* KD - TREE */
1729 /* ====================================================================== */
1730 /* ====================================================================== */
1731
1732 /*****************************************************************************
1733
1734 FUNCTION
1735
1736 swapPhotons
1737
1738 swaps two photons
1739
1740 Precondition:
1741 photon memory initialized
1742 static map_s points to the photon map
1743 'a' and 'b' are indexes within the range of photons in map_s
1744 (NO ERROR CHECKING IS DONE)
1745
1746 Postconditions:
1747 the photons indexed by 'a' and 'b' are swapped
1748
1749 *****************************************************************************/
1750
swapPhotons(int a,int b)1751 static void swapPhotons(int a, int b)
1752 {
1753 int ai,aj,bi,bj;
1754 PHOTON tmp;
1755
1756 /* !!!!!!!!!!! warning
1757 This code does the same function as the macro PHOTON_AMF
1758 It is done here separatly instead of using the macro for
1759 speed reasons (to avoid duplicate operations). If the
1760 macro is changed, this MUST ALSO BE CHANGED!
1761 */
1762 ai = a & PHOTON_BLOCK_MASK;
1763 aj = a >> PHOTON_BLOCK_POWER;
1764 bi = b & PHOTON_BLOCK_MASK;
1765 bj = b >> PHOTON_BLOCK_POWER;
1766
1767 tmp = map_s[aj][ai];
1768 map_s[aj][ai] = map_s[bj][bi];
1769 map_s[bj][bi] = tmp;
1770 }
1771
1772 /*****************************************************************************
1773
1774 FUNCTION
1775
1776 insertSort
1777 (modified from Data Structures textbook)
1778
1779 Preconditions:
1780 photon memory initialized
1781 static map_s points to the photon map
1782 'start' is the index of the first photon
1783 'end' is the index of the last photon
1784 'd' is the dimension to sort on (X, Y, or Z)
1785
1786 Postconditions:
1787 photons from 'start' to 'end' in map_s are sorted in
1788 ascending order on dimension d
1789 ******************************************************************************/
insertSort(int start,int end,int d)1790 static void insertSort(int start, int end, int d)
1791 {
1792 int j,k;
1793 PHOTON tmp;
1794
1795 for(k=end-1; k>=start; k--)
1796 {
1797 j=k+1;
1798 tmp = PHOTON_AMF(map_s, k);
1799 while ( (tmp.Loc[d] > PHOTON_AMF(map_s,j).Loc[d]) )
1800 {
1801 PHOTON_AMF(map_s,j-1) = PHOTON_AMF(map_s,j);
1802 j++;
1803 if (j>end) break;
1804 }
1805 PHOTON_AMF(map_s,j-1) = tmp;
1806 }
1807 }
1808
1809 /*****************************************************************************
1810
1811 FUNCTION
1812
1813 quickSortRec
1814 (modified from Data Structures textbook)
1815
1816 Recursive part of the quicksort routine
1817 This does not sort all the way. once this is done, insertSort
1818 should be called to finish the sorting process!
1819
1820 Preconditions:
1821 photon memory initialized
1822 static map_s points to the photon map
1823 'left' is the index of the first photon
1824 'right' is the index of the last photon
1825 'd' is the dimension to sort on (X, Y, or Z)
1826
1827 Postconditions:
1828 photons from 'left' to 'right' in map_s are MOSTLY sorted in
1829 ascending order on dimension d
1830 ******************************************************************************/
quickSortRec(int left,int right,int d)1831 static void quickSortRec(int left, int right, int d)
1832 {
1833 int j,k;
1834 if(left<right)
1835 {
1836 swapPhotons(((left+right)>>1), left+1);
1837 if(PHOTON_AMF(map_s,left+1).Loc[d] > PHOTON_AMF(map_s,right).Loc[d])
1838 swapPhotons(left+1,right);
1839 if(PHOTON_AMF(map_s,left).Loc[d] > PHOTON_AMF(map_s,right).Loc[d])
1840 swapPhotons(left,right);
1841 if(PHOTON_AMF(map_s,left+1).Loc[d] > PHOTON_AMF(map_s,left).Loc[d])
1842 swapPhotons(left+1,left);
1843
1844 j=left+1; k=right;
1845 while(j<=k)
1846 {
1847 for(j++; ((j<=right)&&(PHOTON_AMF(map_s,j).Loc[d]<PHOTON_AMF(map_s,left).Loc[d])); j++)
1848 #ifdef WHITE_SPACE_BEFORE_SEMICOLON_PATCH
1849 // Some compilers needs white space between for(), if(), while() and ';'
1850 // Otherwise mistake can be reported as warning
1851 #endif
1852 ;
1853 for(k--; ((k>=left)&&(PHOTON_AMF(map_s,k).Loc[d]>PHOTON_AMF(map_s,left).Loc[d])); k--)
1854 #ifdef WHITE_SPACE_BEFORE_SEMICOLON_PATCH
1855 // Some compilers needs white space between for(), if(), while() and ';'
1856 // Otherwise mistake can be reported as warning
1857 #endif
1858 ;
1859
1860 if(j<k)
1861 swapPhotons(j,k);
1862 }
1863
1864 swapPhotons(left,k);
1865 if(k-left > 10)
1866 {
1867 quickSortRec(left,k-1,d);
1868 }
1869 if(right-k > 10)
1870 {
1871 quickSortRec(k+1,right,d);
1872 }
1873 /* leave the rest for insertSort */
1874 }
1875 }
1876
1877 /*****************************************************************************
1878
1879 FUNCTION
1880
1881 halfSortRec
1882 (modified quicksort algorithm)
1883
1884 Recursive part of the quicksort routine, but it only does half
1885 the quicksort. It only recurses one branch - the branch that contains
1886 the midpoint (median).
1887
1888 Preconditions:
1889 photon memory initialized
1890 static map_s points to the photon map
1891 'left' is the index of the first photon
1892 'right' is the index of the last photon
1893 'd' is the dimension to sort on (X, Y, or Z)
1894 'mid' is the index where the median will end up
1895
1896 Postconditions:
1897 the photon at the midpoint (mid) is the median of the photons
1898 when sorted on dimention d.
1899 ******************************************************************************/
halfSortRec(int left,int right,int d,int mid)1900 static void halfSortRec(int left, int right, int d, int mid)
1901 {
1902 int j,k;
1903 if(left<right)
1904 {
1905 swapPhotons(((left+right)>>1), left+1);
1906 if(PHOTON_AMF(map_s,left+1).Loc[d] > PHOTON_AMF(map_s,right).Loc[d])
1907 swapPhotons(left+1,right);
1908 if(PHOTON_AMF(map_s,left).Loc[d] > PHOTON_AMF(map_s,right).Loc[d])
1909 swapPhotons(left,right);
1910 if(PHOTON_AMF(map_s,left+1).Loc[d] > PHOTON_AMF(map_s,left).Loc[d])
1911 swapPhotons(left+1,left);
1912
1913 j=left+1; k=right;
1914 while(j<=k)
1915 {
1916 for(j++; ((j<=right)&&(PHOTON_AMF(map_s,j).Loc[d]<PHOTON_AMF(map_s,left).Loc[d])); j++)
1917 #ifdef WHITE_SPACE_BEFORE_SEMICOLON_PATCH
1918 // Some compilers needs white space between for(), if(), while() and ';'
1919 // Otherwise mistake can be reported as warning
1920 #endif
1921 ;
1922 for(k--; ((k>=left)&&(PHOTON_AMF(map_s,k).Loc[d]>PHOTON_AMF(map_s,left).Loc[d])); k--)
1923 #ifdef WHITE_SPACE_BEFORE_SEMICOLON_PATCH
1924 // Some compilers needs white space between for(), if(), while() and ';'
1925 // Otherwise mistake can be reported as warning
1926 #endif
1927 ;
1928
1929 if(j<k)
1930 swapPhotons(j,k);
1931 }
1932
1933 /* put the pivot into its position */
1934 swapPhotons(left,k);
1935
1936 /* only go down the side that contains the midpoint */
1937 /* don't do anything if the midpoint=k (the pivot, which is
1938 now in the correct position */
1939 if(k-left > 0 && (mid>=left) && (mid<k))
1940 {
1941 halfSortRec(left,k-1,d,mid);
1942 }
1943 else if(right-k > 0 && (mid>k) && (mid<=right))
1944 {
1945 halfSortRec(k+1,right,d,mid);
1946 }
1947 }
1948 }
1949
1950 /*****************************************************************************
1951
1952 FUNCTION
1953
1954 sortAndSubdivide
1955
1956 Finds the dimension with the greates range, sorts the photons on that
1957 dimension. Then it recurses on the left and right halves (keeping
1958 the median photon as a pivot). This produces a balanced kd-tree.
1959
1960 Preconditions:
1961 photon memory initialized
1962 static map_s points to the photon map
1963 'start' is the index of the first photon
1964 'end' is the index of the last photon
1965 'sorted' is the dimension that was last sorted (so we don't sort again)
1966
1967 Postconditions:
1968 photons from 'start' to 'end' in map_s are in a valid kd-tree format
1969 ******************************************************************************/
sortAndSubdivide(int start,int end,int)1970 static void sortAndSubdivide(int start, int end, int /*sorted*/)
1971 {
1972 int i,j; /* counters */
1973 SNGL_VECT min,max; /* min/max vectors for finding range */
1974 int DimToUse; /* which dimesion has the greatest range */
1975 int mid; /* index of median (middle) */
1976 int len; /* length of the array we're sorting */
1977
1978 if (end==start)
1979 {
1980 PHOTON_AMF(map_s, start).info = 0;
1981 return;
1982 }
1983
1984 if(end<start) return;
1985
1986 // this seems to be a good place for this call [trf]
1987 Do_Cooperate(1);
1988
1989 /*
1990 * loop and find greatest range
1991 */
1992 Make_Vector(min, 1/EPSILON, 1/EPSILON, 1/EPSILON);
1993 Make_Vector(max, -1/EPSILON, -1/EPSILON, -1/EPSILON);
1994
1995 for(i=start; i<=end; i++)
1996 {
1997 for(j=X; j<=Z; j++)
1998 {
1999 PHOTON *ph = &(PHOTON_AMF(map_s,i));
2000
2001 if (ph->Loc[j] < min[j])
2002 min[j]=ph->Loc[j];
2003 if (ph->Loc[j] > max[j])
2004 max[j]=ph->Loc[j];
2005 }
2006 }
2007
2008 /* choose which dimension to use */
2009 DimToUse = X;
2010 if((max[Y]-min[Y])>(max[DimToUse]-min[DimToUse]))
2011 DimToUse=Y;
2012 if((max[Z]-min[Z])>(max[DimToUse]-min[DimToUse]))
2013 DimToUse=Z;
2014
2015 /* find midpoint */
2016 mid = (end+start)>>1;
2017
2018 /* use half of a quicksort to find the median */
2019 len = end-start;
2020 if (len>=2)
2021 {
2022 /* only display status every so often */
2023 if(len > 1000)
2024 {
2025 gPhotonStat_end = end;
2026 Send_ProgressUpdate(PROGRESS_SORTING_PHOTONS);
2027 #ifdef FASTER_PROGRESS_DISPLAY_PATCH
2028 Do_Cooperate(0);
2029 #endif
2030 }
2031
2032 halfSortRec(start, end, DimToUse, mid);
2033 //quickSortRec(start, end, DimToUse);
2034 }
2035
2036 /* set DimToUse for the midpoint */
2037 PHOTON_AMF(map_s, mid).info = DimToUse;
2038
2039 /* now recurse to continue building the kd-tree */
2040 sortAndSubdivide(start, mid - 1, DimToUse);
2041 sortAndSubdivide(mid + 1, end, DimToUse);
2042 }
2043
2044 /*****************************************************************************
2045
2046 FUNCTION
2047
2048 buildTree
2049
2050 Builds the kd-tree by calling sortAndSubdivide(). Sets the static
2051 variable map_s.
2052
2053 Preconditions:
2054 photon memory initialized
2055 'map' is a pointer to a photon map containing an array of unsorted
2056 photons
2057
2058 Postconditions:
2059 photons are in a valid kd-tree format
2060 ******************************************************************************/
buildTree(PHOTON_MAP * map)2061 static void buildTree(PHOTON_MAP *map)
2062 {
2063 Send_Progress("Sorting photons", PROGRESS_SORTING_PHOTONS);
2064 map_s = map->head;
2065 sortAndSubdivide(0,map->numPhotons-1,X+Y+Z/*this is not X, Y, or Z*/);
2066 }
2067
2068 /*****************************************************************************
2069
2070 FUNCTION
2071
2072 setGatherOptions
2073
2074 determines gather options
2075
2076 Preconditions:
2077 photon memory initialized
2078 'map' points to an already-built (and sorted) photon map
2079 'mediaMap' is true if 'map' contians media photons, and false if
2080 'map' contains surface photons
2081
2082 Postconditions:
2083 gather gather options are set for this map
2084 ******************************************************************************/
setGatherOptions(PHOTON_MAP * map,int mediaMap)2085 static void setGatherOptions(PHOTON_MAP *map, int mediaMap)
2086 {
2087 DBL r;
2088 DBL density;
2089 VECTOR Point;
2090 int numToSample;
2091 int n,i,j;
2092 DBL mind,maxd,avgd;
2093 DBL sum,sum2;
2094 DBL saveDensity;
2095 //int greaterThan;
2096 int lessThan;
2097
2098 /* if user did not set minimum gather radius,
2099 then we should calculate it */
2100 if (map->minGatherRad <= 0.0)
2101 {
2102 mind=10000000.0;
2103 maxd=avgd=sum=sum2=0.0;
2104
2105 /* use 5% of photons, min 100, max 10000 */
2106 numToSample = map->numPhotons/20;
2107 if (numToSample>1000) numToSample = 1000;
2108 if (numToSample<100) numToSample = 100;
2109
2110 for(i=0; i<numToSample; i++)
2111 {
2112 j = rand() % map->numPhotons;
2113
2114 Assign_Vector(Point,(PHOTON_AMF(map->head, j)).Loc);
2115
2116 n=gatherPhotons(Point, 10000000.0, &r, NULL, false, map);
2117
2118 if(mediaMap)
2119 density = 3.0 * n / (4.0*M_PI*r*r*r); /* should that be 4/3? */
2120 else
2121 density = n / (M_PI*r*r);
2122
2123
2124 if (density>maxd) maxd=density;
2125 if (density<mind) mind=density;
2126 sum+=density;
2127 sum2+=density*density;
2128 }
2129 avgd = sum/numToSample;
2130
2131 /* try using some standard deviation stuff instead of this */
2132 saveDensity = avgd;
2133 /*
2134 greaterThan = 0;
2135 for(i=0; i<numToSample; i++)
2136 {
2137 j = rand() % map->numPhotons;
2138
2139 Assign_Vector(Point,(PHOTON_AMF(map->head, j)).Loc);
2140
2141 n=gatherPhotons(Point, 10000000.0, &r, NULL, false, map);
2142
2143 if(mediaMap)
2144 density = 3.0 * n / (4.0*M_PI*r*r*r); // should that be 4/3?
2145 else
2146 density = n / (M_PI*r*r);
2147
2148 if (density>saveDensity)
2149 greaterThan++;
2150 }
2151
2152 density = saveDensity * (DBL)greaterThan / numToSample;
2153 */
2154 density = saveDensity;
2155
2156 if(mediaMap)
2157 {
2158 map->minGatherRad = pow(3.0 * photonOptions.maxGatherCount / (density*M_PI*4.0), 0.3333);
2159 }
2160 else
2161 map->minGatherRad = sqrt(photonOptions.maxGatherCount / (density*M_PI));
2162
2163 lessThan = 0;
2164 for(i=0; i<numToSample; i++)
2165 {
2166 j = rand() % map->numPhotons;
2167
2168 Assign_Vector(Point,(PHOTON_AMF(map->head, j)).Loc);
2169
2170 n=gatherPhotons(Point, map->minGatherRad, &r, NULL, false, map);
2171
2172 if(mediaMap)
2173 density = 3.0 * n / (4.0*M_PI*r*r*r); // should that be 4/3?
2174 else
2175 density = n / (M_PI*r*r);
2176
2177 // count as "lessThan" if the density is below 70% of the average,
2178 // and if it is at least above 5% of the average.
2179 if (density<(saveDensity*0.7) && density>(saveDensity*0.05))
2180 lessThan++;
2181 }
2182
2183 // the 30.0 is a bit of a fudge-factor.
2184 map->minGatherRad*=(1.0+20.0*((DBL)lessThan/(numToSample)));
2185
2186 }
2187
2188 // Now apply the user-defined multiplier, so that the user can tell
2189 // POV to take shortcuts which will improve speed at the expensive of
2190 // quality. Alternatively, the user could tell POV to use a bigger
2191 // radius to improve quality.
2192 map->minGatherRad *= map->minGatherRadMult;
2193
2194 if(mediaMap)
2195 {
2196 /* double the radius if it is a media map */
2197 map->minGatherRad *= 2;
2198 }
2199
2200 /* always do this! - use 6 times the area */
2201 map->gatherRadStep = map->minGatherRad*2;
2202
2203 /* somehow we could automatically determine the number of steps */
2204 }
2205
2206
2207 /**************************************************************
2208
2209 =========== PRIORITY QUEUES ===============
2210 Each priority stores its data in the static variables below (such as
2211 numfound_s) and in the global variables
2212
2213 Preconditions:
2214
2215 static DBL size_sq_s; - maximum radius given squared
2216 static DBL Size_s; - radius
2217 static DBL dmax_s; - square of radius used so far
2218 static int TargetNum_s; - target number
2219 static DBL *pt_s; - center point
2220 static numfound_s; - number of photons in priority queue
2221
2222 these must be allocated:
2223 photonOptions.photonGatherList - array of photons in priority queue
2224 photonOptions.photonDistances - corresponding priorities(distances)
2225
2226 *Each priority queue has the following functions:
2227
2228 function PQInsert(PHOTON *photon, DBL d)
2229
2230 Inserts 'photon' into the priority queue with priority (distance
2231 from target point) 'd'.
2232
2233 void PQDelMax()
2234
2235 Removes the photon with the greates distance (highest priority)
2236 from the queue.
2237
2238 ********************************************************************/
2239
2240 /* try different priority queues */
2241
2242 #define ORDERED 0
2243 #define UNORDERED 1
2244 #define HEAP 2
2245
2246 #define PRI_QUE HEAP
2247
2248 /* -------------- ordered list implementation ----------------- */
2249 #if (PRI_QUE == ORDERED)
PQInsert(PHOTON * photon,DBL d)2250 static void PQInsert(PHOTON *photon, DBL d)
2251 {
2252 int i,j;
2253
2254 Increase_Counter(stats[Priority_Queue_Add]);
2255 /* save this in order, remove maximum, save new dmax_s */
2256
2257 /* store in array and shift - assumption is that we don't have
2258 to shift often */
2259 for (i=0; photonOptions.photonDistances[i]<d && i<(numfound_s); i++);
2260 for (j=numfound_s; j>i; j--)
2261 {
2262 photonOptions.photonGatherList[j] = photonOptions.photonGatherList[j-1];
2263 photonOptions.photonDistances[j] = photonOptions.photonDistances[j-1];
2264 }
2265
2266 numfound_s++;
2267 photonOptions.photonGatherList[i] = photon;
2268 photonOptions.photonDistances[i] = d;
2269 if (numfound_s==TargetNum_s)
2270 dmax_s=photonOptions.photonDistances[numfound_s-1];
2271
2272 }
2273
2274 #ifndef REMOVE_NOT_USED_CODE_WARNINGS_PATCH
PQDelMax()2275 static void PQDelMax()
2276 {
2277 Increase_Counter(stats[Priority_Queue_Remove]);
2278 numfound_s--;
2279 }
2280 #endif
2281 #endif
2282
2283 /* -------------- unordered list implementation ----------------- */
2284 #if (PRI_QUE == UNORDERED)
PQInsert(PHOTON * photon,DBL d)2285 static void PQInsert(PHOTON *photon, DBL d)
2286 {
2287 Increase_Counter(stats[Priority_Queue_Add]);
2288
2289 photonOptions.photonGatherList[numfound_s] = photon;
2290 photonOptions.photonDistances[numfound_s] = d;
2291
2292 if (d>dmax_s)
2293 dmax_s=d;
2294
2295 numfound_s++;
2296 }
2297
2298 #ifndef REMOVE_NOT_USED_CODE_WARNINGS_PATCH
PQDelMax()2299 static void PQDelMax()
2300 {
2301 int i,max;
2302
2303 Increase_Counter(stats[Priority_Queue_Remove]);
2304
2305 max=0;
2306 /* find max */
2307 for(i=1; i<numfound_s; i++)
2308 if (photonOptions.photonDistances[i]>photonOptions.photonDistances[max]) max = i;
2309
2310 /* remove it, shifting the photons */
2311 for(i=max+1; i<numfound_s; i++)
2312 {
2313 photonOptions.photonGatherList[i-1] = photonOptions.photonGatherList[i];
2314 photonOptions.photonDistances[i-1] = photonOptions.photonDistances[i];
2315 }
2316
2317 numfound_s--;
2318
2319 /* find a new dmax_s */
2320 dmax_s=photonOptions.photonDistances[0];
2321 for(i=1; i<numfound_s; i++)
2322 if (photonOptions.photonDistances[i]>dmax_s) dmax_s = photonOptions.photonDistances[i];
2323 }
2324 #endif
2325 #endif
2326
2327 /* -------------- heap implementation ----------------- */
2328 /* this is from Sejwick (spelling?) */
2329 #if (PRI_QUE == HEAP)
2330
PQInsert(PHOTON * photon,DBL d)2331 static void PQInsert(PHOTON *photon, DBL d)
2332 {
2333 Increase_Counter(stats[Priority_Queue_Add]);
2334 DBL* Distances = photonOptions.photonDistances;
2335 PHOTON** Photons = photonOptions.photonGatherList;
2336
2337 unsigned int k = ++numfound_s;
2338
2339 while (k > 1)
2340 {
2341 unsigned int half_k = k / 2;
2342 DBL d_half_k_m1 = Distances[half_k - 1];
2343
2344 if (d < d_half_k_m1)
2345 break;
2346
2347 Distances [k - 1] = d_half_k_m1;
2348 Photons[k - 1] = Photons[half_k - 1];
2349
2350 k = half_k;
2351 }
2352
2353 Distances [k - 1] = d;
2354 Photons[k - 1] = photon;
2355 }
2356
FullPQInsert(PHOTON * photon,DBL d)2357 static void FullPQInsert(PHOTON *photon, DBL d)
2358 {
2359 Increase_Counter(stats[Priority_Queue_Remove]);
2360 DBL* Distances = photonOptions.photonDistances;
2361 PHOTON** Photons = photonOptions.photonGatherList;
2362
2363 int k = 1, k2 = 2;
2364 for (; k2 < numfound_s; k = k2, k2 += k)
2365 {
2366 DBL d_k2_m1 = Distances[k2 - 1],
2367 d_k2 = Distances[k2];
2368
2369 if (d_k2 > d_k2_m1)
2370 {
2371 d_k2_m1 = d_k2;
2372 ++k2;
2373 }
2374
2375 if (!(d_k2_m1 > d))
2376 break;
2377
2378 Distances [k - 1] = d_k2_m1;
2379 Photons[k - 1] = Photons[k2 - 1];
2380 }
2381
2382 if (k2 == numfound_s) {
2383 DBL d_k2_m1 = Distances[k2 - 1];
2384 if (d_k2_m1 > d) {
2385 Distances [k - 1] = d_k2_m1;
2386 Photons[k - 1] = Photons[k2 - 1];
2387 k = k2;
2388 }
2389 }
2390
2391 Distances [k - 1] = d;
2392 Photons[k - 1] = photon;
2393
2394 /* find a new dmax_s */
2395 dmax_s = Distances[0];
2396 }
2397
2398 #endif
2399
2400 /*****************************************************************************
2401
2402 FUNCTION
2403
2404 gatherPhotonsRec()
2405
2406 Recursive part of gatherPhotons
2407 Searches the kd-tree with range start..end (midpoint is pivot)
2408
2409 Preconditions:
2410 same preconditions as priority queue functions
2411 static variable map_s points to the map to use
2412 'start' is the first photon in this range
2413 'end' is the last photon in this range
2414
2415 the range 'start..end' must have been used in building photon map!!!
2416
2417 Postconditions:
2418 photons within the range of start..end are added to the priority
2419 queue (photons may be delted from the queue to make room for photons
2420 of lower priority)
2421
2422 ******************************************************************************/
2423
gatherPhotonsRec(int start,int end)2424 static void gatherPhotonsRec(int start, int end)
2425 {
2426 DBL delta;
2427 int DimToUse;
2428 DBL d,dx,dy,dz;
2429 int mid;
2430 PHOTON *photon;
2431 VECTOR ptToPhoton;
2432 DBL discFix; /* use disc(ellipsoid) for gathering instead of sphere */
2433
2434 /* find midpoint */
2435 mid = (end+start)>>1;
2436 photon = &(PHOTON_AMF(map_s, mid));
2437
2438 DimToUse = photon->info & PH_MASK_XYZ;
2439
2440 /*
2441 * check this photon
2442 */
2443
2444 /* find distance from pt */
2445 ptToPhoton[X] = - pt_s[X] + photon->Loc[X];
2446 ptToPhoton[Y] = - pt_s[Y] + photon->Loc[Y];
2447 ptToPhoton[Z] = - pt_s[Z] + photon->Loc[Z];
2448 /* all distances are squared */
2449 dx = ptToPhoton[X]*ptToPhoton[X];
2450 dy = ptToPhoton[Y]*ptToPhoton[Y];
2451 dz = ptToPhoton[Z]*ptToPhoton[Z];
2452
2453 if (!( ((dx>dmax_s) && ((DimToUse)==X)) ||
2454 ((dy>dmax_s) && ((DimToUse)==Y)) ||
2455 ((dz>dmax_s) && ((DimToUse)==Z)) ))
2456 {
2457 /* it fits manhatten distance - maybe we can use this photon */
2458
2459 /* find euclidian distance (squared) */
2460 d = dx + dy + dz;
2461
2462 /* now fix this distance so that we gather using an ellipsoid
2463 alligned with the surface normal instead of a sphere. This
2464 minimizes false bleeding of photons at sharp corners
2465
2466 dmax_s is square of radius of major axis
2467 dmax_s/16 is " " " " minor " (1/6 of major axis)
2468 */
2469 /*
2470 VDot(discFix,norm_s,ptToPhoton);
2471 discFix*=discFix*(dmax_s/1000.0-dmax_s);
2472 */
2473
2474 if (flattenFactor!=0.0)
2475 {
2476 VDot(discFix,norm_s,ptToPhoton);
2477 discFix = fabs(discFix);
2478 d += flattenFactor*(discFix)*d*16;
2479 }
2480 /* this will add zero if on the plane, and will double distance from
2481 point to photon if it is ptToPhoton is perpendicular to the surface */
2482
2483 if(d < dmax_s)
2484 {
2485 if (numfound_s+1>TargetNum_s)
2486 {
2487 FullPQInsert(photon, d);
2488 sqrt_dmax_s = sqrt(dmax_s);
2489 }
2490 else
2491 PQInsert(photon, d);
2492 }
2493 }
2494
2495 /* now go left & right if appropriate - if going left or right goes out
2496 the current range, then don't go that way. */
2497 /*
2498 delta=pt_s[DimToUse]-photon->Loc[DimToUse];
2499 if(delta<0)
2500 {
2501 if (end>=mid+1) gatherPhotonsRec(start, mid - 1);
2502 if (delta*delta < dmax_s )
2503 if (mid-1>=start) gatherPhotonsRec(mid + 1, end);
2504 }
2505 else
2506 {
2507 if (mid-1>=start) gatherPhotonsRec(mid+1,end);
2508 if (delta*delta < dmax_s )
2509 if (end>=mid+1) gatherPhotonsRec(start, mid - 1);
2510 }
2511 */
2512 delta=pt_s[DimToUse]-photon->Loc[DimToUse];
2513 if(delta<0)
2514 {
2515 /* on left - go left first */
2516 if (pt_s[DimToUse]-sqrt_dmax_s < photon->Loc[DimToUse])
2517 {
2518 if (mid-1>=start)
2519 gatherPhotonsRec(start, mid - 1);
2520 }
2521 if (pt_s[DimToUse]+sqrt_dmax_s > photon->Loc[DimToUse])
2522 {
2523 if(end>=mid+1)
2524 gatherPhotonsRec(mid + 1, end);
2525 }
2526 }
2527 else
2528 {
2529 /* on right - go right first */
2530 if (pt_s[DimToUse]+sqrt_dmax_s > photon->Loc[DimToUse])
2531 {
2532 if(end>=mid+1)
2533 gatherPhotonsRec(mid + 1, end);
2534 }
2535 if (pt_s[DimToUse]-sqrt_dmax_s < photon->Loc[DimToUse])
2536 {
2537 if (mid-1>=start)
2538 gatherPhotonsRec(start, mid - 1);
2539 }
2540 }
2541 }
2542
2543 /*****************************************************************************
2544
2545 FUNCTION
2546
2547 gatherPhotons()
2548
2549 gathers photons from the global photon map
2550
2551 Preconditons:
2552
2553 photonOptions.photonGatherList and photonOptions.photonDistances
2554 are allocated and are each maxGatherCount in length
2555
2556 'Size' - maximum search radius
2557 'r' points to a double
2558
2559 BuildPhotonMaps() has been called for this scene.
2560
2561 Postconditions:
2562
2563 *r is radius actually used for gathereing (maximum value is 'Size')
2564 photonOptions.photonGatherList and photonOptions.photonDistances
2565 contain the gathered photons
2566 returns number of photons gathered
2567
2568 ******************************************************************************/
2569
gatherPhotons(VECTOR pt,DBL Size,DBL * r,VECTOR norm,int flatten,PHOTON_MAP * map)2570 int gatherPhotons(VECTOR pt, DBL Size, DBL *r, VECTOR norm, int flatten, PHOTON_MAP *map)
2571 {
2572 if (map->numPhotons<=0) return 0; /* no crashes, please... */
2573
2574 /* set the static variables */
2575 numfound_s=0;
2576 size_sq_s = Size*Size;
2577 dmax_s = size_sq_s;
2578 sqrt_dmax_s = Size;
2579 norm_s = norm;
2580
2581 if(flatten)
2582 {
2583 flattenFactor = 1.0;
2584 }
2585 else
2586 {
2587 flattenFactor = 0.0;
2588 }
2589
2590 Size_s = Size;
2591 TargetNum_s = photonOptions.maxGatherCount;
2592 pt_s = pt;
2593
2594 map_s = map->head;
2595
2596 /* now search the kd-tree recursively */
2597 gatherPhotonsRec(0, map->numPhotons-1);
2598
2599 /* set the radius variable */
2600 *r = sqrt_dmax_s;
2601
2602 /* return the number of photons found */
2603 return(numfound_s);
2604 }
2605
2606
2607 /******************************************************************
2608 stuff grabbed from radiosit.h & radiosit.c
2609 ******************************************************************/
2610
2611 extern const BYTE_XYZ rad_samples[];
2612
VUnpack(VECTOR dest_vec,const BYTE_XYZ * pack_vec)2613 static void VUnpack(VECTOR dest_vec, const BYTE_XYZ * pack_vec)
2614 {
2615 dest_vec[X] = ((double)pack_vec->x * (1./ 255.))*2.-1.;
2616 dest_vec[Y] = ((double)pack_vec->y * (1./ 255.))*2.-1.;
2617 dest_vec[Z] = ((double)pack_vec->z * (1./ 255.));
2618
2619 VNormalizeEq(dest_vec); /* already good to about 1%, but we can do better */
2620 }
2621
2622 /******************************************************************
2623 ******************************************************************/
ChooseRay(RAY * NewRay,VECTOR Normal,RAY *,VECTOR Raw_Normal,int WhichRay)2624 void ChooseRay(RAY *NewRay, VECTOR Normal, RAY * /*Ray*/, VECTOR Raw_Normal, int WhichRay)
2625 {
2626 VECTOR random_vec, up, n2, n3;
2627 int i;
2628 DBL /*n,*/ NRay_Direction;
2629
2630 #define REFLECT_FOR_RADIANCE 0
2631 #if (REFLECT_FOR_RADIANCE)
2632 /* Get direction of reflected ray. */
2633 n = -2.0 * (Ray->Direction[X] * Normal[X] + Ray->Direction[Y] * Normal[Y] + Ray->Direction[Z] * Normal[Z]);
2634
2635 VLinComb2(NewRay->Direction, n, Normal, 1.0, Ray->Direction);
2636
2637 VDot(NRay_Direction, NewRay->Direction, Raw_Normal);
2638 if (NRay_Direction < 0.0)
2639 {
2640 /* subtract 2*(projection of NRay.Direction onto Raw_Normal)
2641 from NRay.Direction */
2642 DBL Proj;
2643 Proj = NRay_Direction * -2;
2644 VAddScaledEq(NewRay->Direction, Proj, Raw_Normal);
2645 }
2646 return;
2647 #else
2648 Assign_Vector(NewRay->Direction, Normal);
2649 #endif
2650
2651 if ( fabs(fabs(NewRay->Direction[Z])- 1.) < .1 ) {
2652 /* too close to vertical for comfort, so use cross product with horizon */
2653 up[X] = 0.; up[Y] = 1.; up[Z] = 0.;
2654 }
2655 else
2656 {
2657 up[X] = 0.; up[Y] = 0.; up[Z] = 1.;
2658 }
2659
2660 VCross(n2, NewRay->Direction, up); VNormalizeEq(n2);
2661 VCross(n3, NewRay->Direction, n2); VNormalizeEq(n3);
2662
2663 /*i = (int)(FRAND()*1600);*/
2664 i = WhichRay;
2665 WhichRay = (WhichRay + 1) % 1600;
2666
2667 VUnpack(random_vec, &rad_samples[i]);
2668
2669 if ( fabs(NewRay->Direction[Z] - 1.) < .001 ) /* pretty well straight Z, folks */
2670 {
2671 /* we are within 1/20 degree of pointing in the Z axis. */
2672 /* use all vectors as is--they're precomputed this way */
2673 Assign_Vector(NewRay->Direction, random_vec);
2674 }
2675 else
2676 {
2677 NewRay->Direction[X] = n2[X]*random_vec[X] + n3[X]*random_vec[Y] + NewRay->Direction[X]*random_vec[Z];
2678 NewRay->Direction[Y] = n2[Y]*random_vec[X] + n3[Y]*random_vec[Y] + NewRay->Direction[Y]*random_vec[Z];
2679 NewRay->Direction[Z] = n2[Z]*random_vec[X] + n3[Z]*random_vec[Y] + NewRay->Direction[Z]*random_vec[Z];
2680 }
2681
2682 /* if our new ray goes through, flip it back across raw_normal */
2683
2684 VDot(NRay_Direction, NewRay->Direction, Raw_Normal);
2685 if (NRay_Direction < 0.0)
2686 {
2687 /* subtract 2*(projection of NRay.Direction onto Raw_Normal)
2688 from NRay.Direction */
2689 DBL Proj;
2690 Proj = NRay_Direction * -2;
2691 VAddScaledEq(NewRay->Direction, Proj, Raw_Normal);
2692 }
2693
2694 VNormalizeEq(NewRay->Direction);
2695 }
2696
GetPhotonStat(POVMSType a)2697 int GetPhotonStat(POVMSType a)
2698 {
2699 switch(a)
2700 {
2701 case kPOVAttrib_ObjectPhotonCount:
2702 return gPhotonStat_i;
2703 case kPOVAttrib_TotalPhotonCount:
2704 return photonOptions.photonMap.numPhotons;
2705 case kPOVAttrib_MediaPhotonCount:
2706 return photonOptions.mediaPhotonMap.numPhotons;
2707 case kPOVAttrib_PhotonXSamples:
2708 return gPhotonStat_x_samples;
2709 case kPOVAttrib_PhotonYSamples:
2710 return gPhotonStat_y_samples;
2711 case kPOVAttrib_CurrentPhotonCount:
2712 return gPhotonStat_end;
2713 }
2714
2715 return 0;
2716 }
2717
2718 END_POV_NAMESPACE
2719