ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapcontrib.c
Revision: 2.1
Committed: Tue Feb 24 19:39:26 2015 UTC (9 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Initial check-in of photon map addition by Roland Schregle

File Contents

# User Rev Content
1 greg 2.1 /*
2     ==================================================================
3     Photon map support for light source contributions
4    
5     Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
6     (c) Fraunhofer Institute for Solar Energy Systems,
7     Lucerne University of Applied Sciences & Arts
8     ==================================================================
9    
10     $Id: pmapcontrib.c,v 4.27 2015/02/12 09:57:21 taschreg Exp taschreg $
11     */
12    
13    
14     #include "pmapcontrib.h"
15     #include "pmap.h"
16     #include "pmapmat.h"
17     #include "pmapsrc.h"
18     #include "pmaprand.h"
19     #include "pmapio.h"
20     #include "pmapdiag.h"
21     #include "rcontrib.h"
22     #include "otypes.h"
23     #include <signal.h>
24    
25    
26    
27     extern int contrib; /* coeff/contrib flag */
28    
29    
30    
31     static void setPmapContribParams (PhotonMap *pmap, LUTAB *srcContrib)
32     /* Set parameters for light source contributions */
33     {
34     /* Set light source modifier list and appropriate callback to extract
35     * their contributions from the photon map */
36     if (pmap) {
37     pmap -> srcContrib = srcContrib;
38     pmap -> lookup = photonContrib;
39     /* Ensure we get all requested photon contribs during lookups */
40     pmap -> gatherTolerance = 1.0;
41     }
42     }
43    
44    
45    
46     static void checkPmapContribs (const PhotonMap *pmap, LUTAB *srcContrib)
47     /* Check modifiers for light source contributions */
48     {
49     const PhotonPrimary *primary = pmap -> primary;
50     OBJREC *srcMod;
51     unsigned long i, found = 0;
52    
53     /* Make sure at least one of the modifiers is actually in the pmap,
54     * otherwise findPhotons() winds up in an infinite loop! */
55     for (i = pmap -> primarySize; i; --i, ++primary) {
56     if (primary -> srcIdx < 0 || primary -> srcIdx >= nsources)
57     error(INTERNAL, "invalid light source index in photon map");
58    
59     srcMod = objptr(source [primary -> srcIdx].so -> omod);
60     if ((MODCONT*)lu_find(srcContrib, srcMod -> oname) -> data)
61     ++found;
62     }
63    
64     if (!found)
65     error(USER, "modifiers not in photon map");
66     }
67    
68    
69    
70     void initPmapContrib (LUTAB *srcContrib, unsigned numSrcContrib)
71     {
72     unsigned t;
73    
74     for (t = 0; t < NUM_PMAP_TYPES; t++)
75     if (photonMaps [t] && t != PMAP_TYPE_CONTRIB) {
76     sprintf(errmsg, "%s photon map does not support contributions",
77     pmapName [t]);
78     error(USER, errmsg);
79     }
80    
81     /* Get params */
82     setPmapContribParams(contribPmap, srcContrib);
83    
84     if (contribPhotonMapping) {
85     if (contribPmap -> maxGather < numSrcContrib) {
86     /* Adjust density estimate bandwidth if lower than modifier
87     * count, otherwise contributions are missing */
88     error(WARNING, "contrib density estimate bandwidth too low, "
89     "adjusting to modifier count");
90     contribPmap -> maxGather = numSrcContrib;
91     }
92    
93     /* Sanity check */
94     checkPmapContribs(contribPmap, srcContrib);
95     }
96     }
97    
98    
99    
100     void photonContrib (PhotonMap *pmap, RAY *ray, COLOR irrad)
101     /* Sum up light source contributions from photons in pmap->srcContrib */
102     {
103     unsigned i;
104     PhotonSQNode *sq;
105     float r, invArea;
106     RREAL rayCoeff [3];
107     FVECT rdir, rop;
108    
109     setcolor(irrad, 0, 0, 0);
110    
111     if (!pmap -> maxGather)
112     return;
113    
114     /* Ignore sources */
115     if (ray -> ro)
116     if (islight(objptr(ray -> ro -> omod) -> otype))
117     return;
118    
119     /* Set context for binning function evaluation and get cumulative path
120     * coefficient up to photon lookup point */
121     worldfunc(RCCONTEXT, ray);
122     raycontrib(rayCoeff, ray, PRIMARY);
123    
124     /* Save incident ray's direction and hitpoint */
125     VCOPY(rdir, ray -> rdir);
126     VCOPY(rop, ray -> rop);
127    
128     /* Lookup photons */
129     pmap -> squeueEnd = 0;
130     findPhotons(pmap, ray);
131    
132     /* Need at least 2 photons */
133     if (pmap -> squeueEnd < 2) {
134     #ifdef PMAP_NONEFOUND
135     sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)",
136     ray -> ro ? ray -> ro -> oname : "<null>",
137     ray -> rop [0], ray -> rop [1], ray -> rop [2]);
138     error(WARNING, errmsg);
139     #endif
140    
141     return;
142     }
143    
144     /* Average (squared) radius between furthest two photons to improve
145     * accuracy and get inverse search area 1 / (PI * r^2), with extra
146     * normalisation factor 1 / PI for ambient calculation */
147     sq = pmap -> squeue + 1;
148     r = max(sq -> dist, (sq + 1) -> dist);
149     r = 0.25 * (pmap -> maxDist + r + 2 * sqrt(pmap -> maxDist * r));
150     invArea = 1 / (PI * PI * r);
151    
152     /* Skip the extra photon */
153     for (i = 1 ; i < pmap -> squeueEnd; i++, sq++) {
154     COLOR flux;
155    
156     /* Get photon's contribution to density estimate */
157     getPhotonFlux(sq -> photon, flux);
158     scalecolor(flux, invArea);
159     #ifdef PMAP_EPANECHNIKOV
160     /* Apply Epanechnikov kernel to photon flux (dists are squared) */
161     scalecolor(flux, 2 * (1 - sq -> dist / r));
162     #endif
163     addcolor(irrad, flux);
164    
165     if (pmap -> srcContrib) {
166     const PhotonPrimary *primary = pmap -> primary +
167     sq -> photon -> primary;
168     OBJREC *srcMod = objptr(source [primary -> srcIdx].so -> omod);
169     MODCONT *srcContrib = (MODCONT*)lu_find(pmap -> srcContrib,
170     srcMod -> oname) -> data;
171    
172     if (srcContrib) {
173     /* Photon's emitting light source has modifier whose
174     * contributions are sought */
175     int srcBin;
176    
177     /* Set incident dir and origin of photon's primary ray on
178     * light source for dummy shadow ray, and evaluate binning
179     * function */
180     VCOPY(ray -> rdir, primary -> dir);
181     VCOPY(ray -> rop, primary -> org);
182     srcBin = evalue(srcContrib -> binv) + .5;
183    
184     if (srcBin < 0 || srcBin >= srcContrib -> nbins) {
185     error(WARNING, "bad bin number (ignored)");
186     continue;
187     }
188    
189     if (!contrib) {
190     /* Ray coefficient mode; normalise by light source radiance
191     * after applying distrib pattern */
192     int j;
193     raytexture(ray, srcMod -> omod);
194     setcolor(ray -> rcol, srcMod -> oargs.farg [0],
195     srcMod -> oargs.farg [1], srcMod -> oargs.farg [2]);
196     multcolor(ray -> rcol, ray -> pcol);
197     for (j = 0; j < 3; j++)
198     flux [j] = ray -> rcol [j] ? flux [j] / ray -> rcol [j]
199     : 0;
200     }
201    
202     multcolor(flux, rayCoeff);
203     addcolor(srcContrib -> cbin [srcBin], flux);
204     }
205     else fprintf(stderr, "Skipped contrib from %s\n", srcMod -> oname);
206     }
207     }
208    
209     /* Restore incident ray's direction and hitpoint */
210     VCOPY(ray -> rdir, rdir);
211     VCOPY(ray -> rop, rop);
212    
213     return;
214     }
215    
216    
217    
218     void distribPhotonContrib (PhotonMap* pm)
219     {
220     EmissionMap emap;
221     char errmsg2 [128];
222     unsigned srcIdx;
223     double *srcFlux; /* Emitted flux per light source */
224     const double srcDistribTarget = /* Target photon count per source */
225     nsources ? (double)pm -> distribTarget / nsources : 0;
226    
227     if (!pm)
228     error(USER, "no photon map defined");
229    
230     if (!nsources)
231     error(USER, "no light sources");
232    
233     /* Allocate photon flux per light source; this differs for every
234     * source as all sources contribute the same number of distributed
235     * photons (srcDistribTarget), hence the number of photons emitted per
236     * source does not correlate with its emitted flux. The resulting flux
237     * per photon is therefore adjusted individually for each source. */
238     if (!(srcFlux = calloc(nsources, sizeof(double))))
239     error(SYSTEM, "cannot allocate source flux");
240    
241     /* ================================================================
242     * INITIALISASHUNN - Set up emisshunn and scattering funcs
243     * ================================================================ */
244     emap.samples = NULL;
245     emap.src = NULL;
246     emap.maxPartitions = MAXSPART;
247     emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1);
248     if (!emap.partitions)
249     error(USER, "can't allocate source partitions");
250    
251     initPhotonMap(pm, PMAP_TYPE_CONTRIB);
252     initPhotonEmissionFuncs();
253     initPhotonScatterFuncs();
254    
255     /* Get photon ports if specified */
256     if (ambincl == 1)
257     getPhotonPorts();
258    
259     /* Get photon sensor modifiers */
260     getPhotonSensors(photonSensorList);
261    
262     /* Seed RNGs for photon distribution */
263     pmapSeed(randSeed, partState);
264     pmapSeed(randSeed, emitState);
265     pmapSeed(randSeed, cntState);
266     pmapSeed(randSeed, mediumState);
267     pmapSeed(randSeed, scatterState);
268     pmapSeed(randSeed, rouletteState);
269    
270     /* Record start time and enable progress report signal handler */
271     repStartTime = time(NULL);
272     signal(SIGCONT, pmapDistribReport);
273    
274     for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
275     unsigned portCnt = 0, passCnt = 0, prePassCnt = 0;
276     double srcNumEmit = 0; /* # photons to emit from source */
277     unsigned long srcNumDistrib = pm -> heapEnd; /* # photons stored */
278    
279     srcFlux [srcIdx] = repProgress = 0;
280     emap.src = source + srcIdx;
281    
282     if (photonRepTime)
283     eputs("\n");
284    
285     /* =============================================================
286     * FLUX INTEGRATION - Get total flux emitted from light source
287     * ============================================================= */
288     do {
289     emap.port = emap.src -> sflags & SDISTANT
290     ? photonPorts + portCnt : NULL;
291     photonPartition [emap.src -> so -> otype] (&emap);
292    
293     if (photonRepTime) {
294     sprintf(errmsg, "Integrating flux from source %s (mod %s) ",
295     source [srcIdx].so -> oname,
296     objptr(source [srcIdx].so -> omod) -> oname);
297    
298     if (emap.port) {
299     sprintf(errmsg2, "via port %s ",
300     photonPorts [portCnt].so -> oname);
301     strcat(errmsg, errmsg2);
302     }
303    
304     sprintf(errmsg2, "(%lu partitions)...\n",
305     emap.numPartitions);
306     strcat(errmsg, errmsg2);
307     eputs(errmsg);
308     fflush(stderr);
309     }
310    
311     for (emap.partitionCnt = 0;
312     emap.partitionCnt < emap.numPartitions;
313     emap.partitionCnt++) {
314     initPhotonEmission(&emap, pdfSamples);
315     srcFlux [srcIdx] += colorAvg(emap.partFlux);
316     }
317    
318     portCnt++;
319     } while (portCnt < numPhotonPorts);
320    
321     if (srcFlux [srcIdx] < FTINY) {
322     sprintf(errmsg, "source %s has zero emission",
323     source [srcIdx].so -> oname);
324     error(WARNING, errmsg);
325     }
326     else {
327     /* ==========================================================
328     * 2-PASS PHOTON DISTRIBUTION
329     * Pass 1 (pre): emit fraction of target photon count
330     * Pass 2 (main): based on outcome of pass 1, estimate
331     * remaining number of photons to emit to
332     * approximate target count
333     * ========================================================== */
334     do {
335     if (!passCnt) {
336     /* INIT PASS 1 */
337     if (++prePassCnt > maxPreDistrib) {
338     /* Warn if no photons contributed after sufficient
339     * iterations */
340     sprintf(errmsg, "too many prepasses, no photons "
341     "from source %s", source [srcIdx].so -> oname);
342     error(WARNING, errmsg);
343     break;
344     }
345    
346     /* Num to emit is fraction of target count */
347     srcNumEmit = preDistrib * srcDistribTarget;
348     }
349    
350     else {
351     /* INIT PASS 2 */
352     /* Based on the outcome of the predistribution we can now
353     * figure out how many more photons we have to emit from
354     * the current source to meet the target count,
355     * srcDistribTarget. This value is clamped to 0 in case
356     * the target has already been exceeded in pass 1.
357     * srcNumEmit and srcNumDistrib is the number of photons
358     * emitted and distributed (stored) from the current
359     * source in pass 1, respectively. */
360     srcNumDistrib = pm -> heapEnd - srcNumDistrib;
361     srcNumEmit *= srcNumDistrib
362     ? max(srcDistribTarget/srcNumDistrib, 1) - 1
363     : 0;
364    
365     if (!srcNumEmit)
366     /* No photons left to distribute in main pass */
367     break;
368     }
369    
370     /* Set completion count for progress report */
371     repComplete = srcNumEmit + repProgress;
372     portCnt = 0;
373    
374     do {
375     emap.port = emap.src -> sflags & SDISTANT
376     ? photonPorts + portCnt : NULL;
377     photonPartition [emap.src -> so -> otype] (&emap);
378    
379     if (photonRepTime) {
380     if (!passCnt)
381     sprintf(errmsg, "PREPASS %d on source %s (mod %s) ",
382     prePassCnt, source [srcIdx].so -> oname,
383     objptr(source[srcIdx].so->omod) -> oname);
384     else
385     sprintf(errmsg, "MAIN PASS on source %s (mod %s) ",
386     source [srcIdx].so -> oname,
387     objptr(source[srcIdx].so->omod) -> oname);
388    
389     if (emap.port) {
390     sprintf(errmsg2, "via port %s ",
391     photonPorts [portCnt].so -> oname);
392     strcat(errmsg, errmsg2);
393     }
394    
395     sprintf(errmsg2, "(%lu partitions)...\n",
396     emap.numPartitions);
397     strcat(errmsg, errmsg2);
398     eputs(errmsg);
399     fflush(stderr);
400     }
401    
402     for (emap.partitionCnt = 0;
403     emap.partitionCnt < emap.numPartitions;
404     emap.partitionCnt++) {
405     double partNumEmit;
406     unsigned long partEmitCnt;
407    
408     /* Get photon origin within current source partishunn
409     * and build emission map */
410     photonOrigin [emap.src -> so -> otype] (&emap);
411     initPhotonEmission(&emap, pdfSamples);
412    
413     /* Number of photons to emit from ziss partishunn;
414     * scale according to its normalised contribushunn to
415     * the emitted source flux */
416     partNumEmit = srcNumEmit * colorAvg(emap.partFlux) /
417     srcFlux [srcIdx];
418     partEmitCnt = (unsigned long)partNumEmit;
419    
420     /* Probabilistically account for fractional photons */
421     if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
422     partEmitCnt++;
423    
424     /* Integer counter avoids FP rounding errors */
425     while (partEmitCnt--) {
426     RAY photonRay;
427    
428     /* Emit photon according to PDF (if any), allocate
429     * associated primary ray, and trace through scene
430     * until absorbed/leaked */
431     emitPhoton(&emap, &photonRay);
432     addPhotonPrimary(pm, &photonRay);
433     tracePhoton(&photonRay);
434    
435     /* Record progress */
436     repProgress++;
437    
438     if (photonRepTime > 0 &&
439     time(NULL) >= repLastTime + photonRepTime)
440     pmapDistribReport();
441     #ifndef BSD
442     else signal(SIGCONT, pmapDistribReport);
443     #endif
444     }
445     }
446    
447     portCnt++;
448     } while (portCnt < numPhotonPorts);
449    
450     if (pm -> heapEnd == srcNumDistrib)
451     /* Double preDistrib in case no photons were stored
452     * for this source and redo pass 1 */
453     preDistrib *= 2;
454     else {
455     /* Now do pass 2 */
456     passCnt++;
457     if (photonRepTime)
458     eputs("\n");
459     }
460     } while (passCnt < 2);
461    
462     /* Flux per photon emitted from this source; repProgress is the
463     * number of emitted photons after both passes */
464     srcFlux [srcIdx] = repProgress ? srcFlux [srcIdx] / repProgress
465     : 0;
466     }
467     }
468    
469     /* ================================================================
470     * POST-DISTRIBUTION - Set photon flux and build kd-tree, etc.
471     * ================================================================ */
472     signal(SIGCONT, SIG_DFL);
473     free(emap.samples);
474    
475     if (!pm -> heapEnd)
476     error(USER, "empty photon map");
477    
478     /* Check for valid primary photon rays */
479     if (!pm -> primary)
480     error(INTERNAL, "no primary rays in contribution photon map");
481    
482     if (pm -> primary [pm -> primaryEnd].srcIdx < 0)
483     /* Last primary ray is unused, so decrement counter */
484     pm -> primaryEnd--;
485    
486     if (photonRepTime) {
487     eputs("\nBuilding contrib photon heap...\n");
488     fflush(stderr);
489     }
490    
491     balancePhotons(pm, srcFlux);
492     }