Homec4science
Diffusion RADIANCE Photon Map (pmapmerge) Browse

RadiancePmap/pmapsrc.c
49fd88e2dbf9pmapmerge

Last Change
Raw File

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);
}
c4science · Help