ViewVC Help

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

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

Comparing ray/src/rt/pmapcontrib.c (file contents):
Revision 2.5 by greg, Wed May 20 12:58:31 2015 UTC vs.
Revision 2.19 by greg, Thu Nov 8 00:54:07 2018 UTC

#	Line 1 | Line 1
1	+	#ifndef lint
2	+	static const char RCSid[] = "$Id$";
3	+	#endif
4	+
5		/*
6	<	==================================================================
7	<	Photon map support for light source contributions
6	>	======================================================================
7	>	Photon map for light source contributions
8
9		Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
10		(c) Lucerne University of Applied Sciences and Arts,
11	<	supported by the Swiss National Science Foundation (SNSF, #147053)
12	<	==================================================================
11	>	supported by the Swiss National Science Foundation (SNSF, #147053)
12	>	======================================================================
13
14		$Id$
15		*/
16
17
18		#include "pmapcontrib.h"
15	–	#include "pmap.h"
19		#include "pmapmat.h"
20		#include "pmapsrc.h"
21		#include "pmaprand.h"
#	Line 20 | Line 23
23		#include "pmapdiag.h"
24		#include "rcontrib.h"
25		#include "otypes.h"
26	+	#include "otspecial.h"
27	+	#if NIX
28	+	#include <sys/mman.h>
29	+	#include <sys/wait.h>
30	+	#endif
31
32
33	<
34	<	static void setPmapContribParams (PhotonMap *pmap, LUTAB *srcContrib)
35	<	/* Set parameters for light source contributions */
33	>	static PhotonPrimaryIdx newPhotonPrimary (PhotonMap *pmap,
34	>	const RAY *primRay,
35	>	FILE *primHeap)
36	>	/* Add primary ray for emitted photon and save light source index, origin on
37	>	* source, and emitted direction; used by contrib photons. The current
38	>	* primary is stored in pmap -> lastPrimary.  If the previous primary
39	>	* contributed photons (has srcIdx >= 0), it's appended to primHeap.  If
40	>	* primRay == NULL, the current primary is still flushed, but no new primary
41	>	* is set.  Returns updated primary counter pmap -> numPrimary.  */
42		{
43	<	/* Set light source modifier list and appropriate callback to extract
44	<	* their contributions from the photon map */
45	<	if (pmap) {
46	<	pmap -> srcContrib = srcContrib;
47	<	pmap -> lookup = photonContrib;
48	<	/* Ensure we get all requested photon contribs during lookups */
49	<	pmap -> gatherTolerance = 1.0;
43	>	if (!pmap || !primHeap)
44	>	return 0;
45	>
46	>	/* Check if last primary ray has spawned photons (srcIdx >= 0, see
47	>	* newPhoton()), in which case we save it to the primary heap file
48	>	* before clobbering it */
49	>	if (pmap -> lastPrimary.srcIdx >= 0) {
50	>	if (!fwrite(&pmap -> lastPrimary, sizeof(PhotonPrimary), 1, primHeap))
51	>	error(SYSTEM, "failed writing photon primary in newPhotonPrimary");
52	>
53	>	pmap -> numPrimary++;
54	>	if (pmap -> numPrimary > PMAP_MAXPRIMARY)
55	>	error(INTERNAL, "photon primary overflow in newPhotonPrimary");
56		}
37	–	}
57
58	+	/* Mark unused with negative source index until path spawns a photon (see
59	+	* newPhoton()) */
60	+	pmap -> lastPrimary.srcIdx = -1;
61	+
62	+	if (primRay) {
63	+	FVECT dvec;
64
65	<
66	<	static void checkPmapContribs (const PhotonMap *pmap, LUTAB *srcContrib)
67	<	/* Check modifiers for light source contributions */
68	<	{
69	<	const PhotonPrimary *primary = pmap -> primary;
70	<	OBJREC *srcMod;
71	<	unsigned long i, found = 0;
72	<
73	<	/* Make sure at least one of the modifiers is actually in the pmap,
74	<	* 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	<	srcMod = objptr(source [primary -> srcIdx].so -> omod);
55	<	if ((MODCONT*)lu_find(srcContrib, srcMod -> oname) -> data)
56	<	++found;
65	>	#ifdef PMAP_PRIMARYDIR
66	>	/* Reverse incident direction to point to light source */
67	>	dvec [0] = -primRay -> rdir [0];
68	>	dvec [1] = -primRay -> rdir [1];
69	>	dvec [2] = -primRay -> rdir [2];
70	>	pmap -> lastPrimary.dir = encodedir(dvec);
71	>	#endif
72	>	#ifdef PMAP_PRIMARYPOS
73	>	VCOPY(pmap -> lastPrimary.pos, primRay -> rop);
74	>	#endif
75		}
76
77	<	if (!found)
60	<	error(USER, "modifiers not in photon map");
77	>	return pmap -> numPrimary;
78		}
62	–
63	–
79
80	<	void initPmapContrib (LUTAB *srcContrib, unsigned numSrcContrib)
80	>
81	>
82	>	#ifdef DEBUG_PMAP
83	>	static int checkPrimaryHeap (FILE *file)
84	>	/* Check heap for ordered primaries */
85		{
86	<	unsigned t;
86	>	Photon   p, lastp;
87	>	int      i, dup;
88
89	<	for (t = 0; t < NUM_PMAP_TYPES; t++)
90	<	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	<	}
89	>	rewind(file);
90	>	memset(&lastp, 0, sizeof(lastp));
91
92	<	/* Get params */
93	<	setPmapContribParams(contribPmap, srcContrib);
94	<
95	<	if (contribPhotonMapping) {
96	<	if (contribPmap -> maxGather < numSrcContrib) {
97	<	/* Adjust density estimate bandwidth if lower than modifier
98	<	* count, otherwise contributions are missing */
99	<	error(WARNING, "contrib density estimate bandwidth too low, "
100	<	"adjusting to modifier count");
101	<	contribPmap -> maxGather = numSrcContrib;
92	>	while (fread(&p, sizeof(p), 1, file)) {
93	>	dup = 1;
94	>
95	>	for (i = 0; i <= 2; i++) {
96	>	if (p.pos [i] < thescene.cuorg [i] ||
97	>	p.pos [i] > thescene.cuorg [i] + thescene.cusize) {
98	>
99	>	sprintf(errmsg, "corrupt photon in heap at [%f, %f, %f]\n",
100	>	p.pos [0], p.pos [1], p.pos [2]);
101	>	error(WARNING, errmsg);
102	>	}
103	>
104	>	dup &= p.pos [i] == lastp.pos [i];
105		}
106
107	<	/* Sanity check */
108	<	checkPmapContribs(contribPmap, srcContrib);
107	>	if (dup) {
108	>	sprintf(errmsg,
109	>	"consecutive duplicate photon in heap at [%f, %f, %f]\n",
110	>	p.pos [0], p.pos [1], p.pos [2]);
111	>	error(WARNING, errmsg);
112	>	}
113		}
114	+
115	+	return 0;
116		}
117	+	#endif
118
119
120
121	<	void photonContrib (PhotonMap *pmap, RAY *ray, COLOR irrad)
122	<	/* Sum up light source contributions from photons in pmap->srcContrib */
121	>	static PhotonPrimaryIdx buildPrimaries (PhotonMap *pmap, FILE **primaryHeap,
122	>	char **primaryHeapFname,
123	>	PhotonPrimaryIdx *primaryOfs,
124	>	unsigned numHeaps)
125	>	/* Consolidate per-subprocess photon primary heaps into the primary array
126	>	* pmap -> primaries. Returns offset for primary index linearisation in
127	>	* numPrimary. The heap files in primaryHeap are closed on return. */
128		{
129	<	unsigned       i;
130	<	PhotonSQNode   *sq;
131	<	float          r, invArea;
132	<	RREAL          rayCoeff [3];
133	<
103	<	setcolor(irrad, 0, 0, 0);
104	<
105	<	if (!pmap -> maxGather)
106	<	return;
129	>	PhotonPrimaryIdx  heapLen;
130	>	unsigned          heap;
131	>
132	>	if (!pmap || !primaryHeap || !primaryOfs || !numHeaps)
133	>	return 0;
134
135	<	/* Ignore sources */
136	<	if (ray -> ro)
137	<	if (islight(objptr(ray -> ro -> omod) -> otype))
138	<	return;
135	>	pmap -> numPrimary = 0;
136	>
137	>	for (heap = 0; heap < numHeaps; heap++) {
138	>	primaryOfs [heap] = pmap -> numPrimary;
139	>
140	>	if (fseek(primaryHeap [heap], 0, SEEK_END) < 0)
141	>	error(SYSTEM, "failed photon primary seek in buildPrimaries");
142	>	pmap -> numPrimary += heapLen = ftell(primaryHeap [heap]) /
143	>	sizeof(PhotonPrimary);
144
145	<	/* Get cumulative path
146	<	* coefficient up to photon lookup point */
147	<	raycontrib(rayCoeff, ray, PRIMARY);
145	>	pmap -> primaries = realloc(pmap -> primaries,
146	>	pmap -> numPrimary *
147	>	sizeof(PhotonPrimary));
148	>	if (!pmap -> primaries)
149	>	error(SYSTEM, "failed photon primary alloc in buildPrimaries");
150
151	<	/* Lookup photons */
152	<	pmap -> squeueEnd = 0;
153	<	findPhotons(pmap, ray);
154	<
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
151	>	rewind(primaryHeap [heap]);
152	>	if (fread(pmap -> primaries + primaryOfs [heap], sizeof(PhotonPrimary),
153	>	heapLen, primaryHeap [heap]) != heapLen)
154	>	error(SYSTEM, "failed reading photon primaries in buildPrimaries");
155
156	<	return;
156	>	fclose(primaryHeap [heap]);
157	>	unlink(primaryHeapFname [heap]);
158		}
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);
159
160	<	/* Skip the extra photon */
161	<	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	<	SRCREC *sp = &source[primary -> srcIdx];
158	<	OBJREC *srcMod = objptr(sp -> so -> omod);
159	<	MODCONT *srcContrib = (MODCONT*)lu_find(pmap -> srcContrib,
160	<	srcMod -> oname) -> data;
161	<
162	<	if (srcContrib) {
163	<	/* Photon's emitting light source has modifier whose
164	<	* contributions are sought */
165	<	double srcBinReal;
166	<	int srcBin;
167	<	RAY srcRay;
160	>	return pmap -> numPrimary;
161	>	}
162
169	–	if (srcContrib -> binv -> type != NUM) {
170	–	/* Use intersection function to set shadow ray parameters
171	–	*/
172	–	rayorigin(&srcRay, SHADOW, NULL, NULL);
173	–	srcRay.rsrc = primary -> srcIdx;
174	–	VCOPY(srcRay.rorg, primary -> pos);
175	–	VCOPY(srcRay.rdir, primary -> dir);
176	–	if (!(source [primary -> srcIdx].sflags & SDISTANT ?
177	–	sourcehit(&srcRay) :
178	–	(*ofun[sp -> so -> otype].funp)(sp -> so, &srcRay)))
179	–	continue;           /* XXX shouldn't happen! */
180	–	worldfunc(RCCONTEXT, &srcRay);
181	–	set_eparams((char *)srcContrib -> params);
182	–	}
183	–	if ((srcBinReal = evalue(srcContrib -> binv)) < -.5)
184	–	continue;               /* silently ignore negative bins */
185	–
186	–	if ((srcBin = srcBinReal + .5) >= srcContrib -> nbins) {
187	–	error(WARNING, "bad bin number (ignored)");
188	–	continue;
189	–	}
190	–
191	–	if (!contrib) {
192	–	/* Ray coefficient mode; normalise by light source radiance
193	–	* after applying distrib pattern */
194	–	int j;
195	–	raytexture(ray, srcMod -> omod);
196	–	setcolor(ray -> rcol, srcMod -> oargs.farg [0],
197	–	srcMod -> oargs.farg [1], srcMod -> oargs.farg [2]);
198	–	multcolor(ray -> rcol, ray -> pcol);
199	–	for (j = 0; j < 3; j++)
200	–	flux [j] = ray -> rcol [j] ? flux [j] / ray -> rcol [j]
201	–	: 0;
202	–	}
203	–
204	–	multcolor(flux, rayCoeff);
205	–	addcolor(srcContrib -> cbin [srcBin], flux);
206	–	}
207	–	else fprintf(stderr, "Skipped contrib from %s\n", srcMod -> oname);
208	–	}
209	–	}
210	–
211	–	return;
212	–	}
163
164
165	+	/* Defs for photon emission counter array passed by sub-processes to parent
166	+	* via shared memory */
167	+	typedef  unsigned long  PhotonContribCnt;
168
169	<	void distribPhotonContrib (PhotonMap* pm)
170	<	{
171	<	EmissionMap emap;
172	<	char errmsg2 [128];
220	<	unsigned srcIdx;
221	<	double *srcFlux;                 /* Emitted flux per light source */
222	<	const double srcDistribTarget =  /* Target photon count per source */
223	<	nsources ? (double)pm -> distribTarget / nsources : 0;
169	>	/* Indices for photon emission counter array: num photons stored and num
170	>	* emitted per source */
171	>	#define  PHOTONCNT_NUMPHOT    0
172	>	#define  PHOTONCNT_NUMEMIT(n) (1 + n)
173
174	+
175	+
176	+
177	+
178	+
179	+	void distribPhotonContrib (PhotonMap* pm, unsigned numProc)
180	+	{
181	+	EmissionMap       emap;
182	+	char              errmsg2 [128], shmFname [PMAP_TMPFNLEN];
183	+	unsigned          srcIdx, proc;
184	+	int               shmFile, stat, pid;
185	+	double            *srcFlux,         /* Emitted flux per light source */
186	+	srcDistribTarget; /* Target photon count per source */
187	+	PhotonContribCnt  *photonCnt;       /* Photon emission counter array */
188	+	unsigned          photonCntSize = sizeof(PhotonContribCnt) *
189	+	PHOTONCNT_NUMEMIT(nsources);
190	+	FILE              **primaryHeap = NULL;
191	+	char              **primaryHeapFname = NULL;
192	+	PhotonPrimaryIdx  *primaryOfs = NULL;
193	+
194		if (!pm)
195	<	error(USER, "no photon map defined");
195	>	error(USER, "no photon map defined in distribPhotonContrib");
196
197		if (!nsources)
198	<	error(USER, "no light sources");
199	<
198	>	error(USER, "no light sources in distribPhotonContrib");
199	>
200	>	if (nsources > MAXMODLIST)
201	>	error(USER, "too many light sources in distribPhotonContrib");
202	>
203		/* Allocate photon flux per light source; this differs for every
204		* source as all sources contribute the same number of distributed
205		* photons (srcDistribTarget), hence the number of photons emitted per
206		* source does not correlate with its emitted flux. The resulting flux
207		* per photon is therefore adjusted individually for each source. */
208		if (!(srcFlux = calloc(nsources, sizeof(double))))
209	<	error(SYSTEM, "cannot allocate source flux");
209	>	error(SYSTEM, "can't allocate source flux in distribPhotonContrib");
210
211	<	/* ================================================================
212	<	* INITIALISASHUNN - Set up emisshunn and scattering funcs
213	<	* ================================================================ */
211	>	/* ===================================================================
212	>	* INITIALISATION - Set up emission and scattering funcs
213	>	* =================================================================== */
214		emap.samples = NULL;
215		emap.src = NULL;
216		emap.maxPartitions = MAXSPART;
217		emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1);
218		if (!emap.partitions)
219	<	error(USER, "can't allocate source partitions");
219	>	error(USER, "can't allocate source partitions in distribPhotonContrib");
220
221	+	/* Initialise contrib photon map */
222		initPhotonMap(pm, PMAP_TYPE_CONTRIB);
223	+	initPhotonHeap(pm);
224		initPhotonEmissionFuncs();
225		initPhotonScatterFuncs();
226
227	<	/* Get photon ports if specified */
228	<	if (ambincl == 1)
229	<	getPhotonPorts();
227	>	/* Per-subprocess / per-source target counts */
228	>	pm -> distribTarget /= numProc;
229	>	srcDistribTarget = nsources ? (double)pm -> distribTarget / nsources : 0;
230	>
231	>	if (!pm -> distribTarget)
232	>	error(INTERNAL, "no photons to distribute in distribPhotonContrib");
233	>
234	>	/* Get photon ports from modifier list */
235	>	getPhotonPorts(photonPortList);
236
237		/* Get photon sensor modifiers */
238		getPhotonSensors(photonSensorList);
239	+
240	+	#if NIX
241	+	/* Set up shared mem for photon counters (zeroed by ftruncate) */
242	+	strcpy(shmFname, PMAP_TMPFNAME);
243	+	shmFile = mkstemp(shmFname);
244
245	<	/* Seed RNGs for photon distribution */
246	<	pmapSeed(randSeed, partState);
262	<	pmapSeed(randSeed, emitState);
263	<	pmapSeed(randSeed, cntState);
264	<	pmapSeed(randSeed, mediumState);
265	<	pmapSeed(randSeed, scatterState);
266	<	pmapSeed(randSeed, rouletteState);
245	>	if (shmFile < 0 || ftruncate(shmFile, photonCntSize) < 0)
246	>	error(SYSTEM, "failed shared mem init in distribPhotonContrib");
247
248	<	/* Record start time and enable progress report signal handler */
249	<	repStartTime = time(NULL);
250	<	#ifdef SIGCONT
251	<	signal(SIGCONT, pmapDistribReport);
252	<	#endif
248	>	photonCnt = mmap(NULL, photonCntSize, PROT_READ | PROT_WRITE,
249	>	MAP_SHARED, shmFile, 0);
250	>
251	>	if (photonCnt == MAP_FAILED)
252	>	error(SYSTEM, "failed shared mem mapping in distribPhotonContrib");
253	>	#else
254	>	/* Allocate photon counters statically on Windoze */
255	>	if (!(photonCnt = malloc(photonCntSize)))
256	>	error(SYSTEM, "failed trivial malloc in distribPhotonContrib");
257
258	<	for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
259	<	unsigned portCnt = 0, passCnt = 0, prePassCnt = 0;
260	<	double srcNumEmit = 0;          /* # photons to emit from source */
277	<	unsigned long srcNumDistrib = pm -> heapEnd; /* # photons stored */
258	>	for (srcIdx = 0; srcIdx < PHOTONCNT_NUMEMIT(nsources); srcIdx++)
259	>	photonCnt [srcIdx] = 0;
260	>	#endif /* NIX */
261
262	<	srcFlux [srcIdx] = repProgress = 0;
262	>	if (verbose) {
263	>	sprintf(errmsg, "\nIntegrating flux from %d sources", nsources);
264	>
265	>	if (photonPorts) {
266	>	sprintf(errmsg2, " via %d ports", numPhotonPorts);
267	>	strcat(errmsg, errmsg2);
268	>	}
269	>
270	>	strcat(errmsg, "\n");
271	>	eputs(errmsg);
272	>	}
273	>
274	>	/* =============================================================
275	>	* FLUX INTEGRATION - Get total flux emitted from sources/ports
276	>	* ============================================================= */
277	>	for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
278	>	unsigned portCnt = 0;
279	>	srcFlux [srcIdx] = 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;
282	>	do {  /* Need at least one iteration if no ports! */
283	>	emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt
284	>	: NULL;
285		photonPartition [emap.src -> so -> otype] (&emap);
286	<
287	<	if (photonRepTime) {
288	<	sprintf(errmsg, "Integrating flux from source %s (mod %s) ",
289	<	source [srcIdx].so -> oname,
290	<	objptr(source [srcIdx].so -> omod) -> oname);
297	<
286	>
287	>	if (verbose) {
288	>	sprintf(errmsg, "\tIntegrating flux from source %s ",
289	>	source [srcIdx].so -> oname);
290	>
291		if (emap.port) {
292		sprintf(errmsg2, "via port %s ",
293		photonPorts [portCnt].so -> oname);
294		strcat(errmsg, errmsg2);
295		}
296	<
297	<	sprintf(errmsg2, "(%lu partitions)...\n",
305	<	emap.numPartitions);
296	>
297	>	sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions);
298		strcat(errmsg, errmsg2);
299		eputs(errmsg);
300	+	#if NIX
301		fflush(stderr);
302	<	}
302	>	#endif
303	>	}
304
305	<	for (emap.partitionCnt = 0;
312	<	emap.partitionCnt < emap.numPartitions;
305	>	for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
306		emap.partitionCnt++) {
307		initPhotonEmission(&emap, pdfSamples);
308		srcFlux [srcIdx] += colorAvg(emap.partFlux);
309		}
310
311		portCnt++;
312	<	} while (portCnt < numPhotonPorts);
313	<
312	>	} while (portCnt < numPhotonPorts);
313	>
314		if (srcFlux [srcIdx] < FTINY) {
315		sprintf(errmsg, "source %s has zero emission",
316		source [srcIdx].so -> oname);
317		error(WARNING, errmsg);
318		}
319	<	else {
320	<	/* ==========================================================
321	<	* 2-PASS PHOTON DISTRIBUTION
322	<	* Pass 1 (pre):  emit fraction of target photon count
323	<	* Pass 2 (main): based on outcome of pass 1, estimate
324	<	*                remaining number of photons to emit to
325	<	*                approximate target count
326	<	* ========================================================== */
327	<	do {
328	<	if (!passCnt) {
329	<	/* INIT PASS 1 */
330	<	if (++prePassCnt > maxPreDistrib) {
331	<	/* Warn if no photons contributed after sufficient
332	<	* iterations */
333	<	sprintf(errmsg, "too many prepasses, no photons "
341	<	"from source %s", source [srcIdx].so -> oname);
342	<	error(WARNING, errmsg);
343	<	break;
344	<	}
319	>	}
320	>
321	>	/* Allocate & init per-subprocess primary heap files */
322	>	primaryHeap = calloc(numProc, sizeof(FILE*));
323	>	primaryHeapFname = calloc(numProc, sizeof(char*));
324	>	primaryOfs = calloc(numProc, sizeof(PhotonPrimaryIdx));
325	>	if (!primaryHeap || !primaryHeapFname || !primaryOfs)
326	>	error(SYSTEM, "failed primary heap allocation in "
327	>	"distribPhotonContrib");
328	>
329	>	for (proc = 0; proc < numProc; proc++) {
330	>	primaryHeapFname [proc] = malloc(PMAP_TMPFNLEN);
331	>	if (!primaryHeapFname [proc])
332	>	error(SYSTEM, "failed primary heap file allocation in "
333	>	"distribPhotonContrib");
334
335	<	/* Num to emit is fraction of target count */
336	<	srcNumEmit = preDistrib * srcDistribTarget;
337	<	}
335	>	mktemp(strcpy(primaryHeapFname [proc], PMAP_TMPFNAME));
336	>	if (!(primaryHeap [proc] = fopen(primaryHeapFname [proc], "w+b")))
337	>	error(SYSTEM, "failed opening primary heap file in "
338	>	"distribPhotonContrib");
339	>	}
340
341	<	else {
342	<	/* 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;
341	>	/* Record start time for progress reports */
342	>	repStartTime = time(NULL);
343
344	<	if (!srcNumEmit)
345	<	/* No photons left to distribute in main pass */
346	<	break;
347	<	}
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);
344	>	if (verbose) {
345	>	sprintf(errmsg, "\nPhoton distribution @ %d procs\n", numProc);
346	>	eputs(errmsg);
347	>	}
348
349	<	if (photonRepTime) {
350	<	if (!passCnt)
351	<	sprintf(errmsg, "PREPASS %d on source %s (mod %s) ",
352	<	prePassCnt, source [srcIdx].so -> oname,
353	<	objptr(source[srcIdx].so->omod) -> oname);
354	<	else
355	<	sprintf(errmsg, "MAIN PASS on source %s (mod %s) ",
356	<	source [srcIdx].so -> oname,
357	<	objptr(source[srcIdx].so->omod) -> oname);
358	<
359	<	if (emap.port) {
360	<	sprintf(errmsg2, "via port %s ",
361	<	photonPorts [portCnt].so -> oname);
362	<	strcat(errmsg, errmsg2);
349	>	/* MAIN LOOP */
350	>	for (proc = 0; proc < numProc; proc++) {
351	>	#if NIX
352	>	if (!(pid = fork())) {
353	>	/* SUBPROCESS ENTERS HERE; opened and mmapped files inherited */
354	>	#else
355	>	if (1) {
356	>	/* No subprocess under Windoze */
357	>	#endif
358	>	/* Local photon counters for this subprocess */
359	>	unsigned long  lastNumPhotons = 0, localNumEmitted = 0;
360	>	double         photonFluxSum = 0;   /* Accum. photon flux */
361	>
362	>	/* Seed RNGs from PID for decorellated photon distribution */
363	>	pmapSeed(randSeed + proc, partState);
364	>	pmapSeed(randSeed + (proc + 1) % numProc, emitState);
365	>	pmapSeed(randSeed + (proc + 2) % numProc, cntState);
366	>	pmapSeed(randSeed + (proc + 3) % numProc, mediumState);
367	>	pmapSeed(randSeed + (proc + 4) % numProc, scatterState);
368	>	pmapSeed(randSeed + (proc + 5) % numProc, rouletteState);
369	>
370	>	#ifdef PMAP_SIGUSR
371	>	double partNumEmit;
372	>	unsigned long partEmitCnt;
373	>	double srcPhotonFlux, avgPhotonFlux;
374	>	unsigned       portCnt, passCnt, prePassCnt;
375	>	float          srcPreDistrib;
376	>	double         srcNumEmit;     /* # to emit from source */
377	>	unsigned long  srcNumDistrib;  /* # stored */
378	>
379	>	void sigUsrDiags()
380	>	/* Loop diags via SIGUSR1 */
381	>	{
382	>	sprintf(errmsg,
383	>	"********************* Proc %d Diags *********************\n"
384	>	"srcIdx = %d (%s)\nportCnt = %d (%s)\npassCnt = %d\n"
385	>	"srcFlux = %f\nsrcPhotonFlux = %f\navgPhotonFlux = %f\n"
386	>	"partNumEmit = %f\npartEmitCnt = %lu\n\n",
387	>	proc, srcIdx, findmaterial(source [srcIdx].so) -> oname,
388	>	portCnt, photonPorts [portCnt].so -> oname,
389	>	passCnt, srcFlux [srcIdx], srcPhotonFlux, avgPhotonFlux,
390	>	partNumEmit, partEmitCnt);
391	>	eputs(errmsg);
392	>	fflush(stderr);
393	>	}
394	>	#endif
395	>
396	>	#ifdef PMAP_SIGUSR
397	>	signal(SIGUSR1, sigUsrDiags);
398	>	#endif
399	>
400	>	#ifdef DEBUG_PMAP
401	>	/* Output child process PID after random delay to prevent corrupted
402	>	* console output due to race condition */
403	>	usleep(1e6 * pmapRandom(rouletteState));
404	>	fprintf(stderr, "Proc %d: PID = %d "
405	>	"(waiting 10 sec to attach debugger...)\n",
406	>	proc, getpid());
407	>	/* Allow time for debugger to attach to child process */
408	>	sleep(10);
409	>	#endif
410	>
411	>	/* =============================================================
412	>	* 2-PASS PHOTON DISTRIBUTION
413	>	* Pass 1 (pre):  emit fraction of target photon count
414	>	* Pass 2 (main): based on outcome of pass 1, estimate remaining
415	>	*                number of photons to emit to approximate target
416	>	*                count
417	>	* ============================================================= */
418	>	for (srcIdx = 0; srcIdx < nsources; srcIdx++) {
419	>	#ifndef PMAP_SIGUSR
420	>	unsigned       portCnt, passCnt = 0, prePassCnt = 0;
421	>	float          srcPreDistrib = preDistrib;
422	>	double         srcNumEmit = 0;       /* # to emit from source */
423	>	unsigned long  srcNumDistrib = pm -> numPhotons;  /* # stored */
424	>	#else
425	>	passCnt = prePassCnt = 0;
426	>	srcPreDistrib = preDistrib;
427	>	srcNumEmit = 0;       /* # to emit from source */
428	>	srcNumDistrib = pm -> numPhotons;  /* # stored */
429	>	#endif
430	>
431	>	if (srcFlux [srcIdx] < FTINY)
432	>	continue;
433	>
434	>	while (passCnt < 2) {
435	>	if (!passCnt) {
436	>	/* INIT PASS 1 */
437	>	if (++prePassCnt > maxPreDistrib) {
438	>	/* Warn if no photons contributed after sufficient
439	>	* iterations; only output from subprocess 0 to reduce
440	>	* console clutter */
441	>	if (!proc) {
442	>	sprintf(errmsg,
443	>	"source %s: too many prepasses, skipped",
444	>	source [srcIdx].so -> oname);
445	>	error(WARNING, errmsg);
446	>	}
447	>
448	>	break;
449		}
450
451	<	sprintf(errmsg2, "(%lu partitions)...\n",
452	<	emap.numPartitions);
397	<	strcat(errmsg, errmsg2);
398	<	eputs(errmsg);
399	<	fflush(stderr);
451	>	/* Num to emit is fraction of target count */
452	>	srcNumEmit = srcPreDistrib * srcDistribTarget;
453		}
454	<
455	<	for (emap.partitionCnt = 0;
456	<	emap.partitionCnt < emap.numPartitions;
457	<	emap.partitionCnt++) {
458	<	double partNumEmit;
406	<	unsigned long partEmitCnt;
454	>	else {
455	>	/* INIT PASS 2 */
456	>	#ifndef PMAP_SIGUSR
457	>	double srcPhotonFlux, avgPhotonFlux;
458	>	#endif
459
460	<	/* Get photon origin within current source partishunn
461	<	* and build emission map */
462	<	photonOrigin [emap.src -> so -> otype] (&emap);
463	<	initPhotonEmission(&emap, pdfSamples);
464	<
465	<	/* Number of photons to emit from ziss partishunn;
466	<	* scale according to its normalised contribushunn to
467	<	* the emitted source flux */
468	<	partNumEmit = srcNumEmit * colorAvg(emap.partFlux) /
469	<	srcFlux [srcIdx];
470	<	partEmitCnt = (unsigned long)partNumEmit;
471	<
472	<	/* Probabilistically account for fractional photons */
473	<	if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
474	<	partEmitCnt++;
460	>	/* Based on the outcome of the predistribution we can now
461	>	* figure out how many more photons we have to emit from
462	>	* the current source to meet the target count,
463	>	* srcDistribTarget. This value is clamped to 0 in case
464	>	* the target has already been exceeded in pass 1.
465	>	* srcNumEmit and srcNumDistrib is the number of photons
466	>	* emitted and distributed (stored) from the current
467	>	* source in pass 1, respectively. */
468	>	srcNumDistrib = pm -> numPhotons - srcNumDistrib;
469	>	srcNumEmit *= srcNumDistrib
470	>	? max(srcDistribTarget/srcNumDistrib, 1) - 1
471	>	: 0;
472	>
473	>	if (!srcNumEmit)
474	>	/* No photons left to distribute in main pass */
475	>	break;
476
477	<	/* Integer counter avoids FP rounding errors */
478	<	while (partEmitCnt--) {
426	<	RAY photonRay;
477	>	srcPhotonFlux = srcFlux [srcIdx] / srcNumEmit;
478	>	avgPhotonFlux = photonFluxSum / (srcIdx + 1);
479
480	<	/* Emit photon according to PDF (if any), allocate
481	<	* associated primary ray, and trace through scene
482	<	* until absorbed/leaked */
483	<	emitPhoton(&emap, &photonRay);
484	<	addPhotonPrimary(pm, &photonRay);
485	<	tracePhoton(&photonRay);
480	>	if (avgPhotonFlux > FTINY &&
481	>	srcPhotonFlux / avgPhotonFlux < FTINY) {
482	>	/* Skip source if its photon flux is grossly below the
483	>	* running average, indicating negligible contributions
484	>	* at the expense of excessive distribution time; only
485	>	* output from subproc 0 to reduce console clutter */
486	>	if (!proc) {
487	>	sprintf(errmsg,
488	>	"source %s: itsy bitsy photon flux, skipped",
489	>	source [srcIdx].so -> oname);
490	>	error(WARNING, errmsg);
491	>	}
492	>
493	>	srcNumEmit = 0;   /* Or just break??? */
494	>	}
495	>
496	>	/* Update sum of photon flux per light source */
497	>	photonFluxSum += srcPhotonFlux;
498	>	}
499	>
500	>	portCnt = 0;
501	>	do {    /* Need at least one iteration if no ports! */
502	>	emap.src = source + srcIdx;
503	>	emap.port = emap.src -> sflags & SDISTANT
504	>	? photonPorts + portCnt : NULL;
505	>	photonPartition [emap.src -> so -> otype] (&emap);
506	>
507	>	if (verbose && !proc) {
508	>	/* Output from subproc 0 only to avoid race condition
509	>	* on console I/O */
510	>	if (!passCnt)
511	>	sprintf(errmsg, "\tPREPASS %d on source %s ",
512	>	prePassCnt, source [srcIdx].so -> oname);
513	>	else
514	>	sprintf(errmsg, "\tMAIN PASS on source %s ",
515	>	source [srcIdx].so -> oname);
516	>
517	>	if (emap.port) {
518	>	sprintf(errmsg2, "via port %s ",
519	>	photonPorts [portCnt].so -> oname);
520	>	strcat(errmsg, errmsg2);
521	>	}
522	>
523	>	sprintf(errmsg2, "(%lu partitions)\n",
524	>	emap.numPartitions);
525	>	strcat(errmsg, errmsg2);
526	>	eputs(errmsg);
527	>	#if NIX
528	>	fflush(stderr);
529	>	#endif
530	>	}
531	>
532	>	for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions;
533	>	emap.partitionCnt++) {
534	>	#ifndef PMAP_SIGUSR
535	>	double partNumEmit;
536	>	unsigned long partEmitCnt;
537	>	#endif
538
539	<	/* Record progress */
540	<	repProgress++;
539	>	/* Get photon origin within current source partishunn
540	>	* and build emission map */
541	>	photonOrigin [emap.src -> so -> otype] (&emap);
542	>	initPhotonEmission(&emap, pdfSamples);
543
544	<	if (photonRepTime > 0 &&
545	<	time(NULL) >= repLastTime + photonRepTime)
544	>	/* Number of photons to emit from ziss partishunn;
545	>	* scale according to its normalised contribushunn to
546	>	* the emitted source flux */
547	>	partNumEmit = srcNumEmit * colorAvg(emap.partFlux) /
548	>	srcFlux [srcIdx];
549	>	partEmitCnt = (unsigned long)partNumEmit;
550	>
551	>	/* Probabilistically account for fractional photons */
552	>	if (pmapRandom(cntState) < partNumEmit - partEmitCnt)
553	>	partEmitCnt++;
554	>
555	>	/* Update local and shared global emission counter */
556	>	photonCnt [PHOTONCNT_NUMEMIT(srcIdx)] += partEmitCnt;
557	>	localNumEmitted += partEmitCnt;
558	>
559	>	/* Integer counter avoids FP rounding errors during
560	>	* iteration */
561	>	while (partEmitCnt--) {
562	>	RAY photonRay;
563	>
564	>	/* Emit photon according to PDF (if any), allocate
565	>	* associated primary ray, and trace through scene
566	>	* until absorbed/leaked; emitPhoton() sets the
567	>	* emitting light source index in photonRay */
568	>	emitPhoton(&emap, &photonRay);
569	>	#if 1
570	>	if (emap.port)
571	>	/* !!!  PHOTON PORT REJECTION SAMPLING HACK: set
572	>	* !!!  photon port as fake hit object for
573	>	* !!!  primary ray to check for intersection in
574	>	* !!!  tracePhoton() */
575	>	photonRay.ro = emap.port -> so;
576	>	#endif
577	>	newPhotonPrimary(pm, &photonRay, primaryHeap[proc]);
578	>	/* Set subprocess index in photonRay for post-
579	>	* distrib primary index linearisation; this is
580	>	* propagated with the primary index in photonRay
581	>	* and set for photon hits by newPhoton() */
582	>	PMAP_SETRAYPROC(&photonRay, proc);
583	>	tracePhoton(&photonRay);
584	>	}
585	>
586	>	/* Update shared global photon count */
587	>	photonCnt [PHOTONCNT_NUMPHOT] += pm -> numPhotons -
588	>	lastNumPhotons;
589	>	lastNumPhotons = pm -> numPhotons;
590	>	#if !NIX
591	>	/* Synchronous progress report on Windoze */
592	>	if (!proc && photonRepTime > 0 &&
593	>	time(NULL) >= repLastTime + photonRepTime) {
594	>	unsigned s;
595	>	repComplete = pm -> distribTarget * numProc;
596	>	repProgress = photonCnt [PHOTONCNT_NUMPHOT];
597	>
598	>	for (repEmitted = 0, s = 0; s < nsources; s++)
599	>	repEmitted += photonCnt [PHOTONCNT_NUMEMIT(s)];
600	>
601		pmapDistribReport();
602	<	#ifdef SIGCONT
603	<	else signal(SIGCONT, pmapDistribReport);
443	<	#endif
602	>	}
603	>	#endif
604		}
445	–	}
446	–
447	–	portCnt++;
448	–	} while (portCnt < numPhotonPorts);
605
606	<	if (pm -> heapEnd == srcNumDistrib)
607	<	/* Double preDistrib in case no photons were stored
608	<	* for this source and redo pass 1 */
609	<	preDistrib *= 2;
610	<	else {
611	<	/* Now do pass 2 */
612	<	passCnt++;
613	<	if (photonRepTime)
614	<	eputs("\n");
606	>	portCnt++;
607	>	} while (portCnt < numPhotonPorts);
608	>
609	>	if (pm -> numPhotons == srcNumDistrib) {
610	>	/* Double predistrib factor in case no photons were stored
611	>	* for this source and redo pass 1 */
612	>	srcPreDistrib *= 2;
613	>	}
614	>	else {
615	>	/* Now do pass 2 */
616	>	passCnt++;
617	>	}
618		}
619	<	} while (passCnt < 2);
620	<
621	<	/* Flux per photon emitted from this source; repProgress is the
622	<	* number of emitted photons after both passes */
623	<	srcFlux [srcIdx] = repProgress ? srcFlux [srcIdx] / repProgress
624	<	: 0;
619	>	}
620	>
621	>	/* Flush heap buffa one final time to prevent data corruption */
622	>	flushPhotonHeap(pm);
623	>	/* Flush final photon primary to primary heap file */
624	>	newPhotonPrimary(pm, NULL, primaryHeap [proc]);
625	>	/* Heap files closed automatically on exit
626	>	fclose(pm -> heap);
627	>	fclose(primaryHeap [proc]); */
628	>
629	>	#ifdef DEBUG_PMAP
630	>	sprintf(errmsg, "Proc %d total %ld photons\n", proc,
631	>	pm -> numPhotons);
632	>	eputs(errmsg);
633	>	fflush(stderr);
634	>	#endif
635	>
636	>	#ifdef PMAP_SIGUSR
637	>	signal(SIGUSR1, SIG_DFL);
638	>	#endif
639	>
640	>	#if NIX
641	>	/* Terminate subprocess */
642	>	exit(0);
643	>	#endif
644		}
645	+	else if (pid < 0)
646	+	error(SYSTEM, "failed to fork subprocess in distribPhotonContrib");
647		}
648
649	+	#if NIX
650	+	/* PARENT PROCESS CONTINUES HERE */
651	+	#ifdef SIGCONT
652	+	/* Enable progress report signal handler */
653	+	signal(SIGCONT, pmapDistribReport);
654	+	#endif
655	+	/* Wait for subprocesses to complete while reporting progress */
656	+	proc = numProc;
657	+	while (proc) {
658	+	while (waitpid(-1, &stat, WNOHANG) > 0) {
659	+	/* Subprocess exited; check status */
660	+	if (!WIFEXITED(stat) || WEXITSTATUS(stat))
661	+	error(USER, "failed photon distribution");
662	+
663	+	--proc;
664	+	}
665	+
666	+	/* Nod off for a bit and update progress  */
667	+	sleep(1);
668	+
669	+	/* Asynchronous progress report from shared subprocess counters */
670	+	repComplete = pm -> distribTarget * numProc;
671	+	repProgress = photonCnt [PHOTONCNT_NUMPHOT];
672	+
673	+	for (repEmitted = 0, srcIdx = 0; srcIdx < nsources; srcIdx++)
674	+	repEmitted += photonCnt [PHOTONCNT_NUMEMIT(srcIdx)];
675	+
676	+	/* Get global photon count from shmem updated by subprocs */
677	+	pm -> numPhotons = photonCnt [PHOTONCNT_NUMPHOT];
678	+
679	+	if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime)
680	+	pmapDistribReport();
681	+	#ifdef SIGCONT
682	+	else signal(SIGCONT, pmapDistribReport);
683	+	#endif
684	+	}
685	+	#endif /* NIX */
686	+
687		/* ================================================================
688		* POST-DISTRIBUTION - Set photon flux and build kd-tree, etc.
689		* ================================================================ */
690	<	#ifdef SIGCONT
691	<	signal(SIGCONT, SIG_DFL);
692	<	#endif
690	>	#ifdef SIGCONT
691	>	/* Reset signal handler */
692	>	signal(SIGCONT, SIG_DFL);
693	>	#endif
694		free(emap.samples);
695
696	<	if (!pm -> heapEnd)
697	<	error(USER, "empty photon map");
696	>	if (!pm -> numPhotons)
697	>	error(USER, "empty contribution photon map");
698
699	<	/* Check for valid primary photon rays */
700	<	if (!pm -> primary)
699	>	/* Load per-subprocess primary rays into pm -> primary array */
700	>	/* Dumb compilers apparently need the char** cast */
701	>	pm -> numPrimary = buildPrimaries(pm, primaryHeap,
702	>	(char**)primaryHeapFname,
703	>	primaryOfs, numProc);
704	>	if (!pm -> numPrimary)
705		error(INTERNAL, "no primary rays in contribution photon map");
483	–
484	–	if (pm -> primary [pm -> primaryEnd].srcIdx < 0)
485	–	/* Last primary ray is unused, so decrement counter */
486	–	pm -> primaryEnd--;
706
707	<	if (photonRepTime) {
708	<	eputs("\nBuilding contrib photon heap...\n");
707	>	/* Set photon flux per source */
708	>	for (srcIdx = 0; srcIdx < nsources; srcIdx++)
709	>	srcFlux [srcIdx] /= photonCnt [PHOTONCNT_NUMEMIT(srcIdx)];
710	>	#if NIX
711	>	/* Photon counters no longer needed, unmap shared memory */
712	>	munmap(photonCnt, sizeof(*photonCnt));
713	>	close(shmFile);
714	>	unlink(shmFname);
715	>	#else
716	>	free(photonCnt);
717	>	#endif
718	>
719	>	if (verbose) {
720	>	eputs("\nBuilding contribution photon map...\n");
721	>	#if NIX
722		fflush(stderr);
723	+	#endif
724		}
725	+
726	+	/* Build underlying data structure; heap is destroyed */
727	+	buildPhotonMap(pm, srcFlux, primaryOfs, numProc);
728	+
729	+	/* Free per-subprocess primary heap files */
730	+	for (proc = 0; proc < numProc; proc++)
731	+	free(primaryHeapFname [proc]);
732
733	<	balancePhotons(pm, srcFlux);
733	>	free(primaryHeapFname);
734	>	free(primaryHeap);
735	>	free(primaryOfs);
736	>
737	>	if (verbose)
738	>	eputs("\n");
739		}

Diff Legend

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