ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/pmapcontrib.c
Revision: 2.2
Committed: Tue Apr 21 19:16:51 2015 UTC (9 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.1: +1 -2 lines
Log Message:
Fixes for Windows and photon map

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