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

Comparing ray/src/rt/pmapcontrib.c (file contents):
Revision 2.1 by greg, Tue Feb 24 19:39:26 2015 UTC vs.
Revision 2.18 by rschregle, Thu Jun 7 19:26:04 2018 UTC

# Line 1 | Line 1
1 + #ifndef lint
2 + static const char RCSid[] = "$Id$";
3 + #endif
4 +
5   /*
6 <   ==================================================================
7 <   Photon map support for light source contributions
6 >   ======================================================================
7 >   Photon map for light source contributions
8  
9     Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
10 <   (c) Fraunhofer Institute for Solar Energy Systems,
11 <       Lucerne University of Applied Sciences & Arts  
12 <   ==================================================================
10 >   (c) Lucerne University of Applied Sciences and Arts,
11 >       supported by the Swiss National Science Foundation (SNSF, #147053)
12 >   ======================================================================
13    
14     $Id$
15   */
16  
17  
18   #include "pmapcontrib.h"
15 #include "pmap.h"
19   #include "pmapmat.h"
20   #include "pmapsrc.h"
21   #include "pmaprand.h"
# Line 20 | Line 23
23   #include "pmapdiag.h"
24   #include "rcontrib.h"
25   #include "otypes.h"
26 < #include <signal.h>
26 > #if NIX
27 >   #include <sys/mman.h>
28 >   #include <sys/wait.h>  
29 > #endif
30  
31  
32 <
33 < extern int contrib;     /* coeff/contrib flag */
34 <
35 <
36 <
37 < static void setPmapContribParams (PhotonMap *pmap, LUTAB *srcContrib)
38 < /* Set parameters for light source contributions */
32 > static PhotonPrimaryIdx newPhotonPrimary (PhotonMap *pmap,
33 >                                          const RAY *primRay,
34 >                                          FILE *primHeap)
35 > /* Add primary ray for emitted photon and save light source index, origin on
36 > * source, and emitted direction; used by contrib photons. The current
37 > * primary is stored in pmap -> lastPrimary.  If the previous primary
38 > * contributed photons (has srcIdx >= 0), it's appended to primHeap.  If
39 > * primRay == NULL, the current primary is still flushed, but no new primary
40 > * is set.  Returns updated primary counter pmap -> numPrimary.  */
41   {
42 <   /* Set light source modifier list and appropriate callback to extract
43 <    * their contributions from the photon map */
44 <   if (pmap) {
45 <      pmap -> srcContrib = srcContrib;
46 <      pmap -> lookup = photonContrib;
47 <      /* Ensure we get all requested photon contribs during lookups */
48 <      pmap -> gatherTolerance = 1.0;
42 >   if (!pmap || !primHeap)
43 >      return 0;
44 >      
45 >   /* Check if last primary ray has spawned photons (srcIdx >= 0, see
46 >    * newPhoton()), in which case we save it to the primary heap file
47 >    * before clobbering it */
48 >   if (pmap -> lastPrimary.srcIdx >= 0) {
49 >      if (!fwrite(&pmap -> lastPrimary, sizeof(PhotonPrimary), 1, primHeap))
50 >         error(SYSTEM, "failed writing photon primary in newPhotonPrimary");
51 >        
52 >      pmap -> numPrimary++;
53 >      if (pmap -> numPrimary > PMAP_MAXPRIMARY)
54 >         error(INTERNAL, "photon primary overflow in newPhotonPrimary");
55     }
42 }
56  
57 +   /* Mark unused with negative source index until path spawns a photon (see
58 +    * newPhoton()) */
59 +   pmap -> lastPrimary.srcIdx = -1;
60 +    
61 +   if (primRay) {
62 +      FVECT dvec;
63  
64 <
65 < static void checkPmapContribs (const PhotonMap *pmap, LUTAB *srcContrib)
66 < /* Check modifiers for light source contributions */
67 < {
68 <   const PhotonPrimary *primary = pmap -> primary;
69 <   OBJREC *srcMod;
70 <   unsigned long i, found = 0;
71 <  
72 <   /* Make sure at least one of the modifiers is actually in the pmap,
73 <    * otherwise findPhotons() winds up in an infinite loop! */
55 <   for (i = pmap -> primarySize; i; --i, ++primary) {
56 <      if (primary -> srcIdx < 0 || primary -> srcIdx >= nsources)
57 <         error(INTERNAL, "invalid light source index in photon map");
58 <        
59 <      srcMod = objptr(source [primary -> srcIdx].so -> omod);
60 <      if ((MODCONT*)lu_find(srcContrib, srcMod -> oname) -> data)
61 <         ++found;
64 > #ifdef PMAP_PRIMARYDIR            
65 >      /* Reverse incident direction to point to light source */
66 >      dvec [0] = -primRay -> rdir [0];
67 >      dvec [1] = -primRay -> rdir [1];
68 >      dvec [2] = -primRay -> rdir [2];
69 >      pmap -> lastPrimary.dir = encodedir(dvec);
70 > #endif      
71 > #ifdef PMAP_PRIMARYPOS      
72 >      VCOPY(pmap -> lastPrimary.pos, primRay -> rop);
73 > #endif      
74     }
75    
76 <   if (!found)
65 <      error(USER, "modifiers not in photon map");
76 >   return pmap -> numPrimary;
77   }
67
68  
78  
79 < void initPmapContrib (LUTAB *srcContrib, unsigned numSrcContrib)
79 >
80 >
81 > #ifdef DEBUG_PMAP
82 > static int checkPrimaryHeap (FILE *file)
83 > /* Check heap for ordered primaries */
84   {
85 <   unsigned t;
85 >   Photon   p, lastp;
86 >   int      i, dup;
87    
88 <   for (t = 0; t < NUM_PMAP_TYPES; t++)
89 <      if (photonMaps [t] && t != PMAP_TYPE_CONTRIB) {
76 <         sprintf(errmsg, "%s photon map does not support contributions",
77 <                 pmapName [t]);
78 <         error(USER, errmsg);
79 <      }
88 >   rewind(file);
89 >   memset(&lastp, 0, sizeof(lastp));
90    
91 <   /* Get params */
92 <   setPmapContribParams(contribPmap, srcContrib);
93 <  
94 <   if (contribPhotonMapping) {
95 <      if (contribPmap -> maxGather < numSrcContrib) {
96 <         /* Adjust density estimate bandwidth if lower than modifier
97 <          * count, otherwise contributions are missing */
98 <         error(WARNING, "contrib density estimate bandwidth too low, "
99 <                        "adjusting to modifier count");
100 <         contribPmap -> maxGather = numSrcContrib;
91 >   while (fread(&p, sizeof(p), 1, file)) {
92 >      dup = 1;
93 >      
94 >      for (i = 0; i <= 2; i++) {
95 >         if (p.pos [i] < thescene.cuorg [i] ||
96 >             p.pos [i] > thescene.cuorg [i] + thescene.cusize) {
97 >            
98 >            sprintf(errmsg, "corrupt photon in heap at [%f, %f, %f]\n",
99 >                    p.pos [0], p.pos [1], p.pos [2]);
100 >            error(WARNING, errmsg);
101 >         }
102 >        
103 >         dup &= p.pos [i] == lastp.pos [i];
104        }
105        
106 <      /* Sanity check */
107 <      checkPmapContribs(contribPmap, srcContrib);
106 >      if (dup) {
107 >         sprintf(errmsg,
108 >                 "consecutive duplicate photon in heap at [%f, %f, %f]\n",
109 >                 p.pos [0], p.pos [1], p.pos [2]);
110 >         error(WARNING, errmsg);
111 >      }
112     }
113 +
114 +   return 0;
115   }
116 + #endif
117  
118  
119  
120 < void photonContrib (PhotonMap *pmap, RAY *ray, COLOR irrad)
121 < /* Sum up light source contributions from photons in pmap->srcContrib */
120 > static PhotonPrimaryIdx buildPrimaries (PhotonMap *pmap, FILE **primaryHeap,
121 >                                        char **primaryHeapFname,
122 >                                        PhotonPrimaryIdx *primaryOfs,
123 >                                        unsigned numHeaps)
124 > /* Consolidate per-subprocess photon primary heaps into the primary array
125 > * pmap -> primaries. Returns offset for primary index linearisation in
126 > * numPrimary. The heap files in primaryHeap are closed on return. */
127   {
128 <   unsigned       i;
129 <   PhotonSQNode   *sq;
105 <   float          r, invArea;
106 <   RREAL          rayCoeff [3];
107 <   FVECT          rdir, rop;
108 <
109 <   setcolor(irrad, 0, 0, 0);
110 <
111 <   if (!pmap -> maxGather)
112 <      return;
113 <      
114 <   /* Ignore sources */
115 <   if (ray -> ro)
116 <      if (islight(objptr(ray -> ro -> omod) -> otype))
117 <         return;
118 <
119 <   /* Set context for binning function evaluation and get cumulative path
120 <    * coefficient up to photon lookup point */
121 <   worldfunc(RCCONTEXT, ray);
122 <   raycontrib(rayCoeff, ray, PRIMARY);
123 <
124 <   /* Save incident ray's direction and hitpoint */
125 <   VCOPY(rdir, ray -> rdir);
126 <   VCOPY(rop, ray -> rop);
127 <
128 <   /* Lookup photons */
129 <   pmap -> squeueEnd = 0;
130 <   findPhotons(pmap, ray);
128 >   PhotonPrimaryIdx  heapLen;
129 >   unsigned          heap;
130    
131 <   /* Need at least 2 photons */
132 <   if (pmap -> squeueEnd < 2) {
134 <      #ifdef PMAP_NONEFOUND
135 <         sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)",
136 <                 ray -> ro ? ray -> ro -> oname : "<null>",
137 <                 ray -> rop [0], ray -> rop [1], ray -> rop [2]);
138 <         error(WARNING, errmsg);
139 <      #endif
131 >   if (!pmap || !primaryHeap || !primaryOfs || !numHeaps)
132 >      return 0;
133        
134 <      return;
142 <   }
143 <
144 <   /* Average (squared) radius between furthest two photons to improve
145 <    * accuracy and get inverse search area 1 / (PI * r^2), with extra
146 <    * normalisation factor 1 / PI for ambient calculation */
147 <   sq = pmap -> squeue + 1;
148 <   r = max(sq -> dist, (sq + 1) -> dist);
149 <   r = 0.25 * (pmap -> maxDist + r + 2 * sqrt(pmap -> maxDist * r));  
150 <   invArea = 1 / (PI * PI * r);
134 >   pmap -> numPrimary = 0;
135    
136 <   /* Skip the extra photon */
137 <   for (i = 1 ; i < pmap -> squeueEnd; i++, sq++) {        
154 <      COLOR flux;
136 >   for (heap = 0; heap < numHeaps; heap++) {
137 >      primaryOfs [heap] = pmap -> numPrimary;
138        
139 <      /* Get photon's contribution to density estimate */
140 <      getPhotonFlux(sq -> photon, flux);
141 <      scalecolor(flux, invArea);
142 < #ifdef PMAP_EPANECHNIKOV
160 <      /* Apply Epanechnikov kernel to photon flux (dists are squared) */
161 <      scalecolor(flux, 2 * (1 - sq -> dist / r));
162 < #endif
163 <      addcolor(irrad, flux);
164 <      
165 <      if (pmap -> srcContrib) {
166 <         const PhotonPrimary *primary = pmap -> primary +
167 <                                        sq -> photon -> primary;
168 <         OBJREC *srcMod = objptr(source [primary -> srcIdx].so -> omod);
169 <         MODCONT *srcContrib = (MODCONT*)lu_find(pmap -> srcContrib,
170 <                                                 srcMod -> oname) -> data;
171 <        
172 <         if (srcContrib) {
173 <            /* Photon's emitting light source has modifier whose
174 <             * contributions are sought */
175 <            int srcBin;
139 >      if (fseek(primaryHeap [heap], 0, SEEK_END) < 0)
140 >         error(SYSTEM, "failed photon primary seek in buildPrimaries");
141 >      pmap -> numPrimary += heapLen = ftell(primaryHeap [heap]) /
142 >                                      sizeof(PhotonPrimary);      
143  
144 <            /* Set incident dir and origin of photon's primary ray on
145 <             * light source for dummy shadow ray, and evaluate binning
146 <             * function */
147 <            VCOPY(ray -> rdir, primary -> dir);
148 <            VCOPY(ray -> rop, primary -> org);
182 <            srcBin = evalue(srcContrib -> binv) + .5;
144 >      pmap -> primaries = realloc(pmap -> primaries,
145 >                                  pmap -> numPrimary *
146 >                                  sizeof(PhotonPrimary));  
147 >      if (!pmap -> primaries)
148 >         error(SYSTEM, "failed photon primary alloc in buildPrimaries");
149  
150 <            if (srcBin < 0 || srcBin >= srcContrib -> nbins) {
151 <               error(WARNING, "bad bin number (ignored)");
152 <               continue;
153 <            }
154 <            
155 <            if (!contrib) {
156 <               /* Ray coefficient mode; normalise by light source radiance
191 <                * after applying distrib pattern */
192 <               int j;
193 <               raytexture(ray, srcMod -> omod);
194 <               setcolor(ray -> rcol, srcMod -> oargs.farg [0],
195 <                        srcMod -> oargs.farg [1], srcMod -> oargs.farg [2]);
196 <               multcolor(ray -> rcol, ray -> pcol);
197 <               for (j = 0; j < 3; j++)
198 <                  flux [j] = ray -> rcol [j] ? flux [j] / ray -> rcol [j]
199 <                                             : 0;
200 <            }
201 <                    
202 <            multcolor(flux, rayCoeff);
203 <            addcolor(srcContrib -> cbin [srcBin], flux);
204 <         }
205 <         else fprintf(stderr, "Skipped contrib from %s\n", srcMod -> oname);
206 <      }
150 >      rewind(primaryHeap [heap]);
151 >      if (fread(pmap -> primaries + primaryOfs [heap], sizeof(PhotonPrimary),
152 >                heapLen, primaryHeap [heap]) != heapLen)
153 >         error(SYSTEM, "failed reading photon primaries in buildPrimaries");
154 >      
155 >      fclose(primaryHeap [heap]);
156 >      unlink(primaryHeapFname [heap]);
157     }
158    
159 <   /* Restore incident ray's direction and hitpoint */
160 <   VCOPY(ray -> rdir, rdir);
211 <   VCOPY(ray -> rop, rop);
212 <    
213 <   return;
214 < }
159 >   return pmap -> numPrimary;
160 > }      
161  
162  
163  
164 < void distribPhotonContrib (PhotonMap* pm)
165 < {
166 <   EmissionMap emap;
221 <   char errmsg2 [128];
222 <   unsigned srcIdx;
223 <   double *srcFlux;                 /* Emitted flux per light source */
224 <   const double srcDistribTarget =  /* Target photon count per source */
225 <      nsources ? (double)pm -> distribTarget / nsources : 0;  
164 > /* Defs for photon emission counter array passed by sub-processes to parent
165 > * via shared memory */
166 > typedef  unsigned long  PhotonContribCnt;
167  
168 + /* Indices for photon emission counter array: num photons stored and num
169 + * emitted per source */
170 + #define  PHOTONCNT_NUMPHOT    0
171 + #define  PHOTONCNT_NUMEMIT(n) (1 + n)
172 +
173 +
174 +
175 +
176 +
177 +
178 + void distribPhotonContrib (PhotonMap* pm, unsigned numProc)
179 + {
180 +   EmissionMap       emap;
181 +   char              errmsg2 [128], shmFname [PMAP_TMPFNLEN];
182 +   unsigned          srcIdx, proc;
183 +   int               shmFile, stat, pid;
184 +   double            *srcFlux,         /* Emitted flux per light source */
185 +                     srcDistribTarget; /* Target photon count per source */
186 +   PhotonContribCnt  *photonCnt;       /* Photon emission counter array */
187 +   unsigned          photonCntSize = sizeof(PhotonContribCnt) *
188 +                                     PHOTONCNT_NUMEMIT(nsources);
189 +   FILE              **primaryHeap = NULL;
190 +   char              **primaryHeapFname = NULL;
191 +   PhotonPrimaryIdx  *primaryOfs = NULL;
192 +                                    
193     if (!pm)
194 <      error(USER, "no photon map defined");
194 >      error(USER, "no photon map defined in distribPhotonContrib");
195        
196     if (!nsources)
197 <      error(USER, "no light sources");
198 <  
197 >      error(USER, "no light sources in distribPhotonContrib");
198 >
199 >   if (nsources > MAXMODLIST)
200 >      error(USER, "too many light sources in distribPhotonContrib");
201 >      
202     /* Allocate photon flux per light source; this differs for every
203      * source as all sources contribute the same number of distributed
204      * photons (srcDistribTarget), hence the number of photons emitted per
205      * source does not correlate with its emitted flux. The resulting flux
206      * per photon is therefore adjusted individually for each source. */
207     if (!(srcFlux = calloc(nsources, sizeof(double))))
208 <      error(SYSTEM, "cannot allocate source flux");
208 >      error(SYSTEM, "can't allocate source flux in distribPhotonContrib");
209  
210 <   /* ================================================================
211 <    * INITIALISASHUNN - Set up emisshunn and scattering funcs
212 <    * ================================================================ */
210 >   /* ===================================================================
211 >    * INITIALISATION - Set up emission and scattering funcs
212 >    * =================================================================== */
213     emap.samples = NULL;
214     emap.src = NULL;
215     emap.maxPartitions = MAXSPART;
216     emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1);
217     if (!emap.partitions)
218 <      error(USER, "can't allocate source partitions");
218 >      error(USER, "can't allocate source partitions in distribPhotonContrib");
219  
220 +   /* Initialise contrib photon map */  
221     initPhotonMap(pm, PMAP_TYPE_CONTRIB);
222 +   initPhotonHeap(pm);
223     initPhotonEmissionFuncs();
224     initPhotonScatterFuncs();
225    
226 <   /* Get photon ports if specified */
227 <   if (ambincl == 1)
228 <      getPhotonPorts();
226 >   /* Per-subprocess / per-source target counts */
227 >   pm -> distribTarget /= numProc;
228 >   srcDistribTarget = nsources ? (double)pm -> distribTarget / nsources : 0;  
229 >  
230 >   if (!pm -> distribTarget)
231 >      error(INTERNAL, "no photons to distribute in distribPhotonContrib");
232 >  
233 >   /* Get photon ports from modifier list */
234 >   getPhotonPorts(photonPortList);
235        
236     /* Get photon sensor modifiers */
237     getPhotonSensors(photonSensorList);      
238 +
239 + #if NIX  
240 +   /* Set up shared mem for photon counters (zeroed by ftruncate) */
241 +   strcpy(shmFname, PMAP_TMPFNAME);
242 +   shmFile = mkstemp(shmFname);
243    
244 <   /* Seed RNGs for photon distribution */
245 <   pmapSeed(randSeed, partState);
264 <   pmapSeed(randSeed, emitState);
265 <   pmapSeed(randSeed, cntState);
266 <   pmapSeed(randSeed, mediumState);
267 <   pmapSeed(randSeed, scatterState);
268 <   pmapSeed(randSeed, rouletteState);
244 >   if (shmFile < 0 || ftruncate(shmFile, photonCntSize) < 0)
245 >      error(SYSTEM, "failed shared mem init in distribPhotonContrib");
246  
247 <   /* Record start time and enable progress report signal handler */
248 <   repStartTime = time(NULL);
249 <   signal(SIGCONT, pmapDistribReport);
247 >   photonCnt = mmap(NULL, photonCntSize, PROT_READ | PROT_WRITE,
248 >                    MAP_SHARED, shmFile, 0);
249 >                    
250 >   if (photonCnt == MAP_FAILED)
251 >      error(SYSTEM, "failed shared mem mapping in distribPhotonContrib");
252 > #else
253 >   /* Allocate photon counters statically on Windoze */
254 >   if (!(photonCnt = malloc(photonCntSize)))
255 >      error(SYSTEM, "failed trivial malloc in distribPhotonContrib");
256    
257 <   for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
258 <      unsigned portCnt = 0, passCnt = 0, prePassCnt = 0;
259 <      double srcNumEmit = 0;          /* # photons to emit from source */
277 <      unsigned long srcNumDistrib = pm -> heapEnd; /* # photons stored */
257 >   for (srcIdx = 0; srcIdx < PHOTONCNT_NUMEMIT(nsources); srcIdx++)
258 >      photonCnt [srcIdx] = 0;
259 > #endif /* NIX */
260  
261 <      srcFlux [srcIdx] = repProgress = 0;
261 >   if (verbose) {
262 >      sprintf(errmsg, "\nIntegrating flux from %d sources", nsources);
263 >
264 >      if (photonPorts) {
265 >         sprintf(errmsg2, " via %d ports", numPhotonPorts);
266 >         strcat(errmsg, errmsg2);
267 >      }
268 >
269 >      strcat(errmsg, "\n");
270 >      eputs(errmsg);
271 >   }
272 >
273 >   /* =============================================================
274 >    * FLUX INTEGRATION - Get total flux emitted from sources/ports
275 >    * ============================================================= */  
276 >   for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
277 >      unsigned portCnt = 0;      
278 >      srcFlux [srcIdx] = 0;
279        emap.src = source + srcIdx;
280        
281 <      if (photonRepTime)
282 <         eputs("\n");
283 <      
285 <      /* =============================================================
286 <       * FLUX INTEGRATION - Get total flux emitted from light source
287 <       * ============================================================= */
288 <      do {
289 <         emap.port = emap.src -> sflags & SDISTANT
290 <                     ? photonPorts + portCnt : NULL;
281 >      do {  /* Need at least one iteration if no ports! */      
282 >         emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt
283 >                                                   : NULL;
284           photonPartition [emap.src -> so -> otype] (&emap);
285 <        
286 <         if (photonRepTime) {
287 <            sprintf(errmsg, "Integrating flux from source %s (mod %s) ",
288 <                    source [srcIdx].so -> oname,
289 <                    objptr(source [srcIdx].so -> omod) -> oname);
297 <                    
285 >
286 >         if (verbose) {
287 >            sprintf(errmsg, "\tIntegrating flux from source %s ",
288 >                    source [srcIdx].so -> oname);
289 >
290              if (emap.port) {
291                 sprintf(errmsg2, "via port %s ",
292                         photonPorts [portCnt].so -> oname);
293                 strcat(errmsg, errmsg2);
294              }
295 <            
296 <            sprintf(errmsg2, "(%lu partitions)...\n",
305 <                    emap.numPartitions);
295 >
296 >            sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions);
297              strcat(errmsg, errmsg2);
298              eputs(errmsg);
299 + #if NIX            
300              fflush(stderr);
301 <         }
301 > #endif            
302 >         }                    
303          
304 <         for (emap.partitionCnt = 0;
312 <              emap.partitionCnt < emap.numPartitions;
304 >         for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
305                emap.partitionCnt++) {
306              initPhotonEmission(&emap, pdfSamples);
307              srcFlux [srcIdx] += colorAvg(emap.partFlux);
308           }
309          
310           portCnt++;
311 <      } while (portCnt < numPhotonPorts);
312 <
311 >      } while (portCnt < numPhotonPorts);        
312 >      
313        if (srcFlux [srcIdx] < FTINY) {
314           sprintf(errmsg, "source %s has zero emission",
315                   source [srcIdx].so -> oname);
316           error(WARNING, errmsg);
317        }
318 <      else {
319 <         /* ==========================================================
320 <          * 2-PASS PHOTON DISTRIBUTION
321 <          * Pass 1 (pre):  emit fraction of target photon count
322 <          * Pass 2 (main): based on outcome of pass 1, estimate
323 <          *                remaining number of photons to emit to
324 <          *                approximate target count
325 <          * ========================================================== */
326 <         do {
327 <            if (!passCnt) {  
328 <               /* INIT PASS 1 */
329 <               if (++prePassCnt > maxPreDistrib) {
330 <                  /* Warn if no photons contributed after sufficient
331 <                   * iterations */
332 <                  sprintf(errmsg, "too many prepasses, no photons "
341 <                          "from source %s", source [srcIdx].so -> oname);
342 <                  error(WARNING, errmsg);
343 <                  break;
344 <               }
318 >   }  
319 >  
320 >   /* Allocate & init per-subprocess primary heap files */
321 >   primaryHeap = calloc(numProc, sizeof(FILE*));
322 >   primaryHeapFname = calloc(numProc, sizeof(char*));
323 >   primaryOfs = calloc(numProc, sizeof(PhotonPrimaryIdx));
324 >   if (!primaryHeap || !primaryHeapFname || !primaryOfs)
325 >      error(SYSTEM, "failed primary heap allocation in "
326 >            "distribPhotonContrib");
327 >      
328 >   for (proc = 0; proc < numProc; proc++) {
329 >      primaryHeapFname [proc] = malloc(PMAP_TMPFNLEN);
330 >      if (!primaryHeapFname [proc])
331 >         error(SYSTEM, "failed primary heap file allocation in "
332 >               "distribPhotonContrib");
333                
334 <               /* Num to emit is fraction of target count */
335 <               srcNumEmit = preDistrib * srcDistribTarget;
336 <            }
334 >      mktemp(strcpy(primaryHeapFname [proc], PMAP_TMPFNAME));
335 >      if (!(primaryHeap [proc] = fopen(primaryHeapFname [proc], "w+b")))
336 >         error(SYSTEM, "failed opening primary heap file in "
337 >               "distribPhotonContrib");
338 >   }              
339  
340 <            else {            
341 <               /* INIT PASS 2 */                          
352 <               /* Based on the outcome of the predistribution we can now
353 <                * figure out how many more photons we have to emit from
354 <                * the current source to meet the target count,
355 <                * srcDistribTarget. This value is clamped to 0 in case
356 <                * the target has already been exceeded in pass 1.
357 <                * srcNumEmit and srcNumDistrib is the number of photons
358 <                * emitted and distributed (stored) from the current
359 <                * source in pass 1, respectively. */
360 <               srcNumDistrib = pm -> heapEnd - srcNumDistrib;
361 <               srcNumEmit *= srcNumDistrib
362 <                             ? max(srcDistribTarget/srcNumDistrib, 1) - 1
363 <                             : 0;
340 >   /* Record start time for progress reports */
341 >   repStartTime = time(NULL);
342  
343 <               if (!srcNumEmit)
344 <                  /* No photons left to distribute in main pass */
345 <                  break;
346 <            }
369 <        
370 <            /* Set completion count for progress report */
371 <            repComplete = srcNumEmit + repProgress;
372 <            portCnt = 0;
373 <                    
374 <            do {
375 <               emap.port = emap.src -> sflags & SDISTANT
376 <                           ? photonPorts + portCnt : NULL;
377 <               photonPartition [emap.src -> so -> otype] (&emap);
343 >   if (verbose) {
344 >      sprintf(errmsg, "\nPhoton distribution @ %d procs\n", numProc);
345 >      eputs(errmsg);
346 >   }
347                
348 <               if (photonRepTime) {
349 <                  if (!passCnt)
350 <                     sprintf(errmsg, "PREPASS %d on source %s (mod %s) ",
351 <                             prePassCnt, source [srcIdx].so -> oname,
352 <                             objptr(source[srcIdx].so->omod) -> oname);
353 <                  else
354 <                     sprintf(errmsg, "MAIN PASS on source %s (mod %s) ",
355 <                             source [srcIdx].so -> oname,
356 <                             objptr(source[srcIdx].so->omod) -> oname);
357 <                          
358 <                  if (emap.port) {
359 <                     sprintf(errmsg2, "via port %s ",
360 <                             photonPorts [portCnt].so -> oname);
361 <                     strcat(errmsg, errmsg2);
348 >   /* MAIN LOOP */
349 >   for (proc = 0; proc < numProc; proc++) {
350 > #if NIX          
351 >      if (!(pid = fork())) {
352 >         /* SUBPROCESS ENTERS HERE; opened and mmapped files inherited */
353 > #else
354 >      if (1) {
355 >         /* No subprocess under Windoze */
356 > #endif  
357 >         /* Local photon counters for this subprocess */
358 >         unsigned long  lastNumPhotons = 0, localNumEmitted = 0;
359 >         double         photonFluxSum = 0;   /* Accum. photon flux */
360 >
361 >         /* Seed RNGs from PID for decorellated photon distribution */
362 >         pmapSeed(randSeed + proc, partState);
363 >         pmapSeed(randSeed + (proc + 1) % numProc, emitState);
364 >         pmapSeed(randSeed + (proc + 2) % numProc, cntState);
365 >         pmapSeed(randSeed + (proc + 3) % numProc, mediumState);
366 >         pmapSeed(randSeed + (proc + 4) % numProc, scatterState);
367 >         pmapSeed(randSeed + (proc + 5) % numProc, rouletteState);
368 >
369 > #ifdef PMAP_SIGUSR                      
370 >   double partNumEmit;
371 >   unsigned long partEmitCnt;
372 >   double srcPhotonFlux, avgPhotonFlux;
373 >   unsigned       portCnt, passCnt, prePassCnt;
374 >   float          srcPreDistrib;
375 >   double         srcNumEmit;     /* # to emit from source */
376 >   unsigned long  srcNumDistrib;  /* # stored */
377 >
378 >   void sigUsrDiags()
379 >   /* Loop diags via SIGUSR1 */
380 >   {
381 >      sprintf(errmsg,
382 >              "********************* Proc %d Diags *********************\n"
383 >              "srcIdx = %d (%s)\nportCnt = %d (%s)\npassCnt = %d\n"
384 >              "srcFlux = %f\nsrcPhotonFlux = %f\navgPhotonFlux = %f\n"
385 >              "partNumEmit = %f\npartEmitCnt = %lu\n\n",
386 >              proc, srcIdx, findmaterial(source [srcIdx].so) -> oname,
387 >              portCnt, photonPorts [portCnt].so -> oname,
388 >              passCnt, srcFlux [srcIdx], srcPhotonFlux, avgPhotonFlux,
389 >              partNumEmit, partEmitCnt);
390 >      eputs(errmsg);
391 >      fflush(stderr);
392 >   }
393 > #endif
394 >        
395 > #ifdef PMAP_SIGUSR
396 >         signal(SIGUSR1, sigUsrDiags);
397 > #endif        
398 >
399 > #ifdef DEBUG_PMAP          
400 >         /* Output child process PID after random delay to prevent corrupted
401 >          * console output due to race condition */
402 >         usleep(1e6 * pmapRandom(rouletteState));
403 >         fprintf(stderr, "Proc %d: PID = %d "
404 >                 "(waiting 10 sec to attach debugger...)\n",
405 >                 proc, getpid());
406 >         /* Allow time for debugger to attach to child process */
407 >         sleep(10);
408 > #endif
409 >
410 >         /* =============================================================
411 >          * 2-PASS PHOTON DISTRIBUTION
412 >          * Pass 1 (pre):  emit fraction of target photon count
413 >          * Pass 2 (main): based on outcome of pass 1, estimate remaining
414 >          *                number of photons to emit to approximate target
415 >          *                count
416 >          * ============================================================= */
417 >         for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
418 > #ifndef PMAP_SIGUSR        
419 >            unsigned       portCnt, passCnt = 0, prePassCnt = 0;
420 >            float          srcPreDistrib = preDistrib;
421 >            double         srcNumEmit = 0;       /* # to emit from source */
422 >            unsigned long  srcNumDistrib = pm -> numPhotons;  /* # stored */
423 > #else
424 >            passCnt = prePassCnt = 0;
425 >            srcPreDistrib = preDistrib;
426 >            srcNumEmit = 0;       /* # to emit from source */
427 >            srcNumDistrib = pm -> numPhotons;  /* # stored */
428 > #endif            
429 >
430 >            if (srcFlux [srcIdx] < FTINY)
431 >               continue;
432 >                        
433 >            while (passCnt < 2) {
434 >               if (!passCnt) {  
435 >                  /* INIT PASS 1 */
436 >                  if (++prePassCnt > maxPreDistrib) {
437 >                     /* Warn if no photons contributed after sufficient
438 >                      * iterations; only output from subprocess 0 to reduce
439 >                      * console clutter */
440 >                     if (!proc) {
441 >                        sprintf(errmsg,
442 >                                "source %s: too many prepasses, skipped",
443 >                                source [srcIdx].so -> oname);
444 >                        error(WARNING, errmsg);
445 >                     }
446 >
447 >                     break;
448                    }
449                    
450 <                  sprintf(errmsg2, "(%lu partitions)...\n",
451 <                          emap.numPartitions);
397 <                  strcat(errmsg, errmsg2);
398 <                  eputs(errmsg);
399 <                  fflush(stderr);
450 >                  /* Num to emit is fraction of target count */
451 >                  srcNumEmit = srcPreDistrib * srcDistribTarget;
452                 }
453 <              
454 <               for (emap.partitionCnt = 0;
455 <                    emap.partitionCnt < emap.numPartitions;
456 <                    emap.partitionCnt++) {
457 <                  double partNumEmit;
406 <                  unsigned long partEmitCnt;
453 >               else {
454 >                  /* INIT PASS 2 */
455 > #ifndef PMAP_SIGUSR
456 >                  double srcPhotonFlux, avgPhotonFlux;
457 > #endif
458                    
459 <                  /* Get photon origin within current source partishunn
460 <                   * and build emission map */
461 <                  photonOrigin [emap.src -> so -> otype] (&emap);
462 <                  initPhotonEmission(&emap, pdfSamples);
463 <                  
464 <                  /* Number of photons to emit from ziss partishunn;
465 <                   * scale according to its normalised contribushunn to
466 <                   * the emitted source flux */
467 <                  partNumEmit = srcNumEmit * colorAvg(emap.partFlux) /
468 <                                srcFlux [srcIdx];
469 <                  partEmitCnt = (unsigned long)partNumEmit;
470 <                                      
471 <                  /* Probabilistically account for fractional photons */
472 <                  if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
473 <                     partEmitCnt++;
459 >                  /* Based on the outcome of the predistribution we can now
460 >                   * figure out how many more photons we have to emit from
461 >                   * the current source to meet the target count,
462 >                   * srcDistribTarget. This value is clamped to 0 in case
463 >                   * the target has already been exceeded in pass 1.
464 >                   * srcNumEmit and srcNumDistrib is the number of photons
465 >                   * emitted and distributed (stored) from the current
466 >                   * source in pass 1, respectively. */
467 >                  srcNumDistrib = pm -> numPhotons - srcNumDistrib;
468 >                  srcNumEmit *= srcNumDistrib
469 >                                ? max(srcDistribTarget/srcNumDistrib, 1) - 1
470 >                                : 0;
471 >
472 >                  if (!srcNumEmit)
473 >                     /* No photons left to distribute in main pass */
474 >                     break;
475                      
476 <                  /* Integer counter avoids FP rounding errors */
477 <                  while (partEmitCnt--) {                  
426 <                     RAY photonRay;
476 >                  srcPhotonFlux = srcFlux [srcIdx] / srcNumEmit;
477 >                  avgPhotonFlux = photonFluxSum / (srcIdx + 1);
478                    
479 <                     /* Emit photon according to PDF (if any), allocate
480 <                      * associated primary ray, and trace through scene
481 <                      * until absorbed/leaked */
482 <                     emitPhoton(&emap, &photonRay);
483 <                     addPhotonPrimary(pm, &photonRay);
484 <                     tracePhoton(&photonRay);
479 >                  if (avgPhotonFlux > FTINY &&
480 >                      srcPhotonFlux / avgPhotonFlux < FTINY) {
481 >                     /* Skip source if its photon flux is grossly below the
482 >                      * running average, indicating negligible contributions
483 >                      * at the expense of excessive distribution time; only
484 >                      * output from subproc 0 to reduce console clutter */
485 >                     if (!proc) {
486 >                        sprintf(errmsg,
487 >                                "source %s: itsy bitsy photon flux, skipped",
488 >                                source [srcIdx].so -> oname);                    
489 >                        error(WARNING, errmsg);
490 >                     }
491 >
492 >                     srcNumEmit = 0;   /* Or just break??? */
493 >                  }
494 >                        
495 >                  /* Update sum of photon flux per light source */
496 >                  photonFluxSum += srcPhotonFlux;
497 >               }
498 >                              
499 >               portCnt = 0;
500 >               do {    /* Need at least one iteration if no ports! */
501 >                  emap.src = source + srcIdx;
502 >                  emap.port = emap.src -> sflags & SDISTANT
503 >                              ? photonPorts + portCnt : NULL;
504 >                  photonPartition [emap.src -> so -> otype] (&emap);
505 >
506 >                  if (verbose && !proc) {
507 >                     /* Output from subproc 0 only to avoid race condition
508 >                      * on console I/O */
509 >                     if (!passCnt)
510 >                        sprintf(errmsg, "\tPREPASS %d on source %s ",
511 >                                prePassCnt, source [srcIdx].so -> oname);
512 >                     else
513 >                        sprintf(errmsg, "\tMAIN PASS on source %s ",
514 >                                source [srcIdx].so -> oname);
515 >
516 >                     if (emap.port) {
517 >                        sprintf(errmsg2, "via port %s ",
518 >                                photonPorts [portCnt].so -> oname);
519 >                        strcat(errmsg, errmsg2);
520 >                     }
521 >
522 >                     sprintf(errmsg2, "(%lu partitions)\n",
523 >                             emap.numPartitions);
524 >                     strcat(errmsg, errmsg2);                    
525 >                     eputs(errmsg);
526 > #if NIX                    
527 >                     fflush(stderr);
528 > #endif                    
529 >                  }                
530 >                  
531 >                  for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
532 >                       emap.partitionCnt++) {
533 > #ifndef PMAP_SIGUSR                      
534 >                     double partNumEmit;
535 >                     unsigned long partEmitCnt;
536 > #endif
537                      
538 <                     /* Record progress */
539 <                     repProgress++;
538 >                     /* Get photon origin within current source partishunn
539 >                      * and build emission map */
540 >                     photonOrigin [emap.src -> so -> otype] (&emap);
541 >                     initPhotonEmission(&emap, pdfSamples);
542                      
543 <                     if (photonRepTime > 0 &&
544 <                         time(NULL) >= repLastTime + photonRepTime)
543 >                     /* Number of photons to emit from ziss partishunn;
544 >                      * scale according to its normalised contribushunn to
545 >                      * the emitted source flux */                    
546 >                     partNumEmit = srcNumEmit * colorAvg(emap.partFlux) /
547 >                                   srcFlux [srcIdx];                    
548 >                     partEmitCnt = (unsigned long)partNumEmit;
549 >                                                              
550 >                     /* Probabilistically account for fractional photons */
551 >                     if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
552 >                        partEmitCnt++;
553 >                        
554 >                     /* Update local and shared global emission counter */
555 >                     photonCnt [PHOTONCNT_NUMEMIT(srcIdx)] += partEmitCnt;
556 >                     localNumEmitted += partEmitCnt;                                    
557 >                    
558 >                     /* Integer counter avoids FP rounding errors during
559 >                      * iteration */
560 >                     while (partEmitCnt--) {
561 >                        RAY photonRay;
562 >                    
563 >                        /* Emit photon according to PDF (if any), allocate
564 >                         * associated primary ray, and trace through scene
565 >                         * until absorbed/leaked; emitPhoton() sets the
566 >                         * emitting light source index in photonRay */
567 >                        emitPhoton(&emap, &photonRay);
568 > #if 1
569 >                        if (emap.port)
570 >                           /* !!!  PHOTON PORT REJECTION SAMPLING HACK: set
571 >                            * !!!  photon port as fake hit object for
572 >                            * !!!  primary ray to check for intersection in
573 >                            * !!!  tracePhoton() */                        
574 >                           photonRay.ro = emap.port -> so;
575 > #endif
576 >                        newPhotonPrimary(pm, &photonRay, primaryHeap[proc]);
577 >                        /* Set subprocess index in photonRay for post-
578 >                         * distrib primary index linearisation; this is
579 >                         * propagated with the primary index in photonRay
580 >                         * and set for photon hits by newPhoton() */
581 >                        PMAP_SETRAYPROC(&photonRay, proc);
582 >                        tracePhoton(&photonRay);
583 >                     }
584 >                    
585 >                     /* Update shared global photon count */                    
586 >                     photonCnt [PHOTONCNT_NUMPHOT] += pm -> numPhotons -
587 >                                                      lastNumPhotons;
588 >                     lastNumPhotons = pm -> numPhotons;
589 > #if !NIX
590 >                     /* Synchronous progress report on Windoze */
591 >                     if (!proc && photonRepTime > 0 &&
592 >                           time(NULL) >= repLastTime + photonRepTime) {
593 >                        unsigned s;                        
594 >                        repComplete = pm -> distribTarget * numProc;
595 >                        repProgress = photonCnt [PHOTONCNT_NUMPHOT];
596 >                        
597 >                        for (repEmitted = 0, s = 0; s < nsources; s++)
598 >                           repEmitted += photonCnt [PHOTONCNT_NUMEMIT(s)];
599 >
600                          pmapDistribReport();
601 <                     #ifndef BSD
602 <                        else signal(SIGCONT, pmapDistribReport);
443 <                     #endif
601 >                     }
602 > #endif
603                    }
445               }
446                          
447               portCnt++;
448            } while (portCnt < numPhotonPorts);
604  
605 <            if (pm -> heapEnd == srcNumDistrib)
606 <               /* Double preDistrib in case no photons were stored
607 <                * for this source and redo pass 1 */
608 <               preDistrib *= 2;
609 <            else {
610 <               /* Now do pass 2 */
611 <               passCnt++;
612 <               if (photonRepTime)
613 <                  eputs("\n");
605 >                  portCnt++;
606 >               } while (portCnt < numPhotonPorts);                  
607 >
608 >               if (pm -> numPhotons == srcNumDistrib) {
609 >                  /* Double predistrib factor in case no photons were stored
610 >                   * for this source and redo pass 1 */
611 >                  srcPreDistrib *= 2;
612 >               }
613 >               else {
614 >                  /* Now do pass 2 */
615 >                  passCnt++;
616 >               }
617              }
618 <         } while (passCnt < 2);
619 <        
620 <         /* Flux per photon emitted from this source; repProgress is the
621 <          * number of emitted photons after both passes */
622 <         srcFlux [srcIdx] = repProgress ? srcFlux [srcIdx] / repProgress
623 <                                        : 0;
618 >         }
619 >                        
620 >         /* Flush heap buffa one final time to prevent data corruption */
621 >         flushPhotonHeap(pm);        
622 >         /* Flush final photon primary to primary heap file */
623 >         newPhotonPrimary(pm, NULL, primaryHeap [proc]);
624 >         /* Heap files closed automatically on exit
625 >            fclose(pm -> heap);
626 >            fclose(primaryHeap [proc]); */
627 >                  
628 > #ifdef DEBUG_PMAP
629 >         sprintf(errmsg, "Proc %d total %ld photons\n", proc,
630 >                 pm -> numPhotons);
631 >         eputs(errmsg);
632 >         fflush(stderr);
633 > #endif
634 >
635 > #ifdef PMAP_SIGUSR
636 >         signal(SIGUSR1, SIG_DFL);
637 > #endif
638 >
639 > #if NIX
640 >         /* Terminate subprocess */
641 >         exit(0);
642 > #endif
643        }
644 +      else if (pid < 0)
645 +         error(SYSTEM, "failed to fork subprocess in distribPhotonContrib");
646     }
647  
648 + #if NIX
649 +   /* PARENT PROCESS CONTINUES HERE */
650 + #ifdef SIGCONT
651 +   /* Enable progress report signal handler */
652 +   signal(SIGCONT, pmapDistribReport);
653 + #endif
654 +   /* Wait for subprocesses to complete while reporting progress */
655 +   proc = numProc;
656 +   while (proc) {
657 +      while (waitpid(-1, &stat, WNOHANG) > 0) {
658 +         /* Subprocess exited; check status */
659 +         if (!WIFEXITED(stat) || WEXITSTATUS(stat))
660 +            error(USER, "failed photon distribution");
661 +      
662 +         --proc;
663 +      }
664 +      
665 +      /* Nod off for a bit and update progress  */
666 +      sleep(1);
667 +
668 +      /* Asynchronous progress report from shared subprocess counters */      
669 +      repComplete = pm -> distribTarget * numProc;
670 +      repProgress = photonCnt [PHOTONCNT_NUMPHOT];      
671 +      
672 +      for (repEmitted = 0, srcIdx = 0; srcIdx < nsources; srcIdx++)
673 +         repEmitted += photonCnt [PHOTONCNT_NUMEMIT(srcIdx)];
674 +
675 +      /* Get global photon count from shmem updated by subprocs */
676 +      pm -> numPhotons = photonCnt [PHOTONCNT_NUMPHOT];
677 +
678 +      if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
679 +         pmapDistribReport();
680 + #ifdef SIGCONT
681 +      else signal(SIGCONT, pmapDistribReport);
682 + #endif
683 +   }
684 + #endif /* NIX */
685 +
686     /* ================================================================
687      * POST-DISTRIBUTION - Set photon flux and build kd-tree, etc.
688      * ================================================================ */
689 + #ifdef SIGCONT    
690 +   /* Reset signal handler */
691     signal(SIGCONT, SIG_DFL);
692 + #endif  
693     free(emap.samples);
694  
695 <   if (!pm -> heapEnd)
696 <      error(USER, "empty photon map");
695 >   if (!pm -> numPhotons)
696 >      error(USER, "empty contribution photon map");
697  
698 <   /* Check for valid primary photon rays */
699 <   if (!pm -> primary)
698 >   /* Load per-subprocess primary rays into pm -> primary array */
699 >   /* Dumb compilers apparently need the char** cast */
700 >   pm -> numPrimary = buildPrimaries(pm, primaryHeap,
701 >                                     (char**)primaryHeapFname,
702 >                                     primaryOfs, numProc);
703 >   if (!pm -> numPrimary)
704        error(INTERNAL, "no primary rays in contribution photon map");
481      
482   if (pm -> primary [pm -> primaryEnd].srcIdx < 0)
483      /* Last primary ray is unused, so decrement counter */
484      pm -> primaryEnd--;
705    
706 <   if (photonRepTime) {
707 <      eputs("\nBuilding contrib photon heap...\n");
706 >   /* Set photon flux per source */
707 >   for (srcIdx = 0; srcIdx < nsources; srcIdx++)
708 >      srcFlux [srcIdx] /= photonCnt [PHOTONCNT_NUMEMIT(srcIdx)];
709 > #if NIX
710 >   /* Photon counters no longer needed, unmap shared memory */
711 >   munmap(photonCnt, sizeof(*photonCnt));
712 >   close(shmFile);
713 >   unlink(shmFname);
714 > #else
715 >   free(photonCnt);  
716 > #endif      
717 >  
718 >   if (verbose) {
719 >      eputs("\nBuilding contribution photon map...\n");
720 > #if NIX      
721        fflush(stderr);
722 + #endif      
723     }
724 +  
725 +   /* Build underlying data structure; heap is destroyed */
726 +   buildPhotonMap(pm, srcFlux, primaryOfs, numProc);
727 +  
728 +   /* Free per-subprocess primary heap files */
729 +   for (proc = 0; proc < numProc; proc++)
730 +      free(primaryHeapFname [proc]);
731        
732 <   balancePhotons(pm, srcFlux);
732 >   free(primaryHeapFname);
733 >   free(primaryHeap);
734 >   free(primaryOfs);
735 >  
736 >   if (verbose)
737 >      eputs("\n");
738   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines