ViewVC Help

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

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

Comparing ray/src/rt/pmapdata.c (file contents):
Revision 2.11 by greg, Tue Aug 18 18:45:55 2015 UTC vs.
Revision 2.15 by rschregle, Tue May 17 17:39:47 2016 UTC

#	Line 1 | Line 1
1		#ifndef lint
2		static const char RCSid[] = "$Id$";
3		#endif
4	+
5		/*
6	<	==================================================================
7	<	Photon map data structures and kd-tree handling
6	>	=========================================================================
7	>	Photon map types and interface to nearest neighbour lookups in underlying
8	>	point cloud data structure.
9	>
10	>	The default data structure is an in-core kd-tree (see pmapkdt.{h,c}).
11	>	This can be overriden with the PMAP_OOC compiletime switch, which enables
12	>	an out-of-core octree (see oococt.{h,c}).
13
14		Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
15		(c) Fraunhofer Institute for Solar Energy Systems,
16		(c) Lucerne University of Applied Sciences and Arts,
17	<	supported by the Swiss National Science Foundation (SNSF, #147053)
18	<	==================================================================
17	>	supported by the Swiss National Science Foundation (SNSF, #147053)
18	>	==========================================================================
19
20		$Id$
21		*/
#	Line 32 | Line 38 | PhotonMap *photonMaps [NUM_PMAP_TYPES] = {
38
39
40
41	+	/* Include routines to handle underlying point cloud data structure */
42	+	#ifdef PMAP_OOC
43	+	#include "pmapooc.c"
44	+	#else
45	+	#include "pmapkdt.c"
46	+	#endif
47	+
48	+
49	+
50		void initPhotonMap (PhotonMap *pmap, PhotonMapType t)
51		/* Init photon map 'n' stuff... */
52		{
53		if (!pmap)
54		return;
55
56	<	pmap -> heapSize = pmap -> heapEnd = 0;
42	<	pmap -> heap = NULL;
43	<	pmap -> squeue = NULL;
56	>	pmap -> numPhotons = 0;
57		pmap -> biasCompHist = NULL;
58		pmap -> maxPos [0] = pmap -> maxPos [1] = pmap -> maxPos [2] = -FHUGE;
59		pmap -> minPos [0] = pmap -> minPos [1] = pmap -> minPos [2] = FHUGE;
#	Line 50 | Line 63 | void initPhotonMap (PhotonMap *pmap, PhotonMapType t)
63		pmap -> numDensity = 0;
64		pmap -> distribRatio = 1;
65		pmap -> type = t;
66	+	pmap -> squeue.node = NULL;
67	+	pmap -> squeue.len = 0;
68
69		/* Init local RNG state */
70		pmap -> randState [0] = 10243;
#	Line 61 | Line 76 | void initPhotonMap (PhotonMap *pmap, PhotonMapType t)
76		/* Set up type-specific photon lookup callback */
77		pmap -> lookup = pmapLookup [t];
78
79	<	pmap -> primary = NULL;
80	<	pmap -> primarySize = pmap -> primaryEnd = 0;
79	>	/* Mark primary photon ray as unused */
80	>	pmap -> lastPrimary.srcIdx = -1;
81	>	pmap -> numPrimary = 0;
82	>	pmap -> primaries = NULL;
83	>
84	>	/* Init storage */
85	>	pmap -> heap = NULL;
86	>	pmap -> heapBuf = NULL;
87	>	pmap -> heapBufLen = 0;
88	>	#ifdef PMAP_OOC
89	>	OOC_Null(&pmap -> store);
90	>	#else
91	>	kdT_Null(&pmap -> store);
92	>	#endif
93		}
94
95
96
97	<	const PhotonPrimary* addPhotonPrimary (PhotonMap *pmap, const RAY *ray)
97	>	void initPhotonHeap (PhotonMap *pmap)
98		{
99	<	PhotonPrimary *prim = NULL;
73	<	FVECT dvec;
99	>	int fdFlags;
100
101	<	if (!pmap || !ray)
102	<	return NULL;
101	>	if (!pmap)
102	>	error(INTERNAL, "undefined photon map in initPhotonHeap");
103
104	<	/* Check if last primary ray has spawned photons (srcIdx >= 0, see
105	<	* addPhoton()), in which case we keep it and allocate a new one;
106	<	* otherwise we overwrite the unused entry */
107	<	if (pmap -> primary && pmap -> primary [pmap -> primaryEnd].srcIdx >= 0)
108	<	pmap -> primaryEnd++;
109	<
110	<	if (!pmap -> primarySize || pmap -> primaryEnd >= pmap -> primarySize) {
85	<	/* Allocate/enlarge array */
86	<	pmap -> primarySize += pmap -> heapSizeInc;
87	<
88	<	/* Counter wraparound? */
89	<	if (pmap -> primarySize < pmap -> heapSizeInc)
90	<	error(INTERNAL, "photon primary overflow");
91	<
92	<	pmap -> primary = (PhotonPrimary *)realloc(pmap -> primary,
93	<	sizeof(PhotonPrimary) *
94	<	pmap -> primarySize);
95	<	if (!pmap -> primary)
96	<	error(USER, "can't allocate photon primaries");
104	>	if (!pmap -> heap) {
105	>	/* Open heap file */
106	>	if (!(pmap -> heap = tmpfile()))
107	>	error(SYSTEM, "failed opening heap file in initPhotonHeap");
108	>	fdFlags = fcntl(fileno(pmap -> heap), F_GETFL);
109	>	fcntl(fileno(pmap -> heap), F_SETFL, fdFlags | O_APPEND);
110	>	/*      ftruncate(fileno(pmap -> heap), 0); */
111		}
112	+	}
113	+
114	+
115	+
116	+	void flushPhotonHeap (PhotonMap *pmap)
117	+	{
118	+	int                  fd;
119	+	const unsigned long  len = pmap -> heapBufLen * sizeof(Photon);
120
121	<	prim = pmap -> primary + pmap -> primaryEnd;
121	>	if (!pmap)
122	>	error(INTERNAL, "undefined photon map in flushPhotonHeap");
123	>
124	>	if (!pmap -> heap || !pmap -> heapBuf)
125	>	error(INTERNAL, "undefined heap in flushPhotonHeap");
126	>
127	>	/* Atomically seek and write block to heap */
128	>	/* !!! Unbuffered I/O via pwrite() avoids potential race conditions
129	>	* !!! and buffer corruption which can occur with lseek()/fseek()
130	>	* !!! followed by write()/fwrite(). */
131	>	fd = fileno(pmap -> heap);
132
133	<	/* Mark unused with negative source index until path spawns a photon (see
134	<	* addPhoton()) */
135	<	prim -> srcIdx = -1;
133	>	#ifdef DEBUG_PMAP
134	>	sprintf(errmsg, "Proc %d: flushing %ld photons from pos %ld\n", getpid(),
135	>	pmap -> heapBufLen, lseek(fd, 0, SEEK_END) / sizeof(Photon));
136	>	eputs(errmsg);
137	>	#endif
138	>
139	>	/*if (pwrite(fd, pmap -> heapBuf, len, lseek(fd, 0, SEEK_END)) != len) */
140	>	if (write(fd, pmap -> heapBuf, len) != len)
141	>	error(SYSTEM, "failed append to heap file in flushPhotonHeap");
142	>
143	>	if (fsync(fd))
144	>	error(SYSTEM, "failed fsync in flushPhotonHeap");
145
146	<	/* Reverse incident direction to point to light source */
147	<	dvec [0] = -ray -> rdir [0];
107	<	dvec [1] = -ray -> rdir [1];
108	<	dvec [2] = -ray -> rdir [2];
109	<	prim -> dir = encodedir(dvec);
146	>	pmap -> heapBufLen = 0;
147	>	}
148
149	<	VCOPY(prim -> pos, ray -> rop);
149	>
150	>
151	>	#ifdef DEBUG_OOC
152	>	static int checkPhotonHeap (FILE *file)
153	>	/* Check heap for nonsensical or duplicate photons */
154	>	{
155	>	Photon   p, lastp;
156	>	int      i, dup;
157
158	<	return prim;
158	>	rewind(file);
159	>	memset(&lastp, 0, sizeof(lastp));
160	>
161	>	while (fread(&p, sizeof(p), 1, file)) {
162	>	dup = 1;
163	>
164	>	for (i = 0; i <= 2; i++) {
165	>	if (p.pos [i] < thescene.cuorg [i] ||
166	>	p.pos [i] > thescene.cuorg [i] + thescene.cusize) {
167	>
168	>	sprintf(errmsg, "corrupt photon in heap at [%f, %f, %f]\n",
169	>	p.pos [0], p.pos [1], p.pos [2]);
170	>	error(WARNING, errmsg);
171	>	}
172	>
173	>	dup &= p.pos [i] == lastp.pos [i];
174	>	}
175	>
176	>	if (dup) {
177	>	sprintf(errmsg,
178	>	"consecutive duplicate photon in heap at [%f, %f, %f]\n",
179	>	p.pos [0], p.pos [1], p.pos [2]);
180	>	error(WARNING, errmsg);
181	>	}
182	>	}
183	>
184	>	return 0;
185		}
186	+	#endif
187
188
189
190	<	const Photon* addPhoton (PhotonMap* pmap, const RAY* ray)
190	>	int newPhoton (PhotonMap* pmap, const RAY* ray)
191		{
192		unsigned i;
193	<	Photon* photon = NULL;
193	>	Photon photon;
194		COLOR photonFlux;
195
196		/* Account for distribution ratio */
197		if (!pmap || pmapRandom(pmap -> randState) > pmap -> distribRatio)
198	<	return NULL;
198	>	return -1;
199
200		/* Don't store on sources */
201		if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype))
202	<	return NULL;
202	>	return -1;
203
132	–	#if 0
133	–	if (inContribPmap(pmap)) {
134	–	/* Adding contribution photon */
135	–	if (ray -> parent && ray -> parent -> rtype & PRIMARY)
136	–	/* Add primary photon ray if parent is primary; by putting this
137	–	* here and checking the ray's immediate parent, we only add
138	–	* primaries that actually contribute photons, and we only add them
139	–	* once.  */
140	–	addPhotonPrimary(pmap, ray -> parent);
141	–
142	–	/* Save index to primary ray (remains unchanged if primary already in
143	–	* array) */
144	–	primary = pmap -> primaryEnd;
145	–	}
146	–	#endif
147	–
204		#ifdef PMAP_ROI
205	<	/* Store photon if within region of interest -- for egg-spurtz only! */
205	>	/* Store photon if within region of interest -- for Ze Eckspertz only! */
206		if (ray -> rop [0] >= pmapROI [0] && ray -> rop [0] <= pmapROI [1] &&
207		ray -> rop [1] >= pmapROI [2] && ray -> rop [1] <= pmapROI [3] &&
208		ray -> rop [2] >= pmapROI [4] && ray -> rop [2] <= pmapROI [5])
209		#endif
210		{
155	–	if (pmap -> heapEnd >= pmap -> heapSize) {
156	–	/* Enlarge heap */
157	–	pmap -> heapSize += pmap -> heapSizeInc;
158	–
159	–	/* Counter wraparound? */
160	–	if (pmap -> heapSize < pmap -> heapSizeInc)
161	–	error(INTERNAL, "photon heap overflow");
162	–
163	–	pmap -> heap = (Photon *)realloc(pmap -> heap,
164	–	sizeof(Photon) * pmap -> heapSize);
165	–	if (!pmap -> heap)
166	–	error(USER, "can't allocate photon heap");
167	–	}
168	–
169	–	photon = pmap -> heap + pmap -> heapEnd++;
170	–
211		/* Adjust flux according to distribution ratio and ray weight */
212		copycolor(photonFlux, ray -> rcol);
213		scalecolor(photonFlux,
214		ray -> rweight / (pmap -> distribRatio ? pmap -> distribRatio
215		: 1));
216	<	setPhotonFlux(photon, photonFlux);
216	>	setPhotonFlux(&photon, photonFlux);
217
218		/* Set photon position and flags */
219	<	VCOPY(photon -> pos, ray -> rop);
220	<	photon -> flags = PMAP_CAUSTICRAY(ray) ? PMAP_CAUST_BIT : 0;
219	>	VCOPY(photon.pos, ray -> rop);
220	>	photon.flags = 0;
221	>	photon.caustic = PMAP_CAUSTICRAY(ray);
222
223	<	/* Set primary ray index and mark as used for contrib photons */
223	>	/* Set contrib photon's primary ray and subprocess index (the latter
224	>	* to linearise the primary ray indices after photon distribution is
225	>	* complete). Also set primary ray's source index, thereby marking it
226	>	* as used. */
227		if (isContribPmap(pmap)) {
228	<	photon -> primary = pmap -> primaryEnd;
229	<	pmap -> primary [pmap -> primaryEnd].srcIdx = ray -> rsrc;
228	>	photon.primary = pmap -> numPrimary;
229	>	photon.proc = PMAP_GETRAYPROC(ray);
230	>	pmap -> lastPrimary.srcIdx = ray -> rsrc;
231		}
232	<	else photon -> primary = 0;
232	>	else photon.primary = 0;
233
234	<	/* Update min and max positions & set normal */
235	<	for (i = 0; i <= 2; i++) {
236	<	if (photon -> pos [i] < pmap -> minPos [i])
237	<	pmap -> minPos [i] = photon -> pos [i];
238	<	if (photon -> pos [i] > pmap -> maxPos [i])
239	<	pmap -> maxPos [i] = photon -> pos [i];
240	<	photon -> norm [i] = 127.0 * (isVolumePmap(pmap) ? ray -> rdir [i]
241	<	: ray -> ron [i]);
234	>	/* Set normal */
235	>	for (i = 0; i <= 2; i++)
236	>	photon.norm [i] = 127.0 * (isVolumePmap(pmap) ? ray -> rdir [i]
237	>	: ray -> ron [i]);
238	>
239	>	if (!pmap -> heapBuf) {
240	>	/* Lazily allocate heap buffa */
241	>	#if 1
242	>	/* Randomise buffa size to temporally decorellate buffa flushes */
243	>	srandom(randSeed + getpid());
244	>	pmap -> heapBufSize = PMAP_HEAPBUFSIZE * (0.5 + frandom());
245	>	#else
246	>	/* Randomisation disabled for reproducability during debugging */
247	>	pmap -> heapBufSize = PMAP_HEAPBUFSIZE;
248	>	#endif
249	>	if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon))))
250	>	error(SYSTEM, "failed heap buffer allocation in newPhoton");
251	>	pmap -> heapBufLen = 0;
252		}
253	+
254	+	/* Photon initialised; now append to heap buffa */
255	+	memcpy(pmap -> heapBuf + pmap -> heapBufLen, &photon, sizeof(Photon));
256	+
257	+	if (++pmap -> heapBufLen >= pmap -> heapBufSize)
258	+	/* Heap buffa full, flush to heap file */
259	+	flushPhotonHeap(pmap);
260	+
261	+	pmap -> numPhotons++;
262		}
263
264	<	return photon;
264	>	return 0;
265		}
266
267
268
269	<	static void nearestNeighbours (PhotonMap* pmap, const float pos [3],
270	<	const float norm [3], unsigned long node)
207	<	/* Recursive part of findPhotons(..).
208	<	Note that all heap and priority queue index handling is 1-based, but
209	<	accesses to the arrays are 0-based! */
269	>	void buildPhotonMap (PhotonMap *pmap, double *photonFlux,
270	>	PhotonPrimaryIdx *primaryOfs, unsigned nproc)
271		{
272	<	Photon* p = &pmap -> heap [node - 1];
273	<	unsigned i, j;
274	<	/* Signed distance to current photon's splitting plane */
275	<	float d = pos [photonDiscr(*p)] - p -> pos [photonDiscr(*p)],
276	<	d2 = d * d;
277	<	PhotonSQNode* sq = pmap -> squeue;
278	<	const unsigned sqSize = pmap -> squeueSize;
279	<	float dv [3];
272	>	unsigned long  n, nCheck = 0;
273	>	unsigned       i;
274	>	Photon         *p;
275	>	COLOR          flux;
276	>	FILE           *nuHeap;
277	>	/* Need double here to reduce summation errors */
278	>	double         avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0;
279	>	FVECT          d;
280
281	<	/* Search subtree closer to pos first; exclude other subtree if the
282	<	distance to the splitting plane is greater than maxDist */
222	<	if (d < 0) {
223	<	if (node << 1 <= pmap -> heapSize)
224	<	nearestNeighbours(pmap, pos, norm, node << 1);
225	<	if (d2 < pmap -> maxDist && node << 1 < pmap -> heapSize)
226	<	nearestNeighbours(pmap, pos, norm, (node << 1) + 1);
227	<	}
228	<	else {
229	<	if (node << 1 < pmap -> heapSize)
230	<	nearestNeighbours(pmap, pos, norm, (node << 1) + 1);
231	<	if (d2 < pmap -> maxDist && node << 1 <= pmap -> heapSize)
232	<	nearestNeighbours(pmap, pos, norm, node << 1);
233	<	}
234	<
235	<	/* Reject photon if normal faces away (ignored for volume photons) */
236	<	if (norm && DOT(norm, p -> norm) <= 0.5 * frandom())
237	<	return;
281	>	if (!pmap)
282	>	error(INTERNAL, "undefined photon map in buildPhotonMap");
283
284	<	if (isContribPmap(pmap) && pmap -> srcContrib) {
285	<	/* Lookup in contribution photon map */
286	<	OBJREC *srcMod;
287	<	const int srcIdx = photonSrcIdx(pmap, p);
288	<
289	<	if (srcIdx < 0 || srcIdx >= nsources)
245	<	error(INTERNAL, "invalid light source index in photon map");
246	<
247	<	srcMod = findmaterial(source [srcIdx].so);
284	>	/* Get number of photons from heapfile size */
285	>	fseek(pmap -> heap, 0, SEEK_END);
286	>	pmap -> numPhotons = ftell(pmap -> heap) / sizeof(Photon);
287	>
288	>	if (!pmap -> numPhotons)
289	>	error(INTERNAL, "empty photon map in buildPhotonMap");
290
291	<	/* Reject photon if contributions from light source which emitted it
292	<	* are not sought */
251	<	if (!lu_find(pmap -> srcContrib, srcMod -> oname) -> data)
252	<	return;
291	>	if (!pmap -> heap)
292	>	error(INTERNAL, "no heap in buildPhotonMap");
293
294	<	/* Reject non-caustic photon if lookup for caustic contribs */
295	<	if (pmap -> lookupFlags & PMAP_CAUST_BIT & ~p -> flags)
296	<	return;
294	>	#ifdef DEBUG_PMAP
295	>	eputs("Checking photon heap consistency...\n");
296	>	checkPhotonHeap(pmap -> heap);
297	>
298	>	sprintf(errmsg, "Heap contains %ld photons\n", pmap -> numPhotons);
299	>	eputs(errmsg);
300	>	#endif
301	>
302	>	/* Allocate heap buffa */
303	>	if (!pmap -> heapBuf) {
304	>	pmap -> heapBufSize = PMAP_HEAPBUFSIZE;
305	>	pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon));
306	>	if (!pmap -> heapBuf)
307	>	error(SYSTEM, "failed to allocate postprocessed photon heap in"
308	>	"buildPhotonMap");
309		}
310	+
311	+	/* We REALLY don't need yet another @%&*! heap just to hold the scaled
312	+	* photons, but can't think of a quicker fix... */
313	+	if (!(nuHeap = tmpfile()))
314	+	error(SYSTEM, "failed to open postprocessed photon heap in "
315	+	"buildPhotonMap");
316	+
317	+	rewind(pmap -> heap);
318	+
319	+	#ifdef DEBUG_PMAP
320	+	eputs("Postprocessing photons...\n");
321	+	#endif
322
323	<	/* Squared distance to current photon */
324	<	dv [0] = pos [0] - p -> pos [0];
325	<	dv [1] = pos [1] - p -> pos [1];
326	<	dv [2] = pos [2] - p -> pos [2];
327	<	d2 = DOT(dv, dv);
323	>	while (!feof(pmap -> heap)) {
324	>	pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon),
325	>	PMAP_HEAPBUFSIZE, pmap -> heap);
326	>
327	>	if (pmap -> heapBufLen) {
328	>	for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) {
329	>	/* Update min and max pos and set photon flux */
330	>	for (i = 0; i <= 2; i++) {
331	>	if (p -> pos [i] < pmap -> minPos [i])
332	>	pmap -> minPos [i] = p -> pos [i];
333	>	else if (p -> pos [i] > pmap -> maxPos [i])
334	>	pmap -> maxPos [i] = p -> pos [i];
335
336	<	/* Accept photon if closer than current max dist & add to priority queue */
337	<	if (d2 < pmap -> maxDist) {
338	<	if (pmap -> squeueEnd < sqSize) {
339	<	/* Priority queue not full; append photon and restore heap */
340	<	i = ++pmap -> squeueEnd;
341	<
342	<	while (i > 1 && sq [(i >> 1) - 1].dist <= d2) {
343	<	sq [i - 1].photon = sq [(i >> 1) - 1].photon;
344	<	sq [i - 1].dist = sq [(i >> 1) - 1].dist;
345	<	i >>= 1;
336	>	/* Update centre of gravity with photon position */
337	>	CoG [i] += p -> pos [i];
338	>	}
339	>
340	>	if (primaryOfs)
341	>	/* Linearise photon primary index from subprocess index using the
342	>	* per-subprocess offsets in primaryOfs */
343	>	p -> primary += primaryOfs [p -> proc];
344	>
345	>	/* Scale photon's flux (hitherto normalised to 1 over RGB); in
346	>	* case of a contrib photon map, this is done per light source,
347	>	* and photonFlux is assumed to be an array */
348	>	getPhotonFlux(p, flux);
349	>
350	>	if (photonFlux) {
351	>	scalecolor(flux, photonFlux [isContribPmap(pmap) ?
352	>	photonSrcIdx(pmap, p) : 0]);
353	>	setPhotonFlux(p, flux);
354	>	}
355	>
356	>	/* Update average photon flux; need a double here */
357	>	addcolor(avgFlux, flux);
358		}
359
360	<	sq [--i].photon = p;
361	<	sq [i].dist = d2;
362	<	/* Update maxDist if we've just filled the queue */
363	<	if (pmap -> squeueEnd >= pmap -> squeueSize)
364	<	pmap -> maxDist = sq [0].dist;
360	>	/* Write modified photons to new heap */
361	>	fwrite(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufLen, nuHeap);
362	>
363	>	if (ferror(nuHeap))
364	>	error(SYSTEM, "failed postprocessing photon flux in "
365	>	"buildPhotonMap");
366		}
367	<	else {
368	<	/* Priority queue full; replace maximum, restore heap, and
285	<	update maxDist */
286	<	i = 1;
287	<
288	<	while (i <= sqSize >> 1) {
289	<	j = i << 1;
290	<	if (j < sqSize && sq [j - 1].dist < sq [j].dist)
291	<	j++;
292	<	if (d2 >= sq [j - 1].dist)
293	<	break;
294	<	sq [i - 1].photon = sq [j - 1].photon;
295	<	sq [i - 1].dist = sq [j - 1].dist;
296	<	i = j;
297	<	}
298	<
299	<	sq [--i].photon = p;
300	<	sq [i].dist = d2;
301	<	pmap -> maxDist = sq [0].dist;
302	<	}
367	>
368	>	nCheck += pmap -> heapBufLen;
369		}
370	+
371	+	#ifdef DEBUG_PMAP
372	+	if (nCheck < pmap -> numPhotons)
373	+	error(INTERNAL, "truncated photon heap in buildPhotonMap");
374	+	#endif
375	+
376	+	/* Finalise average photon flux */
377	+	scalecolor(avgFlux, 1.0 / pmap -> numPhotons);
378	+	copycolor(pmap -> photonFlux, avgFlux);
379	+
380	+	/* Average photon positions to get centre of gravity */
381	+	for (i = 0; i < 3; i++)
382	+	pmap -> CoG [i] = CoG [i] /= pmap -> numPhotons;
383	+
384	+	rewind(pmap -> heap);
385	+
386	+	/* Compute average photon distance to centre of gravity */
387	+	while (!feof(pmap -> heap)) {
388	+	pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon),
389	+	PMAP_HEAPBUFSIZE, pmap -> heap);
390	+
391	+	if (pmap -> heapBufLen)
392	+	for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) {
393	+	VSUB(d, p -> pos, CoG);
394	+	CoGdist += DOT(d, d);
395	+	}
396	+	}
397	+
398	+	pmap -> CoGdist = CoGdist /= pmap -> numPhotons;
399	+
400	+	/* Swap heaps */
401	+	fclose(pmap -> heap);
402	+	pmap -> heap = nuHeap;
403	+
404	+	#ifdef PMAP_OOC
405	+	OOC_BuildPhotonMap(pmap, nproc);
406	+	#else
407	+	/* kd-tree not parallelised */
408	+	kdT_BuildPhotonMap(pmap);
409	+	#endif
410	+
411	+	/* Trash heap and its buffa */
412	+	free(pmap -> heapBuf);
413	+	fclose(pmap -> heap);
414	+	pmap -> heap = NULL;
415	+	pmap -> heapBuf = NULL;
416		}
417
418
#	Line 318 | Line 430 | static void nearestNeighbours (PhotonMap* pmap, const
430		/* Coefficient for adaptive maximum search radius */
431		#define PMAP_MAXDIST_COEFF 100
432
321	–
433		void findPhotons (PhotonMap* pmap, const RAY* ray)
434		{
324	–	float pos [3], norm [3];
435		int redo = 0;
436
437	<	if (!pmap -> squeue) {
437	>	if (!pmap -> squeue.len) {
438		/* Lazy init priority queue */
439	<	pmap -> squeueSize = pmap -> maxGather + 1;
440	<	pmap -> squeue = (PhotonSQNode*)malloc(pmap -> squeueSize *
441	<	sizeof(PhotonSQNode));
442	<	if (!pmap -> squeue)
443	<	error(USER, "can't allocate photon priority queue");
334	<
439	>	#ifdef PMAP_OOC
440	>	OOC_InitFindPhotons(pmap);
441	>	#else
442	>	kdT_InitFindPhotons(pmap);
443	>	#endif
444		pmap -> minGathered = pmap -> maxGather;
445		pmap -> maxGathered = pmap -> minGather;
446		pmap -> totalGathered = 0;
#	Line 342 | Line 451 | void findPhotons (PhotonMap* pmap, const RAY* ray)
451		pmap -> rmsError = 0;
452		/* SQUARED max search radius limit is based on avg photon distance to
453		* centre of gravity, unless fixed by user (maxDistFix > 0) */
454	<	pmap -> maxDist0 = pmap -> maxDistLimit =
454	>	pmap -> maxDist0 = pmap -> maxDist2Limit =
455		maxDistFix > 0 ? maxDistFix * maxDistFix
456	<	: PMAP_MAXDIST_COEFF * pmap -> squeueSize *
457	<	pmap -> CoGdist / pmap -> heapSize;
456	>	: PMAP_MAXDIST_COEFF * pmap -> squeue.len *
457	>	pmap -> CoGdist / pmap -> numPhotons;
458		}
459
460		do {
461	<	pmap -> squeueEnd = 0;
462	<	pmap -> maxDist = pmap -> maxDist0;
461	>	pmap -> squeue.tail = 0;
462	>	pmap -> maxDist2 = pmap -> maxDist0;
463
464		/* Search position is ray -> rorg for volume photons, since we have no
465		intersection point. Normals are ignored -- these are incident
466		directions). */
467		if (isVolumePmap(pmap)) {
468	<	VCOPY(pos, ray -> rorg);
469	<	nearestNeighbours(pmap, pos, NULL, 1);
468	>	#ifdef PMAP_OOC
469	>	OOC_FindPhotons(pmap, ray -> rorg, NULL);
470	>	#else
471	>	kdT_FindPhotons(pmap, ray -> rorg, NULL);
472	>	#endif
473		}
474		else {
475	<	VCOPY(pos, ray -> rop);
476	<	VCOPY(norm, ray -> ron);
477	<	nearestNeighbours(pmap, pos, norm, 1);
475	>	#ifdef PMAP_OOC
476	>	OOC_FindPhotons(pmap, ray -> rop, ray -> ron);
477	>	#else
478	>	kdT_FindPhotons(pmap, ray -> rop, ray -> ron);
479	>	#endif
480		}
481
482	<	if (pmap -> maxDist < FTINY) {
483	<	sprintf(errmsg, "itsy bitsy teeny weeny photon search radius %e",
484	<	sqrt(pmap -> maxDist));
485	<	error(WARNING, errmsg);
486	<	}
487	<
488	<	if (pmap -> squeueEnd < pmap -> squeueSize * pmap -> gatherTolerance) {
482	>	#ifdef PMAP_LOOKUP_INFO
483	>	fprintf(stderr, "%d/%d %s photons found within radius %.3f "
484	>	"at (%.2f,%.2f,%.2f) on %s\n", pmap -> squeue.tail,
485	>	pmap -> squeue.len, pmapName [pmap -> type], sqrt(pmap -> maxDist2),
486	>	ray -> rop [0], ray -> rop [1], ray -> rop [2],
487	>	ray -> ro ? ray -> ro -> oname : "<null>");
488	>	#endif
489	>
490	>	if (pmap -> squeue.tail < pmap -> squeue.len * pmap -> gatherTolerance) {
491		/* Short lookup; too few photons found */
492	<	if (pmap -> squeueEnd > PMAP_SHORT_LOOKUP_THRESH) {
492	>	if (pmap -> squeue.tail > PMAP_SHORT_LOOKUP_THRESH) {
493		/* Ignore short lookups which return fewer than
494		* PMAP_SHORT_LOOKUP_THRESH photons under the assumption there
495		* really are no photons in the vicinity, and increasing the max
#	Line 381 | Line 497 | void findPhotons (PhotonMap* pmap, const RAY* ray)
497		#ifdef PMAP_LOOKUP_WARN
498		sprintf(errmsg,
499		"%d/%d %s photons found at (%.2f,%.2f,%.2f) on %s",
500	<	pmap -> squeueEnd, pmap -> squeueSize,
501	<	pmapName [pmap -> type], pos [0], pos [1], pos [2],
500	>	pmap -> squeue.tail, pmap -> squeue.len,
501	>	pmapName [pmap -> type],
502	>	ray -> rop [0], ray -> rop [1], ray -> rop [2],
503		ray -> ro ? ray -> ro -> oname : "<null>");
504		error(WARNING, errmsg);
505		#endif
#	Line 391 | Line 508 | void findPhotons (PhotonMap* pmap, const RAY* ray)
508		if (maxDistFix > 0)
509		return;
510
511	<	if (pmap -> maxDist0 < pmap -> maxDistLimit) {
511	>	if (pmap -> maxDist0 < pmap -> maxDist2Limit) {
512		/* Increase max search radius if below limit & redo search */
513		pmap -> maxDist0 *= PMAP_MAXDIST_INC;
514		#ifdef PMAP_LOOKUP_REDO
#	Line 401 | Line 518 | void findPhotons (PhotonMap* pmap, const RAY* ray)
518		sprintf(errmsg,
519		redo ? "restarting photon lookup with max radius %.1e"
520		: "max photon lookup radius adjusted to %.1e",
521	<	sqrt(pmap -> maxDist0));
521	>	pmap -> maxDist0);
522		error(WARNING, errmsg);
523		#endif
524		}
525		#ifdef PMAP_LOOKUP_REDO
526		else {
527		sprintf(errmsg, "max photon lookup radius clamped to %.1e",
528	<	sqrt(pmap -> maxDist0));
528	>	pmap -> maxDist0);
529		error(WARNING, errmsg);
530		}
531		#endif
#	Line 436 | Line 553 | void findPhotons (PhotonMap* pmap, const RAY* ray)
553
554
555
556	<	static void nearest1Neighbour (PhotonMap *pmap, const float pos [3],
440	<	const float norm [3], Photon **photon,
441	<	unsigned long node)
442	<	/* Recursive part of find1Photon(..).
443	<	Note that all heap index handling is 1-based, but accesses to the
444	<	arrays are 0-based! */
556	>	void find1Photon (PhotonMap *pmap, const RAY* ray, Photon *photon)
557		{
558	<	Photon *p = pmap -> heap + node - 1;
447	<	/* Signed distance to current photon's splitting plane */
448	<	float d = pos [photonDiscr(*p)] - p -> pos [photonDiscr(*p)],
449	<	d2 = d * d;
450	<	float dv [3];
558	>	pmap -> maxDist2 = thescene.cusize;  /* ? */
559
560	<	/* Search subtree closer to pos first; exclude other subtree if the
561	<	distance to the splitting plane is greater than maxDist */
562	<	if (d < 0) {
563	<	if (node << 1 <= pmap -> heapSize)
564	<	nearest1Neighbour(pmap, pos, norm, photon, node << 1);
457	<	if (d2 < pmap -> maxDist && node << 1 < pmap -> heapSize)
458	<	nearest1Neighbour(pmap, pos, norm, photon, (node << 1) + 1);
459	<	}
460	<	else {
461	<	if (node << 1 < pmap -> heapSize)
462	<	nearest1Neighbour(pmap, pos, norm, photon, (node << 1) + 1);
463	<	if (d2 < pmap -> maxDist && node << 1 <= pmap -> heapSize)
464	<	nearest1Neighbour(pmap, pos, norm, photon, node << 1);
465	<	}
466	<
467	<	/* Squared distance to current photon */
468	<	dv [0] = pos [0] - p -> pos [0];
469	<	dv [1] = pos [1] - p -> pos [1];
470	<	dv [2] = pos [2] - p -> pos [2];
471	<	d2 = DOT(dv, dv);
472	<
473	<	if (d2 < pmap -> maxDist && DOT(norm, p -> norm) > 0.5 * frandom()) {
474	<	/* Closest photon so far with similar normal */
475	<	pmap -> maxDist = d2;
476	<	*photon = p;
477	<	}
560	>	#ifdef PMAP_OOC
561	>	OOC_Find1Photon(pmap, ray -> rop, ray -> ron, photon);
562	>	#else
563	>	kdT_Find1Photon(pmap, ray -> rop, ray -> ron, photon);
564	>	#endif
565		}
566
567
568
569	<	Photon* find1Photon (PhotonMap *pmap, const RAY* ray)
569	>	void getPhoton (PhotonMap *pmap, PhotonIdx idx, Photon *photon)
570		{
571	<	float fpos [3], norm [3];
572	<	Photon* photon = NULL;
486	<
487	<	VCOPY(fpos, ray -> rop);
488	<	VCOPY(norm, ray -> ron);
489	<	pmap -> maxDist = thescene.cusize;
490	<	nearest1Neighbour(pmap, fpos, norm, &photon, 1);
491	<
492	<	return photon;
493	<	}
494	<
495	<
496	<
497	<	static unsigned long medianPartition (const Photon* heap,
498	<	unsigned long* heapIdx,
499	<	unsigned long* heapXdi,
500	<	unsigned long left,
501	<	unsigned long right, unsigned dim)
502	<	/* Returns index to median in heap from indices left to right
503	<	(inclusive) in dimension dim. The heap is partitioned relative to
504	<	median using a quicksort algorithm. The heap indices in heapIdx are
505	<	sorted rather than the heap itself. */
506	<	{
507	<	register const float* p;
508	<	const unsigned long n = right - left + 1;
509	<	register unsigned long l, r, lg2, n2, m;
510	<	register unsigned d;
511	<
512	<	/* Round down n to nearest power of 2 */
513	<	for (lg2 = 0, n2 = n; n2 > 1; n2 >>= 1, ++lg2);
514	<	n2 = 1 << lg2;
515	<
516	<	/* Determine median position; this takes into account the fact that
517	<	only the last level in the heap can be partially empty, and that
518	<	it fills from left to right */
519	<	m = left + ((n - n2) > (n2 >> 1) - 1 ? n2 - 1 : n - (n2 >> 1));
520	<
521	<	while (right > left) {
522	<	/* Pivot node */
523	<	p = heap [heapIdx [right]].pos;
524	<	l = left;
525	<	r = right - 1;
571	>	#ifdef PMAP_OOC
572	>	if (OOC_GetPhoton(pmap, idx, photon))
573
574	<	/* l & r converge, swapping elements out of order with respect to
575	<	pivot node. Identical keys are resolved by cycling through
576	<	dim. The convergence point is then the pivot's position. */
577	<	do {
531	<	while (l <= r) {
532	<	d = dim;
533	<
534	<	while (heap [heapIdx [l]].pos [d] == p [d]) {
535	<	d = (d + 1) % 3;
536	<
537	<	if (d == dim) {
538	<	/* Ignore dupes? */
539	<	error(WARNING, "duplicate keys in photon heap");
540	<	l++;
541	<	break;
542	<	}
543	<	}
544	<
545	<	if (heap [heapIdx [l]].pos [d] < p [d])
546	<	l++;
547	<	else break;
548	<	}
549	<
550	<	while (r > l) {
551	<	d = dim;
552	<
553	<	while (heap [heapIdx [r]].pos [d] == p [d]) {
554	<	d = (d + 1) % 3;
555	<
556	<	if (d == dim) {
557	<	/* Ignore dupes? */
558	<	error(WARNING, "duplicate keys in photon heap");
559	<	r--;
560	<	break;
561	<	}
562	<	}
563	<
564	<	if (heap [heapIdx [r]].pos [d] > p [d])
565	<	r--;
566	<	else break;
567	<	}
568	<
569	<	/* Swap indices (not the nodes they point to) */
570	<	n2 = heapIdx [l];
571	<	heapIdx [l] = heapIdx [r];
572	<	heapIdx [r] = n2;
573	<	/* Update reverse indices */
574	<	heapXdi [heapIdx [l]] = l;
575	<	heapXdi [n2] = r;
576	<	} while (l < r);
577	<
578	<	/* Swap indices of convergence and pivot nodes */
579	<	heapIdx [r] = heapIdx [l];
580	<	heapIdx [l] = heapIdx [right];
581	<	heapIdx [right] = n2;
582	<	/* Update reverse indices */
583	<	heapXdi [heapIdx [r]] = r;
584	<	heapXdi [heapIdx [l]] = l;
585	<	heapXdi [n2] = right;
586	<	if (l >= m) right = l - 1;
587	<	if (l <= m) left = l + 1;
588	<	}
589	<
590	<	/* Once left & right have converged at m, we have found the median */
591	<	return m;
574	>	#else
575	>	if (kdT_GetPhoton(pmap, idx, photon))
576	>	#endif
577	>	error(INTERNAL, "failed photon lookup");
578		}
579
580
581
582	<	void buildHeap (Photon* heap, unsigned long* heapIdx,
597	<	unsigned long* heapXdi, const float min [3],
598	<	const float max [3], unsigned long left,
599	<	unsigned long right, unsigned long root)
600	<	/* Recursive part of balancePhotons(..). Builds heap from subarray
601	<	defined by indices left and right. min and max are the minimum resp.
602	<	maximum photon positions in the array. root is the index of the
603	<	current subtree's root, which corresponds to the median's 1-based
604	<	index in the heap. heapIdx are the balanced heap indices. The heap
605	<	is accessed indirectly through these. heapXdi are the reverse indices
606	<	from the heap to heapIdx so that heapXdi [heapIdx [i]] = i. */
582	>	Photon *getNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx)
583		{
584	<	float maxLeft [3], minRight [3];
585	<	Photon rootNode;
586	<	unsigned d;
587	<
588	<	/* Choose median for dimension with largest spread and partition
613	<	accordingly */
614	<	const float d0 = max [0] - min [0],
615	<	d1 = max [1] - min [1],
616	<	d2 = max [2] - min [2];
617	<	const unsigned char dim = d0 > d1 ? d0 > d2 ? 0 : 2
618	<	: d1 > d2 ? 1 : 2;
619	<	const unsigned long median =
620	<	left == right ? left
621	<	: medianPartition(heap, heapIdx, heapXdi, left, right, dim);
622	<
623	<	/* Place median at root of current subtree. This consists of swapping
624	<	the median and the root nodes and updating the heap indices */
625	<	memcpy(&rootNode, heap + heapIdx [median], sizeof(Photon));
626	<	memcpy(heap + heapIdx [median], heap + root - 1, sizeof(Photon));
627	<	setPhotonDiscr(rootNode, dim);
628	<	memcpy(heap + root - 1, &rootNode, sizeof(Photon));
629	<	heapIdx [heapXdi [root - 1]] = heapIdx [median];
630	<	heapXdi [heapIdx [median]] = heapXdi [root - 1];
631	<	heapIdx [median] = root - 1;
632	<	heapXdi [root - 1] = median;
633	<
634	<	/* Update bounds for left and right subtrees and recurse on them */
635	<	for (d = 0; d <= 2; d++)
636	<	if (d == dim)
637	<	maxLeft [d] = minRight [d] = rootNode.pos [d];
638	<	else {
639	<	maxLeft [d] = max [d];
640	<	minRight [d] = min [d];
641	<	}
642	<
643	<	if (left < median)
644	<	buildHeap(heap, heapIdx, heapXdi, min, maxLeft,
645	<	left, median - 1, root << 1);
646	<
647	<	if (right > median)
648	<	buildHeap(heap, heapIdx, heapXdi, minRight, max,
649	<	median + 1, right, (root << 1) + 1);
584	>	#ifdef PMAP_OOC
585	>	return OOC_GetNearestPhoton(squeue, idx);
586	>	#else
587	>	return kdT_GetNearestPhoton(squeue, idx);
588	>	#endif
589		}
590
591
592
593	<	void balancePhotons (PhotonMap* pmap, double *photonFlux)
593	>	PhotonIdx firstPhoton (const PhotonMap *pmap)
594		{
595	<	Photon *heap = pmap -> heap;
596	<	unsigned long i;
597	<	unsigned long *heapIdx;        /* Photon index array */
598	<	unsigned long *heapXdi;        /* Reverse index to heapIdx */
599	<	unsigned j;
661	<	COLOR flux;
662	<	/* Need doubles here to reduce errors from increment */
663	<	double avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0;
664	<	FVECT d;
665	<
666	<	if (pmap -> heapEnd) {
667	<	pmap -> heapSize = pmap -> heapEnd;
668	<	heapIdx = (unsigned long*)malloc(pmap -> heapSize *
669	<	sizeof(unsigned long));
670	<	heapXdi = (unsigned long*)malloc(pmap -> heapSize *
671	<	sizeof(unsigned long));
672	<	if (!heapIdx || !heapXdi)
673	<	error(USER, "can't allocate heap index");
674	<
675	<	for (i = 0; i < pmap -> heapSize; i++) {
676	<	/* Initialize index arrays */
677	<	heapXdi [i] = heapIdx [i] = i;
678	<	getPhotonFlux(heap + i, flux);
679	<
680	<	/* Scale photon's flux (hitherto normalised to 1 over RGB); in case
681	<	* of a contrib photon map, this is done per light source, and
682	<	* photonFlux is assumed to be an array */
683	<	if (photonFlux) {
684	<	scalecolor(flux, photonFlux [isContribPmap(pmap) ?
685	<	photonSrcIdx(pmap, heap + i) : 0]);
686	<	setPhotonFlux(heap + i, flux);
687	<	}
688	<
689	<	/* Need a double here */
690	<	addcolor(avgFlux, flux);
691	<
692	<	/* Add photon position to centre of gravity */
693	<	for (j = 0; j < 3; j++)
694	<	CoG [j] += heap [i].pos [j];
695	<	}
696	<
697	<	/* Average photon positions to get centre of gravity */
698	<	for (j = 0; j < 3; j++)
699	<	pmap -> CoG [j] = CoG [j] /= pmap -> heapSize;
700	<
701	<	/* Compute average photon distance to CoG */
702	<	for (i = 0; i < pmap -> heapSize; i++) {
703	<	VSUB(d, heap [i].pos, CoG);
704	<	CoGdist += DOT(d, d);
705	<	}
706	<
707	<	pmap -> CoGdist = CoGdist /= pmap -> heapSize;
708	<
709	<	/* Average photon flux based on RGBE representation */
710	<	scalecolor(avgFlux, 1.0 / pmap -> heapSize);
711	<	copycolor(pmap -> photonFlux, avgFlux);
712	<
713	<	/* Build kd-tree */
714	<	buildHeap(pmap -> heap, heapIdx, heapXdi, pmap -> minPos,
715	<	pmap -> maxPos, 0, pmap -> heapSize - 1, 1);
716	<
717	<	free(heapIdx);
718	<	free(heapXdi);
719	<	}
595	>	#ifdef PMAP_OOC
596	>	return OOC_FirstPhoton(pmap);
597	>	#else
598	>	return kdT_FirstPhoton(pmap);
599	>	#endif
600		}
601
602
603
604		void deletePhotons (PhotonMap* pmap)
605		{
606	<	free(pmap -> heap);
607	<	free(pmap -> squeue);
606	>	#ifdef PMAP_OOC
607	>	OOC_Delete(&pmap -> store);
608	>	#else
609	>	kdT_Delete(&pmap -> store);
610	>	#endif
611	>
612	>	free(pmap -> squeue.node);
613		free(pmap -> biasCompHist);
614
615	<	pmap -> heapSize = 0;
616	<	pmap -> minGather = pmap -> maxGather =
732	<	pmap -> squeueSize = pmap -> squeueEnd = 0;
615	>	pmap -> numPhotons = pmap -> minGather = pmap -> maxGather =
616	>	pmap -> squeue.len = pmap -> squeue.tail = 0;
617		}

Diff Legend

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