ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmap.c
(Generate patch)

Comparing ray/src/rt/pmap.c (file contents):
Revision 2.4 by rschregle, Thu Apr 23 20:02:04 2015 UTC vs.
Revision 2.16 by rschregle, Wed Nov 21 19:33:30 2018 UTC

# Line 1 | Line 1
1 + #ifndef lint
2 + static const char RCSid[] = "$Id$";
3 + #endif
4 +
5 +
6   /*
7 <   ==================================================================
7 >   ======================================================================
8     Photon map main module
9  
10     Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
11     (c) Fraunhofer Institute for Solar Energy Systems,
12     (c) Lucerne University of Applied Sciences and Arts,
13 <   supported by the Swiss National Science Foundation (SNSF, #147053)
14 <   ==================================================================
13 >       supported by the Swiss National Science Foundation (SNSF, #147053)
14 >   ======================================================================
15    
16     $Id$
17   */
18  
19  
15 –
20   #include "pmap.h"
21   #include "pmapmat.h"
22   #include "pmapsrc.h"
# Line 22 | Line 26
26   #include "pmapdiag.h"
27   #include "otypes.h"
28   #include <time.h>
29 < #include <sys/stat.h>
29 > #if NIX
30 >   #include <sys/stat.h>
31 >   #include <sys/mman.h>
32 >   #include <sys/wait.h>
33 > #endif
34  
35  
28 –
29 – extern char *octname;
30 –
31 – static char PmapRevision [] = "$Revision$";
32 –
33 –
34 –
35 – /* Photon map lookup functions per type */
36 – void (*pmapLookup [NUM_PMAP_TYPES])(PhotonMap*, RAY*, COLOR) = {
37 –   photonDensity, photonPreCompDensity, photonDensity, volumePhotonDensity,
38 –   photonDensity, NULL
39 – };
40 –
41 –
42 –
43 – void colorNorm (COLOR c)
44 – /* Normalise colour channels to average of 1 */
45 – {
46 –   const float avg = colorAvg(c);
47 –  
48 –   if (!avg)
49 –      return;
50 –      
51 –   c [0] /= avg;
52 –   c [1] /= avg;
53 –   c [2] /= avg;
54 – }
55 –
56 –
57 –
58 – void loadPmaps (PhotonMap **pmaps, const PhotonMapParams *parm)
59 – {
60 –   unsigned t;
61 –   struct stat octstat, pmstat;
62 –   PhotonMap *pm;
63 –   PhotonMapType type;
64 –  
65 –   for (t = 0; t < NUM_PMAP_TYPES; t++)
66 –      if (setPmapParam(&pm, parm + t)) {        
67 –         /* Check if photon map newer than octree */
68 –         if (pm -> fileName && octname &&
69 –             !stat(pm -> fileName, &pmstat) && !stat(octname, &octstat) &&
70 –             octstat.st_mtime > pmstat.st_mtime) {
71 –            sprintf(errmsg, "photon map in file %s may be stale",
72 –                    pm -> fileName);
73 –            error(USER, errmsg);
74 –         }
75 –        
76 –         /* Load photon map from file and get its type */
77 –         if ((type = loadPhotonMap(pm, pm -> fileName)) == PMAP_TYPE_NONE)
78 –            error(USER, "failed loading photon map");
79 –            
80 –         /* Assign to appropriate photon map type (deleting previously
81 –          * loaded photon map of same type if necessary) */
82 –         if (pmaps [type]) {
83 –            deletePhotons(pmaps [type]);
84 –            free(pmaps [type]);
85 –         }
86 –         pmaps [type] = pm;
87 –        
88 –         /* Check for invalid density estimate bandwidth */                            
89 –         if (pm -> maxGather > pm -> heapSize) {
90 –            error(WARNING, "adjusting density estimate bandwidth");
91 –            pm -> minGather = pm -> maxGather = pm -> heapSize;
92 –         }
93 –      }
94 – }
95 –
96 –
97 –
36   void savePmaps (const PhotonMap **pmaps, int argc, char **argv)
37   {
38     unsigned t;
39    
40     for (t = 0; t < NUM_PMAP_TYPES; t++) {
41        if (pmaps [t])
42 <         savePhotonMap(pmaps [t], pmaps [t] -> fileName, t, argc, argv);
42 >         savePhotonMap(pmaps [t], pmaps [t] -> fileName, argc, argv);
43     }
44   }                  
45  
46  
109 –  
110 – void cleanUpPmaps (PhotonMap **pmaps)
111 – {
112 –   unsigned t;
113 –  
114 –   for (t = 0; t < NUM_PMAP_TYPES; t++) {
115 –      if (pmaps [t]) {
116 –         deletePhotons(pmaps [t]);
117 –         free(pmaps [t]);
118 –      }
119 –   }
120 – }
121 –
122 –
47      
48   static int photonParticipate (RAY *ray)
49   /* Trace photon through participating medium. Returns 1 if passed through,
# Line 128 | Line 52 | static int photonParticipate (RAY *ray)
52     int i;
53     RREAL cosTheta, cosPhi, du, dv;
54     const float cext = colorAvg(ray -> cext),
55 <               albedo = colorAvg(ray -> albedo);
55 >               albedo = colorAvg(ray -> albedo),
56 >               gecc2 = ray -> gecc * ray -> gecc;
57     FVECT u, v;
58     COLOR cvext;
59  
# Line 136 | Line 61 | static int photonParticipate (RAY *ray)
61     ray -> rmax = -log(pmapRandom(mediumState)) / cext;
62    
63     while (!localhit(ray, &thescene)) {
64 +      if (!incube(&thescene, ray -> rop)) {
65 +         /* Terminate photon if it has leaked from the scene */
66 + #ifdef DEBUG_PMAP
67 +         fprintf(stderr,
68 +                 "Volume photon leaked from scene at [%.3f %.3f %.3f]\n",
69 +                 ray -> rop [0], ray -> rop [1], ray -> rop [2]);
70 + #endif                
71 +         return 0;
72 +      }
73 +        
74        setcolor(cvext, exp(-ray -> rmax * ray -> cext [0]),
75                        exp(-ray -> rmax * ray -> cext [1]),
76                        exp(-ray -> rmax * ray -> cext [2]));
# Line 145 | Line 80 | static int photonParticipate (RAY *ray)
80        colorNorm(ray -> rcol);
81        VCOPY(ray -> rorg, ray -> rop);
82        
83 + #if 0
84 +      if (albedo > FTINY && ray -> rlvl > 0)
85 + #else
86 +      /* Store volume photons unconditionally in mist to also account for
87 +         direct inscattering from sources */
88        if (albedo > FTINY)
89 + #endif      
90           /* Add to volume photon map */
91 <         if (ray -> rlvl > 0) addPhoton(volumePmap, ray);
91 >         newPhoton(volumePmap, ray);
92          
93        /* Absorbed? */
94 <      if (pmapRandom(rouletteState) > albedo) return 0;
94 >      if (pmapRandom(rouletteState) > albedo)
95 >         return 0;
96        
97        /* Colour bleeding without attenuation (?) */
98        multcolor(ray -> rcol, ray -> albedo);
99        scalecolor(ray -> rcol, 1 / albedo);    
100        
101        /* Scatter photon */
102 <      cosTheta = ray -> gecc <= FTINY ? 2 * pmapRandom(scatterState) - 1
103 <                                      : 1 / (2 * ray -> gecc) *
104 <                                            (1 + ray -> gecc * ray -> gecc -
105 <                                               (1 - ray -> gecc * ray -> gecc) /
106 <                                               (1 - ray -> gecc + 2 * ray -> gecc *
165 <                                                  pmapRandom(scatterState)));
102 >      cosTheta = ray -> gecc <= FTINY
103 >                    ? 2 * pmapRandom(scatterState) - 1
104 >                    : 0.5 * (1 + gecc2 -
105 >                         (1 - gecc2) / (1 - ray -> gecc + 2 * ray -> gecc *
106 >                         pmapRandom(scatterState))) / ray -> gecc;
107                                                    
108        cosPhi = cos(2 * PI * pmapRandom(scatterState));
109        du = dv = sqrt(1 - cosTheta * cosTheta);   /* sin(theta) */
# Line 184 | Line 125 | static int photonParticipate (RAY *ray)
125        ray -> rlvl++;
126        ray -> rmax = -log(pmapRandom(mediumState)) / cext;
127     }  
128 <    
128 >  
129 >   /* Passed through medium until intersecting local object */  
130     setcolor(cvext, exp(-ray -> rot * ray -> cext [0]),
131                     exp(-ray -> rot * ray -> cext [1]),
132                     exp(-ray -> rot * ray -> cext [2]));
133                    
134     /* Modify ray color and normalise */
135     multcolor(ray -> rcol, cvext);
136 <   colorNorm(ray -> rcol);
137 <  
196 <   /* Passed through medium */  
136 >   colorNorm(ray -> rcol);  
137 >
138     return 1;
139   }
140  
# Line 203 | Line 144 | void tracePhoton (RAY *ray)
144   /* Follow photon as it bounces around... */
145   {
146     long mod;
147 <   OBJREC* mat;
147 >   OBJREC *mat, *port = NULL;
148 >  
149 >   if (!ray -> parent) {
150 >      /* !!!  PHOTON PORT REJECTION SAMPLING HACK: get photon port for
151 >       * !!!  primary ray from ray -> ro, then reset the latter to NULL so
152 >       * !!!  as not to interfere with localhit() */
153 >      port = ray -> ro;
154 >      ray -> ro = NULL;
155 >   }
156  
157     if (ray -> rlvl > photonMaxBounce) {
158 + #ifdef PMAP_RUNAWAY_WARN  
159        error(WARNING, "runaway photon!");
160 + #endif      
161        return;
162     }
163 <  
163 >  
164     if (colorAvg(ray -> cext) > FTINY && !photonParticipate(ray))
165        return;
166 <      
166 >
167     if (localhit(ray, &thescene)) {
168        mod = ray -> ro -> omod;
169 <      
169 >
170 >      if (port && ray -> ro != port) {
171 >         /* !!! PHOTON PORT REJECTION SAMPLING HACK !!!
172 >          * Terminate photon if emitted from port without intersecting it;
173 >          * this can happen when the port's partitions extend beyond its
174 >          * actual geometry, e.g.  with polygons.  Since the total flux
175 >          * relayed by the port is based on the (in this case) larger
176 >          * partition area, it is overestimated; terminating these photons
177 >          * constitutes rejection sampling and thereby compensates any bias
178 >          * incurred by the overestimated flux.  */
179 > #ifdef PMAP_PORTREJECT_WARN
180 >         sprintf(errmsg, "photon outside port %s", ray -> ro -> oname);
181 >         error(WARNING, errmsg);
182 > #endif        
183 >         return;
184 >      }
185 >
186        if ((ray -> clipset && inset(ray -> clipset, mod)) || mod == OVOID) {
187           /* Transfer ray if modifier is VOID or clipped within antimatta */
188           RAY tray;
# Line 233 | Line 200 | void tracePhoton (RAY *ray)
200  
201  
202   static void preComputeGlobal (PhotonMap *pmap)
203 < /* Precompute irradiance from global photons for final gathering using
204 <   the first finalGather * pmap -> heapSize photons in the heap. Returns
205 <   new heap with precomputed photons. */
203 > /* Precompute irradiance from global photons for final gathering for  
204 >   a random subset of finalGather * pmap -> numPhotons photons, and builds
205 >   the photon map, discarding the original photons. */
206 > /* !!! NOTE: PRECOMPUTATION WITH OOC CURRENTLY WITHOUT CACHE !!! */  
207   {
208 <   unsigned long i, nuHeapSize;
209 <   unsigned j;
210 <   Photon *nuHeap, *p;
211 <   COLOR irrad;
212 <   RAY ray;
213 <   float nuMinPos [3], nuMaxPos [3];
208 >   unsigned long  i, numPreComp;
209 >   unsigned       j;
210 >   PhotonIdx      pIdx;
211 >   Photon         photon;
212 >   RAY            ray;
213 >   PhotonMap      nuPmap;
214  
215 <   repComplete = nuHeapSize = finalGather * pmap -> heapSize;
215 >   repComplete = numPreComp = finalGather * pmap -> numPhotons;
216    
217 <   if (photonRepTime) {
217 >   if (verbose) {
218        sprintf(errmsg,
219 <              "Precomputing irradiance for %ld global photons...\n",
220 <              nuHeapSize);
219 >              "\nPrecomputing irradiance for %ld global photons\n",
220 >              numPreComp);
221        eputs(errmsg);
222 + #if NIX      
223        fflush(stderr);
224 + #endif      
225     }
226    
227 <   p = nuHeap = (Photon*)malloc(nuHeapSize * sizeof(Photon));
228 <   if (!nuHeap)
229 <      error(USER, "can't allocate photon heap");
230 <      
231 <   for (j = 0; j <= 2; j++) {
232 <      nuMinPos [j] = FHUGE;
263 <      nuMaxPos [j] = -FHUGE;
264 <   }
227 >   /* Copy photon map for precomputed photons */
228 >   memcpy(&nuPmap, pmap, sizeof(PhotonMap));
229 >
230 >   /* Zero counters, init new heap and extents */  
231 >   nuPmap.numPhotons = 0;  
232 >   initPhotonHeap(&nuPmap);
233    
234 +   for (j = 0; j < 3; j++) {  
235 +      nuPmap.minPos [j] = FHUGE;
236 +      nuPmap.maxPos [j] = -FHUGE;
237 +   }
238 +
239     /* Record start time, baby */
240     repStartTime = time(NULL);
241 <   #ifdef SIGCONT
242 <      signal(SIGCONT, pmapPreCompReport);
243 <   #endif
241 > #ifdef SIGCONT
242 >   signal(SIGCONT, pmapPreCompReport);
243 > #endif
244     repProgress = 0;
272 –   memcpy(nuHeap, pmap -> heap, nuHeapSize * sizeof(Photon));
245    
246 <   for (i = 0, p = nuHeap; i < nuHeapSize; i++, p++) {
247 <      ray.ro = NULL;
248 <      VCOPY(ray.rop, p -> pos);
246 >   photonRay(NULL, &ray, PRIMARY, NULL);
247 >   ray.ro = NULL;
248 >  
249 >   for (i = 0; i < numPreComp; i++) {
250 >      /* Get random photon from stratified distribution in source heap to
251 >       * avoid duplicates and clustering */
252 >      pIdx = firstPhoton(pmap) +
253 >             (unsigned long)((i + pmapRandom(pmap -> randState)) /
254 >                             finalGather);
255 >      getPhoton(pmap, pIdx, &photon);
256        
257 <      /* Update min and max positions & set ray normal */
258 <      for (j = 0; j < 3; j++) {
259 <         if (p -> pos [j] < nuMinPos [j]) nuMinPos [j] = p -> pos [j];
260 <         if (p -> pos [j] > nuMaxPos [j]) nuMaxPos [j] = p -> pos [j];
282 <         ray.ron [j] = p -> norm [j] / 127.0;
283 <      }
257 >      /* Init dummy photon ray with intersection at photon position */
258 >      VCOPY(ray.rop, photon.pos);
259 >      for (j = 0; j < 3; j++)
260 >         ray.ron [j] = photon.norm [j] / 127.0;
261        
262 <      photonDensity(pmap, &ray, irrad);
263 <      setPhotonFlux(p, irrad);
262 >      /* Get density estimate at photon position */
263 >      photonDensity(pmap, &ray, ray.rcol);
264 >                  
265 >      /* Append photon to new heap from ray */
266 >      newPhoton(&nuPmap, &ray);
267 >      
268 >      /* Update progress */
269        repProgress++;
270        
271        if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
272           pmapPreCompReport();
273 <      #ifdef SIGCONT
274 <         else signal(SIGCONT, pmapPreCompReport);
275 <      #endif
273 > #ifdef SIGCONT
274 >      else signal(SIGCONT, pmapPreCompReport);
275 > #endif
276     }
277    
278 <   #ifdef SIGCONT  
279 <      signal(SIGCONT, SIG_DFL);
298 <   #endif
278 >   /* Flush heap */
279 >   flushPhotonHeap(&nuPmap);
280    
281 <   /* Replace & rebuild heap */
282 <   free(pmap -> heap);
283 <   pmap -> heap = nuHeap;
303 <   pmap -> heapSize = pmap -> heapEnd = nuHeapSize;
304 <   VCOPY(pmap -> minPos, nuMinPos);
305 <   VCOPY(pmap -> maxPos, nuMaxPos);
281 > #ifdef SIGCONT  
282 >   signal(SIGCONT, SIG_DFL);
283 > #endif
284    
285 <   if (photonRepTime) {
286 <      eputs("Rebuilding global photon heap...\n");
285 >   /* Trash original pmap, replace with precomputed one */
286 >   deletePhotons(pmap);
287 >   memcpy(pmap, &nuPmap, sizeof(PhotonMap));
288 >  
289 >   if (verbose) {
290 >      eputs("\nRebuilding precomputed photon map\n");
291 > #if NIX      
292        fflush(stderr);
293 + #endif      
294     }
295 <  
296 <   balancePhotons(pmap, NULL);
295 >
296 >   /* Rebuild underlying data structure, destroying heap */  
297 >   buildPhotonMap(pmap, NULL, NULL, 1);
298   }
299  
300  
301  
302 < void distribPhotons (PhotonMap **pmaps)
302 > typedef struct {
303 >   unsigned long  numPhotons [NUM_PMAP_TYPES],
304 >                  numEmitted, numComplete;
305 > } PhotonCnt;
306 >
307 >
308 >
309 > void distribPhotons (PhotonMap **pmaps, unsigned numProc)
310   {
311 <   EmissionMap emap;
312 <   char errmsg2 [128];
313 <   unsigned t, srcIdx, passCnt = 0, prePassCnt = 0;
314 <   double totalFlux = 0;
315 <   PhotonMap *pm;
311 >   EmissionMap    emap;
312 >   char           errmsg2 [128], shmFname [PMAP_TMPFNLEN];
313 >   unsigned       t, srcIdx, proc;
314 >   double         totalFlux = 0;
315 >   int            shmFile, stat, pid;
316 >   PhotonMap      *pm;
317 >   PhotonCnt      *photonCnt;
318    
319 <   for (t = 0; t < NUM_PMAP_TYPES && !photonMaps [t]; t++);
319 >   for (t = 0; t < NUM_PMAP_TYPES && !pmaps [t]; t++);
320 >  
321     if (t >= NUM_PMAP_TYPES)
322 <      error(USER, "no photon maps defined");
322 >      error(USER, "no photon maps defined in distribPhotons");
323        
324     if (!nsources)
325 <      error(USER, "no light sources");
325 >      error(USER, "no light sources in distribPhotons");
326  
327     /* ===================================================================
328      * INITIALISATION - Set up emission and scattering funcs
# Line 336 | Line 331 | void distribPhotons (PhotonMap **pmaps)
331     emap.maxPartitions = MAXSPART;
332     emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1);
333     if (!emap.partitions)
334 <      error(INTERNAL, "can't allocate source partitions");
334 >      error(INTERNAL, "can't allocate source partitions in distribPhotons");
335        
336     /* Initialise all defined photon maps */
337     for (t = 0; t < NUM_PMAP_TYPES; t++)
338 <      initPhotonMap(photonMaps [t], t);
338 >      if (pmaps [t]) {
339 >         initPhotonMap(pmaps [t], t);
340 >         /* Open photon heapfile */
341 >         initPhotonHeap(pmaps [t]);
342 >         /* Per-subprocess target count */
343 >         pmaps [t] -> distribTarget /= numProc;
344 >        
345 >         if (!pmaps [t] -> distribTarget)
346 >            error(INTERNAL, "no photons to distribute in distribPhotons");
347 >      }
348  
349     initPhotonEmissionFuncs();
350     initPhotonScatterFuncs();
351    
352 <   /* Get photon ports if specified */
353 <   if (ambincl == 1)
350 <      getPhotonPorts();
352 >   /* Get photon ports from modifier list */
353 >   getPhotonPorts(photonPortList);
354  
355     /* Get photon sensor modifiers */
356     getPhotonSensors(photonSensorList);
357    
358 <   /* Seed RNGs for photon distribution */
359 <   pmapSeed(randSeed, partState);
360 <   pmapSeed(randSeed, emitState);
361 <   pmapSeed(randSeed, cntState);
362 <   pmapSeed(randSeed, mediumState);
363 <   pmapSeed(randSeed, scatterState);
364 <   pmapSeed(randSeed, rouletteState);
358 > #if NIX
359 >   /* Set up shared mem for photon counters (zeroed by ftruncate) */
360 >   strcpy(shmFname, PMAP_TMPFNAME);
361 >   shmFile = mkstemp(shmFname);
362 >
363 >   if (shmFile < 0 || ftruncate(shmFile, sizeof(*photonCnt)) < 0)
364 >      error(SYSTEM, "failed shared mem init in distribPhotons");
365 >
366 >   photonCnt = mmap(NULL, sizeof(*photonCnt), PROT_READ | PROT_WRITE,
367 >                    MAP_SHARED, shmFile, 0);
368 >                    
369 >   if (photonCnt == MAP_FAILED)
370 >      error(SYSTEM, "failed mapping shared memory in distribPhotons");
371 > #else
372 >   /* Allocate photon counters statically on Windoze */
373 >   if (!(photonCnt = malloc(sizeof(PhotonCnt))))
374 >      error(SYSTEM, "failed trivial malloc in distribPhotons");
375 >   photonCnt -> numEmitted = photonCnt -> numComplete = 0;      
376 > #endif /* NIX */
377 >
378 >   if (verbose) {
379 >      sprintf(errmsg, "\nIntegrating flux from %d sources", nsources);
380 >      
381 >      if (photonPorts) {
382 >         sprintf(errmsg2, " via %d ports", numPhotonPorts);
383 >         strcat(errmsg, errmsg2);
384 >      }
385 >      
386 >      strcat(errmsg, "\n");
387 >      eputs(errmsg);
388 >   }  
389    
363 –   if (photonRepTime)
364 –      eputs("\n");
365 –  
390     /* ===================================================================
391      * FLUX INTEGRATION - Get total photon flux from light sources
392      * =================================================================== */
393 <   for (srcIdx = 0; srcIdx < nsources; srcIdx++) {        
393 >   for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
394        unsigned portCnt = 0;
395        emap.src = source + srcIdx;
396        
397 <      do {
397 >      do {  /* Need at least one iteration if no ports! */
398           emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt
399                                                     : NULL;
400           photonPartition [emap.src -> so -> otype] (&emap);
401          
402 <         if (photonRepTime) {
403 <            sprintf(errmsg, "Integrating flux from source %s ",
402 >         if (verbose) {
403 >            sprintf(errmsg, "\tIntegrating flux from source %s ",
404                      source [srcIdx].so -> oname);
405 <                    
405 >
406              if (emap.port) {
407                 sprintf(errmsg2, "via port %s ",
408                         photonPorts [portCnt].so -> oname);
409                 strcat(errmsg, errmsg2);
410              }
411 <            
412 <            sprintf(errmsg2, "(%lu partitions)...\n", emap.numPartitions);
411 >
412 >            sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions);
413              strcat(errmsg, errmsg2);
414              eputs(errmsg);
415 + #if NIX            
416              fflush(stderr);
417 + #endif            
418           }
419          
420           for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
# Line 403 | Line 429 | void distribPhotons (PhotonMap **pmaps)
429  
430     if (totalFlux < FTINY)
431        error(USER, "zero flux from light sources");
406 –
407 –   /* Record start time and enable progress report signal handler */
408 –   repStartTime = time(NULL);
409 –   #ifdef SIGCONT
410 –      signal(SIGCONT, pmapDistribReport);
411 –   #endif
412 –   repProgress = prePassCnt = 0;
413 –  
414 –   if (photonRepTime)
415 –      eputs("\n");
416 –  
417 –   /* ===================================================================
418 –    * 2-PASS PHOTON DISTRIBUTION
419 –    * Pass 1 (pre):  emit fraction of target photon count
420 –    * Pass 2 (main): based on outcome of pass 1, estimate remaining number
421 –    *                of photons to emit to approximate target count
422 –    * =================================================================== */
423 –   do {
424 –      double numEmit;
432        
433 <      if (!passCnt) {  
434 <         /* INIT PASS 1 */
428 <         /* Skip if no photons contributed after sufficient iterations; make
429 <          * it clear to user which photon maps are missing so (s)he can
430 <          * check the scene geometry and materials */
431 <         if (++prePassCnt > maxPreDistrib) {
432 <            sprintf(errmsg, "too many prepasses");
433 >   /* Record start time for progress reports */
434 >   repStartTime = time(NULL);
435  
436 <            for (t = 0; t < NUM_PMAP_TYPES; t++)
437 <               if (photonMaps [t] && !photonMaps [t] -> heapEnd) {
438 <                  sprintf(errmsg2, ", no %s photons stored", pmapName [t]);
439 <                  strcat(errmsg, errmsg2);
438 <               }
439 <            
440 <            error(USER, errmsg);
441 <            break;
442 <         }
436 >   if (verbose) {
437 >      sprintf(errmsg, "\nPhoton distribution @ %d procs\n", numProc);
438 >      eputs(errmsg);
439 >   }
440  
441 <         /* Num to emit is fraction of minimum target count */
442 <         numEmit = FHUGE;
441 >   /* MAIN LOOP */  
442 >   for (proc = 0; proc < numProc; proc++) {
443 > #if NIX          
444 >      if (!(pid = fork())) {
445 >         /* SUBPROCESS ENTERS HERE; open and mmapped files inherited */
446 > #else
447 >      if (1) {
448 >         /* No subprocess under Windoze */
449 > #endif
450 >         /* Local photon counters for this subprocess */
451 >         unsigned       passCnt = 0, prePassCnt = 0;
452 >         unsigned long  lastNumPhotons [NUM_PMAP_TYPES];
453 >         unsigned long  localNumEmitted = 0; /* Num photons emitted by this
454 >                                                subprocess alone */
455          
456 <         for (t = 0; t < NUM_PMAP_TYPES; t++)
457 <            if (photonMaps [t])
458 <               numEmit = min(photonMaps [t] -> distribTarget, numEmit);
456 >         /* Seed RNGs from PID for decorellated photon distribution */
457 >         pmapSeed(randSeed + proc, partState);
458 >         pmapSeed(randSeed + (proc + 1) % numProc, emitState);
459 >         pmapSeed(randSeed + (proc + 2) % numProc, cntState);
460 >         pmapSeed(randSeed + (proc + 3) % numProc, mediumState);
461 >         pmapSeed(randSeed + (proc + 4) % numProc, scatterState);
462 >         pmapSeed(randSeed + (proc + 5) % numProc, rouletteState);
463                
464 <         numEmit *= preDistrib;
465 <      }
464 > #ifdef DEBUG_PMAP          
465 >         /* Output child process PID after random delay to prevent corrupted
466 >          * console output due to race condition */
467 >         usleep(1e6 * pmapRandom(rouletteState));
468 >         fprintf(stderr, "Proc %d: PID = %d "
469 >                 "(waiting 10 sec to attach debugger...)\n",
470 >                 proc, getpid());
471 >         /* Allow time for debugger to attach to child process */
472 >         sleep(10);
473 > #endif            
474  
454 –      else {            
455 –         /* INIT PASS 2 */
456 –         /* Based on the outcome of the predistribution we can now estimate
457 –          * how many more photons we have to emit for each photon map to
458 –          * meet its respective target count. This value is clamped to 0 in
459 –          * case the target has already been exceeded in the pass 1. Note
460 –          * repProgress is the number of photons emitted thus far, while
461 –          * heapEnd is the number of photons stored in each photon map. */
462 –         double maxDistribRatio = 0;
463 –
464 –         /* Set the distribution ratio for each map; this indicates how many
465 –          * photons of each respective type are stored per emitted photon,
466 –          * and is used as probability for storing a photon by addPhoton().
467 –          * Since this biases the photon density, addPhoton() promotes the
468 –          * flux of stored photons to compensate. */
475           for (t = 0; t < NUM_PMAP_TYPES; t++)
476 <            if ((pm = photonMaps [t])) {
477 <               pm -> distribRatio = (double)pm -> distribTarget /
478 <                                    pm -> heapEnd - 1;
479 <
480 <               /* Check if photon map "overflowed", i.e. exceeded its target
481 <                * count in the prepass; correcting the photon flux via the
482 <                * distribution ratio is no longer possible, as no more
483 <                * photons of this type will be stored, so notify the user
484 <                * rather than deliver incorrect results.
479 <                * In future we should handle this more intelligently by
480 <                * using the photonFlux in each photon map to individually
481 <                * correct the flux after distribution. */
482 <               if (pm -> distribRatio <= FTINY) {
483 <                  sprintf(errmsg,
484 <                          "%s photon map overflow in prepass, reduce -apD",
485 <                          pmapName [t]);
486 <                  error(INTERNAL, errmsg);
487 <               }
488 <                  
489 <               maxDistribRatio = max(pm -> distribRatio, maxDistribRatio);
490 <            }
491 <        
492 <         /* Normalise distribution ratios and calculate number of photons to
493 <          * emit in main pass */
494 <         for (t = 0; t < NUM_PMAP_TYPES; t++)
495 <            if ((pm = photonMaps [t]))
496 <               pm -> distribRatio /= maxDistribRatio;
497 <              
498 <         if ((numEmit = repProgress * maxDistribRatio) < FTINY)
499 <            /* No photons left to distribute in main pass */
500 <            break;
501 <      }
502 <      
503 <      /* Set completion count for progress report */
504 <      repComplete = numEmit + repProgress;                            
505 <      
506 <      /* PHOTON DISTRIBUTION LOOP */
507 <      for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
508 <         unsigned portCnt = 0;
509 <         emap.src = source + srcIdx;
510 <                  
476 >            lastNumPhotons [t] = 0;
477 >            
478 >         /* =============================================================
479 >          * 2-PASS PHOTON DISTRIBUTION
480 >          * Pass 1 (pre):  emit fraction of target photon count
481 >          * Pass 2 (main): based on outcome of pass 1, estimate remaining
482 >          *                number of photons to emit to approximate target
483 >          *                count
484 >          * ============================================================= */
485           do {
486 <            emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt
513 <                                                      : NULL;
514 <            photonPartition [emap.src -> so -> otype] (&emap);
486 >            double numEmit;
487              
488 <            if (photonRepTime) {
489 <               if (!passCnt)
490 <                  sprintf(errmsg, "PREPASS %d on source %s ",
491 <                          prePassCnt, source [srcIdx].so -> oname);
492 <               else
493 <                  sprintf(errmsg, "MAIN PASS on source %s ",
494 <                          source [srcIdx].so -> oname);
495 <                      
496 <               if (emap.port) {
497 <                  sprintf(errmsg2, "via port %s ",
498 <                          photonPorts [portCnt].so -> oname);
499 <                  strcat(errmsg, errmsg2);
488 >            if (!passCnt) {  
489 >               /* INIT PASS 1 */              
490 >               /* Skip if no photons contributed after sufficient
491 >                * iterations; make it clear to user which photon maps are
492 >                * missing so (s)he can check geometry and materials */
493 >               if (++prePassCnt > maxPreDistrib) {
494 >                  sprintf(errmsg, "proc %d: too many prepasses", proc);
495 >
496 >                  for (t = 0; t < NUM_PMAP_TYPES; t++)
497 >                     if (pmaps [t] && !pmaps [t] -> numPhotons) {
498 >                        sprintf(errmsg2, ", no %s photons stored",
499 >                                pmapName [t]);
500 >                        strcat(errmsg, errmsg2);
501 >                     }
502 >                  
503 >                  error(USER, errmsg);
504 >                  break;
505                 }
506 +
507 +               /* Num to emit is fraction of minimum target count */
508 +               numEmit = FHUGE;
509                
510 <               sprintf(errmsg2, "(%lu partitions)...\n", emap.numPartitions);
511 <               strcat(errmsg, errmsg2);
512 <               eputs(errmsg);
513 <               fflush(stderr);
510 >               for (t = 0; t < NUM_PMAP_TYPES; t++)
511 >                  if (pmaps [t])
512 >                     numEmit = min(pmaps [t] -> distribTarget, numEmit);
513 >                    
514 >               numEmit *= preDistrib;
515              }
516 <            
517 <            for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
518 <                 emap.partitionCnt++) {            
519 <               double partNumEmit;
520 <               unsigned long partEmitCnt;
516 >            else {            
517 >               /* INIT PASS 2 */
518 >               /* Based on the outcome of the predistribution we can now
519 >                * estimate how many more photons we have to emit for each
520 >                * photon map to meet its respective target count.  This
521 >                * value is clamped to 0 in case the target has already been
522 >                * exceeded in the pass 1. */
523 >               double maxDistribRatio = 0;
524 >
525 >               /* Set the distribution ratio for each map; this indicates
526 >                * how many photons of each respective type are stored per
527 >                * emitted photon, and is used as probability for storing a
528 >                * photon by newPhoton().  Since this biases the photon
529 >                * density, newPhoton() promotes the flux of stored photons
530 >                * to compensate.  */
531 >               for (t = 0; t < NUM_PMAP_TYPES; t++)
532 >                  if ((pm = pmaps [t])) {
533 >                     pm -> distribRatio = (double)pm -> distribTarget /
534 >                                          pm -> numPhotons - 1;
535 >
536 >                     /* Check if photon map "overflowed", i.e. exceeded its
537 >                      * target count in the prepass; correcting the photon
538 >                      * flux via the distribution ratio is no longer
539 >                      * possible, as no more photons of this type will be
540 >                      * stored, so notify the user rather than deliver
541 >                      * incorrect results.  In future we should handle this
542 >                      * more intelligently by using the photonFlux in each
543 >                      * photon map to individually correct the flux after
544 >                      * distribution.  */
545 >                     if (pm -> distribRatio <= FTINY) {
546 >                        sprintf(errmsg, "%s photon map overflow in "
547 >                                "prepass, reduce -apD", pmapName [t]);
548 >                        error(INTERNAL, errmsg);
549 >                     }
550 >                        
551 >                     maxDistribRatio = max(pm -> distribRatio,
552 >                                           maxDistribRatio);
553 >                  }
554                
555 <               /* Get photon origin within current source partishunn and
556 <                * build emission map */
557 <               photonOrigin [emap.src -> so -> otype] (&emap);
558 <               initPhotonEmission(&emap, pdfSamples);
559 <              
560 <               /* Number of photons to emit from ziss partishunn --
561 <                * proportional to flux; photon ray weight and scalar flux
562 <                * are uniform (the latter only varying in RGB). */
563 <               partNumEmit = numEmit * colorAvg(emap.partFlux) / totalFlux;
564 <               partEmitCnt = (unsigned long)partNumEmit;
551 <              
552 <               /* Probabilistically account for fractional photons */
553 <               if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
554 <                  partEmitCnt++;
555 >               /* Normalise distribution ratios and calculate number of
556 >                * photons to emit in main pass */
557 >               for (t = 0; t < NUM_PMAP_TYPES; t++)
558 >                  if ((pm = pmaps [t]))
559 >                     pm -> distribRatio /= maxDistribRatio;
560 >                    
561 >               if ((numEmit = localNumEmitted * maxDistribRatio) < FTINY)
562 >                  /* No photons left to distribute in main pass */
563 >                  break;
564 >            }
565  
566 <               /* Integer counter avoids FP rounding errors */
567 <               while (partEmitCnt--) {
568 <                  RAY photonRay;
566 >            /* Update shared completion counter for progreport by parent */
567 >            photonCnt -> numComplete += numEmit;                            
568 >
569 >            /* PHOTON DISTRIBUTION LOOP */
570 >            for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
571 >               unsigned portCnt = 0;
572 >               emap.src = source + srcIdx;
573 >
574 >               do {  /* Need at least one iteration if no ports! */
575 >                  emap.port = emap.src -> sflags & SDISTANT
576 >                              ? photonPorts + portCnt : NULL;
577 >                  photonPartition [emap.src -> so -> otype] (&emap);
578 >
579 >                  if (verbose && !proc) {
580 >                     /* Output from subproc 0 only to avoid race condition
581 >                      * on console I/O */
582 >                     if (!passCnt)
583 >                        sprintf(errmsg, "\tPREPASS %d on source %s ",
584 >                                prePassCnt, source [srcIdx].so -> oname);
585 >                     else
586 >                        sprintf(errmsg, "\tMAIN PASS on source %s ",
587 >                                source [srcIdx].so -> oname);
588 >
589 >                     if (emap.port) {
590 >                        sprintf(errmsg2, "via port %s ",
591 >                                photonPorts [portCnt].so -> oname);
592 >                        strcat(errmsg, errmsg2);
593 >                     }
594 >
595 >                     sprintf(errmsg2, "(%lu partitions)\n",
596 >                             emap.numPartitions);
597 >                     strcat(errmsg, errmsg2);
598 >                     eputs(errmsg);
599 > #if NIX                    
600 >                     fflush(stderr);
601 > #endif                    
602 >                  }
603                    
604 <                  /* Emit photon based on PDF and trace through scene until
605 <                   * absorbed/leaked */
606 <                  emitPhoton(&emap, &photonRay);
607 <                  tracePhoton(&photonRay);
604 >                  for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
605 >                       emap.partitionCnt++) {
606 >                     double partNumEmit;
607 >                     unsigned long partEmitCnt;
608 >                    
609 >                     /* Get photon origin within current source partishunn
610 >                      * and build emission map */
611 >                     photonOrigin [emap.src -> so -> otype] (&emap);
612 >                     initPhotonEmission(&emap, pdfSamples);
613 >                    
614 >                     /* Number of photons to emit from ziss partishunn --
615 >                      * proportional to flux; photon ray weight and scalar
616 >                      * flux are uniform (latter only varying in RGB). */
617 >                     partNumEmit = numEmit * colorAvg(emap.partFlux) /
618 >                                   totalFlux;
619 >                     partEmitCnt = (unsigned long)partNumEmit;
620 >                    
621 >                     /* Probabilistically account for fractional photons */
622 >                     if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
623 >                        partEmitCnt++;
624 >
625 >                     /* Update local and shared (global) emission counter */
626 >                     photonCnt -> numEmitted += partEmitCnt;
627 >                     localNumEmitted += partEmitCnt;
628 >
629 >                     /* Integer counter avoids FP rounding errors during
630 >                      * iteration */
631 >                     while (partEmitCnt--) {
632 >                        RAY photonRay;
633 >                        
634 >                        /* Emit photon based on PDF and trace through scene
635 >                         * until absorbed/leaked */
636 >                        emitPhoton(&emap, &photonRay);
637 > #if 1
638 >                        if (emap.port)
639 >                           /* !!!  PHOTON PORT REJECTION SAMPLING HACK: set
640 >                            * !!!  photon port as fake hit object for
641 >                            * !!!  primary ray to check for intersection in
642 >                            * !!!  tracePhoton() */
643 >                           photonRay.ro = emap.port -> so;
644 > #endif
645 >                        tracePhoton(&photonRay);
646 >                     }                                          
647 >
648 >                     /* Update shared global photon count for each pmap */
649 >                     for (t = 0; t < NUM_PMAP_TYPES; t++)
650 >                        if (pmaps [t]) {
651 >                           photonCnt -> numPhotons [t] +=
652 >                              pmaps [t] -> numPhotons - lastNumPhotons [t];
653 >                           lastNumPhotons [t] = pmaps [t] -> numPhotons;
654 >                        }
655 > #if !NIX
656 >                     /* Synchronous progress report on Windoze */
657 >                     if (!proc && photonRepTime > 0 &&
658 >                           time(NULL) >= repLastTime + photonRepTime) {
659 >                        repEmitted = repProgress = photonCnt -> numEmitted;
660 >                        repComplete = photonCnt -> numComplete;                          
661 >                        pmapDistribReport();
662 >                     }
663 > #endif
664 >                  }
665                    
666 <                  /* Record progress */
667 <                  repProgress++;
668 <                  
669 <                  if (photonRepTime > 0 &&
670 <                      time(NULL) >= repLastTime + photonRepTime)
671 <                     pmapDistribReport();
672 <                  #ifdef SIGCONT
673 <                     else signal(SIGCONT, pmapDistribReport);
674 <                  #endif
666 >                  portCnt++;
667 >               } while (portCnt < numPhotonPorts);
668 >            }
669 >            
670 >            for (t = 0; t < NUM_PMAP_TYPES; t++)
671 >               if (pmaps [t] && !pmaps [t] -> numPhotons) {
672 >                  /* Double preDistrib in case a photon map is empty and
673 >                   * redo pass 1 --> possibility of infinite loop for
674 >                   * pathological scenes (e.g.  absorbing materials) */
675 >                  preDistrib *= 2;
676 >                  break;
677                 }
678 +            
679 +            if (t >= NUM_PMAP_TYPES)
680 +               /* No empty photon maps found; now do pass 2 */
681 +               passCnt++;
682 +         } while (passCnt < 2);
683 +        
684 +         /* Flush heap buffa for every pmap one final time;
685 +          * avoids potential data corruption! */
686 +         for (t = 0; t < NUM_PMAP_TYPES; t++)
687 +            if (pmaps [t]) {
688 +               flushPhotonHeap(pmaps [t]);
689 +               /* Heap file closed automatically on exit
690 +                  fclose(pmaps [t] -> heap); */
691 + #ifdef DEBUG_PMAP              
692 +               sprintf(errmsg, "Proc %d: total %ld photons\n", proc,
693 +                       pmaps [t] -> numPhotons);
694 +               eputs(errmsg);
695 + #endif              
696              }
697 <                        
698 <            portCnt++;
699 <         } while (portCnt < numPhotonPorts);
697 > #if NIX
698 >         /* Terminate subprocess */
699 >         exit(0);
700 > #endif
701        }
702 +      else if (pid < 0)
703 +         error(SYSTEM, "failed to fork subprocess in distribPhotons");        
704 +   }
705 +
706 + #if NIX
707 +   /* PARENT PROCESS CONTINUES HERE */
708 + #ifdef SIGCONT
709 +   /* Enable progress report signal handler */
710 +   signal(SIGCONT, pmapDistribReport);
711 + #endif  
712 +   /* Wait for subprocesses complete while reporting progress */
713 +   proc = numProc;
714 +   while (proc) {
715 +      while (waitpid(-1, &stat, WNOHANG) > 0) {
716 +         /* Subprocess exited; check status */
717 +         if (!WIFEXITED(stat) || WEXITSTATUS(stat))
718 +            error(USER, "failed photon distribution");
719        
720 +         --proc;
721 +      }
722 +      
723 +      /* Nod off for a bit and update progress  */
724 +      sleep(1);
725 +
726 +      /* Asynchronous progress report from shared subprocess counters */  
727 +      repEmitted = repProgress = photonCnt -> numEmitted;
728 +      repComplete = photonCnt -> numComplete;      
729 +
730 +      repProgress = repComplete = 0;
731        for (t = 0; t < NUM_PMAP_TYPES; t++)
732 <         if (photonMaps [t] && !photonMaps [t] -> heapEnd) {
733 <            /* Double preDistrib in case a photon map is empty and redo
734 <             * pass 1 --> possibility of infinite loop for pathological
735 <             * scenes (e.g. absorbing materials) */
586 <            preDistrib *= 2;
587 <            break;
732 >         if ((pm = pmaps [t])) {
733 >            /* Get global photon count from shmem updated by subprocs */
734 >            repProgress += pm -> numPhotons = photonCnt -> numPhotons [t];
735 >            repComplete += pm -> distribTarget;
736           }
737 <      
590 <      if (t >= NUM_PMAP_TYPES) {
591 <         /* No empty photon maps found; now do pass 2 */
592 <         passCnt++;
593 <         if (photonRepTime)
594 <            eputs("\n");
595 <      }
596 <   } while (passCnt < 2);
737 >      repComplete *= numProc;
738  
739 +      if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
740 +         pmapDistribReport();
741 + #ifdef SIGCONT
742 +      else signal(SIGCONT, pmapDistribReport);
743 + #endif
744 +   }
745 + #endif /* NIX */
746 +
747     /* ===================================================================
748 <    * POST-DISTRIBUTION - Set photon flux and build kd-tree, etc.
748 >    * POST-DISTRIBUTION - Set photon flux and build data struct for photon
749 >    * storage, etc.
750      * =================================================================== */
751 <   #ifdef SIGCONT    
752 <      signal(SIGCONT, SIG_DFL);
753 <   #endif
751 > #ifdef SIGCONT    
752 >   /* Reset signal handler */
753 >   signal(SIGCONT, SIG_DFL);
754 > #endif
755     free(emap.samples);
756    
757 <   /* Set photon flux (repProgress is total num emitted) */
758 <   totalFlux /= repProgress;
759 <  
757 >   /* Set photon flux */
758 >   totalFlux /= photonCnt -> numEmitted;
759 > #if NIX  
760 >   /* Photon counters no longer needed, unmap shared memory */
761 >   munmap(photonCnt, sizeof(*photonCnt));
762 >   close(shmFile);
763 >   unlink(shmFname);
764 > #else
765 >   free(photonCnt);  
766 > #endif      
767 >   if (verbose)
768 >      eputs("\n");
769 >      
770     for (t = 0; t < NUM_PMAP_TYPES; t++)
771 <      if (photonMaps [t]) {
772 <         if (photonRepTime) {
773 <            sprintf(errmsg, "\nBuilding %s photon map...\n", pmapName [t]);
771 >      if (pmaps [t]) {
772 >         if (verbose) {
773 >            sprintf(errmsg, "Building %s photon map\n", pmapName [t]);
774              eputs(errmsg);
775 + #if NIX            
776              fflush(stderr);
777 + #endif            
778           }
779 <      
780 <         balancePhotons(photonMaps [t], &totalFlux);
779 >        
780 >         /* Build underlying data structure; heap is destroyed */
781 >         buildPhotonMap(pmaps [t], &totalFlux, NULL, numProc);
782        }
783        
784     /* Precompute photon irradiance if necessary */
785 <   if (preCompPmap)
785 >   if (preCompPmap) {
786 >      if (verbose)
787 >         eputs("\n");
788        preComputeGlobal(preCompPmap);
789 < }
624 <
625 <
626 <
627 < void photonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
628 < /* Photon density estimate. Returns irradiance at ray -> rop. */
629 < {
630 <   unsigned i;
631 <   PhotonSQNode *sq;
632 <   float r;
633 <   COLOR flux;
634 <
635 <   setcolor(irrad, 0, 0, 0);
636 <
637 <   if (!pmap -> maxGather)
638 <      return;
639 <      
640 <   /* Ignore sources */
641 <   if (ray -> ro)
642 <      if (islight(objptr(ray -> ro -> omod) -> otype))
643 <         return;
644 <        
645 <   pmap -> squeueEnd = 0;
646 <   findPhotons(pmap, ray);
789 >   }      
790    
791 <   /* Need at least 2 photons */
792 <   if (pmap -> squeueEnd < 2) {
650 <      #ifdef PMAP_NONEFOUND  
651 <         sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)",
652 <                 ray -> ro ? ray -> ro -> oname : "<null>",
653 <                 ray -> rop [0], ray -> rop [1], ray -> rop [2]);
654 <         error(WARNING, errmsg);
655 <      #endif      
656 <
657 <      return;
658 <   }
659 <      
660 <   if (pmap -> minGather == pmap -> maxGather) {
661 <      /* No bias compensation. Just do a plain vanilla estimate */
662 <      sq = pmap -> squeue + 1;
663 <      
664 <      /* Average radius between furthest two photons to improve accuracy */      
665 <      r = max(sq -> dist, (sq + 1) -> dist);
666 <      r = 0.25 * (pmap -> maxDist + r + 2 * sqrt(pmap -> maxDist * r));  
667 <      
668 <      /* Skip the extra photon */
669 <      for (i = 1 ; i < pmap -> squeueEnd; i++, sq++) {
670 <         getPhotonFlux(sq -> photon, flux);        
671 < #ifdef PMAP_EPANECHNIKOV
672 <         /* Apply Epanechnikov kernel to photon flux (dists are squared) */
673 <         scalecolor(flux, 2 * (1 - sq -> dist / r));
674 < #endif        
675 <         addcolor(irrad, flux);
676 <      }
677 <      
678 <      /* Divide by search area PI * r^2, 1 / PI required as ambient
679 <         normalisation factor */        
680 <      scalecolor(irrad, 1 / (PI * PI * r));
681 <      
682 <      return;
683 <   }
684 <   else
685 <      /* Apply bias compensation to density estimate */
686 <      biasComp(pmap, irrad);
687 < }
688 <
689 <
690 <
691 < void photonPreCompDensity (PhotonMap *pmap, RAY *r, COLOR irrad)
692 < /* Returns precomputed photon density estimate at ray -> rop. */
693 < {
694 <   Photon *p;
695 <  
696 <   setcolor(irrad, 0, 0, 0);
697 <
698 <   /* Ignore sources */
699 <   if (r -> ro && islight(objptr(r -> ro -> omod) -> otype))
700 <      return;
701 <      
702 <   if ((p = find1Photon(preCompPmap, r)))
703 <      getPhotonFlux(p, irrad);
704 < }
705 <
706 <
707 <
708 < void volumePhotonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
709 < /* Photon volume density estimate. Returns irradiance at ray -> rop. */
710 < {
711 <   unsigned i;
712 <   PhotonSQNode *sq;
713 <   float gecc2, r, ph;
714 <   COLOR flux;
715 <
716 <   setcolor(irrad, 0, 0, 0);
717 <  
718 <   if (!pmap -> maxGather)
719 <      return;
720 <      
721 <   pmap -> squeueEnd = 0;
722 <   findPhotons(pmap, ray);
723 <  
724 <   /* Need at least 2 photons */
725 <   if (pmap -> squeueEnd < 2)
726 <      return;
727 <      
728 <   if (pmap -> minGather == pmap -> maxGather) {
729 <      /* No bias compensation. Just do a plain vanilla estimate */
730 <      gecc2 = ray -> gecc * ray -> gecc;
731 <      sq = pmap -> squeue + 1;
732 <      
733 <      /* Average radius between furthest two photons to improve accuracy */      
734 <      r = max(sq -> dist, (sq + 1) -> dist);
735 <      r = 0.25 * (pmap -> maxDist + r + 2 * sqrt(pmap -> maxDist * r));  
736 <      
737 <      /* Skip the extra photon */
738 <      for (i = 1 ; i < pmap -> squeueEnd; i++, sq++) {
739 <         /* Compute phase function for inscattering from photon */
740 <         if (gecc2 <= FTINY)
741 <            ph = 1;
742 <         else {
743 <            ph = DOT(ray -> rdir, sq -> photon -> norm) / 127;
744 <            ph = 1 + gecc2 - 2 * ray -> gecc * ph;
745 <            ph = (1 - gecc2) / (ph * sqrt(ph));
746 <         }
747 <        
748 <         getPhotonFlux(sq -> photon, flux);
749 <         scalecolor(flux, ph);
750 <         addcolor(irrad, flux);
751 <      }
752 <      
753 <      /* Divide by search volume 4 / 3 * PI * r^3 and phase function
754 <         normalization factor 1 / (4 * PI) */
755 <      scalecolor(irrad, 3 / (16 * PI * PI * r * sqrt(r)));
756 <      
757 <      return;
758 <   }
759 <  
760 <   else
761 <      /* Apply bias compensation to density estimate */
762 <      volumeBiasComp(pmap, ray, irrad);
791 >   if (verbose)
792 >      eputs("\n");
793   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines