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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines