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

Diff Legend

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