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.11 by rschregle, Tue May 17 17:39:47 2016 UTC vs.
Revision 2.17 by rschregle, Tue Dec 18 22:14:04 2018 UTC

# Line 2 | Line 2
2   static const char RCSid[] = "$Id$";
3   #endif
4  
5 +
6   /*
7     ======================================================================
8     Photon map main module
# Line 16 | Line 17 | static const char RCSid[] = "$Id$";
17   */
18  
19  
19
20   #include "pmap.h"
21   #include "pmapmat.h"
22   #include "pmapsrc.h"
# Line 26 | Line 26 | static const char RCSid[] = "$Id$";
26   #include "pmapdiag.h"
27   #include "otypes.h"
28   #include <time.h>
29 < #include <sys/stat.h>
30 < #include <sys/mman.h>
31 < #include <sys/wait.h>
29 > #if NIX
30 >   #include <sys/stat.h>
31 >   #include <sys/mman.h>
32 >   #include <sys/wait.h>
33 > #endif
34  
33 #define PMAP_REV  "$Revision$"
35  
35
36 extern char *octname;
37
38
39
40 /* Photon map lookup functions per type */
41 void (*pmapLookup [NUM_PMAP_TYPES])(PhotonMap*, RAY*, COLOR) = {
42   photonDensity, photonPreCompDensity, photonDensity, volumePhotonDensity,
43   photonDensity, NULL
44 };
45
46
47
48 void colorNorm (COLOR c)
49 /* Normalise colour channels to average of 1 */
50 {
51   const float avg = colorAvg(c);
52  
53   if (!avg)
54      return;
55      
56   c [0] /= avg;
57   c [1] /= avg;
58   c [2] /= avg;
59 }
60
61
62
63 void loadPmaps (PhotonMap **pmaps, const PhotonMapParams *parm)
64 {
65   unsigned t;
66   struct stat octstat, pmstat;
67   PhotonMap *pm;
68   PhotonMapType type;
69  
70   for (t = 0; t < NUM_PMAP_TYPES; t++)
71      if (setPmapParam(&pm, parm + t)) {        
72         /* Check if photon map newer than octree */
73         if (pm -> fileName && octname &&
74             !stat(pm -> fileName, &pmstat) && !stat(octname, &octstat) &&
75             octstat.st_mtime > pmstat.st_mtime) {
76            sprintf(errmsg, "photon map in file %s may be stale",
77                    pm -> fileName);
78            error(USER, errmsg);
79         }
80        
81         /* Load photon map from file and get its type */
82         if ((type = loadPhotonMap(pm, pm -> fileName)) == PMAP_TYPE_NONE)
83            error(USER, "failed loading photon map");
84            
85         /* Assign to appropriate photon map type (deleting previously
86          * loaded photon map of same type if necessary) */
87         if (pmaps [type]) {
88            deletePhotons(pmaps [type]);
89            free(pmaps [type]);
90         }
91         pmaps [type] = pm;
92        
93         /* Check for invalid density estimate bandwidth */                            
94         if (pm -> maxGather > pm -> numPhotons) {
95            error(WARNING, "adjusting density estimate bandwidth");
96            pm -> minGather = pm -> maxGather = pm -> numPhotons;
97         }
98      }
99 }
100
101
102
36   void savePmaps (const PhotonMap **pmaps, int argc, char **argv)
37   {
38     unsigned t;
# Line 111 | Line 44 | void savePmaps (const PhotonMap **pmaps, int argc, cha
44   }                  
45  
46  
114  
115 void cleanUpPmaps (PhotonMap **pmaps)
116 {
117   unsigned t;
118  
119   for (t = 0; t < NUM_PMAP_TYPES; t++) {
120      if (pmaps [t]) {
121         deletePhotons(pmaps [t]);
122         free(pmaps [t]);
123      }
124   }
125 }
126
127
47      
48   static int photonParticipate (RAY *ray)
49   /* Trace photon through participating medium. Returns 1 if passed through,
50     or 0 if absorbed and $*%&ed. Analogon to rayparticipate(). */
51   {
52     int i;
53 <   RREAL cosTheta, cosPhi, du, dv;
53 >   RREAL xi1, cosTheta, phi, du, dv;
54     const float cext = colorAvg(ray -> cext),
55 <               albedo = colorAvg(ray -> albedo);
55 >               albedo = colorAvg(ray -> albedo),
56 >               gecc = ray -> gecc, gecc2 = sqr(gecc);
57     FVECT u, v;
58     COLOR cvext;
59  
# Line 141 | Line 61 | static int photonParticipate (RAY *ray)
61     ray -> rmax = -log(pmapRandom(mediumState)) / cext;
62    
63     while (!localhit(ray, &thescene)) {
64 +      if (!incube(&thescene, ray -> rop)) {
65 +         /* Terminate photon if it has leaked from the scene */
66 + #ifdef DEBUG_PMAP
67 +         fprintf(stderr,
68 +                 "Volume photon leaked from scene at [%.3f %.3f %.3f]\n",
69 +                 ray -> rop [0], ray -> rop [1], ray -> rop [2]);
70 + #endif                
71 +         return 0;
72 +      }
73 +        
74        setcolor(cvext, exp(-ray -> rmax * ray -> cext [0]),
75                        exp(-ray -> rmax * ray -> cext [1]),
76                        exp(-ray -> rmax * ray -> cext [2]));
# Line 150 | Line 80 | static int photonParticipate (RAY *ray)
80        colorNorm(ray -> rcol);
81        VCOPY(ray -> rorg, ray -> rop);
82        
83 + #if 0
84        if (albedo > FTINY && ray -> rlvl > 0)
85 + #else
86 +      /* Store volume photons unconditionally in mist to also account for
87 +         direct inscattering from sources */
88 +      if (albedo > FTINY)
89 + #endif
90           /* Add to volume photon map */
91           newPhoton(volumePmap, ray);
92          
# Line 163 | Line 99 | static int photonParticipate (RAY *ray)
99        scalecolor(ray -> rcol, 1 / albedo);    
100        
101        /* Scatter photon */
102 <      cosTheta = ray -> gecc <= FTINY ? 2 * pmapRandom(scatterState) - 1
103 <                                      : 1 / (2 * ray -> gecc) *
104 <                                            (1 + ray -> gecc * ray -> gecc -
105 <                                               (1 - ray -> gecc * ray -> gecc) /
106 <                                               (1 - ray -> gecc + 2 * ray -> gecc *
107 <                                                  pmapRandom(scatterState)));
108 <                                                  
109 <      cosPhi = cos(2 * PI * pmapRandom(scatterState));
110 <      du = dv = sqrt(1 - cosTheta * cosTheta);   /* sin(theta) */
111 <      du *= cosPhi;
112 <      dv *= sqrt(1 - cosPhi * cosPhi);           /* sin(phi) */
102 >      xi1 = pmapRandom(scatterState);
103 >      cosTheta = ray -> gecc <= FTINY
104 >                    ? 2 * xi1 - 1
105 >                    : 0.5 / gecc *
106 >                      (1 + gecc2 - sqr((1 - gecc2) /
107 >                                       (1 + gecc * (2 * xi1 - 1))));
108 >
109 >      phi = 2 * PI * pmapRandom(scatterState);
110 >      du = dv = sqrt(1 - sqr(cosTheta));     /* sin(theta) */
111 >      du *= cos(phi);
112 >      dv *= sin(phi);
113        
114        /* Get axes u & v perpendicular to photon direction */
115        i = 0;
# Line 187 | Line 123 | static int photonParticipate (RAY *ray)
123        for (i = 0; i < 3; i++)
124           ray -> rdir [i] = du * u [i] + dv * v [i] +
125                             cosTheta * ray -> rdir [i];
126 +
127        ray -> rlvl++;
128        ray -> rmax = -log(pmapRandom(mediumState)) / cext;
129     }  
130 <    
130 >  
131 >   /* Passed through medium until intersecting local object */  
132     setcolor(cvext, exp(-ray -> rot * ray -> cext [0]),
133                     exp(-ray -> rot * ray -> cext [1]),
134                     exp(-ray -> rot * ray -> cext [2]));
135                    
136     /* Modify ray color and normalise */
137     multcolor(ray -> rcol, cvext);
138 <   colorNorm(ray -> rcol);
139 <  
202 <   /* Passed through medium */  
138 >   colorNorm(ray -> rcol);  
139 >
140     return 1;
141   }
142  
# Line 209 | Line 146 | void tracePhoton (RAY *ray)
146   /* Follow photon as it bounces around... */
147   {
148     long mod;
149 <   OBJREC* mat;
149 >   OBJREC *mat, *port = NULL;
150 >  
151 >   if (!ray -> parent) {
152 >      /* !!!  PHOTON PORT REJECTION SAMPLING HACK: get photon port for
153 >       * !!!  primary ray from ray -> ro, then reset the latter to NULL so
154 >       * !!!  as not to interfere with localhit() */
155 >      port = ray -> ro;
156 >      ray -> ro = NULL;
157 >   }
158  
159     if (ray -> rlvl > photonMaxBounce) {
160   #ifdef PMAP_RUNAWAY_WARN  
# Line 220 | Line 165 | void tracePhoton (RAY *ray)
165    
166     if (colorAvg(ray -> cext) > FTINY && !photonParticipate(ray))
167        return;
168 <      
168 >
169     if (localhit(ray, &thescene)) {
170        mod = ray -> ro -> omod;
171 <      
171 >
172 >      if (port && ray -> ro != port) {
173 >         /* !!! PHOTON PORT REJECTION SAMPLING HACK !!!
174 >          * Terminate photon if emitted from port without intersecting it;
175 >          * this can happen when the port's partitions extend beyond its
176 >          * actual geometry, e.g.  with polygons.  Since the total flux
177 >          * relayed by the port is based on the (in this case) larger
178 >          * partition area, it is overestimated; terminating these photons
179 >          * constitutes rejection sampling and thereby compensates any bias
180 >          * incurred by the overestimated flux.  */
181 > #ifdef PMAP_PORTREJECT_WARN
182 >         sprintf(errmsg, "photon outside port %s", ray -> ro -> oname);
183 >         error(WARNING, errmsg);
184 > #endif        
185 >         return;
186 >      }
187 >
188        if ((ray -> clipset && inset(ray -> clipset, mod)) || mod == OVOID) {
189           /* Transfer ray if modifier is VOID or clipped within antimatta */
190           RAY tray;
# Line 255 | Line 216 | static void preComputeGlobal (PhotonMap *pmap)
216  
217     repComplete = numPreComp = finalGather * pmap -> numPhotons;
218    
219 <   if (photonRepTime) {
220 <      sprintf(errmsg, "Precomputing irradiance for %ld global photons...\n",
219 >   if (verbose) {
220 >      sprintf(errmsg,
221 >              "\nPrecomputing irradiance for %ld global photons\n",
222                numPreComp);
223        eputs(errmsg);
224 + #if NIX      
225        fflush(stderr);
226 + #endif      
227     }
228    
229     /* Copy photon map for precomputed photons */
# Line 286 | Line 250 | static void preComputeGlobal (PhotonMap *pmap)
250    
251     for (i = 0; i < numPreComp; i++) {
252        /* Get random photon from stratified distribution in source heap to
253 <       * avoid duplicates and clutering */
253 >       * avoid duplicates and clustering */
254        pIdx = firstPhoton(pmap) +
255               (unsigned long)((i + pmapRandom(pmap -> randState)) /
256                               finalGather);
# Line 324 | Line 288 | static void preComputeGlobal (PhotonMap *pmap)
288     deletePhotons(pmap);
289     memcpy(pmap, &nuPmap, sizeof(PhotonMap));
290    
291 <   if (photonRepTime) {
292 <      eputs("Rebuilding precomputed photon map...\n");
291 >   if (verbose) {
292 >      eputs("\nRebuilding precomputed photon map\n");
293 > #if NIX      
294        fflush(stderr);
295 + #endif      
296     }
297  
298     /* Rebuild underlying data structure, destroying heap */  
# Line 345 | Line 311 | typedef struct {
311   void distribPhotons (PhotonMap **pmaps, unsigned numProc)
312   {
313     EmissionMap    emap;
314 <   char           errmsg2 [128], shmFname [255];
314 >   char           errmsg2 [128], shmFname [PMAP_TMPFNLEN];
315     unsigned       t, srcIdx, proc;
316     double         totalFlux = 0;
317     int            shmFile, stat, pid;
# Line 377 | Line 343 | void distribPhotons (PhotonMap **pmaps, unsigned numPr
343           initPhotonHeap(pmaps [t]);
344           /* Per-subprocess target count */
345           pmaps [t] -> distribTarget /= numProc;
346 +        
347 +         if (!pmaps [t] -> distribTarget)
348 +            error(INTERNAL, "no photons to distribute in distribPhotons");
349        }
350  
351     initPhotonEmissionFuncs();
352     initPhotonScatterFuncs();
353    
354 <   /* Get photon ports if specified */
355 <   if (ambincl == 1)
387 <      getPhotonPorts();
354 >   /* Get photon ports from modifier list */
355 >   getPhotonPorts(photonPortList);
356  
357     /* Get photon sensor modifiers */
358     getPhotonSensors(photonSensorList);
359    
360 + #if NIX
361     /* Set up shared mem for photon counters (zeroed by ftruncate) */
362 < #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);
362 >   strcpy(shmFname, PMAP_TMPFNAME);
363     shmFile = mkstemp(shmFname);
399 #endif      
364  
365 <   if (shmFile < 0)
366 <      error(SYSTEM, "failed opening shared memory file in distribPhotons");
365 >   if (shmFile < 0 || ftruncate(shmFile, sizeof(*photonCnt)) < 0)
366 >      error(SYSTEM, "failed shared mem init in distribPhotons");
367  
404   if (ftruncate(shmFile, sizeof(*photonCnt)) < 0)
405      error(SYSTEM, "failed setting shared memory size in distribPhotons");
406
368     photonCnt = mmap(NULL, sizeof(*photonCnt), PROT_READ | PROT_WRITE,
369                      MAP_SHARED, shmFile, 0);
370                      
371     if (photonCnt == MAP_FAILED)
372 <      error(SYSTEM, "failed mapping shared memory in distribPhotons");
372 >      error(SYSTEM, "failed mapping shared memory in distribPhotons");
373 > #else
374 >   /* Allocate photon counters statically on Windoze */
375 >   if (!(photonCnt = malloc(sizeof(PhotonCnt))))
376 >      error(SYSTEM, "failed trivial malloc in distribPhotons");
377 >   photonCnt -> numEmitted = photonCnt -> numComplete = 0;      
378 > #endif /* NIX */
379  
380 <   if (photonRepTime)
381 <      eputs("\n");
380 >   if (verbose) {
381 >      sprintf(errmsg, "\nIntegrating flux from %d sources", nsources);
382 >      
383 >      if (photonPorts) {
384 >         sprintf(errmsg2, " via %d ports", numPhotonPorts);
385 >         strcat(errmsg, errmsg2);
386 >      }
387 >      
388 >      strcat(errmsg, "\n");
389 >      eputs(errmsg);
390 >   }  
391    
392     /* ===================================================================
393      * FLUX INTEGRATION - Get total photon flux from light sources
# Line 425 | Line 401 | void distribPhotons (PhotonMap **pmaps, unsigned numPr
401                                                     : NULL;
402           photonPartition [emap.src -> so -> otype] (&emap);
403          
404 <         if (photonRepTime) {
405 <            sprintf(errmsg, "Integrating flux from source %s ",
404 >         if (verbose) {
405 >            sprintf(errmsg, "\tIntegrating flux from source %s ",
406                      source [srcIdx].so -> oname);
407 <                    
407 >
408              if (emap.port) {
409                 sprintf(errmsg2, "via port %s ",
410                         photonPorts [portCnt].so -> oname);
411                 strcat(errmsg, errmsg2);
412              }
413 <            
414 <            sprintf(errmsg2, "(%lu partitions)...\n", emap.numPartitions);
413 >
414 >            sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions);
415              strcat(errmsg, errmsg2);
416              eputs(errmsg);
417 + #if NIX            
418              fflush(stderr);
419 + #endif            
420           }
421          
422           for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
# Line 453 | Line 431 | void distribPhotons (PhotonMap **pmaps, unsigned numPr
431  
432     if (totalFlux < FTINY)
433        error(USER, "zero flux from light sources");
434 +      
435 +   /* Record start time for progress reports */
436 +   repStartTime = time(NULL);
437  
438 +   if (verbose) {
439 +      sprintf(errmsg, "\nPhoton distribution @ %d procs\n", numProc);
440 +      eputs(errmsg);
441 +   }
442 +
443     /* MAIN LOOP */  
444     for (proc = 0; proc < numProc; proc++) {
445 + #if NIX          
446        if (!(pid = fork())) {
447 <         /* SUBPROCESS ENTERS HERE.
448 <            All opened and memory mapped files are inherited */
447 >         /* SUBPROCESS ENTERS HERE; open and mmapped files inherited */
448 > #else
449 >      if (1) {
450 >         /* No subprocess under Windoze */
451 > #endif
452 >         /* Local photon counters for this subprocess */
453           unsigned       passCnt = 0, prePassCnt = 0;
454           unsigned long  lastNumPhotons [NUM_PMAP_TYPES];
455           unsigned long  localNumEmitted = 0; /* Num photons emitted by this
# Line 466 | Line 457 | void distribPhotons (PhotonMap **pmaps, unsigned numPr
457          
458           /* Seed RNGs from PID for decorellated photon distribution */
459           pmapSeed(randSeed + proc, partState);
460 <         pmapSeed(randSeed + proc, emitState);
461 <         pmapSeed(randSeed + proc, cntState);
462 <         pmapSeed(randSeed + proc, mediumState);
463 <         pmapSeed(randSeed + proc, scatterState);
464 <         pmapSeed(randSeed + proc, rouletteState);
465 <                  
460 >         pmapSeed(randSeed + (proc + 1) % numProc, emitState);
461 >         pmapSeed(randSeed + (proc + 2) % numProc, cntState);
462 >         pmapSeed(randSeed + (proc + 3) % numProc, mediumState);
463 >         pmapSeed(randSeed + (proc + 4) % numProc, scatterState);
464 >         pmapSeed(randSeed + (proc + 5) % numProc, rouletteState);
465 >              
466 > #ifdef DEBUG_PMAP          
467 >         /* Output child process PID after random delay to prevent corrupted
468 >          * console output due to race condition */
469 >         usleep(1e6 * pmapRandom(rouletteState));
470 >         fprintf(stderr, "Proc %d: PID = %d "
471 >                 "(waiting 10 sec to attach debugger...)\n",
472 >                 proc, getpid());
473 >         /* Allow time for debugger to attach to child process */
474 >         sleep(10);
475 > #endif            
476 >
477           for (t = 0; t < NUM_PMAP_TYPES; t++)
478              lastNumPhotons [t] = 0;
479              
# Line 491 | Line 493 | void distribPhotons (PhotonMap **pmaps, unsigned numPr
493                  * iterations; make it clear to user which photon maps are
494                  * missing so (s)he can check geometry and materials */
495                 if (++prePassCnt > maxPreDistrib) {
496 <                  sprintf(errmsg,
495 <                          "proc %d, source %s: too many prepasses",
496 <                          proc, source [srcIdx].so -> oname);              
496 >                  sprintf(errmsg, "proc %d: too many prepasses", proc);
497  
498                    for (t = 0; t < NUM_PMAP_TYPES; t++)
499                       if (pmaps [t] && !pmaps [t] -> numPhotons) {
# Line 565 | Line 565 | void distribPhotons (PhotonMap **pmaps, unsigned numPr
565                    break;
566              }
567  
568 <            /* Update shared completion counter for prog.report by parent */
568 >            /* Update shared completion counter for progreport by parent */
569              photonCnt -> numComplete += numEmit;                            
570  
571              /* PHOTON DISTRIBUTION LOOP */
# Line 578 | Line 578 | void distribPhotons (PhotonMap **pmaps, unsigned numPr
578                                ? photonPorts + portCnt : NULL;
579                    photonPartition [emap.src -> so -> otype] (&emap);
580  
581 <                  if (photonRepTime && !proc) {
581 >                  if (verbose && !proc) {
582 >                     /* Output from subproc 0 only to avoid race condition
583 >                      * on console I/O */
584                       if (!passCnt)
585 <                        sprintf(errmsg, "PREPASS %d on source %s ",
585 >                        sprintf(errmsg, "\tPREPASS %d on source %s ",
586                                  prePassCnt, source [srcIdx].so -> oname);
587                       else
588 <                        sprintf(errmsg, "MAIN PASS on source %s ",
588 >                        sprintf(errmsg, "\tMAIN PASS on source %s ",
589                                  source [srcIdx].so -> oname);
590 <                            
590 >
591                       if (emap.port) {
592                          sprintf(errmsg2, "via port %s ",
593                                  photonPorts [portCnt].so -> oname);
594                          strcat(errmsg, errmsg2);
595                       }
596 <                    
596 >
597                       sprintf(errmsg2, "(%lu partitions)\n",
598                               emap.numPartitions);
599                       strcat(errmsg, errmsg2);
600                       eputs(errmsg);
601 + #if NIX                    
602                       fflush(stderr);
603 + #endif                    
604                    }
605                    
606                    for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
# Line 611 | Line 615 | void distribPhotons (PhotonMap **pmaps, unsigned numPr
615                      
616                       /* Number of photons to emit from ziss partishunn --
617                        * proportional to flux; photon ray weight and scalar
618 <                      * flux are uniform (the latter only varying in RGB).
615 <                      * */
618 >                      * flux are uniform (latter only varying in RGB). */
619                       partNumEmit = numEmit * colorAvg(emap.partFlux) /
620                                     totalFlux;
621                       partEmitCnt = (unsigned long)partNumEmit;
# Line 633 | Line 636 | void distribPhotons (PhotonMap **pmaps, unsigned numPr
636                          /* Emit photon based on PDF and trace through scene
637                           * until absorbed/leaked */
638                          emitPhoton(&emap, &photonRay);
639 + #if 1
640 +                        if (emap.port)
641 +                           /* !!!  PHOTON PORT REJECTION SAMPLING HACK: set
642 +                            * !!!  photon port as fake hit object for
643 +                            * !!!  primary ray to check for intersection in
644 +                            * !!!  tracePhoton() */
645 +                           photonRay.ro = emap.port -> so;
646 + #endif
647                          tracePhoton(&photonRay);
648                       }                                          
649 <                    
649 >
650                       /* Update shared global photon count for each pmap */
651                       for (t = 0; t < NUM_PMAP_TYPES; t++)
652                          if (pmaps [t]) {
# Line 643 | Line 654 | void distribPhotons (PhotonMap **pmaps, unsigned numPr
654                                pmaps [t] -> numPhotons - lastNumPhotons [t];
655                             lastNumPhotons [t] = pmaps [t] -> numPhotons;
656                          }
657 + #if !NIX
658 +                     /* Synchronous progress report on Windoze */
659 +                     if (!proc && photonRepTime > 0 &&
660 +                           time(NULL) >= repLastTime + photonRepTime) {
661 +                        repEmitted = repProgress = photonCnt -> numEmitted;
662 +                        repComplete = photonCnt -> numComplete;                          
663 +                        pmapDistribReport();
664 +                     }
665 + #endif
666                    }
667                    
668                    portCnt++;
# Line 658 | Line 678 | void distribPhotons (PhotonMap **pmaps, unsigned numPr
678                    break;
679                 }
680              
681 <            if (t >= NUM_PMAP_TYPES) {
681 >            if (t >= NUM_PMAP_TYPES)
682                 /* No empty photon maps found; now do pass 2 */
683                 passCnt++;
664 #if 0
665               if (photonRepTime)
666                  eputs("\n");
667 #endif
668            }
684           } while (passCnt < 2);
685          
686 <         /* Unmap shared photon counters */
687 < #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! */
686 >         /* Flush heap buffa for every pmap one final time;
687 >          * avoids potential data corruption! */
688           for (t = 0; t < NUM_PMAP_TYPES; t++)
689              if (pmaps [t]) {
681 #if 0            
682               eputs("Final flush\n");
683 #endif              
690                 flushPhotonHeap(pmaps [t]);
691 <               fclose(pmaps [t] -> heap);
691 >               /* Heap file closed automatically on exit
692 >                  fclose(pmaps [t] -> heap); */
693   #ifdef DEBUG_PMAP              
694 <               sprintf(errmsg, "Proc %d: total %ld photons\n", getpid(),
694 >               sprintf(errmsg, "Proc %d: total %ld photons\n", proc,
695                         pmaps [t] -> numPhotons);
696                 eputs(errmsg);
697   #endif              
698              }
699 <
699 > #if NIX
700 >         /* Terminate subprocess */
701           exit(0);
702 + #endif
703        }
704        else if (pid < 0)
705           error(SYSTEM, "failed to fork subprocess in distribPhotons");        
706     }
707  
708 + #if NIX
709     /* PARENT PROCESS CONTINUES HERE */
700   /* Record start time and enable progress report signal handler */
701   repStartTime = time(NULL);
710   #ifdef SIGCONT
711 +   /* Enable progress report signal handler */
712     signal(SIGCONT, pmapDistribReport);
713 < #endif
714 <  
706 <   if (photonRepTime)
707 <      eputs("\n");
708 <  
709 <   /* Wait for subprocesses to complete while reporting progress */
713 > #endif  
714 >   /* Wait for subprocesses complete while reporting progress */
715     proc = numProc;
716     while (proc) {
717        while (waitpid(-1, &stat, WNOHANG) > 0) {
# Line 719 | Line 724 | void distribPhotons (PhotonMap **pmaps, unsigned numPr
724        
725        /* Nod off for a bit and update progress  */
726        sleep(1);
727 <      /* Update progress report from shared subprocess counters */
727 >
728 >      /* Asynchronous progress report from shared subprocess counters */  
729        repEmitted = repProgress = photonCnt -> numEmitted;
730 <      repComplete = photonCnt -> numComplete;
730 >      repComplete = photonCnt -> numComplete;      
731  
732 +      repProgress = repComplete = 0;
733        for (t = 0; t < NUM_PMAP_TYPES; t++)
734           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            
735              /* Get global photon count from shmem updated by subprocs */
736 <            pm -> numPhotons = photonCnt -> numPhotons [t];
737 < #endif            
736 >            repProgress += pm -> numPhotons = photonCnt -> numPhotons [t];
737 >            repComplete += pm -> distribTarget;
738           }
739 +      repComplete *= numProc;
740  
741        if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
742           pmapDistribReport();
# Line 741 | Line 744 | void distribPhotons (PhotonMap **pmaps, unsigned numPr
744        else signal(SIGCONT, pmapDistribReport);
745   #endif
746     }
747 + #endif /* NIX */
748  
749     /* ===================================================================
750      * POST-DISTRIBUTION - Set photon flux and build data struct for photon
751      * storage, etc.
752      * =================================================================== */
753   #ifdef SIGCONT    
754 +   /* Reset signal handler */
755     signal(SIGCONT, SIG_DFL);
756   #endif
757     free(emap.samples);
758    
759 <   /* Set photon flux (repProgress is total num emitted) */
759 >   /* Set photon flux */
760     totalFlux /= photonCnt -> numEmitted;
761 <  
761 > #if NIX  
762     /* Photon counters no longer needed, unmap shared memory */
763     munmap(photonCnt, sizeof(*photonCnt));
764     close(shmFile);
760 #if 0  
761   shm_unlink(shmFname);
762 #else
765     unlink(shmFname);
766 + #else
767 +   free(photonCnt);  
768   #endif      
769 <  
769 >   if (verbose)
770 >      eputs("\n");
771 >      
772     for (t = 0; t < NUM_PMAP_TYPES; t++)
773        if (pmaps [t]) {
774 <         if (photonRepTime) {
775 <            sprintf(errmsg, "\nBuilding %s photon map...\n", pmapName [t]);
774 >         if (verbose) {
775 >            sprintf(errmsg, "Building %s photon map\n", pmapName [t]);
776              eputs(errmsg);
777 + #if NIX            
778              fflush(stderr);
779 + #endif            
780           }
781          
782           /* Build underlying data structure; heap is destroyed */
783           buildPhotonMap(pmaps [t], &totalFlux, NULL, numProc);
784        }
785 <
785 >      
786     /* Precompute photon irradiance if necessary */
787 <   if (preCompPmap)
787 >   if (preCompPmap) {
788 >      if (verbose)
789 >         eputs("\n");
790        preComputeGlobal(preCompPmap);
791 < }
782 <
783 <
784 <
785 < void photonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
786 < /* Photon density estimate. Returns irradiance at ray -> rop. */
787 < {
788 <   unsigned                      i;
789 <   float                         r;
790 <   COLOR                         flux;
791 <   Photon                        *photon;
792 <   const PhotonSearchQueueNode   *sqn;
793 <
794 <   setcolor(irrad, 0, 0, 0);
795 <
796 <   if (!pmap -> maxGather)
797 <      return;
798 <      
799 <   /* Ignore sources */
800 <   if (ray -> ro && islight(objptr(ray -> ro -> omod) -> otype))
801 <      return;
802 <        
803 <   findPhotons(pmap, ray);
791 >   }      
792    
793 <   /* Need at least 2 photons */
794 <   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 <
817 <   if (pmap -> minGather == pmap -> maxGather) {
818 <      /* No bias compensation. Just do a plain vanilla estimate */
819 <      sqn = pmap -> squeue.node + 1;
820 <      
821 <      /* Average radius between furthest two photons to improve accuracy */      
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 -> 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 based on photon dist */
831 <         scalecolor(flux, 2 * (1 - sqn -> dist2 / r));
832 < #endif  
833 <         addcolor(irrad, flux);
834 <      }
835 <      
836 <      /* Divide by search area PI * r^2, 1 / PI required as ambient
837 <         normalisation factor */        
838 <      scalecolor(irrad, 1 / (PI * PI * r));
839 <      
840 <      return;
841 <   }
842 <   else
843 <      /* Apply bias compensation to density estimate */
844 <      biasComp(pmap, irrad);
845 < }
846 <
847 <
848 <
849 < void photonPreCompDensity (PhotonMap *pmap, RAY *r, COLOR irrad)
850 < /* Returns precomputed photon density estimate at ray -> rop. */
851 < {
852 <   Photon p;
853 <  
854 <   setcolor(irrad, 0, 0, 0);
855 <
856 <   /* Ignore sources */
857 <   if (r -> ro && islight(objptr(r -> ro -> omod) -> otype))
858 <      return;
859 <      
860 <   find1Photon(preCompPmap, r, &p);
861 <   getPhotonFlux(&p, irrad);
862 < }
863 <
864 <
865 <
866 < void volumePhotonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
867 < /* Photon volume density estimate. Returns irradiance at ray -> rop. */
868 < {
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 <      
880 <   findPhotons(pmap, ray);
881 <  
882 <   /* Need at least 2 photons */
883 <   if (pmap -> squeue.tail < 2)
884 <      return;
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 <      sqn = pmap -> squeue.node + 1;
894 <      
895 <      /* Average radius between furthest two photons to improve accuracy */      
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 -> 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, photon -> norm) / 127;
908 <            ph = 1 + gecc2 - 2 * ray -> gecc * ph;
909 <            ph = (1 - gecc2) / (ph * sqrt(ph));
910 <         }
911 <        
912 <         getPhotonFlux(photon, flux);
913 <         scalecolor(flux, ph);
914 <         addcolor(irrad, flux);
915 <      }
916 <      
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)));
920 <      return;
921 <   }
922 < #if 0
923 <   else
924 <      /* Apply bias compensation to density estimate */
925 <      volumeBiasComp(pmap, ray, irrad);
926 < #endif      
793 >   if (verbose)
794 >      eputs("\n");
795   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines