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.13 by greg, Thu Sep 29 21:51:58 2016 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 support for building 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  
31 < static void setPmapContribParams (PhotonMap *pmap, LUTAB *srcContrib)
32 < /* Set parameters for light source contributions */
31 > static PhotonPrimaryIdx newPhotonPrimary (PhotonMap *pmap,
32 >                                          const RAY *primRay,
33 >                                          FILE *primHeap)
34 > /* Add primary ray for emitted photon and save light source index, origin on
35 > * source, and emitted direction; used by contrib photons. The current
36 > * primary is stored in pmap -> lastPrimary.  If the previous primary
37 > * contributed photons (has srcIdx >= 0), it's appended to primHeap.  If
38 > * primRay == NULL, the current primary is still flushed, but no new primary
39 > * is set.  Returns updated primary counter pmap -> numPrimary.  */
40   {
41 <   /* Set light source modifier list and appropriate callback to extract
42 <    * their contributions from the photon map */
43 <   if (pmap) {
44 <      pmap -> srcContrib = srcContrib;
45 <      pmap -> lookup = photonContrib;
46 <      /* Ensure we get all requested photon contribs during lookups */
47 <      pmap -> gatherTolerance = 1.0;
41 >   if (!pmap || !primHeap)
42 >      return 0;
43 >      
44 >   /* Check if last primary ray has spawned photons (srcIdx >= 0, see
45 >    * newPhoton()), in which case we write it to the primary heap file
46 >    * before overwriting it */
47 >   if (pmap -> lastPrimary.srcIdx >= 0) {
48 >      if (!fwrite(&pmap -> lastPrimary, sizeof(PhotonPrimary), 1, primHeap))
49 >         error(SYSTEM, "failed writing photon primary in newPhotonPrimary");
50 >        
51 >      pmap -> numPrimary++;
52 >      if (pmap -> numPrimary > PMAP_MAXPRIMARY)
53 >         error(INTERNAL, "photon primary overflow in newPhotonPrimary");
54     }
37 }
55  
56 <
57 <
58 < static void checkPmapContribs (const PhotonMap *pmap, LUTAB *srcContrib)
59 < /* Check modifiers for light source contributions */
60 < {
61 <   const PhotonPrimary *primary = pmap -> primary;
62 <   OBJREC *srcMod;
63 <   unsigned long i, found = 0;
64 <  
65 <   /* Make sure at least one of the modifiers is actually in the pmap,
66 <    * otherwise findPhotons() winds up in an infinite loop! */
67 <   for (i = pmap -> primarySize; i; --i, ++primary) {
68 <      if (primary -> srcIdx < 0 || primary -> srcIdx >= nsources)
69 <         error(INTERNAL, "invalid light source index in photon map");
70 <        
54 <      srcMod = objptr(source [primary -> srcIdx].so -> omod);
55 <      if ((MODCONT*)lu_find(srcContrib, srcMod -> oname) -> data)
56 <         ++found;
56 >   /* Mark unused with negative source index until path spawns a photon (see
57 >    * newPhoton()) */
58 >   pmap -> lastPrimary.srcIdx = -1;
59 >    
60 >   if (primRay) {
61 >      FVECT dvec;
62 >      
63 >      /* Reverse incident direction to point to light source */
64 >      dvec [0] = -primRay -> rdir [0];
65 >      dvec [1] = -primRay -> rdir [1];
66 >      dvec [2] = -primRay -> rdir [2];
67 >      pmap -> lastPrimary.dir = encodedir(dvec);
68 > #ifdef PMAP_PRIMARYPOS      
69 >      VCOPY(pmap -> lastPrimary.pos, primRay -> rop);
70 > #endif      
71     }
72    
73 <   if (!found)
60 <      error(USER, "modifiers not in photon map");
73 >   return pmap -> numPrimary;
74   }
62
63  
75  
76 < void initPmapContrib (LUTAB *srcContrib, unsigned numSrcContrib)
76 >
77 >
78 > #ifdef DEBUG_PMAP_CONTRIB
79 > static int checkPrimaryHeap (FILE *file)
80 > /* Check heap for ordered primaries */
81   {
82 <   unsigned t;
82 >   Photon   p, lastp;
83 >   int      i, dup;
84    
85 <   for (t = 0; t < NUM_PMAP_TYPES; t++)
86 <      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 <      }
85 >   rewind(file);
86 >   memset(&lastp, 0, sizeof(lastp));
87    
88 <   /* Get params */
89 <   setPmapContribParams(contribPmap, srcContrib);
90 <  
91 <   if (contribPhotonMapping) {
92 <      if (contribPmap -> maxGather < numSrcContrib) {
93 <         /* Adjust density estimate bandwidth if lower than modifier
94 <          * count, otherwise contributions are missing */
95 <         error(WARNING, "contrib density estimate bandwidth too low, "
96 <                        "adjusting to modifier count");
97 <         contribPmap -> maxGather = numSrcContrib;
88 >   while (fread(&p, sizeof(p), 1, file)) {
89 >      dup = 1;
90 >      
91 >      for (i = 0; i <= 2; i++) {
92 >         if (p.pos [i] < thescene.cuorg [i] ||
93 >             p.pos [i] > thescene.cuorg [i] + thescene.cusize) {
94 >            
95 >            sprintf(errmsg, "corrupt photon in heap at [%f, %f, %f]\n",
96 >                    p.pos [0], p.pos [1], p.pos [2]);
97 >            error(WARNING, errmsg);
98 >         }
99 >        
100 >         dup &= p.pos [i] == lastp.pos [i];
101        }
102        
103 <      /* Sanity check */
104 <      checkPmapContribs(contribPmap, srcContrib);
103 >      if (dup) {
104 >         sprintf(errmsg,
105 >                 "consecutive duplicate photon in heap at [%f, %f, %f]\n",
106 >                 p.pos [0], p.pos [1], p.pos [2]);
107 >         error(WARNING, errmsg);
108 >      }
109     }
110 +
111 +   return 0;
112   }
113 + #endif
114  
115  
116  
117 < void photonContrib (PhotonMap *pmap, RAY *ray, COLOR irrad)
118 < /* Sum up light source contributions from photons in pmap->srcContrib */
117 > static PhotonPrimaryIdx buildPrimaries (PhotonMap *pmap, FILE **primaryHeap,
118 >                                        PhotonPrimaryIdx *primaryOfs,
119 >                                        unsigned numHeaps)
120 > /* Consolidate per-subprocess photon primary heaps into the primary array
121 > * pmap -> primaries. Returns offset for primary index linearisation in
122 > * numPrimary. The heap files in primaryHeap are closed on return. */
123   {
124 <   unsigned       i;
125 <   PhotonSQNode   *sq;
126 <   float          r, invArea;
127 <   RREAL          rayCoeff [3];
128 <
103 <   setcolor(irrad, 0, 0, 0);
104 <
105 <   if (!pmap -> maxGather)
106 <      return;
124 >   PhotonPrimaryIdx  heapLen;
125 >   unsigned          heap;
126 >  
127 >   if (!pmap || !primaryHeap || !primaryOfs || !numHeaps)
128 >      return 0;
129        
130 <   /* Ignore sources */
131 <   if (ray -> ro)
132 <      if (islight(objptr(ray -> ro -> omod) -> otype))
133 <         return;
130 >   pmap -> numPrimary = 0;
131 >  
132 >   for (heap = 0; heap < numHeaps; heap++) {
133 >      primaryOfs [heap] = pmap -> numPrimary;
134 >      
135 >      if (fseek(primaryHeap [heap], 0, SEEK_END))
136 >         error(SYSTEM, "failed photon primary seek in buildPrimaries");
137 >      pmap -> numPrimary += heapLen = ftell(primaryHeap [heap]) /
138 >                                      sizeof(PhotonPrimary);      
139  
140 <   /* Get cumulative path
141 <    * coefficient up to photon lookup point */
142 <   raycontrib(rayCoeff, ray, PRIMARY);
140 >      pmap -> primaries = realloc(pmap -> primaries,
141 >                                  pmap -> numPrimary *
142 >                                  sizeof(PhotonPrimary));  
143 >      if (!pmap -> primaries)
144 >         error(SYSTEM, "failed photon primary alloc in buildPrimaries");
145  
146 <   /* Lookup photons */
147 <   pmap -> squeueEnd = 0;
148 <   findPhotons(pmap, ray);
149 <  
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
146 >      rewind(primaryHeap [heap]);
147 >      if (fread(pmap -> primaries + primaryOfs [heap], sizeof(PhotonPrimary),
148 >                heapLen, primaryHeap [heap]) != heapLen)
149 >         error(SYSTEM, "failed reading photon primaries in buildPrimaries");
150        
151 <      return;
151 >      fclose(primaryHeap [heap]);      
152     }
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);
153    
154 <   /* Skip the extra photon */
155 <   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;
154 >   return pmap -> numPrimary;
155 > }      
156  
164         /* Photon's emitting light source has modifier whose
165          * contributions are sought */
166         double srcBinReal;
167         int srcBin;
168         RAY srcRay;
157  
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         }
158  
159 <         if ((srcBinReal = evalue(srcContrib -> binv)) < -.5)
160 <             continue;          /* silently ignore negative bins */
161 <  
189 <         if ((srcBin = srcBinReal + .5) >= srcContrib -> nbins) {
190 <             error(WARNING, "bad bin number (ignored)");
191 <             continue;
192 <         }
193 <            
194 <         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);
209 <      }
210 <   }
211 <        
212 <   return;
213 < }
159 > /* Defs for photon emission counter array passed by sub-processes to parent
160 > * via shared memory */
161 > typedef  unsigned long  PhotonContribCnt;
162  
163 + /* Indices for photon emission counter array: num photons stored and num
164 + * emitted per source */
165 + #define  PHOTONCNT_NUMPHOT    0
166 + #define  PHOTONCNT_NUMEMIT(n) (1 + n)
167  
168  
217 void distribPhotonContrib (PhotonMap* pm)
218 {
219   EmissionMap emap;
220   char errmsg2 [128];
221   unsigned srcIdx;
222   double *srcFlux;                 /* Emitted flux per light source */
223   const double srcDistribTarget =  /* Target photon count per source */
224      nsources ? (double)pm -> distribTarget / nsources : 0;  
169  
170 + void distribPhotonContrib (PhotonMap* pm, unsigned numProc)
171 + {
172 +   EmissionMap       emap;
173 +   char              errmsg2 [128], shmFname [255];
174 +   unsigned          srcIdx, proc;
175 +   int               shmFile, stat, pid;
176 +   double            *srcFlux,         /* Emitted flux per light source */
177 +                     srcDistribTarget; /* Target photon count per source */
178 +   PhotonContribCnt  *photonCnt;       /* Photon emission counter array */
179 +   const unsigned    photonCntSize = sizeof(PhotonContribCnt) *
180 +                                     PHOTONCNT_NUMEMIT(nsources);
181 +   FILE              *primaryHeap [numProc];
182 +   PhotonPrimaryIdx  primaryOfs [numProc];
183 +                                    
184     if (!pm)
185 <      error(USER, "no photon map defined");
185 >      error(USER, "no photon map defined in distribPhotonContrib");
186        
187     if (!nsources)
188 <      error(USER, "no light sources");
189 <  
188 >      error(USER, "no light sources in distribPhotonContrib");
189 >
190 >   if (nsources > MAXMODLIST)
191 >      error(USER, "too many light sources in distribPhotonContrib");
192 >      
193     /* Allocate photon flux per light source; this differs for every
194      * source as all sources contribute the same number of distributed
195      * photons (srcDistribTarget), hence the number of photons emitted per
196      * source does not correlate with its emitted flux. The resulting flux
197      * per photon is therefore adjusted individually for each source. */
198     if (!(srcFlux = calloc(nsources, sizeof(double))))
199 <      error(SYSTEM, "cannot allocate source flux");
199 >      error(SYSTEM, "can't allocate source flux in distribPhotonContrib");
200  
201 <   /* ================================================================
202 <    * INITIALISASHUNN - Set up emisshunn and scattering funcs
203 <    * ================================================================ */
201 >   /* ===================================================================
202 >    * INITIALISATION - Set up emission and scattering funcs
203 >    * =================================================================== */
204     emap.samples = NULL;
205     emap.src = NULL;
206     emap.maxPartitions = MAXSPART;
207     emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1);
208     if (!emap.partitions)
209 <      error(USER, "can't allocate source partitions");
209 >      error(USER, "can't allocate source partitions in distribPhotonContrib");
210  
211 +   /* Initialise contrib photon map */  
212     initPhotonMap(pm, PMAP_TYPE_CONTRIB);
213 +   initPhotonHeap(pm);
214     initPhotonEmissionFuncs();
215     initPhotonScatterFuncs();
216    
217 +   /* Per-subprocess / per-source target counts */
218 +   pm -> distribTarget /= numProc;
219 +   srcDistribTarget = nsources ? (double)pm -> distribTarget / nsources : 0;
220 +  
221     /* Get photon ports if specified */
222     if (ambincl == 1)
223        getPhotonPorts();
# Line 258 | Line 225 | void distribPhotonContrib (PhotonMap* pm)
225     /* Get photon sensor modifiers */
226     getPhotonSensors(photonSensorList);      
227    
228 <   /* Seed RNGs for photon distribution */
229 <   pmapSeed(randSeed, partState);
230 <   pmapSeed(randSeed, emitState);
231 <   pmapSeed(randSeed, cntState);
232 <   pmapSeed(randSeed, mediumState);
233 <   pmapSeed(randSeed, scatterState);
234 <   pmapSeed(randSeed, rouletteState);
235 <
269 <   /* Record start time and enable progress report signal handler */
270 <   repStartTime = time(NULL);
271 <   #ifdef SIGCONT  
272 <      signal(SIGCONT, pmapDistribReport);
273 <   #endif  
228 >   /* Set up shared mem for photon counters (zeroed by ftruncate) */
229 > #if 0  
230 >   snprintf(shmFname, 255, PMAP_SHMFNAME, getpid());
231 >   shmFile = shm_open(shmFname, O_CREAT | O_RDWR, S_IRUSR | S_IWUSR);
232 > #else
233 >   strcpy(shmFname, PMAP_SHMFNAME);
234 >   shmFile = mkstemp(shmFname);
235 > #endif        
236    
237 <   for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
238 <      unsigned portCnt = 0, passCnt = 0, prePassCnt = 0;
277 <      double srcNumEmit = 0;          /* # photons to emit from source */
278 <      unsigned long srcNumDistrib = pm -> heapEnd; /* # photons stored */
237 >   if (shmFile < 0 || ftruncate(shmFile, photonCntSize) < 0)
238 >      error(SYSTEM, "failed shared mem init in distribPhotonContrib");
239  
240 <      srcFlux [srcIdx] = repProgress = 0;
240 >   photonCnt = mmap(NULL, photonCntSize, PROT_READ | PROT_WRITE,
241 >                    MAP_SHARED, shmFile, 0);
242 >                    
243 >   if (photonCnt == MAP_FAILED)
244 >      error(SYSTEM, "failed shared mem mapping in distribPhotonContrib");
245 >
246 >   /* =============================================================
247 >    * FLUX INTEGRATION - Get total flux emitted from light source
248 >    * ============================================================= */  
249 >   for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
250 >      unsigned portCnt = 0;
251 >      
252 >      srcFlux [srcIdx] = 0;
253        emap.src = source + srcIdx;
254        
255        if (photonRepTime)
256           eputs("\n");
257 <      
258 <      /* =============================================================
259 <       * FLUX INTEGRATION - Get total flux emitted from light source
260 <       * ============================================================= */
289 <      do {
290 <         emap.port = emap.src -> sflags & SDISTANT
291 <                     ? photonPorts + portCnt : NULL;
257 >
258 >      do {  /* Need at least one iteration if no ports! */      
259 >         emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt
260 >                                                   : NULL;
261           photonPartition [emap.src -> so -> otype] (&emap);
262          
263           if (photonRepTime) {
# Line 302 | Line 271 | void distribPhotonContrib (PhotonMap* pm)
271                 strcat(errmsg, errmsg2);
272              }
273              
274 <            sprintf(errmsg2, "(%lu partitions)...\n",
306 <                    emap.numPartitions);
274 >            sprintf(errmsg2, "(%lu partitions)...\n", emap.numPartitions);
275              strcat(errmsg, errmsg2);
276              eputs(errmsg);
277              fflush(stderr);
278           }
279          
280 <         for (emap.partitionCnt = 0;
313 <              emap.partitionCnt < emap.numPartitions;
280 >         for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
281                emap.partitionCnt++) {
282              initPhotonEmission(&emap, pdfSamples);
283              srcFlux [srcIdx] += colorAvg(emap.partFlux);
284           }
285          
286           portCnt++;
287 <      } while (portCnt < numPhotonPorts);
288 <
287 >      } while (portCnt < numPhotonPorts);        
288 >      
289        if (srcFlux [srcIdx] < FTINY) {
290           sprintf(errmsg, "source %s has zero emission",
291                   source [srcIdx].so -> oname);
292           error(WARNING, errmsg);
293        }
294 <      else {
295 <         /* ==========================================================
294 >   }  
295 >
296 >   if (photonRepTime)
297 >      eputs("\n");
298 >
299 >   /* Init per-subprocess primary heap files */
300 >   for (proc = 0; proc < numProc; proc++)
301 >      if (!(primaryHeap [proc] = tmpfile()))
302 >         error(SYSTEM, "failed opening primary heap file in "
303 >               "distribPhotonContrib");
304 >              
305 >   /* MAIN LOOP */
306 >   for (proc = 0; proc < numProc; proc++) {
307 >      if (!(pid = fork())) {
308 >         /* SUBPROCESS ENTERS HERE;
309 >          * all opened and memory mapped files are inherited */
310 >          
311 >         /* Local photon counters for this subprocess */
312 >         unsigned long  lastNumPhotons = 0, localNumEmitted = 0;
313 >         double         photonFluxSum = 0;   /* Running photon flux sum */
314 >
315 >         /* Seed RNGs from PID for decorellated photon distribution */
316 >         pmapSeed(randSeed + proc, partState);
317 >         pmapSeed(randSeed + proc, emitState);
318 >         pmapSeed(randSeed + proc, cntState);
319 >         pmapSeed(randSeed + proc, mediumState);
320 >         pmapSeed(randSeed + proc, scatterState);
321 >         pmapSeed(randSeed + proc, rouletteState);
322 >        
323 >         /* =============================================================
324            * 2-PASS PHOTON DISTRIBUTION
325            * Pass 1 (pre):  emit fraction of target photon count
326 <          * Pass 2 (main): based on outcome of pass 1, estimate
327 <          *                remaining number of photons to emit to
328 <          *                approximate target count
329 <          * ========================================================== */
330 <         do {
331 <            if (!passCnt) {  
332 <               /* INIT PASS 1 */
333 <               if (++prePassCnt > maxPreDistrib) {
334 <                  /* Warn if no photons contributed after sufficient
335 <                   * iterations */
336 <                  sprintf(errmsg, "too many prepasses, no photons "
337 <                          "from source %s", source [srcIdx].so -> oname);
338 <                  error(WARNING, errmsg);
339 <                  break;
326 >          * Pass 2 (main): based on outcome of pass 1, estimate remaining
327 >          *                number of photons to emit to approximate target
328 >          *                count
329 >          * ============================================================= */        
330 >         for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
331 >            unsigned       portCnt, passCnt = 0, prePassCnt = 0;
332 >            float          srcPreDistrib = preDistrib;
333 >            double         srcNumEmit = 0;       /* # to emit from source */
334 >            unsigned long  srcNumDistrib = pm -> numPhotons;  /* # stored */
335 >
336 >            if (srcFlux [srcIdx] < FTINY)
337 >               continue;
338 >                        
339 >            while (passCnt < 2) {
340 >               if (!passCnt) {  
341 >                  /* INIT PASS 1 */
342 >                  if (++prePassCnt > maxPreDistrib) {
343 >                     /* Warn if no photons contributed after sufficient
344 >                      * iterations */
345 >                     sprintf(errmsg, "proc %d, source %s: "
346 >                             "too many prepasses, skipped",
347 >                             proc, source [srcIdx].so -> oname);
348 >                     error(WARNING, errmsg);
349 >                     break;
350 >                  }
351 >                  
352 >                  /* Num to emit is fraction of target count */
353 >                  srcNumEmit = srcPreDistrib * srcDistribTarget;
354                 }
355 <              
356 <               /* Num to emit is fraction of target count */
357 <               srcNumEmit = preDistrib * srcDistribTarget;
358 <            }
355 >               else {
356 >                  /* INIT PASS 2 */
357 >                  double srcPhotonFlux, avgPhotonFlux;
358 >                  
359 >                  /* Based on the outcome of the predistribution we can now
360 >                   * figure out how many more photons we have to emit from
361 >                   * the current source to meet the target count,
362 >                   * srcDistribTarget. This value is clamped to 0 in case
363 >                   * the target has already been exceeded in pass 1.
364 >                   * srcNumEmit and srcNumDistrib is the number of photons
365 >                   * emitted and distributed (stored) from the current
366 >                   * source in pass 1, respectively. */
367 >                  srcNumDistrib = pm -> numPhotons - srcNumDistrib;
368 >                  srcNumEmit *= srcNumDistrib
369 >                                ? max(srcDistribTarget/srcNumDistrib, 1) - 1
370 >                                : 0;
371  
372 <            else {            
373 <               /* INIT PASS 2 */                          
374 <               /* 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;
372 >                  if (!srcNumEmit)
373 >                     /* No photons left to distribute in main pass */
374 >                     break;
375                      
376 <            do {
377 <               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 <                  }
376 >                  srcPhotonFlux = srcFlux [srcIdx] / srcNumEmit;
377 >                  avgPhotonFlux = photonFluxSum / (srcIdx + 1);
378                    
379 <                  sprintf(errmsg2, "(%lu partitions)...\n",
380 <                          emap.numPartitions);
381 <                  strcat(errmsg, errmsg2);
382 <                  eputs(errmsg);
383 <                  fflush(stderr);
379 >                  if (avgPhotonFlux > 0 &&
380 >                      srcPhotonFlux / avgPhotonFlux < FTINY) {
381 >                     /* Skip source if its photon flux is grossly below the
382 >                      * running average, indicating negligible contribs at
383 >                      * the expense of excessive distribution time */
384 >                     sprintf(errmsg, "proc %d, source %s: "
385 >                             "itsy bitsy photon flux, skipped",
386 >                             proc, source [srcIdx].so -> oname);
387 >                     error(WARNING, errmsg);
388 >                     srcNumEmit = 0;
389 >                  }
390 >                        
391 >                  /* Update sum of photon flux per light source */
392 >                  photonFluxSum += srcPhotonFlux;
393                 }
394                
395 <               for (emap.partitionCnt = 0;
396 <                    emap.partitionCnt < emap.numPartitions;
397 <                    emap.partitionCnt++) {
398 <                  double partNumEmit;
399 <                  unsigned long partEmitCnt;
395 >               portCnt = 0;
396 >               do {    /* Need at least one iteration if no ports! */
397 >                  emap.src = source + srcIdx;
398 >                  emap.port = emap.src -> sflags & SDISTANT
399 >                              ? photonPorts + portCnt : NULL;
400 >                  photonPartition [emap.src -> so -> otype] (&emap);
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 >                  if (photonRepTime && !proc) {
403 >                     if (!passCnt)
404 >                        sprintf(errmsg, "PREPASS %d on source %s (mod %s) ",
405 >                                prePassCnt, source [srcIdx].so -> oname,
406 >                                objptr(source[srcIdx].so->omod) -> oname);
407 >                     else
408 >                        sprintf(errmsg, "MAIN PASS on source %s (mod %s) ",
409 >                                source [srcIdx].so -> oname,
410 >                                objptr(source[srcIdx].so->omod) -> oname);
411 >                            
412 >                     if (emap.port) {
413 >                        sprintf(errmsg2, "via port %s ",
414 >                                photonPorts [portCnt].so -> oname);
415 >                        strcat(errmsg, errmsg2);
416 >                     }
417                      
418 <                  /* Integer counter avoids FP rounding errors */
419 <                  while (partEmitCnt--) {                  
420 <                     RAY photonRay;
418 >                     sprintf(errmsg2, "(%lu partitions)\n",
419 >                             emap.numPartitions);
420 >                     strcat(errmsg, errmsg2);
421 >                     eputs(errmsg);
422 >                     fflush(stderr);
423 >                  }
424                    
425 <                     /* Emit photon according to PDF (if any), allocate
426 <                      * associated primary ray, and trace through scene
427 <                      * until absorbed/leaked */
428 <                     emitPhoton(&emap, &photonRay);
433 <                     addPhotonPrimary(pm, &photonRay);
434 <                     tracePhoton(&photonRay);
425 >                  for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
426 >                       emap.partitionCnt++) {
427 >                     double partNumEmit;
428 >                     unsigned long partEmitCnt;
429                      
430 <                     /* Record progress */
431 <                     repProgress++;
430 >                     /* Get photon origin within current source partishunn
431 >                      * and build emission map */
432 >                     photonOrigin [emap.src -> so -> otype] (&emap);
433 >                     initPhotonEmission(&emap, pdfSamples);
434                      
435 <                     if (photonRepTime > 0 &&
436 <                         time(NULL) >= repLastTime + photonRepTime)
437 <                        pmapDistribReport();
438 <                     #ifdef SIGCONT
439 <                        else signal(SIGCONT, pmapDistribReport);
440 <                     #endif
435 >                     /* Number of photons to emit from ziss partishunn;
436 >                      * scale according to its normalised contribushunn to
437 >                      * the emitted source flux */                    
438 >                     partNumEmit = srcNumEmit * colorAvg(emap.partFlux) /
439 >                                   srcFlux [srcIdx];                    
440 >                     partEmitCnt = (unsigned long)partNumEmit;
441 >                                                              
442 >                     /* Probabilistically account for fractional photons */
443 >                     if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
444 >                        partEmitCnt++;
445 >                        
446 >                     /* Update local and shared global emission counter */
447 >                     localNumEmitted += partEmitCnt;                                    
448 >                     photonCnt [PHOTONCNT_NUMEMIT(srcIdx)] += partEmitCnt;
449 >                    
450 >                     /* Integer counter avoids FP rounding errors */
451 >                     while (partEmitCnt--) {
452 >                        RAY photonRay;
453 >                    
454 >                        /* Emit photon according to PDF (if any), allocate
455 >                         * associated primary ray, and trace through scene
456 >                         * until absorbed/leaked; emitPhoton() sets the
457 >                         * emitting light source index in photonRay */
458 >                        emitPhoton(&emap, &photonRay);
459 >                        newPhotonPrimary(pm, &photonRay, primaryHeap[proc]);
460 >                        /* Set subprocess index in photonRay for post-
461 >                         * distrib primary index linearisation; this is
462 >                         * propagated with the primary index in photonRay
463 >                         * and set for photon hits by newPhoton() */
464 >                        PMAP_SETRAYPROC(&photonRay, proc);
465 >                        tracePhoton(&photonRay);
466 >                     }
467 >                    
468 >                     /* Update shared global photon count */                    
469 >                     photonCnt [PHOTONCNT_NUMPHOT] += pm -> numPhotons -
470 >                                                      lastNumPhotons;
471 >                     lastNumPhotons = pm -> numPhotons;
472                    }
446               }
447                          
448               portCnt++;
449            } while (portCnt < numPhotonPorts);
473  
474 <            if (pm -> heapEnd == srcNumDistrib)
475 <               /* Double preDistrib in case no photons were stored
476 <                * for this source and redo pass 1 */
477 <               preDistrib *= 2;
478 <            else {
479 <               /* Now do pass 2 */
480 <               passCnt++;
481 <               if (photonRepTime)
482 <                  eputs("\n");
474 >                  portCnt++;
475 >               } while (portCnt < numPhotonPorts);                  
476 >
477 >               if (pm -> numPhotons == srcNumDistrib)
478 >                  /* Double predistrib factor in case no photons were stored
479 >                   * for this source and redo pass 1 */
480 >                  srcPreDistrib *= 2;
481 >               else {
482 >                  /* Now do pass 2 */
483 >                  passCnt++;
484 > /*                if (photonRepTime)
485 >                     eputs("\n"); */
486 >               }
487              }
488 <         } while (passCnt < 2);
488 >         }
489 >                        
490 >         /* Flush heap buffa one final time to prevent data corruption */
491 >         flushPhotonHeap(pm);
492 >         fclose(pm -> heap);
493          
494 <         /* Flux per photon emitted from this source; repProgress is the
495 <          * number of emitted photons after both passes */
496 <         srcFlux [srcIdx] = repProgress ? srcFlux [srcIdx] / repProgress
497 <                                        : 0;
494 >         /* Flush final photon primary to primary heap file */
495 >         newPhotonPrimary(pm, NULL, primaryHeap [proc]);
496 >         fclose(primaryHeap [proc]);
497 >                  
498 > #ifdef DEBUG_PMAP
499 >         sprintf(errmsg, "Proc %d exited with total %ld photons\n", proc,
500 >                 pm -> numPhotons);
501 >         eputs(errmsg);
502 > #endif
503 >
504 >         exit(0);
505        }
506 +      else if (pid < 0)
507 +         error(SYSTEM, "failed to fork subprocess in distribPhotonContrib");
508     }
509  
510 +   /* PARENT PROCESS CONTINUES HERE */
511 +   /* Record start time and enable progress report signal handler */
512 +   repStartTime = time(NULL);
513 + #ifdef SIGCONT
514 +   signal(SIGCONT, pmapDistribReport);
515 + #endif
516 + /*
517 +   if (photonRepTime)
518 +      eputs("\n"); */
519 +  
520 +   /* Wait for subprocesses to complete while reporting progress */
521 +   proc = numProc;
522 +   while (proc) {
523 +      while (waitpid(-1, &stat, WNOHANG) > 0) {
524 +         /* Subprocess exited; check status */
525 +         if (!WIFEXITED(stat) || WEXITSTATUS(stat))
526 +            error(USER, "failed photon distribution");
527 +      
528 +         --proc;
529 +      }
530 +      
531 +      /* Nod off for a bit and update progress  */
532 +      sleep(1);
533 +      
534 +      /* Update progress report from shared subprocess counters */
535 +      repComplete = pm -> distribTarget * numProc;
536 +      repProgress = photonCnt [PHOTONCNT_NUMPHOT];
537 +      for (repEmitted = 0, srcIdx = 0; srcIdx < nsources; srcIdx++)
538 +         repEmitted += photonCnt [PHOTONCNT_NUMEMIT(srcIdx)];
539 +
540 +      /* Get global photon count from shmem updated by subprocs */
541 +      pm -> numPhotons = photonCnt [PHOTONCNT_NUMPHOT];
542 +
543 +      if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
544 +         pmapDistribReport();
545 + #ifdef SIGCONT
546 +      else signal(SIGCONT, pmapDistribReport);
547 + #endif
548 +   }
549 +
550     /* ================================================================
551      * POST-DISTRIBUTION - Set photon flux and build kd-tree, etc.
552      * ================================================================ */
553 <   #ifdef SIGCONT    
554 <      signal(SIGCONT, SIG_DFL);
555 <   #endif  
553 > #ifdef SIGCONT    
554 >   signal(SIGCONT, SIG_DFL);
555 > #endif  
556     free(emap.samples);
557  
558 <   if (!pm -> heapEnd)
558 >   if (!pm -> numPhotons)
559        error(USER, "empty photon map");
560  
561 <   /* Check for valid primary photon rays */
562 <   if (!pm -> primary)
561 >   /* Load per-subprocess primary rays into pm -> primary array */
562 >   pm -> numPrimary = buildPrimaries(pm, primaryHeap, primaryOfs, numProc);
563 >   if (!pm -> numPrimary)
564        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--;
565    
566 +   /* Set photon flux per source */
567 +   for (srcIdx = 0; srcIdx < nsources; srcIdx++)
568 +      srcFlux [srcIdx] /= photonCnt [PHOTONCNT_NUMEMIT(srcIdx)];
569 +
570 +   /* Photon counters no longer needed, unmap shared memory */
571 +   munmap(photonCnt, sizeof(*photonCnt));
572 +   close(shmFile);
573 + #if 0  
574 +   shm_unlink(shmFname);
575 + #else
576 +   unlink(shmFname);
577 + #endif      
578 +  
579     if (photonRepTime) {
580 <      eputs("\nBuilding contrib photon heap...\n");
580 >      eputs("\nBuilding contrib photon map...\n");
581        fflush(stderr);
582     }
583 <      
584 <   balancePhotons(pm, srcFlux);
583 >  
584 >   /* Build underlying data structure; heap is destroyed */
585 >   buildPhotonMap(pm, srcFlux, primaryOfs, numProc);  
586   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines