ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/Development/ray/src/rt/pmap.c

Comparing ray/src/rt/pmap.c (file contents):
Revision 2.10 by greg, Tue Sep 1 16:27:52 2015 UTC vs.
Revision 2.12 by greg, Mon Sep 26 20:19:30 2016 UTC

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

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines (old)
>	Changed lines (new)