diff --git a/pmapsrc.c b/pmapsrc.c index b42f51b..e51cf5d 100644 --- a/pmapsrc.c +++ b/pmapsrc.c @@ -1,977 +1,983 @@ #ifndef lint static const char RCSid[] = "$Id: pmapsrc.c,v 2.19 2020/08/10 19:51:20 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.19 2020/08/10 19:51:20 rschregle Exp $" */ #include "pmapsrc.h" #include "pmap.h" #include "pmaprand.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; void (*photonPartition [NUMOTYPE]) (EmissionMap*); void (*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; int mLen; char **lp; /* Init photon port objects */ photonPorts = NULL; if (!portList [0]) return; for (i = numPhotonPorts = 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 (issurface(obj -> otype) && mat && !islight(mat -> otype)) { mLen = strlen(mat -> oname); for (lp = portList; *lp && strncmp(mat -> oname, *lp, mLen); lp++); if (*lp) { /* 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) 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 void 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; } static void 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; } } static void 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; } 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 void 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; } /* PHOTON ORIGIN ROUTINES ------------------------------------------------ */ static void 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 */ { int i, cent[3], size[3], parr[2]; FVECT vpos; 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 flatPhotonOrigin"); vpos [0] = (1 - 2 * pmapRandom(partState)) * size [0] / emap -> maxPartitions; vpos [1] = (1 - 2 * pmapRandom(partState)) * size [1] / emap -> maxPartitions; vpos [2] = 0; for (i = 0; i < 3; i++) vpos [i] += (double)cent [i] / emap -> maxPartitions; /* Get origin */ 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]; /* Get surface axes */ VCOPY(emap -> us, emap -> src -> ss [SU]); normalize(emap -> us); VCOPY(emap -> ws, emap -> 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 (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; } } static void 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; } } static void 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; if (emap -> port) { /* Relay to photon port; get origin on its surface */ SRCREC *src = emap -> src; emap -> src = emap -> port; 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)); } static void 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; } } /* PHOTON EMISSION ROUTINES ---------------------------------------------- */ static void defaultEmissionFunc (EmissionMap* emap) /* Default behaviour when no emission funcs defined for this source type */ { objerror(emap -> src -> so, INTERNAL, "undefined photon emission function"); } 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, thetaScale; + double phi, cosTheta, sinTheta, du, dv, dOmega; EmissionSample* sample; const OBJREC* mod = findmaterial(emap -> src -> so); static RAY r; photonOrigin [emap -> src -> so -> otype] (emap); setcolor(emap -> partFlux, 0, 0, 0); emap -> cdf = 0; emap -> numSamples = 0; - cosTheta = DOT(emap -> ws, emap -> wh); + 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 */ - thetaScale = (1 - emap -> cosThetaMax); - + --> 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 */ - du = sqrt(pdfSamples * 2 * thetaScale); + 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); - - dOmega = 2 * PI * thetaScale / (emap -> numTheta * emap -> numPhi); - thetaScale /= emap -> numTheta; + 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)) * thetaScale; + cosTheta = 1 - (t + pmapRandom(emitState)) * emap -> dCosTheta; sinTheta = sqrt(1 - sqr(cosTheta)); - phi = 2 * PI * (p + pmapRandom(emitState)) / emap -> numPhi; + 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 cos(theta_surface) */ + + /* Multiply by projection cos(theta_surface) */ scalecolor(r.rcol, r.rod); - - /* Add PDF sample if nonzero; importance info for photon emission - * could go here... */ + + /* 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 dOmega * dA */ + /* Multiply by differential angle & area, dOmega * dA */ scalecolor(emap -> partFlux, dOmega * emap -> partArea); } } #define vomitPhoton emitPhoton #define bluarrrghPhoton vomitPhoton void emitPhoton (const EmissionMap* emap, RAY* ray) /* Emit photon from current partition emap -> partitionCnt based on emission distribution. Returns new photon ray. */ { 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 */ photonOrigin [emap -> src -> so -> otype] ((EmissionMap*)emap); /* 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 photonMaxDist (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 = photonMaxDist; rayorigin(ray, PRIMARY, NULL, NULL); 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)) * - (1 - emap -> cosThetaMax) / emap -> numTheta + cosTheta = ( + 1 - (sample -> theta + pmapRandom(emitState)) * emap -> dCosTheta ); sinTheta = sqrt(1 - sqr(cosTheta)); - phi = 2 * PI * (sample -> phi + pmapRandom(emitState)) / emap -> numPhi; + 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); + 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]; - + 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. */ + /* 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]; - + ray -> rorg [i] -= 2 * FTINY * ray -> rdir [i]; + /* Assign emitting light source index */ ray -> rsrc = emap -> src - source; } /* 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 (photonFlow) if (volumePhotonMapping) /* Photon flow mode: evaluate volume photon map as lightfield * (including direct component), ignoring participating medium. */ volumePhotonSphIrrad(volumePmap, r, photonIrrad); else error(USER, "photon flow mode requires a volume photon map"); else #endif /* 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); } diff --git a/pmapsrc.h b/pmapsrc.h index f26ab7f..54f6113 100644 --- a/pmapsrc.h +++ b/pmapsrc.h @@ -1,107 +1,108 @@ /* RCSid $Id: pmapsrc.h,v 2.7 2020/08/07 01:21:13 rschregle Exp $ */ /* ====================================================================== 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.h,v 2.7 2020/08/07 01:21:13 rschregle Exp $ */ #ifndef PMAPSRC_H #define PMAPSRC_H #include "ray.h" #include "source.h" /* Data structures for photon emission */ typedef struct { unsigned theta, phi; /* Direction indices */ COLOR pdf; /* Probability of emission in this direction per RGB */ float cdf; /* Cumulative probability up to this sample */ } EmissionSample; typedef struct { unsigned numTheta, numPhi; /* Num divisions in theta & phi */ + double dCosTheta, dPhi; /* Increments for above */ RREAL cosThetaMax; /* cos(source aperture) */ FVECT uh, vh, wh; /* Emission aperture axes at origin, w is normal*/ FVECT us, vs, ws; /* Source surface axes at origin, w is normal */ FVECT photonOrg; /* Current photon origin */ EmissionSample *samples; /* PDF samples for photon emission directions */ unsigned long numPartitions; /* Number of surface partitions */ RREAL partArea; /* Area covered by each partition */ SRCREC *src, *port; /* Current source and port */ unsigned long partitionCnt, /* Current surface partition */ maxPartitions, /* Max allocated partitions */ numSamples; /* Number of PDF samples */ unsigned char* partitions; /* Source surface partitions */ COLOR partFlux; /* Total flux emitted by partition */ float cdf; /* Cumulative probability */ } EmissionMap; /* Photon port flags (orientation relative to surface normal): * Forward, backward, both (bidirectional). */ #define PMAP_PORTFWD 1 #define PMAP_PORTBWD 2 #define PMAP_PORTBI (PMAP_PORTFWD | PMAP_PORTBWD) /* Photon ports for emission from geometry en lieu of light sources */ extern char *photonPortList [MAXSET + 1]; extern SRCREC *photonPorts; extern unsigned numPhotonPorts; /* Dispatch tables for partitioning light source surfaces and generating an origin for photon emission */ extern void (*photonPartition []) (EmissionMap*); extern void (*photonOrigin []) (EmissionMap*); void getPhotonPorts (char **portList); /* Find geometry declared as photon ports from modifiers in portList */ void initPhotonEmissionFuncs (); /* Init photonPartition[] and photonOrigin[] dispatch tables with source specific emission routines */ void initPhotonEmission (EmissionMap *emap, float numPdfSamples); /* Initialize 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 */ void emitPhoton (const EmissionMap*, RAY* ray); /* Emit photon from current source partition based on emission distribution and return new photon ray */ void multDirectPmap (RAY *r); /* Factor irradiance from direct photons into r -> rcol; interface to * direct() */ void inscatterVolumePmap (RAY *r, COLOR inscatter); /* Add inscattering from volume photon map; interface to srcscatter() */ #endif