diff --git a/README-photonflow.txt b/README-photonflow.txt index 39ef716..ab343e3 100644 --- a/README-photonflow.txt +++ b/README-photonflow.txt @@ -1,261 +1,275 @@ This is the README file for the RADIANCE photon flow implementation. This software was developed as part of the project "Three-Dimensional Light Flow", hosted by Tokyo University of Science under the supervision of Prof. Nozomu Yoshizawa, and supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI) under the grant number JP19KK0115. $Revision: 1.3 $ $Date: 2020/11/16 22:41:09 $ INTRODUCTION These source files implement a modified volume photon density estimate in RADIANCE to evaluate the illuminance of a physical lightfield using photon flow. This code is an experimental variant of the RADIANCE photon map extension intended for testing prior to committing it to the RADIANCE CVS as part of the official distribution. Whether this (admittedly specialised) functionality should be made available to the RADIANCE community as a default build option is debatable and may be decided in due course. This code is essentially a C implementation of the prototype Python script volPhotonIllum.py. As such, it is orders of magnitude faster and more convenient, which comes to bear with large photon maps and many sensor positions. One crucial difference between the RADIANCE and Python implementations is that the former performs photon lookups in terms of a fixed number of photons around the sensor position (nearest neighbour search), while the latter does so in terms of a fixed radius (range search), and therefore does not adapt to the local photon density, which is suboptimal. The illuminance evaluated from photon flow approximates Cuttle's cubic illuminance for a virtual sensor position in a physical light field. Photon flow constitutes volume photons deposited by mkpmap in a participating medium (mist); due to the interaction of a photon with the medium, it will be deposited at multiple locations along its path; the density of these photons per path is governed by the participating medium's extinction. Evaluating cubic illuminance from photon flow entails performing volume photon lookups for each of the 6 (virtual) cube faces corresponding to 3 orthogonal axes in positive and negative directions. The resulting search volume is a hemisphere, and photons incident from the respective cube face (identified by it is normal or axis) are filtered on the fly based on their stored incident directions. This results in six independent hemispherical lookups with individual radii adapted to the density of the incident photons for each axis. Although this incurs more overhead, it was found that this approach yields consistently better results than a single lookup for all axes with a posteriori filtering, particularly with asymmetric illuminance. It was also found that the illuminance scales with the extinction of the participating medium, and is therefore corrected for the same. Generally, smaller values yield more accurate results but incur longer photon distribution times with mkpmap. The photon mapping code was modified to generate a variant of the existing volume photon map, referred to as a "lightflow photon map", which corrects the photon flux according to the extinction during it generation with mkpmap. When loaded with rtrace, a lightflow photon map triggers the cubic illuminance evaluation via the function volumePhotonCubicIllum(). Unlike the existing evaluation function for volume photons, volumePhotonDensity(), the lightflow function does not apply the Henyey-Greenstein phase function to compute outscattering (thus ignoring the participating medium), but instead performs the per-axis lookups with a volume density estimate as described above. Because photon flow not only includes photons that have been reflected off one or more surfaces (indirect contributions), but also photons that have just been emitted from light sources (direct contributions), volumePhotonCubicIllum() is is called by multiDirectPmap() in pmapsrc.c, which normally calls the direct photon map density estimate (typically only used for debugging). To prevent overcounting, the ambient calculation is therefore disabled in photon flow mode (see ampPmap() in pmapamb.c). COMPILING The photon flow code is not compiled by default as part of the RADIANCE suite, and must be enabled at compile time by defining PMAP_PHOTONFLOW. This may be accelerated by only (re)compiling in ray/src/rt if a RADIANCE build already exists: rmake OPT+=-DPMAP_PHOTONFLOW install Note that these options are independent of the underlying data structure used to store the photon map. The default options build the legacy kd-tree. If a larger, out-of-core photon map is desired (to accommodate >1 billion photons, for example), this build option must be enabled by defining PMAP_OOC: rmake OPT+="-DPMAP_OOC -DPMAP_PHOTONFLOW" install For details, see the RADIANCE Out-Of-Core Photon Map Technical Report at http://dx.doi.org/10.13140/RG.2.1.2158.9363 Similar to stanard (planar) photon mapping, photon flow supports kernel density estimates by weighting individual photon contributions by their inverse distance to the lookup point. This filter reduces bias, and can be enabled with the following: rmake OPT+=-DPMAP_EPANECHNIKOV -DPMAP_OOC -DPMAP_PHOTONFLOW" install This command enables the 3D Epanechnikov kernel (1 - (d_i)^2), where d_i is the normalised distance between the lookup point and photon (i). The bias reduction is moderate, however, and ineffective with excessive lookup bandwidths. Finally, the following checks verify that the photon flow code was successfully built: 1. 'mkpmap -help' lists the '-apV' option to generate lightflow photon maps, 2. 'rtrace -help' lists the '-aS' option to toggle between planar and spherical illuminance evaluation from lightflow photons. USAGE A lightflow photon map representing the physical lightfield is generated with mkpmap by specifying the extinction of a global participating medium either on the command line, or in a RADIANCE scene description using a mist material with boundary geometry (e.g. to limit photonflow to a specific region of interest). The participating medium must not absorb nor scatter, so that a photon's trajectory is not altered as it traverses the medium. Therefore as a convenience, mkpmap implicitly sets the albedo and scattering eccentricity to 1.0, which is equivalent to specifying the options -ma 1 1 1 (no absorption) and -mg 1 (forward scattering): mkpmap -apV lightflow.vpm 1M [-apo port1 ... -apo portN] -apD 0.001 \ -me .1 .1 .1 scene.oct The extinction -me is a free parameter here and determines the photon density per path. Tests have shown that lower values improve accuracy but deposit fewer photons per path, requiring mkpmap to emit more photon paths in order to satisfy the given target photon count. There is a caveat with using higher values of -me, namely that the two-pass photon distribution algorithm used by mkpmap to achieve the target number of photons may fail. The fraction of photons emitted in the 1st pass (prepass) is specified with -apD (predistribution factor). Its default value of 0.25 is tuned for "well behaved" geometry and materials using global or caustic photons. This value may be too high for volume photons in a dense participating medium, and the photon distribution overflows in the prepass with a "photon map overflow" error, meaning that the target photon count was already exceeded (technically this does not yield an incorrect solution, but it's not what the user wanted). When using photon flow, this value may need to be significantly reduced from its default of the prepass overflows. Once the lightflow photon map has been generated, the RGB irradiance can be evaluated using rtrace. This evaluation is enabled when loading the lightflow photon map, which disables the normal ambient calculation. Because the ambient calculation is entirely circumvented, the -ab parameter is now irrelevant; the spherical volume photon density is evaluated immediately at each sensor position. Each normal vector supplied with each sensor point corresponds to the axis for each of the 6 cube faces, while the position is (or should be!) identical for all 6 input records. Though redundant, this parametrisation is consistent with rtrace convention. sensors.dat: <-u> ... <-w> rtrace -h -I -ap lightflow.vpm 1000 scene.oct < sensors.dat Output: ... In this example, the RGB irradiance E_i from photon flow is evaluated with a bandwidth of 1000 photons for each axis i denoting one of the orthogonal vectors {u, v, w} or its reverse. The cubic scalar illuminance Esr can be calculated from these values in a separate script as a postprocess. The advantage of this discrete approach is that it affords access to the RGB components of each axis irradiance E_i, which would be lost if the scalar illuminance were calculated in rtrace directly. In addition, the mean spherical irradiance can be calculated with rtrace to compare with the cubic scalar illuminance Esr. Rtrace accepts a boolean option -aS to toggle between the planar irradiance evaluation for E_i as above (corresponding to the default -aS-), and the mean spherical irradiance evaluation using -aS+ (or just -aS): rtrace -h -I -ap lightflow.vpm 1000 -aS+ scene.oct < sph_sensor.dat Note that the normal vector specified for the sensor position is irrelevant and ignored for this evaluation. +As a further option, rtrace accepts an -aH option to toggle between spherical +and hemispherical lookups for lightflow photons.. In the latter case (-aH+), +photons located behind the plane are not considered, which can reduce bias +near "solid" boundaries (scene cube, geometry), behind which no photons are +deposited. However, this increases the search radius compared to a spherical +lookup, and can lead to proximity bias if the illumination gradient is high. +For lookup points primarily located in free space, it is recommended to use +spherical lookups (-aH-), which is the default. For example, the command + + rtrace -h -I -ap lightflow.vpm 1000 -aH+ scene.oct < wall_sensor.dat + +will locate lightflow photons on the near side of a solid wall. Note that +the -aH+ and -aS+ options are mutually exclusive! + Finally, the lightflow photon distribution may be visualised by either dumping it as a RADIANCE scene description, and rendering the photons as spheres, e.g. pmapdump -n 1m lightflow.vpm | objview or by dumping it as a point list for use in a point cloud viewer: pmapdump -a -f -n 1m lightflow.vpm > photonFlow.xyz In the latter example, pmapdump also outputs each photon's radiant flux (in W) for further processing. See the pmapdump man page or the RADIANCE photon map manual at http://dx.doi.org/10.13140/RG.2.1.4330.8405/1 for details. TRANSIENT PHOTON FLOW An experimental variant of the light flow photon map that supports transient (i.e. time-dependent) rendering is currently being tested. This variant tags each photon with its corresponding path length, which corresponds to a time stamp. In order for the time stamp to be physically valid, the user must specify the speed of light in units of scene geometry (i.e. a multiple of approx. 3e8 m/s) when generating the photon map with mkpmap. This constant is specified in addition to the photon map size with the type-specific -apT option: mkpmap ... -apT lightflow.tpm 1M 3e8 ... scene.oct Transient lightflow can then be evaluated as a spatio-temporal function by specifying a point in time in addition to a position: rtrace -h -I -ap lightflow.tpm 250 1.6e-9 scene.oct < sensor.dat In the above example, the irradiance/illuminance will be evaluated at the given sensor position AND the given time: 1.6 nsec after photon emission begins. This spatio-temporal distribution can also be visualised by passing the parameters into rpict, which leads to surprising effects. To this end, the Python script transpmap.py is being developed to render spatio-temporal light flow sequences and merge these into a video. For example, the commands transpmap.py -t 0 -T 8e-9 -d 1e-10 -p 10m -b 250 -n 8 cornell.vf cornell.oct will render transient lightflow from 10 million photons with a bandwidth of 250 photons on 8 cores within the time interval [0, 8 nsec] in increments of 100 psec per frame. See 'transpmap.py -h' for the supported options. CONCLUSION Having reassessed the approach to evaluate cubic illuminance from volume photons, revised the code accordingly, and tested it with asymmetric illuminance distributions, I feel this release is easier to use, better integrated into the existing RADIANCE framework, and more stable and mature. However, the code, like the approach it implements, is still experimental and needs some more testing in "real world" scenarious, so I welcome your feedback in terms of its results, but also its usage. Once we are confident with the results, we should make the code public and upload it to the RADIANCE CVS repository. Initially, though, it would be prudent to still leave photon flow an option for advanced users only. Many thanks in advance for your cooperation and feedback. Happy photon flowing! --Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) diff --git a/pmap.h b/pmap.h index 5b8ebd3..1e85bb3 100644 --- a/pmap.h +++ b/pmap.h @@ -1,115 +1,112 @@ /* RCSid $Id: pmap.h,v 2.9 2018/01/24 19:39:05 rschregle Exp $ */ /* ====================================================================== Photon map main header 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: pmap.h,v 2.9 2018/01/24 19:39:05 rschregle Exp $ */ /* TODO: Combine type-specific properties of pmap types as binary flags, e.g. PMAP_TYPE_TRANSLIGHTFLOW = PMAP_TYPE_VOLUME | PMAP_TYPE_TRANSIENT | PMAP_TYPE_LIGHTFLOW ??? */ #ifndef PMAP_H #define PMAP_H #ifndef NIX #if defined(_WIN32) || defined(_WIN64) #define NIX 0 #else #define NIX 1 #endif #endif #include "pmapparm.h" #include "pmapdata.h" #ifndef min #define min(a, b) ((a) < (b) ? (a) : (b)) #endif #ifndef max #define max(a, b) ((a) > (b) ? (a) : (b)) #endif #define sqr(a) ((a) * (a)) /* Average over colour channels */ #define colorAvg(col) ((col [0] + col [1] + col [2]) / 3) /* Macros to test for enabled photon map (also see pmapdata.h) */ #define photonMapping (\ globalPmap || preCompPmap || causticPmap || contribPmap || \ transientPmap \ ) #define causticPhotonMapping (causticPmap) #define directPhotonMapping (directPmap || transientPhotonMapping) #define volumePhotonMapping (volumePmap) /* #define contribPhotonMapping (contribPmap && contribPmap -> srcContrib) */ #define contribPhotonMapping (contribPmap) #define transientPhotonMapping (transientPmap) #ifdef PMAP_PHOTONFLOW - /* Enable hemispherical lookups in photon flow mode; this reduces - bias near solid surfaces */ - #define PMAP_PHOTONFLOW_HEMI #define lightFlowPhotonMapping (lightFlowPmap) #endif extern void (*pmapLookup [])(PhotonMap*, RAY*, COLOR); /* Photon map lookup functions per type */ void loadPmaps (PhotonMap **pmaps, const PhotonMapParams *params); /* Load photon and set their respective parameters, checking timestamps * relative to octree for possible staleness */ void savePmaps (const PhotonMap **pmaps, int argc, char **argv); /* Save all defined photon maps with specified command line */ void cleanUpPmaps (PhotonMap **pmaps); /* Trash all photon maps after processing is complete */ void distribPhotons (PhotonMap **pmaps, unsigned numProc); /* Emit photons from light sources and build photon maps for non-NULL * entries in photon map array */ void tracePhoton (RAY*); /* Follow photon as it bounces around the scene. Analogon to * raytrace(). */ void photonDensity (PhotonMap*, RAY*, COLOR irrad); /* Perform surface density estimate from incoming photon flux at ray's intersection point. Returns irradiance from found photons. */ void photonPreCompDensity (PhotonMap *pmap, RAY *r, COLOR irrad); /* Returns precomputed photon density estimate at ray -> rop. */ void volumePhotonDensity (PhotonMap*, RAY*, COLOR); /* Perform volume density estimate from incoming photon flux at ray's intersection point. Returns irradiance. */ void colorNorm (COLOR); /* Normalise colour channels to average of 1 */ #endif diff --git a/pmapdata.c b/pmapdata.c index 2540df8..8c3e64f 100644 --- a/pmapdata.c +++ b/pmapdata.c @@ -1,840 +1,836 @@ #ifndef lint static const char RCSid[] = "$Id: pmapdata.c,v 2.23 2021/04/14 11:26:25 rschregle Exp $"; #endif /* ========================================================================= Photon map types and high level interface to nearest neighbour lookups in underlying point cloud data structure. The default data structure is an in-core kd-tree (see pmapkdt.{h,c}). This can be overriden with the PMAP_OOC compiletime switch, which enables an out-of-core octree (see oococt.{h,c}). 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", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ========================================================================= $Id: pmapdata.c,v 2.23 2021/04/14 11:26:25 rschregle Exp $ */ #include "pmapdata.h" #include "pmapdens.h" #include "pmaprand.h" #include "pmapmat.h" #include "otypes.h" #include "random.h" PhotonMap *photonMaps [NUM_PMAP_TYPES] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL #ifdef PMAP_PHOTONFLOW , NULL, NULL #endif }; /* Include routines to handle underlying point cloud data structure */ #ifdef PMAP_OOC #include "pmapooc.c" #else #include "pmapkdt.c" #include "pmaptkdt.c" #endif /* Ambient include/exclude set (from ambient.c) */ #ifndef MAXASET #define MAXASET 4095 #endif extern OBJECT ambset [MAXASET+1]; /* Callback to print photon attributes acc. to user defined format */ int (*printPhoton)(RAY *r, Photon *p, PhotonMap *pm); /* PHOTON MAP BUILD ROUTINES ------------------------------------------ */ void initPhotonMap (PhotonMap *pmap, PhotonMapType t) /* Init photon map 'n' stuff... */ { if (!pmap) return; pmap -> numPhotons = 0; pmap -> biasCompHist = NULL; pmap -> maxPos [0] = pmap -> maxPos [1] = pmap -> maxPos [2] = -FHUGE; pmap -> minPos [0] = pmap -> minPos [1] = pmap -> minPos [2] = FHUGE; pmap -> minGathered = pmap -> maxGathered = pmap -> totalGathered = 0; pmap -> gatherTolerance = gatherTolerance; pmap -> minError = pmap -> maxError = pmap -> rmsError = 0; pmap -> numDensity = 0; pmap -> distribRatio = 1; pmap -> type = t; pmap -> squeue.node = NULL; pmap -> squeue.len = 0; #ifdef PMAP_PATHFILT /* Init path lookup table and its key buffer */ pmap -> pathLUT = NULL; pmap -> pathLUTKeys = NULL; pmap -> numPathLUTKeys = 0; #endif /* Init transient photon mapping stuff */ pmap -> velocity = photonVelocity; pmap -> minPathLen = pmap -> maxPathLen = pmap -> avgPathLen = 0; /* Init local RNG state */ pmap -> randState [0] = 10243; pmap -> randState [1] = 39829; pmap -> randState [2] = 9433; pmapSeed(randSeed, pmap -> randState); /* Set up type-specific photon lookup callback */ pmap -> lookup = pmapLookup [t]; /* Mark photon contribution source as unused */ pmap -> lastContribSrc.srcIdx = pmap -> lastContribSrc.srcBin = -1; pmap -> numContribSrc = 0; pmap -> contribSrc = NULL; /* Init storage */ pmap -> heap = NULL; pmap -> heapBuf = NULL; pmap -> heapBufLen = 0; #ifdef PMAP_OOC OOC_Null(&pmap -> store); #else kdT_Null(&pmap -> store); #endif } void initPhotonHeap (PhotonMap *pmap) { int fdFlags; if (!pmap) error(INTERNAL, "undefined photon map in initPhotonHeap()"); if (!pmap -> heap) { /* Open heap file */ mktemp(strcpy(pmap -> heapFname, PMAP_TMPFNAME)); if (!(pmap -> heap = fopen(pmap -> heapFname, "w+b"))) error(SYSTEM, "failed opening heap file in initPhotonHeap()"); #ifdef F_SETFL /* XXX is there an alternative needed for Windows? */ fdFlags = fcntl(fileno(pmap -> heap), F_GETFL); fcntl(fileno(pmap -> heap), F_SETFL, fdFlags | O_APPEND); #endif/* ftruncate(fileno(pmap -> heap), 0); */ } } void flushPhotonHeap (PhotonMap *pmap) { int fd; const unsigned long len = pmap -> heapBufLen * sizeof(Photon); if (!pmap) error(INTERNAL, "undefined photon map in flushPhotonHeap()"); if (!pmap -> heap || !pmap -> heapBuf) { /* Silently ignore undefined heap error(INTERNAL, "undefined heap in flushPhotonHeap()"); */ return; } /* Atomically seek and write block to heap */ /* !!! Unbuffered I/O via pwrite() avoids potential race conditions * !!! and buffer corruption which can occur with lseek()/fseek() * !!! followed by write()/fwrite(). */ fd = fileno(pmap -> heap); #ifdef DEBUG_PMAP sprintf(errmsg, "Proc %d: flushing %ld photons from pos %ld\n", getpid(), pmap -> heapBufLen, lseek(fd, 0, SEEK_END) / sizeof(Photon) ); eputs(errmsg); #endif /*if (pwrite(fd, pmap -> heapBuf, len, lseek(fd, 0, SEEK_END)) != len) */ if (write(fd, pmap -> heapBuf, len) != len) { error(SYSTEM, "failed append to heap file in flushPhotonHeap()"); /* Clean up temp file */ fclose(pmap -> heap); unlink(pmap -> heapFname); } #if NIX if (fsync(fd)) error(SYSTEM, "failed fsync in flushPhotonHeap()"); #endif pmap -> heapBufLen = 0; } #ifdef DEBUG_PMAP static int checkPhotonHeap (FILE *file) /* Check heap for nonsensical or duplicate photons */ { Photon p, lastp; int i, dup; rewind(file); memset(&lastp, 0, sizeof(lastp)); while (fread(&p, sizeof(p), 1, file)) { dup = 1; for (i = 0; i <= 2; i++) { if (p.pos [i] < thescene.cuorg [i] || p.pos [i] > thescene.cuorg [i] + thescene.cusize ) { sprintf(errmsg, "corrupt photon in heap at [%f, %f, %f]\n", p.pos [0], p.pos [1], p.pos [2] ); error(WARNING, errmsg); } dup &= p.pos [i] == lastp.pos [i]; } if (dup) { sprintf(errmsg, "consecutive duplicate photon in heap " "at [%f, %f, %f]\n", p.pos [0], p.pos [1], p.pos [2] ); error(WARNING, errmsg); } } return 0; } #endif int newPhoton (PhotonMap* pmap, const RAY* ray) { unsigned i; Photon photon; COLOR photonFlux; /* Account for distribution ratio */ if (!pmap || pmapRandom(pmap -> randState) > pmap -> distribRatio) return -1; /* Don't store on sources */ if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype)) return -1; /* Ignore photon if modifier in/outside exclude/include set */ if (ambincl != -1 && ray -> ro && ambincl != inset(ambset, ray->ro->omod) ) return -1; if (pmapNumROI && pmapROI) { unsigned inROI = 0; FVECT photonDist; /* Store photon if within a region of interest (for ze Ecksperts!) Note size of spherical ROI is squared. */ for (i = 0; !inROI && i < pmapNumROI; i++) { VSUB(photonDist, ray -> rop, pmapROI [i].pos); inROI = (PMAP_ROI_ISSPHERE(pmapROI + i) ? DOT(photonDist, photonDist) <= pmapROI [i].siz [0] : fabs(photonDist [0]) <= pmapROI [i].siz [0] && fabs(photonDist [1]) <= pmapROI [i].siz [1] && fabs(photonDist [2]) <= pmapROI [i].siz [2] ); } if (!inROI) return -1; } /* Adjust flux according to distribution ratio and ray weight */ copycolor(photonFlux, ray -> rcol); #if 0 /* Factored out ray -> rweight as deprecated (?) for pmap, and infact erroneously attenuates volume photon flux based on extinction, which is already factored in by photonParticipate() */ scalecolor( photonFlux, ray -> rweight / (pmap -> distribRatio ? pmap -> distribRatio : 1) ); #else scalecolor(photonFlux, 1.0 / (pmap -> distribRatio ? pmap -> distribRatio : 1) ); #endif setPhotonFlux(&photon, photonFlux); /* Set photon position and flags */ VCOPY(photon.pos, ray -> rop); photon.flags = 0; photon.caustic = PMAP_CAUSTICRAY(ray); /* Set photon's subprocess index */ photon.proc = PMAP_GETRAYPROC(ray); switch (pmap -> type) { case PMAP_TYPE_CONTRIB: /* Set contrib photon's origin and subprocess index (the latter to * linearise the origin array after photon distribution completes). * Also set contribution source index, thereby marking it as used. * Note the contribution source bin has already been set by * newPhotonContribSrc(). */ photon.aux.contribSrc = pmap -> numContribSrc; pmap -> lastContribSrc.srcIdx = ray -> rsrc; break; case PMAP_TYPE_TRANSIENT: #ifdef PMAP_PHOTONFLOW case PMAP_TYPE_TRANSLIGHTFLOW: #endif /* Set cumulative path length for transient photon, corresponding to its time of flight. The cumulative path length is obtained relative to the maximum path length photonMaxPathDist. The last intersection distance is subtracted as this isn't factored in until photonRay() generates a new photon ray. */ - #if 1 photon.aux.pathLen = photonMaxPathDist - ray -> rmax + ray -> rot; - #else - photon.aux.pathLen = photonMaxPathDist - ray -> rmax; - #endif break; default: /* Set photon's path ID from ray serial num (supplemented by photon.proc below) */ photon.aux.pathID = ray -> rno; } /* Set normal */ for (i = 0; i <= 2; i++) photon.norm [i] = 127.0 * (isVolumePmap(pmap) #ifdef PMAP_PHOTONFLOW || isLightFlowPmap(pmap) #endif ? ray -> rdir [i] : ray -> ron [i] ); if (!pmap -> heapBuf) { /* Lazily allocate heap buffa */ #if NIX /* Randomise buffa size to temporally decorellate flushes in * multiprocessing mode */ srandom(randSeed + getpid()); pmap -> heapBufSize = PMAP_HEAPBUFSIZE * (0.5 + frandom()); #else /* Randomisation disabled for single processes on WIN; also useful * for reproducability during debugging */ pmap -> heapBufSize = PMAP_HEAPBUFSIZE; #endif if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)))) error(SYSTEM, "failed heap buffer allocation in newPhoton"); pmap -> heapBufLen = 0; } /* Photon initialised; now append to heap buffa */ memcpy(pmap -> heapBuf + pmap -> heapBufLen, &photon, sizeof(Photon)); if (++pmap -> heapBufLen >= pmap -> heapBufSize) /* Heap buffa full, flush to heap file */ flushPhotonHeap(pmap); pmap -> numPhotons++; /* Print photon attributes */ if (printPhoton) /* Non-const kludge */ printPhoton((RAY*)ray, &photon, pmap); return 0; } void buildPhotonMap (PhotonMap *pmap, double *photonFlux, PhotonContribSourceIdx *contribSrcOfs, unsigned nproc ) { unsigned long n, nCheck = 0; unsigned i; Photon *p; COLOR flux; char nuHeapFname [sizeof(PMAP_TMPFNAME)]; FILE *nuHeap; /* Need double here to reduce summation errors */ double avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0, avgPathLen = 0; FVECT d; if (!pmap) error(INTERNAL, "undefined photon map in buildPhotonMap()"); /* Get number of photons from heapfile size */ if (fseek(pmap -> heap, 0, SEEK_END) < 0) error(SYSTEM, "failed seek to end of photon heap in buildPhotonMap()"); pmap -> numPhotons = ftell(pmap -> heap) / sizeof(Photon); if (!pmap -> numPhotons) error(INTERNAL, "empty photon map in buildPhotonMap()"); if (!pmap -> heap) error(INTERNAL, "no heap in buildPhotonMap()"); #ifdef DEBUG_PMAP eputs("Checking photon heap consistency...\n"); checkPhotonHeap(pmap -> heap); sprintf(errmsg, "Heap contains %ld photons\n", pmap -> numPhotons); eputs(errmsg); #endif /* Allocate heap buffa */ if (!pmap -> heapBuf) { pmap -> heapBufSize = PMAP_HEAPBUFSIZE; pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)); if (!pmap -> heapBuf) error(SYSTEM, "failed to allocate postprocessed photon heap in buildPhotonMap()" ); } /* We REALLY don't need yet another @%&*! heap just to hold the scaled * photons, but can't think of a quicker fix... */ mktemp(strcpy(nuHeapFname, PMAP_TMPFNAME)); if (!(nuHeap = fopen(nuHeapFname, "w+b"))) error(SYSTEM, "failed to open postprocessed photon heap in buildPhotonMap()" ); rewind(pmap -> heap); #ifdef DEBUG_PMAP eputs("Postprocessing photons...\n"); #endif while (!feof(pmap -> heap)) { #ifdef DEBUG_PMAP printf("Reading %lu at %lu... ", pmap -> heapBufSize, ftell(pmap->heap) ); #endif pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufSize, pmap -> heap ); #ifdef DEBUG_PMAP printf("Got %lu\n", pmap -> heapBufLen); #endif if (ferror(pmap -> heap)) error(SYSTEM, "failed to read photon heap in buildPhotonMap()"); for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) { /* Update min and max pos and set photon flux */ for (i = 0; i <= 2; i++) { if (p -> pos [i] < pmap -> minPos [i]) pmap -> minPos [i] = p -> pos [i]; else if (p -> pos [i] > pmap -> maxPos [i]) pmap -> maxPos [i] = p -> pos [i]; /* Update centre of gravity with photon position */ CoG [i] += p -> pos [i]; } if (isContribPmap(pmap) && contribSrcOfs) /* Linearise contrib photon origin index from subprocess index * using the per-subprocess offsets in contribSrcOfs */ p -> aux.contribSrc += contribSrcOfs [p -> proc]; /* Update min and max path length */ if (p -> aux.pathLen < pmap -> minPathLen) pmap -> minPathLen = p -> aux.pathLen; else if (p -> aux.pathLen > pmap -> maxPathLen) pmap -> maxPathLen = p -> aux.pathLen; /* Update avg photon path length */ avgPathLen += p -> aux.pathLen; /* Scale photon's flux (hitherto normalised to 1 over RGB); in * case of a contrib photon map, this is done per light source, * and photonFlux is assumed to be an array */ getPhotonFlux(p, flux); if (photonFlux) { scalecolor(flux, photonFlux [isContribPmap(pmap) ? photonSrcIdx(pmap, p) : 0] ); setPhotonFlux(p, flux); } /* Update average photon flux; need a double here */ addcolor(avgFlux, flux); } /* Write modified photons to new heap */ fwrite(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufLen, nuHeap); if (ferror(nuHeap)) { error(SYSTEM, "failed postprocessing photon flux in buildPhotonMap()" ); /* Clean up temp files */ fclose(pmap -> heap); fclose(nuHeap); unlink(pmap -> heapFname); unlink(nuHeapFname); } nCheck += pmap -> heapBufLen; } #ifdef DEBUG_PMAP if (nCheck < pmap -> numPhotons) error(INTERNAL, "truncated photon heap in buildPhotonMap"); #endif /* Finalise average photon flux */ scalecolor(avgFlux, 1.0 / pmap -> numPhotons); copycolor(pmap -> photonFlux, avgFlux); /* Average photon positions to get centre of gravity */ for (i = 0; i < 3; i++) pmap -> CoG [i] = CoG [i] /= pmap -> numPhotons; /* Average photon path lengths */ pmap -> avgPathLen = avgPathLen /= pmap -> numPhotons; rewind(pmap -> heap); /* Compute average photon distance to centre of gravity */ while (!feof(pmap -> heap)) { pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufSize, pmap -> heap ); for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) { VSUB(d, p -> pos, CoG); CoGdist += DOT(d, d); if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) { /* Add temporal distance (in units of photon path length) for transient photons */ d [0] = p -> aux.pathLen - avgPathLen; CoGdist += d [0] * d [0]; } } } pmap -> CoGdist = CoGdist /= pmap -> numPhotons; /* Swap heaps, discarding unscaled photons */ fclose(pmap -> heap); unlink(pmap -> heapFname); pmap -> heap = nuHeap; strcpy(pmap -> heapFname, nuHeapFname); #ifdef PMAP_OOC OOC_BuildPhotonMap(pmap, nproc); #else if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) kdT_TransBuildPhotonMap(pmap); else kdT_BuildPhotonMap(pmap); #endif /* Trash heap and its buffa */ free(pmap -> heapBuf); fclose(pmap -> heap); unlink(pmap -> heapFname); pmap -> heap = NULL; pmap -> heapBuf = NULL; } /* PHOTON MAP SEARCH ROUTINES ------------------------------------------ */ /* Dynamic max photon search radius increase and reduction factors */ #define PMAP_MAXDIST_INC 4 #define PMAP_MAXDIST_DEC 0.9 /* Num successful lookups before reducing in max search radius */ #define PMAP_MAXDIST_CNT 1000 /* Threshold below which we assume increasing max radius won't help */ #define PMAP_SHORT_LOOKUP_THRESH 1 /* Coefficient for adaptive maximum search radius */ #define PMAP_MAXDIST_COEFF 100 void findPhotons (PhotonMap* pmap, const RAY* ray) { int redo = 0; const RREAL *norm = ray -> ron, *pos = ray -> rop; if (!pmap -> squeue.len) { /* Lazy init priority queue */ #ifdef PMAP_OOC OOC_InitFindPhotons(pmap); #else kdT_InitFindPhotons(pmap); #endif pmap -> minGathered = pmap -> maxGather; pmap -> maxGathered = pmap -> minGather; pmap -> totalGathered = 0; pmap -> numLookups = pmap -> numShortLookups = 0; pmap -> shortLookupPct = 0; pmap -> minError = FHUGE; pmap -> maxError = -FHUGE; pmap -> rmsError = 0; /* SQUARED max search radius limit is based on avg photon distance to * centre of gravity, unless fixed by user (maxDistFix > 0) */ pmap -> maxDist0 = pmap -> maxDist2Limit = maxDistFix > 0 ? maxDistFix * maxDistFix : PMAP_MAXDIST_COEFF * pmap -> squeue.len * pmap -> CoGdist / pmap -> numPhotons; } do { pmap -> squeue.tail = 0; pmap -> maxDist2 = pmap -> maxDist0; if (isVolumePmap(pmap) #ifdef PMAP_PHOTONFLOW || isLightFlowPmap(pmap) && sphericalIrrad #endif ) { /* Volume photons are not associated with a normal or intersection point; null the normal and set origin as lookup point. */ /* If a lightflow photon map is used and sphericalIrrad = 1, the mean spherical irradiance is evaluated, so normals are irrelevant and not filtered. */ norm = NULL; pos = ray -> rorg; } /* XXX/NOTE: status returned by XXX_FindPhotons() is currently ignored; if no photons are found, an empty queue is returned under the assumption all photons are too distant to contribute significant flux. */ #ifdef PMAP_OOC OOC_FindPhotons(pmap, pos, norm); #else if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) kdT_TransFindPhotons(pmap, pos, norm); else kdT_FindPhotons(pmap, pos, norm); #endif #ifdef PMAP_LOOKUP_INFO fprintf(stderr, "%d/%d %s photons found within radius %.3f " "at (%.2f,%.2f,%.2f) on %s\n", pmap -> squeue.tail, pmap -> squeue.len, pmapName [pmap -> type], sqrt(pmap -> maxDist2), ray -> rop [0], ray -> rop [1], ray -> rop [2], ray -> ro ? ray -> ro -> oname : "" ); #endif if (pmap -> squeue.tail < pmap -> squeue.len * pmap->gatherTolerance) { /* Short lookup; too few photons found */ if (pmap -> squeue.tail > PMAP_SHORT_LOOKUP_THRESH) { /* Ignore short lookups which return fewer than * PMAP_SHORT_LOOKUP_THRESH photons under the assumption there * really are no photons in the vicinity, and increasing the max * search radius therefore won't help */ #ifdef PMAP_LOOKUP_WARN sprintf(errmsg, "%d/%d %s photons found at (%.2f,%.2f,%.2f) on %s", pmap -> squeue.tail, pmap->squeue.len, pmapName [pmap->type], ray -> rop [0], ray -> rop [1], ray -> rop [2], ray -> ro ? ray -> ro -> oname : "" ); error(WARNING, errmsg); #endif /* Bail out after warning if maxDist is fixed */ if (maxDistFix > 0) return; if (pmap -> maxDist0 < pmap -> maxDist2Limit) { /* Increase max search radius if below limit & redo search */ pmap -> maxDist0 *= PMAP_MAXDIST_INC; #ifdef PMAP_LOOKUP_REDO redo = 1; #endif #ifdef PMAP_LOOKUP_WARN sprintf(errmsg, redo ? "restarting photon lookup with max radius %.1e" : "max photon lookup radius adjusted to %.1e", pmap -> maxDist0 ); error(WARNING, errmsg); #endif } #ifdef PMAP_LOOKUP_REDO else { sprintf(errmsg, "max photon lookup radius clamped to %.1e", pmap -> maxDist0 ); error(WARNING, errmsg); } #endif } /* Reset successful lookup counter */ pmap -> numLookups = 0; } else { /* Skip search radius reduction if maxDist is fixed */ if (maxDistFix > 0) return; /* Increment successful lookup counter and reduce max search radius if * wraparound */ pmap -> numLookups = (pmap -> numLookups + 1) % PMAP_MAXDIST_CNT; if (!pmap -> numLookups) pmap -> maxDist0 *= PMAP_MAXDIST_DEC; redo = 0; } } while (redo); } Photon *find1Photon (PhotonMap *pmap, const RAY* ray, Photon *photon) { /* Init (squared) search radius to avg photon dist to centre of gravity */ float maxDist2_0 = pmap -> CoGdist; int found = 0; #ifdef PMAP_LOOKUP_REDO #define REDO 1 #else #define REDO 0 #endif do { pmap -> maxDist2 = maxDist2_0; #ifdef PMAP_OOC found = OOC_Find1Photon(pmap, ray -> rop, ray -> ron, photon); #else found = kdT_Find1Photon(pmap, ray -> rop, ray -> ron, photon); #endif if (found < 0) { /* Expand search radius to retry */ maxDist2_0 *= 2; #ifdef PMAP_LOOKUP_WARN sprintf(errmsg, "failed 1-NN photon lookup" #ifdef PMAP_LOOKUP_REDO ", retrying with search radius %.2f", maxDist2_0 #endif ); error(WARNING, errmsg); #endif } } while (REDO && found < 0); /* Return photon buffer containing valid photon, else NULL */ return found < 0 ? NULL : photon; } void getPhoton (PhotonMap *pmap, PhotonIdx idx, Photon *photon) { #ifdef PMAP_OOC if (OOC_GetPhoton(pmap, idx, photon)) #else if (kdT_GetPhoton(pmap, idx, photon)) #endif error(INTERNAL, "failed photon lookup"); } Photon *getNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx) { #ifdef PMAP_OOC return OOC_GetNearestPhoton(squeue, idx); #else return kdT_GetNearestPhoton(squeue, idx); #endif } PhotonIdx firstPhoton (const PhotonMap *pmap) { #ifdef PMAP_OOC return OOC_FirstPhoton(pmap); #else return kdT_FirstPhoton(pmap); #endif } /* PHOTON MAP CLEANUP ROUTINES ------------------------------------------ */ void deletePhotons (PhotonMap* pmap) { #ifdef PMAP_OOC OOC_Delete(&pmap -> store); #else kdT_Delete(&pmap -> store); #endif free(pmap -> squeue.node); free(pmap -> biasCompHist); pmap -> numPhotons = pmap -> minGather = pmap -> maxGather = pmap -> squeue.len = pmap -> squeue.tail = 0; #ifdef PMAP_PATHFILT if (pmap -> pathLUT) { lu_done(pmap -> pathLUT); free(pmap -> pathLUT); } if (pmap -> pathLUTKeys) { unsigned k; for (k = 0; k < pmap->squeue.len + 1; free(pmap->pathLUTKeys [k++])); free(pmap -> pathLUTKeys); } #endif } diff --git a/pmapdens.c b/pmapdens.c index 6d60966..9ca0cce 100644 --- a/pmapdens.c +++ b/pmapdens.c @@ -1,319 +1,320 @@ #ifndef lint static const char RCSid[] = "$Id$"; #endif /* ====================================================================== Photon map density estimation routines 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$ */ #include "pmapdens.h" #include "pmapbias.h" #include "otypes.h" #include "random.h" /* Photon map lookup functions per type */ void (*pmapLookup [NUM_PMAP_TYPES])(PhotonMap*, RAY*, COLOR) = { photonDensity, photonPreCompDensity, photonDensity, volumePhotonDensity, photonDensity, photonDensity, photonDensity #ifdef PMAP_PHOTONFLOW , lightFlowDensity, lightFlowDensity #endif }; void photonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad) /* Surface-bound photon density estimate for global and caustic photons. Returns irradiance at ray -> rop. */ { unsigned i; float r2; COLOR flux; Photon *photon; const PhotonSearchQueueNode *sqn; setcolor(irrad, 0, 0, 0); if (!pmap -> maxGather) return; /* Ignore sources */ if (ray -> ro && islight(objptr(ray -> ro -> omod) -> otype)) return; findPhotons(pmap, ray); /* Need at least 2 photons */ if (pmap -> squeue.tail < 2) { #ifdef PMAP_NONEFOUND sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)", ray -> ro ? ray -> ro -> oname : "", ray -> rop [0], ray -> rop [1], ray -> rop [2] ); error(WARNING, errmsg); #endif return; } if (pmap -> minGather == pmap -> maxGather) { /* No bias compensation. Just do a plain vanilla estimate */ sqn = pmap -> squeue.node + 1; /* Average radius^2 between furthest two photons to improve accuracy */ r2 = max(sqn -> dist2, (sqn + 1) -> dist2); r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2)); /* Skip the extra photon */ for (i = 1 ; i < pmap -> squeue.tail; i++, sqn++) { photon = getNearestPhoton(&pmap -> squeue, sqn -> idx); getPhotonFlux(photon, flux); #ifdef PMAP_EPANECHNIKOV /* Apply Epanechnikov kernel to photon flux based on distance. (see Eq. 4.1, p.76 in Silverman, "Density Estimation for Statistics and Data Analysis", 1st Ed., 1986, and Wann Jensen, "Realistic Image Synthesis using Photon Mapping") */ scalecolor(flux, 1 - sqn -> dist2 / r2); #endif addcolor(irrad, flux); } #ifdef PMAP_EPANECHNIKOV /* Normalise accumulated flux by Epanechnikov kernel integral in 2D (see Eq. 4.4, p.76 in Silverman, "Density Estimation for Statistics and Data Analysis", 1st Ed., 1986, and Wann Jensen, "Realistic Image Synthesis using Photon Mapping"), include RADIANCE-specific lambertian factor PI */ scalecolor(irrad, 2 / (PI * PI * r2)); #else /* Normalise accumulated flux by search area PI * r^2, including RADIANCE-specific lambertian factor PI */ scalecolor(irrad, 1 / (PI * PI * r2)); #endif return; } else /* Apply bias compensation to density estimate */ biasComp(pmap, irrad); } void photonPreCompDensity (PhotonMap *pmap, RAY *r, COLOR irrad) /* Surface-bound photon density estimate for (single) precomputed photon. Returns precomputed irradiance at ray -> rop. */ { Photon p; setcolor(irrad, 0, 0, 0); /* Ignore sources */ if (r -> ro && islight(objptr(r -> ro -> omod) -> otype)) return; if (find1Photon(preCompPmap, r, &p)) /* p contains a found photon, so get its irradiance, otherwise it * remains zero under the assumption all photons are too distant * to contribute significantly */ getPhotonFlux(&p, irrad); } void volumePhotonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad) /* Volume photon density estimate using Henyey-Greenstein phase function. Returns inscattered irradiance at ray -> rop. */ { unsigned i; float r2, gecc2, ph; COLOR flux; Photon *photon; const PhotonSearchQueueNode *sqn; setcolor(irrad, 0, 0, 0); if (!pmap -> maxGather) return; findPhotons(pmap, ray); /* Need at least 2 photons */ if (pmap -> squeue.tail < 2) return; gecc2 = ray -> gecc * ray -> gecc; sqn = pmap -> squeue.node + 1; /* Average radius^2 between furthest two photons to improve accuracy */ r2 = max(sqn -> dist2, (sqn + 1) -> dist2); r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2)); /* Skip the extra photon */ for (i = 1; i < pmap -> squeue.tail; i++, sqn++) { photon = getNearestPhoton(&pmap -> squeue, sqn -> idx); /* Compute phase function for inscattering from photon */ if (gecc2 <= FTINY) ph = 1; else { ph = DOT(ray -> rdir, photon -> norm) / 127; ph = 1 + gecc2 - 2 * ray -> gecc * ph; ph = (1 - gecc2) / (ph * sqrt(ph)); } getPhotonFlux(photon, flux); scalecolor(flux, ph); addcolor(irrad, flux); } /* Divide by search volume 4 / 3 * PI * r^3 and phase function normalization factor 4 * PI */ scalecolor(irrad, 3 / (16 * PI * PI * r2 * sqrt(r2))); } #ifdef PMAP_PHOTONFLOW void lightFlowDensity (PhotonMap *pmap, RAY *ray, COLOR irrad) /* Evaluate (static) irradiance from volume photons incident on plane defined by normal ray -> ron at point ray -> rop, ignoring participating medium. This evaluation interprets the volume photon map as a representation of the "flow" of light, i.e. a physical light field. If sphericalIrrad == 0, the planar irradiance is evaluated with respect to the plane defined by the ray normal ray -> ron. If sphericalIrrad == 1, the mean spherical irradiance is evaluated, - and the normal ray -> ron is ignored. */ + and the normal ray -> ron is ignored. + If hemiSearch == 0, photons are searched in a spherical volume. + If hemiSearch == 1, photons are */ { unsigned p, i; float r2, cosNorm; Photon *photon; const PhotonSearchQueueNode *sqn; COLOR photonFlux; FVECT photonDir; setcolor(irrad, 0, 0, 0); if (!pmap -> maxGather) return; /* Flip normal since volume photon directions point towards plane of incidence (has no effect if evaluating mean spherical irradiance) */ flipsurface(ray); /* Find photons with incident direction filtered by dot product with normal (note that found photons may also be located _behind_ the plane, as if having passed through it. */ findPhotons(pmap, ray); if (pmap -> squeue.tail < pmap -> squeue.len && !isTransLightFlowPmap(pmap) ) { sprintf(errmsg, "short lookup in photon flow; consider -am %.3f or greater", 2 * sqrt(pmap -> maxDist2Limit) ); error(WARNING, errmsg); } /* Need at least 2 photons, else bail out */ if (pmap -> squeue.tail < 2) return; sqn = pmap -> squeue.node + 1; /* Avg radius^2 between furthest two photons to improve accuracy */ r2 = max(sqn -> dist2, (sqn + 1) -> dist2); r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2)); #ifdef PMAP_DEBUGPATHS printf("N = %d\tr = %f\n", pmap -> squeue.tail, sqrt(r2)); #endif /* Skip the extra photon */ for (p = 1; p < pmap -> squeue.tail; p++, sqn++) { photon = getNearestPhoton(&pmap -> squeue, sqn -> idx); getPhotonFlux(photon, photonFlux); #ifdef PMAP_EPANECHNIKOV /* Apply Epanechnikov kernel to photon flux based on distance (see Eq. 4.1, p.76 in Silverman, "Density Estimation for Statistics and Data Analysis", 1st Ed., 1986) */ scalecolor(photonFlux, 1 - sqn -> dist2 / r2); #endif if (!sphericalIrrad) { /* Associate photons with plane defined by ray normal */ for (i = 0; i < 3; i++) { /* Jitter photon direction since stored as 8-bit signed int */ photonDir [i] = photon -> norm [i]; if (photon -> norm [i] && photon -> norm [i] & 127 < 127) photonDir [i] += photon -> norm [i] & 0x80 ? -frandom() : frandom(); } /* Project photon direction onto normal */ cosNorm = DOT(ray -> ron, photonDir) / 127; scalecolor(photonFlux, cosNorm); } addcolor(irrad, photonFlux); #ifdef PMAP_DEBUGPATHS printf("%f\t%f\t%f\t%f\t%f\t%d:%010d\n", sqrt(sqn -> dist2), photon -> pos [0], photon -> pos [1], photon -> pos [2], photonFlux [0], photon -> proc, photon -> primary ); #endif } #ifdef PMAP_EPANECHNIKOV /* Normalise accumulated flux by Epanechnikov kernel integral in 3D (see Eq. 4.4, p.76 in Silverman, "Density Estimation for Statistics and Data Analysis", 1st Ed., 1986), include RADIANCE-specific lambertian factor PI */ scalecolor(irrad, 15 / (8 * PI * PI * r2 * sqrt(r2))); #else /* Normalise accumulated flux by spherical volume to obtain irradiance, include RADIANCE-specific irradiance factor PI */ scalecolor(irrad, 3 / (4 * PI * PI * r2 * sqrt(r2))); #endif - if (sphericalIrrad) - /* Scale irradiance by 1/4; this corresponds to the "true" scalar - illuminance as defined in Eq. 13 in: - Mangkuto R.A., "Research note: The accuracy of the mean spherical - semi-cubic illuminance approach for determining scalar illuminance", - Lighting Research & Technology, 2020; 52: 151-158. */ - scalecolor(irrad, 0.25); - #ifdef PMAP_PHOTONFLOW_HEMI - else + if (hemiLightFlow) /* Hemispherical mode; lightflow photons which lie behind the plane were filtered out (see filterPhoton()), so renormalise for hemisphere (3/(2 PI)). This reduces bias near "hard" boundaries (scene cube limits or actual geometry with opaque surfaces), in which case no photons may be present behind the plane. */ scalecolor(irrad, 2.0); - #endif + else if (sphericalIrrad) + /* Mean spherical irradian mode; scale irradiance by 1/4. This + corresponds to the "true" scalar illuminance as defined in + Eq. 13 in: Mangkuto R.A., "Research note: The accuracy of the + mean spherical semi-cubic illuminance approach for determining + scalar illuminance", Lighting Research & Technology, 2020; 52: + 151-158. */ + scalecolor(irrad, 0.25); /* Restore ray normal */ flipsurface(ray); } #endif /* PMAP_PHOTONFLOW */ diff --git a/pmapfilt.c b/pmapfilt.c index 832833c..3fbadcb 100644 --- a/pmapfilt.c +++ b/pmapfilt.c @@ -1,359 +1,358 @@ /* ====================================================================== Photon map filter callbacks for nearest neighbour search Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ====================================================================== $Id: pmapfilt.c,v 1.2 2020/11/16 22:28:23 u-no-hoo Exp u-no-hoo $ */ #include "pmapfilt.h" -/* pmap.h needed for PMAP_PHOTONFLOW_HEMI; also inludes pmapdata.h */ +#include "pmapdata.h" #include "pmap.h" #include "random.h" #include "source.h" #include "otspecial.h" int filterPhoton (const void *rec, const void *state) /* Filter callback for photon kNN search */ { const Photon *photon = rec; const PhotonSearchFilterState *filtState = state; const PhotonMap *pmap = filtState -> pmap; /* Reject photon if normal faces away (ignored for volume photons) with * tolerance to account for perturbation; note photon normal is coded * in range [-127,127], hence we factor this in */ if (filtState -> norm && DOT(filtState->norm, photon->norm) <= PMAP_NORM_TOL * 127 * frandom() ) return 0; -#if defined(PMAP_PHOTONFLOW) && defined(PMAP_PHOTONFLOW_HEMI) - /* Omit lightflow photons which pass normal test above but lie behind - the plane (= hemispherical instead of spherical lookup). - This reduces bias near solid boundaries (scene cube limits or - actual geometry), in which case no photons may be present behind - the plane. */ - /* NOTE: filtstate->norm is REVERSED in this mode since photon->norm +#ifdef PMAP_PHOTONFLOW + if (isLightFlowPmap(pmap) && hemiLightFlow && !sphericalIrrad) { + /* Hemispherical lookup mode; omit lightflow photons which pass normal + test above but lie behind the plane. This reduces bias near solid + boundaries (scene cube limits or actual geometry), in which case no + photons may be present behind the plane. + NOTE: filtstate->norm is REVERSED in this mode since photon->norm (=photon direction here) points TOWARDS the receiving surface! */ - if (isLightFlowPmap(pmap) && !sphericalIrrad) { FVECT photonVec; VSUB(photonVec, photon -> pos, filtState -> pos); if (DOT(filtState -> norm, photonVec) > FTINY) return 0; } #endif if (isContribPmap(pmap)) { /* Lookup in contribution photon map; accept photon only if its emitting light source modifier matches that of the filter state (if defined). */ if (filtState -> srcMod) { const int srcIdx = photonSrcIdx(pmap, photon); if (srcIdx < 0 || srcIdx >= nsources) error(INTERNAL, "invalid light source index in photon map"); if (filtState -> srcMod != findmaterial(source [srcIdx].so)) return 0; } /* Reject non-caustic photon if lookup for caustic contribs */ if (pmap -> lookupCaustic && !photon -> caustic) return 0; } /* Accept photon */ return 1; } #ifdef PMAP_PATHFILT /* Length of key strings for photon path LUT (plus terminator) */ #define PMAP_PATHKEYLEN ((sizeof(PhotonOrigin) + 1 << 1) + 1) void initPhotonPaths (void *state) /* Lazily allocate and initialise photon path ID lookup table and * associate key buffer */ { PhotonMap *pmap = state; if (!pmap -> pathLUT) { /* Allocate and init photon path ID lookup table of same size as * search queue */ LUTAB initLUT = LU_SINIT(NULL, NULL); if (!(pmap -> pathLUT = malloc(sizeof(LUTAB)))) error(SYSTEM, "failed photon path lookup table allocation"); memcpy(pmap -> pathLUT, &initLUT, sizeof(initLUT)); if (!lu_init(pmap -> pathLUT, pmap -> squeue.len)) error(INTERNAL, "failed photon path lookup table init"); } if (!pmap -> pathLUTKeys) { /* Allocate key buffer of same size as search queue (+1 for * scratchpad at tail) */ unsigned k, maxKeys = pmap -> squeue.len + 1; pmap -> pathLUTKeys = (char**)calloc(maxKeys, sizeof(char*)); for (k = 0; k < maxKeys; k++) if ( pmap -> pathLUTKeys || !(pmap -> pathLUTKeys [k] = (char*)malloc(PMAP_PATHKEYLEN)) ) error(SYSTEM, "failed photon path LUT key buffer allocation"); } /* Reset photon path ID key buffer */ pmap -> numPathLUTKeys = 0; } static char *photonPathKey (const Photon *photon, char *key) { /* Generate LUT key from photon path ID. Writes transposed hex string into specified key buffer. */ PhotonOrigin pathId1 = photon -> org; unsigned char pathId2 = photon -> proc; char *kp = key; unsigned i; for (i = sizeof(PhotonOrigin) << 1; i; *kp++ = (pathId1 & 0x0f) + 'A', pathId1 >>= 4, i-- ); for (i = sizeof(unsigned char) << 1; i; *kp++ = (pathId2 & 0x0f) + 'A', pathId2 >>= 4, i-- ); *kp++ = '\0'; return key; } int checkPhotonPaths (const void *state) /* Run sanity check on photon path ID lookup table */ { const PhotonMap *pmap = state; const PhotonSearchQueue *squeue = &pmap -> squeue; const PhotonSearchQueueNode *squeueNode; const Photon *photon; const LUENT *lutEntry; unsigned squeueIdx; char key [PMAP_PATHKEYLEN]; if (!pmap -> pathLUT || !pmap -> pathLUTKeys) error(INTERNAL, "photon path lookup table not initialised"); for (squeueIdx = 0, squeueNode = squeue -> node; squeueIdx < squeue -> tail; squeueIdx++, squeueNode++ ) { photon = getNearestPhoton(squeue, squeueNode -> idx); photonPathKey(photon, key); lutEntry = lu_find(pmap -> pathLUT, key); if (!lutEntry || !lutEntry -> key || strcmp(lutEntry -> key, key)) error(CONSISTENCY, "invalid key in photon path LUT check"); if (!lutEntry -> data || (PhotonSearchQueueNode*)lutEntry -> data != squeueNode ) error(CONSISTENCY, "invalid data in photon path LUT check"); } return 0; } /* TODO: REPLACE pathLUTKeys SCRATCHPAD WITH LOCAL ARRAYS??? */ void **findPhotonPath (const void *rec, const void *state) /* Photon path callback for k-NN search to map photon path IDs to search * queue nodes via a lookup table. State points to the associated photon * map containing the path lookup table. Returns modifiable pointer to data * field of found LUT entry, which in turn points to the search queue node * for the corresponding record. If no entry exists for the path ID, NULL * is returned. */ { const Photon *photon = rec; const PhotonMap *pmap = state; char *key; LUENT *lutEntry; if (!pmap -> pathLUT || !pmap -> pathLUTKeys) error(INTERNAL, "photon path lookup table not initialised"); /* Generate key for photon's path and write into next free slot in key buffer (using it as scratchpad) */ key = photonPathKey(photon, pmap -> pathLUTKeys [pmap -> numPathLUTKeys] ); #ifdef PMAP_DEBUGPATHS printf("%s ", key); #endif /* Look up attribute for generated key, set last LUT entry in state */ if (!(lutEntry = lu_find(pmap -> pathLUT, key))) /* Outta mammaries */ error(SYSTEM, "failed photon path lookup entry allocation"); return (void**)(lutEntry -> data ? &lutEntry -> data : NULL); } void addPhotonPath (const void *rec, void *data, void *state) /* Photon path callback for k-NN search to add photon path IDs to lookup table entries. State points to the associated photon map containing the path lookup table. The new lookup table entry is initialised with the specified data pointer */ { const Photon *photon = rec; PhotonMap *pmap = state; char *key; LUENT *lutEntry; if (!pmap -> pathLUT || !pmap -> pathLUTKeys) error(INTERNAL, "photon path lookup table not initialised"); /* Generate key for photon's path and write into next free slot in key buffer */ key = photonPathKey(photon, pmap -> pathLUTKeys [pmap -> numPathLUTKeys] ); /* Find free LUT entry for generated key */ if (!(lutEntry = lu_find(pmap -> pathLUT, key))) /* Outta mammaries */ error(SYSTEM, "failed photon path lookup entry allocation"); if (lutEntry -> data) /* Occupied */ error(INTERNAL, "attempt to overwrite photon path lookup entry"); #ifdef PMAP_DEBUGPATHS printf("%s add\n", key); #endif /* Set data and key fields, advance key buffer tail */ lutEntry -> data = data; lutEntry -> key = pmap -> pathLUTKeys [pmap -> numPathLUTKeys]; if (pmap -> numPathLUTKeys > pmap -> squeue.tail) /* Number of keys (and distinct path IDs) can never exceed the number of photons currently in the search queue (+1 for scratchpad) */ error(INTERNAL, "photon path lookup table key overflow"); else pmap -> numPathLUTKeys++; } static char emptyKey = '\0'; void deletePhotonPath (const void *rec, void *state) /* Photon path callback for k-NN search to delete photon path IDs from lookup table. State points to the associated photon map containing the path ID lookup table. */ { const Photon *photon = rec; PhotonMap *pmap = state; LUENT *lutEntry; char *delKey, *relocKey; if (!pmap -> pathLUT || !pmap -> pathLUTKeys) error(INTERNAL, "photon path lookup table not initialised"); /* Generate key for photon's path and write into scratchpad at key buffer tail (will be overwritten) */ delKey = photonPathKey(photon, pmap -> pathLUTKeys [pmap -> numPathLUTKeys] ); #ifdef PMAP_DEBUGPATHS printf("%s del\n", delKey); #endif /* Find photon's LUT entry, check if valid, and get its key buffer slot */ lutEntry = lu_find(pmap -> pathLUT, delKey); if (!lutEntry || !lutEntry -> key || !lutEntry -> data) error(CONSISTENCY, "attempt to delete empty photon path LUT entry"); /* Delete photon LUT entry */ lu_delete(pmap -> pathLUT, delKey); delKey = lutEntry -> key; lutEntry -> key = &emptyKey; /* ??? */ /* Decrement key buffer tail and get last valid key */ relocKey = pmap -> pathLUTKeys [--pmap -> numPathLUTKeys]; if (delKey != relocKey) { /* Reclaim deleted photon's key buffer slot by overwriting with last valid key */ memcpy(delKey, relocKey, PMAP_PATHKEYLEN); relocKey = delKey; #ifdef PMAP_DEBUGPATHS printf("%s reloc\n", relocKey); #endif /* Update the LUT entry owning the relocated key to point to its new position in the overwritten slot, again checking for validity */ lutEntry = lu_find(pmap -> pathLUT, relocKey); if (!lutEntry || !lutEntry -> key || !lutEntry -> data) error(CONSISTENCY, "attempt to delete empty photon path LUT entry" ); lutEntry -> key = relocKey; } } int clearPhotonPath (const LUENT *lutEntry, void *lut) /* Clear entry in photon path ID lookup table */ { #ifdef PMAP_DEBUGPATHS printf("%s clr\n", lutEntry -> key); #endif lu_delete(lut, lutEntry -> key); /* Apologies to Don Gregorio for abusing his LUT here... */ ((LUENT*)lutEntry) -> key = &emptyKey; return 0; } void resetPhotonPaths (void *state) /* Reset (clear) all entries in photon path ID lookup table */ { PhotonMap *pmap = state; if (!pmap -> pathLUT) error(INTERNAL, "photon path lookup table not initialised"); lu_doall(pmap -> pathLUT, clearPhotonPath, pmap -> pathLUT); /* Reset photon path ID key buffer */ pmap -> numPathLUTKeys = 0; } #endif /* PMAP_PATHFILT */ diff --git a/pmapopt.c b/pmapopt.c index e29799d..e60941e 100644 --- a/pmapopt.c +++ b/pmapopt.c @@ -1,173 +1,192 @@ #ifndef lint static const char RCSid[] = "$Id: pmapopt.c,v 2.10 2020/06/15 22:18:57 rschregle Exp $"; #endif /* ================================================================== Photon map interface to RADIANCE render options 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: pmapopt.c,v 2.10 2020/06/15 22:18:57 rschregle Exp $ */ #include "ray.h" #include "pmapparm.h" int getPmapRenderOpt (int ac, char *av []) /* Parse next render option for photon map; interface to getrenderopt(); * return -1 if parsing failed, else number of parameters consumed */ { #define check(ol,al) (av[0][ol] || badarg(ac-1,av+1,al)) /* Evaluate boolean option, setting var accordingly */ #define check_bool(olen, var) switch (av [0][olen]) { \ case '\0': \ var = !var; break; \ case 'y': case 'Y': case 't': case 'T': case '+': case '1': \ var = 1; break; \ case 'n': case 'N': case 'f': case 'F': case '-': case '0': \ var = 0; break; \ default: \ return -1; \ } static int t = -1; /* pmap parameter index */ if (ac < 1 || !av [0] || av [0][0] != '-') return -1; switch (av [0][1]) { case 'a': switch (av [0][2]) { case 'p': /* photon map */ /* Photon map bumps up ambounce from its default 0; SHOULD THIS REMAIN THE DEFAULT BEHAVIOUR? In some cases it's desirable to render photons directly via ab -1 */ ambounce += (ambounce == 0); if (!check(3, "s")) { /* File -> assume bwidth = 1 or precomputed pmap */ if (++t >= NUM_PMAP_TYPES) error(USER, "too many photon maps"); pmapParams [t].fileName = savqstr(av [1]); pmapParams [t].minGather = pmapParams [t].maxGather = defaultGather; } else return -1; if (!check(3, "si")) { /* File + bandwidth */ pmapParams [t].minGather = pmapParams [t].maxGather = atoi(av [2]); if (!pmapParams [t].minGather) return -1; } else { sprintf(errmsg, "missing photon lookup bandwidth, defaulting to %d", defaultGather ); error(WARNING, errmsg); return 1; } if (!check(3, "sii")) { /* File + min bwidth + max bwidth -> bias compensation */ pmapParams [t].maxGather = atoi(av [3]); if (pmapParams [t].minGather >= pmapParams [t].maxGather) return -1; return 3; } /* Transient pmap currently only supported by kd-tree */ #ifndef PMAP_OOC if (!check(3, "sif")) { /* File + bwidth + time -> transient pmap */ pmapParams [t].time = atof(av [3]); if (pmapParams [t].time < 0) error(USER, "invalid time %g -- " "we're not Steven bloody Hawking, you know!" ); return 3; } #endif else return 2; case 'm': /* Fixed max photon search radius */ if (check(3, "f") || (maxDistFix = atof(av [1])) <= 0) error(USER, "invalid max photon search radius"); return 1; #ifdef PMAP_OOC case 'c': /* OOC pmap cache page size ratio */ if (check(3, "f") || (pmapCachePageSize = atof(av [1])) <= 0) error(USER, "invalid photon cache page size ratio"); return 1; case 'C': /* OOC pmap cache size (in photons); 0 = no caching */ if (check(3, "s")) error(USER, "invalid number of cached photons"); /* Parsing failure sets to zero and disables caching */ pmapCacheSize = parseMultiplier(av [1]); return 1; #endif #ifdef PMAP_PHOTONFLOW + case 'H': + /* Hemispherical lookup switch for lightflow pmap */ + check_bool(3, hemiLightFlow); + if (hemiLightFlow && sphericalIrrad) + error(USER, + "hemispherical lookups and mean spherical irradiance " + "are mutually exclusive" + ); + else return 0; + case 'S': /* Mean spherical / planar irrad switch for lightflow pmap */ check_bool(3, sphericalIrrad); - return 0; + if (sphericalIrrad && hemiLightFlow) + error(USER, + "spherical irradiance and hemispherical lookups " + "are mutually exclusive" + ); + else return 0; #endif } } #undef check /* Unknown option */ return -1; } void printPmapDefaults () /* Print defaults for photon map render options */ { printf("-am %.1f\t\t\t\t# max photon search radius\n", maxDistFix); #ifdef PMAP_OOC printf("-ac %.1f\t\t\t\t# photon cache page size ratio\n", pmapCachePageSize); printf("-aC %ld\t\t\t# num cached photons\n", pmapCacheSize); #endif #ifdef PMAP_PHOTONFLOW + printf(hemiLightFlow + ? "-aH+\t\t\t\t# Hemispherical lookups for lightflow photons\n" + : "-aH-\t\t\t\t# Spherical lookups for lightflow photons\n" + ); printf(sphericalIrrad ? "-aS+\t\t\t\t# spherical irradiance from lightflow photons\n" : "-aS-\t\t\t\t# planar irradiance from lightflow photons\n" ); #endif } diff --git a/pmapparm.c b/pmapparm.c index 88cd396..9dd7789 100644 --- a/pmapparm.c +++ b/pmapparm.c @@ -1,129 +1,131 @@ #ifndef lint static const char RCSid[] = "$Id: pmapparm.c,v 2.9 2018/02/02 19:47:55 rschregle Exp $"; #endif /* ====================================================================== Parameters for photon map generation and rendering; used by mkpmap and rpict/rvu/rtrace. 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", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ====================================================================== $Id: pmapparm.c,v 2.9 2018/02/02 19:47:55 rschregle Exp $ */ #include "pmapparm.h" #include "pmapdata.h" #include float pdfSamples = 1000, /* PDF samples per steradian */ finalGather = 0.25, /* fraction of global photons for irradiance precomputation */ preDistrib = 0.25, /* fraction of num photons for distribution prepass */ gatherTolerance = 0.5, /* Photon map lookup tolerance; lookups returning fewer than this fraction of minGather/maxGather are restarted with a larger search radius */ maxDistFix = 0, /* Static maximum photon search radius (radius is adaptive if this is zero) */ compressRatio = 0.9, /* Compression ratio for precomputed contribution photon map */ photonMaxPathDist = 0; /* Maximum cumulative distance of photon path */ double photonVelocity = 299792458.0; /* Photon velocity in units of scene geometry (default m/s) */ #ifdef PMAP_OOC float pmapCachePageSize = 8; /* OOC cache pagesize as multiple * of maxGather */ unsigned long pmapCacheSize = 1e6; /* OOC cache size in photons */ #endif #ifdef PMAP_PHOTONFLOW - int sphericalIrrad = 0; /* Toggle for spherical / planar + int sphericalIrrad = 0, /* Toggle for spherical / planar irrad from lightflow pmap */ + hemiLightFlow = 0; /* Toggle hemispherical/spherical + lookups for lightflow pmap */ #endif /* Regions of interest */ unsigned pmapNumROI = 0; PhotonMapROI *pmapROI = NULL; unsigned verbose = 0; /* Verbose console output */ unsigned long photonMaxBounce = 5000; /* Runaway photon bounce limit */ unsigned photonRepTime = 0, /* Seconds between reports */ maxPreDistrib = 4, /* Max predistrib passes */ defaultGather = 40; /* Default num photons for lookup */ /* Per photon map lookup params */ PhotonMapParams pmapParams [NUM_PMAP_TYPES] = { {NULL, 0, 0, 0, 0.}, {NULL, 0, 0, 0, 0.}, {NULL, 0, 0, 0, 0.}, {NULL, 0, 0, 0, 0.}, {NULL, 0, 0, 0, 0.}, {NULL, 0, 0, 0, 0.}, {NULL, 0, 0, 0, 0.} #ifdef PMAP_PHOTONFLOW , {NULL, 0, 0, 0}, {NULL, 0, 0, 0} #endif }; int setPmapParam (PhotonMap** pm, const PhotonMapParams* parm) { if (parm && parm -> fileName) { if (!(*pm = (PhotonMap*)malloc(sizeof(PhotonMap)))) error(INTERNAL, "failed to allocate photon map"); (*pm) -> fileName = parm -> fileName; (*pm) -> minGather = parm -> minGather; (*pm) -> maxGather = parm -> maxGather; (*pm) -> distribTarget = parm -> distribTarget; (*pm) -> maxDist0 = FHUGE; (*pm) -> contribTab = NULL; (*pm) -> time = parm -> time; return 1; } return 0; } unsigned long parseMultiplier (const char *num) { unsigned long mult = 1; const int strEnd = strlen(num) - 1; if (strEnd <= 0) return 0; if (!isdigit(num [strEnd])) switch (toupper(num [strEnd])) { case 'G': mult *= 1000; case 'M': mult *= 1000; case 'K': mult *= 1000; break; default : error(USER, "unknown multiplier"); } return (unsigned long)(mult * atof(num)); } diff --git a/pmapparm.h b/pmapparm.h index cb3d0a7..5bcd48a 100644 --- a/pmapparm.h +++ b/pmapparm.h @@ -1,95 +1,95 @@ /* RCSid $Id: pmapparm.h,v 2.10 2021/04/14 11:26:25 rschregle Exp $ */ /* ====================================================================== Parameters for photon map generation and rendering; used by mkpmap and rpict/rvu/rtrace. 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", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ====================================================================== $Id: pmapparm.h,v 2.10 2021/04/14 11:26:25 rschregle Exp $ */ #ifndef PMAPPARAMS_H #define PMAPPARAMS_H #include "pmaptype.h" /* Struct for passing lookup params per photon map for rpict/rtrace/rvu */ typedef struct { char *fileName; /* Photon map file */ unsigned minGather, maxGather; /* Num photons to gather */ unsigned long distribTarget; /* Num photons to store */ double time; /* Transient photon time */ } PhotonMapParams; /* Region of interest */ typedef struct { /* siz [1], siz [2] <= 0 --> sphere, else rectangle */ float pos [3], siz [3]; } PhotonMapROI; #define PMAP_ROI_ISSPHERE(roi) ((roi)->siz[1] <= 0 && (roi)->siz[2] <= 0) #define PMAP_ROI_SETSPHERE(roi) ((roi)->siz[1] = (roi)->siz[2] = -1) extern PhotonMapParams pmapParams [NUM_PMAP_TYPES]; /* Macros for type specific photon map parameters */ #define globalPmapParams (pmapParams [PMAP_TYPE_GLOBAL]) #define preCompPmapParams (pmapParams [PMAP_TYPE_PRECOMP]) #define causticPmapParams (pmapParams [PMAP_TYPE_CAUSTIC]) #define volumePmapParams (pmapParams [PMAP_TYPE_VOLUME]) #define directPmapParams (pmapParams [PMAP_TYPE_DIRECT]) #define contribPmapParams (pmapParams [PMAP_TYPE_CONTRIB]) #define transientPmapParams (pmapParams [PMAP_TYPE_TRANSIENT]) #ifdef PMAP_PHOTONFLOW #define lightFlowParams (pmapParams [PMAP_TYPE_LIGHTFLOW]) #define transLightFlowParams (pmapParams [PMAP_TYPE_TRANSLIGHTFLOW]) #endif extern float pdfSamples, preDistrib, finalGather, gatherTolerance, maxDistFix, pmapMaxDist, photonMaxPathDist, compressRatio; extern double photonVelocity; extern unsigned long photonHeapSizeInc, photonMaxBounce; extern unsigned photonRepTime, maxPreDistrib, defaultGather, verbose; extern unsigned pmapNumROI; extern PhotonMapROI *pmapROI; #ifdef PMAP_OOC extern float pmapCachePageSize; extern unsigned long pmapCacheSize; #endif #ifdef PMAP_PHOTONFLOW - extern int sphericalIrrad; + extern int sphericalIrrad, hemiLightFlow; #endif struct PhotonMap; int setPmapParam (struct PhotonMap **pm, const PhotonMapParams *parm); /* Allocate photon map and set its parameters from parm */ unsigned long parseMultiplier (const char *num); /* Evaluate numeric parameter string with optional multiplier suffix (G = 10^9, M = 10^6, K = 10^3). Returns 0 if parsing fails. */ #endif