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.10 by greg, Tue Sep 1 16:27:52 2015 UTC vs.
Revision 2.11 by rschregle, Tue May 17 17:39:47 2016 UTC

# Line 1 | Line 1
1   #ifndef lint
2   static const char RCSid[] = "$Id$";
3   #endif
4 +
5   /*
6 <   ==================================================================
6 >   ======================================================================
7     Photon map main module
8  
9     Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
10     (c) Fraunhofer Institute for Solar Energy Systems,
11     (c) Lucerne University of Applied Sciences and Arts,
12 <   supported by the Swiss National Science Foundation (SNSF, #147053)
13 <   ==================================================================
12 >       supported by the Swiss National Science Foundation (SNSF, #147053)
13 >   ======================================================================
14    
15 +   $Id$
16   */
17  
18  
# Line 25 | Line 27 | static const char RCSid[] = "$Id$";
27   #include "otypes.h"
28   #include <time.h>
29   #include <sys/stat.h>
30 + #include <sys/mman.h>
31 + #include <sys/wait.h>
32  
33 + #define PMAP_REV  "$Revision$"
34  
35  
36   extern char *octname;
37  
33 – static char PmapRevision [] = "$Revision$";
38  
39  
36 –
40   /* Photon map lookup functions per type */
41   void (*pmapLookup [NUM_PMAP_TYPES])(PhotonMap*, RAY*, COLOR) = {
42     photonDensity, photonPreCompDensity, photonDensity, volumePhotonDensity,
# Line 88 | Line 91 | void loadPmaps (PhotonMap **pmaps, const PhotonMapPara
91           pmaps [type] = pm;
92          
93           /* Check for invalid density estimate bandwidth */                            
94 <         if (pm -> maxGather > pm -> heapSize) {
94 >         if (pm -> maxGather > pm -> numPhotons) {
95              error(WARNING, "adjusting density estimate bandwidth");
96 <            pm -> minGather = pm -> maxGather = pm -> heapSize;
96 >            pm -> minGather = pm -> maxGather = pm -> numPhotons;
97           }
98        }
99   }
# Line 147 | Line 150 | static int photonParticipate (RAY *ray)
150        colorNorm(ray -> rcol);
151        VCOPY(ray -> rorg, ray -> rop);
152        
153 <      if (albedo > FTINY)
153 >      if (albedo > FTINY && ray -> rlvl > 0)
154           /* Add to volume photon map */
155 <         if (ray -> rlvl > 0) addPhoton(volumePmap, ray);
155 >         newPhoton(volumePmap, ray);
156          
157        /* Absorbed? */
158 <      if (pmapRandom(rouletteState) > albedo) return 0;
158 >      if (pmapRandom(rouletteState) > albedo)
159 >         return 0;
160        
161        /* Colour bleeding without attenuation (?) */
162        multcolor(ray -> rcol, ray -> albedo);
# Line 237 | Line 241 | void tracePhoton (RAY *ray)
241  
242  
243   static void preComputeGlobal (PhotonMap *pmap)
244 < /* Precompute irradiance from global photons for final gathering using
245 <   the first finalGather * pmap -> heapSize photons in the heap. Returns
246 <   new heap with precomputed photons. */
244 > /* Precompute irradiance from global photons for final gathering for  
245 >   a random subset of finalGather * pmap -> numPhotons photons, and builds
246 >   the photon map, discarding the original photons. */
247 > /* !!! NOTE: PRECOMPUTATION WITH OOC CURRENTLY WITHOUT CACHE !!! */  
248   {
249 <   unsigned long i, nuHeapSize;
250 <   unsigned j;
251 <   Photon *nuHeap, *p;
252 <   COLOR irrad;
253 <   RAY ray;
254 <   float nuMinPos [3], nuMaxPos [3];
249 >   unsigned long  i, numPreComp;
250 >   unsigned       j;
251 >   PhotonIdx      pIdx;
252 >   Photon         photon;
253 >   RAY            ray;
254 >   PhotonMap      nuPmap;
255  
256 <   repComplete = nuHeapSize = finalGather * pmap -> heapSize;
256 >   repComplete = numPreComp = finalGather * pmap -> numPhotons;
257    
258     if (photonRepTime) {
259 <      sprintf(errmsg,
260 <              "Precomputing irradiance for %ld global photons...\n",
256 <              nuHeapSize);
259 >      sprintf(errmsg, "Precomputing irradiance for %ld global photons...\n",
260 >              numPreComp);
261        eputs(errmsg);
262        fflush(stderr);
263     }
264    
265 <   p = nuHeap = (Photon*)malloc(nuHeapSize * sizeof(Photon));
266 <   if (!nuHeap)
267 <      error(USER, "can't allocate photon heap");
268 <      
269 <   for (j = 0; j <= 2; j++) {
270 <      nuMinPos [j] = FHUGE;
267 <      nuMaxPos [j] = -FHUGE;
268 <   }
265 >   /* Copy photon map for precomputed photons */
266 >   memcpy(&nuPmap, pmap, sizeof(PhotonMap));
267 >
268 >   /* Zero counters, init new heap and extents */  
269 >   nuPmap.numPhotons = 0;  
270 >   initPhotonHeap(&nuPmap);
271    
272 +   for (j = 0; j < 3; j++) {  
273 +      nuPmap.minPos [j] = FHUGE;
274 +      nuPmap.maxPos [j] = -FHUGE;
275 +   }
276 +
277     /* Record start time, baby */
278     repStartTime = time(NULL);
279 <   #ifdef SIGCONT
280 <      signal(SIGCONT, pmapPreCompReport);
281 <   #endif
279 > #ifdef SIGCONT
280 >   signal(SIGCONT, pmapPreCompReport);
281 > #endif
282     repProgress = 0;
276 –   memcpy(nuHeap, pmap -> heap, nuHeapSize * sizeof(Photon));
283    
284 <   for (i = 0, p = nuHeap; i < nuHeapSize; i++, p++) {
285 <      ray.ro = NULL;
286 <      VCOPY(ray.rop, p -> pos);
284 >   photonRay(NULL, &ray, PRIMARY, NULL);
285 >   ray.ro = NULL;
286 >  
287 >   for (i = 0; i < numPreComp; i++) {
288 >      /* Get random photon from stratified distribution in source heap to
289 >       * avoid duplicates and clutering */
290 >      pIdx = firstPhoton(pmap) +
291 >             (unsigned long)((i + pmapRandom(pmap -> randState)) /
292 >                             finalGather);
293 >      getPhoton(pmap, pIdx, &photon);
294        
295 <      /* Update min and max positions & set ray normal */
296 <      for (j = 0; j < 3; j++) {
297 <         if (p -> pos [j] < nuMinPos [j]) nuMinPos [j] = p -> pos [j];
298 <         if (p -> pos [j] > nuMaxPos [j]) nuMaxPos [j] = p -> pos [j];
286 <         ray.ron [j] = p -> norm [j] / 127.0;
287 <      }
295 >      /* Init dummy photon ray with intersection at photon position */
296 >      VCOPY(ray.rop, photon.pos);
297 >      for (j = 0; j < 3; j++)
298 >         ray.ron [j] = photon.norm [j] / 127.0;
299        
300 <      photonDensity(pmap, &ray, irrad);
301 <      setPhotonFlux(p, irrad);
300 >      /* Get density estimate at photon position */
301 >      photonDensity(pmap, &ray, ray.rcol);
302 >                  
303 >      /* Append photon to new heap from ray */
304 >      newPhoton(&nuPmap, &ray);
305 >      
306 >      /* Update progress */
307        repProgress++;
308        
309        if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
310           pmapPreCompReport();
311 <      #ifdef SIGCONT
312 <         else signal(SIGCONT, pmapPreCompReport);
313 <      #endif
311 > #ifdef SIGCONT
312 >      else signal(SIGCONT, pmapPreCompReport);
313 > #endif
314     }
315    
316 <   #ifdef SIGCONT  
317 <      signal(SIGCONT, SIG_DFL);
302 <   #endif
316 >   /* Flush heap */
317 >   flushPhotonHeap(&nuPmap);
318    
319 <   /* Replace & rebuild heap */
320 <   free(pmap -> heap);
321 <   pmap -> heap = nuHeap;
307 <   pmap -> heapSize = pmap -> heapEnd = nuHeapSize;
308 <   VCOPY(pmap -> minPos, nuMinPos);
309 <   VCOPY(pmap -> maxPos, nuMaxPos);
319 > #ifdef SIGCONT  
320 >   signal(SIGCONT, SIG_DFL);
321 > #endif
322    
323 +   /* Trash original pmap, replace with precomputed one */
324 +   deletePhotons(pmap);
325 +   memcpy(pmap, &nuPmap, sizeof(PhotonMap));
326 +  
327     if (photonRepTime) {
328 <      eputs("Rebuilding global photon heap...\n");
328 >      eputs("Rebuilding precomputed photon map...\n");
329        fflush(stderr);
330     }
331 <  
332 <   balancePhotons(pmap, NULL);
331 >
332 >   /* Rebuild underlying data structure, destroying heap */  
333 >   buildPhotonMap(pmap, NULL, NULL, 1);
334   }
335  
336  
337  
338 < void distribPhotons (PhotonMap **pmaps)
338 > typedef struct {
339 >   unsigned long  numPhotons [NUM_PMAP_TYPES],
340 >                  numEmitted, numComplete;
341 > } PhotonCnt;
342 >
343 >
344 >
345 > void distribPhotons (PhotonMap **pmaps, unsigned numProc)
346   {
347 <   EmissionMap emap;
348 <   char errmsg2 [128];
349 <   unsigned t, srcIdx, passCnt = 0, prePassCnt = 0;
350 <   double totalFlux = 0;
351 <   PhotonMap *pm;
347 >   EmissionMap    emap;
348 >   char           errmsg2 [128], shmFname [255];
349 >   unsigned       t, srcIdx, proc;
350 >   double         totalFlux = 0;
351 >   int            shmFile, stat, pid;
352 >   PhotonMap      *pm;
353 >   PhotonCnt      *photonCnt;
354    
355     for (t = 0; t < NUM_PMAP_TYPES && !pmaps [t]; t++);
356 +  
357     if (t >= NUM_PMAP_TYPES)
358 <      error(USER, "no photon maps defined");
358 >      error(USER, "no photon maps defined in distribPhotons");
359        
360     if (!nsources)
361 <      error(USER, "no light sources");
361 >      error(USER, "no light sources in distribPhotons");
362  
363     /* ===================================================================
364      * INITIALISATION - Set up emission and scattering funcs
# Line 340 | Line 367 | void distribPhotons (PhotonMap **pmaps)
367     emap.maxPartitions = MAXSPART;
368     emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1);
369     if (!emap.partitions)
370 <      error(INTERNAL, "can't allocate source partitions");
370 >      error(INTERNAL, "can't allocate source partitions in distribPhotons");
371        
372     /* Initialise all defined photon maps */
373     for (t = 0; t < NUM_PMAP_TYPES; t++)
374 <      initPhotonMap(pmaps [t], t);
374 >      if (pmaps [t]) {
375 >         initPhotonMap(pmaps [t], t);
376 >         /* Open photon heapfile */
377 >         initPhotonHeap(pmaps [t]);
378 >         /* Per-subprocess target count */
379 >         pmaps [t] -> distribTarget /= numProc;
380 >      }
381  
382     initPhotonEmissionFuncs();
383     initPhotonScatterFuncs();
# Line 356 | Line 389 | void distribPhotons (PhotonMap **pmaps)
389     /* Get photon sensor modifiers */
390     getPhotonSensors(photonSensorList);
391    
392 <   /* Seed RNGs for photon distribution */
393 <   pmapSeed(randSeed, partState);
394 <   pmapSeed(randSeed, emitState);
395 <   pmapSeed(randSeed, cntState);
396 <   pmapSeed(randSeed, mediumState);
397 <   pmapSeed(randSeed, scatterState);
398 <   pmapSeed(randSeed, rouletteState);
399 <  
392 >   /* Set up shared mem for photon counters (zeroed by ftruncate) */
393 > #if 0  
394 >   snprintf(shmFname, 255, PMAP_SHMFNAME, getpid());
395 >   shmFile = shm_open(shmFname, O_CREAT | O_RDWR, S_IRUSR | S_IWUSR);
396 > #else
397 >   strcpy(shmFname, PMAP_SHMFNAME);
398 >   shmFile = mkstemp(shmFname);
399 > #endif      
400 >
401 >   if (shmFile < 0)
402 >      error(SYSTEM, "failed opening shared memory file in distribPhotons");
403 >
404 >   if (ftruncate(shmFile, sizeof(*photonCnt)) < 0)
405 >      error(SYSTEM, "failed setting shared memory size in distribPhotons");
406 >
407 >   photonCnt = mmap(NULL, sizeof(*photonCnt), PROT_READ | PROT_WRITE,
408 >                    MAP_SHARED, shmFile, 0);
409 >                    
410 >   if (photonCnt == MAP_FAILED)
411 >      error(SYSTEM, "failed mapping shared memory in distribPhotons");
412 >
413     if (photonRepTime)
414        eputs("\n");
415    
416     /* ===================================================================
417      * FLUX INTEGRATION - Get total photon flux from light sources
418      * =================================================================== */
419 <   for (srcIdx = 0; srcIdx < nsources; srcIdx++) {        
419 >   for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
420        unsigned portCnt = 0;
421        emap.src = source + srcIdx;
422        
423 <      do {
423 >      do {  /* Need at least one iteration if no ports! */
424           emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt
425                                                     : NULL;
426           photonPartition [emap.src -> so -> otype] (&emap);
# Line 408 | Line 454 | void distribPhotons (PhotonMap **pmaps)
454     if (totalFlux < FTINY)
455        error(USER, "zero flux from light sources");
456  
457 <   /* Record start time and enable progress report signal handler */
458 <   repStartTime = time(NULL);
459 <   #ifdef SIGCONT
460 <      signal(SIGCONT, pmapDistribReport);
461 <   #endif
462 <   repProgress = prePassCnt = 0;
463 <  
464 <   if (photonRepTime)
465 <      eputs("\n");
420 <  
421 <   /* ===================================================================
422 <    * 2-PASS PHOTON DISTRIBUTION
423 <    * Pass 1 (pre):  emit fraction of target photon count
424 <    * Pass 2 (main): based on outcome of pass 1, estimate remaining number
425 <    *                of photons to emit to approximate target count
426 <    * =================================================================== */
427 <   do {
428 <      double numEmit;
429 <      
430 <      if (!passCnt) {  
431 <         /* INIT PASS 1 */
432 <         /* Skip if no photons contributed after sufficient iterations; make
433 <          * it clear to user which photon maps are missing so (s)he can
434 <          * check the scene geometry and materials */
435 <         if (++prePassCnt > maxPreDistrib) {
436 <            sprintf(errmsg, "too many prepasses");
437 <
438 <            for (t = 0; t < NUM_PMAP_TYPES; t++)
439 <               if (pmaps [t] && !pmaps [t] -> heapEnd) {
440 <                  sprintf(errmsg2, ", no %s photons stored", pmapName [t]);
441 <                  strcat(errmsg, errmsg2);
442 <               }
443 <            
444 <            error(USER, errmsg);
445 <            break;
446 <         }
447 <
448 <         /* Num to emit is fraction of minimum target count */
449 <         numEmit = FHUGE;
457 >   /* MAIN LOOP */  
458 >   for (proc = 0; proc < numProc; proc++) {
459 >      if (!(pid = fork())) {
460 >         /* SUBPROCESS ENTERS HERE.
461 >            All opened and memory mapped files are inherited */
462 >         unsigned       passCnt = 0, prePassCnt = 0;
463 >         unsigned long  lastNumPhotons [NUM_PMAP_TYPES];
464 >         unsigned long  localNumEmitted = 0; /* Num photons emitted by this
465 >                                                subprocess alone */
466          
467 <         for (t = 0; t < NUM_PMAP_TYPES; t++)
468 <            if (pmaps [t])
469 <               numEmit = min(pmaps [t] -> distribTarget, numEmit);
470 <              
471 <         numEmit *= preDistrib;
472 <      }
473 <
458 <      else {            
459 <         /* INIT PASS 2 */
460 <         /* Based on the outcome of the predistribution we can now estimate
461 <          * how many more photons we have to emit for each photon map to
462 <          * meet its respective target count. This value is clamped to 0 in
463 <          * case the target has already been exceeded in the pass 1. Note
464 <          * repProgress is the number of photons emitted thus far, while
465 <          * heapEnd is the number of photons stored in each photon map. */
466 <         double maxDistribRatio = 0;
467 <
468 <         /* Set the distribution ratio for each map; this indicates how many
469 <          * photons of each respective type are stored per emitted photon,
470 <          * and is used as probability for storing a photon by addPhoton().
471 <          * Since this biases the photon density, addPhoton() promotes the
472 <          * flux of stored photons to compensate. */
473 <         for (t = 0; t < NUM_PMAP_TYPES; t++)
474 <            if ((pm = pmaps [t])) {
475 <               pm -> distribRatio = (double)pm -> distribTarget /
476 <                                    pm -> heapEnd - 1;
477 <
478 <               /* Check if photon map "overflowed", i.e. exceeded its target
479 <                * count in the prepass; correcting the photon flux via the
480 <                * distribution ratio is no longer possible, as no more
481 <                * photons of this type will be stored, so notify the user
482 <                * rather than deliver incorrect results.
483 <                * In future we should handle this more intelligently by
484 <                * using the photonFlux in each photon map to individually
485 <                * correct the flux after distribution. */
486 <               if (pm -> distribRatio <= FTINY) {
487 <                  sprintf(errmsg,
488 <                          "%s photon map overflow in prepass, reduce -apD",
489 <                          pmapName [t]);
490 <                  error(INTERNAL, errmsg);
491 <               }
467 >         /* Seed RNGs from PID for decorellated photon distribution */
468 >         pmapSeed(randSeed + proc, partState);
469 >         pmapSeed(randSeed + proc, emitState);
470 >         pmapSeed(randSeed + proc, cntState);
471 >         pmapSeed(randSeed + proc, mediumState);
472 >         pmapSeed(randSeed + proc, scatterState);
473 >         pmapSeed(randSeed + proc, rouletteState);
474                    
493 –               maxDistribRatio = max(pm -> distribRatio, maxDistribRatio);
494 –            }
495 –        
496 –         /* Normalise distribution ratios and calculate number of photons to
497 –          * emit in main pass */
475           for (t = 0; t < NUM_PMAP_TYPES; t++)
476 <            if ((pm = pmaps [t]))
477 <               pm -> distribRatio /= maxDistribRatio;
478 <              
479 <         if ((numEmit = repProgress * maxDistribRatio) < FTINY)
480 <            /* No photons left to distribute in main pass */
481 <            break;
482 <      }
483 <      
484 <      /* Set completion count for progress report */
508 <      repComplete = numEmit + repProgress;                            
509 <      
510 <      /* PHOTON DISTRIBUTION LOOP */
511 <      for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
512 <         unsigned portCnt = 0;
513 <         emap.src = source + srcIdx;
514 <
476 >            lastNumPhotons [t] = 0;
477 >            
478 >         /* =============================================================
479 >          * 2-PASS PHOTON DISTRIBUTION
480 >          * Pass 1 (pre):  emit fraction of target photon count
481 >          * Pass 2 (main): based on outcome of pass 1, estimate remaining
482 >          *                number of photons to emit to approximate target
483 >          *                count
484 >          * ============================================================= */
485           do {
486 <            emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt
517 <                                                      : NULL;
518 <            photonPartition [emap.src -> so -> otype] (&emap);
486 >            double numEmit;
487              
488 <            if (photonRepTime) {
489 <               if (!passCnt)
490 <                  sprintf(errmsg, "PREPASS %d on source %s ",
491 <                          prePassCnt, source [srcIdx].so -> oname);
492 <               else
493 <                  sprintf(errmsg, "MAIN PASS on source %s ",
494 <                          source [srcIdx].so -> oname);
495 <                      
496 <               if (emap.port) {
497 <                  sprintf(errmsg2, "via port %s ",
498 <                          photonPorts [portCnt].so -> oname);
499 <                  strcat(errmsg, errmsg2);
488 >            if (!passCnt) {  
489 >               /* INIT PASS 1 */              
490 >               /* Skip if no photons contributed after sufficient
491 >                * iterations; make it clear to user which photon maps are
492 >                * missing so (s)he can check geometry and materials */
493 >               if (++prePassCnt > maxPreDistrib) {
494 >                  sprintf(errmsg,
495 >                          "proc %d, source %s: too many prepasses",
496 >                          proc, source [srcIdx].so -> oname);              
497 >
498 >                  for (t = 0; t < NUM_PMAP_TYPES; t++)
499 >                     if (pmaps [t] && !pmaps [t] -> numPhotons) {
500 >                        sprintf(errmsg2, ", no %s photons stored",
501 >                                pmapName [t]);
502 >                        strcat(errmsg, errmsg2);
503 >                     }
504 >                  
505 >                  error(USER, errmsg);
506 >                  break;
507                 }
508 +
509 +               /* Num to emit is fraction of minimum target count */
510 +               numEmit = FHUGE;
511                
512 <               sprintf(errmsg2, "(%lu partitions)...\n", emap.numPartitions);
513 <               strcat(errmsg, errmsg2);
514 <               eputs(errmsg);
515 <               fflush(stderr);
512 >               for (t = 0; t < NUM_PMAP_TYPES; t++)
513 >                  if (pmaps [t])
514 >                     numEmit = min(pmaps [t] -> distribTarget, numEmit);
515 >                    
516 >               numEmit *= preDistrib;
517              }
518 <            
519 <            for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
520 <                 emap.partitionCnt++) {            
521 <               double partNumEmit;
522 <               unsigned long partEmitCnt;
518 >            else {            
519 >               /* INIT PASS 2 */
520 >               /* Based on the outcome of the predistribution we can now
521 >                * estimate how many more photons we have to emit for each
522 >                * photon map to meet its respective target count.  This
523 >                * value is clamped to 0 in case the target has already been
524 >                * exceeded in the pass 1. */
525 >               double maxDistribRatio = 0;
526 >
527 >               /* Set the distribution ratio for each map; this indicates
528 >                * how many photons of each respective type are stored per
529 >                * emitted photon, and is used as probability for storing a
530 >                * photon by newPhoton().  Since this biases the photon
531 >                * density, newPhoton() promotes the flux of stored photons
532 >                * to compensate.  */
533 >               for (t = 0; t < NUM_PMAP_TYPES; t++)
534 >                  if ((pm = pmaps [t])) {
535 >                     pm -> distribRatio = (double)pm -> distribTarget /
536 >                                          pm -> numPhotons - 1;
537 >
538 >                     /* Check if photon map "overflowed", i.e. exceeded its
539 >                      * target count in the prepass; correcting the photon
540 >                      * flux via the distribution ratio is no longer
541 >                      * possible, as no more photons of this type will be
542 >                      * stored, so notify the user rather than deliver
543 >                      * incorrect results.  In future we should handle this
544 >                      * more intelligently by using the photonFlux in each
545 >                      * photon map to individually correct the flux after
546 >                      * distribution.  */
547 >                     if (pm -> distribRatio <= FTINY) {
548 >                        sprintf(errmsg, "%s photon map overflow in "
549 >                                "prepass, reduce -apD", pmapName [t]);
550 >                        error(INTERNAL, errmsg);
551 >                     }
552 >                        
553 >                     maxDistribRatio = max(pm -> distribRatio,
554 >                                           maxDistribRatio);
555 >                  }
556                
557 <               /* Get photon origin within current source partishunn and
558 <                * build emission map */
559 <               photonOrigin [emap.src -> so -> otype] (&emap);
560 <               initPhotonEmission(&emap, pdfSamples);
561 <              
562 <               /* Number of photons to emit from ziss partishunn --
563 <                * proportional to flux; photon ray weight and scalar flux
564 <                * are uniform (the latter only varying in RGB). */
565 <               partNumEmit = numEmit * colorAvg(emap.partFlux) / totalFlux;
566 <               partEmitCnt = (unsigned long)partNumEmit;
555 <              
556 <               /* Probabilistically account for fractional photons */
557 <               if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
558 <                  partEmitCnt++;
557 >               /* Normalise distribution ratios and calculate number of
558 >                * photons to emit in main pass */
559 >               for (t = 0; t < NUM_PMAP_TYPES; t++)
560 >                  if ((pm = pmaps [t]))
561 >                     pm -> distribRatio /= maxDistribRatio;
562 >                    
563 >               if ((numEmit = localNumEmitted * maxDistribRatio) < FTINY)
564 >                  /* No photons left to distribute in main pass */
565 >                  break;
566 >            }
567  
568 <               /* Integer counter avoids FP rounding errors */
569 <               while (partEmitCnt--) {
570 <                  RAY photonRay;
568 >            /* Update shared completion counter for prog.report by parent */
569 >            photonCnt -> numComplete += numEmit;                            
570 >
571 >            /* PHOTON DISTRIBUTION LOOP */
572 >            for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
573 >               unsigned portCnt = 0;
574 >               emap.src = source + srcIdx;
575 >
576 >               do {  /* Need at least one iteration if no ports! */
577 >                  emap.port = emap.src -> sflags & SDISTANT
578 >                              ? photonPorts + portCnt : NULL;
579 >                  photonPartition [emap.src -> so -> otype] (&emap);
580 >
581 >                  if (photonRepTime && !proc) {
582 >                     if (!passCnt)
583 >                        sprintf(errmsg, "PREPASS %d on source %s ",
584 >                                prePassCnt, source [srcIdx].so -> oname);
585 >                     else
586 >                        sprintf(errmsg, "MAIN PASS on source %s ",
587 >                                source [srcIdx].so -> oname);
588 >                            
589 >                     if (emap.port) {
590 >                        sprintf(errmsg2, "via port %s ",
591 >                                photonPorts [portCnt].so -> oname);
592 >                        strcat(errmsg, errmsg2);
593 >                     }
594 >                    
595 >                     sprintf(errmsg2, "(%lu partitions)\n",
596 >                             emap.numPartitions);
597 >                     strcat(errmsg, errmsg2);
598 >                     eputs(errmsg);
599 >                     fflush(stderr);
600 >                  }
601                    
602 <                  /* Emit photon based on PDF and trace through scene until
603 <                   * absorbed/leaked */
604 <                  emitPhoton(&emap, &photonRay);
605 <                  tracePhoton(&photonRay);
602 >                  for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
603 >                       emap.partitionCnt++) {
604 >                     double partNumEmit;
605 >                     unsigned long partEmitCnt;
606 >                    
607 >                     /* Get photon origin within current source partishunn
608 >                      * and build emission map */
609 >                     photonOrigin [emap.src -> so -> otype] (&emap);
610 >                     initPhotonEmission(&emap, pdfSamples);
611 >                    
612 >                     /* Number of photons to emit from ziss partishunn --
613 >                      * proportional to flux; photon ray weight and scalar
614 >                      * flux are uniform (the latter only varying in RGB).
615 >                      * */
616 >                     partNumEmit = numEmit * colorAvg(emap.partFlux) /
617 >                                   totalFlux;
618 >                     partEmitCnt = (unsigned long)partNumEmit;
619 >                    
620 >                     /* Probabilistically account for fractional photons */
621 >                     if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
622 >                        partEmitCnt++;
623 >
624 >                     /* Update local and shared (global) emission counter */
625 >                     photonCnt -> numEmitted += partEmitCnt;
626 >                     localNumEmitted += partEmitCnt;
627 >
628 >                     /* Integer counter avoids FP rounding errors during
629 >                      * iteration */
630 >                     while (partEmitCnt--) {
631 >                        RAY photonRay;
632 >                        
633 >                        /* Emit photon based on PDF and trace through scene
634 >                         * until absorbed/leaked */
635 >                        emitPhoton(&emap, &photonRay);
636 >                        tracePhoton(&photonRay);
637 >                     }                                          
638 >                    
639 >                     /* Update shared global photon count for each pmap */
640 >                     for (t = 0; t < NUM_PMAP_TYPES; t++)
641 >                        if (pmaps [t]) {
642 >                           photonCnt -> numPhotons [t] +=
643 >                              pmaps [t] -> numPhotons - lastNumPhotons [t];
644 >                           lastNumPhotons [t] = pmaps [t] -> numPhotons;
645 >                        }
646 >                  }
647                    
648 <                  /* Record progress */
649 <                  repProgress++;
650 <                  
651 <                  if (photonRepTime > 0 &&
652 <                      time(NULL) >= repLastTime + photonRepTime)
653 <                     pmapDistribReport();
654 <                  #ifdef SIGCONT
655 <                     else signal(SIGCONT, pmapDistribReport);
656 <                  #endif
648 >                  portCnt++;
649 >               } while (portCnt < numPhotonPorts);
650 >            }
651 >            
652 >            for (t = 0; t < NUM_PMAP_TYPES; t++)
653 >               if (pmaps [t] && !pmaps [t] -> numPhotons) {
654 >                  /* Double preDistrib in case a photon map is empty and
655 >                   * redo pass 1 --> possibility of infinite loop for
656 >                   * pathological scenes (e.g.  absorbing materials) */
657 >                  preDistrib *= 2;
658 >                  break;
659                 }
660 +            
661 +            if (t >= NUM_PMAP_TYPES) {
662 +               /* No empty photon maps found; now do pass 2 */
663 +               passCnt++;
664 + #if 0
665 +               if (photonRepTime)
666 +                  eputs("\n");
667 + #endif
668              }
669 <                        
670 <            portCnt++;
671 <         } while (portCnt < numPhotonPorts);
669 >         } while (passCnt < 2);
670 >        
671 >         /* Unmap shared photon counters */
672 > #if 0        
673 >         munmap(photonCnt, sizeof(*photonCnt));
674 >         close(shmFile);
675 > #endif
676 >        
677 >         /* Flush heap buffa for every pmap one final time; this is required
678 >          * to prevent data corruption! */
679 >         for (t = 0; t < NUM_PMAP_TYPES; t++)
680 >            if (pmaps [t]) {
681 > #if 0            
682 >               eputs("Final flush\n");
683 > #endif              
684 >               flushPhotonHeap(pmaps [t]);
685 >               fclose(pmaps [t] -> heap);
686 > #ifdef DEBUG_PMAP              
687 >               sprintf(errmsg, "Proc %d: total %ld photons\n", getpid(),
688 >                       pmaps [t] -> numPhotons);
689 >               eputs(errmsg);
690 > #endif              
691 >            }
692 >
693 >         exit(0);
694        }
695 +      else if (pid < 0)
696 +         error(SYSTEM, "failed to fork subprocess in distribPhotons");        
697 +   }
698 +
699 +   /* PARENT PROCESS CONTINUES HERE */
700 +   /* Record start time and enable progress report signal handler */
701 +   repStartTime = time(NULL);
702 + #ifdef SIGCONT
703 +   signal(SIGCONT, pmapDistribReport);
704 + #endif
705 +  
706 +   if (photonRepTime)
707 +      eputs("\n");
708 +  
709 +   /* Wait for subprocesses to complete while reporting progress */
710 +   proc = numProc;
711 +   while (proc) {
712 +      while (waitpid(-1, &stat, WNOHANG) > 0) {
713 +         /* Subprocess exited; check status */
714 +         if (!WIFEXITED(stat) || WEXITSTATUS(stat))
715 +            error(USER, "failed photon distribution");
716        
717 +         --proc;
718 +      }
719 +      
720 +      /* Nod off for a bit and update progress  */
721 +      sleep(1);
722 +      /* Update progress report from shared subprocess counters */
723 +      repEmitted = repProgress = photonCnt -> numEmitted;
724 +      repComplete = photonCnt -> numComplete;
725 +
726        for (t = 0; t < NUM_PMAP_TYPES; t++)
727 <         if (pmaps [t] && !pmaps [t] -> heapEnd) {
728 <            /* Double preDistrib in case a photon map is empty and redo
729 <             * pass 1 --> possibility of infinite loop for pathological
730 <             * scenes (e.g. absorbing materials) */
731 <            preDistrib *= 2;
732 <            break;
727 >         if ((pm = pmaps [t])) {
728 > #if 0        
729 >            /* Get photon count from heapfile size for progress update */
730 >            fseek(pm -> heap, 0, SEEK_END);
731 >            pm -> numPhotons = ftell(pm -> heap) / sizeof(Photon); */
732 > #else            
733 >            /* Get global photon count from shmem updated by subprocs */
734 >            pm -> numPhotons = photonCnt -> numPhotons [t];
735 > #endif            
736           }
593 –      
594 –      if (t >= NUM_PMAP_TYPES) {
595 –         /* No empty photon maps found; now do pass 2 */
596 –         passCnt++;
597 –         if (photonRepTime)
598 –            eputs("\n");
599 –      }
600 –   } while (passCnt < 2);
737  
738 +      if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
739 +         pmapDistribReport();
740 + #ifdef SIGCONT
741 +      else signal(SIGCONT, pmapDistribReport);
742 + #endif
743 +   }
744 +
745     /* ===================================================================
746 <    * POST-DISTRIBUTION - Set photon flux and build kd-tree, etc.
746 >    * POST-DISTRIBUTION - Set photon flux and build data struct for photon
747 >    * storage, etc.
748      * =================================================================== */
749 <   #ifdef SIGCONT    
750 <      signal(SIGCONT, SIG_DFL);
751 <   #endif
749 > #ifdef SIGCONT    
750 >   signal(SIGCONT, SIG_DFL);
751 > #endif
752     free(emap.samples);
753    
754     /* Set photon flux (repProgress is total num emitted) */
755 <   totalFlux /= repProgress;
755 >   totalFlux /= photonCnt -> numEmitted;
756    
757 +   /* Photon counters no longer needed, unmap shared memory */
758 +   munmap(photonCnt, sizeof(*photonCnt));
759 +   close(shmFile);
760 + #if 0  
761 +   shm_unlink(shmFname);
762 + #else
763 +   unlink(shmFname);
764 + #endif      
765 +  
766     for (t = 0; t < NUM_PMAP_TYPES; t++)
767        if (pmaps [t]) {
768           if (photonRepTime) {
# Line 617 | Line 770 | void distribPhotons (PhotonMap **pmaps)
770              eputs(errmsg);
771              fflush(stderr);
772           }
773 <      
774 <         balancePhotons(pmaps [t], &totalFlux);
773 >        
774 >         /* Build underlying data structure; heap is destroyed */
775 >         buildPhotonMap(pmaps [t], &totalFlux, NULL, numProc);
776        }
777 <      
777 >
778     /* Precompute photon irradiance if necessary */
779     if (preCompPmap)
780        preComputeGlobal(preCompPmap);
# Line 631 | Line 785 | void distribPhotons (PhotonMap **pmaps)
785   void photonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
786   /* Photon density estimate. Returns irradiance at ray -> rop. */
787   {
788 <   unsigned i;
789 <   PhotonSQNode *sq;
790 <   float r;
791 <   COLOR flux;
788 >   unsigned                      i;
789 >   float                         r;
790 >   COLOR                         flux;
791 >   Photon                        *photon;
792 >   const PhotonSearchQueueNode   *sqn;
793  
794     setcolor(irrad, 0, 0, 0);
795  
# Line 642 | Line 797 | void photonDensity (PhotonMap *pmap, RAY *ray, COLOR i
797        return;
798        
799     /* Ignore sources */
800 <   if (ray -> ro)
801 <      if (islight(objptr(ray -> ro -> omod) -> otype))
647 <         return;
800 >   if (ray -> ro && islight(objptr(ray -> ro -> omod) -> otype))
801 >      return;
802          
649 –   pmap -> squeueEnd = 0;
803     findPhotons(pmap, ray);
804    
805     /* Need at least 2 photons */
806 <   if (pmap -> squeueEnd < 2) {
807 <      #ifdef PMAP_NONEFOUND  
808 <         sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)",
809 <                 ray -> ro ? ray -> ro -> oname : "<null>",
810 <                 ray -> rop [0], ray -> rop [1], ray -> rop [2]);
811 <         error(WARNING, errmsg);
812 <      #endif      
806 >   if (pmap -> squeue.tail < 2) {
807 > #ifdef PMAP_NONEFOUND  
808 >      sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)",
809 >              ray -> ro ? ray -> ro -> oname : "<null>",
810 >              ray -> rop [0], ray -> rop [1], ray -> rop [2]);
811 >      error(WARNING, errmsg);
812 > #endif      
813  
814        return;
815     }
816 <      
816 >
817     if (pmap -> minGather == pmap -> maxGather) {
818        /* No bias compensation. Just do a plain vanilla estimate */
819 <      sq = pmap -> squeue + 1;
819 >      sqn = pmap -> squeue.node + 1;
820        
821        /* Average radius between furthest two photons to improve accuracy */      
822 <      r = max(sq -> dist, (sq + 1) -> dist);
823 <      r = 0.25 * (pmap -> maxDist + r + 2 * sqrt(pmap -> maxDist * r));  
822 >      r = max(sqn -> dist2, (sqn + 1) -> dist2);
823 >      r = 0.25 * (pmap -> maxDist2 + r + 2 * sqrt(pmap -> maxDist2 * r));  
824        
825        /* Skip the extra photon */
826 <      for (i = 1 ; i < pmap -> squeueEnd; i++, sq++) {
827 <         getPhotonFlux(sq -> photon, flux);        
826 >      for (i = 1 ; i < pmap -> squeue.tail; i++, sqn++) {
827 >         photon = getNearestPhoton(&pmap -> squeue, sqn -> idx);
828 >         getPhotonFlux(photon, flux);        
829   #ifdef PMAP_EPANECHNIKOV
830 <         /* Apply Epanechnikov kernel to photon flux (dists are squared) */
831 <         scalecolor(flux, 2 * (1 - sq -> dist / r));
832 < #endif        
830 >         /* Apply Epanechnikov kernel to photon flux based on photon dist */
831 >         scalecolor(flux, 2 * (1 - sqn -> dist2 / r));
832 > #endif  
833           addcolor(irrad, flux);
834        }
835        
# Line 695 | Line 849 | void photonDensity (PhotonMap *pmap, RAY *ray, COLOR i
849   void photonPreCompDensity (PhotonMap *pmap, RAY *r, COLOR irrad)
850   /* Returns precomputed photon density estimate at ray -> rop. */
851   {
852 <   Photon *p;
852 >   Photon p;
853    
854     setcolor(irrad, 0, 0, 0);
855  
# Line 703 | Line 857 | void photonPreCompDensity (PhotonMap *pmap, RAY *r, CO
857     if (r -> ro && islight(objptr(r -> ro -> omod) -> otype))
858        return;
859        
860 <   if ((p = find1Photon(preCompPmap, r)))
861 <      getPhotonFlux(p, irrad);
860 >   find1Photon(preCompPmap, r, &p);
861 >   getPhotonFlux(&p, irrad);
862   }
863  
864  
# Line 712 | Line 866 | void photonPreCompDensity (PhotonMap *pmap, RAY *r, CO
866   void volumePhotonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
867   /* Photon volume density estimate. Returns irradiance at ray -> rop. */
868   {
869 <   unsigned i;
870 <   PhotonSQNode *sq;
871 <   float gecc2, r, ph;
872 <   COLOR flux;
869 >   unsigned                      i;
870 >   float                         r, gecc2, ph;
871 >   COLOR                         flux;
872 >   Photon                        *photon;
873 >   const PhotonSearchQueueNode   *sqn;
874  
875     setcolor(irrad, 0, 0, 0);
876    
877     if (!pmap -> maxGather)
878        return;
879        
725 –   pmap -> squeueEnd = 0;
880     findPhotons(pmap, ray);
881    
882     /* Need at least 2 photons */
883 <   if (pmap -> squeueEnd < 2)
883 >   if (pmap -> squeue.tail < 2)
884        return;
885 <      
886 <   if (pmap -> minGather == pmap -> maxGather) {
885 >
886 > #if 0      
887 >   /* Volume biascomp disabled (probably redundant) */
888 >   if (pmap -> minGather == pmap -> maxGather)
889 > #endif  
890 >   {
891        /* No bias compensation. Just do a plain vanilla estimate */
892        gecc2 = ray -> gecc * ray -> gecc;
893 <      sq = pmap -> squeue + 1;
893 >      sqn = pmap -> squeue.node + 1;
894        
895        /* Average radius between furthest two photons to improve accuracy */      
896 <      r = max(sq -> dist, (sq + 1) -> dist);
897 <      r = 0.25 * (pmap -> maxDist + r + 2 * sqrt(pmap -> maxDist * r));  
896 >      r = max(sqn -> dist2, (sqn + 1) -> dist2);
897 >      r = 0.25 * (pmap -> maxDist2 + r + 2 * sqrt(pmap -> maxDist2 * r));  
898        
899        /* Skip the extra photon */
900 <      for (i = 1 ; i < pmap -> squeueEnd; i++, sq++) {
900 >      for (i = 1; i < pmap -> squeue.tail; i++, sqn++) {
901 >         photon = getNearestPhoton(&pmap -> squeue, sqn -> idx);
902 >        
903           /* Compute phase function for inscattering from photon */
904           if (gecc2 <= FTINY)
905              ph = 1;
906           else {
907 <            ph = DOT(ray -> rdir, sq -> photon -> norm) / 127;
907 >            ph = DOT(ray -> rdir, photon -> norm) / 127;
908              ph = 1 + gecc2 - 2 * ray -> gecc * ph;
909              ph = (1 - gecc2) / (ph * sqrt(ph));
910           }
911          
912 <         getPhotonFlux(sq -> photon, flux);
912 >         getPhotonFlux(photon, flux);
913           scalecolor(flux, ph);
914           addcolor(irrad, flux);
915        }
# Line 757 | Line 917 | void volumePhotonDensity (PhotonMap *pmap, RAY *ray, C
917        /* Divide by search volume 4 / 3 * PI * r^3 and phase function
918           normalization factor 1 / (4 * PI) */
919        scalecolor(irrad, 3 / (16 * PI * PI * r * sqrt(r)));
760 –      
920        return;
921     }
922 <  
922 > #if 0
923     else
924        /* Apply bias compensation to density estimate */
925        volumeBiasComp(pmap, ray, irrad);
926 + #endif      
927   }

Diff Legend

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