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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines