diff --git a/mkpmap.c b/mkpmap.c index 2315b95..7879ebd 100644 --- a/mkpmap.c +++ b/mkpmap.c @@ -1,1001 +1,1001 @@ #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 #ifndef PMAP_OOC /* Transient photon map only supported by kd-tree */ #define PMAP_TRANSIENT #endif 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 #ifdef PMAP_CONTRIB /* 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); unsigned numContribs = 0; int contribMode = 0; #endif /* 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_CONTRIB printf("-am %.1f\t\t\t\t# max photon search radius\n", maxDistFix); #endif #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"); #ifdef PMAP_CONTRIB puts("-apC file nPhotons bwidth compression" "\t# precomputed contribution photon map" ); #endif 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"); #ifdef PMAP_TRANSIENT puts("-apt file nPhotons velocity\t\t# transient photon map"); #ifdef PMAP_PHOTONFLOW puts("-apT file nPhotons velocity\t\t# transient light flow photon map"); #endif #endif #ifdef PMAP_PHOTONFLOW 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"); #ifdef PMAP_CONTRIB puts("-bn\t\t\t\t\t# number of contribution bins"); #endif 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); #ifdef PMAP_CONTRIB puts("-e expr\t\t\t\t# init expression for contrib binning"); #endif 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# limit photon distance\n", photonMaxPathDist); printf("-lr %ld\t\t\t\t# limit photon bounces\n", photonMaxBounce); #endif #ifdef PMAP_CONTRIB puts("-m mod\t\t\t\t# contribution modifier"); #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; stuff the wInD00zE Weenies! */ printf("-n %d\t\t\t\t\t# number of parallel processes\n", nproc); #endif #ifdef PMAP_CONTRIB puts("-p\t\t\t\t\t# per-modifier contrib binning params"); #endif printf("-t %-9d\t\t\t\t# time between reports\n", photonRepTime); printf(verbose ? "-v+\t\t\t\t\t# verbose console output\n" : "-v-\t\t\t\t\t# terse console output\n" ); #ifdef PMAP_CONTRIB /* Toggle to precompute contributions or coefficients */ printf("-V%c\t\t\t\t\t# precompute %s\n", contribMode ? '+' : '-', contribMode ? "contributions" : "coefficients" ); #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, nBins = 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; #ifdef PMAP_CONTRIB case 'm': /* Fixed max photon search radius for precomputed * contribution photons */ check(3, "f"); if ((maxDistFix = atof(argv [++i])) <= 0) error(USER, "invalid max photon search radius"); break; #endif 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; #ifdef PMAP_CONTRIB case 'C': /* Precomputed contribution photon map */ checkNumPhotons(i + 2); check(4, "ssif"); contribPmapParams.fileName = argv [++i]; contribPmapParams.distribTarget = parseMultiplier(argv [++i]); if (!contribPmapParams.distribTarget) goto badopt; contribPmapParams.minGather = contribPmapParams.maxGather = atoi(argv [++i]); if (!contribPmapParams.maxGather) goto badopt; compressRatio = atof(argv [++i]); if (compressRatio < 0 || compressRatio > 1 ) error(USER, "contribution photon compression ratio must " "be in range [0, 1]" ); break; #endif #ifdef PMAP_TRANSIENT 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; #ifndef PMAP_TRANSIENT 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; #ifdef PMAP_CONTRIB case 'b': switch (argv [i][2]) { case 'v': /* Back face visibility */ check_bool(3, backvis); break; case 'n': /* Number of bins for precomp contrib pmap */ check(3, "i"); nBins = atoi(argv [++i]); /* Round to nearest integer square */ nBins = sqrt(nBins) + 0.5; nBins *= nBins; - if (nBins <= 1 || nBins > PMAP_CONTRIB_MAXBINS) + if (!nBins) goto badopt; break; default: goto badopt; } break; #endif 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: #ifdef PMAP_CONTRIB /* Functional expression for precomputed contrib. pmap, e.g. to set constants for disk2square.cal */ check(2, "s"); scompile(argv [++i], NULL, 0); #else goto badopt; #endif } 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: #ifdef PMAP_CONTRIB /* Modifier name for precomputed contrib. photon map */ check(2, "s"); addContribModifier(&contribTab, &numContribs, argv [++i], binParm, nBins ); break; #else goto badopt; #endif } break; #ifdef PMAP_CONTRIB case 'M': /* Modifier file for precomputed contribution photon map */ check(2, "s"); addContribModfile(&contribTab, &numContribs, argv [++i], binParm, nBins ); break; #endif #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 #ifdef PMAP_CONTIRB case 'p': /* Parameter setting(s) for precomputed contrib. photon map */ check(2, "s"); set_eparams(binParm = argv [++i]); break; #endif case 't': /* Timer */ check(2, "i"); photonRepTime = atoi(argv [++i]); break; case 'v': /* Verbosity */ check_bool(2, verbose); break; #ifdef PMAP_CONTRIB case 'V': /* Precomputed contributions/coefficients */ check_bool(2, contribMode); break; #endif #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 */ #ifdef PMAP_CONTRIB if (contribPmap) /* Just build contrib pmap, ignore others */ distribPhotonContrib(contribPmap, &contribTab, numContribs, &contribMode, nproc ); else #endif distribPhotons(photonMaps, nproc); /* 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/mrgbe.c b/mrgbe.c index 08afa73..4740b5d 100644 --- a/mrgbe.c +++ b/mrgbe.c @@ -1,423 +1,423 @@ /* ========================================================================= mRGBE (miniRGBE): reduced RGBE encoding of normalised RGB floats plus associated integer payload data. See mrgbe.h for encoding details. Compile with -DmRGBE_TEST to enable unit test (requires libm) Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") ========================================================================= $Id$ */ #include "mrgbe.h" #include #include #include #define MIN(a,b) ((a) < (b) ? (a) : (b)) #define MAX(a,b) ((a) > (b) ? (a) : (b)) #define ABS(a) ((a) >= 0 ? (a) : -(a)) mRGBERange *mRGBEinit (mRGBERange *range, double rgbMin [3], double rgbMax [3] ) /* Initialise range state for mRGBEencode()/mRGBEdecode() with the specified * float input ranges [rgbMin, rgbMax] in _absolute_ values per colour * channel. This sets the offset and normalisation for the encoding. */ { unsigned i; double orgMin; for (i = 0; i < 3; i++) { if (rgbMin [i] >= rgbMax [i]) { fprintf(stderr, "ERROR - invalid mRGBE range\n"); return NULL; } range -> min [i] = orgMin = fabs(rgbMin [i]); range -> max [i] = fabs(rgbMax [i]); /* Iteratively reduce minimum and set normalisation factor until * normalising the original minimum exceeds mRGBE_MIN (and therefore * won't be clamped to zero). This avoids sign reversals when * encoding negative values close or equal to the original minimum. */ do { range -> min [i] = MAX(0, range -> min [i] - mRGBE_MIN); range -> norm [i] = mRGBE_RANGE / (range -> max [i] - range -> min [i]); } while (range -> min [i] > 0 && (orgMin - range -> min [i]) * range -> norm [i] <= mRGBE_MIN ); if (range -> norm [i] < mRGBE_MIN) { fprintf(stderr, "ERROR - Zero mRGBE normalisation / range overflow\n" ); return NULL; } } return range; } mRGBE mRGBEencode (double rgb [3], const mRGBERange *range, unsigned data) /* Encode signed float RGB tuple to mRGBE along with payload data. The * supplied range state must be previously initialised with mRGBEinit(). * A warning is issued if rgb lies outside the initialised range. */ { #define SGN(a, b) ((a) >= 0 ? (b) : -(b)) mRGBE mrgbe = {0, 0, 0, 0}; double rgbRange [3], rgbNorm [3], mantNorm; int i, exp; /* Check payload data */ if (data > mRGBE_DATAMAX) { - fprintf(stderr, "ERROR - mRGBE data overflow"); + fprintf(stderr, "ERROR - mRGBE data overflow\n"); return mrgbe; } if (!range) { fprintf(stderr, "ERROR - mRGBE range not initialised\n"); return mrgbe; } /* Offset and normalise RGB according to initialised range; note rgbNorm * is in ABSOLUTE VALUES */ for (i = 0; i < 3; i++) { rgbNorm [i] = fabs(rgb [i]); /* Check range */ if (rgbNorm [i] < range -> min [i] || rgbNorm [i] > range -> max [i]) { fprintf(stderr, "ERROR - mRGBE input outside range\n"); return mrgbe; } rgbNorm [i] = (rgbNorm [i] - range -> min [i]) * range -> norm [i]; if (fabs(rgbNorm [i] - mRGBE_MAX) < mRGBE_MIN) /* Normalisation to maximum. This must be avoided by offsetting by minimum, otherwise frexp() returns a positive exponent. Normalisation to >mRGBE_MAX is caught by overflow check below */ rgbNorm [i] = mRGBE_RANGE; } /* Mantissa normalisation factor (= RGB maximum) */ mantNorm = MAX(MAX(rgbNorm [0], rgbNorm [1]), rgbNorm [2]); /* Special handling for zero */ if (mantNorm < mRGBE_MIN) mantNorm = exp = 0; else /* Get normalised mantissa and exponent to base 2 */ mantNorm = frexp(mantNorm, &exp) * mRGBE_MANTMAX / mantNorm; /* Set RGB integer mantissas, adding excess (zero offset) and rounding up or down, depending on sign. Set sign according to original value. Clamp to mRGBE_MANTOVF to prevent rounding beyond overflow into negative! */ mrgbe.red = MIN((int)( SGN(rgb [0], rgbNorm [0] * mantNorm + mRGBE_ROUND) + mRGBE_MANTMAX ), mRGBE_MANTOVF ); mrgbe.grn = MIN((int)( SGN(rgb [1], rgbNorm [1] * mantNorm + mRGBE_ROUND) + mRGBE_MANTMAX ), mRGBE_MANTOVF ); mrgbe.blu = MIN((int)( SGN(rgb [2], rgbNorm [2] * mantNorm + mRGBE_ROUND) + mRGBE_MANTMAX ), mRGBE_MANTOVF ); /* Drop the exponent sign, since implicitly negative */ mrgbe.exp = abs(exp); /* Set payload data */ mrgbe.dat = data; return mrgbe; #undef SGN } unsigned mRGBEdecode (mRGBE mrgbe, const mRGBERange *range, double rgb [3]) /* Decode mRGBE, returning signed float RGB tuple in array rgb. The * associated payload data is decoded as return value. The supplied range * state must be previously initialised with mRGBEinit(). */ { #define SGN(a, b) ((a) >= mRGBE_MANTMAX ? (b) : -(b)) #ifdef mRGBE_ZEROJITTER #define JITTER(m, e) (mRGBE_JITTER * drand48()) #else #define JITTER(m, e) ((m) == mRGBE_MANTMAX && !(e) \ ? 0 : mRGBE_JITTER * drand48() \ ) #endif double mantNorm; if (!range) { fprintf(stderr, "ERROR - mRGBE range not initialised\n"); return 0; } /* Mantissa normalisation factor; exponent is implicitly negative */ mantNorm = ldexp(1.0, -mrgbe.exp) / mRGBE_MANTMAX; /* Decode and set RGB components, subtracting excess (zero offset) and * adding jitter to break up quantisation artefacts. Clamp to -mRGBE_MAX * in case jitter underflows. (Note jitter is added/subtracted for * negative/positive values to compensate rounding down/up during * encoding). Finally denormalise and offset to original range with * saved state. */ rgb [0] = SGN(mrgbe.red, 1) * ( ABS(mrgbe.red - mRGBE_MANTMAX + JITTER(mrgbe.red, mrgbe.exp)) * mantNorm / range -> norm [0] + range -> min [0] ); rgb [1] = SGN(mrgbe.grn, 1) * ( ABS(mrgbe.grn - mRGBE_MANTMAX + JITTER(mrgbe.grn, mrgbe.exp)) * mantNorm / range -> norm [1] + range -> min [1] ); rgb [2] = SGN(mrgbe.blu, 1) * ( ABS(mrgbe.blu - mRGBE_MANTMAX + JITTER(mrgbe.blu, mrgbe.exp)) * mantNorm / range -> norm [2] + range -> min [2] ); /* Return payload data */ return mrgbe.dat; #undef SGN #undef JITTER } #ifdef mRGBE_TEST /* Build standalone to run unit test */ #define RGBMIN 0 #define RGBMAX 1 int main (int argc, char **argv) { double rgb [3] = {0, 0, 0}, rgbDec [3], rgbMin [3] = {0, 0, 0}, rgbMax [3] = {1, 1, 1}, rgbAvg [3] = {0, 0, 0}, rgbDecAvg [3] = {0, 0, 0}, rgbAbs, relDiff, rmse, rmseSum; unsigned numTests, i, j, k, data, dataDec; const short rndState0 [3] = {11, 23, 31}; short encFlag = 0, rndState [3]; mRGBERange mrgbeRange; mRGBE mrgbe; #if 0 rgb [0] = -1.481; rgb [1] = -0.437; rgb [2] = -1.392; rgbMin [0] = 1.003; rgbMin [1] = 0.436; rgbMin [2] = 0.299; rgbMax [0] = 14.863; rgbMax [1] = 13.505; rgbMax [2] = 13.359; mRGBEinit(&mrgbeRange, rgbMin, rgbMax); mrgbe = mRGBEencode(rgb, &mrgbeRange, 0); printf("% 7.4f\t% 7.4f\t% 7.4f\t%4d\t-->\t" "%02d:%02d:%02d:%02d\t-->\t", rgb [0], rgb [1], rgb [2], 0, mrgbe.red, mrgbe.grn, mrgbe.blu, mrgbe.exp ); /* Decode */ dataDec = mRGBEdecode(mrgbe, &mrgbeRange, rgbDec); for (j = rmse = 0; j < 3; j++) { relDiff = fabs(rgb [j]) > 0 ? 1 - rgbDec [j] / rgb [j] : 0; rmse += relDiff * relDiff; } rmseSum += rmse = sqrt(rmse / 3); printf("% 7.4f\t% 7.4f\t% 7.4f\t%4d\tRMSE = %.1f%%\n", rgbDec [0], rgbDec [1], rgbDec [2], dataDec, 100 * rmse ); return 0; #endif if (argc > 1 && (numTests = atoi(argv [1])) > 0) { /* Test range checks */ mRGBEinit(&mrgbeRange, rgbMin, rgbMax); puts("Check zero"); rgb [0] = rgb [1] = rgb [2] = 0; mrgbe = mRGBEencode(rgb, &mrgbeRange, 0); mRGBEdecode(mrgbe, &mrgbeRange, rgbDec); printf("% 7.3f\t% 7.3f\t% 7.3f\t-->\t%02d:%02d:%02d:%02d\t-->" "\t% 7.3f\t% 7.3f\t%7.3f\n", rgb [0], rgb [1], rgb [2], mrgbe.red, mrgbe.grn, mrgbe.blu, mrgbe.exp, rgbDec [0], rgbDec [1], rgbDec [2] ); puts("Check data overflow"); rgb [0] = rgb [1] = rgb [2] = mRGBE_MAX / 2; mRGBEencode(rgb, &mrgbeRange, mRGBE_DATAMAX << 1); puts("\nCheck underflow"); rgb [0] = rgb [1] = rgb [2] = mRGBE_MIN / 2; mRGBEencode(rgb, &mrgbeRange, 0); puts("\nCheck overflow (positive)"); rgb [0] = rgb [1] = rgb [2] = 2 * mRGBE_MAX; mRGBEencode(rgb, &mrgbeRange, 0); puts("\nCheck overflow (negative)"); rgb [0] = rgb [1] = rgb [2] = -2 * mRGBE_MAX; mRGBEencode(rgb, &mrgbeRange, 0); puts("\nCheck empty range"); rgbMax [0] = rgbMax [1] = rgbMax [2] = 0; mRGBEinit(&mrgbeRange, rgbMin, rgbMax); putchar('\n'); rgbMin [0] = rgbMin [1] = rgbMin [2] = RGBMAX; rgbMax [0] = rgbMax [1] = rgbMax [2] = RGBMIN; /* Iterate twice with same random RGB tuples; gather minimum and * maximum in 1st pass, then initialise and run encode/decode tests * in 2nd pass. */ do { /* Initialise RNG state so this can be restored to obtain the same RGB tuple sequence in both passes. This state needs to be kept separate from the RNG used for jittering in the encode/decode routines. */ for (i = 0; i < 3; i++) rndState [i] = rndState0 [i]; /* Random encode vs. decode */ for (k = 0, rmseSum = 0; k < numTests; k++) { data = (unsigned)(mRGBE_DATAMAX * erand48(rndState)); for (i = 0; i < 3; i++) { rgb [i] = (erand48(rndState) > 0.5 ? 1 : -1) * ((RGBMAX - RGBMIN) * erand48(rndState) + RGBMIN); #if 1 /* Correlate RGB */ if (i > 0) rgb [i] = rgb [0] * (0.1 + 0.9 * erand48(rndState)); #endif #if 0 /* Check zeros */ rgb [i] = 0; #endif } if (encFlag) { printf("% 4d:\t", k); /* Encode */ mrgbe = mRGBEencode(rgb, &mrgbeRange, data); printf("% 7.3f\t% 7.3f\t% 7.3f\t%4d\t-->\t" "%02d:%02d:%02d:%02d\t-->\t", rgb [0], rgb [1], rgb [2], data, mrgbe.red, mrgbe.grn, mrgbe.blu, mrgbe.exp ); /* Decode */ dataDec = mRGBEdecode(mrgbe, &mrgbeRange, rgbDec); for (j = rmse = 0; j < 3; j++) { relDiff = fabs(rgb [j]) > 0 ? 1 - rgbDec [j] / rgb [j] : 0; rmse += relDiff * relDiff; /* Update original and decoded averages */ rgbAvg [j] += rgb [j]; rgbDecAvg [j] += rgbDec [j]; } rmseSum += rmse = sqrt(rmse / 3); printf("% 7.3f\t% 7.3f\t% 7.3f\t%4d\tRMSE = %.1f%% %s\n", rgbDec [0], rgbDec [1], rgbDec [2], dataDec, 100 * rmse, rmse > 0.1 ? "!" : "" ); #if 1 /* Stop at sign reversal; this error isn't acceptable * under any circumstances! */ for (j = 0; j < 3; j++) if (rgb [j] / rgbDec [j] < 0) { puts("**** SIGN REVERSAL!!! ****\n"); return -1; } #endif } else { /* Update RGB minimum and maximum */ for (i = 0; i < 3; i++) { rgbAbs = fabs(rgb [i]); if (rgbAbs > rgbMax [i]) rgbMax [i] = rgbAbs; if (rgbAbs < rgbMin [i]) rgbMin [i] = rgbAbs; #if 0 /* NOTE: Clamping range minimum to 0 increases RMSE if RGBMIN > 0 */ rgbMin [i] = 0; #endif } } } /* Dump results */ if (!encFlag) { /* Initialise mRGBE encoding state with min/max */ if (!mRGBEinit(&mrgbeRange, rgbMin, rgbMax)) return -1; printf( "RGB Range = [%.3f, %.3f, %.3f] - [%.3f, %.3f, %.3f]\n\n", rgbMin [0], rgbMin [1], rgbMin [2], rgbMax [0], rgbMax [1], rgbMax [2] ); } else { putchar(' '); for (i = 0; i < 118; i++) putchar('-'); printf("\n Avg:\t% 7.3f\t% 7.3f\t% 7.3f\t\t\t\t\t\t" "% 7.3f\t% 7.3f\t% 7.3f\n", rgbAvg [0] / numTests, rgbAvg [1] / numTests, rgbAvg [2] / numTests, rgbDecAvg [0] / numTests, rgbDecAvg [1] / numTests, rgbDecAvg [2] / numTests ); printf("\nAvg RMSE = %.2f%%\n", 100 * rmseSum / numTests); } } while (encFlag = !encFlag); return 0; } else { printf("Usage: %s \n", argv [0]); return -1; } } #endif diff --git a/pmapcontrib.c b/pmapcontrib.c index 48eb083..32df644 100644 --- a/pmapcontrib.c +++ b/pmapcontrib.c @@ -1,1153 +1,1165 @@ #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 handle contribution binning, compression and encoding, and are used by 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 "pmutil.h" #include "otspecial.h" #include "otypes.h" #include "lookup.h" #ifdef PMAP_CONTRIB /* The following are convenient placeholders interfacing to mkpmap and rcontrib. They can be externally set via initPmapContribTab() and referenced within the contrib pmap modules. These variables can then be safely ignored by rtrace/rpict/rvu, without annoying linking errors. */ /* Global pointer to rcontrib's contribution binning LUT */ LUTAB *pmapContribTab = NULL; /* Contribution/coefficient mode flag */ int *pmapContribMode; extern void SDdisk2square(double sq[2], double diskx, double disky); #if 0 int xy2lin (unsigned scDim, int x, int y) /* Serialise 2D contribution coords in (scDim x scDim) Shirley-Chiu square to 1D index. Returns -1 if coordinates invalid */ { return x < 0 || y < 0 ? -1 : (x * scDim) + y; } void lin2xy (unsigned scDim, int bin, int *x, int *y) /* Deserialise 1D contribution index to 2D coordinates in (scDim x scDim) * Shirley-Chiu square. Returns -1 if bin invalid */ { *x = bin < 0 ? -1 : bin / scDim; *y = bin < 0 ? -1 : bin % scDim; } #endif static int ray2bin (const RAY *ray, unsigned scDim) /* Map ray dir (pointing away from origin) to its 1D bin in an (scDim x scDim) Shirley-Chiu square. Returns -1 if mapped bin is invalid (e.g. behind plane) */ { 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 PMAP_CONTRIB_XY2LIN(scDim, scBin [0], scBin [1]); } else return -1; } /* ------------------ CONTRIBSRC STUFF --------------------- */ #ifndef PMAP_CONTRIB_TEST MODCONT *addContribModifier (LUTAB *contribTab, unsigned *numContribs, 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; - - /* Check for bin overflow in mRGBE encoding BEFORE photon distrib! */ - if (binCnt > PMAP_CONTRIB_MAXBINS) { - sprintf(errmsg, "too many bins for modifier %s", mod); + unsigned numCoeffs; + + /* Check for overflow in mRGBE encoding BEFORE photon distrib! + This requires getting the number of wavelet coefficients generated + by the transform a priori. */ + if (!(numCoeffs = padWaveletXform2(NULL, NULL, sqrt(binCnt), NULL))) { + sprintf(errmsg, "can't determine number of wavelet coefficients " + "for modifier %s", mod + ); + error(INTERNAL, errmsg); + } + + if (numCoeffs * numCoeffs > PMAP_CONTRIB_MAXCOEFFS) { + sprintf(errmsg, "too many coefficients (max %d) for modifier %s; " + "reduce -bn", PMAP_CONTRIB_MAXCOEFFS, mod + ); error(USER, errmsg); - } + } + 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 (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, 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, 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; #if 0 worldfunc(RCCONTEXT, &srcRay); #endif set_eparams((char*)srcCont -> params); if ((bin = ray2bin(&srcRay, sqrt(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 ); } 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( pmapContribTab, 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 ); error(WARNING, errmsg); } } return pmap -> numContribSrc; } PhotonContribSourceIdx buildContribSources (PhotonMap *pmap, FILE **contribSrcHeap, char **contribSrcHeapFname, PhotonContribSourceIdx *contribSrcOfs, unsigned numHeaps ) /* 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; } /* ----------------- CONTRIB PMAP ENCODING STUFF -------------------- */ static int coeffCompare (const void *c1, const void *c2) /* Comparison function to REVERSE sort thresholded coefficients */ { const PreComputedContribCoeff *tcoeff1 = c1, *tcoeff2 = c2; /* Use dot product as magnitude to compare _absolute_ values */ const WAVELET_COEFF v1 = DOT(tcoeff1 -> coeff, tcoeff1 -> coeff), v2 = DOT(tcoeff2 -> coeff, tcoeff2 -> coeff); if (v1 < v2) return 1; else if (v1 > v2) return -1; else return 0; } static int thresholdContribs (PreComputedContrib *preCompContrib) /* Threshold wavelet detail coefficients in preCompContrib -> waveletMatrix [1..coeffDim-1] [1..coeffDim-1] by keeping the (preCompContrib -> nCompressedCoeffs) largest of these and returning them in preCompContrib -> compressedCoeffs along with their original order. NOTE: preCompContrib -> waveletMatrix [0][0] is the average wavelet coefficient and excluded from thresholding. Returns 0 on success. */ { unsigned i, j, coeffDim, coeffIdx; WaveletMatrix2 waveletMatrix; PreComputedContribCoeff *threshCoeffs; if (!preCompContrib || !(coeffDim = preCompContrib -> coeffDim) || !(threshCoeffs = preCompContrib -> compressedCoeffs) || !(waveletMatrix = preCompContrib -> waveletMatrix) ) /* Should be initialised by caller! */ return -1; /* Set up coefficient buffer for compression (thresholding), skipping * coefficient at waveletMatrix [0][0] since it's the average and * therefore not thresholded. The 2D coefficient indices (matrix * coords) are linearised to 1D. */ for (i = 0; i < coeffDim; i++) for (j = 0; j < coeffDim; j++) { coeffIdx = PMAP_CONTRIB_XY2LIN(coeffDim, i, j); if (coeffIdx) { /* Set up pointers to coeffs (sorted faster than 3 doubles) and remember original (linearised) index prior to sort */ threshCoeffs [coeffIdx - 1].idx = coeffIdx; threshCoeffs [coeffIdx - 1].coeff = (WAVELET_COEFF*)&( waveletMatrix [i] [j] ); } } /* REVERSE sort coeffs by magnitude; the non-thresholded coeffs will then start at threshCoeffs [0] (note nCoeffs == coeffDim^2) */ /* XXX: PARTITIONING WOULD ACTUALLY BE SUFFICIENT AND FASTER! */ qsort(threshCoeffs, preCompContrib -> nCoeffs - 1, sizeof(PreComputedContribCoeff), coeffCompare ); return 0; } static int encodeContribs (PreComputedContrib *preCompContrib, float compressRatio ) /* Apply wavelet transform to input matrix preCompContrib -> waveletMatrix and compress according to compressRatio, storing thresholded and mRGBE-encoded coefficients in preCompContrib -> mrgbeCoeffs. Note that the average coefficient is not encoded, and returned in preCompContrib -> waveletMatrix [0][0]. Returns 0 on success. */ { unsigned i, j, k, scDim; WaveletMatrix2 waveletMatrix, tWaveletMatrix; PreComputedContribCoeff *threshCoeffs; mRGBERange *mrgbeRange; mRGBE *mrgbeCoeffs; WAVELET_COEFF absCoeff; #ifdef PMAP_CONTRIB_DBG WaveletCoeff3 decCoeff; unsigned decIdx; #endif if (!preCompContrib || !preCompContrib -> mrgbeCoeffs || !preCompContrib -> compressedCoeffs || !(waveletMatrix = preCompContrib -> waveletMatrix) || !(tWaveletMatrix = preCompContrib -> tWaveletMatrix) || !(scDim = preCompContrib -> scDim) ) /* Should be initialised by the caller! */ return -1; #ifdef PMAP_CONTRIB_DBG for (i = 0; i < scDim; i++) { for (j = 0; j < scDim; j++) { for (k = 0; k < 3; k++) { #if 0 /* Set contributions to bins for debugging */ waveletMatrix [i] [j] [k] = PMAP_CONTRIB_XY2LIN(scDim, i, j); #else /* Replace contribs with "bump" function */ waveletMatrix [i] [j] [k] = (1. - fabs(i - scDim/2. + 0.5) / (scDim/2. - 0.5)) * (1. - fabs(j - scDim/2. + 0.5) / (scDim/2. - 0.5)); #endif } #if 0 /* Dump contribs prior to encoding for debugging */ printf("% 7.3g\t", colorAvg(waveletMatrix [i] [j])); } putchar('\n'); } putchar('\n'); #else } } #endif #endif /* Do 2D wavelet transform on preCompContrib -> waveletMatrix */ if (padWaveletXform2(waveletMatrix, tWaveletMatrix, scDim, NULL) != preCompContrib -> coeffDim ) error(INTERNAL, "failed wavelet transform in encodeContribs()"); /* Compress wavelet detail coeffs by thresholding */ if (thresholdContribs(preCompContrib) < 0) error(INTERNAL, "failed wavelet compression in encodeContribs()"); threshCoeffs = preCompContrib -> compressedCoeffs; /* Init per-channel coefficient range for mRGBE encoding */ mrgbeRange = &preCompContrib -> mrgbeRange; setcolor(mrgbeRange -> min, FHUGE, FHUGE, FHUGE); setcolor(mrgbeRange -> max, 0, 0, 0); /* Update per-channel coefficient range */ for (i = 0; i < preCompContrib -> nCompressedCoeffs; i++) { for (k = 0; k < 3; k++) { #ifdef PMAP_CONTRIB_DBG #if 0 /* Replace wavelet coefficient with bin for debugging, ordering bins linearly */ threshCoeffs [i].coeff [k] = threshCoeffs [i].idx = i + 1; #endif #endif absCoeff = fabs(threshCoeffs [i].coeff [k]); if (absCoeff < mrgbeRange -> min [k]) mrgbeRange -> min [k] = absCoeff; if (absCoeff > mrgbeRange -> max [k]) mrgbeRange -> max [k] = absCoeff; } } if (preCompContrib -> nCompressedCoeffs == 1) /* Maximum compression with just 1 (!) compressed detail coeff (is * this even useful?), so mRGBE range is undefined since min & max * coincide; set minimum to 0, maximum to the single remaining * coefficient */ setcolor(mrgbeRange -> min, 0, 0, 0); /* Init mRGBE coefficient normalisation from range */ if (!mRGBEinit(mrgbeRange, mrgbeRange -> min, mrgbeRange -> max)) return -1; mrgbeCoeffs = preCompContrib -> mrgbeCoeffs; /* Encode wavelet detail coefficients to mRGBE */ for (i = 0; i < preCompContrib -> nCompressedCoeffs; i++) { /* XXX: threshCoeffs [i].idx could overflow the mRGBE data field, * since there are more padded coefficients than bins!!! * SORT BY IDX AND USE AN INCREMENTAL INDEXING HERE INSTEAD??? */ mrgbeCoeffs [i] = mRGBEencode(threshCoeffs [i].coeff, mrgbeRange, threshCoeffs [i].idx ); + if (!mrgbeCoeffs [i].all) + error(INTERNAL, "failed mRGBE encoding in encodeContribs()"); #ifdef PMAP_CONTRIB_DBG /* Encoding sanity check */ decIdx = mRGBEdecode(mrgbeCoeffs [i], mrgbeRange, decCoeff); #if 0 if (decCoeff != threshCoeffs [i].idx || sqrt(dist2(decCoeff, threshCoeffs [i].coeff)) > 1.7 ) error(CONSISTENCY, "failed sanity check in encodeContribs()"); #endif for (k = 0; k < 3; k++) if (decCoeff [k] / threshCoeffs [i].coeff [k] < 0) error(CONSISTENCY, "coefficient sign reversal in encodeContribs()" ); #endif } return 0; } 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 initContribHeap()" ); #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( pmapContribTab, 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 distance */ scalecolor(flux, 1 - sqn -> dist2 / r2); #endif addcolor(irrad, flux); addcolor(contrib -> cbin [photonSrcBin(pmap, photon)], flux); } return contrib; } void freePreCompContribNode (void *p) /* Clean up precomputed contribution LUT entry */ { PhotonMap *preCompContribPmap = (PhotonMap*)p; PreComputedContrib *preCompContrib; if (preCompContribPmap) { preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib ); if (preCompContrib) { /* Free primary and transposed wavelet matrices */ freeWaveletMatrix2(preCompContrib -> waveletMatrix, preCompContrib -> coeffDim ); freeWaveletMatrix2(preCompContrib -> tWaveletMatrix, preCompContrib -> coeffDim ); /* Free thresholded coefficients */ free(preCompContrib -> compressedCoeffs); /* Free mRGBE encoded coefficients */ free(preCompContrib -> mrgbeCoeffs); free(preCompContrib -> waveletFname); if (preCompContrib -> cache) { /* Free contribution cache */ OOC_DeleteCache(preCompContrib -> cache); free(preCompContrib -> cache); } } /* Clean up precomputed contrib photon map */ deletePhotons(preCompContribPmap); free(preCompContribPmap); } } void initPreComputedContrib (PreComputedContrib *preCompContrib) /* Initialise precomputed contribution container in photon map */ { preCompContrib -> waveletFname = NULL; preCompContrib -> waveletFile = NULL; preCompContrib -> waveletMatrix = preCompContrib -> tWaveletMatrix = NULL; preCompContrib -> compressedCoeffs = NULL; setcolor(preCompContrib -> mrgbeRange.min, 0, 0, 0); setcolor(preCompContrib -> mrgbeRange.max, 0, 0, 0); setcolor(preCompContrib -> mrgbeRange.norm, 0, 0, 0); preCompContrib -> mrgbeCoeffs = NULL; preCompContrib -> scDim = preCompContrib -> nBins = preCompContrib -> coeffDim = preCompContrib -> nCoeffs = preCompContrib -> nNonZeroCoeffs = preCompContrib -> nCompressedCoeffs = preCompContrib -> contribSize = 0; preCompContrib -> cache = NULL; } void preComputeContrib (PhotonMap *pmap) /* Precompute contributions for a random subset of (finalGather * pmap -> numPhotons) photons, and return the per-modifier precomputed contribution photon maps in LUT preCompContribTab, discarding the original photons. */ { unsigned long p, numPreComp; unsigned i, scDim, nBins, coeffDim, nCoeffs, nCompressedCoeffs; PhotonIdx pIdx; Photon photon; RAY ray; MODCONT *binnedContribs; COLR mrgbeRange32 [2]; LUENT *preCompContribNode; PhotonMap *preCompContribPmap, nuPmap; PreComputedContrib *preCompContrib; LUTAB lutInit = LU_SINIT(NULL, freePreCompContribNode); WaveletMatrix2 waveletMatrix; #ifdef PMAP_CONTRIB_DBG RAY dbgRay; #endif if (verbose) { sprintf(errmsg, "\nPrecomputing contributions for %ld photons\n", numPreComp ); eputs(errmsg); #if NIX fflush(stderr); #endif } /* Init new parent photon map and set output filename */ initPhotonMap(&nuPmap, pmap -> type); nuPmap.fileName = pmap -> fileName; /* Set contrib/coeff mode in new parent photon map */ nuPmap.contribMode = pmap -> contribMode; /* Allocate and init LUT containing per-modifier child photon maps */ if (!(nuPmap.preCompContribTab = malloc(sizeof(LUTAB)))) error(SYSTEM, "out of memory allocating LUT in preComputeContrib()" ); memcpy(nuPmap.preCompContribTab, &lutInit, sizeof(LUTAB)); /* 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 (p = 0; p < numPreComp; p++) { /* Get random photon from stratified distribution in source heap * to avoid duplicates and clustering */ pIdx = firstPhoton(pmap) + (unsigned long)( (p + 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 (i = 0; i < 3; i++) ray.ron [i] = photon.norm [i] / 127.0; #ifndef PMAP_CONTRIB_DBG /* Get contributions at photon position; retry with another photon if no contribs found */ binnedContribs = getPhotonContrib(pmap, &ray, ray.rcol); #else /* Get all contribs from same photon for debugging */ /* Position must differ otherwise too many identical photon keys * will cause ooC octree to overflow! */ VCOPY(dbgRay.rop, photon.pos); getPhoton(pmap, 0, &photon); memcpy(&dbgRay, &ray, sizeof(RAY)); dbgRay.rsrc = photonSrcIdx(pmap, &photon); for (i = 0; i < 3; i++) dbgRay.ron [i] = photon.norm [i] / 127.0; binnedContribs = getPhotonContrib(pmap, &dbgRay, dbgRay.rcol); #endif if (binnedContribs) { #if 0 if (!(p & PMAP_CNTUPDATE)) printf("Precomputing contribution photon %lu / %lu\n", p, numPreComp ); #endif /* Found contributions at photon position, so generate * precomputed photon */ if (!(preCompContribNode = lu_find(nuPmap.preCompContribTab, binnedContribs -> modname) ) ) error(SYSTEM, "out of memory allocating LUT entry in " "preComputeContrib()" ); if (!preCompContribNode -> key) { /* New LUT node for precomputed contribs for this modifier */ preCompContribNode -> key = (char*)binnedContribs -> modname; preCompContribNode -> data = (char*)( preCompContribPmap = malloc(sizeof(PhotonMap)) ); if (preCompContribPmap) { /* Init new child photon map and its contributions */ initPhotonMap(preCompContribPmap, PMAP_TYPE_CONTRIB_CHILD ); initPhotonHeap(preCompContribPmap); initContribHeap(preCompContribPmap); preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib = malloc(sizeof(PreComputedContrib)) ); initPreComputedContrib(preCompContrib); } if (!preCompContribPmap || !preCompContrib) error(SYSTEM, "out of memory allocating new photon map " "in preComputeContrib()" ); /* Set output filename from parent photon map */ preCompContribPmap -> fileName = nuPmap.fileName; /* Set contrib/coeff mode from parent photon map */ preCompContribPmap -> contribMode = nuPmap.contribMode; /* Set Shirley-Chiu square & wavelet matrix dimensions (number of bins resp. coefficients). */ preCompContrib -> scDim = scDim = sqrt( preCompContrib -> nBins = nBins = binnedContribs -> nbins ); if (scDim * scDim != nBins) /* Mkpmap shoulda took care of dis */ error(INTERNAL, "contribution bin count not " "integer square in preComputeContrib()" ); if (nBins > 1) { /* BINNING ENABLED; set up wavelet xform & compression. Get dimensions of wavelet coefficient matrix after boundary padding, and number of nonzero coefficients. The number of compressed coefficients is fixed and based on the number of nonzero coefficients (minus one as the average coefficient is not thresholded), since zero coeffs are already thresholded by default. */ preCompContrib -> coeffDim = coeffDim = padWaveletXform2( NULL, NULL, scDim, &preCompContrib -> nNonZeroCoeffs ); preCompContrib -> nCoeffs = nCoeffs = coeffDim * coeffDim; nCompressedCoeffs = (1 - compressRatio) * (preCompContrib -> nNonZeroCoeffs - 1); if (nCompressedCoeffs < 1) { error(WARNING, "maximum contribution compression, clamping number " "of wavelet coefficients in preComputeContrib()" ); nCompressedCoeffs = 1; } if (nCompressedCoeffs >= preCompContrib -> nCoeffs) { error(WARNING, "minimum contribution compression, clamping number " "of wavelet coefficients in preComputeContrib()" ); nCompressedCoeffs = nCoeffs - 1; } preCompContrib -> nCompressedCoeffs = nCompressedCoeffs; /* Lazily allocate primary and transposed wavelet * coefficient matrix */ preCompContrib -> waveletMatrix = waveletMatrix = allocWaveletMatrix2(coeffDim); preCompContrib -> tWaveletMatrix = allocWaveletMatrix2( coeffDim ); if (!waveletMatrix || !preCompContrib -> tWaveletMatrix) error(SYSTEM, "out of memory allocating wavelet " "coefficient matrix in preComputeContrib()" ); /* Lazily allocate thresholded detail coefficients (minus the average coeff) */ preCompContrib -> compressedCoeffs = calloc( nCoeffs - 1, sizeof(PreComputedContribCoeff) ); /* Lazily alloc mRGBE-encoded compressed wavelet coeffs */ preCompContrib -> mrgbeCoeffs = calloc( nCompressedCoeffs, sizeof(mRGBE) ); if (!preCompContrib -> compressedCoeffs || !preCompContrib -> mrgbeCoeffs ) error(SYSTEM, "out of memory allocating compressed " "contribution buffer in preComputeContrib()" ); /* Set size of compressed contributions, in bytes */ preCompContrib -> contribSize = PMAP_CONTRIB_ENCSIZE( nCompressedCoeffs ); } } else { /* LUT node already exists for this modifier; use its data */ preCompContribPmap = (PhotonMap*)preCompContribNode -> data; preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib ); scDim = preCompContrib -> scDim; nBins = preCompContrib -> nBins; nCoeffs = preCompContrib -> nCoeffs; nCompressedCoeffs = preCompContrib -> nCompressedCoeffs; waveletMatrix = preCompContrib -> waveletMatrix; } if (nBins > 1) { /* Binning enabled; copy binned contribs to wavelet xform * input matrix, then compress & encode */ for (i = 0; i < scDim; i++) memcpy(waveletMatrix [i], binnedContribs->cbin + PMAP_CONTRIB_XY2LIN(scDim, i, 0), scDim * WAVELET_COEFFSIZE ); if (encodeContribs(preCompContrib, compressRatio)) error(INTERNAL, "failed contribution " "compression/encoding in preComputeContrib()" ); /* Dump encoded wavelet coeffs to contrib heap file, prepended by mRGBE range in 32-bit RGBE. NOTE: minimum and maximum are reversed in the file to facilitate handling the special (but pointless) case of a single compressed coeff (see below) */ setcolr(mrgbeRange32 [0], preCompContrib -> mrgbeRange.max [0], preCompContrib -> mrgbeRange.max [1], preCompContrib -> mrgbeRange.max [2] ); setcolr(mrgbeRange32 [1], preCompContrib -> mrgbeRange.min [0], preCompContrib -> mrgbeRange.min [1], preCompContrib -> mrgbeRange.min [2] ); /* Dump only mRGBE maximum if 1 compressed coeff (maximum * compression, pointless as it may be!), since minimum is * implicitly zero and can be omitted to save space */ putbinary(mrgbeRange32, sizeof(COLR), 1 + (nCompressedCoeffs > 1), preCompContribPmap -> contribHeap ); if (putbinary(preCompContrib -> mrgbeCoeffs, sizeof(mRGBE), nCompressedCoeffs, preCompContribPmap -> contribHeap ) != nCompressedCoeffs ) error(SYSTEM, "failed writing to contribution heap " "in preComputeContrib()" ); } - /* Set photon flux to coarsest average wavelet coefficient - waveletMatrix [0] [0]. IF BINNING IS DISABLED, this is the - total contribution from all directions */ - copycolor(ray.rcol, waveletMatrix [0] [0]); - /* HACK: signal newPhoton() to set precomputed photon's contribution source from ray -> rsrc */ /* XXX: DO WE STILL NEED THIS CRAP AFTER PRECOMP, SINCE CHILD PMAPS ARE PER-MODIFIER ANYWAY? */ preCompContribPmap -> lastContribSrc.srcIdx = -2; /* Append photon to new heap from ray and increment total count for all modifiers in parent photon map */ newPhoton(preCompContribPmap, &ray); nuPmap.numPhotons++; } /* Update progress */ repProgress++; if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) pmapPreCompReport(); #ifdef SIGCONT else signal(SIGCONT, pmapPreCompReport); #endif } /* Trash original pmap and its leaf file, and replace with new parent, which now acts as container for the per-modifier child pmaps */ unlink(pmap -> store.leafFname); deletePhotons(pmap); memcpy(pmap, &nuPmap, sizeof(PhotonMap)); #ifdef SIGCONT signal(SIGCONT, SIG_DFL); #endif /* Bail out if no photons could be precomputed */ if (!pmap -> numPhotons) error(USER, "no contribution photons precomputed; try increasing -am" ); /* Build per-modifier precomputed photon maps from their contribution heaps */ lu_doall(pmap -> preCompContribTab, buildPreCompContribPmap, NULL); } #else /* ------------------------------------------------------------------- U N I T T E S T S ------------------------------------------------------------------- */ #include #include #include "func.h" int main (int argc, char *argv []) { unsigned i, scdim, 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 Shirley-Chiu dimension scDim>1, " "number of samples nsamp\n", stderr ); return -1; } if (!(scdim = atoi(argv [1]))) { fputs("Invalid Shirley-Chiu dimension\n", stderr); return -1; } 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, scdim); printf("%.3f\t%.3f\t%.3f\t-->\t%d\n", ray.rdir [0], ray.rdir [1], ray.rdir [2], bin ); } return 0; } #endif #endif /* PMAP_CONTRIB */ diff --git a/pmapcontrib.h b/pmapcontrib.h index 046145c..fa7aa36 100644 --- a/pmapcontrib.h +++ b/pmapcontrib.h @@ -1,224 +1,224 @@ /* 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 #if defined(PMAP_OOC) && !defined(PMAP_CONTRIB) #define PMAP_CONTRIB #endif #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 per-modifier contrib photon maps and wavelet coeffs; these are held in a separate subdirectory PMAP_CONTRIB_DIR */ #define PMAP_CONTRIB_DIR "%s.rc" #define PMAP_CONTRIB_FILE "%s/%s.pm" #define PMAP_CONTRIB_WAVELETFILE "%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 bin dimension */ - #define PMAP_CONTRIB_MAXBINS (mRGBE_DATAMAX) + /* Maximum number of mRGBE encodable wavelet coefficients */ + #define PMAP_CONTRIB_MAXCOEFFS (mRGBE_DATAMAX + 1) /* Size of encoded coefficients in wavelet file (in bytes) as a function * of the number of compressed coefficients nComp. * NOTE: the coefficient range minimum is omitted from the mRGBE range if * nComp == 1, since it is implicitly zero! */ #define PMAP_CONTRIB_ENCSIZE(nComp) ( \ sizeof(COLR) * (1 + ((nComp) > 1)) + sizeof(mRGBE) * (nComp) \ ) /* Serialise 2D coordinates (x,y) in (dim x dim) Shirley-Chiu square to 1D index. Returns -1 if coordinates invalid */ #define PMAP_CONTRIB_XY2LIN(dim, x, y) ( \ (x) < 0 || (y) < 0 ? -1 : (x) * (dim) + (y) \ ) /* Deserialise 1D index idx to 2D coordinates (x,y) in (dim x dim) Shirley-Chiu square. Returns -1 if bin invalid */ #define PMAP_CONTRIB_LIN2X(dim, idx) ((idx) < 0 ? -1 : (idx) / (dim)) #define PMAP_CONTRIB_LIN2Y(dim, idx) ((idx) < 0 ? -1 : (idx) % (dim)) /* Struct for wavelet coeff thresholding; saves original coeff index */ typedef struct { WAVELET_COEFF *coeff; unsigned idx; } PreComputedContribCoeff; typedef struct { char *waveletFname; FILE *waveletFile; WaveletMatrix2 waveletMatrix, tWaveletMatrix; /* Wavelet coeff compression/encoding buffer */ PreComputedContribCoeff *compressedCoeffs; mRGBERange mrgbeRange; mRGBE *mrgbeCoeffs; unsigned scDim, nBins, coeffDim, nCoeffs, nCompressedCoeffs, nNonZeroCoeffs; unsigned long contribSize; OOC_Cache *cache; } PreComputedContrib; /* The following are convenient placeholders interfacing to mkpmap and rcontrib. They can be externally set via initPmapContribTab() and referenced within the contrib pmap modules. These variables can then be safely ignored by rtrace/rpict/rvu, without annoying linking errors. */ /* Global pointer to rcontrib's contribution binning LUT */ extern LUTAB *pmapContribTab; /* Contribution/coefficient mode flag */ extern int *pmapContribMode; #if 0 int xy2lin (unsigned scDim, int x, int y); /* Serialise 2D contribution coords in (scDim x scDim) Shirley-Chiu square to 1D index. Returns -1 if coordinates invalid */ void lin2xy (unsigned scDim, int bin, int *x, int *y); /* Deserialise 1D contribution index to 2D coordinates in (scDim x scDim) Shirley-Chiu square. Returns -1 if bin invalid */ #endif MODCONT *addContribModifier(LUTAB *contribTab, unsigned *numContribs, 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, int binCnt ); /* Add light source modifiers from file modFile to contribution lookup * table contribTab, and update numContribs. */ void freePreCompContribNode (void *p); /* Clean up precomputed contribution LUT entry */ void initPmapContribTab (LUTAB *contribTab, int *contribMode); /* Set contribution table and contrib/coeff mode flag (interface to rcmain.c, see also pmapContribTab above) */ void initPmapContrib (PhotonMap *pmap); /* Initialise contribution photon map and check modifiers. NOTE: pmapContribTab must be set before via initPmapContribTab() */ 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 initPreComputedContrib (PreComputedContrib *preCompContrib); /* Initialise precomputed contribution container in photon map */ void preComputeContrib (PhotonMap *pmap); /* Precompute contributions for a random subset of (finalGather * pmap -> numPhotons) photons, and init the per-modifier precomputed contribution photon maps in LUT pmap -> preCompContribTab, discarding the original photons. */ void distribPhotonContrib (PhotonMap *pmap, LUTAB *contribTab, unsigned numContribs, int *contribMode, unsigned numProc ); /* Emit photons with binned light source contributions, precompute their contributions and build photon map */ int buildPreCompContribPmap (const LUENT *preCompContribNode, void *p); /* Build per-modifier precomputed photon map. Returns 0 on success. */ void saveContribPhotonMap (const PhotonMap *pmap, const char *fname, int argc, char **argv ); /* Save contribution pmap and its per-modifier precomputed children */ void loadContribPhotonMap (PhotonMap *pmap, const char *fname); /* Load contribution pmap and its per-modifier precomputed children */ void getPreCompPhotonContrib (PhotonMap *pmap, RAY *ray, COLOR totalContrib ); /* Get precomputed light source contributions at ray -> rop for all specified modifiers and accumulate them in pmapContribTab (which maps to rcontrib's binning LUT). Also return total precomputed contribs. */ #endif diff --git a/pmapdata.h b/pmapdata.h index 9efb605..2d7e32d 100644 --- a/pmapdata.h +++ b/pmapdata.h @@ -1,433 +1,431 @@ /* 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 #else #define PMAP_PROCBITS 5 #endif #define PMAP_MAXPROC (1 << PMAP_PROCBITS) #define PMAP_GETRAYPROC(r) ((r) -> crtype >> 8) #define PMAP_SETRAYPROC(r,p) ((r) -> crtype |= p << 8) typedef struct { float pos [3]; /* Photon position */ signed char norm [3]; /* Encoded normal / incident direction [volume photons] */ union { struct { #ifndef PMAP_OOC unsigned char discr : 2; /* kd-tree discriminator axis */ #endif unsigned char caustic : 1; /* Specularly scattered (=caustic) */ /* Photon's generating subprocess index, used for primary ray * index linearisation when building contrib pmaps; note this is * reduced for kd-tree to accommodate the discriminator field */ unsigned char proc : PMAP_PROCBITS; }; unsigned char flags; }; - /* Photon flux in watts or lumen / - photon contribution [contrib photons] / - average wavelet coefficient [precomputed contrib photons] */ + /* Photon flux in watts or lumen */ #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) /* Forward declarations */ struct PhotonMap; struct PreComputedContrib; struct PhotonBiasCompNode; 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 *preCompContribTab; /* LUT for per-modifier precomp contrib child photon maps (in parent) or NULL (in child) */ struct PreComputedContrib *preCompContrib; /* Precomputed contribs (in child) or NULL (in parent) */ FILE *contribHeap; /* Out-of-core heap containing unsorted precomputed contrib photon bins prior to construction of store */ char contribHeapFname [sizeof(PMAP_TMPFNAME)]; int contribMode; /* If 0, photon map contains precomputed coefficients, else contributions */ /* ================================================================ * BIAS COMPENSATION STUFF * ================================================================ */ struct 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]) #else #define lightflowPmap NULL #define transLightFlowPmap NULL #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 isContribChildPmap(p) ((p) -> type == PMAP_TYPE_CONTRIB_CHILD) #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, void *photonIdx ); /* Find single closest photon to ray -> rop with similar normal. Return NULL if none found, else the supplied Photon* buffer and the photon's index/ID (if photonIdx != NULL), 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/wavelet2.h b/wavelet2.h index d10deb7..7cb6eb9 100644 --- a/wavelet2.h +++ b/wavelet2.h @@ -1,245 +1,245 @@ /* ========================================================================= Definitions for 2D wavelet transforms on arrays of 3-tuples. Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") ========================================================================= $Id$ */ #ifndef _WAVELET2_H #define _WAVELET2_H /* Wavelet coefficient type; defaults to double if not already defined. (may be overridden in compiler command line, for example). NOTE: single precision float may improve cache performance during 2D transform with large arrays. */ #ifndef WAVELET_COEFF #define WAVELET_COEFF double #endif /* Perform final Haar wavelet transform on coarsest two detail levels. * This has the advantage of resulting in a single coarsest level * approximation coefficient, possibly at the expense of poorer * compression */ /* #define WAVELET_FINAL_HAAR */ /* Boundary extension modes for padded Daubechies D4 wavelet transform */ #define WAVELET_EXTEND_CONST 0 /* Constant (default if undef) */ #define WAVELET_EXTEND_GRAD1 1 /* 1st order gradient */ #define WAVELET_EXTEND_GRAD2 2 /* 2nd order gradient */ #define WAVELET_EXTEND_CYCL 3 /* Periodic (wraparound) */ #define WAVELET_EXTEND_SYMM 4 /* Symmetric (reflection) */ #ifndef WAVELET_EXTEND - #define WAVELET_EXTEND_MODE WAVELET_EXTEND_SYMM + #define WAVELET_EXTEND_MODE WAVELET_EXTEND_GRAD1 #endif /* Convenience macros for handling coefficients */ #define min(a, b) ((a) < (b) ? (a) : (b)) #define max(a, b) ((a) > (b) ? (a) : (b)) #define coeffAvg(a) (((a) [0] + (a) [1] + (a) [2]) / 3) #define coeffIsEmpty(c) ( \ isnan((c) [0]) || isnan((c) [1]) || isnan((c) [2]) \ ) #define coeffIsZero(c) ((c) [0] == 0 && (c) [1] == 0 && (c) [2] == 0) /* Determine if coeff tuple c is thresholded, also including empty coeffs * in debugging mode, since these will be unconditionally dropped. Note * that the 3x3 coarsest approximation coefficients in the upper left of * matrix c are not thresholded, since they are essential for a * reasonably accurate reconstruction! */ /* #ifdef WAVELET_FINAL_HAAR */ #if 0 #define coeffThresh(c, i, j, t) ( \ fabs(t) > 0 && ((i) || (j)) && ( \ coeffIsEmpty((c) [i][j]) || fabs((c) [i][j][0]) < (t) || \ fabs((c) [i][j][1]) < (t) || fabs((c) [i][j][2]) < (t) \ ) \ ) #else #define coeffThresh(c, i, j, t) ( \ fabs(t) > 0 && ((i) > 2 || (j) > 2) && ( \ coeffIsEmpty((c) [i][j]) || fabs((c) [i][j][0]) < (t) || \ fabs((c) [i][j][1]) < (t) || fabs((c) [i][j][2]) < (t) \ ) \ ) #endif /* Haar / Daubechies D4 coefficients */ #define SQRT2 1.41421356 #define SQRT3 1.73205081 #define H4NORM (0.25 / SQRT2) /* Haar wavelet coeffs */ extern const WAVELET_COEFF h2; /* Daubechies D4 wavelet coeffs */ extern const WAVELET_COEFF h4 [4]; extern const WAVELET_COEFF g4 [4]; /* Wavelet matrix defs; these are stored as arrays of pointers (a.k.a. * Iliffe vectors), as this proved to be more performant than a flattened * representation. Encapsulating the coefficient's triplets in a struct * prevents the compiler from introducing another layer of indirection. * */ typedef WAVELET_COEFF WaveletCoeff3 [3]; #define WAVELET_COEFFSIZE (sizeof(WaveletCoeff3)) typedef WaveletCoeff3 **WaveletMatrix2; /* ----------- 2D COEFFICIENT MATRIX ALLOC/INIT ROUTINES ----------- */ WaveletMatrix2 allocWaveletMatrix2 (unsigned len); /* Allocate and return a 2D coefficient array of size (len x len). Returns NULL if allocation failed. */ void freeWaveletMatrix2 (WaveletMatrix2 y, unsigned len); /* Free previously allocated 2D coeff array y of size (len x len) */ void zeroCoeffs2 (WaveletMatrix2 y, unsigned len); /* Set 2D array coefficients to zero */ /* ---------- 2D WAVELET TRANSFORM OPTIMISED FOR SIZES 2^l ---------- ------ THESE USE PERIODIC (WRAP-AROUND) BOUNDARY EXTENSION ------- */ int waveletXform2 (WaveletMatrix2 y, WaveletMatrix2 yt, unsigned l); /* Perform full 2D multiresolution forward wavelet transform on array y of size (2^l) x (2^l) containing original signal as 3-tuples, where l >= 1. Note no intra-tuple transform occurs. The wavelet coefficients are returned in array y, containing the coarsest approximation in y [0][0] followed by horizontal/vertical details in order of increasing resolution/frequency. A preallocated array yt of identical dimensions to y can be supplied as buffer for intermediate results. If yt == NULL, a buffer is automatically allocated and freed on demand, but this is inefficient for frequent calls. It is recommended to preallocate yt to the maximum expected size. The dimensions of yt are not checked; this is the caller's responsibility. Returns 0 on success, else -1. */ int waveletInvXform2 (WaveletMatrix2 y, WaveletMatrix2 yt, unsigned l); /* Perform full 2D multiresolution inverse wavelet transform on array y of size (2^l) x (2^l) containing wavelet coefficients as 3-tuples, where l >= 1. Note no intra-tuple transform occurs. A preallocated array yt of identical dimensions to y can be supplied as buffer for intermediate results. If yt == NULL, a buffer is automatically allocated and freed on demand, but this is inefficient for frequent calls. It is recommended to preallocate yt to the maximum expected size. The dimensions of yt are not checked; this is the caller's responsibility. The reconstructed signal is returned in array y. Returns 0 on success, else -1. */ /* ------------ 2D WAVELET TRANSFORM FOR ARBITRARY SIZES ------------- ---- THESE USE CONSTANT BOUNDARY EXTENSION WITH PADDING COEFFS----- */ unsigned padWaveletXform2 (WaveletMatrix2 y, WaveletMatrix2 yt, unsigned yLen0, unsigned *numNonZero ); /* Perform full 2D multiresolution forward wavelet transform with padded boundary coefficients on array y containing (yLen0 x yLen0) samples of 3-tuples. Note no intra-tuple transform occurs. An additional array yt of identical dimension to y is required as buffer for intermediate results. The wavelet coefficients are returned in array y, containing the coarsest approximation coefficients in y [0][0] followed by horizontal/vertical detail coefficients in order of increasing resolution/frequency. NOTE: Both y and yt must be large enough to accommodate the additional padding coefficients introduced at the array boundaries during the transform. It is the caller's responsibility to preallocate and dimension these arrays accordingly with allocWaveletMatrix2(). The dimension of the output array (which generally exceeds the input dimension yLen0) must be interrogated beforehand by calling the function with y or yt set to NULL. The returned value is the output dimension, or 0 if the transform failed. In addition, the number of nonzero coefficients is returned in numNonZero (if not NULL), which is a subset of the number of output coefficients, i.e. (output dimension)^2. */ unsigned padWaveletInvXform2 (WaveletMatrix2 y, WaveletMatrix2 yt, unsigned yLen0, unsigned yLenOrg ); /* Perform full 2D multiresolution inverse wavelet transform on array y containing (yLen0 x yLen0) boundary padded wavelet coefficients as 3-tuples. The original number of samples (yLenOrg x yLenOrg) is required for the reconstruction. An additional array yt of identical dimension to y is required as buffer for intermediate results. The reconstructed samples are returned in y. This is smaller than (yLen0 x yLen0) since the boundary coefficients from the forward transform are collapsed. The returned value is the reconstructed number of samples per dimension (yLenOrg), or 0 if the inverse transform failed. */ #ifdef WAVELET_DBG /* ----------------- SUPPORT ROUTINES FOR DEBUGGING ----------------- */ void clearCoeffs2 (WaveletMatrix2 y, unsigned len); /* Clear 2D array to facilitate debugging */ void dumpCoeffs2 (const WaveletMatrix2 y1, const WaveletMatrix2 y2, unsigned y1Len, unsigned y2Len, float thresh ); /* Dump 2D arrays y1 and y2 of dims y1Len resp. y2Len side-by-side to stdout (skip y2 if NULL). Coefficients with absolute magnitude less than thresh are marked with brackets ('[]') as thresholded. */ float rmseCoeffs2 (const WaveletMatrix2 y1, const WaveletMatrix2 y2, unsigned len ); /* Calculate RMSE between 2D matrices y1 and y2 */ #endif #endif