ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapcontrib.c
Revision: 2.8
Committed: Thu May 21 13:54:59 2015 UTC (9 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.7: +3 -3 lines
Log Message:
Added calls to findmaterial() and simplified photon map shadow check

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 rschregle 2.4 (c) Lucerne University of Applied Sciences and Arts,
7     supported by the Swiss National Science Foundation (SNSF, #147053)
8 greg 2.1 ==================================================================
9    
10 greg 2.8 $Id: pmapcontrib.c,v 2.7 2015/05/20 14:44:12 greg Exp $
11 greg 2.1 */
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    
24    
25    
26     static void setPmapContribParams (PhotonMap *pmap, LUTAB *srcContrib)
27     /* Set parameters for light source contributions */
28     {
29     /* Set light source modifier list and appropriate callback to extract
30     * their contributions from the photon map */
31     if (pmap) {
32     pmap -> srcContrib = srcContrib;
33     pmap -> lookup = photonContrib;
34     /* Ensure we get all requested photon contribs during lookups */
35     pmap -> gatherTolerance = 1.0;
36     }
37     }
38    
39    
40    
41     static void checkPmapContribs (const PhotonMap *pmap, LUTAB *srcContrib)
42     /* Check modifiers for light source contributions */
43     {
44     const PhotonPrimary *primary = pmap -> primary;
45     OBJREC *srcMod;
46     unsigned long i, found = 0;
47    
48     /* Make sure at least one of the modifiers is actually in the pmap,
49     * otherwise findPhotons() winds up in an infinite loop! */
50     for (i = pmap -> primarySize; i; --i, ++primary) {
51     if (primary -> srcIdx < 0 || primary -> srcIdx >= nsources)
52     error(INTERNAL, "invalid light source index in photon map");
53    
54 greg 2.8 srcMod = findmaterial(source [primary -> srcIdx].so);
55 greg 2.1 if ((MODCONT*)lu_find(srcContrib, srcMod -> oname) -> data)
56     ++found;
57     }
58    
59     if (!found)
60     error(USER, "modifiers not in photon map");
61     }
62    
63    
64    
65     void initPmapContrib (LUTAB *srcContrib, unsigned numSrcContrib)
66     {
67     unsigned t;
68    
69     for (t = 0; t < NUM_PMAP_TYPES; t++)
70     if (photonMaps [t] && t != PMAP_TYPE_CONTRIB) {
71     sprintf(errmsg, "%s photon map does not support contributions",
72     pmapName [t]);
73     error(USER, errmsg);
74     }
75    
76     /* Get params */
77     setPmapContribParams(contribPmap, srcContrib);
78    
79     if (contribPhotonMapping) {
80     if (contribPmap -> maxGather < numSrcContrib) {
81     /* Adjust density estimate bandwidth if lower than modifier
82     * count, otherwise contributions are missing */
83     error(WARNING, "contrib density estimate bandwidth too low, "
84     "adjusting to modifier count");
85     contribPmap -> maxGather = numSrcContrib;
86     }
87    
88     /* Sanity check */
89     checkPmapContribs(contribPmap, srcContrib);
90     }
91     }
92    
93    
94    
95     void photonContrib (PhotonMap *pmap, RAY *ray, COLOR irrad)
96     /* Sum up light source contributions from photons in pmap->srcContrib */
97     {
98     unsigned i;
99     PhotonSQNode *sq;
100     float r, invArea;
101     RREAL rayCoeff [3];
102 greg 2.5
103 greg 2.1 setcolor(irrad, 0, 0, 0);
104    
105     if (!pmap -> maxGather)
106     return;
107    
108     /* Ignore sources */
109     if (ray -> ro)
110     if (islight(objptr(ray -> ro -> omod) -> otype))
111     return;
112    
113 greg 2.5 /* Get cumulative path
114 greg 2.1 * coefficient up to photon lookup point */
115     raycontrib(rayCoeff, ray, PRIMARY);
116    
117     /* Lookup photons */
118     pmap -> squeueEnd = 0;
119     findPhotons(pmap, ray);
120    
121     /* Need at least 2 photons */
122     if (pmap -> squeueEnd < 2) {
123     #ifdef PMAP_NONEFOUND
124     sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)",
125     ray -> ro ? ray -> ro -> oname : "<null>",
126     ray -> rop [0], ray -> rop [1], ray -> rop [2]);
127     error(WARNING, errmsg);
128     #endif
129    
130     return;
131     }
132    
133     /* Average (squared) radius between furthest two photons to improve
134     * accuracy and get inverse search area 1 / (PI * r^2), with extra
135     * normalisation factor 1 / PI for ambient calculation */
136     sq = pmap -> squeue + 1;
137     r = max(sq -> dist, (sq + 1) -> dist);
138     r = 0.25 * (pmap -> maxDist + r + 2 * sqrt(pmap -> maxDist * r));
139     invArea = 1 / (PI * PI * r);
140    
141     /* Skip the extra photon */
142     for (i = 1 ; i < pmap -> squeueEnd; i++, sq++) {
143     COLOR flux;
144    
145     /* Get photon's contribution to density estimate */
146     getPhotonFlux(sq -> photon, flux);
147     scalecolor(flux, invArea);
148     #ifdef PMAP_EPANECHNIKOV
149     /* Apply Epanechnikov kernel to photon flux (dists are squared) */
150     scalecolor(flux, 2 * (1 - sq -> dist / r));
151     #endif
152     addcolor(irrad, flux);
153    
154     if (pmap -> srcContrib) {
155     const PhotonPrimary *primary = pmap -> primary +
156     sq -> photon -> primary;
157 greg 2.7 const SRCREC *sp = &source[primary -> srcIdx];
158 greg 2.8 OBJREC *srcMod = findmaterial(sp -> so);
159 greg 2.1 MODCONT *srcContrib = (MODCONT*)lu_find(pmap -> srcContrib,
160     srcMod -> oname) -> data;
161 greg 2.6 if (!srcContrib)
162     continue;
163    
164     /* Photon's emitting light source has modifier whose
165     * contributions are sought */
166     double srcBinReal;
167     int srcBin;
168     RAY srcRay;
169    
170     if (srcContrib -> binv -> type != NUM) {
171     /* Use intersection function to set shadow ray parameters
172     * if it's not simply a constant
173     */
174     rayorigin(&srcRay, SHADOW, NULL, NULL);
175     srcRay.rsrc = primary -> srcIdx;
176     VCOPY(srcRay.rorg, primary -> pos);
177 greg 2.7 decodedir(srcRay.rdir, primary -> dir);
178    
179     if (!(sp->sflags & SDISTANT ? sourcehit(&srcRay)
180     : (*ofun[sp -> so -> otype].funp)(sp -> so, &srcRay)))
181 greg 2.6 continue; /* XXX shouldn't happen! */
182 greg 2.7
183 greg 2.6 worldfunc(RCCONTEXT, &srcRay);
184     set_eparams((char *)srcContrib -> params);
185     }
186    
187     if ((srcBinReal = evalue(srcContrib -> binv)) < -.5)
188     continue; /* silently ignore negative bins */
189 greg 2.5
190 greg 2.6 if ((srcBin = srcBinReal + .5) >= srcContrib -> nbins) {
191     error(WARNING, "bad bin number (ignored)");
192     continue;
193     }
194 greg 2.1
195 greg 2.6 if (!contrib) {
196     /* Ray coefficient mode; normalise by light source radiance
197     * after applying distrib pattern */
198     int j;
199     raytexture(ray, srcMod -> omod);
200     setcolor(ray -> rcol, srcMod -> oargs.farg [0],
201 greg 2.1 srcMod -> oargs.farg [1], srcMod -> oargs.farg [2]);
202 greg 2.6 multcolor(ray -> rcol, ray -> pcol);
203     for (j = 0; j < 3; j++)
204     flux [j] = ray -> rcol [j] ? flux [j] / ray -> rcol [j]
205 greg 2.1 : 0;
206 greg 2.6 }
207 greg 2.1
208 greg 2.6 multcolor(flux, rayCoeff);
209     addcolor(srcContrib -> cbin [srcBin], flux);
210 greg 2.1 }
211     }
212 greg 2.5
213 greg 2.1 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 rschregle 2.3 #ifdef SIGCONT
273     signal(SIGCONT, pmapDistribReport);
274     #endif
275 greg 2.1
276     for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
277     unsigned portCnt = 0, passCnt = 0, prePassCnt = 0;
278     double srcNumEmit = 0; /* # photons to emit from source */
279     unsigned long srcNumDistrib = pm -> heapEnd; /* # photons stored */
280    
281     srcFlux [srcIdx] = repProgress = 0;
282     emap.src = source + srcIdx;
283    
284     if (photonRepTime)
285     eputs("\n");
286    
287     /* =============================================================
288     * FLUX INTEGRATION - Get total flux emitted from light source
289     * ============================================================= */
290     do {
291     emap.port = emap.src -> sflags & SDISTANT
292     ? photonPorts + portCnt : NULL;
293     photonPartition [emap.src -> so -> otype] (&emap);
294    
295     if (photonRepTime) {
296     sprintf(errmsg, "Integrating flux from source %s (mod %s) ",
297     source [srcIdx].so -> oname,
298     objptr(source [srcIdx].so -> omod) -> oname);
299    
300     if (emap.port) {
301     sprintf(errmsg2, "via port %s ",
302     photonPorts [portCnt].so -> oname);
303     strcat(errmsg, errmsg2);
304     }
305    
306     sprintf(errmsg2, "(%lu partitions)...\n",
307     emap.numPartitions);
308     strcat(errmsg, errmsg2);
309     eputs(errmsg);
310     fflush(stderr);
311     }
312    
313     for (emap.partitionCnt = 0;
314     emap.partitionCnt < emap.numPartitions;
315     emap.partitionCnt++) {
316     initPhotonEmission(&emap, pdfSamples);
317     srcFlux [srcIdx] += colorAvg(emap.partFlux);
318     }
319    
320     portCnt++;
321     } while (portCnt < numPhotonPorts);
322    
323     if (srcFlux [srcIdx] < FTINY) {
324     sprintf(errmsg, "source %s has zero emission",
325     source [srcIdx].so -> oname);
326     error(WARNING, errmsg);
327     }
328     else {
329     /* ==========================================================
330     * 2-PASS PHOTON DISTRIBUTION
331     * Pass 1 (pre): emit fraction of target photon count
332     * Pass 2 (main): based on outcome of pass 1, estimate
333     * remaining number of photons to emit to
334     * approximate target count
335     * ========================================================== */
336     do {
337     if (!passCnt) {
338     /* INIT PASS 1 */
339     if (++prePassCnt > maxPreDistrib) {
340     /* Warn if no photons contributed after sufficient
341     * iterations */
342     sprintf(errmsg, "too many prepasses, no photons "
343     "from source %s", source [srcIdx].so -> oname);
344     error(WARNING, errmsg);
345     break;
346     }
347    
348     /* Num to emit is fraction of target count */
349     srcNumEmit = preDistrib * srcDistribTarget;
350     }
351    
352     else {
353     /* INIT PASS 2 */
354     /* Based on the outcome of the predistribution we can now
355     * figure out how many more photons we have to emit from
356     * the current source to meet the target count,
357     * srcDistribTarget. This value is clamped to 0 in case
358     * the target has already been exceeded in pass 1.
359     * srcNumEmit and srcNumDistrib is the number of photons
360     * emitted and distributed (stored) from the current
361     * source in pass 1, respectively. */
362     srcNumDistrib = pm -> heapEnd - srcNumDistrib;
363     srcNumEmit *= srcNumDistrib
364     ? max(srcDistribTarget/srcNumDistrib, 1) - 1
365     : 0;
366    
367     if (!srcNumEmit)
368     /* No photons left to distribute in main pass */
369     break;
370     }
371    
372     /* Set completion count for progress report */
373     repComplete = srcNumEmit + repProgress;
374     portCnt = 0;
375    
376     do {
377     emap.port = emap.src -> sflags & SDISTANT
378     ? photonPorts + portCnt : NULL;
379     photonPartition [emap.src -> so -> otype] (&emap);
380    
381     if (photonRepTime) {
382     if (!passCnt)
383     sprintf(errmsg, "PREPASS %d on source %s (mod %s) ",
384     prePassCnt, source [srcIdx].so -> oname,
385     objptr(source[srcIdx].so->omod) -> oname);
386     else
387     sprintf(errmsg, "MAIN PASS on source %s (mod %s) ",
388     source [srcIdx].so -> oname,
389     objptr(source[srcIdx].so->omod) -> oname);
390    
391     if (emap.port) {
392     sprintf(errmsg2, "via port %s ",
393     photonPorts [portCnt].so -> oname);
394     strcat(errmsg, errmsg2);
395     }
396    
397     sprintf(errmsg2, "(%lu partitions)...\n",
398     emap.numPartitions);
399     strcat(errmsg, errmsg2);
400     eputs(errmsg);
401     fflush(stderr);
402     }
403    
404     for (emap.partitionCnt = 0;
405     emap.partitionCnt < emap.numPartitions;
406     emap.partitionCnt++) {
407     double partNumEmit;
408     unsigned long partEmitCnt;
409    
410     /* Get photon origin within current source partishunn
411     * and build emission map */
412     photonOrigin [emap.src -> so -> otype] (&emap);
413     initPhotonEmission(&emap, pdfSamples);
414    
415     /* Number of photons to emit from ziss partishunn;
416     * scale according to its normalised contribushunn to
417     * the emitted source flux */
418     partNumEmit = srcNumEmit * colorAvg(emap.partFlux) /
419     srcFlux [srcIdx];
420     partEmitCnt = (unsigned long)partNumEmit;
421    
422     /* Probabilistically account for fractional photons */
423     if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
424     partEmitCnt++;
425    
426     /* Integer counter avoids FP rounding errors */
427     while (partEmitCnt--) {
428     RAY photonRay;
429    
430     /* Emit photon according to PDF (if any), allocate
431     * associated primary ray, and trace through scene
432     * until absorbed/leaked */
433     emitPhoton(&emap, &photonRay);
434     addPhotonPrimary(pm, &photonRay);
435     tracePhoton(&photonRay);
436    
437     /* Record progress */
438     repProgress++;
439    
440     if (photonRepTime > 0 &&
441     time(NULL) >= repLastTime + photonRepTime)
442     pmapDistribReport();
443 rschregle 2.3 #ifdef SIGCONT
444 greg 2.1 else signal(SIGCONT, pmapDistribReport);
445     #endif
446     }
447     }
448    
449     portCnt++;
450     } while (portCnt < numPhotonPorts);
451    
452     if (pm -> heapEnd == srcNumDistrib)
453     /* Double preDistrib in case no photons were stored
454     * for this source and redo pass 1 */
455     preDistrib *= 2;
456     else {
457     /* Now do pass 2 */
458     passCnt++;
459     if (photonRepTime)
460     eputs("\n");
461     }
462     } while (passCnt < 2);
463    
464     /* Flux per photon emitted from this source; repProgress is the
465     * number of emitted photons after both passes */
466     srcFlux [srcIdx] = repProgress ? srcFlux [srcIdx] / repProgress
467     : 0;
468     }
469     }
470    
471     /* ================================================================
472     * POST-DISTRIBUTION - Set photon flux and build kd-tree, etc.
473     * ================================================================ */
474 rschregle 2.3 #ifdef SIGCONT
475     signal(SIGCONT, SIG_DFL);
476     #endif
477 greg 2.1 free(emap.samples);
478    
479     if (!pm -> heapEnd)
480     error(USER, "empty photon map");
481    
482     /* Check for valid primary photon rays */
483     if (!pm -> primary)
484     error(INTERNAL, "no primary rays in contribution photon map");
485    
486     if (pm -> primary [pm -> primaryEnd].srcIdx < 0)
487     /* Last primary ray is unused, so decrement counter */
488     pm -> primaryEnd--;
489    
490     if (photonRepTime) {
491     eputs("\nBuilding contrib photon heap...\n");
492     fflush(stderr);
493     }
494    
495     balancePhotons(pm, srcFlux);
496     }