diff --git a/mkpmap.c b/mkpmap.c
index 76661c9..b210ae3 100644
--- a/mkpmap.c
+++ b/mkpmap.c
@@ -1,876 +1,879 @@
#ifndef lint
static const char RCSid[] = "$Id: mkpmap.c,v 2.11 2021/04/14 11:26:25 rschregle Exp $";
#endif
/*
=========================================================================
Photon map generator
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: mkpmap.c,v 2.11 2021/04/14 11:26:25 rschregle Exp $
*/
#include "pmap.h"
#include "pmapmat.h"
#include "pmapsrc.h"
#include "pmapcontrib.h"
#include "pmaprand.h"
#include "paths.h"
#include "ambient.h"
#include "resolu.h"
#include "source.h"
#include <ctype.h>
#include <string.h>
#include <sys/stat.h>
/* Enable options for Ze Ekspertz only! */
#define PMAP_EKSPERTZ
extern char VersionID [];
char* progname; /* argv[0] */
int dimlist [MAXDIM]; /* sampling dimensions */
int ndims = 0; /* number of sampling dimensions */
char* octname = NULL; /* octree name */
CUBE thescene; /* scene top-level octree */
OBJECT nsceneobjs; /* number of objects in scene */
double srcsizerat = 0.01; /* source partition size ratio */
int backvis = 1; /* back face visibility */
int clobber = 0; /* overwrite output */
COLOR cextinction = BLKCOLOR; /* global extinction coefficient */
COLOR salbedo = BLKCOLOR; /* global scattering albedo */
double seccg = 0; /* global scattering eccentricity */
char *amblist [AMBLLEN + 1]; /* ambient include/exclude list */
int ambincl = -1; /* include == 1, exclude == 0 */
char *diagFile = NULL; /* diagnostics output file */
int rand_samp = 1; /* uncorrelated random sampling */
int nproc = 1; /* number of parallel processes */
#ifdef EVALDRC_HACK
char *angsrcfile = NULL; /* angular source file for EvalDRC */
#endif
/* Contribution photon map stuff: light source modifier lookup table and
associated cleanup func */
static void chickenMcFreeNuggetz(void *contrib)
{
epfree((*(MODCONT *)contrib).binv);
free(contrib);
}
LUTAB srcContrib = LU_SINIT(NULL, chickenMcFreeNuggetz);
#if 0
char srcMod [MAXMODLIST];
#endif
unsigned numSrcContrib = 0;
/* The variables below interface with the RADIANCE suite merely as dummies
to resolve external references and make the linker happy; most of them
are unused. */
COLOR ambval = BLKCOLOR;
double shadthresh = .05, ambacc = 0.2, shadcert = .5, minweight = 5e-3,
ssampdist = 0, dstrsrc = 0.0, specthresh = 0.15, specjitter = 1.0,
avgrefl = 0.5;
int ambvwt = 0, ambssamp = 0, ambres = 32, ambounce = 0, directrelay = 1,
directvis = 1, samplendx, do_irrad = 0, ambdiv = 128, vspretest = 512,
maxdepth = 6;
char *shm_boundary = NULL, *ambfile = NULL,
RCCONTEXT [] = "RC."; /* Evaluation context for contrib pmap */
void (*trace)() = NULL, (*addobjnotify [])() = {ambnotify, NULL};
void printdefaults()
/* print default values to stdout */
{
#ifdef EVALDRC_HACK
/* EvalDRC support */
puts("-A\t\t\t\t# angular source file");
#endif
puts("-ae mod\t\t\t\t# exclude modifier");
puts("-aE file\t\t\t\t# exclude modifiers from file");
puts("-ai mod\t\t\t\t# include modifier");
puts("-aI file\t\t\t\t# include modifiers from file");
#ifdef PMAP_EKSPERTZ
puts("-api xmin ymin zmin xmax ymax zmax"
"\t# rectangular region of interest"
);
puts("-apI xpos ypos zpos radius\t\t# spherical region of interest");
#endif
puts("-apg file nPhotons\t\t\t# global photon map");
puts("-apc file nPhotons\t\t\t# caustic photon map");
puts("-apd file nPhotons\t\t\t# direct photon map");
puts("-app file nPhotons bwidth\t\t# precomputed global photon map");
puts("-apv file nPhotons\t\t\t# volume photon map");
puts("-apC file nPhotons bwidth compRatio"
"\t# precomputed contribution photon map"
);
#ifdef PMAP_PHOTONFLOW
- puts("-apV file nPhotons\t\t\t# light flow volume photon map");
+ puts("-apV file nPhotons\t\t\t# light flow photon map");
+ puts("-apT file nPhotons\t\t\t# temporal light flow photon map");
#endif
printf("-apD %f\t\t\t\t# predistribution factor\n", preDistrib);
printf("-apM %d\t\t\t\t\t# max predistrib passes\n", maxPreDistrib);
#if 1
/* Kept for backwards compat, will be gradually phased out by -ld, -lr */
printf("-apm %ld\t\t\t\t# limit photon bounces\n", photonMaxBounce);
#endif
puts("-apo+ mod\t\t\t\t# photon port modifier");
puts("-apO+ file\t\t\t\t# photon ports from file");
printf("-apP %f\t\t\t\t# precomputation factor\n", finalGather);
printf("-apr %d\t\t\t\t\t# random seed\n", randSeed);
puts("-aps mod\t\t\t\t# antimatter sensor modifier");
puts("-apS file\t\t\t\t# antimatter sensors from file");
printf(backvis
? "-bv+\t\t\t\t\t# back face visibility on\n"
: "-bv-\t\t\t\t\t# back face visibility off\n"
);
printf("-dp %.1f\t\t\t\t# PDF samples / sr\n", pdfSamples);
printf("-ds %f\t\t\t\t# source partition size ratio\n", srcsizerat);
printf("-e %s\t\t\t\t# diagnostics output file\n", diagFile);
printf(clobber
? "-fo+\t\t\t\t\t# force overwrite\n"
: "-fo-\t\t\t\t\t# do not overwrite\n"
);
#ifdef PMAP_EKSPERTZ
/* (Formerly) NU STUFF for Ze Exspertz! */
printf("-ld %.1f\t\t\t\t\t# limit photon distance\n", photonMaxDist);
printf("-lr %ld\t\t\t\t# limit photon bounces\n", photonMaxBounce);
#endif
printf("-ma %.2f %.2f %.2f\t\t\t# scattering albedo\n",
colval(salbedo,RED), colval(salbedo,GRN), colval(salbedo,BLU)
);
printf("-me %.2e %.2e %.2e\t\t# extinction coefficient\n",
colval(cextinction,RED), colval(cextinction,GRN),
colval(cextinction,BLU)
);
printf("-mg %.2f\t\t\t\t# scattering eccentricity\n", seccg);
#if NIX
/* Multiprocessing on NIX only; so tuff luck, Windoze Weenies! */
printf("-n %d\t\t\t\t\t# number of parallel processes\n", nproc);
#endif
printf("-t %-9d\t\t\t\t# time between reports\n", photonRepTime);
printf(verbose
? "-v+\t\t\t\t\t# verbose console output\n"
: "-v-\t\t\t\t\t# terse console output\n"
);
#if 0
/* TODO: Do contrib/coeffs here or in rcontrib? */
printf(
"-V%c\t\t\t\t# precompute %s\n",
contrib ? '+' : '-',
contrib ? "contributions" : "coefficients"
);
#endif
}
#if 0
/* Disabled for now; will attempt to extract bins from photon map header */
static char *appendParams (char *accParams, char **params, int n)
/* Accumulate n next parameter strings from params to string accParams.
Returns new accmulated string. */
{
int i;
if (params && n)
for (i = 0; i < n; i++)
if (params [i] && strlen(params [i]))
/* Reallocate and factor in delimiter (space) and terminator */
if (accParams = realloc(
accParams, strlen(accParams) + strlen(params [i]) + 2
)) {
strcat(accParams, " ");
strcat(accParams, params [i]);
}
else
error(
SYSTEM,
"cannot append binning parameters for contrib photon map"
);
return accParams;
}
#endif
int main (int argc, char* argv [])
{
#define check(ol, al) if ( \
argv [i][ol] || badarg(argc - i - 1,argv + i + 1, al) \
) goto badopt
/* Evaluate boolean option, setting var accordingly */
#define check_bool(olen, var) switch (argv [i][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: \
goto badopt; \
}
/* Evaluate trinary option, setting bits v1 and v2 in var accordingly */
#define check_tri(olen, v1, v2, var) switch (argv [i][olen]) { \
case '\0': case '+': \
var = v1; break; \
case '-': \
var = v2; break;\
case '0': \
var = v1 | v2; break; \
default: \
goto badopt; \
}
/* Check target number of photons with optional multiplier suffix
begins with a digit */
#define checkNumPhotons(i) if (!isdigit(argv [i][0])) \
if (argv [i][0] == '-') \
goto badopt; \
else { \
sprintf(errmsg, "invalid number of photons %s", argv [i]); \
error(USER, errmsg); \
}
int loadflags = IO_CHECK | IO_SCENE | IO_TREE | IO_BOUNDS,
rval, i, j, n, binCnt = 0;
char **portLp = photonPortList, **sensLp = photonSensorList,
**amblp = NULL, sbuf [MAXSTR], portFlags [2] = "\0\0",
*binVal = NULL, *binParm = NULL;
struct stat pmstat;
/* Global program name */
progname = fixargv0(argv [0]);
/* Initialize object types */
initotypes();
/* Initialize cal function routines for contrib photon binning */
initfunc();
setcontext(RCCONTEXT);
#if defined(_WIN32) || defined(_WIN64)
/* Increase file limit to maximum */
/* XXX: DO WE NEED THIS FOR CONTRIB PHOTONS?
_setmaxstdio(2048); */
#endif
/* Parse options */
for (i = 1; i < argc; i++) {
/* Eggs-pand arguments */
while ((rval = expandarg(&argc, &argv, i)))
if (rval < 0) {
sprintf(errmsg, "cannot eggs-pand '%s'", argv [i]);
error(SYSTEM, errmsg);
}
if (argv[i] == NULL)
break;
if (!strcmp(argv [i], "-version")) {
puts(VersionID);
quit(0);
}
if (!strcmp(argv [i], "-defaults") || !strcmp(argv [i], "-help")) {
printdefaults();
quit(0);
}
/* Get octree */
if (i == argc - 1) {
octname = argv [i];
break;
}
switch (argv [i][1]) {
case 'a': /* Ambient */
switch (argv [i][2]) {
case 'i': /* Ambient include */
case 'I':
check(3, "s");
if (ambincl != 1) {
ambincl = 1;
amblp = amblist;
}
if (isupper(argv [i][2])) {
/* Add modifiers from file */
rval = wordfile(
amblp, AMBLLEN - (amblp - amblist),
getpath(argv [++i], getrlibpath(), R_OK)
);
if (rval < 0) {
sprintf(
errmsg, "cannot open ambient include file \"%s\"",
argv [i]
);
error(SYSTEM, errmsg);
}
amblp += rval;
}
else {
/* Add modifier from next arg */
*amblp++ = savqstr(argv [++i]);
*amblp = NULL;
}
break;
case 'e': /* Ambient exclude */
case 'E':
check(3, "s");
if (ambincl != 0) {
ambincl = 0;
amblp = amblist;
}
if (isupper(argv [i][2])) {
/* Add modifiers from file */
rval = wordfile(
amblp, AMBLLEN - (amblp - amblist),
getpath(argv [++i], getrlibpath(), R_OK)
);
if (rval < 0) {
sprintf(errmsg,
"cannot open ambient exclude file \"%s\"", argv [i]
);
error(SYSTEM, errmsg);
}
amblp += rval;
}
else {
/* Add modifier from next arg */
*amblp++ = savqstr(argv [++i]);
*amblp = NULL;
}
break;
case 'p': /* Pmap-specific */
switch (argv [i][3]) {
case 'g': /* Global photon map */
checkNumPhotons(i + 2);
check(4, "ss");
globalPmapParams.fileName = argv [++i];
globalPmapParams.distribTarget =
parseMultiplier(argv [++i]);
if (!globalPmapParams.distribTarget)
goto badopt;
globalPmapParams.minGather =
globalPmapParams.maxGather = 0;
break;
case 'p': /* Precomputed global photon map */
checkNumPhotons(i + 2);
check(4, "ssi");
preCompPmapParams.fileName = argv [++i];
preCompPmapParams.distribTarget =
parseMultiplier(argv [++i]);
if (!preCompPmapParams.distribTarget)
goto badopt;
preCompPmapParams.minGather =
preCompPmapParams.maxGather = atoi(argv [++i]);
if (!preCompPmapParams.maxGather)
goto badopt;
break;
case 'c': /* Caustic photon map */
checkNumPhotons(i + 2);
check(4, "ss");
causticPmapParams.fileName = argv [++i];
causticPmapParams.distribTarget =
parseMultiplier(argv [++i]);
if (!causticPmapParams.distribTarget)
goto badopt;
break;
case 'v': /* Volume photon map */
checkNumPhotons(i + 2);
check(4, "ss");
volumePmapParams.fileName = argv [++i];
volumePmapParams.distribTarget =
parseMultiplier(argv [++i]);
if (!volumePmapParams.distribTarget)
goto badopt;
break;
case 'd': /* Direct photon map */
checkNumPhotons(i + 2);
check(4, "ss");
directPmapParams.fileName = argv [++i];
directPmapParams.distribTarget =
parseMultiplier(argv [++i]);
if (!directPmapParams.distribTarget)
goto badopt;
break;
case 'C': /* Precomputed contribution photon map */
checkNumPhotons(i + 2);
check(4, "ssif");
contribPmapParams.fileName = argv [++i];
contribPmapParams.distribTarget =
parseMultiplier(argv [++i]);
if (!contribPmapParams.distribTarget)
goto badopt;
contribPmapParams.minGather =
contribPmapParams.maxGather = atoi(argv [++i]);
if (!contribPmapParams.maxGather)
goto badopt;
compressRatio = atof(argv [++i]);
if (compressRatio < FTINY ||
compressRatio > 1 - FTINY
) {
error(USER, "contribution photon "
"compression ratio must be in range (0, 1)"
);
goto badopt;
}
break;
#ifdef PMAP_PHOTONFLOW
- case 'V': /* Light flow volume photon map */
+ /* Light flow is a variant of volume photon map */
+ case 'V': /* Light flow photon map */
+ case 'T': /* Temporal light flow photon map */
check(4, "ss");
lightFlowPmapParams.fileName = argv [++i];
lightFlowPmapParams.distribTarget =
parseMultiplier(argv [++i]);
if (!lightFlowPmapParams.distribTarget)
goto badopt;
/* Set zero absorption and forward scattering for
global mist; extinction up to user of local
mist defined in octree */
/* setcolor(cextinction, 1, 1, 1); */
setcolor(salbedo, 1, 1, 1);
seccg = 1;
break;
#endif
case 'D': /* Predistribution factor */
check(4, "f");
preDistrib = atof(argv [++i]);
if (preDistrib <= 0)
error(USER, "predistrib factor must be > 0");
break;
case 'M': /* Max predistribution passes */
check(4, "i");
maxPreDistrib = atoi(argv [++i]);
if (maxPreDistrib <= 0)
error(USER, "max predistrib passes must be > 0");
break;
#if 1
/* Kept for backwards compat, to be phased out by -lr */
case 'm': /* Max photon bounces */
check(4, "i");
photonMaxBounce = atol(argv [++i]);
if (photonMaxBounce <= 0)
error(USER, "max photon bounces must be > 0");
break;
#endif
#ifdef PMAP_EKSPERTZ
case 'i': /* Add rectangular region of interest */
case 'I': /* Add spherical region of interest */
check(4, isupper(argv [j=i][3]) ? "ffff" : "ffffff");
n = pmapNumROI;
pmapROI = realloc(pmapROI,
++pmapNumROI * sizeof(PhotonMapROI)
);
if (!pmapROI)
error(SYSTEM, "failed to allocate ROI");
pmapROI [n].pos [0] = atof(argv [++i]);
pmapROI [n].pos [1] = atof(argv [++i]);
pmapROI [n].pos [2] = atof(argv [++i]);
pmapROI [n].siz [0] = atof(argv [++i]);
if (isupper(argv [j][3])) {
/* Spherical ROI; radius^2 */
pmapROI [n].siz [0] *= pmapROI [n].siz [0];
PMAP_ROI_SETSPHERE(pmapROI + n);
if (pmapROI [n].siz [0] <= FTINY)
error(USER,
"region of interest has invalid radius"
);
}
else {
/* Rectangular ROI */
pmapROI [n].siz [1] = atof(argv [++i]);
pmapROI [n].siz [2] = atof(argv [++i]);
for (j = 0; j < 3; j++) {
/* Pos at rectangle centre, siz symmetric */
pmapROI [n].pos [j] = 0.5 * (
pmapROI [n].pos [j] + pmapROI [n].siz [j]
);
pmapROI [n].siz [j] = fabs(
pmapROI [n].siz [j] - pmapROI [n].pos [j]
);
if (pmapROI [n].siz [j] <= FTINY)
error(USER,
"region of interest has invalid size"
);
}
}
break;
#endif
case 'P': /* Global photon precomp ratio */
check(4, "f");
finalGather = atof(argv [++i]);
if (finalGather <= 0 || finalGather > 1)
error(USER, "global photon precomputation ratio "
"must be in range ]0, 1]"
);
break;
case 'o': /* Photon port */
case 'O':
/* Check for bad arg and length, taking into account
* default forward orientation if none specified, in
* order to maintain previous behaviour */
check(argv [i][4] ? 5 : 4, "s");
/* Get port orientation flags */
check_tri(4, PMAP_PORTFWD, PMAP_PORTBWD,
portFlags [0]
);
if (isupper(argv [i][3])) {
/* Add port modifiers from file */
rval = wordfile(portLp,
MAXSET - (portLp - photonPortList),
getpath(argv [++i], getrlibpath(), R_OK)
);
if (rval < 0) {
sprintf(errmsg,
"cannot open photon port file %s", argv [i]
);
error(SYSTEM, errmsg);
}
/* HACK: append port orientation flags to every
* modifier; note this requires reallocation */
for (; rval--; portLp++) {
j = strlen(*portLp);
if (!(*portLp = realloc(*portLp, j + 2))) {
sprintf(errmsg,
"cannot allocate photon port modifiers"
" from file %s", argv [i]
);
error(SYSTEM, errmsg);
}
strcat(*portLp, portFlags);
}
}
else {
/* Append port flags to port modifier, add to
* port list and mark of end list with NULL */
strcpy(sbuf, argv [++i]);
strcat(sbuf, portFlags);
*portLp++ = savqstr(sbuf);
*portLp = NULL;
}
break;
case 'r': /* Random seed */
check(4, "i");
randSeed = atoi(argv [++i]);
break;
case 's': /* Antimatter sensor */
case 'S':
check(4, "s");
if (isupper(argv[i][3])) {
/* Add sensor modifiers from file */
rval = wordfile(sensLp,
MAXSET - (sensLp - photonSensorList),
getpath(argv [++i], getrlibpath(), R_OK)
);
if (rval < 0) {
sprintf(errmsg,
"cannot open antimatter sensor file %s",
argv [i]
);
error(SYSTEM, errmsg);
}
sensLp += rval;
}
else {
/* Append modifier to sensor list, mark end with
* NULL */
*sensLp++ = savqstr(argv [++i]);
*sensLp = NULL;
}
break;
default:
goto badopt;
}
break;
default:
goto badopt;
}
break;
case 'b':
switch (argv [i][2]) {
case 'v':
/* Back face visibility */
check_bool(3, backvis);
break;
case 'n':
/* Binning count for precomputed contrib. photon map.
Save params for export to file for rcontrib. */
check(3, "s");
/* binParams = appendParams(binParams, argv + i, 2); */
binCnt = (int)(eval(argv [++i]) + .5);
break;
default:
/* Binning expression for precomp. contrib. photon map.
Save params for export to file for rcontrib. */
check(2,"s");
/* binParams = appendParams(binParams, argv + i, 2); */
binVal = argv [++i];
}
break;
case 'd': /* Direct */
switch (argv [i][2]) {
case 'p': /* PDF samples */
check(3, "f");
pdfSamples = atof(argv [++i]);
break;
case 's': /* Source partition size ratio */
check(3, "f");
srcsizerat = atof(argv [++i]);
break;
default:
goto badopt;
}
break;
case 'e':
switch (argv [i][2]) {
case 'f':
/* Diagnostics file */
check(2, "s");
diagFile = argv [++i];
break;
default:
/* Functional expression for precomputed contrib. pmap */
check(2, "s");
/* binParams = appendParams(binParams, argv + i, 2); */
scompile(argv [++i], NULL, 0);
}
break;
case 'f':
switch (argv [i][2]) {
case 'o':
/* Force overwrite */
check_bool(3, clobber);
break;
default:
/* Function file for precomputed contribution photon map */
check(2, "s");
/* binParams = appendParams(binParams, argv + i, 2); */
loadfunc(argv[++i]);
}
break;
#ifdef PMAP_EKSPERTZ
case 'l': /* Limits */
switch (argv [i][2]) {
case 'd': /* Limit photon path distance */
check(3, "f");
photonMaxDist = atof(argv [++i]);
if (photonMaxDist <= 0)
error(USER, "max photon distance must be > 0");
break;
case 'r': /* Limit photon bounces */
check(3, "i");
photonMaxBounce = atol(argv [++i]);
if (photonMaxBounce <= 0)
error(USER, "max photon bounces must be > 0");
break;
default:
goto badopt;
}
break;
#endif
case 'm':
switch (argv[i][2]) {
case 'e': /* Mist Eggs-tinction coefficient */
check(3, "fff");
setcolor(cextinction, atof(argv [i + 1]),
atof(argv [i + 2]), atof(argv [i + 3])
);
i += 3;
break;
case 'a':
/* Mist albedo */
check(3, "fff");
setcolor(salbedo, atof(argv [i + 1]),
atof(argv [i + 2]), atof(argv [i + 3])
);
i += 3;
break;
case 'g':
/* Mist scattering eccentricity */
check(3, "f");
seccg = atof(argv [++i]);
break;
default:
/* Modifier name for precomputed contrib. photon map */
check(2, "s");
/* binParams = appendParams(binParams, argv + i, 2); */
addSourceModifier(&srcContrib, &numSrcContrib,
argv [++i], binParm, binVal, binCnt
);
break;
}
break;
case 'M':
/* Modifier file for precomputed contribution photon map */
check(2, "s");
/* binParams = appendParams(binParams, argv + i, 2); */
addSourceModfile(&srcContrib, &numSrcContrib,
argv [++i], binParm, binVal, binCnt
);
break;
#if NIX
case 'n': /* Num parallel processes (NIX only) */
check(2, "i");
nproc = atoi(argv [++i]);
if (nproc > PMAP_MAXPROC) {
nproc = PMAP_MAXPROC;
sprintf(errmsg, "too many parallel processes, "
"clamping to %d\n", nproc
);
error(WARNING, errmsg);
}
break;
#endif
case 'p':
/* Parameter setting(s) for precomputed contrib. photon map */
check(2, "s");
/* binParams = appendParams(binParams, argv + i, 2); */
set_eparams(binParm = argv [++i]);
break;
case 't': /* Timer */
check(2, "i");
photonRepTime = atoi(argv [++i]);
break;
case 'v': /* Verbosity */
check_bool(2, verbose);
break;
#ifdef EVALDRC_HACK
case 'A': /* Capt. B's (now historic) angular source file hack */
check(2,"s");
angsrcfile = argv[++i];
break;
#endif
default:
goto badopt;
}
}
/* Open diagnostics file */
if (diagFile) {
if (!freopen(diagFile, "a", stderr))
quit(2);
fprintf(stderr, "**************\n*** PID %5d: ", getpid());
printargs(argc, argv, stderr);
putc('\n', stderr);
fflush(stderr);
}
#ifdef NICE
/* Lower priority */
nice(NICE);
#endif
if (octname == NULL)
error(USER, "missing octree argument");
/* Allocate photon maps and set parameters */
for (i = 0; i < NUM_PMAP_TYPES; i++) {
setPmapParam(photonMaps + i, pmapParams + i);
/* Don't overwrite existing photon map unless clobbering enabled */
if (photonMaps [i] && !stat(photonMaps [i] -> fileName, &pmstat) &&
!clobber
) {
sprintf(errmsg, "photon map file %s exists, not overwritten",
photonMaps [i] -> fileName
);
error(USER, errmsg);
}
}
for (i = 0; i < NUM_PMAP_TYPES && !photonMaps [i]; i++);
if (i >= NUM_PMAP_TYPES)
error(USER, "no photon maps specified");
readoct(octname, loadflags, &thescene, NULL);
#ifdef EVALDRC_HACK
if (angsrcfile)
readobj(angsrcfile); /* load angular sources */
#endif
nsceneobjs = nobjects;
/* Get sources */
marksources();
/* Do forward pass and build photon maps */
if (contribPmap)
/* Just build contrib pmap, ignore others */
distribPhotonContrib(contribPmap, &srcContrib, numSrcContrib, nproc);
else
distribPhotons(photonMaps, nproc);
/* TODO: Where do we save binParams and wavelet coeffs? */
/* Save photon maps; no idea why GCC needs an explicit cast here... */
savePmaps((const PhotonMap**)photonMaps, argc, argv);
cleanUpPmaps(photonMaps);
quit(0);
badopt:
sprintf(errmsg, "command line error at '%s'", argv[i]);
error(USER, errmsg);
#undef check
#undef check_bool
return 0;
}
diff --git a/pmapdata.c b/pmapdata.c
index 91df123..1687d15 100644
--- a/pmapdata.c
+++ b/pmapdata.c
@@ -1,768 +1,774 @@
#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,
#ifdef PMAP_PHOTONFLOW
NULL
#endif
};
/* Include routines to handle underlying point cloud data structure */
#ifdef PMAP_OOC
#include "pmapooc.c"
#else
#include "pmapkdt.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_PHOTONFLOW
/* Init path lookup table and its key buffer */
pmap -> pathLUT = NULL;
pmap -> pathLUTKeys = NULL;
pmap -> numPathLUTKeys = 0;
#endif
/* 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 origin as unused */
pmap -> lastContribOrg.srcIdx = pmap -> lastContribOrg.srcBin = -1;
pmap -> numContribOrg = 0;
pmap -> contribOrg = 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");
#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);
if (isContribPmap(pmap)) {
/* Set contrib photon's origin and subprocess index (the latter to
* linearise the origin array after photon distribution is complete).
* Also set origin's source index, thereby marking it as used. Note
* the origin's bin has already been set by newPhotonOrigin(). */
photon.aux.contribOrgIdx = pmap -> numContribOrg;
photon.proc = PMAP_GETRAYPROC(ray);
pmap -> lastContribOrg.srcIdx = ray -> rsrc;
}
else {
/* Set photon's path ID from ray serial num and subprocess index */
photon.aux.pathID = ray -> rno;
photon.proc = PMAP_GETRAYPROC(ray);
}
/* Set normal */
for (i = 0; i <= 2; i++)
photon.norm [i] = 127.0 * (isVolumePmap(pmap)
? 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,
PhotonContribOriginIdx *originOfs, 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;
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 (originOfs)
/* Linearise contrib photon origin index from subprocess index
* using the per-subprocess offsets in originOfs */
p -> aux.contribOrgIdx += originOfs [p -> proc];
/* 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"
);
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;
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);
}
}
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
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;
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;
/* 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_PHOTONFLOW
- if (isVolumePmap(pmap) && !isLightFlowPmap(pmap)) {
+ /* If sphericalIrrad = 0, planar irradiance is evaluated from
+ lightflow photons, hence these are associated with the plane
+ defined by the ray normal, and filtered accordingly by their
+ incident directions.
+ If spherical Irrad = 1, mean spherical irradiance is evaluated,
+ the normal is passed as NULL, which disables normal filtering. */
+ if (isVolumePmap(pmap) && sphericalIrrad) {
#else
if (isVolumePmap(pmap)) {
#endif
/* Search position is ray -> rorg for volume photons, since we have
no intersection point. Photon normals are ignored; these are
incident directions for the volume scattering func). */
#ifdef PMAP_OOC
OOC_FindPhotons(pmap, ray -> rorg, NULL);
#else
kdT_FindPhotons(pmap, ray -> rorg, NULL);
#endif
}
else {
#ifdef PMAP_OOC
OOC_FindPhotons(pmap, ray -> rop, ray -> ron);
#else
kdT_FindPhotons(pmap, ray -> rop, ray -> ron);
#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 : "<null>"
);
#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 : "<null>"
);
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_PHOTONFLOW
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/pmapdata.h b/pmapdata.h
index f0e36af..d15e83e 100644
--- a/pmapdata.h
+++ b/pmapdata.h
@@ -1,392 +1,392 @@
/* RCSid $Id: pmapdata.h,v 2.14 2020/04/08 15:14:21 rschregle Exp $ */
/*
=========================================================================
Photon map types and 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}).
Defining PMAP_FLOAT_FLUX stores photon flux as floats rather than packed
RGBE for greater precision; this may be necessary when the flux differs
significantly in individual colour channels, e.g. with highly saturated
colours.
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.h,v 2.14 2020/04/08 15:14:21 rschregle Exp $
*/
#ifndef PMAPDATA_H
#define PMAPDATA_H
#ifndef NIX
#if defined(_WIN32) || defined(_WIN64)
#define NIX 0
#else
#define NIX 1
#endif
#endif
#if (defined(PMAP_OOC) && !NIX)
#error "OOC currently only supported on NIX -- tuff luck."
#endif
#include "ray.h"
#include "pmaptype.h"
#include "paths.h"
#include "lookup.h"
#include <stdint.h>
/* Photon origin for light source contributions; these are only used
to precompute contribution photons. They are referenced by contribution
photons (see origin field in struct Photon below) in a surjective
mapping, since multiple photons may share the same origin if they
lie along its associated path. For this reason it is more efficient
to factor this data out of the photons themselves and consolidate
it here until the photons have been precomputed, after which it is
no longer needed. */
typedef struct {
int16 srcIdx, /* Index of emitting light source */
srcBin; /* Bin number corresponding to incident dir */
} PhotonContribOrigin;
typedef uint32 PhotonPathID;
typedef uint32 PhotonContribOriginIdx;
#define PMAP_MAXCONTRIBORG UINT32_MAX
#define photonSrcIdx(pm, p) ((pm) -> contribOrg \
? (pm) -> contribOrg [(p) -> aux.contribOrgIdx].srcIdx \
: (p) -> aux.pathID\
)
#define photonSrcBin(pm, p) ( \
(pm) -> contribOrg [(p) -> aux.contribOrgIdx].srcBin \
)
#define photonSrcMod(pm, p) findmaterial(source [photonSrcIdx(pm, p)].so)
- /* Multipurpose auxiliary photon attribute type and limit */
+ /* Multipurpose auxiliary photon attribute type */
typedef union {
/* Distance travelled by photon (= path length / time of flight)
for temporal light flow */
float pathLen;
/* Index into PhotonContribOrigin array for contrib photons */
PhotonContribOriginIdx contribOrgIdx;
/* Unique path ID for all other photon types */
PhotonPathID pathID;
} PhotonAuxAttrib;
/* Macros for photon's generating subprocess field */
#ifdef PMAP_OOC
#define PMAP_PROCBITS 7
#else
#define PMAP_PROCBITS 5
#endif
#define PMAP_MAXPROC (1 << PMAP_PROCBITS)
#define PMAP_GETRAYPROC(r) ((r) -> crtype >> 8)
#define PMAP_SETRAYPROC(r,p) ((r) -> crtype |= p << 8)
typedef struct {
float pos [3]; /* Photon position */
signed char norm [3]; /* Encoded normal / incident
direction [volume photons] */
union {
struct {
#ifndef PMAP_OOC
unsigned char discr : 2; /* kd-tree discriminator axis */
#endif
unsigned char caustic : 1; /* Specularly scattered (=caustic) */
/* Photon's generating subprocess index, used for primary ray
* index linearisation when building contrib pmaps; note this is
* reduced for kd-tree to accommodate the discriminator field */
unsigned char proc : PMAP_PROCBITS;
};
unsigned char flags;
};
/* Photon flux in watts or lumen /
photon contribution [contrib photons] /
average wavelet coefficient [precomputed contrib photons] */
#ifdef PMAP_FLOAT_FLUX
COLOR flux;
#else
COLR flux;
#endif
/* Auxiliary field; this is a multipurpose, type-specific field used
by the following photon types (as identified by enum PhotonMapType
in pmaptype.h):
PMAP_TYPE_CONTRIB: Index into photon map's contrib
origin array.
PMAP_TYPE_TEMPLIGHTFLOW: Distance travelled by photon / time
of flight
All others: Photon path ID.
*/
PhotonAuxAttrib aux;
} Photon;
/* Define PMAP_FLOAT_FLUX to store photon flux as floats instead of
* compact RGBE, which was found to improve accuracy in analytical
* validation. */
#ifdef PMAP_FLOAT_FLUX
#define setPhotonFlux(p,f) copycolor((p) -> flux, f)
#define getPhotonFlux(p,f) copycolor(f, (p) -> flux)
#else
#define setPhotonFlux(p,f) setcolr((p)->flux, (f)[0], (f)[1], (f)[2])
#define getPhotonFlux(p,f) colr_color(f, (p) -> flux)
#endif
/* Bias compensation history node */
typedef struct {
COLOR irrad;
float weight;
} PhotonBiasCompNode;
/* Forward declaration */
struct PhotonMap;
/* Define search queue and underlying data struct types */
#ifdef PMAP_OOC
#include "pmapooc.h"
#else
#include "pmapkdt.h"
#endif
/* Mean size of heapfile write buffer, in number of photons */
#define PMAP_HEAPBUFSIZE 1e6
/* Mean idle time between heap locking attempts, in usec */
#define PMAP_HEAPBUFSLEEP 2e6
/* Temporary filename for photon heaps */
#define PMAP_TMPFNAME TEMPLATE
#define PMAP_TMPFNLEN (TEMPLEN + 1)
typedef struct PhotonMap {
PhotonMapType type; /* See pmaptype.h */
char *fileName; /* Photon map file */
/* ================================================================
* PRE/POST-BUILD STORAGE
* ================================================================ */
FILE *heap; /* Unsorted photon heap prior to
construction of store */
char heapFname [sizeof(PMAP_TMPFNAME)];
Photon *heapBuf; /* Write buffer for above */
unsigned long heapBufLen, /* Current & max size of heapBuf */
heapBufSize;
PhotonStorage store; /* Photon storage in space
subdividing data struct */
/* ================================================================
* PHOTON DISTRIBUTION STUFF
* ================================================================ */
unsigned long distribTarget, /* Num stored specified by user */
numPhotons; /* Num actually stored */
float distribRatio; /* Probability of photon storage */
COLOR photonFlux; /* Average photon flux */
unsigned short randState [3]; /* Local RNG state */
/* ================================================================
* PHOTON LOOKUP STUFF
* ================================================================ */
union { /* Flags passed to findPhotons() */
char lookupCaustic : 1;
char lookupFlags;
};
PhotonSearchQueue squeue; /* Search queue for photon lookups */
unsigned minGather, /* Specified min/max photons per */
maxGather; /* density estimate */
/* NOTE: All distances are SQUARED */
float maxDist2, /* Max search radius */
maxDist0, /* Initial value for above */
maxDist2Limit, /* Hard limit for above */
gatherTolerance; /* Fractional deviation from minGather/
maxGather for short lookup */
void (*lookup)(
struct PhotonMap*, RAY*, COLOR
); /* Callback for type-specific photon
* lookup (usually density estimate) */
#ifdef PMAP_PHOTONFLOW
LUTAB *pathLUT; /* Photon path lookup table to filter
volume photons */
char **pathLUTKeys; /* Preallocated buffer to store keys
for path lookup table */
unsigned numPathLUTKeys; /* Num keys currently in key buffer
(= next free entry at tail) */
#endif
/* ================================================================
* BIAS COMPENSATION STUFF
* ================================================================ */
PhotonBiasCompNode *biasCompHist; /* Bias compensation history */
/* ================================================================
* STATISTIX
* ================================================================ */
unsigned long totalGathered, /* Total photons gathered */
numDensity, /* Num density estimates */
numLookups, /* Counters for short photon lookups */
numShortLookups;
unsigned minGathered, /* Min/max photons actually gathered */
maxGathered, /* per density estimate */
shortLookupPct; /* % of short lookups for stats */
float minError, /* Min/max/rms density estimate error */
maxError,
rmsError,
CoGdist, /* Avg distance to centre of gravity */
maxPos [3], /* Max & min photon positions */
minPos [3];
FVECT CoG; /* Centre of gravity (avg photon pos) */
/* ================================================================
* CONTRIBUTION PHOTON STUFF
* ================================================================ */
PhotonContribOrigin *contribOrg, /* Contribution origin array */
lastContribOrg; /* Current contrib origin for
contrib photon distribution */
PhotonContribOriginIdx numContribOrg; /* Number of contrib origins */
LUTAB *srcContrib; /* LUT for source contribs */
} PhotonMap;
/* Photon maps by type (see PhotonMapType) */
extern PhotonMap *photonMaps [];
/* Macros for specific photon map types */
#define globalPmap (photonMaps [PMAP_TYPE_GLOBAL])
#define preCompPmap (photonMaps [PMAP_TYPE_PRECOMP])
#define causticPmap (photonMaps [PMAP_TYPE_CAUSTIC])
#define directPmap (photonMaps [PMAP_TYPE_DIRECT])
#define contribPmap (photonMaps [PMAP_TYPE_CONTRIB])
#ifdef PMAP_PHOTONFLOW
#define volumePmap ( \
lightFlowPmap ? lightFlowPmap : photonMaps [PMAP_TYPE_VOLUME] \
)
#define lightFlowPmap (photonMaps [PMAP_TYPE_LIGHTFLOW])
#else
#define volumePmap (photonMaps [PMAP_TYPE_VOLUME])
#endif
/* Photon map type tests */
#define isGlobalPmap(p) ((p) -> type == PMAP_TYPE_GLOBAL)
#define isCausticPmap(p) ((p) -> type == PMAP_TYPE_CAUSTIC)
#define isContribPmap(p) ((p) -> type == PMAP_TYPE_CONTRIB)
#ifdef PMAP_PHOTONFLOW
#define isVolumePmap(p) ( \
(p) -> type == PMAP_TYPE_VOLUME || isLightFlowPmap(p) \
)
#define isLightFlowPmap(p) ((p) -> type == PMAP_TYPE_LIGHTFLOW)
#else
#define isVolumePmap(p) ((p) -> type == PMAP_TYPE_VOLUME)
#endif
void initPhotonMap (PhotonMap *pmap, PhotonMapType t);
/* Initialise empty photon map of specified type */
int newPhoton (PhotonMap *pmap, const RAY *ray);
/* Create new photon with ray's direction, intersection point, and flux,
and append to unsorted photon heap pmap -> heap. The photon is
accepted with probability pmap -> distribRatio for global density
control; if the photon is rejected, -1 is returned, else 0. The flux
is scaled by ray -> rweight and 1 / pmap -> distribRatio. */
void initPhotonHeap (PhotonMap *pmap);
/* Open photon heap file */
void flushPhotonHeap (PhotonMap *pmap);
/* Flush photon heap buffa pmap -> heapBuf to heap file pmap -> heap;
* used by newPhoton() and to finalise heap in distribPhotons(). */
void buildPhotonMap (
PhotonMap *pmap, double *photonFlux,
PhotonContribOriginIdx *originOfs, unsigned nproc
);
/* Postpress unsorted photon heap pmap -> heap and build underlying data
* structure pmap -> store. This is prerequisite to photon lookups with
* findPhotons(). */
/* PhotonFlux is the flux per photon averaged over RGB; this is
* multiplied with each photon's flux during the postprocess. In the
* case of a contribution photon map, this is an array with a separate
* flux specific to each light source due to non-uniform photon emission;
* Otherwise it is referenced as a scalar value. Flux is not scaled if
* photonFlux == NULL. */
/* Photon map construction may be parallelised if nproc > 1, if
* supported. The heap is destroyed on return. */
/* OriginOfs is an array of index offsets for the contribution photon
* origins in pmap->contribOrg generated by each of the nproc subprocesses
* during contrib photon distribution (see distribPhotonContrib()). These
* offsets are used to linearise the photon origin indices in the
* postprocess. This linearisation is skipped if originOfs == NULL,
* e.g. when building a global/caustic/volume photon map, where the
* origins are serial path IDs. */
void findPhotons (PhotonMap* pmap, const RAY *ray);
/* Find pmap -> squeue.len closest photons to ray -> rop with similar
normal. For volume photons ray -> rorg is used and the normal is
ignored (being the incident direction in this case). Found photons
are placed search queue starting with the furthest photon at pmap ->
squeue.node, and pmap -> squeue.tail being the number actually found. */
Photon *find1Photon (PhotonMap *pmap, const RAY *ray, Photon *photon);
/* Find single closest photon to ray -> rop with similar normal.
Return NULL if none found, else the supplied Photon* buffer,
indicating that it contains a valid photon. */
void getPhoton (PhotonMap *pmap, PhotonIdx idx, Photon *photon);
/* Retrieve photon referenced by idx from pmap -> store */
Photon *getNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx);
/* Retrieve photon from NN search queue after calling findPhotons() */
PhotonIdx firstPhoton (const PhotonMap *pmap);
/* Index to first photon, to be passed to getPhoton(). Indices to
* subsequent photons can be optained via increment operator (++) */
void deletePhotons (PhotonMap*);
/* Free dem mammaries... */
#endif
diff --git a/pmapdens.c b/pmapdens.c
index 46a0278..2053763 100644
--- a/pmapdens.c
+++ b/pmapdens.c
@@ -1,574 +1,586 @@
#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,
#ifdef PMAP_PHOTONFLOW
lightFlowPhotonDensity
#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 : "<null>",
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 photon dist */
scalecolor(flux, 2 * (1 - sqn -> dist2 / r2));
#endif
addcolor(irrad, flux);
}
/* Divide by search area PI * r^2, 1 / PI required as ambient
normalisation factor */
scalecolor(irrad, 1 / (PI * PI * r2));
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 1 / (4 * PI) */
scalecolor(irrad, 3 / (16 * PI * PI * r2 * sqrt(r2)));
}
#ifdef PMAP_PHOTONFLOW
void volumePhotonSphIrrad (PhotonMap *pmap, RAY *ray, COLOR irrad)
/* Evaluate mean spherical irradiance from volume photons inscattered at
* ray->rop, ignoring participating medium. This evaluates the scalar
* irradiance of a physical light field represented by "photon flow". */
{
unsigned i;
float r2;
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;
/* Issue warning if few too photons; this will particularly happen
if PMAP_PHOTONFILT is enabled, since photons from duplicate paths
are pruned during the search, and this is not accounted for by
the automatic search radius adjustment. As a workaround for now,
it's up to the user provide a sufficiently large fixed search
radius with the -ma option to rtrace */
if (pmap -> squeue.tail < pmap -> squeue.len) {
- sprintf(
- errmsg,
+ sprintf(errmsg,
"short lookup in photon flow; consider -am %.3f or greater",
2 * sqrt(pmap -> maxDist2)
);
error(WARNING, errmsg);
}
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 photon dist */
/* XXX: THIS KERNEL IS UNTESTED!!! */
scalecolor(flux, 1 - sqn -> dist2 / r2);
#endif
addcolor(irrad, flux);
#ifdef PMAP_DEBUGPATHS
- printf(
- "%f\t%f\t%f\t%f\t%f\t%d:%010d\n",
- sqrt(sqn -> dist2),
+ 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],
flux [0], photon -> proc, photon -> primary
);
#endif
}
#ifdef PMAP_DEBUGPATHS
printf("N = %d\tr = %f\n", pmap -> squeue.tail, sqrt(r2));
#endif
#ifdef PMAP_PATHFILT
/* Divide accumulated flux by sphere surface area and ambient
normalisation factor PI */
scalecolor(irrad, 1 / (4 * PI * PI * r2));
#else
/* Divide accumulated flux by sphere volume and ambient normalisation
factor PI */
scalecolor(irrad, 3 / (4 * PI * PI * r2 * sqrt(r2)));
#endif
}
void volumePhotonCubicIllum (PhotonMap *pmap, RAY *ray, COLOR illum)
/* Evaluate cubic illuminance according to Cuttle from volume photons
* inscattered at ray->rop, ignoring participating medium. The evaluation
* interprets the volume photon map as a representation of the physical
* light field ("photon flow"). */
{
unsigned ax, p, i ,j;
float r2, cosNorm, volNorm;
Photon *photon;
const PhotonSearchQueueNode *sqn;
COLOR photonFlux, E [6],
Evec [3], Esymm [3],
EvecLum, EsymmSum = {0, 0, 0};
FVECT norm0;
const FVECT norm [6] = {
/* Cube face normals (INVERTED since volume photon directions
point TOWARDS their halfspace of incidence */
{-1, 0, 0}, {1, 0, 0}, {0, -1, 0}, {0, 1, 0}, {0, 0, -1}, {0, 0, 1}
};
setcolor(illum, 0, 0, 0);
if (!pmap -> maxGather)
return;
/* Save ray normal */
VCOPY(norm0, ray -> ron);
#if 1
/* Perform separate lookup for each axis. This improves accuracy
particularly for low bandwitdhs, as the radius adjusts to the
density of incident photons along each axis. This comes at the
expense of 6 lookups.*/
for (ax = 0; ax < 6; ax++) {
/* Set cube face normal for current axis */
VCOPY(ray -> ron, norm [ax]);
setcolor(E [ax], 0, 0, 0);
/* Find photons incident from this axis, filtering by normal */
findPhotons(pmap, ray);
/* Need at least 2 photons, else skip */
if (pmap -> squeue.tail < 2)
continue;
if (pmap -> squeue.tail < pmap -> squeue.len) {
sprintf(
errmsg,
"short lookup in photon flow; consider -am %.3f or greater",
2 * sqrt(pmap -> maxDist2)
);
error(WARNING, errmsg);
}
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 its dist */
/* XXX: THIS KERNEL IS UNTESTED!!! */
scalecolor(photonFlux, 1 - sqn -> dist2 / r2);
#endif
/* Project photon dir onto axis */
cosNorm = DOT(norm [ax], photon -> norm) / 127;
scalecolor(photonFlux, cosNorm);
addcolor(E [ax], 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
}
/* Normalise illuminance for spherical volume */
scalecolor(E [ax], 3 / (4 * PI * r2 * sqrt(r2)));
printf(
"E [%s%d] = %.1f\n",
ax & 1 ? "-" : "+", ax >> 1, luminance(E [ax])
);
}
#else
/* Perform 1 lookup for all axes, and distribute photon flux
according incident direction. While faster than separate lookups
per axis, the accuracy is poor along those axes with few incident
photons, particularly in conjunction with low bandwidths. */
findPhotons(pmap, ray);
/* Need at least 2 photons */
if (pmap -> squeue.tail < 2)
return;
if (pmap -> squeue.tail < pmap -> squeue.len) {
sprintf(
errmsg,
"short lookup in photon flow; consider -am %.3f or greater",
2 * sqrt(pmap -> maxDist2)
);
error(WARNING, errmsg);
}
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));
/* Volumetric normalisation factor */
#ifdef PMAP_PATHFILT
volNorm = 1 / (2 * PI * r2);
#else
# if 1
volNorm = 3 / (4 * PI * r2 * sqrt(r2));
#else
/* Hemispherical volume */
volNorm = 3 / (2 * PI * r2 * sqrt(r2));
#endif
#endif
#ifdef PMAP_DEBUGPATHS
printf("N = %d\tr = %f\n", pmap -> squeue.tail, sqrt(r2));
#endif
for (ax = 0; ax < 6; ax++) {
/* Set cube face normal for current axis */
setcolor(E [ax], 0, 0, 0);
}
/* 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 its dist */
/* XXX: THIS KERNEL IS UNTESTED!!! */
scalecolor(photonFlux, 1 - sqn -> dist2 / r2);
#endif
for (ax = 0; ax < 6; ax++) {
cosNorm = DOT(norm [ax], photon -> norm) / 127;
if (cosNorm > 0) {
/* Project photon dir onto axis */
scalecolor(photonFlux, cosNorm);
addcolor(E [ax], 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
}
for (ax = 0; ax < 6; ax++) {
/* Normalise illuminance for spherical volume */
scalecolor(E [ax], volNorm);
printf(
"E [%s%d] = %.1f\n",
ax & 1 ? "-" : "+", ax >> 1, luminance(E [ax])
);
}
#endif
/* Calculate illuminance vector and symmetric illuminance. Unlike
Cuttle's original formulation, we defer scalar reduction and
preserve RGB components in chromatic analysis is required. */
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
Evec [i][j] = E [i << 1][j] - E [(i << 1) + 1][j];
Esymm [i][j] = 0.5 * (
E [i << 1][j] + E [(i << 1) + 1][j] - fabs(Evec [i][j])
);
}
/* Reduce Evec to luminance per axis, sum Esymm over axes */
EvecLum [i] = luminance(Evec [i]);
addcolor(EsymmSum, Esymm [i]);
}
scalecolor(EsymmSum, 1.0 / 3);
printf(
"Evec\t= [%.1f\t%.1f\t%.1f]\n",
EvecLum [0], EvecLum [1], EvecLum [2]
);
printf(
"Esym\t= [%.1f\t%.1f\t%.1f]\n",
luminance(Esymm [0]), luminance(Esymm [1]), luminance(Esymm [2])
);
/* Reduce to scalars and return as triplet:
- illuminance vector magnitude |Evec|
- scalar symmetry |Esymm|
- scalar illuminance |Esr|. */
illum [0] = sqrt(DOT(EvecLum, EvecLum));
illum [1] = luminance(EsymmSum);
illum [2] = 0.25 * illum [0] + illum [1];
/* Normalise by RADIANCE-specific irradiance factor PI */
scalecolor(illum, 1 / PI);
/* Restore ray normal */
VCOPY(ray -> ron, norm0);
}
#endif /* PMAP_PHOTONFLOW */
#ifdef PMAP_PHOTONFLOW
void lightFlowPhotonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad)
- /* Evaluate irradiance from volume photons at ray->rop for plane defined
- * by normal ray -> ron, ignoring participating medium. The evaluation
- * interprets the volume photon map as a representation of the "flow"
- * of light, i.e. a physical light field. */
+ /* Evaluate irradiance from volume photons at position 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. */
{
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 */
+ incidence (has no effect if evaluating mean spherical irrad) */
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) {
sprintf(errmsg,
"short lookup in photon flow; consider -am %.3f or greater",
2 * sqrt(pmap -> maxDist2)
);
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 its dist */
/* XXX: THIS KERNEL IS UNTESTED!!! */
scalecolor(photonFlux, 1 - sqn -> dist2 / r2);
#endif
- /* Jitter photon direction since stored as 8-bit signed int */
- for (i = 0; i < 3; i++) {
- photonDir [i] = photon -> norm [i];
+ 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();
+ }
- 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);
}
-
- /* 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
}
- /* Normalise accumulated flux by spherical volume, include
- RADIANCE-specific irradiance factor PI */
+ /* Normalise accumulated flux by spherical volume to obtain irradiance,
+ include RADIANCE-specific irradiance factor PI */
scalecolor(irrad, 3 / (4 * PI * PI * r2 * sqrt(r2)));
+ if (sphericalIrrad)
+ /* Scale irradiance by 1/4; this corresponds to the "true" scalar
+ illuminance as defined in equation (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/pmapopt.c b/pmapopt.c
index dbde40c..5623f9a 100644
--- a/pmapopt.c
+++ b/pmapopt.c
@@ -1,142 +1,159 @@
#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 */
+ /* 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);
+ 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;
}
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 */
+#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 */
+ 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
- }
+#endif
+
+#ifdef PMAP_PHOTONFLOW
+ case 'S':
+ /* Mean spherical / planar irrad switch for lightflow pmap */
+ check_bool(3, sphericalIrrad);
+ 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(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 04eced8..58add9a 100644
--- a/pmapparm.c
+++ b/pmapparm.c
@@ -1,116 +1,121 @@
#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 <ctype.h>
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) */
photonMaxDist = 0, /* Maximum cumulative distance of
photon path */
compressRatio = 0.9; /* Compression ratio for precomputed
contribution photon map */
#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
+ irrad from 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 params */
PhotonMapParams pmapParams [NUM_PMAP_TYPES] = {
{NULL, 0, 0, 0}, {NULL, 0, 0, 0}, {NULL, 0, 0, 0}, {NULL, 0, 0, 0},
{NULL, 0, 0, 0}, {NULL, 0, 0, 0}
#ifdef PMAP_PHOTONFLOW
, {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) -> srcContrib = NULL;
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 62117ea..9dfd6b5 100644
--- a/pmapparm.h
+++ b/pmapparm.h
@@ -1,87 +1,91 @@
/* 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 params per photon map from rpict/rtrace/rvu */
typedef struct {
char *fileName; /* Photon map file */
unsigned minGather, maxGather; /* Num photons to gather */
unsigned long distribTarget; /* Num photons to store */
} 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])
#ifdef PMAP_PHOTONFLOW
#define lightFlowPmapParams (pmapParams [PMAP_TYPE_LIGHTFLOW])
#endif
extern float pdfSamples, preDistrib, finalGather,
gatherTolerance, maxDistFix, pmapMaxDist,
photonMaxDist, compressRatio;
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
+ #endif
+
+ #ifdef PMAP_PHOTONFLOW
+ int sphericalIrrad;
+ #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