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

Diff Legend

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