Diffusion RADIANCE Photon Map (pmapmerge) Browse

RadiancePmap/pmapsrc.c
8c7f4c6bcfb5pmapmerge

pmapsrc.c
View Options

			#ifndef lint
			static const char RCSid[] = "$Id: pmapsrc.c,v 2.20 2021/04/09 17:42:34 rschregle Exp $";
			#endif
			/*
			======================================================================
			Photon map support routines for emission from light sources

			Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
			(c) Fraunhofer Institute for Solar Energy Systems,
			supported by the German Research Foundation
			(DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme" (FARESYS))
			(c) Lucerne University of Applied Sciences and Arts,
			supported by the Swiss National Science Foundation
			(SNSF #147053, "Daylight Redirecting Components")
			(c) Tokyo University of Science,
			supported by the JSPS Grants-in-Aid for Scientific Research
			(KAKENHI JP19KK0115, "Three-Dimensional Light Flow")
			======================================================================

			$Id: pmapsrc.c,v 2.20 2021/04/09 17:42:34 rschregle Exp $"
			*/


			/* TODO: Add option to suppress backscattered photons on non-emitting side
			of ports? */


			#include "pmapsrc.h"
			#include "pmap.h"
			#include "pmaprand.h"
			#include "face.h"
			#include "otypes.h"
			#include "otspecial.h"
			#ifdef PMAP_PHOTONFLOW
			#include "pmapdens.h"
			#endif



			/* List of photon port modifier names */
			char *photonPortList [MAXSET + 1] = {NULL};
			/* Photon port objects (with modifiers in photonPortMods) */
			SRCREC *photonPorts = NULL;
			unsigned numPhotonPorts = 0;

			int (photonPartition [NUMOTYPE]) (EmissionMap);
			int (photonOrigin [NUMOTYPE]) (EmissionMap);



			/* PHOTON PORT SUPPORT ROUTINES ------------------------------------------ */



			/* Get/set photon port orientation flags from/in source flags.
			* HACK: the port orientation flags are embedded in the source flags and
			* shifted so they won't clobber the latter, since these are interpreted
			* by the PhotonPartition() and PhotonOrigin() routines! */
			#define PMAP_SETPORTFLAGS(portdir) ((int)(portdir) << 14)
			#define PMAP_GETPORTFLAGS(sflags) ((sflags) >> 14)

			/* Set number of source partitions.
			* HACK: this is doubled if the source acts as a bidirectionally
			* emitting photon port, resulting in alternating front/backside partitions,
			* although essentially each partition is just used twice with opposing
			* normals. */
			#define PMAP_SETNUMPARTITIONS(emap) ( \
			(emap) -> numPartitions <<= ( \
			(emap) -> port && \
			PMAP_GETPORTFLAGS((emap) -> port -> sflags) == PMAP_PORTBI \
			) \
			)

			/* Get current source partition and numer of partitions
			* HACK: halve the number partitions if the source acts as a bidrectionally
			* emitting photon port, since each partition is used twice with opposing
			* normals. */
			#define PMAP_GETNUMPARTITIONS(emap) (\
			(emap) -> numPartitions >> ( \
			(emap) -> port && \
			PMAP_GETPORTFLAGS((emap) -> port -> sflags) == PMAP_PORTBI \
			) \
			)
			#define PMAP_GETPARTITION(emap) ( \
			(emap) -> partitionCnt >> ( \
			(emap) -> port && \
			PMAP_GETPORTFLAGS((emap) -> port -> sflags) == PMAP_PORTBI \
			) \
			)



			void getPhotonPorts (char **portList)
			/* Find geometry declared as photon ports from modifiers in portList */
			{
			OBJECT i;
			OBJREC obj, mat;
			unsigned lastNumPhotonPorts;
			int mLen;
			char **lp;

			/* Init photon port objects */
			photonPorts = NULL;
			numPhotonPorts = 0;

			if (!portList [0])
			return;

			for (lp = portList; *lp; lp++) {
			lastNumPhotonPorts = numPhotonPorts;

			for (i = 0; i < nobjects; i++) {
			obj = objptr(i);
			mat = findmaterial(obj);

			/* Check if object is a surface and NOT a light source (duh) and
			* resolve its material (if any) via any aliases, then check for
			* inclusion in modifier list; note use of strncmp() to ignore port
			* flags */
			if (mat && !islight(mat -> otype) && issurface(obj -> otype)) {
			mLen = strlen(mat -> oname);

			if (!strncmp(mat -> oname, *lp, mLen)) {
			/* Add photon port */
			photonPorts = (SRCREC*)realloc(
			photonPorts, (numPhotonPorts + 1) * sizeof(SRCREC)
			);
			if (!photonPorts)
			error(USER, "can't allocate photon ports");

			photonPorts [numPhotonPorts].so = obj;
			/* Extract port orientation flags and embed in source flags.
			* Note setsource() combines (i.e. ORs) these with the actual
			* source flags below. */
			photonPorts [numPhotonPorts].sflags =
			PMAP_SETPORTFLAGS((*lp) [mLen]);

			if (!sfun [obj -> otype].of \|\| !sfun[obj -> otype].of -> setsrc)
			objerror(obj, USER, "illegal photon port");

			setsource(photonPorts + numPhotonPorts, obj);
			numPhotonPorts++;
			}
			}
			}

			if (numPhotonPorts == lastNumPhotonPorts) {
			/* No port added for this modifier */
			sprintf(errmsg, "no void ports found for modifier %s", *lp);
			error(USER, errmsg);
			}
			}

			if (!numPhotonPorts)
			/* Now redundant due to check above */
			error(USER, "no valid photon ports found");
			}



			static void setPhotonPortNormal (EmissionMap* emap)
			/* Set normal for current photon port partition (if defined) based on its
			* orientation */
			{

			int i, portFlags;

			if (emap -> port) {
			/* Extract photon port orientation flags, set surface normal as follows:
			-- Port oriented forwards --> flip surface normal to point outwards,
			since normal points inwards per mkillum convention)
			-- Port oriented backwards --> surface normal is NOT flipped, since
			it already points inwards.
			-- Port is bidirectionally/bilaterally oriented --> flip normal based
			on the parity of the current partition emap -> partitionCnt. In
			this case, photon emission alternates between port front/back
			faces for consecutive partitions.
			*/
			portFlags = PMAP_GETPORTFLAGS(emap -> port -> sflags);

			if (
			portFlags == PMAP_PORTFWD \|\|
			portFlags == PMAP_PORTBI && !(emap -> partitionCnt & 1)
			)
			for (i = 0; i < 3; i++)
			emap -> ws [i] = -emap -> ws [i];
			}
			}



			/* SOURCE / PHOTON PORT PARTITIONING ROUTINES----------------------------- */



			static int flatPhotonPartition2 (
			EmissionMap* emap, unsigned long mp, FVECT cent, FVECT u, FVECT v,
			double du2, double dv2
			)
			/* Recursive part of flatPhotonPartition(..) */
			{
			FVECT newct, newax;
			unsigned long npl, npu;

			if (mp > emap -> maxPartitions) {
			/* Enlarge partition array */
			emap -> maxPartitions <<= 1;
			emap -> partitions = (unsigned char*)realloc(
			emap -> partitions, emap -> maxPartitions >> 1
			);

			if (!emap -> partitions)
			error(USER, "can't allocate source partitions");

			memset(
			emap -> partitions + (emap -> maxPartitions >> 2), 0,
			emap -> maxPartitions >> 2
			);
			}

			if (du2 * dv2 <= 1) { /* hit limit? */
			setpart(emap -> partitions, emap -> partitionCnt, S0);
			emap -> partitionCnt++;
			return 1;
			}

			if (du2 > dv2) { /* subdivide in U */
			setpart(emap -> partitions, emap -> partitionCnt, SU);
			emap -> partitionCnt++;
			newax [0] = 0.5 * u [0];
			newax [1] = 0.5 * u [1];
			newax [2] = 0.5 * u [2];
			u = newax;
			du2 *= 0.25;
			}

			else { /* subdivide in V */
			setpart(emap -> partitions, emap -> partitionCnt, SV);
			emap -> partitionCnt++;
			newax [0] = 0.5 * v [0];
			newax [1] = 0.5 * v [1];
			newax [2] = 0.5 * v [2];
			v = newax;
			dv2 *= 0.25;
			}

			/* lower half */
			newct [0] = cent [0] - newax [0];
			newct [1] = cent [1] - newax [1];
			newct [2] = cent [2] - newax [2];
			npl = flatPhotonPartition2(emap, mp << 1, newct, u, v, du2, dv2);

			/* upper half */
			newct [0] = cent [0] + newax [0];
			newct [1] = cent [1] + newax [1];
			newct [2] = cent [2] + newax [2];
			npu = flatPhotonPartition2(emap, mp << 1, newct, u, v, du2, dv2);

			/* return total */
			return npl + npu;
			}



			static int flatPhotonPartition (EmissionMap* emap)
			/* Partition flat source for photon emission */
			{
			RREAL *vp;
			double du2, dv2;

			memset(emap -> partitions, 0, emap -> maxPartitions >> 1);
			emap -> partArea = srcsizerat * thescene.cusize;
			emap -> partArea *= emap -> partArea;
			vp = emap -> src -> ss [SU];
			du2 = DOT(vp, vp) / emap -> partArea;
			vp = emap -> src -> ss [SV];
			dv2 = DOT(vp, vp) / emap -> partArea;
			emap -> partitionCnt = 0;
			emap -> numPartitions = flatPhotonPartition2(
			emap, 1, emap -> src -> sloc,
			emap -> src -> ss [SU], emap -> src -> ss [SV], du2, dv2
			);
			emap -> partitionCnt = 0;
			emap -> partArea = emap -> src -> ss2 / emap -> numPartitions;

			return 1;
			}



			static int sourcePhotonPartition (EmissionMap* emap)
			/* Partition scene cube faces or photon port for photon emission from
			distant source */
			{
			if (emap -> port) {
			/* Relay partitioning to photon port */
			SRCREC *src = emap -> src;
			emap -> src = emap -> port;
			photonPartition [emap -> src -> so -> otype] (emap);
			PMAP_SETNUMPARTITIONS(emap);
			emap -> src = src;
			}

			else {
			/* No photon ports defined; partition scene cube faces (SUBOPTIMAL) */
			memset(emap -> partitions, 0, emap -> maxPartitions >> 1);
			setpart(emap -> partitions, 0, S0);
			emap -> partitionCnt = 0;
			emap -> numPartitions = 1 / srcsizerat;
			emap -> numPartitions *= emap -> numPartitions;
			emap -> partArea = sqr(thescene.cusize) / emap -> numPartitions;
			emap -> numPartitions *= 6;
			}

			return 1;
			}



			static int spherePhotonPartition (EmissionMap* emap)
			/* Partition spherical source into equal solid angles using uniform
			mapping for photon emission */
			{
			unsigned numTheta, numPhi;

			memset(emap -> partitions, 0, emap -> maxPartitions >> 1);
			setpart(emap -> partitions, 0, S0);
			emap -> partArea = 4 * PI * sqr(emap -> src -> srad);
			emap -> numPartitions =
			emap -> partArea / sqr(srcsizerat * thescene.cusize);

			numTheta = max(sqrt(2 * emap -> numPartitions / PI) + 0.5, 1);
			numPhi = 0.5 * PI * numTheta + 0.5;

			emap -> numPartitions = (unsigned long)numTheta * numPhi;
			emap -> partitionCnt = 0;
			emap -> partArea /= emap -> numPartitions;

			return 1;
			}



			static int cylPhotonPartition2 (
			EmissionMap* emap, unsigned long mp, FVECT cent, FVECT axis, double d2
			)
			/* Recursive part of cyPhotonPartition(..) */
			{
			FVECT newct, newax;
			unsigned long npl, npu;

			if (mp > emap -> maxPartitions) {
			/* Enlarge partition array */
			emap -> maxPartitions <<= 1;
			emap -> partitions = (unsigned char*)realloc(
			emap -> partitions, emap -> maxPartitions >> 1
			);
			if (!emap -> partitions)
			error(USER, "can't allocate source partitions");

			memset(
			emap -> partitions + (emap -> maxPartitions >> 2), 0,
			emap -> maxPartitions >> 2
			);
			}

			if (d2 <= 1) {
			/* hit limit? */
			setpart(emap -> partitions, emap -> partitionCnt, S0);
			emap -> partitionCnt++;
			return 1;
			}

			/* subdivide */
			setpart(emap -> partitions, emap -> partitionCnt, SU);
			emap -> partitionCnt++;
			newax [0] = 0.5 * axis [0];
			newax [1] = 0.5 * axis [1];
			newax [2] = 0.5 * axis [2];
			d2 *= 0.25;

			/* lower half */
			newct [0] = cent [0] - newax [0];
			newct [1] = cent [1] - newax [1];
			newct [2] = cent [2] - newax [2];
			npl = cylPhotonPartition2(emap, mp << 1, newct, newax, d2);

			/* upper half */
			newct [0] = cent [0] + newax [0];
			newct [1] = cent [1] + newax [1];
			newct [2] = cent [2] + newax [2];
			npu = cylPhotonPartition2(emap, mp << 1, newct, newax, d2);

			/* return total */
			return npl + npu;
			}



			static int cylPhotonPartition (EmissionMap* emap)
			/* Partition cylindrical source for photon emission */
			{
			double d2;

			memset(emap -> partitions, 0, emap -> maxPartitions >> 1);
			d2 = srcsizerat * thescene.cusize;
			d2 = PI * emap -> src -> ss2 / (2 * emap -> src -> srad * sqr(d2));
			d2 = d2 DOT(emap -> src -> ss [SU], emap -> src -> ss [SU]);

			emap -> partitionCnt = 0;
			emap -> numPartitions = cylPhotonPartition2(
			emap, 1, emap -> src -> sloc, emap -> src -> ss [SU], d2
			);
			emap -> partitionCnt = 0;
			emap -> partArea = PI * emap -> src -> ss2 / emap -> numPartitions;

			return 1;
			}



			/* PHOTON ORIGIN ROUTINES ------------------------------------------------ */



			static int flatPhotonOrigin (EmissionMap* emap)
			/* Init emission map with photon origin and associated surface axes on flat
			light source surface. Also sets source aperture and sampling hemisphere
			axes at origin. Returns 1 if photon lies within source boundary, else 0
			(signalling rejection to trim non-quadrilateral geometry).
			Some code fragments are adapted from srcsamp.c:nextssamp(). */
			{
			int i, cent[3], size[3], parr[2];
			FVECT vpos;
			const SRCREC *src = emap -> src;

			/* Get surface axes */
			VCOPY(emap -> us, src -> ss [SU]);
			normalize(emap -> us);
			VCOPY(emap -> ws, src -> ss [SW]);
			/* Flip normal emap -> ws if port and required by its orientation */
			setPhotonPortNormal(emap);
			fcross(emap -> vs, emap -> ws, emap -> us);

			/* Get hemisphere axes & aperture */
			if (src -> sflags & SSPOT) {
			VCOPY(emap -> wh, src -> sl.s -> aim);
			i = 0;
			do {
			emap -> vh [0] = emap -> vh [1] = emap -> vh [2] = 0;
			emap -> vh [i++] = 1;
			fcross(emap -> uh, emap -> vh, emap -> wh);
			} while (normalize(emap -> uh) < FTINY);

			fcross(emap -> vh, emap -> wh, emap -> uh);
			emap -> cosThetaMax = 1 - src -> sl.s -> siz / (2 * PI);
			}
			else {
			VCOPY(emap -> uh, emap -> us);
			VCOPY(emap -> vh, emap -> vs);
			VCOPY(emap -> wh, emap -> ws);
			emap -> cosThetaMax = 0;
			}

			cent [SU] = cent [SV] = cent [SW] = 0;
			size [SU] = size [SV] = size [SW] = emap -> maxPartitions;
			parr [0] = 0;
			parr [1] = PMAP_GETPARTITION(emap);

			if (!skipparts(cent, size, parr, emap -> partitions))
			error(CONSISTENCY, "bad source partition in flatPhotonOrigin");

			vpos [SU] = (1 - 2 * pmapRandom(partState)) * size [SU] /
			emap -> maxPartitions;
			vpos [SV] = (1 - 2 * pmapRandom(partState)) * size [SV] /
			emap -> maxPartitions;
			vpos [SW] = 0;

			for (i = 0; i < 3; i++)
			vpos [i] += (double)cent [i] / emap -> maxPartitions;

			if (src -> sflags & SCIR) {
			/* Disk source; trim origin by correcting partition extent
			(note implicit trim [SW] = 0) */
			FVECT trim;
			trim [SU] = 1.12837917 * sqrt(1 - 0.5 * vpos [SV] * vpos [SV]);
			trim [SV] = 1.12837917 * sqrt(1 - 0.5 * vpos [SU] * vpos [SU]);

			for (i = 0; i < 3; i++)
			emap -> photonOrg [i] = src -> sloc [i] +
			trim [SU] * vpos [SU] * src -> ss [SU][i] +
			trim [SV] * vpos [SV] * src -> ss [SV][i];

			/* Accept origin unconditionally */
			return 1;
			}
			else {
			/* Polygonal source */
			FACE *inYaFace = getface(src -> so);

			if (inYaFace -> nv == 3) {
			/* Triangle; sample separately; code adapted from rfluxmtx.c */
			vpos [SU] = 0.5 * (vpos [SU] + 1);
			vpos [SV] = 0.5 * (vpos [SV] + 1);
			vpos [SU] *= vpos [SV] = sqrt(vpos [SV]);
			vpos [SV] = 1 - vpos [SV];

			for (i = 0; i < 3; i++) {
			const RREAL v0 = VERTEX(inYaFace, 0) [i];
			emap -> photonOrg [i] = v0 +
			vpos [SU] * (VERTEX(inYaFace, 1) [i] - v0) +
			vpos [SV] * (VERTEX(inYaFace, 2) [i] - v0);
			}

			/* Accept photon origin unconditionally */
			return 1;
			}
			else if (inYaFace -> nv == 4) {
			/* Quadrilateral; use standard sampling from srcsamp.c */
			for (i = 0; i < 3; i++)
			emap -> photonOrg [i] =
			emap -> src -> sloc [i] +
			vpos [SU] * emap -> src -> ss [SU][i] +
			vpos [SV] * emap -> src -> ss [SV][i] +
			vpos [SW] * emap -> src -> ss [SW][i];

			/* Accept photon origin unconditionally */
			return 1;
			}
			else {
			/* Arbitrary polygon */
			for (i = 0; i < 3; i++)
			emap -> photonOrg [i] = emap -> src -> sloc [i] +
			vpos [SU] * src -> srad * emap -> us [i] +
			vpos [SV] * src -> srad * emap -> vs [i];
			/* Accept origin only if within polygon bounds
			(= rejection sampling) */
			return inface(emap -> photonOrg, inYaFace);
			}
			}
			}



			static int spherePhotonOrigin (EmissionMap* emap)
			/* Init emission map with photon origin and associated surface axes on
			spherical light source. Also sets source aperture and sampling
			hemisphere axes at origin */
			{
			int i = 0;
			unsigned numTheta, numPhi, t, p;
			RREAL cosTheta, sinTheta, phi;

			/* Get current partition */
			numTheta = max(sqrt(2 * PMAP_GETNUMPARTITIONS(emap) / PI) + 0.5, 1);
			numPhi = 0.5 * PI * numTheta + 0.5;

			t = PMAP_GETPARTITION(emap) / numPhi;
			p = PMAP_GETPARTITION(emap) - t * numPhi;

			emap -> ws [2] = cosTheta = 1 - 2 * (t + pmapRandom(partState)) / numTheta;
			sinTheta = sqrt(1 - sqr(cosTheta));
			phi = 2 * PI * (p + pmapRandom(partState)) / numPhi;
			emap -> ws [0] = cos(phi) * sinTheta;
			emap -> ws [1] = sin(phi) * sinTheta;
			/* Flip normal emap -> ws if port and required by its orientation */
			setPhotonPortNormal(emap);

			/* Get surface axes us & vs perpendicular to ws */
			do {
			emap -> vs [0] = emap -> vs [1] = emap -> vs [2] = 0;
			emap -> vs [i++] = 1;
			fcross(emap -> us, emap -> vs, emap -> ws);
			} while (normalize(emap -> us) < FTINY);

			fcross(emap -> vs, emap -> ws, emap -> us);

			/* Get origin */
			for (i = 0; i < 3; i++)
			emap -> photonOrg [i] = emap -> src -> sloc [i] +
			emap -> src -> srad * emap -> ws [i];

			/* Get hemisphere axes & aperture */
			if (emap -> src -> sflags & SSPOT) {
			VCOPY(emap -> wh, emap -> src -> sl.s -> aim);
			i = 0;

			do {
			emap -> vh [0] = emap -> vh [1] = emap -> vh [2] = 0;
			emap -> vh [i++] = 1;
			fcross(emap -> uh, emap -> vh, emap -> wh);
			} while (normalize(emap -> uh) < FTINY);

			fcross(emap -> vh, emap -> wh, emap -> uh);
			emap -> cosThetaMax = 1 - emap -> src -> sl.s -> siz / (2 * PI);
			}

			else {
			VCOPY(emap -> uh, emap -> us);
			VCOPY(emap -> vh, emap -> vs);
			VCOPY(emap -> wh, emap -> ws);
			emap -> cosThetaMax = 0;
			}

			return 1;
			}



			static int sourcePhotonOrigin (EmissionMap* emap)
			/* Init emission map with photon origin and associated surface axes
			on scene cube face for distant light source. Also sets source
			aperture (solid angle) and sampling hemisphere axes at origin */
			{
			unsigned long i, partsPerDim, partsPerFace;
			unsigned face;
			RREAL du, dv;
			int acceptOrigin = 1;

			if (emap -> port) {
			/* Relay to photon port; get origin on its surface */
			SRCREC *src = emap -> src;
			emap -> src = emap -> port;
			acceptOrigin = photonOrigin [emap -> src -> so -> otype] (emap);
			emap -> src = src;
			}

			else {
			/* No ports defined, so get origin on scene cube face (SUBOPTIMAL) */
			/* Get current face from partition number */
			partsPerDim = 1 / srcsizerat;
			partsPerFace = sqr(partsPerDim);
			face = emap -> partitionCnt / partsPerFace;
			if (!(emap -> partitionCnt % partsPerFace)) {
			/* Skipped to a new face; get its normal */
			emap -> ws [0] = emap -> ws [1] = emap -> ws [2] = 0;
			emap -> ws [face >> 1] = face & 1 ? 1 : -1;

			/* Get surface axes us & vs perpendicular to ws */
			face >>= 1;
			emap -> vs [0] = emap -> vs [1] = emap -> vs [2] = 0;
			emap -> vs [(face + (emap -> ws [face] > 0 ? 2 : 1)) % 3] = 1;
			fcross(emap -> us, emap -> vs, emap -> ws);
			}

			/* Get jittered offsets within face from partition number
			(in range [-0.5, 0.5]) */
			i = emap -> partitionCnt % partsPerFace;
			du = (i / partsPerDim + pmapRandom(partState)) / partsPerDim - 0.5;
			dv = (i % partsPerDim + pmapRandom(partState)) / partsPerDim - 0.5;

			/* Jittered destination point within partition */
			for (i = 0; i < 3; i++)
			emap -> photonOrg [i] = thescene.cuorg [i] + thescene.cusize * (
			0.5 + du * emap -> us [i] + dv * emap -> vs [i] +
			0.5 * emap -> ws [i]
			);
			}

			/* Get hemisphere axes & aperture */
			VCOPY(emap -> wh, emap -> src -> sloc);
			i = 0;

			do {
			emap -> vh [0] = emap -> vh [1] = emap -> vh [2] = 0;
			emap -> vh [i++] = 1;
			fcross(emap -> uh, emap -> vh, emap -> wh);
			} while (normalize(emap -> uh) < FTINY);

			fcross(emap -> vh, emap -> wh, emap -> uh);

			/* Get aperture */
			emap -> cosThetaMax = 1 - emap -> src -> ss2 / (2 * PI);
			emap -> cosThetaMax = min(1, max(-1, emap -> cosThetaMax));

			return acceptOrigin;
			}



			static int cylPhotonOrigin (EmissionMap* emap)
			/* Init emission map with photon origin and associated surface axes
			on cylindrical light source surface. Also sets source aperture
			and sampling hemisphere axes at origin */
			{
			int i, cent[3], size[3], parr[2];
			FVECT v;

			cent [0] = cent [1] = cent [2] = 0;
			size [0] = size [1] = size [2] = emap -> maxPartitions;
			parr [0] = 0;
			parr [1] = PMAP_GETPARTITION(emap);

			if (!skipparts(cent, size, parr, emap -> partitions))
			error(CONSISTENCY, "bad source partition in cylPhotonOrigin");

			v [SU] = 0;
			v [SV] = (1 - 2 * pmapRandom(partState)) * (double)size [1] /
			emap -> maxPartitions;
			v [SW] = (1 - 2 * pmapRandom(partState)) * (double)size [2] /
			emap -> maxPartitions;
			normalize(v);
			v [SU] = (1 - 2 * pmapRandom(partState)) * (double)size [1] /
			emap -> maxPartitions;

			for (i = 0; i < 3; i++)
			v [i] += (double)cent [i] / emap -> maxPartitions;

			/* Get surface axes */
			for (i = 0; i < 3; i++)
			emap -> photonOrg [i] = emap -> ws [i] = (
			v [SV] * emap -> src -> ss [SV][i] +
			v [SW] * emap -> src -> ss [SW][i]
			) / 0.8559;

			/* Flip normal emap -> ws if port and required by its orientation */
			setPhotonPortNormal(emap);

			normalize(emap -> ws);
			VCOPY(emap -> us, emap -> src -> ss [SU]);
			normalize(emap -> us);
			fcross(emap -> vs, emap -> ws, emap -> us);

			/* Get origin */
			for (i = 0; i < 3; i++)
			emap -> photonOrg [i] +=
			v [SU] * emap -> src -> ss [SU][i] + emap -> src -> sloc [i];

			/* Get hemisphere axes & aperture */
			if (emap -> src -> sflags & SSPOT) {
			VCOPY(emap -> wh, emap -> src -> sl.s -> aim);
			i = 0;

			do {
			emap -> vh [0] = emap -> vh [1] = emap -> vh [2] = 0;
			emap -> vh [i++] = 1;
			fcross(emap -> uh, emap -> vh, emap -> wh);
			} while (normalize(emap -> uh) < FTINY);

			fcross(emap -> vh, emap -> wh, emap -> uh);
			emap -> cosThetaMax = 1 - emap -> src -> sl.s -> siz / (2 * PI);
			}

			else {
			VCOPY(emap -> uh, emap -> us);
			VCOPY(emap -> vh, emap -> vs);
			VCOPY(emap -> wh, emap -> ws);
			emap -> cosThetaMax = 0;
			}

			return 1;
			}



			/* PHOTON EMISSION ROUTINES ---------------------------------------------- */



			static int defaultEmissionFunc (EmissionMap* emap)
			/* Default behaviour when no emission funcs defined for this source type */
			{
			objerror(emap -> src -> so, INTERNAL,
			"undefined photon emission function"
			);

			return 0;
			}



			void initPhotonEmissionFuncs ()
			/* Init photonPartition[] and photonOrigin[] dispatch tables */
			{
			int i;

			for (i = 0; i < NUMOTYPE; i++)
			photonPartition [i] = photonOrigin [i] = defaultEmissionFunc;

			photonPartition [OBJ_FACE] = photonPartition [OBJ_RING] =
			flatPhotonPartition;
			photonPartition [OBJ_SOURCE] = sourcePhotonPartition;
			photonPartition [OBJ_SPHERE] = spherePhotonPartition;
			photonPartition [OBJ_CYLINDER] = cylPhotonPartition;
			photonOrigin [OBJ_FACE] = photonOrigin [OBJ_RING] = flatPhotonOrigin;
			photonOrigin [OBJ_SOURCE] = sourcePhotonOrigin;
			photonOrigin [OBJ_SPHERE] = spherePhotonOrigin;
			photonOrigin [OBJ_CYLINDER] = cylPhotonOrigin;
			}



			void initPhotonEmission (EmissionMap *emap, float numPdfSamples)
			/* Initialise photon emission from partitioned light source emap -> src;
			* this involves integrating the flux emitted from the current photon
			* origin emap -> photonOrg and setting up a PDF to sample the emission
			* distribution with numPdfSamples samples */
			{
			unsigned i, t, p;
			double phi, cosTheta, sinTheta, du, dv, dOmega;
			EmissionSample* sample;
			const OBJREC* mod = findmaterial(emap -> src -> so);
			static RAY r;

			/* Skip virtual sources (we doan' need them with pmap) */
			if (emap -> src -> sflags & SVIRTUAL)
			error(INTERNAL, "got unexpected virtual source");

			photonOrigin [emap -> src -> so -> otype] (emap);
			setcolor(emap -> partFlux, 0, 0, 0);
			emap -> cdf = 0;
			emap -> numSamples = 0;
			cosTheta = DOT(emap -> ws, emap -> wh);

			if (cosTheta <= 0 && sqrt(1-sqr(cosTheta)) <= emap -> cosThetaMax + FTINY)
			/* Aperture completely below surface; no emission from current origin */
			return;

			if (mod -> omod == OVOID && !emap -> port && (cosTheta >= 1 - FTINY \|\| (
			emap -> src -> sflags & SDISTANT &&
			acos(cosTheta) + acos(emap -> cosThetaMax) <= 0.5 * PI
			))) {
			/* Source is unmodified and has no port (which requires testing for
			occlusion), and is either local with normal aligned aperture or
			distant with aperture above surface
			--> get flux & PDF via analytical solution */
			setcolor(emap -> partFlux, mod -> oargs.farg [0],
			mod -> oargs.farg [1], mod -> oargs.farg [2]
			);

			/* Multiply radiance by projected Omega * dA to get flux */
			scalecolor(emap -> partFlux,
			PI * cosTheta * (1 - sqr(max(emap -> cosThetaMax, 0))) *
			emap -> partArea
			);
			}

			else {
			/* Source is either modified, has a port, is local with off-normal
			aperture, or distant with aperture partly below surface
			--> get flux via numerical integration using aperture samples
			with constant differential solid angle */

			/* Figga out numba of aperture samples for integration;
			numTheta / numPhi ratio is 1 / PI */
			emap -> dCosTheta = (1 - emap -> cosThetaMax);
			du = sqrt(pdfSamples * 2 * emap -> dCosTheta);
			emap -> numTheta = max(du + 0.5, 1);
			emap -> numPhi = max(PI * du + 0.5, 1);
			emap -> dCosTheta /= emap -> numTheta;
			emap -> dPhi = 2 * PI / emap -> numPhi;
			/* Differential solid angle is constant for this mapping */
			dOmega = emap -> dCosTheta * emap -> dPhi;

			/* Allocate PDF, baby */
			sample = emap -> samples = (EmissionSample*)realloc(emap -> samples,
			sizeof(EmissionSample) * emap -> numTheta * emap -> numPhi
			);
			if (!emap -> samples)
			error(USER, "can't allocate emission PDF");

			VCOPY(r.rorg, emap -> photonOrg);
			VCOPY(r.rop, emap -> photonOrg);
			r.rmax = 0;

			for (t = 0; t < emap -> numTheta; t++) {
			for (p = 0; p < emap -> numPhi; p++) {
			/* This uniform mapping handles 0 <= thetaMax <= PI */
			cosTheta = 1 - (t + pmapRandom(emitState)) * emap -> dCosTheta;
			sinTheta = sqrt(1 - sqr(cosTheta));
			phi = (p + pmapRandom(emitState)) * emap -> dPhi;
			du = cos(phi) * sinTheta;
			dv = sin(phi) * sinTheta;
			rayorigin(&r, PRIMARY, NULL, NULL);

			for (i = 0; i < 3; i++)
			r.rdir [i] = (du * emap -> uh [i] + dv * emap -> vh [i] +
			cosTheta * emap -> wh [i]
			);

			/* Sample behind surface? */
			VCOPY(r.ron, emap -> ws);
			if ((r.rod = DOT(r.rdir, r.ron)) <= 0)
			continue;

			/* Get radiance emitted in this direction; to get flux we
			multiply by cos(theta_surface), dOmega, and dA. Ray
			is directed towards light source for raytexture(). */
			if (!(emap -> src -> sflags & SDISTANT))
			for (i = 0; i < 3; i++)
			r.rdir [i] = -r.rdir [i];

			/* Port occluded in this direction? */
			if (emap -> port && localhit(&r, &thescene))
			continue;

			raytexture(&r, mod -> omod);
			setcolor(r.rcol, mod -> oargs.farg [0], mod -> oargs.farg [1],
			mod -> oargs.farg [2]
			);
			multcolor(r.rcol, r.pcol);

			/* Multiply by projection cos(theta_surface) */
			scalecolor(r.rcol, r.rod);

			/* Add PDF sample if nonzero.
			TODO: importance info for photon emission could go here.*/
			if (colorAvg(r.rcol)) {
			copycolor(sample -> pdf, r.rcol);
			sample -> cdf = emap -> cdf += colorAvg(sample -> pdf);
			sample -> theta = t;
			sample++ -> phi = p;
			emap -> numSamples++;
			addcolor(emap -> partFlux, r.rcol);
			}
			}
			}

			/* Multiply by differential angle & area, dOmega * dA */
			scalecolor(emap -> partFlux, dOmega * emap -> partArea);
			}
			}



			#define vomitPhoton emitPhoton
			#define bluarrrghPhoton vomitPhoton

			int emitPhoton (const EmissionMap* emap, RAY* ray)
			/* Emit photon from current partition emap -> partitionCnt based on
			emission distribution. Returns new photon ray and an acceptance/rejection
			flag: 1 if origin lies within the emitting source/port, else 0 if it lies
			outside, in which case the emitted photon should be rejected. */
			{
			unsigned long i, lo, hi;
			const EmissionSample* sample = emap -> samples;
			RREAL du, dv, cosTheta, cosThetaSqr, sinTheta, phi;
			const OBJREC* mod = findmaterial(emap -> src -> so);

			/* Choose a new origin within current partition for every
			emitted photon to break up clustering artifacts */
			if (!photonOrigin [emap -> src -> so -> otype] ((EmissionMap*)emap))
			/* Origin lies outside source/port geometry; signal its rejection */
			return 0;

			/* If we have a local glow source with a maximum radius, then
			restrict our photon to the specified distance, otherwise we set
			the limit imposed by photonMaxPathDist (or no limit if 0) */
			if (mod -> otype == MAT_GLOW && !(emap -> src -> sflags & SDISTANT) &&
			mod -> oargs.farg[3] > FTINY
			)
			ray -> rmax = mod -> oargs.farg[3];
			else
			ray -> rmax = photonMaxPathDist;
			rayorigin(ray, PRIMARY, NULL, NULL);

			/* Init directional samples */
			if (!emap -> numSamples) {
			/* Source is unmodified and has no port, and either local with
			normal aligned aperture, or distant with aperture above surface
			--> use cosine weighted distribution */
			cosThetaSqr = (1 -
			pmapRandom(emitState) * (1 - sqr(max(emap -> cosThetaMax, 0)))
			);
			cosTheta = sqrt(cosThetaSqr);
			sinTheta = sqrt(1 - cosThetaSqr);
			phi = 2 * PI * pmapRandom(emitState);
			setcolor(
			ray -> rcol, mod -> oargs.farg [0], mod -> oargs.farg [1],
			mod -> oargs.farg [2]
			);
			}

			else {
			/* Source is either modified, has a port, is local with off-normal
			aperture, or distant with aperture partly below surface
			--> choose direction from constructed cumulative distribution
			function with Monte Carlo inversion using binary search. */
			du = pmapRandom(emitState) * emap -> cdf;
			lo = 1;
			hi = emap -> numSamples;

			while (hi > lo) {
			i = (lo + hi) >> 1;
			sample = emap -> samples + i - 1;

			if (sample -> cdf >= du)
			hi = i;
			if (sample -> cdf < du)
			lo = i + 1;
			}

			/* Finalise found sample */
			i = (lo + hi) >> 1;
			sample = emap -> samples + i - 1;

			/* This is a uniform mapping, mon */
			cosTheta = (
			1 - (sample -> theta + pmapRandom(emitState)) * emap -> dCosTheta
			);
			sinTheta = sqrt(1 - sqr(cosTheta));
			phi = (sample -> phi + pmapRandom(emitState)) * emap -> dPhi;
			copycolor(ray -> rcol, sample -> pdf);
			}

			/* Normalize photon flux so that average over RGB is 1 */
			colorNorm(ray -> rcol);

			VCOPY(ray -> rorg, emap -> photonOrg);
			du = cos(phi) * sinTheta;
			dv = sin(phi) * sinTheta;

			for (i = 0; i < 3; i++)
			ray -> rdir [i] = (
			du * emap -> uh [i] + dv * emap -> vh [i] + cosTheta * emap -> wh [i]
			);

			if (emap -> src -> sflags & SDISTANT)
			/* Distant source; reverse ray direction to point into the scene. */
			for (i = 0; i < 3; i++)
			ray -> rdir [i] = -ray -> rdir [i];

			if (emap -> port) {
			/* Photon emitted from port; move origin just behind port so it
			will be scattered */
			for (i = 0; i < 3; i++)
			ray -> rorg [i] -= 2 * FTINY * ray -> rdir [i];
			/* !!! PHOTON PORT REJECTION SAMPLING HACK: set photon port as fake
			!!! hit object for primary ray to check for its intersection in
			!!! tracePhoton() */
			ray -> ro = emap -> port -> so;
			}

			/* Assign emitting light source index */
			ray -> rsrc = emap -> src - source;

			/* Signal photon's acceptance */
			return 1;
			}



			/* SOURCE CONTRIBS FROM DIRECT / VOLUME PHOTONS -------------------------- */



			void multDirectPmap (RAY *r)
			/* Factor irradiance from direct photons into r -> rcol; interface to
			* direct() */
			{
			COLOR photonIrrad;

			#ifdef PMAP_PHOTONFLOW
			if (lightFlowPhotonMapping)
			/* Photon flow mode: evaluate volume photon map as lightfield
			* (including direct component), ignoring participating medium. */
			(lightFlowPmap -> lookup)(lightFlowPmap, r, photonIrrad);
			else
			#endif
			if (transientPhotonMapping)
			/* Lookup transient photons (direct & indirect) */
			(transientPmap -> lookup)(transientPmap, r, photonIrrad);
			else
			/* Lookup direct photon irradiance */
			(directPmap -> lookup)(directPmap, r, photonIrrad);

			/* Multiply with coefficient in ray */
			multcolor(r -> rcol, photonIrrad);

			return;
			}



			void inscatterVolumePmap (RAY *r, COLOR inscatter)
			/* Add inscattering from volume photon map; interface to srcscatter() */
			{
			/* Add ambient in-scattering via lookup callback */
			(volumePmap -> lookup)(volumePmap, r, inscatter);
			}

RadiancePmap/pmapsrc.c8c7f4c6bcfb5pmapmerge

pmapsrc.cView Options

RadiancePmap/pmapsrc.c
8c7f4c6bcfb5pmapmerge

pmapsrc.c
View Options