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

Diff Legend

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