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.6 by greg, Wed May 20 13:43:28 2015 UTC vs.
Revision 2.12 by rschregle, Tue May 17 17:39:47 2016 UTC

# Line 1 | Line 1
1 + #ifndef lint
2 + static const char RCSid[] = "$Id$";
3 + #endif
4 +
5   /*
6 <   ==================================================================
6 >   ======================================================================
7     Photon map support 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 + #include <sys/mman.h>
27 + #include <sys/wait.h>
28  
29  
30  
# Line 41 | Line 46 | static void setPmapContribParams (PhotonMap *pmap, LUT
46   static void checkPmapContribs (const PhotonMap *pmap, LUTAB *srcContrib)
47   /* Check modifiers for light source contributions */
48   {
49 <   const PhotonPrimary *primary = pmap -> primary;
49 >   const PhotonPrimary  *primary = pmap -> primaries;
50 >   PhotonPrimaryIdx     i, found = 0;
51     OBJREC *srcMod;
46   unsigned long i, found = 0;
52    
53     /* Make sure at least one of the modifiers is actually in the pmap,
54      * otherwise findPhotons() winds up in an infinite loop! */
55 <   for (i = pmap -> primarySize; i; --i, ++primary) {
55 >   for (i = pmap -> numPrimary; 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);
59 >      srcMod = findmaterial(source [primary -> srcIdx].so);
60        if ((MODCONT*)lu_find(srcContrib, srcMod -> oname) -> data)
61           ++found;
62     }
# Line 92 | Line 97 | void initPmapContrib (LUTAB *srcContrib, unsigned numS
97  
98  
99  
100 < void photonContrib (PhotonMap *pmap, RAY *ray, COLOR irrad)
101 < /* Sum up light source contributions from photons in pmap->srcContrib */
100 > static PhotonPrimaryIdx newPhotonPrimary (PhotonMap *pmap,
101 >                                          const RAY *primRay,
102 >                                          FILE *primHeap)
103 > /* Add primary ray for emitted photon and save light source index, origin on
104 > * source, and emitted direction; used by contrib photons. The current
105 > * primary is stored in pmap -> lastPrimary.  If the previous primary
106 > * contributed photons (has srcIdx >= 0), it's appended to primHeap.  If
107 > * primRay == NULL, the current primary is still flushed, but no new primary
108 > * is set.  Returns updated primary counter pmap -> numPrimary.  */
109   {
110 <   unsigned       i;
111 <   PhotonSQNode   *sq;
100 <   float          r, invArea;
101 <   RREAL          rayCoeff [3];
102 <
103 <   setcolor(irrad, 0, 0, 0);
104 <
105 <   if (!pmap -> maxGather)
106 <      return;
110 >   if (!pmap || !primHeap)
111 >      return 0;
112        
113 <   /* Ignore sources */
114 <   if (ray -> ro)
115 <      if (islight(objptr(ray -> ro -> omod) -> otype))
116 <         return;
113 >   /* Check if last primary ray has spawned photons (srcIdx >= 0, see
114 >    * newPhoton()), in which case we write it to the primary heap file
115 >    * before overwriting it */
116 >   if (pmap -> lastPrimary.srcIdx >= 0) {
117 >      if (!fwrite(&pmap -> lastPrimary, sizeof(PhotonPrimary), 1, primHeap))
118 >         error(SYSTEM, "failed writing photon primary in newPhotonPrimary");
119 >        
120 >      pmap -> numPrimary++;
121 >      if (pmap -> numPrimary > PMAP_MAXPRIMARY)
122 >         error(INTERNAL, "photon primary overflow in newPhotonPrimary");
123 >   }
124  
125 <   /* Get cumulative path
126 <    * coefficient up to photon lookup point */
127 <   raycontrib(rayCoeff, ray, PRIMARY);
128 <
129 <   /* Lookup photons */
130 <   pmap -> squeueEnd = 0;
119 <   findPhotons(pmap, ray);
120 <  
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
125 >   /* Mark unused with negative source index until path spawns a photon (see
126 >    * newPhoton()) */
127 >   pmap -> lastPrimary.srcIdx = -1;
128 >    
129 >   if (primRay) {
130 >      FVECT dvec;
131        
132 <      return;
132 >      /* Reverse incident direction to point to light source */
133 >      dvec [0] = -primRay -> rdir [0];
134 >      dvec [1] = -primRay -> rdir [1];
135 >      dvec [2] = -primRay -> rdir [2];
136 >      pmap -> lastPrimary.dir = encodedir(dvec);
137 > #ifdef PMAP_PRIMARYPOS      
138 >      VCOPY(pmap -> lastPrimary.pos, primRay -> rop);
139 > #endif      
140     }
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);
141    
142 <   /* Skip the extra photon */
143 <   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 <         if (!srcContrib)
162 <            continue;
142 >   return pmap -> numPrimary;
143 > }
144  
164         /* Photon's emitting light source has modifier whose
165          * contributions are sought */
166         double srcBinReal;
167         int srcBin;
168         RAY srcRay;
145  
170         if (srcContrib -> binv -> type != NUM) {
171            /* Use intersection function to set shadow ray parameters
172             * if it's not simply a constant
173             */
174            rayorigin(&srcRay, SHADOW, NULL, NULL);
175            srcRay.rsrc = primary -> srcIdx;
176            VCOPY(srcRay.rorg, primary -> pos);
177            VCOPY(srcRay.rdir, primary -> dir);
178            if (!(source [primary -> srcIdx].sflags & SDISTANT ?
179                        sourcehit(&srcRay) :
180                        (*ofun[sp -> so -> otype].funp)(sp -> so, &srcRay)))
181                continue;               /* XXX shouldn't happen! */
182            worldfunc(RCCONTEXT, &srcRay);
183            set_eparams((char *)srcContrib -> params);
184         }
146  
147 <         if ((srcBinReal = evalue(srcContrib -> binv)) < -.5)
148 <             continue;          /* silently ignore negative bins */
149 <  
150 <         if ((srcBin = srcBinReal + .5) >= srcContrib -> nbins) {
151 <             error(WARNING, "bad bin number (ignored)");
152 <             continue;
147 > #ifdef DEBUG_PMAP_CONTRIB
148 > static int checkPrimaryHeap (FILE *file)
149 > /* Check heap for ordered primaries */
150 > {
151 >   Photon   p, lastp;
152 >   int      i, dup;
153 >  
154 >   rewind(file);
155 >   memset(&lastp, 0, sizeof(lastp));
156 >  
157 >   while (fread(&p, sizeof(p), 1, file)) {
158 >      dup = 1;
159 >      
160 >      for (i = 0; i <= 2; i++) {
161 >         if (p.pos [i] < thescene.cuorg [i] ||
162 >             p.pos [i] > thescene.cuorg [i] + thescene.cusize) {
163 >            
164 >            sprintf(errmsg, "corrupt photon in heap at [%f, %f, %f]\n",
165 >                    p.pos [0], p.pos [1], p.pos [2]);
166 >            error(WARNING, errmsg);
167           }
168 <            
169 <         if (!contrib) {
195 <             /* Ray coefficient mode; normalise by light source radiance
196 <              * after applying distrib pattern */
197 <             int j;
198 <             raytexture(ray, srcMod -> omod);
199 <             setcolor(ray -> rcol, srcMod -> oargs.farg [0],
200 <                        srcMod -> oargs.farg [1], srcMod -> oargs.farg [2]);
201 <             multcolor(ray -> rcol, ray -> pcol);
202 <             for (j = 0; j < 3; j++)
203 <                flux [j] = ray -> rcol [j] ? flux [j] / ray -> rcol [j]
204 <                                             : 0;
205 <         }
206 <                    
207 <         multcolor(flux, rayCoeff);
208 <         addcolor(srcContrib -> cbin [srcBin], flux);
168 >        
169 >         dup &= p.pos [i] == lastp.pos [i];
170        }
171 +      
172 +      if (dup) {
173 +         sprintf(errmsg,
174 +                 "consecutive duplicate photon in heap at [%f, %f, %f]\n",
175 +                 p.pos [0], p.pos [1], p.pos [2]);
176 +         error(WARNING, errmsg);
177 +      }
178     }
179 <        
180 <   return;
179 >
180 >   return 0;
181   }
182 + #endif
183  
184  
185  
186 < void distribPhotonContrib (PhotonMap* pm)
186 > static PhotonPrimaryIdx buildPrimaries (PhotonMap *pmap, FILE **primaryHeap,
187 >                                        PhotonPrimaryIdx *primaryOfs,
188 >                                        unsigned numHeaps)
189 > /* Consolidate per-subprocess photon primary heaps into the primary array
190 > * pmap -> primaries. Returns offset for primary index linearisation in
191 > * numPrimary. The heap files in primaryHeap are closed on return. */
192   {
193 <   EmissionMap emap;
194 <   char errmsg2 [128];
195 <   unsigned srcIdx;
196 <   double *srcFlux;                 /* Emitted flux per light source */
197 <   const double srcDistribTarget =  /* Target photon count per source */
198 <      nsources ? (double)pm -> distribTarget / nsources : 0;  
193 >   PhotonPrimaryIdx  heapLen;
194 >   unsigned          heap;
195 >  
196 >   if (!pmap || !primaryHeap || !primaryOfs || !numHeaps)
197 >      return 0;
198 >      
199 >   pmap -> numPrimary = 0;
200 >  
201 >   for (heap = 0; heap < numHeaps; heap++) {
202 >      primaryOfs [heap] = pmap -> numPrimary;
203 >      
204 >      if (fseek(primaryHeap [heap], 0, SEEK_END))
205 >         error(SYSTEM, "failed photon primary seek in buildPrimaries");
206 >      pmap -> numPrimary += heapLen = ftell(primaryHeap [heap]) /
207 >                                      sizeof(PhotonPrimary);      
208  
209 +      pmap -> primaries = realloc(pmap -> primaries,
210 +                                  pmap -> numPrimary *
211 +                                  sizeof(PhotonPrimary));  
212 +      if (!pmap -> primaries)
213 +         error(SYSTEM, "failed photon primary alloc in buildPrimaries");
214 +
215 +      rewind(primaryHeap [heap]);
216 +      if (fread(pmap -> primaries + primaryOfs [heap], sizeof(PhotonPrimary),
217 +                heapLen, primaryHeap [heap]) != heapLen)
218 +         error(SYSTEM, "failed reading photon primaries in buildPrimaries");
219 +      
220 +      fclose(primaryHeap [heap]);      
221 +   }
222 +  
223 +   return pmap -> numPrimary;
224 + }      
225 +
226 +
227 +
228 + /* Defs for photon emission counter array passed by sub-processes to parent
229 + * via shared memory */
230 + typedef  unsigned long  PhotonContribCnt;
231 +
232 + /* Indices for photon emission counter array: num photons stored and num
233 + * emitted per source */
234 + #define  PHOTONCNT_NUMPHOT    0
235 + #define  PHOTONCNT_NUMEMIT(n) (1 + n)
236 +
237 +
238 +
239 + void distribPhotonContrib (PhotonMap* pm, unsigned numProc)
240 + {
241 +   EmissionMap       emap;
242 +   char              errmsg2 [128], shmFname [255];
243 +   unsigned          srcIdx, proc;
244 +   int               shmFile, stat, pid;
245 +   double            *srcFlux,         /* Emitted flux per light source */
246 +                     srcDistribTarget; /* Target photon count per source */
247 +   PhotonContribCnt  *photonCnt;       /* Photon emission counter array */
248 +   const unsigned    photonCntSize = sizeof(PhotonContribCnt) *
249 +                                     PHOTONCNT_NUMEMIT(nsources);
250 +   FILE              *primaryHeap [numProc];
251 +   PhotonPrimaryIdx  primaryOfs [numProc];
252 +                                    
253     if (!pm)
254 <      error(USER, "no photon map defined");
254 >      error(USER, "no photon map defined in distribPhotonContrib");
255        
256     if (!nsources)
257 <      error(USER, "no light sources");
258 <  
257 >      error(USER, "no light sources in distribPhotonContrib");
258 >
259 >   if (nsources > MAXMODLIST)
260 >      error(USER, "too many light sources in distribPhotonContrib");
261 >      
262     /* Allocate photon flux per light source; this differs for every
263      * source as all sources contribute the same number of distributed
264      * photons (srcDistribTarget), hence the number of photons emitted per
265      * source does not correlate with its emitted flux. The resulting flux
266      * per photon is therefore adjusted individually for each source. */
267     if (!(srcFlux = calloc(nsources, sizeof(double))))
268 <      error(SYSTEM, "cannot allocate source flux");
268 >      error(SYSTEM, "can't allocate source flux in distribPhotonContrib");
269  
270 <   /* ================================================================
271 <    * INITIALISASHUNN - Set up emisshunn and scattering funcs
272 <    * ================================================================ */
270 >   /* ===================================================================
271 >    * INITIALISATION - Set up emission and scattering funcs
272 >    * =================================================================== */
273     emap.samples = NULL;
274     emap.src = NULL;
275     emap.maxPartitions = MAXSPART;
276     emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1);
277     if (!emap.partitions)
278 <      error(USER, "can't allocate source partitions");
278 >      error(USER, "can't allocate source partitions in distribPhotonContrib");
279  
280 +   /* Initialise contrib photon map */  
281     initPhotonMap(pm, PMAP_TYPE_CONTRIB);
282 +   initPhotonHeap(pm);
283     initPhotonEmissionFuncs();
284     initPhotonScatterFuncs();
285    
286 +   /* Per-subprocess / per-source target counts */
287 +   pm -> distribTarget /= numProc;
288 +   srcDistribTarget = nsources ? (double)pm -> distribTarget / nsources : 0;
289 +  
290     /* Get photon ports if specified */
291     if (ambincl == 1)
292        getPhotonPorts();
# Line 258 | Line 294 | void distribPhotonContrib (PhotonMap* pm)
294     /* Get photon sensor modifiers */
295     getPhotonSensors(photonSensorList);      
296    
297 <   /* Seed RNGs for photon distribution */
298 <   pmapSeed(randSeed, partState);
299 <   pmapSeed(randSeed, emitState);
300 <   pmapSeed(randSeed, cntState);
301 <   pmapSeed(randSeed, mediumState);
302 <   pmapSeed(randSeed, scatterState);
303 <   pmapSeed(randSeed, rouletteState);
304 <
269 <   /* Record start time and enable progress report signal handler */
270 <   repStartTime = time(NULL);
271 <   #ifdef SIGCONT  
272 <      signal(SIGCONT, pmapDistribReport);
273 <   #endif  
297 >   /* Set up shared mem for photon counters (zeroed by ftruncate) */
298 > #if 0  
299 >   snprintf(shmFname, 255, PMAP_SHMFNAME, getpid());
300 >   shmFile = shm_open(shmFname, O_CREAT | O_RDWR, S_IRUSR | S_IWUSR);
301 > #else
302 >   strcpy(shmFname, PMAP_SHMFNAME);
303 >   shmFile = mkstemp(shmFname);
304 > #endif        
305    
306 <   for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
307 <      unsigned portCnt = 0, passCnt = 0, prePassCnt = 0;
277 <      double srcNumEmit = 0;          /* # photons to emit from source */
278 <      unsigned long srcNumDistrib = pm -> heapEnd; /* # photons stored */
306 >   if (shmFile < 0 || ftruncate(shmFile, photonCntSize) < 0)
307 >      error(SYSTEM, "failed shared mem init in distribPhotonContrib");
308  
309 <      srcFlux [srcIdx] = repProgress = 0;
309 >   photonCnt = mmap(NULL, photonCntSize, PROT_READ | PROT_WRITE,
310 >                    MAP_SHARED, shmFile, 0);
311 >                    
312 >   if (photonCnt == MAP_FAILED)
313 >      error(SYSTEM, "failed shared mem mapping in distribPhotonContrib");
314 >
315 >   /* =============================================================
316 >    * FLUX INTEGRATION - Get total flux emitted from light source
317 >    * ============================================================= */  
318 >   for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
319 >      unsigned portCnt = 0;
320 >      
321 >      srcFlux [srcIdx] = 0;
322        emap.src = source + srcIdx;
323        
324        if (photonRepTime)
325           eputs("\n");
326 <      
327 <      /* =============================================================
328 <       * FLUX INTEGRATION - Get total flux emitted from light source
329 <       * ============================================================= */
289 <      do {
290 <         emap.port = emap.src -> sflags & SDISTANT
291 <                     ? photonPorts + portCnt : NULL;
326 >
327 >      do {  /* Need at least one iteration if no ports! */      
328 >         emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt
329 >                                                   : NULL;
330           photonPartition [emap.src -> so -> otype] (&emap);
331          
332           if (photonRepTime) {
# Line 302 | Line 340 | void distribPhotonContrib (PhotonMap* pm)
340                 strcat(errmsg, errmsg2);
341              }
342              
343 <            sprintf(errmsg2, "(%lu partitions)...\n",
306 <                    emap.numPartitions);
343 >            sprintf(errmsg2, "(%lu partitions)...\n", emap.numPartitions);
344              strcat(errmsg, errmsg2);
345              eputs(errmsg);
346              fflush(stderr);
347           }
348          
349 <         for (emap.partitionCnt = 0;
313 <              emap.partitionCnt < emap.numPartitions;
349 >         for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
350                emap.partitionCnt++) {
351              initPhotonEmission(&emap, pdfSamples);
352              srcFlux [srcIdx] += colorAvg(emap.partFlux);
353           }
354          
355           portCnt++;
356 <      } while (portCnt < numPhotonPorts);
357 <
356 >      } while (portCnt < numPhotonPorts);        
357 >      
358        if (srcFlux [srcIdx] < FTINY) {
359           sprintf(errmsg, "source %s has zero emission",
360                   source [srcIdx].so -> oname);
361           error(WARNING, errmsg);
362        }
363 <      else {
364 <         /* ==========================================================
363 >   }  
364 >
365 >   if (photonRepTime)
366 >      eputs("\n");
367 >
368 >   /* Init per-subprocess primary heap files */
369 >   for (proc = 0; proc < numProc; proc++)
370 >      if (!(primaryHeap [proc] = tmpfile()))
371 >         error(SYSTEM, "failed opening primary heap file in "
372 >               "distribPhotonContrib");
373 >              
374 >   /* MAIN LOOP */
375 >   for (proc = 0; proc < numProc; proc++) {
376 >      if (!(pid = fork())) {
377 >         /* SUBPROCESS ENTERS HERE;
378 >          * all opened and memory mapped files are inherited */
379 >          
380 >         /* Local photon counters for this subprocess */
381 >         unsigned long  lastNumPhotons = 0, localNumEmitted = 0;
382 >         double         photonFluxSum = 0;   /* Running photon flux sum */
383 >
384 >         /* Seed RNGs from PID for decorellated photon distribution */
385 >         pmapSeed(randSeed + proc, partState);
386 >         pmapSeed(randSeed + proc, emitState);
387 >         pmapSeed(randSeed + proc, cntState);
388 >         pmapSeed(randSeed + proc, mediumState);
389 >         pmapSeed(randSeed + proc, scatterState);
390 >         pmapSeed(randSeed + proc, rouletteState);
391 >        
392 >         /* =============================================================
393            * 2-PASS PHOTON DISTRIBUTION
394            * Pass 1 (pre):  emit fraction of target photon count
395 <          * Pass 2 (main): based on outcome of pass 1, estimate
396 <          *                remaining number of photons to emit to
397 <          *                approximate target count
398 <          * ========================================================== */
399 <         do {
400 <            if (!passCnt) {  
401 <               /* INIT PASS 1 */
402 <               if (++prePassCnt > maxPreDistrib) {
403 <                  /* Warn if no photons contributed after sufficient
404 <                   * iterations */
405 <                  sprintf(errmsg, "too many prepasses, no photons "
406 <                          "from source %s", source [srcIdx].so -> oname);
407 <                  error(WARNING, errmsg);
408 <                  break;
395 >          * Pass 2 (main): based on outcome of pass 1, estimate remaining
396 >          *                number of photons to emit to approximate target
397 >          *                count
398 >          * ============================================================= */        
399 >         for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
400 >            unsigned       portCnt, passCnt = 0, prePassCnt = 0;
401 >            float          srcPreDistrib = preDistrib;
402 >            double         srcNumEmit = 0;       /* # to emit from source */
403 >            unsigned long  srcNumDistrib = pm -> numPhotons;  /* # stored */
404 >
405 >            if (srcFlux [srcIdx] < FTINY)
406 >               continue;
407 >                        
408 >            while (passCnt < 2) {
409 >               if (!passCnt) {  
410 >                  /* INIT PASS 1 */
411 >                  if (++prePassCnt > maxPreDistrib) {
412 >                     /* Warn if no photons contributed after sufficient
413 >                      * iterations */
414 >                     sprintf(errmsg, "proc %d, source %s: "
415 >                             "too many prepasses, skipped",
416 >                             proc, source [srcIdx].so -> oname);
417 >                     error(WARNING, errmsg);
418 >                     break;
419 >                  }
420 >                  
421 >                  /* Num to emit is fraction of target count */
422 >                  srcNumEmit = srcPreDistrib * srcDistribTarget;
423                 }
424 <              
425 <               /* Num to emit is fraction of target count */
426 <               srcNumEmit = preDistrib * srcDistribTarget;
427 <            }
424 >               else {
425 >                  /* INIT PASS 2 */
426 >                  double srcPhotonFlux, avgPhotonFlux;
427 >                  
428 >                  /* Based on the outcome of the predistribution we can now
429 >                   * figure out how many more photons we have to emit from
430 >                   * the current source to meet the target count,
431 >                   * srcDistribTarget. This value is clamped to 0 in case
432 >                   * the target has already been exceeded in pass 1.
433 >                   * srcNumEmit and srcNumDistrib is the number of photons
434 >                   * emitted and distributed (stored) from the current
435 >                   * source in pass 1, respectively. */
436 >                  srcNumDistrib = pm -> numPhotons - srcNumDistrib;
437 >                  srcNumEmit *= srcNumDistrib
438 >                                ? max(srcDistribTarget/srcNumDistrib, 1) - 1
439 >                                : 0;
440  
441 <            else {            
442 <               /* INIT PASS 2 */                          
443 <               /* Based on the outcome of the predistribution we can now
354 <                * figure out how many more photons we have to emit from
355 <                * the current source to meet the target count,
356 <                * srcDistribTarget. This value is clamped to 0 in case
357 <                * the target has already been exceeded in pass 1.
358 <                * srcNumEmit and srcNumDistrib is the number of photons
359 <                * emitted and distributed (stored) from the current
360 <                * source in pass 1, respectively. */
361 <               srcNumDistrib = pm -> heapEnd - srcNumDistrib;
362 <               srcNumEmit *= srcNumDistrib
363 <                             ? max(srcDistribTarget/srcNumDistrib, 1) - 1
364 <                             : 0;
365 <
366 <               if (!srcNumEmit)
367 <                  /* No photons left to distribute in main pass */
368 <                  break;
369 <            }
370 <        
371 <            /* Set completion count for progress report */
372 <            repComplete = srcNumEmit + repProgress;
373 <            portCnt = 0;
441 >                  if (!srcNumEmit)
442 >                     /* No photons left to distribute in main pass */
443 >                     break;
444                      
445 <            do {
446 <               emap.port = emap.src -> sflags & SDISTANT
377 <                           ? photonPorts + portCnt : NULL;
378 <               photonPartition [emap.src -> so -> otype] (&emap);
379 <              
380 <               if (photonRepTime) {
381 <                  if (!passCnt)
382 <                     sprintf(errmsg, "PREPASS %d on source %s (mod %s) ",
383 <                             prePassCnt, source [srcIdx].so -> oname,
384 <                             objptr(source[srcIdx].so->omod) -> oname);
385 <                  else
386 <                     sprintf(errmsg, "MAIN PASS on source %s (mod %s) ",
387 <                             source [srcIdx].so -> oname,
388 <                             objptr(source[srcIdx].so->omod) -> oname);
389 <                          
390 <                  if (emap.port) {
391 <                     sprintf(errmsg2, "via port %s ",
392 <                             photonPorts [portCnt].so -> oname);
393 <                     strcat(errmsg, errmsg2);
394 <                  }
445 >                  srcPhotonFlux = srcFlux [srcIdx] / srcNumEmit;
446 >                  avgPhotonFlux = photonFluxSum / (srcIdx + 1);
447                    
448 <                  sprintf(errmsg2, "(%lu partitions)...\n",
449 <                          emap.numPartitions);
450 <                  strcat(errmsg, errmsg2);
451 <                  eputs(errmsg);
452 <                  fflush(stderr);
448 >                  if (avgPhotonFlux > 0 &&
449 >                      srcPhotonFlux / avgPhotonFlux < FTINY) {
450 >                     /* Skip source if its photon flux is grossly below the
451 >                      * running average, indicating negligible contribs at
452 >                      * the expense of excessive distribution time */
453 >                     sprintf(errmsg, "proc %d, source %s: "
454 >                             "itsy bitsy photon flux, skipped",
455 >                             proc, source [srcIdx].so -> oname);
456 >                     error(WARNING, errmsg);
457 >                     srcNumEmit = 0;
458 >                  }
459 >                        
460 >                  /* Update sum of photon flux per light source */
461 >                  photonFluxSum += srcPhotonFlux;
462                 }
463                
464 <               for (emap.partitionCnt = 0;
465 <                    emap.partitionCnt < emap.numPartitions;
466 <                    emap.partitionCnt++) {
467 <                  double partNumEmit;
468 <                  unsigned long partEmitCnt;
464 >               portCnt = 0;
465 >               do {    /* Need at least one iteration if no ports! */
466 >                  emap.src = source + srcIdx;
467 >                  emap.port = emap.src -> sflags & SDISTANT
468 >                              ? photonPorts + portCnt : NULL;
469 >                  photonPartition [emap.src -> so -> otype] (&emap);
470                    
471 <                  /* Get photon origin within current source partishunn
472 <                   * and build emission map */
473 <                  photonOrigin [emap.src -> so -> otype] (&emap);
474 <                  initPhotonEmission(&emap, pdfSamples);
475 <                  
476 <                  /* Number of photons to emit from ziss partishunn;
477 <                   * scale according to its normalised contribushunn to
478 <                   * the emitted source flux */
479 <                  partNumEmit = srcNumEmit * colorAvg(emap.partFlux) /
480 <                                srcFlux [srcIdx];
481 <                  partEmitCnt = (unsigned long)partNumEmit;
482 <                                      
483 <                  /* Probabilistically account for fractional photons */
484 <                  if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
485 <                     partEmitCnt++;
471 >                  if (photonRepTime && !proc) {
472 >                     if (!passCnt)
473 >                        sprintf(errmsg, "PREPASS %d on source %s (mod %s) ",
474 >                                prePassCnt, source [srcIdx].so -> oname,
475 >                                objptr(source[srcIdx].so->omod) -> oname);
476 >                     else
477 >                        sprintf(errmsg, "MAIN PASS on source %s (mod %s) ",
478 >                                source [srcIdx].so -> oname,
479 >                                objptr(source[srcIdx].so->omod) -> oname);
480 >                            
481 >                     if (emap.port) {
482 >                        sprintf(errmsg2, "via port %s ",
483 >                                photonPorts [portCnt].so -> oname);
484 >                        strcat(errmsg, errmsg2);
485 >                     }
486                      
487 <                  /* Integer counter avoids FP rounding errors */
488 <                  while (partEmitCnt--) {                  
489 <                     RAY photonRay;
487 >                     sprintf(errmsg2, "(%lu partitions)\n",
488 >                             emap.numPartitions);
489 >                     strcat(errmsg, errmsg2);
490 >                     eputs(errmsg);
491 >                     fflush(stderr);
492 >                  }
493                    
494 <                     /* Emit photon according to PDF (if any), allocate
495 <                      * associated primary ray, and trace through scene
496 <                      * until absorbed/leaked */
497 <                     emitPhoton(&emap, &photonRay);
433 <                     addPhotonPrimary(pm, &photonRay);
434 <                     tracePhoton(&photonRay);
494 >                  for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
495 >                       emap.partitionCnt++) {
496 >                     double partNumEmit;
497 >                     unsigned long partEmitCnt;
498                      
499 <                     /* Record progress */
500 <                     repProgress++;
499 >                     /* Get photon origin within current source partishunn
500 >                      * and build emission map */
501 >                     photonOrigin [emap.src -> so -> otype] (&emap);
502 >                     initPhotonEmission(&emap, pdfSamples);
503                      
504 <                     if (photonRepTime > 0 &&
505 <                         time(NULL) >= repLastTime + photonRepTime)
506 <                        pmapDistribReport();
507 <                     #ifdef SIGCONT
508 <                        else signal(SIGCONT, pmapDistribReport);
509 <                     #endif
504 >                     /* Number of photons to emit from ziss partishunn;
505 >                      * scale according to its normalised contribushunn to
506 >                      * the emitted source flux */                    
507 >                     partNumEmit = srcNumEmit * colorAvg(emap.partFlux) /
508 >                                   srcFlux [srcIdx];                    
509 >                     partEmitCnt = (unsigned long)partNumEmit;
510 >                                                              
511 >                     /* Probabilistically account for fractional photons */
512 >                     if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
513 >                        partEmitCnt++;
514 >                        
515 >                     /* Update local and shared global emission counter */
516 >                     localNumEmitted += partEmitCnt;                                    
517 >                     photonCnt [PHOTONCNT_NUMEMIT(srcIdx)] += partEmitCnt;
518 >                    
519 >                     /* Integer counter avoids FP rounding errors */
520 >                     while (partEmitCnt--) {
521 >                        RAY photonRay;
522 >                    
523 >                        /* Emit photon according to PDF (if any), allocate
524 >                         * associated primary ray, and trace through scene
525 >                         * until absorbed/leaked; emitPhoton() sets the
526 >                         * emitting light source index in photonRay */
527 >                        emitPhoton(&emap, &photonRay);
528 >                        newPhotonPrimary(pm, &photonRay, primaryHeap[proc]);
529 >                        /* Set subprocess index in photonRay for post-
530 >                         * distrib primary index linearisation; this is
531 >                         * propagated with the primary index in photonRay
532 >                         * and set for photon hits by newPhoton() */
533 >                        PMAP_SETRAYPROC(&photonRay, proc);
534 >                        tracePhoton(&photonRay);
535 >                     }
536 >                    
537 >                     /* Update shared global photon count */                    
538 >                     photonCnt [PHOTONCNT_NUMPHOT] += pm -> numPhotons -
539 >                                                      lastNumPhotons;
540 >                     lastNumPhotons = pm -> numPhotons;
541                    }
446               }
447                          
448               portCnt++;
449            } 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 >               else {
551 >                  /* Now do pass 2 */
552 >                  passCnt++;
553 > /*                if (photonRepTime)
554 >                     eputs("\n"); */
555 >               }
556              }
557 <         } while (passCnt < 2);
557 >         }
558 >                        
559 >         /* Flush heap buffa one final time to prevent data corruption */
560 >         flushPhotonHeap(pm);
561 >         fclose(pm -> heap);
562          
563 <         /* Flux per photon emitted from this source; repProgress is the
564 <          * number of emitted photons after both passes */
565 <         srcFlux [srcIdx] = repProgress ? srcFlux [srcIdx] / repProgress
566 <                                        : 0;
563 >         /* Flush final photon primary to primary heap file */
564 >         newPhotonPrimary(pm, NULL, primaryHeap [proc]);
565 >         fclose(primaryHeap [proc]);
566 >                  
567 > #ifdef DEBUG_PMAP
568 >         sprintf(errmsg, "Proc %d exited with total %ld photons\n", proc,
569 >                 pm -> numPhotons);
570 >         eputs(errmsg);
571 > #endif
572 >
573 >         exit(0);
574        }
575 +      else if (pid < 0)
576 +         error(SYSTEM, "failed to fork subprocess in distribPhotonContrib");
577     }
578  
579 +   /* PARENT PROCESS CONTINUES HERE */
580 +   /* Record start time and enable progress report signal handler */
581 +   repStartTime = time(NULL);
582 + #ifdef SIGCONT
583 +   signal(SIGCONT, pmapDistribReport);
584 + #endif
585 + /*
586 +   if (photonRepTime)
587 +      eputs("\n"); */
588 +  
589 +   /* Wait for subprocesses to complete while reporting progress */
590 +   proc = numProc;
591 +   while (proc) {
592 +      while (waitpid(-1, &stat, WNOHANG) > 0) {
593 +         /* Subprocess exited; check status */
594 +         if (!WIFEXITED(stat) || WEXITSTATUS(stat))
595 +            error(USER, "failed photon distribution");
596 +      
597 +         --proc;
598 +      }
599 +      
600 +      /* Nod off for a bit and update progress  */
601 +      sleep(1);
602 +      
603 +      /* Update progress report from shared subprocess counters */
604 +      repComplete = pm -> distribTarget * numProc;
605 +      repProgress = photonCnt [PHOTONCNT_NUMPHOT];
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 +
619     /* ================================================================
620      * POST-DISTRIBUTION - Set photon flux and build kd-tree, etc.
621      * ================================================================ */
622 <   #ifdef SIGCONT    
623 <      signal(SIGCONT, SIG_DFL);
624 <   #endif  
622 > #ifdef SIGCONT    
623 >   signal(SIGCONT, SIG_DFL);
624 > #endif  
625     free(emap.samples);
626  
627 <   if (!pm -> heapEnd)
627 >   if (!pm -> numPhotons)
628        error(USER, "empty photon map");
629  
630 <   /* Check for valid primary photon rays */
631 <   if (!pm -> primary)
630 >   /* Load per-subprocess primary rays into pm -> primary array */
631 >   pm -> numPrimary = buildPrimaries(pm, primaryHeap, primaryOfs, numProc);
632 >   if (!pm -> numPrimary)
633        error(INTERNAL, "no primary rays in contribution photon map");
484      
485   if (pm -> primary [pm -> primaryEnd].srcIdx < 0)
486      /* Last primary ray is unused, so decrement counter */
487      pm -> primaryEnd--;
634    
635 +   /* Set photon flux per source */
636 +   for (srcIdx = 0; srcIdx < nsources; srcIdx++)
637 +      srcFlux [srcIdx] /= photonCnt [PHOTONCNT_NUMEMIT(srcIdx)];
638 +
639 +   /* Photon counters no longer needed, unmap shared memory */
640 +   munmap(photonCnt, sizeof(*photonCnt));
641 +   close(shmFile);
642 + #if 0  
643 +   shm_unlink(shmFname);
644 + #else
645 +   unlink(shmFname);
646 + #endif      
647 +  
648     if (photonRepTime) {
649 <      eputs("\nBuilding contrib photon heap...\n");
649 >      eputs("\nBuilding contrib photon map...\n");
650        fflush(stderr);
651     }
652 +  
653 +   /* Build underlying data structure; heap is destroyed */
654 +   buildPhotonMap(pm, srcFlux, primaryOfs, numProc);  
655 + }
656 +
657 +
658 +
659 + void photonContrib (PhotonMap *pmap, RAY *ray, COLOR irrad)
660 + /* Sum up light source contributions from photons in pmap->srcContrib */
661 + {
662 +   unsigned                i;
663 +   PhotonSearchQueueNode   *sqn;
664 +   float                   r, invArea;
665 +   RREAL                   rayCoeff [3];
666 +   Photon                  *photon;
667 +   static char             warn = 1;
668 +
669 +   setcolor(irrad, 0, 0, 0);
670 +
671 +   if (!pmap -> maxGather)
672 +      return;
673        
674 <   balancePhotons(pm, srcFlux);
674 >   /* Ignore sources */
675 >   if (ray -> ro && islight(objptr(ray -> ro -> omod) -> otype))
676 >      return;
677 >
678 >   /* Get cumulative path coefficient up to photon lookup point */
679 >   raycontrib(rayCoeff, ray, PRIMARY);
680 >
681 >   /* Lookup photons */
682 >   pmap -> squeue.tail = 0;
683 >   findPhotons(pmap, ray);
684 >  
685 >   /* Need at least 2 photons */
686 >   if (pmap -> squeue.tail < 2) {
687 > #ifdef PMAP_NONEFOUND
688 >      sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)",
689 >              ray -> ro ? ray -> ro -> oname : "<null>",
690 >              ray -> rop [0], ray -> rop [1], ray -> rop [2]);
691 >      error(WARNING, errmsg);
692 > #endif
693 >      
694 >      return;
695 >   }
696 >
697 >   /* Average (squared) radius between furthest two photons to improve
698 >    * accuracy and get inverse search area 1 / (PI * r^2), with extra
699 >    * normalisation factor 1 / PI for ambient calculation */
700 >   sqn = pmap -> squeue.node + 1;
701 >   r = max(sqn -> dist2, (sqn + 1) -> dist2);
702 >   r = 0.25 * (pmap -> maxDist2 + r + 2 * sqrt(pmap -> maxDist2 * r));  
703 >   invArea = 1 / (PI * PI * r);
704 >  
705 >   /* Skip the extra photon */
706 >   for (i = 1 ; i < pmap -> squeue.tail; i++, sqn++) {        
707 >      COLOR flux;
708 >      
709 >      /* Get photon's contribution to density estimate */
710 >      photon = getNearestPhoton(&pmap -> squeue, sqn -> idx);
711 >      getPhotonFlux(photon, flux);
712 >      scalecolor(flux, invArea);      
713 > #ifdef PMAP_EPANECHNIKOV
714 >      /* Apply Epanechnikov kernel to photon flux based on photon distance */
715 >      scalecolor(flux, 2 * (1 - sqn -> dist2 / r));
716 > #endif
717 >      addcolor(irrad, flux);
718 >      
719 >      if (pmap -> srcContrib) {      
720 >         const PhotonPrimary  *primary = pmap -> primaries +
721 >                                         photon -> primary;
722 >         const SRCREC         *sp = &source [primary -> srcIdx];
723 >         OBJREC               *srcMod = findmaterial(sp -> so);
724 >         MODCONT *srcContrib = (MODCONT*)lu_find(pmap -> srcContrib,
725 >                                                 srcMod -> oname) -> data;
726 >         double   srcBinReal;
727 >         int      srcBin;
728 >         RAY      srcRay;                                                
729 >                                                
730 >         if (!srcContrib)
731 >            continue;
732 >
733 >         /* Photon's emitting light source has modifier whose contributions
734 >          * are sought */
735 >         if (srcContrib -> binv -> type != NUM) {
736 >            /* Use intersection function to set shadow ray parameters if
737 >             * it's not simply a constant */
738 >            rayorigin(&srcRay, SHADOW, NULL, NULL);
739 >            srcRay.rsrc = primary -> srcIdx;
740 > #ifdef PMAP_PRIMARYPOS
741 >            VCOPY(srcRay.rorg, primary -> pos);
742 > #else
743 >            /* No primary hitpoints; set dummy ray origin and warn once */
744 >            srcRay.rorg [0] = srcRay.rorg [1] = srcRay.rorg [2] = 0;
745 >            if (warn) {
746 >               error(WARNING, "no photon primary hitpoints for bin evaluation;"
747 >                     " using dummy (0,0,0) !");
748 >               warn = 0;
749 >            }
750 > #endif          
751 >            decodedir(srcRay.rdir, primary -> dir);
752 >
753 >            if (!(sp->sflags & SDISTANT
754 >                  ? sourcehit(&srcRay)
755 >                  : (*ofun[sp -> so -> otype].funp)(sp -> so, &srcRay)))
756 >               continue;                /* XXX shouldn't happen! */
757 >
758 >            worldfunc(RCCONTEXT, &srcRay);
759 >            set_eparams((char *)srcContrib -> params);
760 >         }
761 >
762 >         if ((srcBinReal = evalue(srcContrib -> binv)) < -.5)
763 >            continue;           /* silently ignore negative bins */
764 >  
765 >         if ((srcBin = srcBinReal + .5) >= srcContrib -> nbins) {
766 >            error(WARNING, "bad bin number (ignored)");
767 >            continue;
768 >         }
769 >            
770 >         if (!contrib) {
771 >            /* Ray coefficient mode; normalise by light source radiance
772 >             * after applying distrib pattern */
773 >            int j;
774 >            
775 >            raytexture(ray, srcMod -> omod);
776 >            setcolor(ray -> rcol, srcMod -> oargs.farg [0],
777 >                     srcMod -> oargs.farg [1], srcMod -> oargs.farg [2]);
778 >            multcolor(ray -> rcol, ray -> pcol);
779 >            for (j = 0; j < 3; j++)
780 >               flux [j] = ray -> rcol [j] ? flux [j] / ray -> rcol [j] : 0;
781 >         }
782 >                    
783 >         multcolor(flux, rayCoeff);
784 >         addcolor(srcContrib -> cbin [srcBin], flux);
785 >      }
786 >   }
787 >        
788 >   return;
789   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines