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

Comparing ray/src/rt/pmap.c (file contents):
Revision 2.2 by greg, Tue Apr 21 19:16:51 2015 UTC vs.
Revision 2.17 by rschregle, Tue Dec 18 22:14:04 2018 UTC

# Line 1 | Line 1
1 + #ifndef lint
2 + static const char RCSid[] = "$Id$";
3 + #endif
4 +
5 +
6   /*
7 <   ==================================================================
7 >   ======================================================================
8     Photon map main module
9  
10     Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
11     (c) Fraunhofer Institute for Solar Energy Systems,
12 <       Lucerne University of Applied Sciences & Arts  
13 <   ==================================================================
12 >   (c) Lucerne University of Applied Sciences and Arts,
13 >       supported by the Swiss National Science Foundation (SNSF, #147053)
14 >   ======================================================================
15    
16     $Id$
17   */
18  
19  
14 –
20   #include "pmap.h"
21   #include "pmapmat.h"
22   #include "pmapsrc.h"
# Line 21 | Line 26
26   #include "pmapdiag.h"
27   #include "otypes.h"
28   #include <time.h>
29 < #include <sys/stat.h>
29 > #if NIX
30 >   #include <sys/stat.h>
31 >   #include <sys/mman.h>
32 >   #include <sys/wait.h>
33 > #endif
34  
35  
27 –
28 – extern char *octname;
29 –
30 – static char PmapRevision [] = "$Revision$";
31 –
32 –
33 –
34 – /* Photon map lookup functions per type */
35 – void (*pmapLookup [NUM_PMAP_TYPES])(PhotonMap*, RAY*, COLOR) = {
36 –   photonDensity, photonPreCompDensity, photonDensity, volumePhotonDensity,
37 –   photonDensity, NULL
38 – };
39 –
40 –
41 –
42 – void colorNorm (COLOR c)
43 – /* Normalise colour channels to average of 1 */
44 – {
45 –   const float avg = colorAvg(c);
46 –  
47 –   if (!avg)
48 –      return;
49 –      
50 –   c [0] /= avg;
51 –   c [1] /= avg;
52 –   c [2] /= avg;
53 – }
54 –
55 –
56 –
57 – void loadPmaps (PhotonMap **pmaps, const PhotonMapParams *parm)
58 – {
59 –   unsigned t;
60 –   struct stat octstat, pmstat;
61 –   PhotonMap *pm;
62 –   PhotonMapType type;
63 –  
64 –   for (t = 0; t < NUM_PMAP_TYPES; t++)
65 –      if (setPmapParam(&pm, parm + t)) {        
66 –         /* Check if photon map newer than octree */
67 –         if (!stat(pm -> fileName, &pmstat) && !stat(octname, &octstat) &&
68 –             octstat.st_mtime > pmstat.st_mtime) {
69 –            sprintf(errmsg, "photon map in file %s may be stale",
70 –                    pm -> fileName);
71 –            error(USER, errmsg);
72 –         }
73 –        
74 –         /* Load photon map from file and get its type */
75 –         if ((type = loadPhotonMap(pm, pm -> fileName)) == PMAP_TYPE_NONE)
76 –            error(USER, "failed loading photon map");
77 –            
78 –         /* Assign to appropriate photon map type (deleting previously
79 –          * loaded photon map of same type if necessary) */
80 –         if (pmaps [type]) {
81 –            deletePhotons(pmaps [type]);
82 –            free(pmaps [type]);
83 –         }
84 –         pmaps [type] = pm;
85 –        
86 –         /* Check for invalid density estimate bandwidth */                            
87 –         if (pm -> maxGather > pm -> heapSize) {
88 –            error(WARNING, "adjusting density estimate bandwidth");
89 –            pm -> minGather = pm -> maxGather = pm -> heapSize;
90 –         }
91 –      }
92 – }
93 –
94 –
95 –
36   void savePmaps (const PhotonMap **pmaps, int argc, char **argv)
37   {
38     unsigned t;
39    
40     for (t = 0; t < NUM_PMAP_TYPES; t++) {
41        if (pmaps [t])
42 <         savePhotonMap(pmaps [t], pmaps [t] -> fileName, t, argc, argv);
42 >         savePhotonMap(pmaps [t], pmaps [t] -> fileName, argc, argv);
43     }
44   }                  
45  
46  
107 –  
108 – void cleanUpPmaps (PhotonMap **pmaps)
109 – {
110 –   unsigned t;
111 –  
112 –   for (t = 0; t < NUM_PMAP_TYPES; t++) {
113 –      if (pmaps [t]) {
114 –         deletePhotons(pmaps [t]);
115 –         free(pmaps [t]);
116 –      }
117 –   }
118 – }
119 –
120 –
47      
48   static int photonParticipate (RAY *ray)
49   /* Trace photon through participating medium. Returns 1 if passed through,
50     or 0 if absorbed and $*%&ed. Analogon to rayparticipate(). */
51   {
52     int i;
53 <   RREAL cosTheta, cosPhi, du, dv;
53 >   RREAL xi1, cosTheta, phi, du, dv;
54     const float cext = colorAvg(ray -> cext),
55 <               albedo = colorAvg(ray -> albedo);
55 >               albedo = colorAvg(ray -> albedo),
56 >               gecc = ray -> gecc, gecc2 = sqr(gecc);
57     FVECT u, v;
58     COLOR cvext;
59  
# Line 134 | Line 61 | static int photonParticipate (RAY *ray)
61     ray -> rmax = -log(pmapRandom(mediumState)) / cext;
62    
63     while (!localhit(ray, &thescene)) {
64 +      if (!incube(&thescene, ray -> rop)) {
65 +         /* Terminate photon if it has leaked from the scene */
66 + #ifdef DEBUG_PMAP
67 +         fprintf(stderr,
68 +                 "Volume photon leaked from scene at [%.3f %.3f %.3f]\n",
69 +                 ray -> rop [0], ray -> rop [1], ray -> rop [2]);
70 + #endif                
71 +         return 0;
72 +      }
73 +        
74        setcolor(cvext, exp(-ray -> rmax * ray -> cext [0]),
75                        exp(-ray -> rmax * ray -> cext [1]),
76                        exp(-ray -> rmax * ray -> cext [2]));
# Line 143 | Line 80 | static int photonParticipate (RAY *ray)
80        colorNorm(ray -> rcol);
81        VCOPY(ray -> rorg, ray -> rop);
82        
83 + #if 0
84 +      if (albedo > FTINY && ray -> rlvl > 0)
85 + #else
86 +      /* Store volume photons unconditionally in mist to also account for
87 +         direct inscattering from sources */
88        if (albedo > FTINY)
89 + #endif
90           /* Add to volume photon map */
91 <         if (ray -> rlvl > 0) addPhoton(volumePmap, ray);
91 >         newPhoton(volumePmap, ray);
92          
93        /* Absorbed? */
94 <      if (pmapRandom(rouletteState) > albedo) return 0;
94 >      if (pmapRandom(rouletteState) > albedo)
95 >         return 0;
96        
97        /* Colour bleeding without attenuation (?) */
98        multcolor(ray -> rcol, ray -> albedo);
99        scalecolor(ray -> rcol, 1 / albedo);    
100        
101        /* Scatter photon */
102 <      cosTheta = ray -> gecc <= FTINY ? 2 * pmapRandom(scatterState) - 1
103 <                                      : 1 / (2 * ray -> gecc) *
104 <                                            (1 + ray -> gecc * ray -> gecc -
105 <                                               (1 - ray -> gecc * ray -> gecc) /
106 <                                               (1 - ray -> gecc + 2 * ray -> gecc *
107 <                                                  pmapRandom(scatterState)));
108 <                                                  
109 <      cosPhi = cos(2 * PI * pmapRandom(scatterState));
110 <      du = dv = sqrt(1 - cosTheta * cosTheta);   /* sin(theta) */
111 <      du *= cosPhi;
112 <      dv *= sqrt(1 - cosPhi * cosPhi);           /* sin(phi) */
102 >      xi1 = pmapRandom(scatterState);
103 >      cosTheta = ray -> gecc <= FTINY
104 >                    ? 2 * xi1 - 1
105 >                    : 0.5 / gecc *
106 >                      (1 + gecc2 - sqr((1 - gecc2) /
107 >                                       (1 + gecc * (2 * xi1 - 1))));
108 >
109 >      phi = 2 * PI * pmapRandom(scatterState);
110 >      du = dv = sqrt(1 - sqr(cosTheta));     /* sin(theta) */
111 >      du *= cos(phi);
112 >      dv *= sin(phi);
113        
114        /* Get axes u & v perpendicular to photon direction */
115        i = 0;
# Line 179 | Line 123 | static int photonParticipate (RAY *ray)
123        for (i = 0; i < 3; i++)
124           ray -> rdir [i] = du * u [i] + dv * v [i] +
125                             cosTheta * ray -> rdir [i];
126 +
127        ray -> rlvl++;
128        ray -> rmax = -log(pmapRandom(mediumState)) / cext;
129     }  
130 <    
130 >  
131 >   /* Passed through medium until intersecting local object */  
132     setcolor(cvext, exp(-ray -> rot * ray -> cext [0]),
133                     exp(-ray -> rot * ray -> cext [1]),
134                     exp(-ray -> rot * ray -> cext [2]));
135                    
136     /* Modify ray color and normalise */
137     multcolor(ray -> rcol, cvext);
138 <   colorNorm(ray -> rcol);
139 <  
194 <   /* Passed through medium */  
138 >   colorNorm(ray -> rcol);  
139 >
140     return 1;
141   }
142  
# Line 201 | Line 146 | void tracePhoton (RAY *ray)
146   /* Follow photon as it bounces around... */
147   {
148     long mod;
149 <   OBJREC* mat;
149 >   OBJREC *mat, *port = NULL;
150 >  
151 >   if (!ray -> parent) {
152 >      /* !!!  PHOTON PORT REJECTION SAMPLING HACK: get photon port for
153 >       * !!!  primary ray from ray -> ro, then reset the latter to NULL so
154 >       * !!!  as not to interfere with localhit() */
155 >      port = ray -> ro;
156 >      ray -> ro = NULL;
157 >   }
158  
159     if (ray -> rlvl > photonMaxBounce) {
160 + #ifdef PMAP_RUNAWAY_WARN  
161        error(WARNING, "runaway photon!");
162 + #endif      
163        return;
164     }
165 <  
165 >  
166     if (colorAvg(ray -> cext) > FTINY && !photonParticipate(ray))
167        return;
168 <      
168 >
169     if (localhit(ray, &thescene)) {
170        mod = ray -> ro -> omod;
171 <      
171 >
172 >      if (port && ray -> ro != port) {
173 >         /* !!! PHOTON PORT REJECTION SAMPLING HACK !!!
174 >          * Terminate photon if emitted from port without intersecting it;
175 >          * this can happen when the port's partitions extend beyond its
176 >          * actual geometry, e.g.  with polygons.  Since the total flux
177 >          * relayed by the port is based on the (in this case) larger
178 >          * partition area, it is overestimated; terminating these photons
179 >          * constitutes rejection sampling and thereby compensates any bias
180 >          * incurred by the overestimated flux.  */
181 > #ifdef PMAP_PORTREJECT_WARN
182 >         sprintf(errmsg, "photon outside port %s", ray -> ro -> oname);
183 >         error(WARNING, errmsg);
184 > #endif        
185 >         return;
186 >      }
187 >
188        if ((ray -> clipset && inset(ray -> clipset, mod)) || mod == OVOID) {
189           /* Transfer ray if modifier is VOID or clipped within antimatta */
190           RAY tray;
# Line 231 | Line 202 | void tracePhoton (RAY *ray)
202  
203  
204   static void preComputeGlobal (PhotonMap *pmap)
205 < /* Precompute irradiance from global photons for final gathering using
206 <   the first finalGather * pmap -> heapSize photons in the heap. Returns
207 <   new heap with precomputed photons. */
205 > /* Precompute irradiance from global photons for final gathering for  
206 >   a random subset of finalGather * pmap -> numPhotons photons, and builds
207 >   the photon map, discarding the original photons. */
208 > /* !!! NOTE: PRECOMPUTATION WITH OOC CURRENTLY WITHOUT CACHE !!! */  
209   {
210 <   unsigned long i, nuHeapSize;
211 <   unsigned j;
212 <   Photon *nuHeap, *p;
213 <   COLOR irrad;
214 <   RAY ray;
215 <   float nuMinPos [3], nuMaxPos [3];
210 >   unsigned long  i, numPreComp;
211 >   unsigned       j;
212 >   PhotonIdx      pIdx;
213 >   Photon         photon;
214 >   RAY            ray;
215 >   PhotonMap      nuPmap;
216  
217 <   repComplete = nuHeapSize = finalGather * pmap -> heapSize;
217 >   repComplete = numPreComp = finalGather * pmap -> numPhotons;
218    
219 <   if (photonRepTime) {
219 >   if (verbose) {
220        sprintf(errmsg,
221 <              "Precomputing irradiance for %ld global photons...\n",
222 <              nuHeapSize);
221 >              "\nPrecomputing irradiance for %ld global photons\n",
222 >              numPreComp);
223        eputs(errmsg);
224 + #if NIX      
225        fflush(stderr);
226 + #endif      
227     }
228    
229 <   p = nuHeap = (Photon*)malloc(nuHeapSize * sizeof(Photon));
230 <   if (!nuHeap)
231 <      error(USER, "can't allocate photon heap");
232 <      
233 <   for (j = 0; j <= 2; j++) {
234 <      nuMinPos [j] = FHUGE;
261 <      nuMaxPos [j] = -FHUGE;
262 <   }
229 >   /* Copy photon map for precomputed photons */
230 >   memcpy(&nuPmap, pmap, sizeof(PhotonMap));
231 >
232 >   /* Zero counters, init new heap and extents */  
233 >   nuPmap.numPhotons = 0;  
234 >   initPhotonHeap(&nuPmap);
235    
236 +   for (j = 0; j < 3; j++) {  
237 +      nuPmap.minPos [j] = FHUGE;
238 +      nuPmap.maxPos [j] = -FHUGE;
239 +   }
240 +
241     /* Record start time, baby */
242     repStartTime = time(NULL);
243 + #ifdef SIGCONT
244     signal(SIGCONT, pmapPreCompReport);
245 + #endif
246     repProgress = 0;
268 –   memcpy(nuHeap, pmap -> heap, nuHeapSize * sizeof(Photon));
247    
248 <   for (i = 0, p = nuHeap; i < nuHeapSize; i++, p++) {
249 <      ray.ro = NULL;
250 <      VCOPY(ray.rop, p -> pos);
248 >   photonRay(NULL, &ray, PRIMARY, NULL);
249 >   ray.ro = NULL;
250 >  
251 >   for (i = 0; i < numPreComp; i++) {
252 >      /* Get random photon from stratified distribution in source heap to
253 >       * avoid duplicates and clustering */
254 >      pIdx = firstPhoton(pmap) +
255 >             (unsigned long)((i + pmapRandom(pmap -> randState)) /
256 >                             finalGather);
257 >      getPhoton(pmap, pIdx, &photon);
258        
259 <      /* Update min and max positions & set ray normal */
260 <      for (j = 0; j < 3; j++) {
261 <         if (p -> pos [j] < nuMinPos [j]) nuMinPos [j] = p -> pos [j];
262 <         if (p -> pos [j] > nuMaxPos [j]) nuMaxPos [j] = p -> pos [j];
278 <         ray.ron [j] = p -> norm [j] / 127.0;
279 <      }
259 >      /* Init dummy photon ray with intersection at photon position */
260 >      VCOPY(ray.rop, photon.pos);
261 >      for (j = 0; j < 3; j++)
262 >         ray.ron [j] = photon.norm [j] / 127.0;
263        
264 <      photonDensity(pmap, &ray, irrad);
265 <      setPhotonFlux(p, irrad);
264 >      /* Get density estimate at photon position */
265 >      photonDensity(pmap, &ray, ray.rcol);
266 >                  
267 >      /* Append photon to new heap from ray */
268 >      newPhoton(&nuPmap, &ray);
269 >      
270 >      /* Update progress */
271        repProgress++;
272        
273        if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
274           pmapPreCompReport();
275 < #ifndef BSD
275 > #ifdef SIGCONT
276        else signal(SIGCONT, pmapPreCompReport);
277   #endif
278     }
279    
280 +   /* Flush heap */
281 +   flushPhotonHeap(&nuPmap);
282 +  
283 + #ifdef SIGCONT  
284     signal(SIGCONT, SIG_DFL);
285 + #endif
286    
287 <   /* Replace & rebuild heap */
288 <   free(pmap -> heap);
289 <   pmap -> heap = nuHeap;
297 <   pmap -> heapSize = pmap -> heapEnd = nuHeapSize;
298 <   VCOPY(pmap -> minPos, nuMinPos);
299 <   VCOPY(pmap -> maxPos, nuMaxPos);
287 >   /* Trash original pmap, replace with precomputed one */
288 >   deletePhotons(pmap);
289 >   memcpy(pmap, &nuPmap, sizeof(PhotonMap));
290    
291 <   if (photonRepTime) {
292 <      eputs("Rebuilding global photon heap...\n");
291 >   if (verbose) {
292 >      eputs("\nRebuilding precomputed photon map\n");
293 > #if NIX      
294        fflush(stderr);
295 + #endif      
296     }
297 <  
298 <   balancePhotons(pmap, NULL);
297 >
298 >   /* Rebuild underlying data structure, destroying heap */  
299 >   buildPhotonMap(pmap, NULL, NULL, 1);
300   }
301  
302  
303  
304 < void distribPhotons (PhotonMap **pmaps)
304 > typedef struct {
305 >   unsigned long  numPhotons [NUM_PMAP_TYPES],
306 >                  numEmitted, numComplete;
307 > } PhotonCnt;
308 >
309 >
310 >
311 > void distribPhotons (PhotonMap **pmaps, unsigned numProc)
312   {
313 <   EmissionMap emap;
314 <   char errmsg2 [128];
315 <   unsigned t, srcIdx, passCnt = 0, prePassCnt = 0;
316 <   double totalFlux = 0;
317 <   PhotonMap *pm;
313 >   EmissionMap    emap;
314 >   char           errmsg2 [128], shmFname [PMAP_TMPFNLEN];
315 >   unsigned       t, srcIdx, proc;
316 >   double         totalFlux = 0;
317 >   int            shmFile, stat, pid;
318 >   PhotonMap      *pm;
319 >   PhotonCnt      *photonCnt;
320    
321 <   for (t = 0; t < NUM_PMAP_TYPES && !photonMaps [t]; t++);
321 >   for (t = 0; t < NUM_PMAP_TYPES && !pmaps [t]; t++);
322 >  
323     if (t >= NUM_PMAP_TYPES)
324 <      error(USER, "no photon maps defined");
324 >      error(USER, "no photon maps defined in distribPhotons");
325        
326     if (!nsources)
327 <      error(USER, "no light sources");
327 >      error(USER, "no light sources in distribPhotons");
328  
329     /* ===================================================================
330      * INITIALISATION - Set up emission and scattering funcs
# Line 330 | Line 333 | void distribPhotons (PhotonMap **pmaps)
333     emap.maxPartitions = MAXSPART;
334     emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1);
335     if (!emap.partitions)
336 <      error(INTERNAL, "can't allocate source partitions");
336 >      error(INTERNAL, "can't allocate source partitions in distribPhotons");
337        
338     /* Initialise all defined photon maps */
339     for (t = 0; t < NUM_PMAP_TYPES; t++)
340 <      initPhotonMap(photonMaps [t], t);
340 >      if (pmaps [t]) {
341 >         initPhotonMap(pmaps [t], t);
342 >         /* Open photon heapfile */
343 >         initPhotonHeap(pmaps [t]);
344 >         /* Per-subprocess target count */
345 >         pmaps [t] -> distribTarget /= numProc;
346 >        
347 >         if (!pmaps [t] -> distribTarget)
348 >            error(INTERNAL, "no photons to distribute in distribPhotons");
349 >      }
350  
351     initPhotonEmissionFuncs();
352     initPhotonScatterFuncs();
353    
354 <   /* Get photon ports if specified */
355 <   if (ambincl == 1)
344 <      getPhotonPorts();
354 >   /* Get photon ports from modifier list */
355 >   getPhotonPorts(photonPortList);
356  
357     /* Get photon sensor modifiers */
358     getPhotonSensors(photonSensorList);
359    
360 <   /* Seed RNGs for photon distribution */
361 <   pmapSeed(randSeed, partState);
362 <   pmapSeed(randSeed, emitState);
363 <   pmapSeed(randSeed, cntState);
364 <   pmapSeed(randSeed, mediumState);
365 <   pmapSeed(randSeed, scatterState);
366 <   pmapSeed(randSeed, rouletteState);
360 > #if NIX
361 >   /* Set up shared mem for photon counters (zeroed by ftruncate) */
362 >   strcpy(shmFname, PMAP_TMPFNAME);
363 >   shmFile = mkstemp(shmFname);
364 >
365 >   if (shmFile < 0 || ftruncate(shmFile, sizeof(*photonCnt)) < 0)
366 >      error(SYSTEM, "failed shared mem init in distribPhotons");
367 >
368 >   photonCnt = mmap(NULL, sizeof(*photonCnt), PROT_READ | PROT_WRITE,
369 >                    MAP_SHARED, shmFile, 0);
370 >                    
371 >   if (photonCnt == MAP_FAILED)
372 >      error(SYSTEM, "failed mapping shared memory in distribPhotons");
373 > #else
374 >   /* Allocate photon counters statically on Windoze */
375 >   if (!(photonCnt = malloc(sizeof(PhotonCnt))))
376 >      error(SYSTEM, "failed trivial malloc in distribPhotons");
377 >   photonCnt -> numEmitted = photonCnt -> numComplete = 0;      
378 > #endif /* NIX */
379 >
380 >   if (verbose) {
381 >      sprintf(errmsg, "\nIntegrating flux from %d sources", nsources);
382 >      
383 >      if (photonPorts) {
384 >         sprintf(errmsg2, " via %d ports", numPhotonPorts);
385 >         strcat(errmsg, errmsg2);
386 >      }
387 >      
388 >      strcat(errmsg, "\n");
389 >      eputs(errmsg);
390 >   }  
391    
357 –   if (photonRepTime)
358 –      eputs("\n");
359 –  
392     /* ===================================================================
393      * FLUX INTEGRATION - Get total photon flux from light sources
394      * =================================================================== */
395 <   for (srcIdx = 0; srcIdx < nsources; srcIdx++) {        
395 >   for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
396        unsigned portCnt = 0;
397        emap.src = source + srcIdx;
398        
399 <      do {
399 >      do {  /* Need at least one iteration if no ports! */
400           emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt
401                                                     : NULL;
402           photonPartition [emap.src -> so -> otype] (&emap);
403          
404 <         if (photonRepTime) {
405 <            sprintf(errmsg, "Integrating flux from source %s ",
404 >         if (verbose) {
405 >            sprintf(errmsg, "\tIntegrating flux from source %s ",
406                      source [srcIdx].so -> oname);
407 <                    
407 >
408              if (emap.port) {
409                 sprintf(errmsg2, "via port %s ",
410                         photonPorts [portCnt].so -> oname);
411                 strcat(errmsg, errmsg2);
412              }
413 <            
414 <            sprintf(errmsg2, "(%lu partitions)...\n", emap.numPartitions);
413 >
414 >            sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions);
415              strcat(errmsg, errmsg2);
416              eputs(errmsg);
417 + #if NIX            
418              fflush(stderr);
419 + #endif            
420           }
421          
422           for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
# Line 397 | Line 431 | void distribPhotons (PhotonMap **pmaps)
431  
432     if (totalFlux < FTINY)
433        error(USER, "zero flux from light sources");
400 –
401 –   /* Record start time and enable progress report signal handler */
402 –   repStartTime = time(NULL);
403 –   signal(SIGCONT, pmapDistribReport);
404 –   repProgress = prePassCnt = 0;
405 –  
406 –   if (photonRepTime)
407 –      eputs("\n");
408 –  
409 –   /* ===================================================================
410 –    * 2-PASS PHOTON DISTRIBUTION
411 –    * Pass 1 (pre):  emit fraction of target photon count
412 –    * Pass 2 (main): based on outcome of pass 1, estimate remaining number
413 –    *                of photons to emit to approximate target count
414 –    * =================================================================== */
415 –   do {
416 –      double numEmit;
434        
435 <      if (!passCnt) {  
436 <         /* INIT PASS 1 */
420 <         /* Skip if no photons contributed after sufficient iterations; make
421 <          * it clear to user which photon maps are missing so (s)he can
422 <          * check the scene geometry and materials */
423 <         if (++prePassCnt > maxPreDistrib) {
424 <            sprintf(errmsg, "too many prepasses");
435 >   /* Record start time for progress reports */
436 >   repStartTime = time(NULL);
437  
438 <            for (t = 0; t < NUM_PMAP_TYPES; t++)
439 <               if (photonMaps [t] && !photonMaps [t] -> heapEnd) {
440 <                  sprintf(errmsg2, ", no %s photons stored", pmapName [t]);
441 <                  strcat(errmsg, errmsg2);
430 <               }
431 <            
432 <            error(USER, errmsg);
433 <            break;
434 <         }
438 >   if (verbose) {
439 >      sprintf(errmsg, "\nPhoton distribution @ %d procs\n", numProc);
440 >      eputs(errmsg);
441 >   }
442  
443 <         /* Num to emit is fraction of minimum target count */
444 <         numEmit = FHUGE;
443 >   /* MAIN LOOP */  
444 >   for (proc = 0; proc < numProc; proc++) {
445 > #if NIX          
446 >      if (!(pid = fork())) {
447 >         /* SUBPROCESS ENTERS HERE; open and mmapped files inherited */
448 > #else
449 >      if (1) {
450 >         /* No subprocess under Windoze */
451 > #endif
452 >         /* Local photon counters for this subprocess */
453 >         unsigned       passCnt = 0, prePassCnt = 0;
454 >         unsigned long  lastNumPhotons [NUM_PMAP_TYPES];
455 >         unsigned long  localNumEmitted = 0; /* Num photons emitted by this
456 >                                                subprocess alone */
457          
458 <         for (t = 0; t < NUM_PMAP_TYPES; t++)
459 <            if (photonMaps [t])
460 <               numEmit = min(photonMaps [t] -> distribTarget, numEmit);
458 >         /* Seed RNGs from PID for decorellated photon distribution */
459 >         pmapSeed(randSeed + proc, partState);
460 >         pmapSeed(randSeed + (proc + 1) % numProc, emitState);
461 >         pmapSeed(randSeed + (proc + 2) % numProc, cntState);
462 >         pmapSeed(randSeed + (proc + 3) % numProc, mediumState);
463 >         pmapSeed(randSeed + (proc + 4) % numProc, scatterState);
464 >         pmapSeed(randSeed + (proc + 5) % numProc, rouletteState);
465                
466 <         numEmit *= preDistrib;
467 <      }
466 > #ifdef DEBUG_PMAP          
467 >         /* Output child process PID after random delay to prevent corrupted
468 >          * console output due to race condition */
469 >         usleep(1e6 * pmapRandom(rouletteState));
470 >         fprintf(stderr, "Proc %d: PID = %d "
471 >                 "(waiting 10 sec to attach debugger...)\n",
472 >                 proc, getpid());
473 >         /* Allow time for debugger to attach to child process */
474 >         sleep(10);
475 > #endif            
476  
446 –      else {            
447 –         /* INIT PASS 2 */
448 –         /* Based on the outcome of the predistribution we can now estimate
449 –          * how many more photons we have to emit for each photon map to
450 –          * meet its respective target count. This value is clamped to 0 in
451 –          * case the target has already been exceeded in the pass 1. Note
452 –          * repProgress is the number of photons emitted thus far, while
453 –          * heapEnd is the number of photons stored in each photon map. */
454 –         double maxDistribRatio = 0;
455 –
456 –         /* Set the distribution ratio for each map; this indicates how many
457 –          * photons of each respective type are stored per emitted photon,
458 –          * and is used as probability for storing a photon by addPhoton().
459 –          * Since this biases the photon density, addPhoton() promotes the
460 –          * flux of stored photons to compensate. */
477           for (t = 0; t < NUM_PMAP_TYPES; t++)
478 <            if ((pm = photonMaps [t])) {
479 <               pm -> distribRatio = (double)pm -> distribTarget /
480 <                                    pm -> heapEnd - 1;
481 <
482 <               /* Check if photon map "overflowed", i.e. exceeded its target
483 <                * count in the prepass; correcting the photon flux via the
484 <                * distribution ratio is no longer possible, as no more
485 <                * photons of this type will be stored, so notify the user
486 <                * rather than deliver incorrect results.
471 <                * In future we should handle this more intelligently by
472 <                * using the photonFlux in each photon map to individually
473 <                * correct the flux after distribution. */
474 <               if (pm -> distribRatio <= FTINY) {
475 <                  sprintf(errmsg,
476 <                          "%s photon map overflow in prepass, reduce -apD",
477 <                          pmapName [t]);
478 <                  error(INTERNAL, errmsg);
479 <               }
480 <                  
481 <               maxDistribRatio = max(pm -> distribRatio, maxDistribRatio);
482 <            }
483 <        
484 <         /* Normalise distribution ratios and calculate number of photons to
485 <          * emit in main pass */
486 <         for (t = 0; t < NUM_PMAP_TYPES; t++)
487 <            if ((pm = photonMaps [t]))
488 <               pm -> distribRatio /= maxDistribRatio;
489 <              
490 <         if ((numEmit = repProgress * maxDistribRatio) < FTINY)
491 <            /* No photons left to distribute in main pass */
492 <            break;
493 <      }
494 <      
495 <      /* Set completion count for progress report */
496 <      repComplete = numEmit + repProgress;                            
497 <      
498 <      /* PHOTON DISTRIBUTION LOOP */
499 <      for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
500 <         unsigned portCnt = 0;
501 <         emap.src = source + srcIdx;
502 <                  
478 >            lastNumPhotons [t] = 0;
479 >            
480 >         /* =============================================================
481 >          * 2-PASS PHOTON DISTRIBUTION
482 >          * Pass 1 (pre):  emit fraction of target photon count
483 >          * Pass 2 (main): based on outcome of pass 1, estimate remaining
484 >          *                number of photons to emit to approximate target
485 >          *                count
486 >          * ============================================================= */
487           do {
488 <            emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt
505 <                                                      : NULL;
506 <            photonPartition [emap.src -> so -> otype] (&emap);
488 >            double numEmit;
489              
490 <            if (photonRepTime) {
491 <               if (!passCnt)
492 <                  sprintf(errmsg, "PREPASS %d on source %s ",
493 <                          prePassCnt, source [srcIdx].so -> oname);
494 <               else
495 <                  sprintf(errmsg, "MAIN PASS on source %s ",
496 <                          source [srcIdx].so -> oname);
497 <                      
498 <               if (emap.port) {
499 <                  sprintf(errmsg2, "via port %s ",
500 <                          photonPorts [portCnt].so -> oname);
501 <                  strcat(errmsg, errmsg2);
490 >            if (!passCnt) {  
491 >               /* INIT PASS 1 */              
492 >               /* Skip if no photons contributed after sufficient
493 >                * iterations; make it clear to user which photon maps are
494 >                * missing so (s)he can check geometry and materials */
495 >               if (++prePassCnt > maxPreDistrib) {
496 >                  sprintf(errmsg, "proc %d: too many prepasses", proc);
497 >
498 >                  for (t = 0; t < NUM_PMAP_TYPES; t++)
499 >                     if (pmaps [t] && !pmaps [t] -> numPhotons) {
500 >                        sprintf(errmsg2, ", no %s photons stored",
501 >                                pmapName [t]);
502 >                        strcat(errmsg, errmsg2);
503 >                     }
504 >                  
505 >                  error(USER, errmsg);
506 >                  break;
507                 }
508 +
509 +               /* Num to emit is fraction of minimum target count */
510 +               numEmit = FHUGE;
511                
512 <               sprintf(errmsg2, "(%lu partitions)...\n", emap.numPartitions);
513 <               strcat(errmsg, errmsg2);
514 <               eputs(errmsg);
515 <               fflush(stderr);
512 >               for (t = 0; t < NUM_PMAP_TYPES; t++)
513 >                  if (pmaps [t])
514 >                     numEmit = min(pmaps [t] -> distribTarget, numEmit);
515 >                    
516 >               numEmit *= preDistrib;
517              }
518 <            
519 <            for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
520 <                 emap.partitionCnt++) {            
521 <               double partNumEmit;
522 <               unsigned long partEmitCnt;
518 >            else {            
519 >               /* INIT PASS 2 */
520 >               /* Based on the outcome of the predistribution we can now
521 >                * estimate how many more photons we have to emit for each
522 >                * photon map to meet its respective target count.  This
523 >                * value is clamped to 0 in case the target has already been
524 >                * exceeded in the pass 1. */
525 >               double maxDistribRatio = 0;
526 >
527 >               /* Set the distribution ratio for each map; this indicates
528 >                * how many photons of each respective type are stored per
529 >                * emitted photon, and is used as probability for storing a
530 >                * photon by newPhoton().  Since this biases the photon
531 >                * density, newPhoton() promotes the flux of stored photons
532 >                * to compensate.  */
533 >               for (t = 0; t < NUM_PMAP_TYPES; t++)
534 >                  if ((pm = pmaps [t])) {
535 >                     pm -> distribRatio = (double)pm -> distribTarget /
536 >                                          pm -> numPhotons - 1;
537 >
538 >                     /* Check if photon map "overflowed", i.e. exceeded its
539 >                      * target count in the prepass; correcting the photon
540 >                      * flux via the distribution ratio is no longer
541 >                      * possible, as no more photons of this type will be
542 >                      * stored, so notify the user rather than deliver
543 >                      * incorrect results.  In future we should handle this
544 >                      * more intelligently by using the photonFlux in each
545 >                      * photon map to individually correct the flux after
546 >                      * distribution.  */
547 >                     if (pm -> distribRatio <= FTINY) {
548 >                        sprintf(errmsg, "%s photon map overflow in "
549 >                                "prepass, reduce -apD", pmapName [t]);
550 >                        error(INTERNAL, errmsg);
551 >                     }
552 >                        
553 >                     maxDistribRatio = max(pm -> distribRatio,
554 >                                           maxDistribRatio);
555 >                  }
556                
557 <               /* Get photon origin within current source partishunn and
558 <                * build emission map */
559 <               photonOrigin [emap.src -> so -> otype] (&emap);
560 <               initPhotonEmission(&emap, pdfSamples);
561 <              
562 <               /* Number of photons to emit from ziss partishunn --
563 <                * proportional to flux; photon ray weight and scalar flux
564 <                * are uniform (the latter only varying in RGB). */
565 <               partNumEmit = numEmit * colorAvg(emap.partFlux) / totalFlux;
566 <               partEmitCnt = (unsigned long)partNumEmit;
543 <              
544 <               /* Probabilistically account for fractional photons */
545 <               if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
546 <                  partEmitCnt++;
557 >               /* Normalise distribution ratios and calculate number of
558 >                * photons to emit in main pass */
559 >               for (t = 0; t < NUM_PMAP_TYPES; t++)
560 >                  if ((pm = pmaps [t]))
561 >                     pm -> distribRatio /= maxDistribRatio;
562 >                    
563 >               if ((numEmit = localNumEmitted * maxDistribRatio) < FTINY)
564 >                  /* No photons left to distribute in main pass */
565 >                  break;
566 >            }
567  
568 <               /* Integer counter avoids FP rounding errors */
569 <               while (partEmitCnt--) {
570 <                  RAY photonRay;
568 >            /* Update shared completion counter for progreport by parent */
569 >            photonCnt -> numComplete += numEmit;                            
570 >
571 >            /* PHOTON DISTRIBUTION LOOP */
572 >            for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
573 >               unsigned portCnt = 0;
574 >               emap.src = source + srcIdx;
575 >
576 >               do {  /* Need at least one iteration if no ports! */
577 >                  emap.port = emap.src -> sflags & SDISTANT
578 >                              ? photonPorts + portCnt : NULL;
579 >                  photonPartition [emap.src -> so -> otype] (&emap);
580 >
581 >                  if (verbose && !proc) {
582 >                     /* Output from subproc 0 only to avoid race condition
583 >                      * on console I/O */
584 >                     if (!passCnt)
585 >                        sprintf(errmsg, "\tPREPASS %d on source %s ",
586 >                                prePassCnt, source [srcIdx].so -> oname);
587 >                     else
588 >                        sprintf(errmsg, "\tMAIN PASS on source %s ",
589 >                                source [srcIdx].so -> oname);
590 >
591 >                     if (emap.port) {
592 >                        sprintf(errmsg2, "via port %s ",
593 >                                photonPorts [portCnt].so -> oname);
594 >                        strcat(errmsg, errmsg2);
595 >                     }
596 >
597 >                     sprintf(errmsg2, "(%lu partitions)\n",
598 >                             emap.numPartitions);
599 >                     strcat(errmsg, errmsg2);
600 >                     eputs(errmsg);
601 > #if NIX                    
602 >                     fflush(stderr);
603 > #endif                    
604 >                  }
605                    
606 <                  /* Emit photon based on PDF and trace through scene until
607 <                   * absorbed/leaked */
608 <                  emitPhoton(&emap, &photonRay);
609 <                  tracePhoton(&photonRay);
606 >                  for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
607 >                       emap.partitionCnt++) {
608 >                     double partNumEmit;
609 >                     unsigned long partEmitCnt;
610 >                    
611 >                     /* Get photon origin within current source partishunn
612 >                      * and build emission map */
613 >                     photonOrigin [emap.src -> so -> otype] (&emap);
614 >                     initPhotonEmission(&emap, pdfSamples);
615 >                    
616 >                     /* Number of photons to emit from ziss partishunn --
617 >                      * proportional to flux; photon ray weight and scalar
618 >                      * flux are uniform (latter only varying in RGB). */
619 >                     partNumEmit = numEmit * colorAvg(emap.partFlux) /
620 >                                   totalFlux;
621 >                     partEmitCnt = (unsigned long)partNumEmit;
622 >                    
623 >                     /* Probabilistically account for fractional photons */
624 >                     if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
625 >                        partEmitCnt++;
626 >
627 >                     /* Update local and shared (global) emission counter */
628 >                     photonCnt -> numEmitted += partEmitCnt;
629 >                     localNumEmitted += partEmitCnt;
630 >
631 >                     /* Integer counter avoids FP rounding errors during
632 >                      * iteration */
633 >                     while (partEmitCnt--) {
634 >                        RAY photonRay;
635 >                        
636 >                        /* Emit photon based on PDF and trace through scene
637 >                         * until absorbed/leaked */
638 >                        emitPhoton(&emap, &photonRay);
639 > #if 1
640 >                        if (emap.port)
641 >                           /* !!!  PHOTON PORT REJECTION SAMPLING HACK: set
642 >                            * !!!  photon port as fake hit object for
643 >                            * !!!  primary ray to check for intersection in
644 >                            * !!!  tracePhoton() */
645 >                           photonRay.ro = emap.port -> so;
646 > #endif
647 >                        tracePhoton(&photonRay);
648 >                     }                                          
649 >
650 >                     /* Update shared global photon count for each pmap */
651 >                     for (t = 0; t < NUM_PMAP_TYPES; t++)
652 >                        if (pmaps [t]) {
653 >                           photonCnt -> numPhotons [t] +=
654 >                              pmaps [t] -> numPhotons - lastNumPhotons [t];
655 >                           lastNumPhotons [t] = pmaps [t] -> numPhotons;
656 >                        }
657 > #if !NIX
658 >                     /* Synchronous progress report on Windoze */
659 >                     if (!proc && photonRepTime > 0 &&
660 >                           time(NULL) >= repLastTime + photonRepTime) {
661 >                        repEmitted = repProgress = photonCnt -> numEmitted;
662 >                        repComplete = photonCnt -> numComplete;                          
663 >                        pmapDistribReport();
664 >                     }
665 > #endif
666 >                  }
667                    
668 <                  /* Record progress */
669 <                  repProgress++;
670 <                  
671 <                  if (photonRepTime > 0 &&
672 <                      time(NULL) >= repLastTime + photonRepTime)
673 <                     pmapDistribReport();
674 <                  #ifndef BSD
675 <                     else signal(SIGCONT, pmapDistribReport);
676 <                  #endif
668 >                  portCnt++;
669 >               } while (portCnt < numPhotonPorts);
670 >            }
671 >            
672 >            for (t = 0; t < NUM_PMAP_TYPES; t++)
673 >               if (pmaps [t] && !pmaps [t] -> numPhotons) {
674 >                  /* Double preDistrib in case a photon map is empty and
675 >                   * redo pass 1 --> possibility of infinite loop for
676 >                   * pathological scenes (e.g.  absorbing materials) */
677 >                  preDistrib *= 2;
678 >                  break;
679                 }
680 +            
681 +            if (t >= NUM_PMAP_TYPES)
682 +               /* No empty photon maps found; now do pass 2 */
683 +               passCnt++;
684 +         } while (passCnt < 2);
685 +        
686 +         /* Flush heap buffa for every pmap one final time;
687 +          * avoids potential data corruption! */
688 +         for (t = 0; t < NUM_PMAP_TYPES; t++)
689 +            if (pmaps [t]) {
690 +               flushPhotonHeap(pmaps [t]);
691 +               /* Heap file closed automatically on exit
692 +                  fclose(pmaps [t] -> heap); */
693 + #ifdef DEBUG_PMAP              
694 +               sprintf(errmsg, "Proc %d: total %ld photons\n", proc,
695 +                       pmaps [t] -> numPhotons);
696 +               eputs(errmsg);
697 + #endif              
698              }
699 <                        
700 <            portCnt++;
701 <         } while (portCnt < numPhotonPorts);
699 > #if NIX
700 >         /* Terminate subprocess */
701 >         exit(0);
702 > #endif
703        }
704 +      else if (pid < 0)
705 +         error(SYSTEM, "failed to fork subprocess in distribPhotons");        
706 +   }
707 +
708 + #if NIX
709 +   /* PARENT PROCESS CONTINUES HERE */
710 + #ifdef SIGCONT
711 +   /* Enable progress report signal handler */
712 +   signal(SIGCONT, pmapDistribReport);
713 + #endif  
714 +   /* Wait for subprocesses complete while reporting progress */
715 +   proc = numProc;
716 +   while (proc) {
717 +      while (waitpid(-1, &stat, WNOHANG) > 0) {
718 +         /* Subprocess exited; check status */
719 +         if (!WIFEXITED(stat) || WEXITSTATUS(stat))
720 +            error(USER, "failed photon distribution");
721        
722 +         --proc;
723 +      }
724 +      
725 +      /* Nod off for a bit and update progress  */
726 +      sleep(1);
727 +
728 +      /* Asynchronous progress report from shared subprocess counters */  
729 +      repEmitted = repProgress = photonCnt -> numEmitted;
730 +      repComplete = photonCnt -> numComplete;      
731 +
732 +      repProgress = repComplete = 0;
733        for (t = 0; t < NUM_PMAP_TYPES; t++)
734 <         if (photonMaps [t] && !photonMaps [t] -> heapEnd) {
735 <            /* Double preDistrib in case a photon map is empty and redo
736 <             * pass 1 --> possibility of infinite loop for pathological
737 <             * scenes (e.g. absorbing materials) */
578 <            preDistrib *= 2;
579 <            break;
734 >         if ((pm = pmaps [t])) {
735 >            /* Get global photon count from shmem updated by subprocs */
736 >            repProgress += pm -> numPhotons = photonCnt -> numPhotons [t];
737 >            repComplete += pm -> distribTarget;
738           }
739 <      
582 <      if (t >= NUM_PMAP_TYPES) {
583 <         /* No empty photon maps found; now do pass 2 */
584 <         passCnt++;
585 <         if (photonRepTime)
586 <            eputs("\n");
587 <      }
588 <   } while (passCnt < 2);
739 >      repComplete *= numProc;
740  
741 +      if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
742 +         pmapDistribReport();
743 + #ifdef SIGCONT
744 +      else signal(SIGCONT, pmapDistribReport);
745 + #endif
746 +   }
747 + #endif /* NIX */
748 +
749     /* ===================================================================
750 <    * POST-DISTRIBUTION - Set photon flux and build kd-tree, etc.
750 >    * POST-DISTRIBUTION - Set photon flux and build data struct for photon
751 >    * storage, etc.
752      * =================================================================== */
753 + #ifdef SIGCONT    
754 +   /* Reset signal handler */
755     signal(SIGCONT, SIG_DFL);
756 + #endif
757     free(emap.samples);
758    
759 <   /* Set photon flux (repProgress is total num emitted) */
760 <   totalFlux /= repProgress;
761 <  
759 >   /* Set photon flux */
760 >   totalFlux /= photonCnt -> numEmitted;
761 > #if NIX  
762 >   /* Photon counters no longer needed, unmap shared memory */
763 >   munmap(photonCnt, sizeof(*photonCnt));
764 >   close(shmFile);
765 >   unlink(shmFname);
766 > #else
767 >   free(photonCnt);  
768 > #endif      
769 >   if (verbose)
770 >      eputs("\n");
771 >      
772     for (t = 0; t < NUM_PMAP_TYPES; t++)
773 <      if (photonMaps [t]) {
774 <         if (photonRepTime) {
775 <            sprintf(errmsg, "\nBuilding %s photon map...\n", pmapName [t]);
773 >      if (pmaps [t]) {
774 >         if (verbose) {
775 >            sprintf(errmsg, "Building %s photon map\n", pmapName [t]);
776              eputs(errmsg);
777 + #if NIX            
778              fflush(stderr);
779 + #endif            
780           }
781 <      
782 <         balancePhotons(photonMaps [t], &totalFlux);
781 >        
782 >         /* Build underlying data structure; heap is destroyed */
783 >         buildPhotonMap(pmaps [t], &totalFlux, NULL, numProc);
784        }
785        
786     /* Precompute photon irradiance if necessary */
787 <   if (preCompPmap)
787 >   if (preCompPmap) {
788 >      if (verbose)
789 >         eputs("\n");
790        preComputeGlobal(preCompPmap);
791 < }
614 <
615 <
616 <
617 < void photonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
618 < /* Photon density estimate. Returns irradiance at ray -> rop. */
619 < {
620 <   unsigned i;
621 <   PhotonSQNode *sq;
622 <   float r;
623 <   COLOR flux;
624 <
625 <   setcolor(irrad, 0, 0, 0);
626 <
627 <   if (!pmap -> maxGather)
628 <      return;
629 <      
630 <   /* Ignore sources */
631 <   if (ray -> ro)
632 <      if (islight(objptr(ray -> ro -> omod) -> otype))
633 <         return;
634 <        
635 <   pmap -> squeueEnd = 0;
636 <   findPhotons(pmap, ray);
791 >   }      
792    
793 <   /* Need at least 2 photons */
794 <   if (pmap -> squeueEnd < 2) {
640 <      #ifdef PMAP_NONEFOUND  
641 <         sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)",
642 <                 ray -> ro ? ray -> ro -> oname : "<null>",
643 <                 ray -> rop [0], ray -> rop [1], ray -> rop [2]);
644 <         error(WARNING, errmsg);
645 <      #endif      
646 <
647 <      return;
648 <   }
649 <      
650 <   if (pmap -> minGather == pmap -> maxGather) {
651 <      /* No bias compensation. Just do a plain vanilla estimate */
652 <      sq = pmap -> squeue + 1;
653 <      
654 <      /* Average radius between furthest two photons to improve accuracy */      
655 <      r = max(sq -> dist, (sq + 1) -> dist);
656 <      r = 0.25 * (pmap -> maxDist + r + 2 * sqrt(pmap -> maxDist * r));  
657 <      
658 <      /* Skip the extra photon */
659 <      for (i = 1 ; i < pmap -> squeueEnd; i++, sq++) {
660 <         getPhotonFlux(sq -> photon, flux);        
661 < #ifdef PMAP_EPANECHNIKOV
662 <         /* Apply Epanechnikov kernel to photon flux (dists are squared) */
663 <         scalecolor(flux, 2 * (1 - sq -> dist / r));
664 < #endif        
665 <         addcolor(irrad, flux);
666 <      }
667 <      
668 <      /* Divide by search area PI * r^2, 1 / PI required as ambient
669 <         normalisation factor */        
670 <      scalecolor(irrad, 1 / (PI * PI * r));
671 <      
672 <      return;
673 <   }
674 <   else
675 <      /* Apply bias compensation to density estimate */
676 <      biasComp(pmap, irrad);
677 < }
678 <
679 <
680 <
681 < void photonPreCompDensity (PhotonMap *pmap, RAY *r, COLOR irrad)
682 < /* Returns precomputed photon density estimate at ray -> rop. */
683 < {
684 <   Photon *p;
685 <  
686 <   setcolor(irrad, 0, 0, 0);
687 <
688 <   /* Ignore sources */
689 <   if (r -> ro && islight(objptr(r -> ro -> omod) -> otype))
690 <      return;
691 <      
692 <   if ((p = find1Photon(preCompPmap, r)))
693 <      getPhotonFlux(p, irrad);
694 < }
695 <
696 <
697 <
698 < void volumePhotonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
699 < /* Photon volume density estimate. Returns irradiance at ray -> rop. */
700 < {
701 <   unsigned i;
702 <   PhotonSQNode *sq;
703 <   float gecc2, r, ph;
704 <   COLOR flux;
705 <
706 <   setcolor(irrad, 0, 0, 0);
707 <  
708 <   if (!pmap -> maxGather)
709 <      return;
710 <      
711 <   pmap -> squeueEnd = 0;
712 <   findPhotons(pmap, ray);
713 <  
714 <   /* Need at least 2 photons */
715 <   if (pmap -> squeueEnd < 2)
716 <      return;
717 <      
718 <   if (pmap -> minGather == pmap -> maxGather) {
719 <      /* No bias compensation. Just do a plain vanilla estimate */
720 <      gecc2 = ray -> gecc * ray -> gecc;
721 <      sq = pmap -> squeue + 1;
722 <      
723 <      /* Average radius between furthest two photons to improve accuracy */      
724 <      r = max(sq -> dist, (sq + 1) -> dist);
725 <      r = 0.25 * (pmap -> maxDist + r + 2 * sqrt(pmap -> maxDist * r));  
726 <      
727 <      /* Skip the extra photon */
728 <      for (i = 1 ; i < pmap -> squeueEnd; i++, sq++) {
729 <         /* Compute phase function for inscattering from photon */
730 <         if (gecc2 <= FTINY)
731 <            ph = 1;
732 <         else {
733 <            ph = DOT(ray -> rdir, sq -> photon -> norm) / 127;
734 <            ph = 1 + gecc2 - 2 * ray -> gecc * ph;
735 <            ph = (1 - gecc2) / (ph * sqrt(ph));
736 <         }
737 <        
738 <         getPhotonFlux(sq -> photon, flux);
739 <         scalecolor(flux, ph);
740 <         addcolor(irrad, flux);
741 <      }
742 <      
743 <      /* Divide by search volume 4 / 3 * PI * r^3 and phase function
744 <         normalization factor 1 / (4 * PI) */
745 <      scalecolor(irrad, 3 / (16 * PI * PI * r * sqrt(r)));
746 <      
747 <      return;
748 <   }
749 <  
750 <   else
751 <      /* Apply bias compensation to density estimate */
752 <      volumeBiasComp(pmap, ray, irrad);
793 >   if (verbose)
794 >      eputs("\n");
795   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines