diff --git a/mkpmap.c b/mkpmap.c index a20e1a9..f44e7f0 100644 --- a/mkpmap.c +++ b/mkpmap.c @@ -1,991 +1,991 @@ #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 "pmaproi.h" #include "pmapio.h" #include "paths.h" #include "ambient.h" #include "resolu.h" #include "source.h" #include #include #include /* 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 contribTab = LU_SINIT(NULL, chickenMcFreeNuggetz); #if 0 char srcMod [MAXMODLIST]; #endif unsigned numContribs = 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("-aph mod\t\t\t\t# polyhedral region of interest"); puts("-apH file\t\t\t\t# polyhedral region of interest"); 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("-apc file nPhotons\t\t\t# caustic photon map"); puts("-apC file nPhotons bwidth compression" "\t# precomputed contribution photon map" ); puts("-apd file nPhotons\t\t\t# direct photon map"); puts("-apg file nPhotons\t\t\t# global photon map"); puts("-app file nPhotons bwidth\t\t# precomputed global photon map"); /* Transient pmap currently only supported by kd-tree */ #ifndef PMAP_OOC puts("-apt file nPhotons velocity\t\t# transient photon map"); #endif #ifdef PMAP_PHOTONFLOW /* Transient pmap currently only supported by kd-tree */ #ifndef PMAP_OOC puts("-apT file nPhotons velocity\t\t# transient light flow photon map"); #endif puts("-apV file nPhotons\t\t\t# light flow photon map"); #endif puts("-apv file nPhotons\t\t\t# volume photon map"); 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"); puts("-bn\t\t\t\t\t# contribution binning resolution"); 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); puts("-e expr\t\t\t\t\t# init expression for contrib binning"); printf("-ef %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", photonMaxPathDist); 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 puts("-p\t\t\t\t\t# per-modifier contrib binning params"); 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, binRes = 0; char **portLp = photonPortList, **sensLp = photonSensorList, **roiModLp = pmapROImodList, **amblp = NULL, sbuf [MAXSTR], portFlags [2] = "\0\0", sensFlags [2] = "\0\0", *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(); #if 0 setcontext(RCCONTEXT); #endif /* Compile default orientation variables for contribution binning */ scompile(PMAP_CONTRIB_SCDEFAULTS, NULL, 0); #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 */ #ifdef PMAP_CONTRIB 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[" ); #else error(INTERNAL, "sorry, contribution photon map is currently " "UNDER CONSTRUCTION. The management apologises " "for any inconvenience." ); #endif break; /* Transient pmap currently only supported by kd-tree */ #ifndef PMAP_OOC case 't': /* Transient photon map */ checkNumPhotons(i + 2); check(4, "ssf"); transientPmapParams.fileName = argv [++i]; transientPmapParams.distribTarget = parseMultiplier(argv [++i]); if (!transientPmapParams.distribTarget) goto badopt; photonVelocity = atof(argv [++i]); if (photonVelocity < FTINY) error(USER, "can't halt or reverse time! " "We're not Steven bloody Hawking, you know!" ); break; #endif #ifdef PMAP_PHOTONFLOW /* Light flow is a variant of volume photon map */ case 'V': /* Light flow photon map */ checkNumPhotons(i + 2); check(4, "ss"); lightFlowParams.fileName = argv [++i]; lightFlowParams.distribTarget = parseMultiplier(argv [++i]); if (!lightFlowParams.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; /* Transient pmap currently only supported by kd-tree */ #ifndef PMAP_OOC case 'T': /* Transient light flow photon map */ checkNumPhotons(i + 2); check(4, "ssf"); transLightFlowParams.fileName = argv [++i]; transLightFlowParams.distribTarget = parseMultiplier(argv [++i]); if (!transLightFlowParams.distribTarget) goto badopt; photonVelocity = atof(argv [++i]); if (photonVelocity < FTINY) error(USER, "can't halt or reverse time! " "We're not Steven bloody Hawking, you know!" ); /* Set zero absorption and forward scattering for global mist; extinction is left up to user or defined as local mist boundary in octree */ /* setcolor(cextinction, 1, 1, 1); */ setcolor(salbedo, 1, 1, 1); seccg = 1; break; #endif #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 'h': case 'H': /* Add polyhedral region(s) of interest by modifier(s) */ check(4, "s"); if (isupper(argv[i][3])) { /* Add ROI modifiers from file */ rval = wordfile(roiModLp, MAXSET - (roiModLp - pmapROImodList), getpath(argv [++i], getrlibpath(), R_OK) ); if (rval < 0) { sprintf(errmsg, "cannot open ROI modifier file %s", argv [i] ); error(SYSTEM, errmsg); } roiModLp += rval; } else { /* Add modifier string to list, mark end of list * with NULL */ *roiModLp++ = savqstr(argv [++i]); *roiModLp = NULL; } break; 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': /* Photon precomputation ratio */ check(4, "f"); finalGather = atof(argv [++i]); if (finalGather <= 0 || finalGather > 1) error(USER, "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 orientation flags to all modifier * strings; 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 { /* HACK: append flags to modifier string, add to * port list and mark end of 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 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_SENSFWD, PMAP_SENSBWD, sensFlags [0] ); 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); } /* HACK: append orientation flags to all modifier * strings; note this requires reallocation */ for (; rval--; sensLp++) { j = strlen(*sensLp); if (!(*sensLp = realloc(*sensLp, j + 2))) { sprintf(errmsg, "cannot allocate antimatter sensor " "modifiers from file %s", argv [i] ); error(SYSTEM, errmsg); } strcat(*sensLp, sensFlags); } } else { /* HACK: append flags to modifier string, add to * sensor list and mark of end list with NULL */ strcpy(sbuf, argv [++i]); strcat(sbuf, sensFlags); *sensLp++ = savqstr(sbuf); *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': /* Bin resolution for precomp contrib pmap */ check(3, "i"); binRes = atoi(argv [++i]); if (binRes <= 1 || binRes > PMAP_CONTRIB_MAXBINRES) goto badopt; break; default: goto badopt; } 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(3, "s"); diagFile = argv [++i]; break; default: /* Functional expression for precomputed contrib. pmap, e.g. to set constants for disk2square.cal */ check(2, "s"); scompile(argv [++i], NULL, 0); } break; case 'f': switch (argv [i][2]) { case 'o': /* Force overwrite */ check_bool(3, clobber); break; default: goto badopt; } break; #ifdef PMAP_EKSPERTZ case 'l': /* Limits */ switch (argv [i][2]) { case 'd': /* Limit photon path distance */ check(3, "f"); photonMaxPathDist = atof(argv [++i]); if (photonMaxPathDist <= 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"); addContribModifier(&contribTab, &numContribs, - argv [++i], binParm, NULL, PMAP_CONTRIB_SCBINS(binRes) + argv [++i], binParm, PMAP_CONTRIB_SCBINS(binRes) ); break; } break; case 'M': /* Modifier file for precomputed contribution photon map */ check(2, "s"); addContribModfile(&contribTab, &numContribs, - argv [++i], binParm, NULL, PMAP_CONTRIB_SCBINS(binRes) + argv [++i], binParm, PMAP_CONTRIB_SCBINS(binRes) ); 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"); 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, &contribTab, numContribs, 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/pmap.c b/pmap.c index 1b3301d..1042896 100644 --- a/pmap.c +++ b/pmap.c @@ -1,929 +1,929 @@ #ifndef lint static const char RCSid[] = "$Id: pmap.c,v 2.18 2021/02/19 02:10:35 rschregle Exp $"; #endif /* ====================================================================== Photon map main module 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.c,v 2.18 2021/02/19 02:10:35 rschregle Exp $ */ #include "pmap.h" #include "pmapmat.h" #include "pmapsrc.h" #include "pmaprand.h" #include "pmaproi.h" #include "pmapio.h" #include "pmapdens.h" #include "pmapbias.h" #include "pmapdiag.h" #include "pmutil.h" #include "otypes.h" #include "otspecial.h" #include #if NIX #include #include #include #include #endif static int photonParticipate (RAY *ray) /* Trace photon through participating medium. Returns 1 if passed through, or 0 if absorbed and $*%&ed. Analogon to rayparticipate(). */ { int i; RREAL xi1, cosTheta, phi, du, dv, rmax0; const RREAL cext = colorAvg(ray -> cext), invCext = 1 / cext, albedo = colorAvg(ray -> albedo), invAlbedo = 1 / albedo, gecc = ray -> gecc, gecc2 = sqr(gecc); FVECT u, v; COLOR cvext; /* Save incoming ray's max distance; this is required since rmax is set to the mean free distance for each path segment (see below) */ rmax0 = ray -> rmax; /* Mean free distance until interaction with medium */ ray -> rmax = -log(pmapRandom(mediumState)) / cext; while (!localhit(ray, &thescene)) { /* No geometry hit within mean free dist ray->rmax; photon interacts with medium */ if (!incube(&thescene, ray -> rop)) { /* Terminate photon if it has leaked from the scene */ #ifdef DEBUG_PMAP_LEAKED fprintf( stderr, "Volume photon leaked from scene at [%.3f %.3f %.3f]\n", ray -> rop [0], ray -> rop [1], ray -> rop [2] ); #endif return 0; } /* RGB attenuation over mean free distance */ setcolor(cvext, exp(-ray -> rmax * ray -> cext [0]), exp(-ray -> rmax * ray -> cext [1]), exp(-ray -> rmax * ray -> cext [2]) ); /* Modify ray color and normalise */ multcolor(ray -> rcol, cvext); colorNorm(ray -> rcol); /* Set ray origin for next interaction */ VCOPY(ray -> rorg, ray -> rop); /* Update ray's path length; terminate photon if maximum reached. See also photonRay(). */ rmax0 -= ray -> rot; if (photonMaxPathDist > 0 && rmax0 < 0) return 0; /* Store volume photon if nonzero albedo; this also accounts for direct inscattering from light sources */ if (albedo > FTINY) { /* Set ray's path length to store with photon */ ray -> rmax = rmax0; #ifdef PMAP_PHOTONFLOW if (lightFlowPhotonMapping) { /* Temporarily correct normalised flux for lightflow photons based on extinction (= mean num photons deposited per unit path length). */ scalecolor(ray -> rcol, invCext); newPhoton(lightFlowPmap, ray); /* Undo correction for next iteration after storing photon */ scalecolor(ray -> rcol, cext); } else #endif newPhoton(volumePmap, ray); } /* Colour bleeding without attenuation (?) */ multcolor(ray -> rcol, ray -> albedo); scalecolor(ray -> rcol, invAlbedo); if (pmapRandom(rouletteState) < albedo) { /* Scatter photon */ xi1 = pmapRandom(scatterState); cosTheta = ray -> gecc <= FTINY ? 2 * xi1 - 1 : 0.5 / gecc * ( 1 + gecc2 - sqr((1 - gecc2) / (1 + gecc * (2 * xi1 - 1))) ); phi = 2 * PI * pmapRandom(scatterState); du = dv = sqrt(1 - sqr(cosTheta)); /* sin(theta) */ du *= cos(phi); dv *= sin(phi); /* Get axes u & v perpendicular to photon direction */ i = 0; do { v [0] = v [1] = v [2] = 0; v [i++] = 1; fcross(u, v, ray -> rdir); } while (normalize(u) < FTINY); fcross(v, ray -> rdir, u); for (i = 0; i < 3; i++) ray -> rdir [i] = du * u [i] + dv * v [i] + cosTheta * ray -> rdir [i]; } ray -> rlvl++; ray -> rmax = -log(pmapRandom(mediumState)) / cext; /* fprintf(stderr, "%.3f\n", ray -> rmax); */ } /* Passed through medium until intersecting local object */ setcolor(cvext, exp(-ray -> rot * ray -> cext [0]), exp(-ray -> rot * ray -> cext [1]), exp(-ray -> rot * ray -> cext [2]) ); /* Modify ray color and normalise */ multcolor(ray -> rcol, cvext); colorNorm(ray -> rcol); /* Restore outgoing ray's max distance */ ray -> rmax = rmax0; return 1; } void tracePhoton (RAY *ray) /* Follow photon as it bounces around... */ { long mod; OBJREC *mat, *port = NULL; if (!ray -> parent) { /* !!! PHOTON PORT REJECTION SAMPLING HACK: get photon port for * !!! primary ray from ray -> ro, then reset the latter to NULL so * !!! as not to interfere with localhit() */ port = ray -> ro; ray -> ro = NULL; } if (ray -> rlvl > photonMaxBounce) { #ifdef PMAP_RUNAWAY_WARN error(WARNING, "runaway photon!"); #endif return; } if (colorAvg(ray -> cext) > FTINY && !photonParticipate(ray)) return; /* NOTE: localhit() limits intersections to the aft plane's distance ray->rmax. This mechanism is (ab)used here to limit the photon path length by initialising ray->rmax to photonMaxPathDist; see emitPhoton(). With unlimited path length (photonMaxPathDist=0), ray->rmax becomes negative (see photonRay()), and the aft plane has no effect. */ if (localhit(ray, &thescene)) { mod = ray -> ro -> omod; if (port && ray -> ro != port) { /* Terminate photon if emitted from port without intersecting it */ #ifdef PMAP_PORTREJECT_WARN sprintf(errmsg, "photon outside port %s", ray -> ro -> oname); error(WARNING, errmsg); #endif return; } if ((ray -> clipset && inset(ray -> clipset, mod)) || mod == OVOID) { /* Transfer ray if modifier is VOID or clipped within antimatta */ RAY tray; photonRay(ray, &tray, PMAP_XFER, NULL); tracePhoton(&tray); } else { /* Scatter for modifier material */ mat = objptr(mod); photonScatter [mat -> otype] (mat, ray); } } } static void preComputeGlobal (PhotonMap *pmap) /* Precompute irradiance from global photons for final gathering for a random subset of finalGather * pmap -> numPhotons photons, and builds the photon map, discarding the original photons. */ /* !!! NOTE: PRECOMPUTATION WITH OOC CURRENTLY WITHOUT CACHE !!! */ { unsigned long i, numPreComp; unsigned j; PhotonIdx pIdx; Photon photon; RAY ray; PhotonMap nuPmap; repComplete = numPreComp = finalGather * pmap -> numPhotons; if (verbose) { sprintf(errmsg, "\nPrecomputing irradiance for %ld global photons\n", numPreComp); eputs(errmsg); #if NIX fflush(stderr); #endif } /* Copy photon map for precomputed photons */ memcpy(&nuPmap, pmap, sizeof(PhotonMap)); /* Zero counters, init new heap and extents */ nuPmap.numPhotons = 0; initPhotonHeap(&nuPmap); for (j = 0; j < 3; j++) { nuPmap.minPos [j] = FHUGE; nuPmap.maxPos [j] = -FHUGE; } /* Record start time, baby */ repStartTime = time(NULL); #ifdef SIGCONT signal(SIGCONT, pmapPreCompReport); #endif repProgress = 0; photonRay(NULL, &ray, PRIMARY, NULL); ray.ro = NULL; for (i = 0; i < numPreComp; i++) { /* Get random photon from stratified distribution in source heap to * avoid duplicates and clustering */ pIdx = firstPhoton(pmap) + (unsigned long)((i + pmapRandom(pmap -> randState)) / finalGather); getPhoton(pmap, pIdx, &photon); /* Init dummy photon ray with intersection at photon position */ VCOPY(ray.rop, photon.pos); for (j = 0; j < 3; j++) ray.ron [j] = photon.norm [j] / 127.0; /* Get density estimate at photon position */ photonDensity(pmap, &ray, ray.rcol); /* Append photon to new heap from ray */ newPhoton(&nuPmap, &ray); /* Update progress */ repProgress++; if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) pmapPreCompReport(); #ifdef SIGCONT else signal(SIGCONT, pmapPreCompReport); #endif } /* Flush heap */ flushPhotonHeap(&nuPmap); #ifdef SIGCONT signal(SIGCONT, SIG_DFL); #endif /* Trash original pmap, replace with precomputed one */ deletePhotons(pmap); memcpy(pmap, &nuPmap, sizeof(PhotonMap)); if (verbose) { eputs("\nRebuilding precomputed photon map\n"); #if NIX fflush(stderr); #endif } /* Rebuild underlying data structure, destroying heap */ buildPhotonMap(pmap, NULL, NULL, 1); } typedef struct { unsigned long numPhotons [NUM_PMAP_TYPES], numEmitted, numComplete; } PhotonCnt; /* Photon distribution counter update interval expressed as bitmask; counters shared among subling subprocesses will only be updated in multiples of PMAP_CNTUPDATE in order to reduce contention */ #define PMAP_CNTUPDATE 0xffL void distribPhotons (PhotonMap **pmaps, unsigned numProc) { EmissionMap emap; char errmsg2 [128], shmFname [PMAP_TMPFNLEN]; unsigned t, srcIdx, proc; double totalFlux = 0; int shmFile, stat, pid; PhotonMap *pm; PhotonCnt *photonCnt, lastPhotonCnt; pid_t procPids [PMAP_MAXPROC]; for (t = 0; t < NUM_PMAP_TYPES && !pmaps [t]; t++); if (t >= NUM_PMAP_TYPES) error(USER, "no photon maps defined in distribPhotons"); if (!nsources) error(USER, "no light sources in distribPhotons"); /* =================================================================== * INITIALISATION - Set up emission and scattering funcs * =================================================================== */ emap.samples = NULL; emap.maxPartitions = MAXSPART; emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1); if (!emap.partitions) error(INTERNAL, "can't allocate source partitions in distribPhotons"); /* Initialise all defined photon maps */ for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t]) { initPhotonMap(pmaps [t], t); /* Open photon heapfile */ initPhotonHeap(pmaps [t]); /* Per-subprocess target count */ pmaps [t] -> distribTarget /= numProc; if (!pmaps [t] -> distribTarget) error(INTERNAL, "no photons to distribute in distribPhotons"); } initPhotonEmissionFuncs(); initPhotonScatterFuncs(); /* Get photon ports from modifier list */ getPhotonPorts(photonPortList); /* Get photon sensor modifiers */ getPhotonSensors(photonSensorList); /* Get polyhedral regions of interest */ getPolyROIs(pmapROImodList); #if NIX /* Set up shared mem for photon counters (zeroed by ftruncate) */ strcpy(shmFname, PMAP_TMPFNAME); shmFile = mkstemp(shmFname); if (shmFile < 0 || ftruncate(shmFile, sizeof(*photonCnt)) < 0) error(SYSTEM, "failed shared mem init in distribPhotons"); photonCnt = mmap(NULL, sizeof(*photonCnt), PROT_READ | PROT_WRITE, MAP_SHARED, shmFile, 0 ); if (photonCnt == MAP_FAILED) error(SYSTEM, "failed mapping shared memory in distribPhotons"); #else /* Allocate photon counters statically on Windoze */ if (!(photonCnt = malloc(sizeof(PhotonCnt)))) error(SYSTEM, "failed trivial malloc in distribPhotons"); photonCnt -> numEmitted = photonCnt -> numComplete = 0; #endif /* NIX */ /* Zero overflow tracking counters */ memset(&lastPhotonCnt, 0, sizeof(lastPhotonCnt)); if (verbose) { sprintf(errmsg, "\nIntegrating flux from %d sources", nsources); if (photonPorts) { sprintf(errmsg2, " via %d ports", numPhotonPorts); strcat(errmsg, errmsg2); } strcat(errmsg, "\n"); eputs(errmsg); } /* =================================================================== * FLUX INTEGRATION - Get total photon flux from light sources * =================================================================== */ for (srcIdx = 0; srcIdx < nsources; srcIdx++) { unsigned portCnt = 0; emap.src = source + srcIdx; do { /* Set next photon port if defined; note this loop iterates once if * no ports are defined. */ emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt : NULL; photonPartition [emap.src -> so -> otype] (&emap); if (verbose) { sprintf(errmsg, "\tIntegrating flux from source %s ", source [srcIdx].so -> oname ); if (emap.port) { sprintf(errmsg2, "via port %s ", photonPorts [portCnt].so -> oname ); strcat(errmsg, errmsg2); } sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions); strcat(errmsg, errmsg2); eputs(errmsg); #if NIX fflush(stderr); #endif } for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions; emap.partitionCnt++ ) { initPhotonEmission(&emap, pdfSamples); totalFlux += colorAvg(emap.partFlux); } portCnt++; } while (portCnt < numPhotonPorts); } if (totalFlux < FTINY) error(USER, "zero flux from light sources"); /* Record start time for progress reports */ repStartTime = time(NULL); if (verbose) { sprintf(errmsg, "\nPhoton distribution @ %d procs\n", numProc); eputs(errmsg); } /* MAIN LOOP */ for (proc = 0; proc < numProc; proc++) { #if NIX if (!(pid = fork())) { /* SUBPROCESS ENTERS HERE; open and mmapped files inherited */ #else if (1) { /* No subprocess under Windoze */ #endif /* Local photon counters for this subprocess */ unsigned passCnt = 0, prePassCnt = 0; unsigned long lastNumPhotons [NUM_PMAP_TYPES]; unsigned long localNumEmitted = 0; /* Num photons emitted by this subprocess alone */ /* Seed RNGs from PID for decorellated photon distribution */ pmapSeed(randSeed + proc, partState); pmapSeed(randSeed + (proc + 1) % numProc, emitState); pmapSeed(randSeed + (proc + 2) % numProc, cntState); pmapSeed(randSeed + (proc + 3) % numProc, mediumState); pmapSeed(randSeed + (proc + 4) % numProc, scatterState); pmapSeed(randSeed + (proc + 5) % numProc, rouletteState); #ifdef DEBUG_PMAP /* Output child process PID after random delay to prevent corrupted * console output due to race condition */ usleep(1e6 * pmapRandom(rouletteState)); fprintf(stderr, "Proc %d: PID = %d (waiting 10 sec to attach debugger...)\n", proc, getpid() ); /* Allow time for debugger to attach to child process */ sleep(10); #endif for (t = 0; t < NUM_PMAP_TYPES; t++) lastNumPhotons [t] = 0; /* ============================================================= * 2-PASS PHOTON DISTRIBUTION * Pass 1 (pre): emit fraction of target photon count * Pass 2 (main): based on outcome of pass 1, estimate remaining * number of photons to emit to approximate target * count * ============================================================= */ do { double numEmit; if (!passCnt) { /* INIT PASS 1 */ /* Skip if no photons contributed after sufficient * iterations; make it clear to user which photon maps are * missing so (s)he can check geometry and materials */ if (++prePassCnt > maxPreDistrib) { sprintf(errmsg, "proc %d: too many prepasses", proc); for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t] && !pmaps [t] -> numPhotons) { sprintf(errmsg2, ", no %s photons stored", pmapName [t] ); strcat(errmsg, errmsg2); } sprintf(errmsg2, "; increase -apD"); strcat(errmsg, errmsg2); error(USER, errmsg); break; } /* Num to emit is fraction of minimum target count */ numEmit = FHUGE; for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t]) numEmit = min(pmaps [t] -> distribTarget, numEmit); numEmit *= preDistrib; } else { /* INIT PASS 2 */ /* Based on the outcome of the predistribution we can now * estimate how many more photons we have to emit for each * photon map to meet its respective target count. This * value is clamped to 0 in case the target has already been * exceeded in the pass 1. */ double maxDistribRatio = 0; /* Set the distribution ratio for each map; this indicates * how many photons of each respective type are stored per * emitted photon, and is used as probability for storing a * photon by newPhoton(). Since this biases the photon * density, newPhoton() promotes the flux of stored photons * to compensate. */ for (t = 0; t < NUM_PMAP_TYPES; t++) if ((pm = pmaps [t])) { pm -> distribRatio = (double)pm -> distribTarget / pm -> numPhotons - 1; /* Check if photon map "overflowed", i.e. exceeded its * target count in the prepass; correcting the photon * flux via the distribution ratio is no longer * possible, as no more photons of this type will be * stored, so notify the user rather than deliver * incorrect results. In future we should handle this * more intelligently by using the photonFlux in each * photon map to individually correct the flux after * distribution. */ if (pm -> distribRatio <= FTINY) { sprintf(errmsg, "%s photon map overflow in prepass, reduce -apD ", pmapName [t] ); error(INTERNAL, errmsg); } maxDistribRatio = max(pm -> distribRatio, maxDistribRatio ); } /* Normalise distribution ratios and calculate number of * photons to emit in main pass */ for (t = 0; t < NUM_PMAP_TYPES; t++) if ((pm = pmaps [t])) pm -> distribRatio /= maxDistribRatio; if ((numEmit = localNumEmitted * maxDistribRatio) < FTINY) /* No photons left to distribute in main pass */ break; } /* PHOTON DISTRIBUTION LOOP */ for (srcIdx = 0; srcIdx < nsources; srcIdx++) { unsigned portCnt = 0; emap.src = source + srcIdx; do { /* Set next photon port if defined; note this loop iterates * once if no ports are defined. */ emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt : NULL; photonPartition [emap.src -> so -> otype] (&emap); if (verbose && !proc) { /* Output from subproc 0 only to avoid race condition * on console I/O */ if (!passCnt) sprintf(errmsg, "\tPREPASS %d on source %s ", prePassCnt, source [srcIdx].so -> oname ); else sprintf(errmsg, "\tMAIN PASS on source %s ", source [srcIdx].so -> oname ); if (emap.port) { sprintf(errmsg2, "via port %s ", photonPorts [portCnt].so -> oname ); strcat(errmsg, errmsg2); } sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions ); strcat(errmsg, errmsg2); eputs(errmsg); #if NIX fflush(stderr); #endif } for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions; emap.partitionCnt++ ) { double partNumEmit; unsigned long partEmitCnt; /* Get photon origin within current source partishunn * and build emission map */ photonOrigin [emap.src -> so -> otype] (&emap); initPhotonEmission(&emap, pdfSamples); /* Number of photons to emit from ziss partishunn -- * proportional to flux; photon ray weight and scalar * flux are uniform (latter only varying in RGB). */ partNumEmit = numEmit * colorAvg(emap.partFlux) / totalFlux; partEmitCnt = (unsigned long)partNumEmit; /* Probabilistically account for fractional photons */ if (pmapRandom(cntState) < partNumEmit - partEmitCnt) partEmitCnt++; /* Integer counter avoids FP rounding errors during * iteration */ while (partEmitCnt--) { RAY photonRay; /* Emit photon based on PDF, skip if rejected. */ /* NOTE: rejection sampling skips incrementing the * emission counter (see below), thus compensating * for the rejected photons by increasing the photon * flux in proportion to the lower effective * emission count. * BUG: THIS INTERFERES WITH THE PROGRESS COUNTER * REPORTED TO THE PARENT, AND WITH THE * PREDISTRIBUTION PASS --> PHOTON DISTRIBUTION WILL * FINISH EARLY, WITH FEWER PHOTONS THAN TARGETTED! */ if (!emitPhoton(&emap, &photonRay)) continue; /* Set subprocess ID in photonRay; extends path ID. Used by newPhoton() as photon attributes. */ PMAP_SETRAYPROC(&photonRay, proc); tracePhoton(&photonRay); /* Update local emission counter*/ localNumEmitted++; if (!(localNumEmitted & PMAP_CNTUPDATE)) { /* Update global counters shared with siblings in intervals to reduce overhead and contention */ lastPhotonCnt.numEmitted = photonCnt -> numEmitted; photonCnt -> numEmitted += PMAP_CNTUPDATE; lastPhotonCnt.numComplete = photonCnt -> numComplete; photonCnt -> numComplete += PMAP_CNTUPDATE; /* Check for overflow using previous values */ if (photonCnt -> numEmitted < lastPhotonCnt.numEmitted || photonCnt -> numComplete < lastPhotonCnt.numComplete ) { sprintf(errmsg, "photon emission counter " "overflow (was: %ld, is: %ld)", lastPhotonCnt.numEmitted, photonCnt -> numEmitted ); error(INTERNAL, errmsg); } /* Update global photon count for each pmap */ for (t = 0; t < NUM_PMAP_TYPES; t++) { if (pmaps [t]) { lastPhotonCnt.numPhotons [t] = photonCnt -> numPhotons [t]; - /* Differential increment using local + /* Differentially increment using local counter lastNumPhotons */ photonCnt -> numPhotons [t] += pmaps [t] -> numPhotons - lastNumPhotons [t]; lastNumPhotons [t] = pmaps [t] -> numPhotons; /* Check for photon counter overflow (this * could only happen before an emission * counter overflow if the scene has an * absurdly high albedo and/or very dense * geometry) */ if (photonCnt -> numPhotons [t] < lastPhotonCnt.numPhotons [t] ) { sprintf(errmsg, "photon counter " "overflow (was: %ld, is: %ld)", lastPhotonCnt.numPhotons [t], photonCnt -> numPhotons [t] ); error(INTERNAL, errmsg); } } } } } #if !NIX /* Synchronous progress report on Windoze */ if (!proc && photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) { repEmitted = repProgress = photonCnt -> numEmitted; repComplete = photonCnt -> numComplete; pmapDistribReport(); } #endif } portCnt++; } while (portCnt < numPhotonPorts); } for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t] && !pmaps [t] -> numPhotons) { /* Double preDistrib in case a photon map is empty and * redo pass 1 --> possibility of infinite loop for * pathological scenes (e.g. absorbing materials) */ preDistrib *= 2; break; } if (t >= NUM_PMAP_TYPES) /* No empty photon maps found; now do pass 2 */ passCnt++; } while (passCnt < 2); /* Flush heap buffa for every pmap one final time; * avoids potential data corruption! */ for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t]) { flushPhotonHeap(pmaps [t]); /* Heap file closed automatically on exit fclose(pmaps [t] -> heap); */ #ifdef DEBUG_PMAP sprintf(errmsg, "Proc %d: total %ld photons\n", proc, pmaps [t] -> numPhotons ); eputs(errmsg); #endif } #if NIX /* Terminate subprocess */ exit(0); #endif } else { /* PARENT PROCESS CONTINUES LOOP ITERATION HERE */ if (pid < 0) error(SYSTEM, "failed to fork subprocess in distribPhotons"); else /* Save child process IDs */ procPids [proc] = pid; } } #if NIX /* PARENT PROCESS CONTINUES HERE */ #ifdef SIGCONT /* Enable progress report signal handler */ signal(SIGCONT, pmapDistribReport); #endif /* Wait for subprocesses to complete while reporting progress */ proc = numProc; while (proc) { while (waitpid(-1, &stat, WNOHANG) > 0) { /* Subprocess exited; check status */ if (!WIFEXITED(stat) || WEXITSTATUS(stat)) { /* Exited with error; terminate siblings, clean up & bail out */ for (proc = 0; proc < numProc; proc++) kill(procPids [proc], SIGKILL); /* Unmap shared memory, squish mapped file */ munmap(photonCnt, sizeof(*photonCnt)); close(shmFile); unlink(shmFname); error(USER, "failed photon distribution"); } --proc; } /* Nod off for a bit and update progress */ sleep(1); /* Asynchronous progress report from shared subprocess counters */ repEmitted = repProgress = photonCnt -> numEmitted; repComplete = photonCnt -> numComplete; repProgress = repComplete = 0; for (t = 0; t < NUM_PMAP_TYPES; t++) if ((pm = pmaps [t])) { /* Get global photon count from shmem updated by subprocs */ repProgress += pm -> numPhotons = photonCnt -> numPhotons [t]; repComplete += pm -> distribTarget; } repComplete *= numProc; if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) pmapDistribReport(); #ifdef SIGCONT else signal(SIGCONT, pmapDistribReport); #endif } #endif /* NIX */ /* =================================================================== * POST-DISTRIBUTION - Set photon flux and build data struct for photon * storage, etc. * =================================================================== */ #ifdef SIGCONT /* Reset signal handler */ signal(SIGCONT, SIG_DFL); #endif free(emap.samples); /* Set photon flux */ totalFlux /= photonCnt -> numEmitted; #if NIX /* Photon counters no longer needed, unmap shared memory */ munmap(photonCnt, sizeof(*photonCnt)); close(shmFile); unlink(shmFname); #else free(photonCnt); #endif if (verbose) eputs("\n"); for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t]) { if (verbose) { sprintf(errmsg, "Building %s photon map\n", pmapName [t]); eputs(errmsg); #if NIX fflush(stderr); #endif } /* Build underlying data structure; heap is destroyed */ buildPhotonMap(pmaps [t], &totalFlux, NULL, numProc); } /* Precompute photon irradiance if necessary */ if (preCompPmap) { if (verbose) eputs("\n"); preComputeGlobal(preCompPmap); } if (verbose) eputs("\n"); } diff --git a/pmapcontrib.c b/pmapcontrib.c index ffa6e2a..60415f4 100644 --- a/pmapcontrib.c +++ b/pmapcontrib.c @@ -1,831 +1,907 @@ #ifndef lint static const char RCSid[] = "$Id: pmapcontrib.c,v 2.20 2021/02/16 20:06:06 greg Exp $"; #endif /* ========================================================================= - Photon map routines for precomputed light source contributions; these - routines interface to mkpmap. + Photon map routines for precomputed light source contributions. + These routines interface to mkpmap. 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") ========================================================================= $Id: pmapcontrib.c,v 2.20 2021/02/16 20:06:06 greg Exp $ */ #include "pmapcontrib.h" #include "pmapdiag.h" #include "pmaprand.h" #include "pmapmat.h" #include "pmapsrc.h" #include "pmapio.h" #include "pmutil.h" #include "otspecial.h" #include "otypes.h" #include "lookup.h" extern void SDdisk2square(double sq[2], double diskx, double disky); static int logDim (unsigned size) /* Return log2(sqrt(size)) = l, where size = (2^l)(2^l). Is there a faster way (e.g. binary search)? */ { unsigned i, sz, dim; if (!size) return 0; for (i = 0, dim = 0, sz = size; sz >>= 2; dim++); return dim; } static int xy2bin (unsigned l, int x, int y) /* Serialise 2D coordinates in range (2^l) x (2^l) to 1D bin. Returns -1 if coordinates invalid */ { return x < 0 || y < 0 ? -1 : (x << l) + y; } static void bin2xy (unsigned l, int bin, int *x, int *y) /* Deserialise 1D bin to 2D coordinates in range (2^l) x (2^l). Returns -1 if bin invalid */ { /* Obviously this is faster than integer division/modulo */ *x = bin < 0 ? -1 : bin >> l; *y = bin < 0 ? -1 : bin & PMAP_CONTRIB_SCDIM(l) - 1; } static int ray2bin (const RAY *ray, unsigned nbins) /* Map ray dir (pointing away from origin) to its 1D Shirley-Chiu bin, where nbins = (2^l) x (2^l) for l > 1. Returns -1 if mapped bin is invalid (e.g. behind plane) */ { const unsigned l = logDim(nbins), scDim = PMAP_CONTRIB_SCDIM(l); static int scRHS, varInit = 0; static FVECT scNorm, scUp; unsigned scBin [2]; FVECT diskPlane; RREAL dz, diskx, disky, rad, diskd2, scCoord [2]; if (!varInit) { /* Lazy init shirley-Chiu mapping orientation from function variables */ scRHS = varvalue(PMAP_CONTRIB_SCRHS); scNorm [0] = varvalue(PMAP_CONTRIB_SCNORMX); scNorm [1] = varvalue(PMAP_CONTRIB_SCNORMY); scNorm [2] = varvalue(PMAP_CONTRIB_SCNORMZ); scUp [0] = varvalue(PMAP_CONTRIB_SCUPX); scUp [1] = varvalue(PMAP_CONTRIB_SCUPY); scUp [2] = varvalue(PMAP_CONTRIB_SCUPZ); varInit ^= 1; } /* Map incident direction to disk */ dz = DOT(ray -> rdir, scNorm); /* normal proj */ if (dz > 0) { fcross(diskPlane, scUp, scNorm); /* y-axis in plane, perp to up */ diskPlane [0] *= scRHS; diskx = DOT(ray -> rdir, diskPlane); disky = DOT(ray -> rdir, scUp) - dz * DOT(ray -> rdir, scUp); diskd2 = diskx * diskx + disky * disky; /* in-plane length^2 */ rad = dz>FTINY ? sqrt((1 - dz*dz) / diskd2) : 0; /* radial factor */ disky *= rad; /* diskx, disky now in range [-1, 1] */ /* Apply Shirley-Chiu mapping of (diskx, disky) to square */ SDdisk2square(scCoord, diskx, disky); /* Map Shirley-Chiu square coords to 1D bin */ scBin [0] = scCoord [0] * scDim; scBin [1] = scCoord [1] * scDim; return xy2bin(l, scBin [0], scBin [1]); } else return -1; } #ifndef PMAP_CONTRIB_TEST MODCONT *addContribModifier (LUTAB *contribTab, unsigned *numContribs, - char *mod, char *binParm, char *binVal, int binCnt + char *mod, char *binParm, int binCnt ) /* Add light source modifier mod to contribution lookup table contribTab, and update numContribs. Return initialised contribution data for this modifier. This code adapted from rcontrib.c:addmodifier(). */ { LUENT *lutEntry = lu_find(contribTab, mod); MODCONT *contrib; EPNODE *eBinVal; if (lutEntry -> data) { /* Reject duplicate modifiers */ sprintf(errmsg, "duplicate light source modifier %s", mod); error(USER, errmsg); } if (*numContribs >= MAXMODLIST) { sprintf(errmsg, "too many modifiers (max. %d)", MAXMODLIST); error(INTERNAL, errmsg); } lutEntry -> key = mod; - if (!binVal) - /* Fall back to single bin if no value specified */ - binVal = "0"; - eBinVal = eparse(binVal); - if (eBinVal -> type == NUM) { - /* Bin value is constant */ - binCnt = (int)(evalue(eBinVal) + 1.5); - if (binCnt != 1) { - sprintf(errmsg, "invalid non-zero constant for bin %s", binVal); - error(USER, errmsg); - } - } - else if (binCnt <= 0) { + if (binCnt <= 0) { sprintf(errmsg, "undefined/invalid bin count for modifier %s", mod); error(USER, errmsg); } /* Allocate and init contributions */ contrib = (MODCONT *)malloc( sizeof(MODCONT) + sizeof(DCOLOR) * (binCnt - 1) ); if (!contrib) error(SYSTEM, "out of memory in addContribModifier()"); contrib -> outspec = NULL; contrib -> modname = mod; contrib -> params = binParm; contrib -> nbins = binCnt; contrib -> binv = eBinVal; contrib -> bin0 = 0; memset(contrib -> cbin, 0, sizeof(DCOLOR) * binCnt); lutEntry -> data = lutEntry -> data = (char *)contrib; ++(*numContribs); return(contrib); } void addContribModfile (LUTAB *contribTab, unsigned *numContribs, - char *modFile, char *binParm, char *binVal, int binCnt + char *modFile, char *binParm, int binCnt ) /* Add light source modifiers from file modFile to contribution lookup table * contribTab, and update numContribs. * NOTE: This code is adapted from rcontrib.c */ { char *mods [MAXMODLIST]; int i; /* Find file and store strings */ i = wordfile(mods, MAXMODLIST, getpath(modFile, getrlibpath(), R_OK)); if (i < 0) { sprintf(errmsg, "can't open modifier file %s", modFile); error(SYSTEM, errmsg); } if (*numContribs + i >= MAXMODLIST - 1) { sprintf(errmsg, "too many modifiers (max. %d) in file %s", MAXMODLIST - 1, modFile ); error(INTERNAL, errmsg); } for (i = 0; mods [i]; i++) /* Add each modifier */ - addContribModifier(contribTab, numContribs, mods [i], - binParm, binVal, binCnt + addContribModifier(contribTab, numContribs, + mods [i], binParm, binCnt ); } static int contribSourceBin (LUTAB *contribs, const RAY *ray) /* Map contribution source ray to its bin for light source ray -> rsrc, using Shirley-Chiu disk-to-square mapping. Return invalid bin -1 if mapping failed. */ { const SRCREC *src; const OBJREC *srcMod; const MODCONT *srcCont; RAY srcRay; int bin, i; /* Check we have a valid ray and contribution LUT */ if (!ray || !contribs) return -1; src = &source [ray -> rsrc]; srcMod = findmaterial(src -> so); srcCont = (MODCONT*)lu_find(contribs, srcMod -> oname) -> data; if (!srcCont) /* Not interested in this source (modifier not in contrib LUT) */ return -1; /* Set up shadow ray pointing to source for disk2square mapping */ rayorigin(&srcRay, SHADOW, NULL, NULL); srcRay.rsrc = ray -> rsrc; VCOPY(srcRay.rorg, ray -> rop); for (i = 0; i < 3; i++) srcRay.rdir [i] = -ray -> rdir [i]; if (!(src -> sflags & SDISTANT ? sourcehit(&srcRay) : (*ofun[srcMod -> otype].funp)(srcMod, &srcRay) )) /* (Redundant?) sanity check for valid source ray? */ return -1; worldfunc(RCCONTEXT, &srcRay); set_eparams((char*)srcCont -> params); if ((bin = ray2bin(&srcRay, srcCont -> nbins)) < 0) error(WARNING, "Ignoring invalid bin in contribSourceBin()"); return bin; } PhotonContribSourceIdx newPhotonContribSource (PhotonMap *pmap, const RAY *contribSrcRay, FILE *contribSrcHeap ) /* Add contribution source ray for emitted contribution photon and save * light source index and binned direction. The current contribution * source is stored in pmap -> lastContribSrc. If the previous contrib * source spawned photons (i.e. has srcIdx >= 0), it's appended to * contribSrcHeap. If contribSrcRay == NULL, the current contribution * source is still flushed, but no new source is set. Returns updated * contribution source counter pmap -> numContribSrc. */ { if (!pmap || !contribSrcHeap) return 0; /* Check if last contribution source has spawned photons (srcIdx >= 0, * see newPhoton()), in which case we save it to the heap file before * clobbering it. (Note this is short-term storage, so we doan' need * da portable I/O stuff here). */ if (pmap -> lastContribSrc.srcIdx >= 0) { if (!fwrite(&pmap -> lastContribSrc, sizeof(PhotonContribSource), 1, contribSrcHeap )) error(SYSTEM, "failed writing photon contrib source in " "newPhotonContribSource()" ); pmap -> numContribSrc++; if (!pmap -> numContribSrc || pmap -> numContribSrc > PMAP_MAXCONTRIBSRC ); } - /* Mark this contribution source unused with a negative source index - until its path spawns a photon (see newPhoton()) */ - pmap -> lastContribSrc.srcIdx = -1; - /* Map ray to bin in anticipation that this contribution source will - be used, since the ray will be lost once a photon is spawned */ - pmap -> lastContribSrc.srcBin = contribSourceBin( - pmap -> contribTab, contribSrcRay - ); - - if (pmap -> lastContribSrc.srcBin < 0) { - /* Warn if invalid bin, but trace photon nonetheless. It will - count as emitted to prevent bias, but will not be stored in - newPhoton(), as it contributes zero flux */ - sprintf(errmsg, "invalid bin for light source %s, " - "contribution photons will be discarded", - source [contribSrcRay -> rsrc].so -> oname + if (contribSrcRay) { + /* Mark this contribution source unused with a negative source + index until its path spawns a photon (see newPhoton()) */ + pmap -> lastContribSrc.srcIdx = -1; + + /* Map ray to bin in anticipation that this contrib source will be + used, since the ray will be lost once a photon is spawned */ + pmap -> lastContribSrc.srcBin = contribSourceBin( + pmap -> contribTab, contribSrcRay ); - error(WARNING, errmsg); + + if (pmap -> lastContribSrc.srcBin < 0) { + /* Warn if invalid bin, but trace photon nonetheless. It will + count as emitted to prevent bias, but will not be stored in + newPhoton(), as it contributes zero flux */ + sprintf(errmsg, "invalid bin for light source %s, " + "contribution photons will be discarded", + source [contribSrcRay -> rsrc].so -> oname + ); + error(WARNING, errmsg); + } } return pmap -> numContribSrc; } PhotonContribSourceIdx buildContribSources (PhotonMap *pmap, FILE **contribSrcHeap, char **contribSrcHeapFname, PhotonContribSourceIdx *contribSrcOfs, unsigned numHeaps ) - /* Consolidate per-subprocess contribution sources heaps into array + /* Consolidate per-subprocess contribution source heaps into array * pmap -> contribSrc. Returns offset for contribution source index * linearisation in pmap -> numContribSrc. The heap files in * contribSrcHeap are closed on return. */ { PhotonContribSourceIdx heapLen; unsigned heap; - + if (!pmap || !contribSrcHeap || !contribSrcOfs || !numHeaps) return 0; - + pmap -> numContribSrc = 0; - + for (heap = 0; heap < numHeaps; heap++) { contribSrcOfs [heap] = pmap -> numContribSrc; if (fseek(contribSrcHeap [heap], 0, SEEK_END) < 0) error(SYSTEM, "failed photon contrib source seek " "in buildContribSources()" ); pmap -> numContribSrc += heapLen = ftell(contribSrcHeap [heap]) / sizeof(PhotonContribSource); if (!(pmap -> contribSrc = realloc(pmap -> contribSrc, pmap -> numContribSrc * sizeof(PhotonContribSource) ))) error(SYSTEM, "failed photon contrib source alloc " "in buildContribSources()" ); rewind(contribSrcHeap [heap]); if (fread(pmap -> contribSrc + contribSrcOfs [heap], sizeof(PhotonContribSource), heapLen, contribSrcHeap [heap] ) != heapLen) error(SYSTEM, "failed reading photon contrib source " "in buildContribSources()" ); fclose(contribSrcHeap [heap]); unlink(contribSrcHeapFname [heap]); } return pmap -> numContribSrc; } - static MODCONT *getPhotonContrib (PhotonMap *pmap, RAY *ray, COLOR irrad) - /* Locate photons near ray -> rop which originated from the light source - modifier indexed by ray -> rsrc, and accumulate their contributions - in pmap -> srcContrib. Return total contributions in irrad and - pointer to binned contributions (or NULL if none were found). */ - { - const OBJREC *srcMod = - findmaterial(source [ray -> rsrc].so); - MODCONT *contrib = - (MODCONT*)lu_find(pmap -> contribTab, srcMod -> oname) -> data; - Photon *photon; - COLOR flux; - PhotonSearchQueueNode *sqn; - float r2, norm; - unsigned i; - - /* Zero bins for this source modifier */ - memset(contrib -> cbin, 0, sizeof(DCOLOR) * contrib -> nbins); - setcolor(irrad, 0, 0, 0); - - if (!contrib || !pmap -> maxGather) - /* Modifier not in LUT or zero bandwidth */ - return NULL; - - /* Lookup photons */ - pmap -> squeue.tail = 0; - /* Pass light source index to filter in findPhotons() via - pmap -> lastContribSrc */ - pmap -> lastContribSrc.srcIdx = ray -> rsrc; - 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 NULL; - } - - /* Avg radius^2 between furthest two photons to improve accuracy */ - sqn = pmap -> squeue.node + 1; - r2 = max(sqn -> dist2, (sqn + 1) -> dist2); - r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2)); - - /* XXX: OMIT extra normalisation factor 1/PI for ambient calc? */ - #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 */ - norm = 2 / (PI * PI * r2); - #else - /* Normalise accumulated flux by search area PI * r^2, including - RADIANCE-specific lambertian factor PI */ - norm = 1 / (PI * PI * r2); - #endif - - /* Skip the extra photon */ - for (i = 1 ; i < pmap -> squeue.tail; i++, sqn++) { - /* Get photon's contribution to density estimate */ - photon = getNearestPhoton(&pmap -> squeue, sqn -> idx); - getPhotonFlux(photon, flux); - scalecolor(flux, norm); - #ifdef PMAP_EPANECHNIKOV - /* Apply Epanechnikov kernel to photon flux based on photon distance */ - scalecolor(flux, 1 - sqn -> dist2 / r2); - #endif - addcolor(irrad, flux); - addcolor(contrib -> cbin [photonSrcBin(pmap, photon)], flux); - } - - return contrib; - } - - - - static int CoeffCompare (const void *c1, const void *c2) + static int coeffCompare (const void *c1, const void *c2) /* Comparison function to REVERSE sort thresholded coefficients */ { const ThreshWaveletCoeff *tcoeff1 = c1, *tcoeff2 = c2; /* TODO: USE MAX HERE INSTEAD OF AVERAGE? */ const WAVELET_COEFF v1 = colorAvg(tcoeff1 -> coeff), v2 = colorAvg(tcoeff2 -> coeff); if (v1 < v2) return 1; else if (v1 > v2) return -1; else return 0; } - + static int thresholdContribs (MODCONT *binnedContribs, PreComputedPhotonContrib *pcContribData) /* Threshold binned detail coefficients in &binnedContribs -> cbin [1] by keeping the (pcContribData -> nCompressedBins) largest and returning these in pcContribData -> threshCoeffs along with their original bin order. NOTE: binnedContribs -> cbin [0] is the average wavelet coefficient and excluded from thresholding. */ { - unsigned b, i; + unsigned b; ThreshWaveletCoeff *threshCoeffs; - mRGBERange *mrgbeRange; if (!pcContribData) return -1; - /* Init per-channel coefficient range */ - for (i = 0, mrgbeRange = &pcContribData -> mrgbeRange; i < 3; i++) { - mrgbeRange -> min [i] = FHUGE; - mrgbeRange -> max [i] = 0; - } - /* Skip 1st coeff; it's the average and therefore not thresholded */ for (b = 0, threshCoeffs = pcContribData -> threshCoeffs; - b < binnedContribs -> nbins - 1; b++ + b < binnedContribs -> nbins - 1; + b++ ) { /* Set up pointers to coeffs (sorted faster than 3 floats/doubles) and remember original bin position prior to sorting. - NOTE: The original bin position _excludes_ the 1st coeff, + NOTE: The original bin position *excludes* the 1st coeff, i.e. the average! */ - threshCoeffs [b].coeff = - (WAVELET_COEFF*)&(binnedContribs -> cbin [b + 1]); + threshCoeffs [b].coeff = (WAVELET_COEFF*)( + &(binnedContribs->cbin [b+1]) + ); threshCoeffs [b].bin = b; - - /* Update per-channel coefficient range */ - for (i = 0; i < 3; i++) { - if (threshCoeffs [b].coeff [i] < mrgbeRange -> min [i]) - mrgbeRange -> min [i] = threshCoeffs [b].coeff [i]; - if (threshCoeffs [b].coeff [i] > mrgbeRange -> max [i]) - mrgbeRange -> max [i] = threshCoeffs [b].coeff [i]; - } } /* REVERSE sort coeffs by magnitude; the non-thresholded coeffs will then start at threshCoeffs [0] */ - qsort(threshCoeffs, binnedContribs -> nbins, - sizeof(ThreshWaveletCoeff), CoeffCompare + qsort(threshCoeffs, binnedContribs -> nbins - 1, + sizeof(ThreshWaveletCoeff), coeffCompare ); return 0; } - - static mRGBE *encodeContribs (MODCONT *binnedContribs, - LUTAB *preCompContribLUT, float compressRatio + + static PreComputedPhotonContrib *encodeContribs ( + MODCONT *binnedContribs, LUTAB *preCompContribLUT, float compressRatio ) /* Apply wavelet transform to binnedContribs -> cbin and compress according to compressRatio, storing results in LUT preCompContribLUT - per modifier. Returns the encoded detail coefficients in - preCompCongtribLUT -> data -> mrgbeCoeffs, which is also the return - value. Note that the average coefficient is not encoded, and returned + per modifier. Returns struct with wavelet transformed, thresholded + and encoded detail coefficients. + Note that the average coefficient is not encoded, and returned in binnedContribs -> cbin [0]. Returns NULL on error. */ { - unsigned b; + unsigned b, i; + long nCompressedBins; LUENT *pcContribLUTnode; PreComputedPhotonContrib *pcContribData; ThreshWaveletCoeff *threshCoeffs; mRGBERange *mrgbeRange; mRGBE *mrgbeCoeff; if (!binnedContribs || !preCompContribLUT) return NULL; pcContribLUTnode = lu_find(preCompContribLUT, binnedContribs->modname); if (!pcContribLUTnode) error(SYSTEM, "out of LUT memory in encodeContrib()"); if (!pcContribLUTnode -> key) { /* New LUT entry for precomputed contributions for this modifier */ pcContribLUTnode -> key = (char*)binnedContribs -> modname; pcContribData = (PreComputedPhotonContrib*)( pcContribLUTnode -> data = malloc(sizeof(PreComputedPhotonContrib)) ); if (!pcContribData) error(SYSTEM, "out of memory allocating precomputed photon " "contributions in encodeContrib()" ); /* Get Shirley-Chiu square (= wavelet matrix) resolution, bins per dimension, and fixed number of compressed coeffs/bins (note we subtract one as the average coefficient is _not_ thresholded) */ pcContribData -> l = logDim(binnedContribs -> nbins); pcContribData -> scDim = PMAP_CONTRIB_SCDIM(pcContribData -> l); - pcContribData -> nCompressedBins = - (binnedContribs -> nbins - 1) * (1 - compressRatio); - - /* Lazily "allocate" primary wavelet matrix pointing to + nCompressedBins = (1 - compressRatio) * (binnedContribs->nbins - 1); + if (nCompressedBins > 0) + pcContribData -> nCompressedBins = nCompressedBins; + else + error(USER, "too few bins to compress in encodeContrib()"); + + /* Lazily "allocate" primary wavelet matrix pointing to array binnedContribs -> cbin; this imposes the 2D matrix structure - expected by the wavelet xform routines onto it. */ + expected by the wavelet xform routines onto the 1D array. */ pcContribData -> waveletMatrix = calloc( pcContribData -> scDim, sizeof(WaveletCoeff3*) ); if (!pcContribData -> waveletMatrix) error(SYSTEM, "out of memory allocating primary " "wavelet matrix in encodeContrib()" ); for (b = 0; b < pcContribData -> scDim; b++) - /* Point to each row in existing contrib array */ + /* Point to each row in existing 1D contrib array */ pcContribData -> waveletMatrix [b] = &binnedContribs -> cbin [b * pcContribData -> scDim]; /* Lazily allocate transposed wavelet matrix */ pcContribData -> tWaveletMatrix = allocWaveletMatrix(pcContribData -> l); if (!pcContribData -> tWaveletMatrix) error(SYSTEM, "out of memory allocating transposed " "wavelet matrix in encodeContrib()" ); /* Lazily allocate thresholded detail coefficients (minus the average coeff) */ pcContribData -> threshCoeffs = calloc( binnedContribs -> nbins - 1, sizeof(ThreshWaveletCoeff) ); /* Lazily allocate mRGBE-encoded, compressed wavelet coeffs */ pcContribData -> mrgbeCoeffs = calloc( pcContribData -> nCompressedBins, sizeof(mRGBE) ); if (!pcContribData -> mrgbeCoeffs) error(SYSTEM, "out of memory allocating mRGBE coeffs " "in encodeContrib()" ); } - else pcContribData = - (PreComputedPhotonContrib*)pcContribLUTnode -> data; + else pcContribData = (PreComputedPhotonContrib*)pcContribLUTnode->data; /* Do 2D wavelet transform on pcContribData -> waveletMatrix (which really maps to binnedContribs -> cbin) */ - if (waveletXform2(pcContribData -> waveletMatrix, - pcContribData -> tWaveletMatrix, pcContribData -> l - ) < 0) + if ( + waveletXform2(pcContribData -> waveletMatrix, + pcContribData -> tWaveletMatrix, pcContribData -> l + ) < 0 + ) error(INTERNAL, "failed wavelet transform in encodeContrib()"); /* Compress wavelet detail coeffs by thresholding */ if (thresholdContribs(binnedContribs, pcContribData) < 0) error(INTERNAL, "failed wavelet compression in encodeContrib()"); + threshCoeffs = pcContribData -> threshCoeffs; - /* Init mRGBE coefficient normalisation from range */ + /* Init per-channel coefficient range for mRGBE encoding */ mrgbeRange = &pcContribData -> mrgbeRange; + setcolor(mrgbeRange -> min, FHUGE, FHUGE, FHUGE); + setcolor(mrgbeRange -> max, 0, 0, 0); + + /* Update per-channel coefficient range */ + for (b = 0; b < pcContribData -> nCompressedBins; b++) { + for (i = 0; i < 3; i++) { + if (threshCoeffs [b].coeff [i] < mrgbeRange -> min [i]) + mrgbeRange -> min [i] = threshCoeffs [b].coeff [i]; + if (threshCoeffs [b].coeff [i] > mrgbeRange -> max [i]) + mrgbeRange -> max [i] = threshCoeffs [b].coeff [i]; + } + } + + /* Init mRGBE coefficient normalisation from range */ mRGBEinit(mrgbeRange, mrgbeRange -> min, mrgbeRange -> max); - + mrgbeCoeff = pcContribData -> mrgbeCoeffs; + /* Encode wavelet detail coefficients to mRGBE */ - for (b = 0, threshCoeffs = pcContribData -> threshCoeffs, - mrgbeCoeff = pcContribData -> mrgbeCoeffs; - b < pcContribData -> nCompressedBins; b++ - ) - mrgbeCoeff [b] = mRGBEencode(threshCoeffs [b].coeff, mrgbeRange, - threshCoeffs [b].bin + for (b = 0; b < pcContribData -> nCompressedBins; b++) + mrgbeCoeff [b] = mRGBEencode(threshCoeffs [b].coeff, + mrgbeRange, threshCoeffs [b].bin ); - return mrgbeCoeff; + return pcContribData; + } + + + + static void freePreCompContribLUTentry (void *p) + /* Free data for per-modifier precomputed contribution LUT entry */ + { + PreComputedPhotonContrib *pcContribData = p; + + /* Free primary and transposed wavelet matrices */ + free(pcContribData -> waveletMatrix); + freeWaveletMatrix(pcContribData -> tWaveletMatrix, pcContribData -> l); + + /* Free thresholded coefficients */ + free(pcContribData -> threshCoeffs); + + /* Free mRGBE encoded coefficients */ + free(pcContribData -> mrgbeCoeffs); + } + + + + static void initContribHeap (PhotonMap *pmap) + /* Initialise precomputed contribution heap */ + { + int fdFlags; + + if (!pmap -> contribHeap) { + /* Open heap file for precomputed contributions */ + mktemp(strcpy(pmap -> contribHeapFname, PMAP_TMPFNAME)); + if (!(pmap -> contribHeap = fopen(pmap -> contribHeapFname, "w+b"))) + error(SYSTEM, "failed opening precomputed contribution " + "heap file in initPhotonHeap()" + ); +#ifdef F_SETFL /* XXX is there an alternative needed for Wind0z? */ + fdFlags = fcntl(fileno(pmap -> contribHeap), F_GETFL); + fcntl(fileno(pmap -> contribHeap), F_SETFL, fdFlags | O_APPEND); +#endif /* ftruncate(fileno(pmap -> heap), 0); */ + } + } + + + + static MODCONT *getPhotonContrib (PhotonMap *pmap, RAY *ray, COLOR irrad) + /* Locate photons near ray -> rop which originated from the light source + modifier indexed by ray -> rsrc, and accumulate their contributions + in pmap -> srcContrib. Return total contributions in irrad and + pointer to binned contributions (or NULL if none were found). */ + { + const OBJREC *srcMod = + findmaterial(source [ray -> rsrc].so); + MODCONT *contrib = + (MODCONT*)lu_find(pmap -> contribTab, srcMod -> oname) -> data; + Photon *photon; + COLOR flux; + PhotonSearchQueueNode *sqn; + float r2, norm; + unsigned i; + + /* Zero bins for this source modifier */ + memset(contrib -> cbin, 0, sizeof(DCOLOR) * contrib -> nbins); + setcolor(irrad, 0, 0, 0); + + if (!contrib || !pmap -> maxGather) + /* Modifier not in LUT or zero bandwidth */ + return NULL; + + /* Lookup photons */ + pmap -> squeue.tail = 0; + /* Pass light source index to filter in findPhotons() via + pmap -> lastContribSrc */ + pmap -> lastContribSrc.srcIdx = ray -> rsrc; + 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 NULL; + } + + /* Avg radius^2 between furthest two photons to improve accuracy */ + sqn = pmap -> squeue.node + 1; + r2 = max(sqn -> dist2, (sqn + 1) -> dist2); + r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2)); + + /* XXX: OMIT extra normalisation factor 1/PI for ambient calc? */ + #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 */ + norm = 2 / (PI * PI * r2); + #else + /* Normalise accumulated flux by search area PI * r^2, including + RADIANCE-specific lambertian factor PI */ + norm = 1 / (PI * PI * r2); + #endif + + /* Skip the extra photon */ + for (i = 1 ; i < pmap -> squeue.tail; i++, sqn++) { + /* Get photon's contribution to density estimate */ + photon = getNearestPhoton(&pmap -> squeue, sqn -> idx); + getPhotonFlux(photon, flux); + scalecolor(flux, norm); + #ifdef PMAP_EPANECHNIKOV + /* Apply Epanechnikov kernel to photon flux based on photon distance */ + scalecolor(flux, 1 - sqn -> dist2 / r2); + #endif + addcolor(irrad, flux); + addcolor(contrib -> cbin [photonSrcBin(pmap, photon)], flux); + } + + return contrib; } void preComputeContrib (PhotonMap *pmap) /* Precompute contributions for a random subset of finalGather * - pmap -> numPhotons photons, and build the precomputed contribution photon - map, discarding the original photons. */ + (pmap -> numPhotons) photons, and build the precomputed contribution + photon map, discarding the original photons. */ /* !!! NOTE: PRECOMPUTATION WITH OOC CURRENTLY WITHOUT CACHE !!! */ { - unsigned long i, numPreComp; - unsigned j; - PhotonIdx pIdx; - Photon photon; - RAY ray; - MODCONT *binnedContribs; - PhotonMap nuPmap; + unsigned long i, numPreComp; + unsigned j; + PhotonIdx pIdx; + Photon photon; + RAY ray; + MODCONT *binnedContribs; + PhotonMap nuPmap; + COLR coeffNorm; + PreComputedPhotonContrib *pcContribData; /* Init LUT for per-modifier wavelet matrices and mRGBE coeffs */ - LUTAB preCompContribLUT = LU_SINIT(NULL, NULL); + LUTAB preCompContribLUT = LU_SINIT(NULL, + freePreCompContribLUTentry + ); if (verbose) { sprintf(errmsg, "\nPrecomputing contributions for %ld photons\n", numPreComp ); eputs(errmsg); #if NIX fflush(stderr); #endif } /* Copy photon map for precomputed photons */ memcpy(&nuPmap, pmap, sizeof(PhotonMap)); - /* Zero counters, init new heap and extents */ + /* Zero counters, init new photon and contribution heaps */ nuPmap.numPhotons = 0; + /* ------------------------------------------- */ + /* !!!BUG: NEED SEPARATE HEAPS PER MODIFIER!!! */ + /* ------------------------------------------- */ initPhotonHeap(&nuPmap); + initContribHeap(&nuPmap); + /* Init extents of precomputed photon distribution */ for (j = 0; j < 3; j++) { nuPmap.minPos [j] = FHUGE; nuPmap.maxPos [j] = -FHUGE; } /* Record start time, baby */ repStartTime = time(NULL); #ifdef SIGCONT signal(SIGCONT, pmapPreCompReport); #endif repComplete = numPreComp = finalGather * pmap -> numPhotons; repProgress = 0; photonRay(NULL, &ray, PRIMARY, NULL); ray.ro = NULL; for (i = 0; i < numPreComp; i++) { do { /* Get random photon from stratified distribution in source heap * to avoid duplicates and clustering */ pIdx = firstPhoton(pmap) + (unsigned long)( (i + pmapRandom(pmap -> randState)) / finalGather ); getPhoton(pmap, pIdx, &photon); /* Init dummy photon ray with intersection and normal at photon position. Set emitting light source index from origin. */ VCOPY(ray.rop, photon.pos); ray.rsrc = photonSrcIdx(pmap, &photon); for (j = 0; j < 3; j++) ray.ron [j] = photon.norm [j] / 127.0; - /* Get contributions at photon position; retry with another photon - * if no contribs found */ + /* Get contributions at photon position; retry with another + photon if no contribs found */ binnedContribs = getPhotonContrib(pmap, &ray, ray.rcol); - - /* Compress & encode binned contribs */ - if (!encodeContribs(binnedContribs, &preCompContribLUT, - compressRatio - )) - error(INTERNAL, "failed compression/encoding of photon " - "contributions in preComputeContrib()" - ); - - /* TODO: Dump encoded bins to heap file */ } while (!binnedContribs); + + /* Compress & encode binned contribs */ + if (!(pcContribData = encodeContribs(binnedContribs, + &preCompContribLUT, compressRatio) + ) + ) + error(INTERNAL, "failed compression/encoding of photon " + "contributions in preComputeContrib()" + ); + + /* Set photon flux to average wavelet coefficient */ + copycolor(ray.rcol, binnedContribs -> cbin [0]); + + /* Dump encoded bins to contrib heap file, prepended by + normalisation in 32-bit RGBE */ + setcolr(coeffNorm, + pcContribData -> mrgbeRange.norm [0], + pcContribData -> mrgbeRange.norm [1], + pcContribData -> mrgbeRange.norm [2] + ); + putbinary(coeffNorm, 1, sizeof(coeffNorm), nuPmap.contribHeap); + + if (putbinary(pcContribData -> mrgbeCoeffs, sizeof(mRGBE), + pcContribData -> nCompressedBins, nuPmap.contribHeap + ) != pcContribData -> nCompressedBins + ) + error(SYSTEM, "failed writing precomputed contribution " + "heap in preComputeContrib()" + ); + /* HACK: signal newPhoton() to set precomputed photon's + contribution source from ray -> rsrc */ + nuPmap.lastContribSrc.srcIdx = -2; + /* Append photon to new heap from ray */ newPhoton(&nuPmap, &ray); /* Update progress */ repProgress++; if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) pmapPreCompReport(); #ifdef SIGCONT else signal(SIGCONT, pmapPreCompReport); #endif } - /* Flush heap */ + /* Flush photon/contrib heaps */ flushPhotonHeap(&nuPmap); + fflush(nuPmap.contribHeap); #ifdef SIGCONT signal(SIGCONT, SIG_DFL); #endif + /* Trash original pmap, replace with precomputed one */ deletePhotons(pmap); memcpy(pmap, &nuPmap, sizeof(PhotonMap)); + + /* Trash precomputed contribs LUT along with associated data */ + lu_done(&preCompContribLUT); if (verbose) { eputs("\nRebuilding precomputed contribution photon map\n"); #if NIX fflush(stderr); #endif } - /* Rebuild underlying data structure, destroying heap */ + /* Rebuild underlying data structure, destroying heap */ buildPhotonMap(pmap, NULL, NULL, 1); /* TODO: Save wavelet coeffs here? */ } - -#else +#else + /* + -------------------------------------------------------------------- + U N I T T E S T S + -------------------------------------------------------------------- + */ #include #include #include "func.h" int main (int argc, char *argv []) { unsigned i, l, nbins, nsamp, numTheta = 0, numPhi = 0; double t, p; RAY ray; int bin; if (argc < 3) { fprintf(stderr, "%s [=; ..]\n", argv [0]); fprintf(stderr, "Default variable defs: %s\n", PMAP_CONTRIB_SCDEFAULTS ); fputs("\nMissing resolution l>1, number of samples nsamp\n", stderr ); return -1; } if (!(l = atoi(argv [1]))) { fputs("Invalid resolution\n", stderr); return -1; } else nbins = PMAP_CONTRIB_SCBINS(l); if (!(nsamp = atoi(argv [2]))) { fputs("Invalid num samples\n", stderr); return -1; } else { numTheta = (int)(sqrt(nsamp) / 2); numPhi = 4* numTheta; } /* (Doan' need to) Init cal func routines for binning */ #if 0 initfunc(); setcontext(RCCONTEXT); #endif /* Compile default orientation variables for contribution binning */ scompile(PMAP_CONTRIB_SCDEFAULTS, NULL, 0); /* Compile custom orientation variabls from command line */ for (i = 3; i < argc; i++) scompile(argv [i], NULL, 0); for (i = 0; i < nsamp; i++) { #if 0 /* Random */ t = 0.5 * PI * drand48(); p = 2 * PI * drand48(); #else /* Stratified */ t = 0.5 * PI * ((i % numTheta) + drand48()) / numTheta; p = 2.0 * PI * ((i / numTheta) + drand48()) / numPhi; #endif ray.rdir [0] = sin(t) * cos(p); ray.rdir [1] = sin(t) * sin(p); ray.rdir [2] = cos(t); bin = ray2bin(&ray, nbins); printf("%.3f\t%.3f\t%.3f\t-->\t%d\n", ray.rdir [0], ray.rdir [1], ray.rdir [2], bin ); } return 0; } #endif diff --git a/pmapcontrib.h b/pmapcontrib.h index 31d9a74..3672458 100644 --- a/pmapcontrib.h +++ b/pmapcontrib.h @@ -1,145 +1,145 @@ /* RCSid $Id: pmapcontrib.h,v 2.5 2016/05/17 17:39:47 rschregle Exp $ */ /* ========================================================================= Photon map routines for precomputed light source contributions; these routines interface to mkpmap and rcontrib. 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") ========================================================================= $Id: pmapcontrib.h,v 2.5 2016/05/17 17:39:47 rschregle Exp $ */ #ifndef _PMAPCONTRIB_H #define _PMAPCONTRIB_H #include "pmapdata.h" #include "wavelet2.h" #include "mrgbe.h" #ifndef MAXPROCESS #include "rcontrib.h" #endif #ifndef MAXMODLIST /* Max number of contributing sources */ #define MAXMODLIST 1024 #endif /* Filename templates for precomputed contrib photon map: Wavelet compressed contributions (per modifier) */ #define PMAP_CONTRIB_WAVELETSUFFIX "%s-%s.wvt" /* The following variables can be specified to override the orientation of the Shirley-Chiu mapping (see also disk2square.cal). We use the built-in functions in disk2square.c for efficiency and in order to not depend on an external function file. These variables merely mimic the function file interace. RHS : right-hand-coordinate system (-1 for left-hand) rNx, rNy, rNz : surface normal Ux, Uy, Uz : up vector (defines phi = 0) */ #define PMAP_CONTRIB_SCRHS "RHS" #define PMAP_CONTRIB_SCNORMX "rNx" #define PMAP_CONTRIB_SCNORMY "rNy" #define PMAP_CONTRIB_SCNORMZ "rNz" #define PMAP_CONTRIB_SCUPX "Ux" #define PMAP_CONTRIB_SCUPY "Uy" #define PMAP_CONTRIB_SCUPZ "Uz" #define PMAP_CONTRIB_SCDEFAULTS ( \ "RHS=1; rNx=0; rNy=0; rNz=1; Ux=0; Uy=1; Uz=0;" \ ) /* Maximum Shirley-Chiu binning resolution l*/ #define PMAP_CONTRIB_MAXBINRES (mRGBE_DATABITS >> 1) /* Shirley-Chiu square dimensions and number of bins for resolution l */ #define PMAP_CONTRIB_SCDIM(l) (1 << (l)) #define PMAP_CONTRIB_SCBINS(l) (1 << ((l) << 1)) typedef struct { WAVELET_COEFF *coeff; unsigned bin; } ThreshWaveletCoeff; typedef struct { WaveletMatrix waveletMatrix, tWaveletMatrix; ThreshWaveletCoeff *threshCoeffs; mRGBERange mrgbeRange; mRGBE *mrgbeCoeffs; unsigned l, scDim, nCompressedBins; } PreComputedPhotonContrib; MODCONT *addContribModifier(LUTAB *contribTab, unsigned *numContribs, - char *mod, char *binParm, char *binVal, int binCnt + char *mod, char *binParm, int binCnt ); /* Add light source modifier mod to contribution lookup table contribsTab, and update numContribs. Return initialised contribution data for this modifier. */ void addContribModfile(LUTAB *contribTab, unsigned *numContribs, - char *modFile, char *binParm, char *binVal, int binCnt + char *modFile, char *binParm, int binCnt ); /* Add light source modifiers from file modFile to contribution lookup * table contribTab, and update numContribs. */ void initPmapContrib (LUTAB *contribTab); /* Set up photon map contributions (interface to rcmain.c) */ PhotonContribSourceIdx newPhotonContribSource (PhotonMap *pmap, const RAY *contribSrcRay, FILE *contribSrcHeap ); /* Add contribution source ray for emitted contribution photon and save * light source index and binned direction. The current contribution source * is stored in pmap -> lastContribSrc. If the previous contribution source * spawned photons (i.e. has srcIdx >= 0), it's appended to contribSrcHeap. * If contribSrcRay == NULL, the current contribution source is still * flushed, but no new source is set. Returns updated contribution source * counter pmap -> numContribSrc. */ PhotonContribSourceIdx buildContribSources (PhotonMap *pmap, FILE **contribSrcHeap, char **contribSrcHeapFname, PhotonContribSourceIdx *contribSrcOfs, unsigned numHeaps ); /* Consolidate per-subprocess contribution sources heaps into array * pmap -> contribSrc. Returns offset for contribution source index * linearisation in pmap -> numContribSrc. The heap files in * contribSrcHeap are closed on return. */ void preComputeContrib (PhotonMap *pmap); /* Precompute contributions for a random subset of finalGather * pmap -> numPhotons photons, and build the precomputed contribution photon map, discarding the original photons. */ void distribPhotonContrib (PhotonMap *pmap, LUTAB *contribTab, unsigned numContribs, unsigned numProc ); /* Emit photons with binned light source contributions, precompute their contributions and build photon map */ void getPreCompPhotonContrib (PhotonMap *pmap, RAY *ray, COLOR totalContrib ); /* Fetch and decode precomputed light source contributions from single closest precomputed contribution photon at ray -> rop, and accumulate them in pmap -> contribTab. Also returns total precomputed contribs. */ #endif diff --git a/pmapdata.c b/pmapdata.c index 30f8488..2d8ef7c 100644 --- a/pmapdata.c +++ b/pmapdata.c @@ -1,834 +1,837 @@ #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 "pmaproi.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]; + /* Init precomputed contrib photon stuff */ + pmap -> contribHeap = NULL; /* 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) - ); + 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; /* Store photon if within a region of interest (for ze Ecksperts!) */ if (!photonInROI(ray)) 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 */ + /* Set photon's subprocess index (to consolidate contribution sources + after photon distribution completes) */ 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). - * Note the contribution source bin has already been set by - * newPhotonContribSrc(). */ - photon.aux.contribSrc = pmap -> numContribSrc; -#if 0 /* NOTE: Invalid bins are now handled in distribPhotonContrib() - immediately after returning from this routine */ - if (photonSrcBin(pmap, &photon) < 0) - /* Photon has an invalid bin, and therefore implicitly zero - contributions; don't bother storing! */ - return -1; -#endif - /* Photon will be stored; set contribution source index, - thereby marking it as having spawned photon(s) */ - pmap -> lastContribSrc.srcIdx = ray -> rsrc; + case PMAP_TYPE_CONTRIB: + if (pmap -> lastContribSrc.srcIdx < -1) { + /* HACK: Contrib photon being precomputed; set contribution + source from index passed in ray */ + photon.aux.contribSrc = ray -> rsrc; + } + else { + /* HACK: Contrib photon before precomputation (i.e. in forward + pass); set contribution source from last index in contrib + source array. Note the contribution source bin has already + been set by newPhotonContribSrc(). */ + photon.aux.contribSrc = pmap -> numContribSrc; + /* Photon will be stored; set contribution source index, + * thereby marking it as having spawned photon(s) */ + 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. */ photon.aux.pathLen = photonMaxPathDist - ray -> rmax + ray -> rot; 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]; else if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) { /* 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/pmapdata.h b/pmapdata.h index 198a0a2..5f2aaef 100644 --- a/pmapdata.h +++ b/pmapdata.h @@ -1,419 +1,427 @@ /* 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 /* Source of a contribution photon. This consists of the emitting light source and binned direction. These records are only used to precompute contribution photons. They are referenced by contribution photons (see contribIdx field in struct Photon below) in a surjective mapping, since multiple photons may share the same emitting source and direction 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; /* Binned incident direction */ } PhotonContribSource; typedef uint32 PhotonPathID; typedef uint32 PhotonContribSourceIdx; #define PMAP_MAXCONTRIBSRC UINT32_MAX #define photonSrcIdx(pm, p) ((pm) -> contribSrc \ ? (pm) -> contribSrc [(p) -> aux.contribSrc].srcIdx \ : (p) -> aux.pathID\ ) #define photonSrcBin(pm, p) ( \ (pm) -> contribSrc [(p) -> aux.contribSrc].srcBin \ ) #define photonSrcMod(pm, p) findmaterial(source [photonSrcIdx(pm, p)].so) /* Multipurpose auxiliary photon attribute type */ typedef union { /* Photon's propagation distance (= path length / time of flight) for temporal light flow */ float pathLen; /* Index into contribution photon's emitting source and binned direction; see struct PhotonContribSource above */ PhotonContribSourceIdx contribSrc; /* 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 + #define PMAP_PROCBITS 7 #else - #define PMAP_PROCBITS 5 + #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) + #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 /* Define search queue and underlying data struct types */ #ifdef PMAP_OOC #include "pmapooc.h" #else #include "pmapkdt.h" #include "pmaptkdt.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) /* Bias compensation history node */ typedef struct { COLOR irrad; float weight; } PhotonBiasCompNode; /* Forward declaration */ struct PhotonMap; 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) */ /* ================================================================ * TRANSIENT PHOTON STUFF * ================================================================ */ double velocity, /* Speed of light in units of scene geometry [1/sec] */ time, /* Photons' time of flight for transient lookups */ minPathLen, maxPathLen, avgPathLen; /* Min/max/avg path length */ /* ================================================================ * CONTRIBUTION PHOTON STUFF * ================================================================ */ - PhotonContribSource *contribSrc, /* Contribution source array */ - lastContribSrc; /* Current contrib source */ - PhotonContribSourceIdx numContribSrc; /* Number of contrib sources */ - LUTAB *contribTab; /* LUT for binned contribs */ - + PhotonContribSource + *contribSrc, /* Contribution source array */ + lastContribSrc; /* Current contrib source */ + PhotonContribSourceIdx + numContribSrc; /* Number of contrib sources */ + LUTAB *contribTab; /* LUT for binned contribs */ + + FILE *contribHeap; /* Out-of-core heap containing + unsorted precomputed contribution + photon bins prior to construction + of store */ + char contribHeapFname [sizeof(PMAP_TMPFNAME)]; + /* ================================================================ * 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) */ #ifdef PMAP_PATHFILT /* ================================================================ * PHOTON PATH FILTERING STUFF * ================================================================ */ 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 } 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]) #define volumePmap (photonMaps [PMAP_TYPE_VOLUME]) #define transientPmap (photonMaps [PMAP_TYPE_TRANSIENT]) #ifdef PMAP_PHOTONFLOW /* Transient lightflow has precedence */ #define lightFlowPmap (transLightFlowPmap \ ? transLightFlowPmap : photonMaps [PMAP_TYPE_LIGHTFLOW] \ ) #define transLightFlowPmap (photonMaps [PMAP_TYPE_TRANSLIGHTFLOW]) #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) #define isVolumePmap(p) ((p) -> type == PMAP_TYPE_VOLUME) #define isTransientPmap(p) ((p) -> type == PMAP_TYPE_TRANSIENT) #ifdef PMAP_PHOTONFLOW /* lightflow also implies transient lightflow */ #define isLightFlowPmap(p) ( \ (p) -> type == PMAP_TYPE_LIGHTFLOW || isTransLightFlowPmap(p) \ ) #define isTransLightFlowPmap(p) ( \ (p) -> type == PMAP_TYPE_TRANSLIGHTFLOW \ ) #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, PhotonContribSourceIdx *contribSrcOfs, unsigned nproc ); /* Postprocess 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/pmcontrib2.c b/pmcontrib2.c index 4c7f0e6..8118327 100644 --- a/pmcontrib2.c +++ b/pmcontrib2.c @@ -1,732 +1,732 @@ #ifndef lint static const char RCSid[] = "$Id$"; #endif /* ========================================================================= - Photon map routines for precomputed light source contributions; this - module contains the main photon distribution routine as part of mkpmap. + Photon map routines for precomputed light source contributions. + This module contains the main photon distribution routine in mkpmap. 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") ========================================================================= $Id$ */ #include "pmapcontrib.h" #include "pmapdiag.h" #include "pmaprand.h" #include "pmapmat.h" #include "pmaproi.h" #include "pmapsrc.h" #include "otspecial.h" #if NIX #include #include #endif /* Defs for photon emission counter array passed by sub-processes to parent * via shared memory */ typedef unsigned long PhotonContribCnt; /* Indices for photon emission counter array: num photons stored and num * emitted per source */ #define PHOTONCNT_NUMPHOT 0 #define PHOTONCNT_NUMEMIT(n) (1 + n) /* Photon distribution counter update interval expressed as bitmask; counters shared among subling subprocesses will only be updated in multiples of PMAP_CNTUPDATE in order to reduce contention */ -#define PMAP_CNTUPDATE 0xffL +#define PMAP_CNTUPDATE 0xffL void distribPhotonContrib (PhotonMap* pmap, LUTAB *contribTab, unsigned numContribs, unsigned numProc ) { EmissionMap emap; char errmsg2 [128], shmFname [PMAP_TMPFNLEN]; unsigned srcIdx, proc; int shmFile, stat, pid; double *srcFlux, /* Emitted flux per light source */ srcDistribTarget; /* Target photon count per source */ PhotonContribCnt *photonCnt, /* Photon emission counter array */ lastPhotonCnt [PHOTONCNT_NUMEMIT(nsources)]; unsigned photonCntSize = (sizeof(PhotonContribCnt) * PHOTONCNT_NUMEMIT(nsources) ); FILE **contribSrcHeap = NULL; char **contribSrcHeapFname = NULL; PhotonContribSourceIdx *contribSrcOfs = NULL; pid_t procPids [PMAP_MAXPROC]; if (!pmap) error(USER, "no contribution photon map specified"); if (!nsources) error(USER, "no light sources"); if (nsources > MAXMODLIST) error(USER, "too many light sources"); if (!contribTab || !numContribs) error(USER, "no modifiers specified for contribution photon map"); /* Allocate photon flux per light source; this differs for every * source as all sources contribute the same number of distributed * photons (srcDistribTarget), hence the number of photons emitted per * source does not correlate with its emitted flux. The resulting flux * per photon is therefore adjusted individually for each source. */ if (!(srcFlux = calloc(nsources, sizeof(double)))) error(SYSTEM, "can't allocate source flux in distribPhotonContrib"); /* =================================================================== * INITIALISATION - Set up emission and scattering funcs * =================================================================== */ emap.samples = NULL; emap.src = NULL; emap.maxPartitions = MAXSPART; emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1); if (!emap.partitions) error(USER, "can't allocate source partitions in distribPhotonContrib"); /* Initialise contrib photon map */ initPhotonMap(pmap, PMAP_TYPE_CONTRIB); initPmapContrib(contribTab); initPhotonHeap(pmap); initPhotonEmissionFuncs(); initPhotonScatterFuncs(); /* Per-subprocess / per-source target counts */ pmap -> distribTarget /= numProc; srcDistribTarget = nsources ? (double)pmap -> distribTarget / nsources : 0; if (!pmap -> distribTarget) error(INTERNAL, "no photons to distribute in distribPhotonContrib"); /* Get photon ports from modifier list */ getPhotonPorts(photonPortList); /* Get photon sensor modifiers */ getPhotonSensors(photonSensorList); /* Get polyhedral regions of interest */ getPolyROIs(pmapROImodList); #if NIX /* Set up shared mem for photon counters (zeroed by ftruncate) */ strcpy(shmFname, PMAP_TMPFNAME); shmFile = mkstemp(shmFname); if (shmFile < 0 || ftruncate(shmFile, photonCntSize) < 0) error(SYSTEM, "failed shared mem init in distribPhotonContrib"); photonCnt = mmap(NULL, photonCntSize, PROT_READ | PROT_WRITE, MAP_SHARED, shmFile, 0 ); if (photonCnt == MAP_FAILED) - error(SYSTEM, "failed shsared mem mapping in distribPhotonContrib"); + error(SYSTEM, "failed mapping shared memory in distribPhotonContrib"); #else /* Allocate photon counters statically on Windoze */ if (!(photonCnt = malloc(photonCntSize))) error(SYSTEM, "failed trivial malloc in distribPhotonContrib"); for (srcIdx = 0; srcIdx < PHOTONCNT_NUMEMIT(nsources); srcIdx++) photonCnt [srcIdx] = 0; #endif /* NIX */ /* Zero overflow tracking counters */ memset(&lastPhotonCnt, 0, sizeof(lastPhotonCnt)); if (verbose) { sprintf(errmsg, "\nIntegrating flux from %d sources", nsources); if (photonPorts) { sprintf(errmsg2, " via %d ports", numPhotonPorts); strcat(errmsg, errmsg2); } strcat(errmsg, "\n"); eputs(errmsg); } /* ============================================================= * FLUX INTEGRATION - Get total flux emitted from sources/ports * ============================================================= */ for (srcIdx = 0; srcIdx < nsources; srcIdx++) { unsigned portCnt = 0; const OBJREC *srcMod = findmaterial(source [srcIdx].so); srcFlux [srcIdx] = 0; /* Skip this source if its contributions are not sought */ if (!lu_find(pmap -> contribTab, srcMod -> oname) -> data) { sprintf(errmsg, "ignoring contributions from source %s", source [srcIdx].so -> oname ); error(WARNING, errmsg); continue; } emap.src = source + srcIdx; do { /* Need at least one iteration if no ports! */ emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt : NULL; photonPartition [emap.src -> so -> otype] (&emap); if (verbose) { sprintf(errmsg, "\tIntegrating flux from source %s ", source [srcIdx].so -> oname ); if (emap.port) { sprintf(errmsg2, "via port %s ", photonPorts [portCnt].so -> oname ); strcat(errmsg, errmsg2); } sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions); strcat(errmsg, errmsg2); eputs(errmsg); #if NIX fflush(stderr); #endif } for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions; emap.partitionCnt++ ) { initPhotonEmission(&emap, pdfSamples); srcFlux [srcIdx] += colorAvg(emap.partFlux); } portCnt++; } while (portCnt < numPhotonPorts); if (srcFlux [srcIdx] < FTINY) { sprintf(errmsg, "source %s has zero emission", source [srcIdx].so -> oname ); error(WARNING, errmsg); } } /* Allocate & init per-subprocess contribution source heap files */ contribSrcHeap = calloc(numProc, sizeof(FILE*)); contribSrcHeapFname = calloc(numProc, sizeof(char*)); contribSrcOfs = calloc(numProc, sizeof(PhotonContribSourceIdx)); if (!contribSrcHeap || !contribSrcHeapFname || !contribSrcOfs) error(SYSTEM, "failed contribution source heap allocation " "in distribPhotonContrib()" ); for (proc = 0; proc < numProc; proc++) { contribSrcHeapFname [proc] = malloc(PMAP_TMPFNLEN); if (!contribSrcHeapFname [proc]) error(SYSTEM, "failed contribution source heap file allocation " "in distribPhotonContrib()" ); mktemp(strcpy(contribSrcHeapFname [proc], PMAP_TMPFNAME)); if (!(contribSrcHeap [proc] = fopen(contribSrcHeapFname [proc], "w+b"))) error(SYSTEM, "failed opening contribution source heap file " "in distribPhotonContrib()" ); } /* Record start time for progress reports */ repStartTime = time(NULL); if (verbose) { sprintf(errmsg, "\nPhoton distribution @ %d procs\n", numProc); eputs(errmsg); } /* MAIN LOOP */ for (proc = 0; proc < numProc; proc++) { #if NIX if (!(pid = fork())) { /* SUBPROCESS ENTERS HERE; opened and mmapped files inherited */ #else if (1) { /* No subprocess under Windoze */ #endif /* Local photon counters for this subprocess */ unsigned long lastNumPhotons = 0, localNumEmitted = 0; double photonFluxSum = 0; /* Accum. photon flux */ /* Seed RNGs from PID for decorellated photon distribution */ pmapSeed(randSeed + proc, partState); pmapSeed(randSeed + (proc + 1) % numProc, emitState); pmapSeed(randSeed + (proc + 2) % numProc, cntState); pmapSeed(randSeed + (proc + 3) % numProc, mediumState); pmapSeed(randSeed + (proc + 4) % numProc, scatterState); pmapSeed(randSeed + (proc + 5) % numProc, rouletteState); #ifdef PMAP_SIGUSR { double partNumEmit; unsigned long partEmitCnt; double srcPhotonFlux, avgPhotonFlux; unsigned portCnt, passCnt, prePassCnt; float srcPreDistrib; double srcNumEmit; /* # to emit from source */ unsigned long srcNumDistrib; /* # stored */ void sigUsrDiags() /* Loop diags via SIGUSR1 */ { sprintf(errmsg, "********************* Proc %d Diags *********************\n" "srcIdx = %d (%s)\nportCnt = %d (%s)\npassCnt = %d\n" "srcFlux = %f\nsrcPhotonFlux = %f\navgPhotonFlux = %f\n" "partNumEmit = %f\npartEmitCnt = %lu\n\n", proc, srcIdx, findmaterial(source [srcIdx].so) -> oname, portCnt, photonPorts [portCnt].so -> oname, passCnt, srcFlux [srcIdx], srcPhotonFlux, avgPhotonFlux, partNumEmit, partEmitCnt ); eputs(errmsg); fflush(stderr); } } #endif #ifdef PMAP_SIGUSR signal(SIGUSR1, sigUsrDiags); #endif #ifdef DEBUG_PMAP /* Output child process PID after random delay to prevent corrupted * console output due to race condition */ usleep(1e6 * pmapRandom(rouletteState)); fprintf(stderr, "Proc %d: PID = %d (waiting 10 sec to attach debugger...)\n", proc, getpid() ); /* Allow time for debugger to attach to child process */ sleep(10); #endif /* ============================================================= * 2-PASS PHOTON DISTRIBUTION * Pass 1 (pre): emit fraction of target photon count * Pass 2 (main): based on outcome of pass 1, estimate remaining * number of photons to emit to approximate target * count * ============================================================= */ for (srcIdx = 0; srcIdx < nsources; srcIdx++) { const unsigned numEmitIdx = PHOTONCNT_NUMEMIT(srcIdx); #ifndef PMAP_SIGUSR unsigned portCnt, passCnt = 0, prePassCnt = 0; float srcPreDistrib = preDistrib; double srcNumEmit = 0; /* # to emit from source */ unsigned long srcNumDistrib = pmap -> numPhotons; /* #stored */ #else passCnt = prePassCnt = 0; srcPreDistrib = preDistrib; srcNumEmit = 0; /* # to emit from source */ srcNumDistrib = pmap -> numPhotons; /* # stored */ #endif if (srcFlux [srcIdx] < FTINY) /* Source has zero emission or was skipped in prepass because its contributions are not sought */ continue; while (passCnt < 2) { if (!passCnt) { /* INIT PASS 1 */ if (++prePassCnt > maxPreDistrib) { /* Warn if no photons contributed after sufficient * iterations; only output from subprocess 0 to reduce * console clutter */ if (!proc) { sprintf(errmsg, "source %s: too many prepasses, skipped", source [srcIdx].so -> oname ); error(WARNING, errmsg); } break; } /* Num to emit is fraction of target count */ srcNumEmit = srcPreDistrib * srcDistribTarget; } else { /* INIT PASS 2 */ #ifndef PMAP_SIGUSR double srcPhotonFlux, avgPhotonFlux; #endif /* Based on the outcome of the predistribution we can now * figure out how many more photons we have to emit from * the current source to meet the target count, * srcDistribTarget. This value is clamped to 0 in case * the target has already been exceeded in pass 1. * srcNumEmit and srcNumDistrib is the number of photons * emitted and distributed (stored) from the current * source in pass 1, respectively. */ srcNumDistrib = pmap -> numPhotons - srcNumDistrib; srcNumEmit *= srcNumDistrib ? max(srcDistribTarget/srcNumDistrib, 1) - 1 : 0; if (!srcNumEmit) /* No photons left to distribute in main pass */ break; srcPhotonFlux = srcFlux [srcIdx] / srcNumEmit; avgPhotonFlux = photonFluxSum / (srcIdx + 1); if (avgPhotonFlux > FTINY && srcPhotonFlux / avgPhotonFlux < FTINY ) { /* Skip source if its photon flux is grossly below the * running average, indicating negligible contributions * at the expense of excessive distribution time; only * output from subproc 0 to reduce console clutter */ if (!proc) { sprintf(errmsg, "source %s: itsy bitsy photon flux, skipped", source [srcIdx].so -> oname ); error(WARNING, errmsg); } srcNumEmit = 0; /* Or just break??? */ } /* Update sum of photon flux per light source */ photonFluxSum += srcPhotonFlux; } portCnt = 0; do { /* Need at least one iteration if no ports! */ emap.src = source + srcIdx; emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt : NULL; photonPartition [emap.src -> so -> otype] (&emap); if (verbose && !proc) { /* Output from subproc 0 only to avoid race condition * on console I/O */ if (!passCnt) sprintf(errmsg, "\tPREPASS %d on source %s ", prePassCnt, source [srcIdx].so -> oname ); else sprintf(errmsg, "\tMAIN PASS on source %s ", source [srcIdx].so -> oname ); if (emap.port) { sprintf(errmsg2, "via port %s ", photonPorts [portCnt].so -> oname ); strcat(errmsg, errmsg2); } sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions ); strcat(errmsg, errmsg2); eputs(errmsg); #if NIX fflush(stderr); #endif } for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions; emap.partitionCnt++ ) { #ifndef PMAP_SIGUSR double partNumEmit; unsigned long partEmitCnt; #endif /* Get photon origin within current source partishunn * and build emission map */ photonOrigin [emap.src -> so -> otype] (&emap); initPhotonEmission(&emap, pdfSamples); /* Number of photons to emit from ziss partishunn; * scale according to its normalised contribushunn to * the emitted source flux */ partNumEmit = (srcNumEmit * colorAvg(emap.partFlux) / srcFlux [srcIdx] ); partEmitCnt = (unsigned long)partNumEmit; /* Probabilistically account for fractional photons */ if (pmapRandom(cntState) < partNumEmit - partEmitCnt) partEmitCnt++; /* Integer counter avoids FP rounding errors during * iteration */ while (partEmitCnt--) { RAY photonRay; /* Emit photon according to PDF (if any). If * accepted, allocate associated contribution * origin, and trace through scene until * absorbed/leaked; emitPhoton() sets the emitting * light source index in photonRay */ /* NOTE: rejection sampling skips incrementing the * emission counter (see below), thus compensating * for the rejected photons by increasing the photon * flux in proportion to the lower effective * emission count. * BUG: THIS INTERFERES WITH THE PROGRESS COUNTER * REPORTED TO THE PARENT, AND WITH THE * PREDISTRIBUTION PASS --> PHOTON DISTRIBUTION WILL * FINISH EARLY, WITH FEWER PHOTONS THAN TARGETTED! */ if (!emitPhoton(&emap, &photonRay)) continue; newPhotonContribSource(pmap, &photonRay, contribSrcHeap [proc] ); /* Skip photon if it has an invalid bin. It will implicitly contribute zero flux, so don't bother tracing and storing it. However, it counts as emitted (see enclosing while() loop) to avoid bias. */ if (pmap -> lastContribSrc.srcBin < 0) continue; /* Set subprocess index in photonRay for post- * distrib contribution source index linearisation; * this is propagated with the contrib source index * in photonRay and set for photon hits by * newPhoton() */ PMAP_SETRAYPROC(&photonRay, proc); tracePhoton(&photonRay); /* Update local emission counter */ localNumEmitted++; if (!(localNumEmitted & PMAP_CNTUPDATE)) { /* Update global counters shared with siblings in intervals to reduce overhead/contention */ lastPhotonCnt [numEmitIdx] = photonCnt [numEmitIdx]; photonCnt [numEmitIdx] += PMAP_CNTUPDATE; /* Check for overflow using previous values */ if (photonCnt [numEmitIdx] < lastPhotonCnt [numEmitIdx] ) { sprintf(errmsg, "photon emission counter " "overflow (was: %ld, is: %ld)", lastPhotonCnt [numEmitIdx], photonCnt [numEmitIdx] ); error(INTERNAL, errmsg); } /* Differentially increment photon counter */ - lastNumPhotons = pmap -> numPhotons; + lastNumPhotons = photonCnt [PHOTONCNT_NUMPHOT]; photonCnt [PHOTONCNT_NUMPHOT] += pmap -> numPhotons - lastNumPhotons; /* Check for photon counter overflow (this could only happen before an emission counter overflow if the scene has an absurdly high albedo and/or very dense geometry) */ if (photonCnt [PHOTONCNT_NUMPHOT] < lastNumPhotons ) { - sprintf(errmsg, "photon counter overlow " + sprintf(errmsg, "photon counter overflow " "(was: %ld, is: %ld)", lastNumPhotons, photonCnt [PHOTONCNT_NUMPHOT] ); error(INTERNAL, errmsg); } } } #if !NIX /* Synchronous progress report on Windoze */ if (!proc && photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime ) { unsigned s; repComplete = pmap -> distribTarget * numProc; repProgress = photonCnt [PHOTONCNT_NUMPHOT]; for (repEmitted = 0, s = 0; s < nsources; s++) repEmitted += photonCnt [numEmitIdx]; pmapDistribReport(); } #endif } portCnt++; } while (portCnt < numPhotonPorts); if (pmap -> numPhotons == srcNumDistrib) { /* Double predistrib factor in case no photons were stored * for this source and redo pass 1 */ srcPreDistrib *= 2; } else { /* Now do pass 2 */ passCnt++; } } } /* Flush heap buffa one final time to prevent data corruption */ flushPhotonHeap(pmap); /* Flush last contribution origin to origin heap file */ newPhotonContribSource(pmap, NULL, contribSrcHeap [proc]); /* Heap files closed automatically on exit fclose(pmap -> heap); fclose(orgHeap [proc]); */ #ifdef DEBUG_PMAP sprintf( errmsg, "Proc %d total %ld photons\n", proc, pmap -> numPhotons ); eputs(errmsg); fflush(stderr); #endif #ifdef PMAP_SIGUSR signal(SIGUSR1, SIG_DFL); #endif #if NIX /* Terminate subprocess */ exit(0); #endif } else { /* PARENT PROCESS CONTINUES LOOP ITERATION HERE */ if (pid < 0) error(SYSTEM, "failed to fork subprocess in distribPhotonContrib()" ); else /* Saves child process IDs */ procPids [proc] = pid; } } #if NIX /* PARENT PROCESS CONTINUES HERE */ #ifdef SIGCONT /* Enable progress report signal handler */ signal(SIGCONT, pmapDistribReport); #endif /* Wait for subprocesses to complete while reporting progress */ proc = numProc; while (proc) { while (waitpid(-1, &stat, WNOHANG) > 0) { /* Subprocess exited; check status */ if (!WIFEXITED(stat) || WEXITSTATUS(stat)) { /* Exited with error; terminate siblings, clean up & bail out */ for (proc = 0; proc < numProc; proc++) kill(procPids [proc], SIGKILL); /* Unmap shared memory, squish mapped file */ munmap(photonCnt, sizeof(*photonCnt)); close(shmFile); unlink(shmFname); error(USER, "failed photon distribution"); } --proc; } /* Nod off for a bit and update progress */ sleep(1); /* Asynchronous progress report from shared subprocess counters */ repComplete = pmap -> distribTarget * numProc; repProgress = photonCnt [PHOTONCNT_NUMPHOT]; for (repEmitted = 0, srcIdx = 0; srcIdx < nsources; srcIdx++) repEmitted += photonCnt [PHOTONCNT_NUMEMIT(srcIdx)]; /* Get global photon count from shmem updated by subprocs */ pmap -> numPhotons = photonCnt [PHOTONCNT_NUMPHOT]; if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) pmapDistribReport(); #ifdef SIGCONT else signal(SIGCONT, pmapDistribReport); #endif } #endif /* NIX */ /* ================================================================ * POST-DISTRIBUTION - Set photon flux and build kd-tree, etc. * ================================================================ */ #ifdef SIGCONT /* Reset signal handler */ signal(SIGCONT, SIG_DFL); #endif free(emap.samples); if (!pmap -> numPhotons) error(USER, "empty contribution photon map"); /* Load per-subprocess contribution sources into pmap -> contribSrc */ /* Dumb compilers apparently need the char** cast */ pmap -> numContribSrc = buildContribSources(pmap, contribSrcHeap, (char**)contribSrcHeapFname, contribSrcOfs, numProc ); if (!pmap -> numContribSrc) error(INTERNAL, "no contribution sources in contribution photon map"); /* Set photon flux per source */ /* TODO: HOW/WHERE DO WE HANDLE COEFFICIENT MODE??? */ for (srcIdx = 0; srcIdx < nsources; srcIdx++) srcFlux [srcIdx] /= photonCnt [PHOTONCNT_NUMEMIT(srcIdx)]; #if NIX /* Photon counters no longer needed, unmap shared memory */ munmap(photonCnt, sizeof(*photonCnt)); close(shmFile); unlink(shmFname); #else free(photonCnt); #endif if (verbose) { eputs("\nBuilding contribution photon map...\n"); #if NIX fflush(stderr); #endif } /* Build underlying data structure; heap is destroyed */ buildPhotonMap(pmap, srcFlux, contribSrcOfs, numProc); /* Precompute binned photon contributions */ if (verbose) eputs("\n"); preComputeContrib(contribPmap); /* Free per-subprocess origin heap files */ for (proc = 0; proc < numProc; proc++) free(contribSrcHeapFname [proc]); free(contribSrcHeapFname); free(contribSrcHeap); free(contribSrcOfs); if (verbose) eputs("\n"); } diff --git a/pmcontrib3.c b/pmcontrib3.c index 220c0d3..d22d699 100644 --- a/pmcontrib3.c +++ b/pmcontrib3.c @@ -1,190 +1,190 @@ #ifndef lint static const char RCSid[] = "$Id: pmcontrib2.c,v 2.5 2018/11/08 00:54:07 greg Exp $"; #endif /* ========================================================================= - Photon map routines for precomputed light source contributions; these - routines interface to rcontrib. + Photon map routines for precomputed light source contributions. + These routines interface to rcontrib. 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") ========================================================================= $Id$ */ #include "pmapcontrib.h" #include "pmapmat.h" #include "pmapsrc.h" #include "pmaprand.h" #include "pmapio.h" #include "pmapdiag.h" #include "otypes.h" #include "otspecial.h" static void setPmapContribParams (PhotonMap *pmap, LUTAB *contribTab) /* Set parameters for light source contributions */ { /* Set light source modifier list and appropriate callback to extract * their contributions from the photon map */ if (pmap) { pmap -> contribTab = contribTab; pmap -> lookup = getPreCompPhotonContrib; /* XXX: Need tighter tolerance for contribution photon lookups? */ pmap -> gatherTolerance = 1.0; } } static int checkSourceContrib (const LUENT *lutEntry, void *p) /* Check lookup table entry for light source contributions; Returns 0 on success. */ { char *modName = (char *)((MODCONT *)lutEntry -> data) -> modname; OBJECT modObj = modifier(modName); if (modObj == OVOID) { sprintf(errmsg, "invalid modifier %s", modName); error(USER, errmsg); } else if (!islight(objptr(modObj) -> otype)) { sprintf(errmsg, "%s is not a light source modifier", modName); error(USER, errmsg); } else return 0; } void initPmapContrib (LUTAB *srcContrib) { unsigned t; /* XXX: Redundant check? */ for (t = 0; t < NUM_PMAP_TYPES; t++) if (photonMaps [t] && t != PMAP_TYPE_CONTRIB) { sprintf( errmsg, "%s photon map does not support contributions", pmapName [t] ); error(USER, errmsg); } /* Check for valid modifiers in lookup table; this doesn't return on error. */ lu_doall(srcContrib, checkSourceContrib, NULL); /* Get params */ setPmapContribParams(contribPmap, srcContrib); #if 0 if (contribPhotonMapping) { if (contribPmap -> maxGather < numSrcContrib) { /* Warn if density estimate bandwidth is lower than modifier * count, since some contributions will be missing */ error(WARNING, "photon density estimate bandwidth too low," " contributions may be underestimated"); /* Sanity check */ checkPmapContribs(contribPmap, srcContrib); } #endif } void getPreCompPhotonContrib (PhotonMap *pmap, RAY *ray, COLOR totalContrib) /* Fetch and decode precomputed light source contributions from single closest precomputed contribution photon at ray -> rop, and accumulate them in pmap -> srcContrib. Also returns total precomputed contributions. */ { unsigned b; RREAL rayCoeff [3]; Photon photon; /* Get cumulative path coefficient up to photon lookup point */ raycontrib(rayCoeff, ray, PRIMARY); setcolor(totalContrib, 0, 0, 0); /* Ignore sources */ if (ray -> ro && islight(objptr(ray -> ro -> omod) -> otype)) return; if (find1Photon(pmap, ray, &photon)) { /* Get average wavelet coeff as total contributions */ getPhotonFlux(&photon, totalContrib); if (pmap -> contribTab) { /* TODO: fetch decoded binned contribs here */ const OBJREC *mod = photonSrcMod(pmap, &photon); MODCONT *contrib = (MODCONT*)lu_find( pmap -> contribTab, mod -> oname ) -> data; COLOR binContrib; if (!contrib) { /* Contribs from this light source not sought; this indicates * the photon map was precomputed with different rcontrib * parameters */ sprintf(errmsg, "missing LUT entry for light source modifier %s; " "mismatched contribution photon map parameters?", mod -> oname ); error(CONSISTENCY, errmsg); } #if 0 /* TODO */ decodeBins(&photon, &srcBins); if (srcBins -> numBins != contrib -> nbins) { /* Number of precomputed and allocated bins don't match */ sprintf( errmsg, "wrong number of bins for light source modifier %s; " "mismatched contribution photon map parameters?", srcMod -> oname ); error(CONSISTENCY, errmsg); } for (b = 0; b <= srcBins -> numBins; b++) { /* Accumulate contribs from decoded bin */ copycolor(binContrib, srcBins [b].contrib); if (!contrib) { /* Ray coefficient mode; normalise by light source radiance * after applying pattern */ int i; raytexture(ray, srcMod -> omod); setcolor( ray -> rcol, srcMod -> oargs.farg [0], srcMod -> oargs.farg [1], srcMod -> oargs.farg [2] ); multcolor(ray -> rcol, ray -> pcol); for (i = 0; i < 3; i++) binContrib [i] = ray -> rcol [i] ? binContrib [i] / ray -> rcol [i] : 0; } multcolor(contrib, rayCoeff); addcolor(srcContrib -> cbin [b], binContrib); } #endif } } }