diff --git a/ooccache.h b/ooccache.h index 1ae4a19..4be39e1 100644 --- a/ooccache.h +++ b/ooccache.h @@ -1,97 +1,97 @@ /* ======================================================================= Cache for out-of-core octree. Uses an LRU page replacement scheme. Page table nodes are stored in order of reference in a doubly linked list, with the most recently used (MRU) page at the head, and the least recently used (LRU) at the tail. Nodes are referenced by the page's base address via a hashtable for constant time access (in the absence of hash collisions). Hash collisions can be limited by adjusting the hashtable load factor OOC_CACHE_LOAD; this is the fraction of the number of pages that will actually be cached. 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) ======================================================================= $Id: ooccache.h,v 2.1 2016/05/17 17:39:47 rschregle Exp $ */ #ifndef OOC_CACHE_H #define OOC_CACHE_H #include #include /* Hashtable load factor to limit hash collisions */ #define OOC_CACHE_LOAD 0.75 /* Minimum number of cache pages */ #define OOC_CACHE_MINPAGES 8 /* Hashing function for key -> pagetable mapping (assumes n is prime) */ #define OOC_CACHE_PRIME 19603 #define OOC_CACHE_HASH(k,n) ((k) * OOC_CACHE_PRIME % (n)) /* Next index in pagetable probe sequence for collision resolution */ #define OOC_CACHE_NEXT(i,n) (((i) + 1) % (n)) /* Cache key and index types */ typedef uint32_t OOC_CacheKey; - typedef uint32_t OOC_CacheIdx; + typedef uint32_t OOC_CacheIdx; /* Undefined value for cache index */ #define OOC_CACHEIDX_NULL UINT32_MAX typedef struct { OOC_CacheKey key; /* Key to resolve hashing collisions */ OOC_CacheIdx prev, next; /* Previous/next page in history */ void *data; /* Start of page data */ } OOC_CacheNode; typedef struct { unsigned recSize, /* Data record size in bytes */ recPerPage; /* Page size in number of records */ unsigned long pageSize; /* Page size in bytes (precomp) */ OOC_CacheIdx pageCnt, /* Num pages used */ numPages; /* Pagetable size */ /* NOTE: The hashtable load factor limits the number of pages actually * used, such that pageCnt <= OOC_CACHE_LOAD * numPages */ OOC_CacheNode *pageTable; /* Pagetable with numPages nodes */ OOC_CacheIdx mru, lru; /* Most/least recently used page == head/tail of page history */ OOC_CacheKey lastKey; /* Previous key to detect repeat lookups */ OOC_CacheNode *lastPage; /* Previous page for repeat lookups */ unsigned long numReads, /* Statistics counters */ numHits, numColl, numRept; } OOC_Cache; /* Initialise cache containing up to numPages pages; this is rounded up * the next prime number to reduce hash clustering and collisions. The * page size recPerPage is specified in units of records of size recSize * bytes. Returns 0 on success, else -1. */ int OOC_CacheInit (OOC_Cache *cache, unsigned numPages, unsigned recPerPage, unsigned recSize); /* Return pointer to cached record with index recIdx, loading from file * and possibly evicting the LRU page */ void *OOC_CacheData (OOC_Cache *cache, FILE *file, unsigned long recIdx); /* Delete cache and free allocated pages */ void OOC_DeleteCache (OOC_Cache *cache); #endif diff --git a/oococt.c b/oococt.c index 718e8c4..2348458 100644 --- a/oococt.c +++ b/oococt.c @@ -1,482 +1,484 @@ #ifndef lint static const char RCSid[] = "$Id: oococt.c,v 2.4 2017/08/14 21:12:10 rschregle Exp $"; #endif /* ====================================================================== Out-of-core octree data structure 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) ====================================================================== $Id: oococt.c,v 2.4 2017/08/14 21:12:10 rschregle Exp $ */ #if !defined(_WIN32) && !defined(_WIN64) || defined(PMAP_OOC) /* No Windoze support for now */ #include "oococt.h" #include "rtio.h" #include #include void OOC_Null (OOC_Octree *oct) /* Init empty octree */ { oct -> maxNodes = oct -> numNodes = oct -> leafMax = oct -> maxDepth = oct -> numData = oct -> size = oct -> recSize = oct -> mortonScale = 0; oct -> org [0] = oct -> org [1] = oct -> org [2] = oct -> bound [0] = oct -> bound [1] = oct -> bound [2] = 0; oct -> key = NULL; oct -> nodes = NULL; oct -> leafFile = NULL; + oct -> leafFname = NULL; oct -> cache = NULL; } void OOC_Init (OOC_Octree *oct, unsigned recSize, FVECT org, RREAL size, - RREAL *(*key)(const void*), FILE *leafFile + RREAL *(*key)(const void*), FILE *leafFile, char *leafFname ) { oct -> maxNodes = oct -> numNodes = oct -> leafMax = oct -> maxDepth = oct -> numData = 0; VCOPY(oct -> org, org); oct -> size = oct -> bound[0] = oct -> bound[1] = oct -> bound[2] = size; VADD(oct -> bound, oct -> bound, org); oct -> mortonScale = size > 0 ? MORTON3D_MAX / size : 0; oct -> recSize = recSize; oct -> key = key; oct -> nodes = NULL; oct -> leafFile = leafFile; + oct -> leafFname = leafFname; oct -> cache = NULL; /* Cache currently initialised externally */ } OOC_Node *NewNode (OOC_Octree *oct) /* Allocate new octree node, enlarge array if necessary. Return pointer to new node or NULL if failed. */ { OOC_Node *n = NULL; if (oct -> numNodes >= OOC_NODEIDX_MAX) { /* Node index overflow */ fprintf(stderr, "OOC_NewNode: node index overflow (numNodes = %d)\n", oct -> numNodes ); return NULL; } if (++oct -> numNodes > oct -> maxNodes) { /* Need to allocate root or enlarge array */ oct -> maxNodes += OOC_BLKSIZE / sizeof(OOC_Node); oct -> nodes = realloc(oct -> nodes, oct -> maxNodes * sizeof(OOC_Node) ); if (!oct -> nodes) { perror("OOC_NewNode: couldn't realloc() nodes"); return NULL; } } n = oct -> nodes + oct -> numNodes - 1; n -> node.num = n -> node.kid = n -> node.branch = n -> node.type = 0; return n; } int OOC_GetData (OOC_Octree *oct, OOC_DataIdx idx, void *rec) /* Reads indexed data record from leaf file and copies it to rec, else * returns -1 on failure */ { if (!oct || !rec) return -1; if (idx >= oct -> numData) { fprintf(stderr, "OOC_GetData: invalid data record index %u\n", idx); return -1; } if (oct -> cache) { /* Retrieve record from leaf file via I/O cache */ void *cachedRec; if (!(cachedRec = OOC_CacheData(oct -> cache, oct -> leafFile, idx))) return -1; /* It's safer to copy the record to the caller's local buffer as a * returned pointer may be invalidated by subsequent calls if the * page it points to is swapped out */ memcpy(rec, cachedRec, oct -> recSize); } else { /* No I/O caching; do (SLOW!) per-record I/O from leaf file */ const unsigned long pos = (unsigned long)oct -> recSize * idx; if (pread(fileno(oct -> leafFile), rec, oct -> recSize, pos) != oct -> recSize) { perror("OOC_GetData: failed seek/read in leaf file"); return -1; } } return 0; } int OOC_InBBox (const OOC_Octree *oct, const FVECT org, RREAL size, const FVECT key ) /* Check whether key lies within bounding box defined by org and size. * Compares Morton codes rather than the coordinates directly to account * for dicretisation error introduced by the former. */ { const MortonIdx keyZ = OOC_KEY2MORTON(key, oct); FVECT bbox = {size, size, size}; VADD(bbox, org, bbox); return keyZ > OOC_KEY2MORTON(org, oct) && keyZ < OOC_KEY2MORTON(bbox, oct); } int OOC_Branch (const OOC_Octree *oct, const FVECT org, RREAL size, const FVECT key, FVECT nuOrg ) /* Return index of octant containing key and corresponding origin if nuOrg * != NULL, or -1 if key lies outside all octants. Size is already assumed * to be halved, i.e. corresponding to the length of an octant axis. * Compares Morton codes rather than the coordinates directly to account * for dicretisation error introduced by the former. */ { int o; FVECT octOrg; for (o = 0; o < 8; o++) { VCOPY(octOrg, org); OOC_OCTORIGIN(octOrg, o, size); if (OOC_InBBox(oct, octOrg, size, key)) { if (nuOrg) VCOPY(nuOrg, octOrg); return o; } } return -1; } int OOC_BranchZ (const OOC_Octree *oct, const FVECT org, RREAL size, MortonIdx keyZ, FVECT nuOrg ) /* Optimised version of OOC_Branch() with precomputed key Morton code */ { int o; const FVECT cent = {size, size, size}; FVECT octOrg, bbox; for (o = 0; o < 8; o++) { VCOPY(octOrg, org); OOC_OCTORIGIN(octOrg, o, size); VADD(bbox, octOrg, cent); if (keyZ > OOC_KEY2MORTON(octOrg, oct) && keyZ < OOC_KEY2MORTON(bbox, oct) ) { if (nuOrg) VCOPY(nuOrg, octOrg); return o; } } return -1; } OOC_DataIdx OOC_GetKid (const OOC_Octree *oct, OOC_Node **node, unsigned kid) /* For a given octree node, return the sum of the data counters for kids * [0..k-1]. On return, the node either points to the k'th kid, or * NULL if the kid is nonexistent or the node is a leaf */ { unsigned k; OOC_Node *kidNode = OOC_KID1(oct, *node); /* NULL if leaf */ OOC_DataIdx dataIdx = 0; for (k = 0; k < kid; k++) { if (OOC_ISLEAF(*node)) /* Sum data counters of leaf octants */ dataIdx += (*node) -> leaf.num [k]; else /* Sum data counter of inner node's nonempty kids and advance * pointer to sibling */ if (OOC_ISBRANCH(*node, k)) { dataIdx += OOC_DATACNT(kidNode); kidNode++; } } /* Update node pointer to selected kid (or NULL if nonexistent or node is * a leaf) */ *node = OOC_ISBRANCH(*node, kid) ? kidNode : NULL; return dataIdx; } #ifdef DEBUG_OOC int OOC_Check (OOC_Octree *oct, const OOC_Node *node, FVECT org, RREAL size, OOC_DataIdx dataIdx ) /* Traverse tree & check for valid nodes; oct -> leafFile must be open; * return 0 if ok, otherwise -1 */ { unsigned k; FVECT kidOrg; const RREAL kidSize = size * 0.5; if (!oct || !node) return -1; if (OOC_ISLEAF(node)) { /* Node is a leaf; check records in each octant */ char rec [oct -> recSize]; /* Violates C89! */ OOC_OctCnt d; RREAL *key; for (k = 0; k < 8; k++) { VCOPY(kidOrg, org); OOC_OCTORIGIN(kidOrg, k, kidSize); for (d = node -> leaf.num [k]; d; d--, dataIdx++) { #ifdef DEBUG_OOC_CHECK fprintf(stderr, "OOC_Check: checking record %lu\n", (unsigned long)dataIdx ); #endif if (OOC_GetData(oct, dataIdx, rec)) { fprintf(stderr, "OOC_Check: failed getting record at %lu\n", (unsigned long)dataIdx ); return -1; } key = oct -> key(rec); if (!OOC_InBBox(oct, kidOrg, kidSize, key)) { fprintf(stderr, "OOC_Check: key [%f, %f, %f] (zIdx %020lu) " "outside bbox [[%f-%f], [%f-%f], [%f-%f]] " "in octant %d (should be %d)\n", key [0], key [1], key [2], OOC_KEY2MORTON(key, oct), kidOrg [0], kidOrg [0] + kidSize, kidOrg [1], kidOrg [1] + kidSize, kidOrg [2], kidOrg [2] + kidSize, k, OOC_Branch(oct, org, kidSize, key, NULL) ); /* return -1; */ } } } } else { /* Internal node; recurse on all kids */ const OOC_Node *kid = OOC_KID1(oct, node); OOC_DataIdx numData = 0; if (node -> node.kid >= oct -> numNodes) { fprintf(stderr, "OOC_Check: invalid node index %u\n", node -> node.kid ); return -1; } if (!node -> node.num) { fputs("OOC_Check: empty octree node\n", stderr); return -1; } for (k = 0; k < 8; k++) if (OOC_ISBRANCH(node, k)) { VCOPY(kidOrg, org); OOC_OCTORIGIN(kidOrg, k, kidSize); if (OOC_Check(oct, kid, kidOrg, kidSize, dataIdx)) return -1; dataIdx += OOC_DATACNT(kid); numData += OOC_DATACNT(kid); kid++; } if (node -> node.num != numData) { fprintf(stderr, "Parent/kid data count mismatch; expected %u, found %u\n", node -> node.num, numData ); return -1; } } return 0; } #endif int OOC_SaveOctree (const OOC_Octree *oct, FILE *out) /* Appends octree nodes to specified file along with metadata. Uses * RADIANCE's portable I/O routines. Returns 0 on success, else -1. Note * the internal representation of the nodes is platform specific as they * contain unions, therefore we use the portable I/O routines from the * RADIANCE library */ { int i; OOC_NodeIdx nc; OOC_Node *np = NULL; if (!oct || !oct->nodes || !oct->numData || !oct->numNodes || !out) { fputs("OOC_SaveOctree: empty octree\n", stderr); return -1; } /* Write octree origin, size, data record size, number of records, and * number of nodes */ for (i = 0; i < 3; i++) putflt(oct -> org [i], out); putflt(oct -> size, out); putint(oct -> recSize, sizeof(oct -> recSize), out); putint(oct -> numData, sizeof(oct -> numData), out); putint(oct -> numNodes, sizeof(oct -> numNodes), out); /* Write nodes to file */ for (nc = oct -> numNodes, np = oct -> nodes; nc; nc--, np++) { if (OOC_ISLEAF(np)) { /* Write leaf marker */ putint(-1, 1, out); /* Write leaf octant counters */ for (i = 0; i < 8; i++) putint(np -> leaf.num [i], sizeof(np -> leaf.num [i]), out); } else { /* Write node marker */ putint(0, 1, out); /* Write node, rounding field size up to nearest number of bytes */ putint(np -> node.kid, (OOC_NODEIDX_BITS + 7) >> 3, out); putint(np -> node.num, (OOC_DATAIDX_BITS + 7) >> 3, out); putint(np -> node.branch, 1, out); } if (ferror(out)) { fputs("OOC_SaveOctree: failed writing to node file", stderr); return -1; } } fflush(out); return 0; } int OOC_LoadOctree (OOC_Octree *oct, FILE *nodeFile, RREAL *(*key)(const void*), FILE *leafFile ) /* Load octree nodes and metadata from nodeFile and associate with records * in leafFile. Uses RADIANCE's portable I/O routines. Returns 0 on * success, else -1. */ { int i; OOC_NodeIdx nc; OOC_Node *np = NULL; if (!oct || !nodeFile) return -1; OOC_Null(oct); /* Read octree origin, size, record size, #records and #nodes */ for (i = 0; i < 3; i++) oct -> org [i] = getflt(nodeFile); oct -> size = getflt(nodeFile); oct -> bound [0] = oct -> bound [1] = oct -> bound [2] = oct -> size; VADD(oct -> bound, oct -> bound, oct -> org); oct -> mortonScale = MORTON3D_MAX / oct -> size; oct -> recSize = getint(sizeof(oct -> recSize), nodeFile); oct -> numData = getint(sizeof(oct -> numData), nodeFile); oct -> numNodes = getint(sizeof(oct -> numNodes), nodeFile); oct -> key = key; oct -> leafFile = leafFile; if (feof(nodeFile) || !oct -> numData || !oct -> numNodes) { fputs("OOC_LoadOctree: empty octree\n", stderr); return -1; } /* Allocate octree nodes */ if (!(oct -> nodes = calloc(oct -> numNodes, sizeof(OOC_Node)))) { fputs("OOC_LoadOctree: failed octree allocation\n", stderr); return -1; } /* Read nodes from file */ for (nc = oct -> numNodes, np = oct -> nodes; nc; nc--, np++) { OOC_CLEARNODE(np); /* Read node type marker */ if (getint(1, nodeFile)) { /* Read leaf octant counters and set node type */ for (i = 0; i < 8; i++) np -> leaf.num [i] = getint(sizeof(np -> leaf.num [i]), nodeFile ); OOC_SETLEAF(np); } else { /* Read node, rounding field size up to nearest number of bytes */ np -> node.kid = getint((OOC_NODEIDX_BITS + 7) >> 3, nodeFile); np -> node.num = getint((OOC_DATAIDX_BITS + 7) >> 3, nodeFile); np -> node.branch = getint(1, nodeFile); } if (ferror(nodeFile)) { fputs("OOC_LoadOctree: failed reading from node file\n", stderr); return -1; } } return 0; } void OOC_Delete (OOC_Octree *oct) /* Self-destruct */ { free(oct -> nodes); if (oct -> leafFile) fclose(oct -> leafFile); if (oct -> cache) OOC_DeleteCache(oct -> cache); } #endif /* NIX / PMAP_OOC */ diff --git a/oococt.h b/oococt.h index 360f22e..9105ae1 100644 --- a/oococt.h +++ b/oococt.h @@ -1,256 +1,257 @@ /* ====================================================================== Out-of-core octree data structure 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) ====================================================================== $Id: oococt.h,v 2.3 2016/05/17 17:39:47 rschregle Exp $ */ #ifndef _OOC_OCT_H #define _OOC_OCT_H #include "morton.h" #include "ooccache.h" #include #include /* =================================================================== * CONSTS * =================================================================== */ /* Octree index & data counter types */ typedef uint32_t OOC_NodeIdx; typedef uint32_t OOC_DataIdx; typedef uint8_t OOC_OctCnt; /* Octree node field sizes for above */ #define OOC_NODEIDX_BITS 31 #define OOC_DATAIDX_BITS 32 /* Limits for above */ #define OOC_NODEIDX_MAX ((1L << OOC_NODEIDX_BITS) - 1) #define OOC_DATAIDX_MAX ((1L << OOC_DATAIDX_BITS) - 2) #define OOC_OCTCNT_MAX UINT8_MAX /* Reserved values for above */ #define OOC_DATAIDX_ERR (OOC_DATAIDX_MAX + 1) /* Block size for node allocation (8k) */ #define OOC_BLKSIZE (1L << 13) /* =================================================================== * UTILZ * =================================================================== */ /* Test/set leaf flag for node n */ #define OOC_ISLEAF(n) ((n) -> node.type) #define OOC_SETLEAF(n) ((n) -> node.type = 1) /* Test/set branch bit for octant k in node n */ #define OOC_ISBRANCH(n,k) ( \ !OOC_ISLEAF(n) && (n) -> node.branch & 1 << (k) \ ) #define OOC_SETBRANCH(n,k) ((n) -> node.branch |= 1 << (k)) /* Return index to node's 1st kid */ #define OOC_KID1(o,n) (!OOC_ISLEAF(n) \ ? (o)->nodes + (n)->node.kid \ : NULL \ ) /* Get index to last node in octree; -1 indicates empty octree */ #define OOC_ROOTIDX(o) ((o) -> nodes && (o) -> numNodes > 0 \ ? (o) -> numNodes - 1 \ : -1 \ ) /* Get octree root; this is at the *end* of the node array! */ #define OOC_ROOT(o) (OOC_ROOTIDX(o) >= 0 \ ? (o) -> nodes + OOC_ROOTIDX(o) \ : NULL \ ) /* Copy node to octree root */ #define OOC_SETROOT(o,n) (memcpy(OOC_ROOT(o), (n), sizeof(OOC_Node))) /* Set all node fields to 0 */ #define OOC_CLEARNODE(n) (memset((n), 0, sizeof(OOC_Node))) /* Modify origin o for octant k of size s */ #define OOC_OCTORIGIN(o,k,s) ( \ (o) [0] += (s) * ((k) & 1), \ (o) [1] += (s) * ((k) >> 1 & 1), \ (o) [2] += (s) * ((k) >> 2 & 1) \ ) /* Get num records stored in node (also handles leaves and NULL pointer) */ #define OOC_DATACNT(n) ((n) \ ? OOC_ISLEAF(n) \ ? (n) -> leaf.num [0] + (n) -> leaf.num [1] + \ (n) -> leaf.num [2] + (n) -> leaf.num [3] + \ (n) -> leaf.num [4] + (n) -> leaf.num [5] + \ (n) -> leaf.num [6] + (n) -> leaf.num [7] \ : (n) -> node.num \ : 0 \ ) /* Shortcuts for morton code from key and data record */ #define OOC_KEY2MORTON(k, o) Key2Morton3D( \ (k), (o) -> org, (o) -> mortonScale \ ) #define OOC_REC2MORTON(r, o) Key2Morton3D( \ (o) -> key(r), (o) -> org, (o) -> mortonScale \ ) /* =================================================================== * DATA TYPEZ * =================================================================== */ /* Octree node */ typedef struct { union { struct { /* Interior node (node.type = 0) with: node.kid = index to 1st child node in octree array (a.k.a "Number One Son"), immediately followed by its nonzero siblings as indicated by branch node.num = num records stored in this subtree, node.branch = branch bitmask indicating nonzero kid nodes */ char type : 1; OOC_NodeIdx kid : OOC_NODEIDX_BITS; OOC_DataIdx num : OOC_DATAIDX_BITS; uint8_t branch; /* NOTE that we can't make the kid's index relative to parent * (which would be more compact), since we don't know the * parent's index a priori during the bottom-up construction! */ } node; struct { /* Leaf node (leaf.type = node.type = 1 with: leaf.num [k] = num records stored in octant k */ char type : 1; OOC_OctCnt num [8]; } leaf; }; } OOC_Node; /* Top level octree container */ typedef struct { FVECT org, bound; /* Cube origin (min. XYZ), size, and resulting bounds (max. XYZ) */ RREAL size, *(*key)(const void*); /* Callback for data rec coords */ RREAL mortonScale; /* Scale factor for generating Morton codes from coords */ OOC_Node *nodes; /* Pointer to node array */ /* ************************************* * **** NODES ARE STORED IN POSTFIX **** * **** ORDER, I.E. ROOT IS LAST! **** * *************************************/ OOC_DataIdx numData; /* Num records in leaves */ OOC_NodeIdx maxNodes, /* Num currently allocated nodes */ numNodes; /* Num used nodes (<= maxNodes) */ unsigned recSize, /* Size of data record in bytes */ leafMax, /* Max #records per leaf (for build) */ maxDepth; /* Max subdivision depth (for build) */ FILE *leafFile; /* File containing data in leaves */ + char *leafFname; /* Filename for above */ OOC_Cache *cache; /* I/O cache for leafFile */ } OOC_Octree; /* DOAN' FORGET: NODES IN ***POSTFIX*** ORDAH!!! */ /* =================================================================== * FUNKZ * =================================================================== */ /* Init empty octree */ void OOC_Null (OOC_Octree *oct); /* Initialises octree for records of size recSize with 3D coordinates, * accessed via key() callback. The octree's bounding box is defined by * its origin org and size per axis, and must contain all keys. The * octree is associated with the specified leaf file containing the * records in Morton order */ void OOC_Init (OOC_Octree *oct, unsigned recSize, FVECT org, RREAL size, - RREAL *(*key)(const void*), FILE *leafFile + RREAL *(*key)(const void*), FILE *leafFile, char *leafFname ); /* Allocate new octree node, enlarge array if necessary. Returns pointer to new node or NULL if failed. */ OOC_Node *NewNode (OOC_Octree *oct); /* Reads indexed data record from leaf file and copies it to rec, else * returns -1 on failure */ int OOC_GetData (OOC_Octree *oct, OOC_DataIdx idx, void *rec); /* Appends octree nodes to specified file along with metadata. Uses * RADIANCE's portable I/O routines. Returns 0 on success, else -1. */ int OOC_SaveOctree (const OOC_Octree *oct, FILE *nodeFile); /* Load octree nodes and metadata from nodeFile and associate with * records in leafFile. Uses RADIANCE's portable I/O routines. Returns * 0 on success, else -1. */ int OOC_LoadOctree (OOC_Octree *oct, FILE *nodeFile, RREAL *(*key)(const void*), FILE *leafFile ); #ifdef DEBUG_OOC /* Traverse tree & check for valid nodes; oct -> leafFile must be open; * return 0 if ok, otherwise -1 */ int OOC_Check (OOC_Octree *oct, const OOC_Node *node, FVECT org, RREAL size, OOC_DataIdx dataIdx ); #endif /* Check whether key lies within bounding box defined by org and size */ int OOC_InBBox (const OOC_Octree *oct, const FVECT org, RREAL size, const FVECT key ); /* Return index of octant containing key and corresponding origin if * nuOrg = NULL, or -1 if key lies outside all octants. Size is already * assumed to be halved, i.e. corresponding to the length of an octant * axis. */ int OOC_Branch (const OOC_Octree *oct, const FVECT org, RREAL size, const FVECT key, FVECT nuOrg ); /* Optimised version of OOC_Branch() with precomputed key Morton code */ int OOC_BranchZ (const OOC_Octree *oct, const FVECT org, RREAL size, MortonIdx keyZ, FVECT nuOrg ); /* For a given octree node, return the sum of the data counters for kids * [0..k-1]. On return, the node either points to the k'th kid, or * NULL if the kid is nonexistent or the node is a leaf */ OOC_DataIdx OOC_GetKid (const OOC_Octree *oct, OOC_Node **node, unsigned k) ; /* Self-destruct */ void OOC_Delete (OOC_Octree *oct); /* DID WE MENTION NODES ARE IN POSTFIX ORDAH? */ #endif /* _OOC_OCT_H */ diff --git a/pmap.c b/pmap.c index 1042896..5a1a4da 100644 --- a/pmap.c +++ b/pmap.c @@ -1,929 +1,929 @@ #ifndef lint static const char RCSid[] = "$Id: pmap.c,v 2.18 2021/02/19 02:10:35 rschregle Exp $"; #endif /* ====================================================================== Photon map main module Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Fraunhofer Institute for Solar Energy Systems, supported by the German Research Foundation (DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme FARESYS") (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ====================================================================== $Id: pmap.c,v 2.18 2021/02/19 02:10:35 rschregle Exp $ */ #include "pmap.h" #include "pmapmat.h" #include "pmapsrc.h" #include "pmaprand.h" #include "pmaproi.h" #include "pmapio.h" #include "pmapdens.h" #include "pmapbias.h" #include "pmapdiag.h" #include "pmutil.h" #include "otypes.h" #include "otspecial.h" #include #if NIX #include #include #include #include #endif static int photonParticipate (RAY *ray) /* Trace photon through participating medium. Returns 1 if passed through, or 0 if absorbed and $*%&ed. Analogon to rayparticipate(). */ { int i; RREAL xi1, cosTheta, phi, du, dv, rmax0; const RREAL cext = colorAvg(ray -> cext), invCext = 1 / cext, albedo = colorAvg(ray -> albedo), invAlbedo = 1 / albedo, gecc = ray -> gecc, gecc2 = sqr(gecc); FVECT u, v; COLOR cvext; /* Save incoming ray's max distance; this is required since rmax is set to the mean free distance for each path segment (see below) */ rmax0 = ray -> rmax; /* Mean free distance until interaction with medium */ ray -> rmax = -log(pmapRandom(mediumState)) / cext; while (!localhit(ray, &thescene)) { /* No geometry hit within mean free dist ray->rmax; photon interacts with medium */ if (!incube(&thescene, ray -> rop)) { /* Terminate photon if it has leaked from the scene */ #ifdef DEBUG_PMAP_LEAKED fprintf( stderr, "Volume photon leaked from scene at [%.3f %.3f %.3f]\n", ray -> rop [0], ray -> rop [1], ray -> rop [2] ); #endif return 0; } /* RGB attenuation over mean free distance */ setcolor(cvext, exp(-ray -> rmax * ray -> cext [0]), exp(-ray -> rmax * ray -> cext [1]), exp(-ray -> rmax * ray -> cext [2]) ); /* Modify ray color and normalise */ multcolor(ray -> rcol, cvext); colorNorm(ray -> rcol); /* Set ray origin for next interaction */ VCOPY(ray -> rorg, ray -> rop); /* Update ray's path length; terminate photon if maximum reached. See also photonRay(). */ rmax0 -= ray -> rot; if (photonMaxPathDist > 0 && rmax0 < 0) return 0; /* Store volume photon if nonzero albedo; this also accounts for direct inscattering from light sources */ if (albedo > FTINY) { /* Set ray's path length to store with photon */ ray -> rmax = rmax0; #ifdef PMAP_PHOTONFLOW if (lightFlowPhotonMapping) { /* Temporarily correct normalised flux for lightflow photons based on extinction (= mean num photons deposited per unit path length). */ scalecolor(ray -> rcol, invCext); newPhoton(lightFlowPmap, ray); /* Undo correction for next iteration after storing photon */ scalecolor(ray -> rcol, cext); } else #endif newPhoton(volumePmap, ray); } /* Colour bleeding without attenuation (?) */ multcolor(ray -> rcol, ray -> albedo); scalecolor(ray -> rcol, invAlbedo); if (pmapRandom(rouletteState) < albedo) { /* Scatter photon */ xi1 = pmapRandom(scatterState); cosTheta = ray -> gecc <= FTINY ? 2 * xi1 - 1 : 0.5 / gecc * ( 1 + gecc2 - sqr((1 - gecc2) / (1 + gecc * (2 * xi1 - 1))) ); phi = 2 * PI * pmapRandom(scatterState); du = dv = sqrt(1 - sqr(cosTheta)); /* sin(theta) */ du *= cos(phi); dv *= sin(phi); /* Get axes u & v perpendicular to photon direction */ i = 0; do { v [0] = v [1] = v [2] = 0; v [i++] = 1; fcross(u, v, ray -> rdir); } while (normalize(u) < FTINY); fcross(v, ray -> rdir, u); for (i = 0; i < 3; i++) ray -> rdir [i] = du * u [i] + dv * v [i] + cosTheta * ray -> rdir [i]; } ray -> rlvl++; ray -> rmax = -log(pmapRandom(mediumState)) / cext; /* fprintf(stderr, "%.3f\n", ray -> rmax); */ } /* Passed through medium until intersecting local object */ setcolor(cvext, exp(-ray -> rot * ray -> cext [0]), exp(-ray -> rot * ray -> cext [1]), exp(-ray -> rot * ray -> cext [2]) ); /* Modify ray color and normalise */ multcolor(ray -> rcol, cvext); colorNorm(ray -> rcol); /* Restore outgoing ray's max distance */ ray -> rmax = rmax0; return 1; } void tracePhoton (RAY *ray) /* Follow photon as it bounces around... */ { long mod; OBJREC *mat, *port = NULL; if (!ray -> parent) { /* !!! PHOTON PORT REJECTION SAMPLING HACK: get photon port for * !!! primary ray from ray -> ro, then reset the latter to NULL so * !!! as not to interfere with localhit() */ port = ray -> ro; ray -> ro = NULL; } if (ray -> rlvl > photonMaxBounce) { #ifdef PMAP_RUNAWAY_WARN error(WARNING, "runaway photon!"); #endif return; } if (colorAvg(ray -> cext) > FTINY && !photonParticipate(ray)) return; /* NOTE: localhit() limits intersections to the aft plane's distance ray->rmax. This mechanism is (ab)used here to limit the photon path length by initialising ray->rmax to photonMaxPathDist; see emitPhoton(). With unlimited path length (photonMaxPathDist=0), ray->rmax becomes negative (see photonRay()), and the aft plane has no effect. */ if (localhit(ray, &thescene)) { mod = ray -> ro -> omod; if (port && ray -> ro != port) { /* Terminate photon if emitted from port without intersecting it */ #ifdef PMAP_PORTREJECT_WARN sprintf(errmsg, "photon outside port %s", ray -> ro -> oname); error(WARNING, errmsg); #endif return; } if ((ray -> clipset && inset(ray -> clipset, mod)) || mod == OVOID) { /* Transfer ray if modifier is VOID or clipped within antimatta */ RAY tray; photonRay(ray, &tray, PMAP_XFER, NULL); tracePhoton(&tray); } else { /* Scatter for modifier material */ mat = objptr(mod); photonScatter [mat -> otype] (mat, ray); } } } static void preComputeGlobal (PhotonMap *pmap) /* Precompute irradiance from global photons for final gathering for a random subset of finalGather * pmap -> numPhotons photons, and builds the photon map, discarding the original photons. */ /* !!! NOTE: PRECOMPUTATION WITH OOC CURRENTLY WITHOUT CACHE !!! */ { unsigned long i, numPreComp; unsigned j; PhotonIdx pIdx; Photon photon; RAY ray; PhotonMap nuPmap; repComplete = numPreComp = finalGather * pmap -> numPhotons; if (verbose) { sprintf(errmsg, - "\nPrecomputing irradiance for %ld global photons\n", - numPreComp); + "\nPrecomputing irradiance for %ld global photons\n", numPreComp + ); eputs(errmsg); #if NIX fflush(stderr); #endif } /* Copy photon map for precomputed photons */ memcpy(&nuPmap, pmap, sizeof(PhotonMap)); /* Zero counters, init new heap and extents */ nuPmap.numPhotons = 0; initPhotonHeap(&nuPmap); for (j = 0; j < 3; j++) { nuPmap.minPos [j] = FHUGE; nuPmap.maxPos [j] = -FHUGE; } /* Record start time, baby */ repStartTime = time(NULL); #ifdef SIGCONT signal(SIGCONT, pmapPreCompReport); #endif repProgress = 0; photonRay(NULL, &ray, PRIMARY, NULL); ray.ro = NULL; for (i = 0; i < numPreComp; i++) { /* Get random photon from stratified distribution in source heap to * avoid duplicates and clustering */ pIdx = firstPhoton(pmap) + (unsigned long)((i + pmapRandom(pmap -> randState)) / finalGather); getPhoton(pmap, pIdx, &photon); /* Init dummy photon ray with intersection at photon position */ VCOPY(ray.rop, photon.pos); for (j = 0; j < 3; j++) ray.ron [j] = photon.norm [j] / 127.0; /* Get density estimate at photon position */ photonDensity(pmap, &ray, ray.rcol); /* Append photon to new heap from ray */ newPhoton(&nuPmap, &ray); /* Update progress */ repProgress++; if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) pmapPreCompReport(); #ifdef SIGCONT else signal(SIGCONT, pmapPreCompReport); #endif } /* Flush heap */ flushPhotonHeap(&nuPmap); #ifdef SIGCONT signal(SIGCONT, SIG_DFL); #endif /* Trash original pmap, replace with precomputed one */ deletePhotons(pmap); memcpy(pmap, &nuPmap, sizeof(PhotonMap)); if (verbose) { eputs("\nRebuilding precomputed photon map\n"); #if NIX fflush(stderr); #endif } /* Rebuild underlying data structure, destroying heap */ buildPhotonMap(pmap, NULL, NULL, 1); } typedef struct { unsigned long numPhotons [NUM_PMAP_TYPES], numEmitted, numComplete; } PhotonCnt; /* Photon distribution counter update interval expressed as bitmask; counters shared among subling subprocesses will only be updated in multiples of PMAP_CNTUPDATE in order to reduce contention */ #define PMAP_CNTUPDATE 0xffL void distribPhotons (PhotonMap **pmaps, unsigned numProc) { EmissionMap emap; char errmsg2 [128], shmFname [PMAP_TMPFNLEN]; unsigned t, srcIdx, proc; double totalFlux = 0; int shmFile, stat, pid; PhotonMap *pm; PhotonCnt *photonCnt, lastPhotonCnt; pid_t procPids [PMAP_MAXPROC]; for (t = 0; t < NUM_PMAP_TYPES && !pmaps [t]; t++); if (t >= NUM_PMAP_TYPES) error(USER, "no photon maps defined in distribPhotons"); if (!nsources) error(USER, "no light sources in distribPhotons"); /* =================================================================== * INITIALISATION - Set up emission and scattering funcs * =================================================================== */ emap.samples = NULL; emap.maxPartitions = MAXSPART; emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1); if (!emap.partitions) error(INTERNAL, "can't allocate source partitions in distribPhotons"); /* Initialise all defined photon maps */ for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t]) { initPhotonMap(pmaps [t], t); /* Open photon heapfile */ initPhotonHeap(pmaps [t]); /* Per-subprocess target count */ pmaps [t] -> distribTarget /= numProc; if (!pmaps [t] -> distribTarget) error(INTERNAL, "no photons to distribute in distribPhotons"); } initPhotonEmissionFuncs(); initPhotonScatterFuncs(); /* Get photon ports from modifier list */ getPhotonPorts(photonPortList); /* Get photon sensor modifiers */ getPhotonSensors(photonSensorList); /* Get polyhedral regions of interest */ getPolyROIs(pmapROImodList); #if NIX /* Set up shared mem for photon counters (zeroed by ftruncate) */ strcpy(shmFname, PMAP_TMPFNAME); shmFile = mkstemp(shmFname); if (shmFile < 0 || ftruncate(shmFile, sizeof(*photonCnt)) < 0) error(SYSTEM, "failed shared mem init in distribPhotons"); photonCnt = mmap(NULL, sizeof(*photonCnt), PROT_READ | PROT_WRITE, MAP_SHARED, shmFile, 0 ); if (photonCnt == MAP_FAILED) error(SYSTEM, "failed mapping shared memory in distribPhotons"); #else /* Allocate photon counters statically on Windoze */ if (!(photonCnt = malloc(sizeof(PhotonCnt)))) error(SYSTEM, "failed trivial malloc in distribPhotons"); photonCnt -> numEmitted = photonCnt -> numComplete = 0; #endif /* NIX */ /* Zero overflow tracking counters */ memset(&lastPhotonCnt, 0, sizeof(lastPhotonCnt)); if (verbose) { sprintf(errmsg, "\nIntegrating flux from %d sources", nsources); if (photonPorts) { sprintf(errmsg2, " via %d ports", numPhotonPorts); strcat(errmsg, errmsg2); } strcat(errmsg, "\n"); eputs(errmsg); } /* =================================================================== * FLUX INTEGRATION - Get total photon flux from light sources * =================================================================== */ for (srcIdx = 0; srcIdx < nsources; srcIdx++) { unsigned portCnt = 0; emap.src = source + srcIdx; do { /* Set next photon port if defined; note this loop iterates once if * no ports are defined. */ emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt : NULL; photonPartition [emap.src -> so -> otype] (&emap); if (verbose) { sprintf(errmsg, "\tIntegrating flux from source %s ", source [srcIdx].so -> oname ); if (emap.port) { sprintf(errmsg2, "via port %s ", photonPorts [portCnt].so -> oname ); strcat(errmsg, errmsg2); } sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions); strcat(errmsg, errmsg2); eputs(errmsg); #if NIX fflush(stderr); #endif } for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions; emap.partitionCnt++ ) { initPhotonEmission(&emap, pdfSamples); totalFlux += colorAvg(emap.partFlux); } portCnt++; } while (portCnt < numPhotonPorts); } if (totalFlux < FTINY) error(USER, "zero flux from light sources"); /* Record start time for progress reports */ repStartTime = time(NULL); if (verbose) { sprintf(errmsg, "\nPhoton distribution @ %d procs\n", numProc); eputs(errmsg); } /* MAIN LOOP */ for (proc = 0; proc < numProc; proc++) { #if NIX if (!(pid = fork())) { /* SUBPROCESS ENTERS HERE; open and mmapped files inherited */ #else if (1) { /* No subprocess under Windoze */ #endif /* Local photon counters for this subprocess */ unsigned passCnt = 0, prePassCnt = 0; unsigned long lastNumPhotons [NUM_PMAP_TYPES]; unsigned long localNumEmitted = 0; /* Num photons emitted by this subprocess alone */ /* Seed RNGs from PID for decorellated photon distribution */ pmapSeed(randSeed + proc, partState); pmapSeed(randSeed + (proc + 1) % numProc, emitState); pmapSeed(randSeed + (proc + 2) % numProc, cntState); pmapSeed(randSeed + (proc + 3) % numProc, mediumState); pmapSeed(randSeed + (proc + 4) % numProc, scatterState); pmapSeed(randSeed + (proc + 5) % numProc, rouletteState); #ifdef DEBUG_PMAP /* Output child process PID after random delay to prevent corrupted * console output due to race condition */ usleep(1e6 * pmapRandom(rouletteState)); fprintf(stderr, "Proc %d: PID = %d (waiting 10 sec to attach debugger...)\n", proc, getpid() ); /* Allow time for debugger to attach to child process */ sleep(10); #endif for (t = 0; t < NUM_PMAP_TYPES; t++) lastNumPhotons [t] = 0; /* ============================================================= * 2-PASS PHOTON DISTRIBUTION * Pass 1 (pre): emit fraction of target photon count * Pass 2 (main): based on outcome of pass 1, estimate remaining * number of photons to emit to approximate target * count * ============================================================= */ do { double numEmit; if (!passCnt) { /* INIT PASS 1 */ /* Skip if no photons contributed after sufficient * iterations; make it clear to user which photon maps are * missing so (s)he can check geometry and materials */ if (++prePassCnt > maxPreDistrib) { sprintf(errmsg, "proc %d: too many prepasses", proc); for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t] && !pmaps [t] -> numPhotons) { sprintf(errmsg2, ", no %s photons stored", pmapName [t] ); strcat(errmsg, errmsg2); } sprintf(errmsg2, "; increase -apD"); strcat(errmsg, errmsg2); error(USER, errmsg); break; } /* Num to emit is fraction of minimum target count */ numEmit = FHUGE; for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t]) numEmit = min(pmaps [t] -> distribTarget, numEmit); numEmit *= preDistrib; } else { /* INIT PASS 2 */ /* Based on the outcome of the predistribution we can now * estimate how many more photons we have to emit for each * photon map to meet its respective target count. This * value is clamped to 0 in case the target has already been * exceeded in the pass 1. */ double maxDistribRatio = 0; /* Set the distribution ratio for each map; this indicates * how many photons of each respective type are stored per * emitted photon, and is used as probability for storing a * photon by newPhoton(). Since this biases the photon * density, newPhoton() promotes the flux of stored photons * to compensate. */ for (t = 0; t < NUM_PMAP_TYPES; t++) if ((pm = pmaps [t])) { pm -> distribRatio = (double)pm -> distribTarget / pm -> numPhotons - 1; /* Check if photon map "overflowed", i.e. exceeded its * target count in the prepass; correcting the photon * flux via the distribution ratio is no longer * possible, as no more photons of this type will be * stored, so notify the user rather than deliver * incorrect results. In future we should handle this * more intelligently by using the photonFlux in each * photon map to individually correct the flux after * distribution. */ if (pm -> distribRatio <= FTINY) { sprintf(errmsg, "%s photon map overflow in prepass, reduce -apD ", pmapName [t] ); error(INTERNAL, errmsg); } maxDistribRatio = max(pm -> distribRatio, maxDistribRatio ); } /* Normalise distribution ratios and calculate number of * photons to emit in main pass */ for (t = 0; t < NUM_PMAP_TYPES; t++) if ((pm = pmaps [t])) pm -> distribRatio /= maxDistribRatio; if ((numEmit = localNumEmitted * maxDistribRatio) < FTINY) /* No photons left to distribute in main pass */ break; } /* PHOTON DISTRIBUTION LOOP */ for (srcIdx = 0; srcIdx < nsources; srcIdx++) { unsigned portCnt = 0; emap.src = source + srcIdx; do { /* Set next photon port if defined; note this loop iterates * once if no ports are defined. */ emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt : NULL; photonPartition [emap.src -> so -> otype] (&emap); if (verbose && !proc) { /* Output from subproc 0 only to avoid race condition * on console I/O */ if (!passCnt) sprintf(errmsg, "\tPREPASS %d on source %s ", prePassCnt, source [srcIdx].so -> oname ); else sprintf(errmsg, "\tMAIN PASS on source %s ", source [srcIdx].so -> oname ); if (emap.port) { sprintf(errmsg2, "via port %s ", photonPorts [portCnt].so -> oname ); strcat(errmsg, errmsg2); } sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions ); strcat(errmsg, errmsg2); eputs(errmsg); #if NIX fflush(stderr); #endif } for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions; emap.partitionCnt++ ) { double partNumEmit; unsigned long partEmitCnt; /* Get photon origin within current source partishunn * and build emission map */ photonOrigin [emap.src -> so -> otype] (&emap); initPhotonEmission(&emap, pdfSamples); /* Number of photons to emit from ziss partishunn -- * proportional to flux; photon ray weight and scalar * flux are uniform (latter only varying in RGB). */ partNumEmit = numEmit * colorAvg(emap.partFlux) / totalFlux; partEmitCnt = (unsigned long)partNumEmit; /* Probabilistically account for fractional photons */ if (pmapRandom(cntState) < partNumEmit - partEmitCnt) partEmitCnt++; /* Integer counter avoids FP rounding errors during * iteration */ while (partEmitCnt--) { RAY photonRay; /* Emit photon based on PDF, skip if rejected. */ /* NOTE: rejection sampling skips incrementing the * emission counter (see below), thus compensating * for the rejected photons by increasing the photon * flux in proportion to the lower effective * emission count. * BUG: THIS INTERFERES WITH THE PROGRESS COUNTER * REPORTED TO THE PARENT, AND WITH THE * PREDISTRIBUTION PASS --> PHOTON DISTRIBUTION WILL * FINISH EARLY, WITH FEWER PHOTONS THAN TARGETTED! */ if (!emitPhoton(&emap, &photonRay)) continue; /* Set subprocess ID in photonRay; extends path ID. Used by newPhoton() as photon attributes. */ PMAP_SETRAYPROC(&photonRay, proc); tracePhoton(&photonRay); /* Update local emission counter*/ localNumEmitted++; if (!(localNumEmitted & PMAP_CNTUPDATE)) { /* Update global counters shared with siblings in intervals to reduce overhead and contention */ lastPhotonCnt.numEmitted = photonCnt -> numEmitted; photonCnt -> numEmitted += PMAP_CNTUPDATE; lastPhotonCnt.numComplete = photonCnt -> numComplete; photonCnt -> numComplete += PMAP_CNTUPDATE; /* Check for overflow using previous values */ if (photonCnt -> numEmitted < lastPhotonCnt.numEmitted || photonCnt -> numComplete < lastPhotonCnt.numComplete ) { sprintf(errmsg, "photon emission counter " "overflow (was: %ld, is: %ld)", lastPhotonCnt.numEmitted, photonCnt -> numEmitted ); error(INTERNAL, errmsg); } /* Update global photon count for each pmap */ for (t = 0; t < NUM_PMAP_TYPES; t++) { if (pmaps [t]) { lastPhotonCnt.numPhotons [t] = photonCnt -> numPhotons [t]; /* Differentially increment using local counter lastNumPhotons */ photonCnt -> numPhotons [t] += pmaps [t] -> numPhotons - lastNumPhotons [t]; lastNumPhotons [t] = pmaps [t] -> numPhotons; /* Check for photon counter overflow (this * could only happen before an emission * counter overflow if the scene has an * absurdly high albedo and/or very dense * geometry) */ if (photonCnt -> numPhotons [t] < lastPhotonCnt.numPhotons [t] ) { sprintf(errmsg, "photon counter " "overflow (was: %ld, is: %ld)", lastPhotonCnt.numPhotons [t], photonCnt -> numPhotons [t] ); error(INTERNAL, errmsg); } } } } } #if !NIX /* Synchronous progress report on Windoze */ if (!proc && photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) { repEmitted = repProgress = photonCnt -> numEmitted; repComplete = photonCnt -> numComplete; pmapDistribReport(); } #endif } portCnt++; } while (portCnt < numPhotonPorts); } for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t] && !pmaps [t] -> numPhotons) { /* Double preDistrib in case a photon map is empty and * redo pass 1 --> possibility of infinite loop for * pathological scenes (e.g. absorbing materials) */ preDistrib *= 2; break; } if (t >= NUM_PMAP_TYPES) /* No empty photon maps found; now do pass 2 */ passCnt++; } while (passCnt < 2); /* Flush heap buffa for every pmap one final time; * avoids potential data corruption! */ for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t]) { flushPhotonHeap(pmaps [t]); /* Heap file closed automatically on exit fclose(pmaps [t] -> heap); */ #ifdef DEBUG_PMAP sprintf(errmsg, "Proc %d: total %ld photons\n", proc, pmaps [t] -> numPhotons ); eputs(errmsg); #endif } #if NIX /* Terminate subprocess */ exit(0); #endif } else { /* PARENT PROCESS CONTINUES LOOP ITERATION HERE */ if (pid < 0) error(SYSTEM, "failed to fork subprocess in distribPhotons"); else /* Save child process IDs */ procPids [proc] = pid; } } #if NIX /* PARENT PROCESS CONTINUES HERE */ #ifdef SIGCONT /* Enable progress report signal handler */ signal(SIGCONT, pmapDistribReport); #endif /* Wait for subprocesses to complete while reporting progress */ proc = numProc; while (proc) { while (waitpid(-1, &stat, WNOHANG) > 0) { /* Subprocess exited; check status */ if (!WIFEXITED(stat) || WEXITSTATUS(stat)) { /* Exited with error; terminate siblings, clean up & bail out */ for (proc = 0; proc < numProc; proc++) kill(procPids [proc], SIGKILL); /* Unmap shared memory, squish mapped file */ munmap(photonCnt, sizeof(*photonCnt)); close(shmFile); unlink(shmFname); error(USER, "failed photon distribution"); } --proc; } /* Nod off for a bit and update progress */ sleep(1); /* Asynchronous progress report from shared subprocess counters */ repEmitted = repProgress = photonCnt -> numEmitted; repComplete = photonCnt -> numComplete; repProgress = repComplete = 0; for (t = 0; t < NUM_PMAP_TYPES; t++) if ((pm = pmaps [t])) { /* Get global photon count from shmem updated by subprocs */ repProgress += pm -> numPhotons = photonCnt -> numPhotons [t]; repComplete += pm -> distribTarget; } repComplete *= numProc; if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) pmapDistribReport(); #ifdef SIGCONT else signal(SIGCONT, pmapDistribReport); #endif } #endif /* NIX */ /* =================================================================== * POST-DISTRIBUTION - Set photon flux and build data struct for photon * storage, etc. * =================================================================== */ #ifdef SIGCONT /* Reset signal handler */ signal(SIGCONT, SIG_DFL); #endif free(emap.samples); /* Set photon flux */ totalFlux /= photonCnt -> numEmitted; #if NIX /* Photon counters no longer needed, unmap shared memory */ munmap(photonCnt, sizeof(*photonCnt)); close(shmFile); unlink(shmFname); #else free(photonCnt); #endif if (verbose) eputs("\n"); for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t]) { if (verbose) { sprintf(errmsg, "Building %s photon map\n", pmapName [t]); eputs(errmsg); #if NIX fflush(stderr); #endif } /* Build underlying data structure; heap is destroyed */ buildPhotonMap(pmaps [t], &totalFlux, NULL, numProc); } /* Precompute photon irradiance if necessary */ if (preCompPmap) { if (verbose) eputs("\n"); preComputeGlobal(preCompPmap); } if (verbose) eputs("\n"); } diff --git a/pmapbias.c b/pmapbias.c index 595efc4..6d8e2ee 100644 --- a/pmapbias.c +++ b/pmapbias.c @@ -1,438 +1,442 @@ #ifndef lint static const char RCSid[] = "$Id: pmapbias.c,v 2.5 2016/05/17 17:39:47 rschregle Exp $"; #endif /* ================================================================== Bias compensation for photon density estimates For background see: R. Schregle, "Bias Compensation for Photon Maps", Computer Graphics Forum, v22:n4, pp. 729-742, Dec. 2003. Roland Schregle (roland.schregle@gmail.com) (c) Fraunhofer Institute for Solar Energy Systems ================================================================== $Id: pmapbias.c,v 2.5 2016/05/17 17:39:47 rschregle Exp $ */ #include "pmapbias.h" #include "pmap.h" #include "pmaprand.h" void squeuePartition (PhotonSearchQueueNode* squeue, unsigned lo, unsigned mid, unsigned hi) /* REVERSE Partition squeue such that all photons in squeue-hi+1..squeue-mid are farther than the median at squeue-mid+1, and those in squeue-mid+2..squeue-lo+1 are closer than the median. This means that squeue points to the END of the queue, and the (1-based) indices are offsets relative to it. This convoluted scheme is adopted since the queue is initially a maxheap, so reverse sorting is expected to be faster. */ { unsigned l, h, p; PhotonSearchQueueNode *lp, *hp, *pp, tmp; float pivot; while (hi > lo) { /* Grab pivot node in middle as an educated guess, since our queue is sorta sorted. */ l = lo; h = hi; p = mid; lp = squeue - lo + 1; hp = squeue - hi + 1; pp = squeue - p + 1; pivot = pp -> dist2; /* l & h converge, swapping elements out of order with respect to pivot node. */ while (l < h) { while (lp -> dist2 <= pivot && l <= h && l < hi) ++l, --lp; while (hp -> dist2 >= pivot && h >= l && h > lo) --h, ++hp; if (l < h) { /* Swap */ tmp.idx = lp -> idx; tmp.dist2 = lp -> dist2; lp -> idx = hp -> idx; lp -> dist2 = hp -> dist2; hp -> idx = tmp.idx; hp -> dist2 = tmp.dist2; } } /* Swap convergence and pivot node */ if (p > h) { /* Need this otherwise shit happens! Since lp -> dist2 > hp -> dist2, we swap either l or p depending on whether we're above or below p */ h = l; hp = lp; } tmp.idx = hp -> idx; tmp.dist2 = hp -> dist2; hp -> idx = pp -> idx; hp -> dist2 = pivot; pp -> idx = tmp.idx; pp -> dist2 = tmp.dist2; if (h >= mid) hi = h - 1; if (h <= mid) lo = h + 1; } /* Once lo & hi have converged, we have found the median! */ } void biasComp (PhotonMap* pmap, COLOR irrad) /* Photon density estimate with bias compensation -- czech dis shit out! */ { unsigned i, numLo, numHi, numMid; float r, totalWeight = 0; COLOR fluxLo, fluxMid, irradVar, irradAvg, p, d; Photon *photon; PhotonSearchQueueNode *sqn, *sqEnd; PhotonBiasCompNode *hist, *histEnd; if (!pmap -> biasCompHist) { /* Allocate bias compensation history */ numHi = pmap -> maxGather - pmap -> minGather; for (i = pmap -> minGather + 1; numHi > 1; numHi >>= 1, ++i); pmap -> biasCompHist = calloc(i, sizeof(PhotonBiasCompNode)); if (!pmap -> biasCompHist) error(USER, "can't allocate bias compensation history"); } numLo = min(pmap -> minGather, pmap -> squeue.tail - 1); numHi = min(pmap -> maxGather, pmap -> squeue.tail - 1); sqn = sqEnd = pmap -> squeue.node + pmap -> squeue.tail - 1; - histEnd = pmap -> biasCompHist; + histEnd = (PhotonBiasCompNode*)pmap -> biasCompHist; setcolor(fluxLo, 0, 0, 0); /* Partition to get numLo closest photons starting from END of queue */ squeuePartition(sqEnd, 1, numLo + 1, numHi); /* Get irradiance estimates (ignoring 1 / PI) using 1..numLo photons and chuck in history. Queue is traversed BACKWARDS. */ for (i = 1; i <= numLo; i++, sqn--) { /* Make sure furthest two photons are consecutive w.r.t. distance */ squeuePartition(sqEnd, i, i + 1, numLo + 1); photon = getNearestPhoton(&pmap -> squeue, sqn -> idx); getPhotonFlux(photon, irrad); addcolor(fluxLo, irrad); /* Average radius between furthest two photons to improve accuracy */ r = 0.25 * (sqn -> dist2 + (sqn - 1) -> dist2 + 2 * sqrt(sqn -> dist2 * (sqn - 1) -> dist2)); /* Add irradiance and weight to history. Weights should increase monotonically based on number of photons used for the estimate. */ histEnd -> irrad [0] = fluxLo [0] / r; histEnd -> irrad [1] = fluxLo [1] / r; histEnd -> irrad [2] = fluxLo [2] / r; totalWeight += histEnd++ -> weight = BIASCOMP_WGT((float)i); } /* Compute expected value (average) and variance of irradiance based on history for numLo photons. */ setcolor(irradAvg, 0, 0, 0); setcolor(irradVar, 0, 0, 0); - for (hist = pmap -> biasCompHist; hist < histEnd; ++hist) + for (hist = (PhotonBiasCompNode*)pmap -> biasCompHist; + hist < histEnd; ++hist + ) for (i = 0; i <= 2; ++i) { irradAvg [i] += r = hist -> weight * hist -> irrad [i]; irradVar [i] += r * hist -> irrad [i]; } for (i = 0; i <= 2; ++i) { r = irradAvg [i] /= totalWeight; irradVar [i] = irradVar [i] / totalWeight - r * r; } /* Do binary search within interval [numLo, numHi]. numLo is towards the END of the queue. */ while (numHi - numLo > 1) { numMid = (numLo + numHi) >> 1; /* Split interval to get numMid closest photons starting from the END of the queue */ squeuePartition(sqEnd, numLo, numMid, numHi); /* Make sure furthest two photons are consecutive wrt distance */ squeuePartition(sqEnd, numMid, numMid + 1, numHi); copycolor(fluxMid, fluxLo); sqn = sqEnd - numLo; /* Get irradiance for numMid photons (ignoring 1 / PI) */ for (i = numLo; i < numMid; i++, sqn--) { photon = getNearestPhoton(&pmap -> squeue, sqn -> idx); getPhotonFlux(photon, irrad); addcolor(fluxMid, irrad); } /* Average radius between furthest two photons to improve accuracy */ r = 0.25 * (sqn -> dist2 + (sqn + 1) -> dist2 + 2 * sqrt(sqn -> dist2 * (sqn + 1) -> dist2)); /* Get deviation from current average, and probability that it's due to noise from gaussian distribution based on current variance. Since we are doing this for each colour channel we can also detect chromatic bias. */ for (i = 0; i <= 2; ++i) { d [i] = irradAvg [i] - (irrad [i] = fluxMid [i] / r); p [i] = exp(-0.5 * d [i] * d [i] / irradVar [i]); } if (pmapRandom(pmap -> randState) < colorAvg(p)) { /* Deviation is probably noise, so add mid irradiance to history */ copycolor(histEnd -> irrad, irrad); totalWeight += histEnd++ -> weight = BIASCOMP_WGT((float)numMid); setcolor(irradAvg, 0, 0, 0); setcolor(irradVar, 0, 0, 0); /* Update average and variance */ - for (hist = pmap -> biasCompHist; hist < histEnd; ++hist) + for (hist = (PhotonBiasCompNode*)pmap -> biasCompHist; + hist < histEnd; ++hist + ) for (i = 0; i <= 2; i++) { r = hist -> irrad [i]; irradAvg [i] += hist -> weight * r; irradVar [i] += hist -> weight * r * r; } for (i = 0; i <= 2; i++) { r = irradAvg [i] /= totalWeight; irradVar [i] = irradVar [i] / totalWeight - r * r; } /* Need more photons --> recurse on [numMid, numHi] */ numLo = numMid; copycolor(fluxLo, fluxMid); } else /* Deviation is probably bias --> need fewer photons, so recurse on [numLo, numMid] */ numHi = numMid; } --histEnd; for (i = 0; i <= 2; i++) { /* Estimated relative error */ d [i] = histEnd -> irrad [i] / irradAvg [i] - 1; #ifdef BIASCOMP_BWIDTH /* Return bandwidth instead of irradiance */ irrad [i] = numHi / PI; #else /* 1 / PI required as ambient normalisation factor */ irrad [i] = histEnd -> irrad [i] / (PI * PI); #endif } /* Update statistix */ r = colorAvg(d); pmap -> rmsError += r * r; if (r < pmap -> minError) pmap -> minError = r; if (r > pmap -> maxError) pmap -> maxError = r; if (numHi < pmap -> minGathered) pmap -> minGathered = numHi; if (numHi > pmap -> maxGathered) pmap -> maxGathered = numHi; pmap -> totalGathered += numHi; ++pmap -> numDensity; } /* Volume bias compensation disabled (probably redundant) */ #if 0 void volumeBiasComp (PhotonMap* pmap, const RAY* ray, COLOR irrad) /* Photon volume density estimate with bias compensation -- czech dis shit out! */ { unsigned i, numLo, numHi, numMid = 0; PhotonSQNode *sq; PhotonBCNode *hist; const float gecc2 = ray -> gecc * ray -> gecc; float r, totalWeight = 0; PhotonSQNode *squeueEnd; PhotonBCNode *histEnd; COLOR fluxLo, fluxMid, irradVar, irradAvg, p, d; if (!pmap -> biasCompHist) { /* Allocate bias compensation history */ numHi = pmap -> maxGather - pmap -> minGather; for (i = pmap -> minGather + 1; numHi > 1; numHi >>= 1, ++i); pmap -> biasCompHist = (PhotonBCNode*)malloc(i * sizeof(PhotonBCNode)); if (!pmap -> biasCompHist) error(USER, "can't allocate bias compensation history"); } numLo = min(pmap -> minGather, pmap -> squeueEnd - 1); numHi = min(pmap -> maxGather, pmap -> squeueEnd - 1); sq = squeueEnd = pmap -> squeue + pmap -> squeueEnd - 1; histEnd = pmap -> biasCompHist; setcolor(fluxLo, 0, 0, 0); /* Partition to get numLo closest photons starting from END of queue */ squeuePartition(squeueEnd, 1, numLo, numHi); /* Get irradiance estimates (ignoring constants) using 1..numLo photons and chuck in history. Queue is traversed BACKWARDS. */ for (i = 1; i <= numLo; i++, sq--) { /* Make sure furthest two photons are consecutive wrt distance */ squeuePartition(squeueEnd, i, i + 1, numHi); /* Compute phase function for inscattering from photon */ if (gecc2 <= FTINY) r = 1; else { r = DOT(ray -> rdir, sq -> photon -> norm) / 127; r = 1 + gecc2 - 2 * ray -> gecc * r; r = (1 - gecc2) / (r * sqrt(r)); } getPhotonFlux(sq -> photon, irrad); scalecolor(irrad, r); addcolor(fluxLo, irrad); /* Average radius between furthest two photons to improve accuracy */ r = 0.25 * (sq -> dist + (sq - 1) -> dist + 2 * sqrt(sq -> dist * (sq - 1) -> dist)); r *= sqrt(r); /* Add irradiance and weight to history. Weights should increase monotonically based on number of photons used for the estimate. */ histEnd -> irrad [0] = fluxLo [0] / r; histEnd -> irrad [1] = fluxLo [1] / r; histEnd -> irrad [2] = fluxLo [2] / r; totalWeight += histEnd++ -> weight = BIASCOMP_WGT((float)i); } /* Compute expected value (average) and variance of irradiance based on history for numLo photons. */ setcolor(irradAvg, 0, 0, 0); setcolor(irradVar, 0, 0, 0); for (hist = pmap -> biasCompHist; hist < histEnd; ++hist) for (i = 0; i <= 2; ++i) { irradAvg [i] += r = hist -> weight * hist -> irrad [i]; irradVar [i] += r * hist -> irrad [i]; } for (i = 0; i <= 2; ++i) { r = irradAvg [i] /= totalWeight; irradVar [i] = irradVar [i] / totalWeight - r * r; } /* Do binary search within interval [numLo, numHi]. numLo is towards the END of the queue. */ while (numHi - numLo > 1) { numMid = (numLo + numHi) >> 1; /* Split interval to get numMid closest photons starting from the END of the queue */ squeuePartition(squeueEnd, numLo, numMid, numHi); /* Make sure furthest two photons are consecutive wrt distance */ squeuePartition(squeueEnd, numMid, numMid + 1, numHi); copycolor(fluxMid, fluxLo); sq = squeueEnd - numLo; /* Get irradiance for numMid photons (ignoring constants) */ for (i = numLo; i < numMid; i++, sq--) { /* Compute phase function for inscattering from photon */ if (gecc2 <= FTINY) r = 1; else { r = DOT(ray -> rdir, sq -> photon -> norm) / 127; r = 1 + gecc2 - 2 * ray -> gecc * r; r = (1 - gecc2) / (r * sqrt(r)); } getPhotonFlux(sq -> photon, irrad); scalecolor(irrad, r); addcolor(fluxMid, irrad); } /* Average radius between furthest two photons to improve accuracy */ r = 0.25 * (sq -> dist + (sq + 1) -> dist + 2 * sqrt(sq -> dist * (sq + 1) -> dist)); r *= sqrt(r); /* Get deviation from current average, and probability that it's due to noise from gaussian distribution based on current variance. Since we are doing this for each colour channel we can also detect chromatic bias. */ for (i = 0; i <= 2; ++i) { d [i] = irradAvg [i] - (irrad [i] = fluxMid [i] / r); p [i] = exp(-0.5 * d [i] * d [i] / irradVar [i]); } if (pmapRandom(pmap -> randState) < colorAvg(p)) { /* Deviation is probably noise, so add mid irradiance to history */ copycolor(histEnd -> irrad, irrad); totalWeight += histEnd++ -> weight = BIASCOMP_WGT((float)numMid); setcolor(irradAvg, 0, 0, 0); setcolor(irradVar, 0, 0, 0); /* Update average and variance */ for (hist = pmap -> biasCompHist; hist < histEnd; ++hist) for (i = 0; i <= 2; i++) { r = hist -> irrad [i]; irradAvg [i] += hist -> weight * r; irradVar [i] += hist -> weight * r * r; } for (i = 0; i <= 2; i++) { r = irradAvg [i] /= totalWeight; irradVar [i] = irradVar [i] / totalWeight - r * r; } /* Need more photons -- recurse on [numMid, numHi] */ numLo = numMid; copycolor(fluxLo, fluxMid); } else /* Deviation is probably bias -- need fewer photons, so recurse on [numLo, numMid] */ numHi = numMid; } --histEnd; for (i = 0; i <= 2; i++) { /* Estimated relative error */ d [i] = histEnd -> irrad [i] / irradAvg [i] - 1; /* Divide by 4 / 3 * PI for search volume (r^3 already accounted for) and phase function normalization factor 1 / (4 * PI) */ irrad [i] = histEnd -> irrad [i] * 3 / (16 * PI * PI); } /* Update statistix */ r = colorAvg(d); if (r < pmap -> minError) pmap -> minError = r; if (r > pmap -> maxError) pmap -> maxError = r; pmap -> rmsError += r * r; if (numMid < pmap -> minGathered) pmap -> minGathered = numMid; if (numMid > pmap -> maxGathered) pmap -> maxGathered = numMid; pmap -> totalGathered += numMid; ++pmap -> numDensity; } #endif diff --git a/pmapbias.h b/pmapbias.h index 9939297..091543d 100644 --- a/pmapbias.h +++ b/pmapbias.h @@ -1,44 +1,50 @@ /* RCSid $Id: pmapbias.h,v 2.5 2016/05/17 17:39:47 rschregle Exp $ */ /* ================================================================== Bias compensation for photon density estimates For background see: R. Schregle, "Bias Compensation for Photon Maps", Computer Graphics Forum, v22:n4, pp. 729-742, Dec. 2003. Roland Schregle (roland.schregle@gmail.com) (c) Fraunhofer Institute for Solar Energy Systems ================================================================== $Id: pmapbias.h,v 2.5 2016/05/17 17:39:47 rschregle Exp $ */ #ifndef PMAPBIASCOMP_H #define PMAPBIASCOMP_H #include "pmapdata.h" + /* Bias compensation history node */ + typedef struct { + COLOR irrad; + float weight; + } PhotonBiasCompNode; + /* Bias compensation weighting function */ /* #define BIASCOMP_WGT(n) 1 */ /* #define BIASCOMP_WGT(n) (n) */ #define BIASCOMP_WGT(n) ((n) * (n)) /* #define BIASCOMP_WGT(n) ((n) * (n) * (n)) */ /* #define BIASCOMP_WGT(n) exp(0.003 * (n)) */ /* Dump photon bandwidth for bias compensated density estimates */ /* #define BIASCOMP_BWIDTH */ void biasComp (PhotonMap*, COLOR); /* Photon density estimate with bias compensation, returning irradiance. Expects photons in search queue after a kd-tree lookup. */ void volumeBiasComp (PhotonMap*, const RAY*, COLOR); /* Photon volume density estimate with bias compensation, returning irradiance. Expects photons in search queue after a kd-tree lookup. */ #endif diff --git a/pmapcontrib.c b/pmapcontrib.c index c817a55..b38110f 100644 --- a/pmapcontrib.c +++ b/pmapcontrib.c @@ -1,954 +1,978 @@ #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" extern void SDdisk2square(double sq[2], double diskx, double disky); static int logDim (unsigned size) /* Return log2(sqrt(size)) = l, where size = (2^l)(2^l). Is there a faster way (e.g. binary search)? */ { unsigned i, sz, dim; if (!size) return 0; for (i = 0, dim = 0, sz = size; sz >>= 2; dim++); return dim; } static int xy2bin (unsigned l, int x, int y) /* Serialise 2D coordinates in range (2^l) x (2^l) to 1D bin. Returns -1 if coordinates invalid */ { return x < 0 || y < 0 ? -1 : (x << l) + y; } static void bin2xy (unsigned l, int bin, int *x, int *y) /* Deserialise 1D bin to 2D coordinates in range (2^l) x (2^l). Returns -1 if bin invalid */ { /* Obviously this is faster than integer division/modulo */ *x = bin < 0 ? -1 : bin >> l; *y = bin < 0 ? -1 : bin & PMAP_CONTRIB_SCDIM(l) - 1; } static int ray2bin (const RAY *ray, unsigned nbins) /* Map ray dir (pointing away from origin) to its 1D Shirley-Chiu bin, where nbins = (2^l) x (2^l) for l > 1. Returns -1 if mapped bin is invalid (e.g. behind plane) */ { const unsigned l = logDim(nbins), scDim = PMAP_CONTRIB_SCDIM(l); static int scRHS, varInit = 0; static FVECT scNorm, scUp; unsigned scBin [2]; FVECT diskPlane; RREAL dz, diskx, disky, rad, diskd2, scCoord [2]; if (!varInit) { /* Lazy init shirley-Chiu mapping orientation from function variables */ scRHS = varvalue(PMAP_CONTRIB_SCRHS); scNorm [0] = varvalue(PMAP_CONTRIB_SCNORMX); scNorm [1] = varvalue(PMAP_CONTRIB_SCNORMY); scNorm [2] = varvalue(PMAP_CONTRIB_SCNORMZ); scUp [0] = varvalue(PMAP_CONTRIB_SCUPX); scUp [1] = varvalue(PMAP_CONTRIB_SCUPY); scUp [2] = varvalue(PMAP_CONTRIB_SCUPZ); varInit ^= 1; } /* Map incident direction to disk */ dz = DOT(ray -> rdir, scNorm); /* normal proj */ if (dz > 0) { fcross(diskPlane, scUp, scNorm); /* y-axis in plane, perp to up */ diskPlane [0] *= scRHS; diskx = DOT(ray -> rdir, diskPlane); disky = DOT(ray -> rdir, scUp) - dz * DOT(ray -> rdir, scUp); diskd2 = diskx * diskx + disky * disky; /* in-plane length^2 */ rad = dz>FTINY ? sqrt((1 - dz*dz) / diskd2) : 0; /* radial factor */ disky *= rad; /* diskx, disky now in range [-1, 1] */ /* Apply Shirley-Chiu mapping of (diskx, disky) to square */ SDdisk2square(scCoord, diskx, disky); /* Map Shirley-Chiu square coords to 1D bin */ scBin [0] = scCoord [0] * scDim; scBin [1] = scCoord [1] * scDim; return xy2bin(l, scBin [0], scBin [1]); } else return -1; } #ifndef PMAP_CONTRIB_TEST MODCONT *addContribModifier (LUTAB *contribTab, unsigned *numContribs, char *mod, char *binParm, int binCnt ) /* Add light source modifier mod to contribution lookup table contribTab, and update numContribs. Return initialised contribution data for this modifier. This code adapted from rcontrib.c:addmodifier(). */ { LUENT *lutEntry = lu_find(contribTab, mod); MODCONT *contrib; EPNODE *eBinVal; if (lutEntry -> data) { /* Reject duplicate modifiers */ sprintf(errmsg, "duplicate light source modifier %s", mod); error(USER, errmsg); } if (*numContribs >= MAXMODLIST) { sprintf(errmsg, "too many modifiers (max. %d)", MAXMODLIST); error(INTERNAL, errmsg); } lutEntry -> key = mod; if (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, 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( pmap -> contribTab, contribSrcRay ); if (pmap -> lastContribSrc.srcBin < 0) { /* Warn if invalid bin, but trace photon nonetheless. It will count as emitted to prevent bias, but will not be stored in newPhoton(), as it contributes zero flux */ sprintf(errmsg, "invalid bin for light source %s, " "contribution photons will be discarded", source [contribSrcRay -> rsrc].so -> oname ); 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; } static int coeffCompare (const void *c1, const void *c2) /* Comparison function to REVERSE sort thresholded coefficients */ { const ThreshWaveletCoeff *tcoeff1 = c1, *tcoeff2 = c2; /* TODO: USE MAX HERE INSTEAD OF AVERAGE? */ const WAVELET_COEFF v1 = colorAvg(tcoeff1 -> coeff), v2 = colorAvg(tcoeff2 -> coeff); if (v1 < v2) return 1; else if (v1 > v2) return -1; else return 0; } static int thresholdContribs (MODCONT *binnedContribs, PreComputedContrib *preCompContrib ) /* Threshold binned detail coefficients in binnedContribs -> cbin [1] by keeping the (preCompContrib -> nCompressedBins) largest and returning these in preContrib -> threshCoeffs along with their original bin order. NOTE: binnedContribs -> cbin [0] is the average wavelet coefficient and excluded from thresholding. Returns 0 on success. */ { unsigned b; ThreshWaveletCoeff *threshCoeffs; if (!preCompContrib) return -1; /* Skip 1st coeff; it's the average and therefore not thresholded */ for (b = 0, threshCoeffs = preCompContrib -> threshCoeffs; b < binnedContribs -> nbins - 1; b++ ) { /* Set up pointers to coeffs (sorted faster than 3 floats/doubles) and remember original bin position prior to sorting. NOTE: The original bin position *excludes* the 1st coeff, i.e. the average! */ threshCoeffs [b].coeff = (WAVELET_COEFF*)( &(binnedContribs->cbin [b+1]) ); threshCoeffs [b].bin = b; } /* REVERSE sort coeffs by magnitude; the non-thresholded coeffs will then start at threshCoeffs [0] */ qsort(threshCoeffs, binnedContribs -> nbins - 1, sizeof(ThreshWaveletCoeff), coeffCompare ); return 0; } static int encodeContribs (MODCONT *binnedContribs, PreComputedContrib *preCompContrib, float compressRatio ) /* Apply wavelet transform to binnedContribs -> cbin 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 binnedContribs -> cbin [0]. Returns 0 on success. */ { unsigned b, i; ThreshWaveletCoeff *threshCoeffs; mRGBERange *mrgbeRange; mRGBE *mrgbeCoeff; if (!binnedContribs || !preCompContrib) return -1; /* Do 2D wavelet transform on preCompContrib -> waveletMatrix (which really maps to binnedContribs -> cbin) */ if (waveletXform2(preCompContrib -> waveletMatrix, preCompContrib -> tWaveletMatrix, preCompContrib -> l ) < 0 ) error(INTERNAL, "failed wavelet transform of contributions"); /* Compress wavelet detail coeffs by thresholding */ if (thresholdContribs(binnedContribs, preCompContrib) < 0) error(INTERNAL, "failed wavelet compression of contributions"); threshCoeffs = preCompContrib -> threshCoeffs; /* 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 (b = 0; b < preCompContrib -> nCompressedBins; b++) { for (i = 0; i < 3; i++) { if (threshCoeffs [b].coeff [i] < mrgbeRange -> min [i]) mrgbeRange -> min [i] = threshCoeffs [b].coeff [i]; if (threshCoeffs [b].coeff [i] > mrgbeRange -> max [i]) mrgbeRange -> max [i] = threshCoeffs [b].coeff [i]; } } /* Init mRGBE coefficient normalisation from range */ mRGBEinit(mrgbeRange, mrgbeRange -> min, mrgbeRange -> max); mrgbeCoeff = preCompContrib -> mrgbeCoeffs; /* Encode wavelet detail coefficients to mRGBE */ for (b = 0; b < preCompContrib -> nCompressedBins; b++) mrgbeCoeff [b] = mRGBEencode(threshCoeffs [b].coeff, mrgbeRange, threshCoeffs [b].bin ); 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( pmap -> contribTab, srcMod -> oname ) -> data; Photon *photon; COLOR flux; PhotonSearchQueueNode *sqn; float r2, norm; unsigned i; /* Zero bins for this source modifier */ memset(contrib -> cbin, 0, sizeof(DCOLOR) * contrib -> nbins); setcolor(irrad, 0, 0, 0); if (!contrib || !pmap -> maxGather) /* Modifier not in LUT or zero bandwidth */ return NULL; /* Lookup photons */ pmap -> squeue.tail = 0; /* Pass light source index to filter in findPhotons() via pmap -> lastContribSrc */ pmap -> lastContribSrc.srcIdx = ray -> rsrc; findPhotons(pmap, ray); /* Need at least 2 photons */ if (pmap -> squeue.tail < 2) { #ifdef PMAP_NONEFOUND sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)", ray -> ro ? ray -> ro -> oname : "", ray -> rop [0], ray -> rop [1], ray -> rop [2] ); error(WARNING, errmsg); #endif return NULL; } /* Avg radius^2 between furthest two photons to improve accuracy */ sqn = pmap -> squeue.node + 1; r2 = max(sqn -> dist2, (sqn + 1) -> dist2); r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2)); /* XXX: OMIT extra normalisation factor 1/PI for ambient calc? */ #ifdef PMAP_EPANECHNIKOV /* Normalise accumulated flux by Epanechnikov kernel integral in 2D (see Eq. 4.4, p.76 in Silverman, "Density Estimation for Statistics and Data Analysis", 1st Ed., 1986, and Wann Jensen, "Realistic Image Synthesis using Photon Mapping"), include RADIANCE-specific lambertian factor PI */ norm = 2 / (PI * PI * r2); #else /* Normalise accumulated flux by search area PI * r^2, including RADIANCE-specific lambertian factor PI */ norm = 1 / (PI * PI * r2); #endif /* Skip the extra photon */ for (i = 1 ; i < pmap -> squeue.tail; i++, sqn++) { /* Get photon's contribution to density estimate */ photon = getNearestPhoton(&pmap -> squeue, sqn -> idx); getPhotonFlux(photon, flux); scalecolor(flux, norm); #ifdef PMAP_EPANECHNIKOV /* Apply Epanechnikov kernel to photon flux based on photon distance */ scalecolor(flux, 1 - sqn -> dist2 / r2); #endif addcolor(irrad, flux); addcolor(contrib -> cbin [photonSrcBin(pmap, photon)], flux); } return contrib; } static void freePreCompContribNode (void *p) /* Free per-modifier precomputed contributions LUT entry */ { PhotonMap *preCompContribPmap = (PhotonMap*)p; PreComputedContrib *preCompContrib; if (preCompContribPmap) { preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib ); if (preCompContrib) { /* Free primary and transposed wavelet matrices */ free(preCompContrib -> waveletMatrix); freeWaveletMatrix(preCompContrib -> tWaveletMatrix, preCompContrib -> l ); /* Free thresholded coefficients */ free(preCompContrib -> threshCoeffs); /* Free mRGBE encoded coefficients */ free(preCompContrib -> mrgbeCoeffs); free(preCompContrib -> waveletFname); } /* Clean up precomputed contrib photon map */ deletePhotons(preCompContribPmap); free(preCompContribPmap); } } 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. */ - /* !!! NOTE: PRECOMPUTATION CURRENTLY WITHOUT OOC CACHE !!! */ { unsigned long i, numPreComp; unsigned b, j; long nCompressedBins; PhotonIdx pIdx; Photon photon; RAY ray; MODCONT *binnedContribs; COLR coeffNorm; LUENT *preCompContribNode; PhotonMap *preCompContribPmap, nuPmap; PreComputedContrib *preCompContrib; LUTAB lutInit = LU_SINIT(NULL, freePreCompContribNode ); 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; /* 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 (i = 0; i < numPreComp; i++) { - do { - /* Get random photon from stratified distribution in source heap - * to avoid duplicates and clustering */ - pIdx = firstPhoton(pmap) + (unsigned long)( - (i + pmapRandom(pmap -> randState)) / finalGather - ); - getPhoton(pmap, pIdx, &photon); - - /* Init dummy photon ray with intersection and normal at photon - position. Set emitting light source index from origin. */ - VCOPY(ray.rop, photon.pos); - ray.rsrc = photonSrcIdx(pmap, &photon); - for (j = 0; j < 3; j++) - ray.ron [j] = photon.norm [j] / 127.0; - - /* Get contributions at photon position; retry with another - photon if no contribs found */ - binnedContribs = getPhotonContrib(pmap, &ray, ray.rcol); - } while (!binnedContribs); - - if (!(preCompContribNode = lu_find(nuPmap.preCompContribTab, - binnedContribs -> modname) + /* Get random photon from stratified distribution in source heap + * to avoid duplicates and clustering */ + pIdx = firstPhoton(pmap) + (unsigned long)( + (i + pmapRandom(pmap -> randState)) / finalGather + ); + getPhoton(pmap, pIdx, &photon); + + /* Init dummy photon ray with intersection and normal at photon + position. Set emitting light source index from origin. */ + VCOPY(ray.rop, photon.pos); + ray.rsrc = photonSrcIdx(pmap, &photon); + for (j = 0; j < 3; j++) + ray.ron [j] = photon.norm [j] / 127.0; + + /* Get contributions at photon position; retry with another + photon if no contribs found */ + binnedContribs = getPhotonContrib(pmap, &ray, ray.rcol); + + if (binnedContribs) { + /* 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 contrib heap */ - initPhotonMap(preCompContribPmap, nuPmap.type); - initPhotonHeap(preCompContribPmap); - initContribHeap(preCompContribPmap); - preCompContrib = (PreComputedContrib*)( - preCompContribPmap -> preCompContrib = - malloc(sizeof(PreComputedContrib)) - ); - } - if (!preCompContribPmap || !preCompContrib) - error(SYSTEM, "out of memory allocating new photon map " - "in preComputeContrib()" + error(SYSTEM, + "out of memory allocating LUT entry in preComputeContrib()" ); - /* Set output filename from parent photon map */ - preCompContribPmap -> fileName = nuPmap.fileName; + if (!preCompContribNode -> key) { + /* New LUT node for precomputed contribs for this modifier */ + preCompContribNode -> key = (char*)binnedContribs -> modname; + preCompContribNode -> data = (char*)( + preCompContribPmap = malloc(sizeof(PhotonMap)) + ); - /* Get Shirley-Chiu square & wavelet matrix resolution, bins - per dimension, and fixed number of compressed coeffs/bins - (note we subtract one as the average coefficient is not - thresholded) */ - preCompContrib -> l = logDim(binnedContribs -> nbins); - preCompContrib -> scDim = PMAP_CONTRIB_SCDIM(preCompContrib -> l); + if (preCompContribPmap) { + /* Init new child photon map and its contrib heap */ + initPhotonMap(preCompContribPmap, PMAP_TYPE_CONTRIB_CHILD); + initPhotonHeap(preCompContribPmap); + initContribHeap(preCompContribPmap); + preCompContrib = (PreComputedContrib*)( + preCompContribPmap -> preCompContrib = + malloc(sizeof(PreComputedContrib)) + ); + } + if (!preCompContribPmap || !preCompContrib) + error(SYSTEM, + "out of memory allocating new photon map " + "in preComputeContrib()" + ); - if (binnedContribs -> nbins > 1) { - /* Binning enabled; set up wavelet xform & compression */ - nCompressedBins = (binnedContribs -> nbins - 1) * - (1 - compressRatio); - if (nCompressedBins > 0) - preCompContrib -> nCompressedBins = nCompressedBins; - else error(USER, - "invalid compression ratio in preComputeContrib()" - ); + /* Set output filename from parent photon map */ + preCompContribPmap -> fileName = nuPmap.fileName; - /* Lazily "allocate" primary wavelet matrix pointing to array - binnedContribs -> cbin; this imposes the 2D matrix structure - expected by the wavelet xform routines onto the 1D array. */ - if (!(preCompContrib -> waveletMatrix = calloc( - preCompContrib -> scDim, sizeof(WaveletCoeff3*)) - ) - ) - error(SYSTEM, "out of memory allocating primary wavelet " - "matrix in preComputeContrib()" + /* Get Shirley-Chiu square & wavelet matrix resolution, bins + per dimension, and fixed number of compressed coeffs/bins + (note we subtract one as the average coefficient is not + thresholded) */ + preCompContrib -> l = logDim(binnedContribs -> nbins); + preCompContrib -> scDim = PMAP_CONTRIB_SCDIM(preCompContrib -> l); + + if (binnedContribs -> nbins > 1) { + /* Binning enabled; set up wavelet xform & compression */ + nCompressedBins = (binnedContribs -> nbins - 1) * + (1 - compressRatio); + if (nCompressedBins > 0) + preCompContrib -> nCompressedBins = nCompressedBins; + else error(USER, + "invalid compression ratio in preComputeContrib()" ); - for (b = 0; b < preCompContrib -> scDim; b++) - /* Point to each row in existing 1D contrib array */ - preCompContrib -> waveletMatrix [b] = - &binnedContribs -> cbin [b * preCompContrib -> scDim]; - - /* Lazily allocate transposed wavelet matrix */ - preCompContrib -> tWaveletMatrix = - allocWaveletMatrix(preCompContrib -> l); - if (!preCompContrib -> tWaveletMatrix) - error(SYSTEM, "out of memory allocating transposed wavelet " - "matrix in preComputeContrib()" + /* Lazily "allocate" primary wavelet matrix pointing to array + binnedContribs -> cbin; this imposes the 2D matrix structure + expected by the wavelet xform routines onto the 1D array. */ + if (!(preCompContrib -> waveletMatrix = calloc( + preCompContrib -> scDim, sizeof(WaveletCoeff3*)) + ) + ) + error(SYSTEM, "out of memory allocating primary wavelet " + "matrix in preComputeContrib()" + ); + + for (b = 0; b < preCompContrib -> scDim; b++) + /* Point to each row in existing 1D contrib array */ + preCompContrib -> waveletMatrix [b] = + &binnedContribs -> cbin [b * preCompContrib -> scDim]; + + /* Lazily allocate transposed wavelet matrix */ + preCompContrib -> tWaveletMatrix = + allocWaveletMatrix(preCompContrib -> l); + if (!preCompContrib -> tWaveletMatrix) + error(SYSTEM, + "out of memory allocating transposed wavelet " + "matrix in preComputeContrib()" + ); + + /* Lazily allocate thresholded detail coefficients (minus + the average coeff) */ + preCompContrib -> threshCoeffs = calloc( + binnedContribs -> nbins - 1, sizeof(ThreshWaveletCoeff) ); - - /* Lazily allocate thresholded detail coefficients (minus - the average coeff) */ - preCompContrib -> threshCoeffs = calloc( - binnedContribs -> nbins - 1, sizeof(ThreshWaveletCoeff) - ); - /* Lazily allocate mRGBE-encoded, compressed wavelet coeffs */ - preCompContrib -> mrgbeCoeffs = calloc( - preCompContrib -> nCompressedBins, sizeof(mRGBE) - ); - if (!preCompContrib -> threshCoeffs || - !preCompContrib -> mrgbeCoeffs - ) - error(SYSTEM, "out of memory allocating compressed " - "coefficients in preComputeContrib()" + /* Lazily alloc mRGBE-encoded, compressed wavelet coeffs */ + preCompContrib -> mrgbeCoeffs = calloc( + preCompContrib -> nCompressedBins, sizeof(mRGBE) ); + if (!preCompContrib -> threshCoeffs || + !preCompContrib -> mrgbeCoeffs + ) + error(SYSTEM, "out of memory allocating compressed " + "coefficients in preComputeContrib()" + ); + + /* Set size of compressed contributions, in bytes, + including prepended normalisation */ + preCompContrib -> contribSize = sizeof(coeffNorm) + + sizeof(mRGBE) * preCompContrib -> nCompressedBins; + } + else { + /* No binning; skip wavelet xform & compression */ + preCompContrib -> nCompressedBins = + preCompContrib -> contribSize = 0; + preCompContrib -> waveletMatrix = + preCompContrib -> tWaveletMatrix = NULL; + preCompContrib -> threshCoeffs = NULL; + preCompContrib -> mrgbeCoeffs = NULL; + } } else { + /* LUT node already exists for this modifier; use its data */ preCompContribPmap = (PhotonMap*)preCompContribNode -> data; preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib ); } + + if (binnedContribs -> nbins > 1) { + /* Binning enabled; compress & encode binned contribs */ + if (encodeContribs(binnedContribs, preCompContrib, + compressRatio + ) + ) + error(INTERNAL, + "failed contribution compression/encoding " + "in preComputeContrib()" + ); - /* Compress & encode binned contribs */ - if (encodeContribs(binnedContribs, preCompContrib, - compressRatio - ) < 0 - ) - error(INTERNAL, "failed contribution compression/encoding " - "in preComputeContrib()" + /* Dump encoded bins to contrib heap file, prepended by + normalisation in 32-bit RGBE */ + setcolr(coeffNorm, + preCompContrib -> mrgbeRange.norm [0], + preCompContrib -> mrgbeRange.norm [1], + preCompContrib -> mrgbeRange.norm [2] ); - - /* Dump encoded bins to contrib heap file, prepended by - normalisation in 32-bit RGBE */ - setcolr(coeffNorm, - preCompContrib -> mrgbeRange.norm [0], - preCompContrib -> mrgbeRange.norm [1], - preCompContrib -> mrgbeRange.norm [2] - ); - putbinary(coeffNorm, 1, sizeof(coeffNorm), - preCompContribPmap -> contribHeap - ); - - if (putbinary(preCompContrib -> mrgbeCoeffs, sizeof(mRGBE), - preCompContrib -> nCompressedBins, + #if 0 + /* Replace normalisation with photon's Morton code + (truncated to 32 bits) to check contribPhotonSort() */ + FVECT org = {0,0,0}; + *(uint32*)coeffNorm = Key2Morton3D(ray.rop, org, + MORTON3D_MAX / FHUGE + ) >> 32; + #endif + putbinary(coeffNorm, 1, sizeof(coeffNorm), preCompContribPmap -> contribHeap - ) != preCompContrib -> nCompressedBins - ) - error(SYSTEM, "failed writing to coefficients to " - "contribution heap in preComputeContrib()" ); + + if (putbinary(preCompContrib -> mrgbeCoeffs, sizeof(mRGBE), + preCompContrib -> nCompressedBins, + preCompContribPmap -> contribHeap + ) != preCompContrib -> nCompressedBins + ) + error(SYSTEM, "failed writing to coefficients to " + "contribution heap in preComputeContrib()" + ); + } + + /* Set photon flux to coarsest average wavelet coefficient + NOTE: binnedContrib -> cbin [0] == + preCompContrib -> waveletMatrix [0][0] == + preCompContrib -> tWaveletMatrix [0][0]. + IF BINNING IS DISABLED, this is the total contribution + from all directions */ + copycolor(ray.rcol, binnedContribs -> cbin [0]); + + /* HACK: signal newPhoton() to set precomputed photon's + contribution source from ray -> rsrc */ + /* XXX: DO WE STILL NEED THIS SHIT 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++; } - else { - /* No binning; skip wavelet xform & compression */ - preCompContrib -> nCompressedBins = 0; - preCompContrib -> waveletMatrix = - preCompContrib -> tWaveletMatrix = NULL; - preCompContrib -> threshCoeffs = NULL; - preCompContrib -> mrgbeCoeffs = NULL; - } - - /* Set photon flux to coarsest average wavelet coefficient - NOTE: binnedContrib -> cbin [0] == - preCompContrib -> waveletMatrix [0][0] == - preCompContrib -> tWaveletMatrix [0][0]. - If binning is disabled, this is the total contribution from all - directions */ - copycolor(ray.rcol, binnedContribs -> cbin [0]); - - /* HACK: signal newPhoton() to set precomputed photon's - contribution source from ray -> rsrc */ - 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 replace with new parent, which now acts - * as a for the per-modifier child pmaps */ + /* 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 + #endif /* 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, l, nbins, nsamp, numTheta = 0, numPhi = 0; double t, p; RAY ray; int bin; if (argc < 3) { fprintf(stderr, "%s [=; ..]\n", argv [0]); fprintf(stderr, "Default variable defs: %s\n", PMAP_CONTRIB_SCDEFAULTS ); fputs("\nMissing resolution l>1, number of samples nsamp\n", stderr ); return -1; } if (!(l = atoi(argv [1]))) { fputs("Invalid resolution\n", stderr); return -1; } else nbins = PMAP_CONTRIB_SCBINS(l); if (!(nsamp = atoi(argv [2]))) { fputs("Invalid num samples\n", stderr); return -1; } else { numTheta = (int)(sqrt(nsamp) / 2); numPhi = 4* numTheta; } /* (Doan' need to) Init cal func routines for binning */ #if 0 initfunc(); setcontext(RCCONTEXT); #endif /* Compile default orientation variables for contribution binning */ scompile(PMAP_CONTRIB_SCDEFAULTS, NULL, 0); /* Compile custom orientation variabls from command line */ for (i = 3; i < argc; i++) scompile(argv [i], NULL, 0); for (i = 0; i < nsamp; i++) { #if 0 /* Random */ t = 0.5 * PI * drand48(); p = 2 * PI * drand48(); #else /* Stratified */ t = 0.5 * PI * ((i % numTheta) + drand48()) / numTheta; p = 2.0 * PI * ((i / numTheta) + drand48()) / numPhi; #endif ray.rdir [0] = sin(t) * cos(p); ray.rdir [1] = sin(t) * sin(p); ray.rdir [2] = cos(t); bin = ray2bin(&ray, nbins); printf("%.3f\t%.3f\t%.3f\t-->\t%d\n", ray.rdir [0], ray.rdir [1], ray.rdir [2], bin ); } return 0; } #endif diff --git a/pmapcontrib.h b/pmapcontrib.h index 5db6188..cb7173e 100644 --- a/pmapcontrib.h +++ b/pmapcontrib.h @@ -1,156 +1,157 @@ /* RCSid $Id: pmapcontrib.h,v 2.5 2016/05/17 17:39:47 rschregle Exp $ */ /* ========================================================================= Photon map routines for precomputed light source contributions; these routines interface to mkpmap and rcontrib. Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") ========================================================================= $Id: pmapcontrib.h,v 2.5 2016/05/17 17:39:47 rschregle Exp $ */ #ifndef _PMAPCONTRIB_H #define _PMAPCONTRIB_H #include "pmapdata.h" #include "wavelet2.h" #include "mrgbe.h" #ifndef MAXPROCESS #include "rcontrib.h" #endif #ifndef MAXMODLIST /* Max number of contributing sources */ #define MAXMODLIST 1024 #endif /* Filename templates for 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 binning resolution l*/ #define PMAP_CONTRIB_MAXBINRES (mRGBE_DATABITS >> 1) /* Shirley-Chiu square dimensions and number of bins for resolution l */ #define PMAP_CONTRIB_SCDIM(l) (1 << (l)) #define PMAP_CONTRIB_SCBINS(l) (1 << ((l) << 1)) typedef struct { WAVELET_COEFF *coeff; unsigned bin; } ThreshWaveletCoeff; typedef struct { char *waveletFname; FILE *waveletFile; WaveletMatrix waveletMatrix, tWaveletMatrix; ThreshWaveletCoeff *threshCoeffs; mRGBERange mrgbeRange; mRGBE *mrgbeCoeffs; unsigned l, scDim, nCompressedBins; + unsigned long contribSize; } PreComputedContrib; 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 initPmapContrib (LUTAB *contribTab); /* Set up photon map contributions (interface to rcmain.c) */ PhotonContribSourceIdx newPhotonContribSource (PhotonMap *pmap, const RAY *contribSrcRay, FILE *contribSrcHeap ); /* Add contribution source ray for emitted contribution photon and save * light source index and binned direction. The current contribution source * is stored in pmap -> lastContribSrc. If the previous contribution source * spawned photons (i.e. has srcIdx >= 0), it's appended to contribSrcHeap. * If contribSrcRay == NULL, the current contribution source is still * flushed, but no new source is set. Returns updated contribution source * counter pmap -> numContribSrc. */ PhotonContribSourceIdx buildContribSources (PhotonMap *pmap, FILE **contribSrcHeap, char **contribSrcHeapFname, PhotonContribSourceIdx *contribSrcOfs, unsigned numHeaps ); /* Consolidate per-subprocess contribution sources heaps into array * pmap -> contribSrc. Returns offset for contribution source index * linearisation in pmap -> numContribSrc. The heap files in * contribSrcHeap are closed on return. */ void preComputeContrib (PhotonMap *pmap); /* Precompute contributions for a random subset of (finalGather * pmap -> numPhotons) photons, and init the per-modifier precomputed contribution photon maps in LUT pmap -> preCompContribTab, discarding the original photons. */ void distribPhotonContrib (PhotonMap *pmap, LUTAB *contribTab, unsigned numContribs, unsigned numProc ); /* Emit photons with binned light source contributions, precompute their contributions and build photon map */ int buildPreCompContribPmap (const LUENT *preCompContribNode, void *p); /* Build per-modifier precomputed photon map. Returns 0 on success. */ void saveContribPhotonMap (const PhotonMap *pmap, int argc, char **argv); /* Save contribution photon map and with its constituent per-modifier photon maps in pmap -> preCompContribTab */ void getPreCompPhotonContrib (PhotonMap *pmap, RAY *ray, COLOR totalContrib ); /* Fetch and decode precomputed light source contributions from single closest precomputed contribution photon at ray -> rop, and accumulate them in pmap -> contribTab. Also returns total precomputed contribs. */ #endif diff --git a/pmapdata.c b/pmapdata.c index 91c7c6a..6e52f4a 100644 --- a/pmapdata.c +++ b/pmapdata.c @@ -1,852 +1,858 @@ #ifndef lint static const char RCSid[] = "$Id: pmapdata.c,v 2.23 2021/04/14 11:26:25 rschregle Exp $"; #endif /* ========================================================================= Photon map types and high level interface to nearest neighbour lookups in underlying point cloud data structure. The default data structure is an in-core kd-tree (see pmapkdt.{h,c}). This can be overriden with the PMAP_OOC compiletime switch, which enables an out-of-core octree (see oococt.{h,c}). Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Fraunhofer Institute for Solar Energy Systems, supported by the German Research Foundation (DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme" (FARESYS)) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ========================================================================= $Id: pmapdata.c,v 2.23 2021/04/14 11:26:25 rschregle Exp $ */ #include "pmapdata.h" #include "pmapdens.h" #include "pmaprand.h" #include "pmapmat.h" #include "pmaproi.h" #include "otypes.h" #include "random.h" PhotonMap *photonMaps [NUM_PMAP_TYPES] = { - NULL, NULL, NULL, NULL, NULL, NULL, NULL + NULL, NULL, NULL, NULL, NULL, NULL, NULL, #ifdef PMAP_PHOTONFLOW - , NULL, NULL + NULL, NULL, #endif + NULL }; /* Include routines to handle underlying point cloud data structure */ #ifdef PMAP_OOC #include "pmapooc.c" #else #include "pmapkdt.c" #include "pmaptkdt.c" #endif /* Ambient include/exclude set (from ambient.c) */ #ifndef MAXASET #define MAXASET 4095 #endif extern OBJECT ambset [MAXASET+1]; /* Callback to print photon attributes acc. to user defined format */ int (*printPhoton)(RAY *r, Photon *p, PhotonMap *pm); /* PHOTON MAP BUILD ROUTINES ------------------------------------------ */ void initPhotonMap (PhotonMap *pmap, PhotonMapType t) /* Init photon map 'n' stuff... */ { if (!pmap) return; pmap -> numPhotons = 0; pmap -> biasCompHist = NULL; pmap -> maxPos [0] = pmap -> maxPos [1] = pmap -> maxPos [2] = -FHUGE; pmap -> minPos [0] = pmap -> minPos [1] = pmap -> minPos [2] = FHUGE; pmap -> minGathered = pmap -> maxGathered = pmap -> totalGathered = 0; pmap -> gatherTolerance = gatherTolerance; pmap -> minError = pmap -> maxError = pmap -> rmsError = 0; pmap -> numDensity = 0; pmap -> distribRatio = 1; pmap -> type = t; pmap -> squeue.node = NULL; pmap -> squeue.len = 0; #ifdef PMAP_PATHFILT /* Init path lookup table and its key buffer */ pmap -> pathLUT = NULL; pmap -> pathLUTKeys = NULL; pmap -> numPathLUTKeys = 0; #endif /* Init stats */ setcolor(pmap -> photonFlux, 0, 0, 0); pmap -> CoG [0] = pmap -> CoG [1] = pmap -> CoG [2] = 0; pmap -> CoGdist = 0; pmap -> maxDist2 = 0; /* Init transient photon mapping stuff */ pmap -> velocity = photonVelocity; pmap -> minPathLen = pmap -> maxPathLen = pmap -> avgPathLen = 0; /* Init local RNG state */ pmap -> randState [0] = 10243; pmap -> randState [1] = 39829; pmap -> randState [2] = 9433; pmapSeed(randSeed, pmap -> randState); /* Set up type-specific photon lookup callback */ pmap -> lookup = pmapLookup [t]; /* Init precomputed contrib photon stuff */ pmap -> preCompContribTab = NULL; pmap -> contribHeap = NULL; pmap -> preCompContrib = NULL; /* Mark photon contribution source as unused */ pmap -> lastContribSrc.srcIdx = pmap -> lastContribSrc.srcBin = -1; pmap -> numContribSrc = 0; pmap -> contribSrc = NULL; /* Init storage */ pmap -> heap = NULL; pmap -> heapBuf = NULL; pmap -> heapBufLen = 0; #ifdef PMAP_OOC OOC_Null(&pmap -> store); #else kdT_Null(&pmap -> store); #endif } void initPhotonHeap (PhotonMap *pmap) { int fdFlags; if (!pmap) error(INTERNAL, "undefined photon map in initPhotonHeap()"); if (!pmap -> heap) { /* Open heap file */ mktemp(strcpy(pmap -> heapFname, PMAP_TMPFNAME)); if (!(pmap -> heap = fopen(pmap -> heapFname, "w+b"))) error(SYSTEM, "failed opening heap file in initPhotonHeap()"); #ifdef F_SETFL /* XXX is there an alternative needed for Windows? */ fdFlags = fcntl(fileno(pmap -> heap), F_GETFL); fcntl(fileno(pmap -> heap), F_SETFL, fdFlags | O_APPEND); #endif/* ftruncate(fileno(pmap -> heap), 0); */ } } void flushPhotonHeap (PhotonMap *pmap) { int fd; const unsigned long len = pmap -> heapBufLen * sizeof(Photon); if (!pmap) error(INTERNAL, "undefined photon map in flushPhotonHeap()"); if (!pmap -> heap || !pmap -> heapBuf) { /* Silently ignore undefined heap error(INTERNAL, "undefined heap in flushPhotonHeap()"); */ return; } /* Atomically seek and write block to heap */ /* !!! Unbuffered I/O via pwrite() avoids potential race conditions * !!! and buffer corruption which can occur with lseek()/fseek() * !!! followed by write()/fwrite(). */ fd = fileno(pmap -> heap); #ifdef DEBUG_PMAP sprintf(errmsg, "Proc %d: flushing %ld photons from pos %ld\n", getpid(), pmap -> heapBufLen, lseek(fd, 0, SEEK_END) / sizeof(Photon) ); eputs(errmsg); #endif /*if (pwrite(fd, pmap -> heapBuf, len, lseek(fd, 0, SEEK_END)) != len) */ if (write(fd, pmap -> heapBuf, len) != len) { error(SYSTEM, "failed append to heap file in flushPhotonHeap()"); /* Clean up temp file */ fclose(pmap -> heap); unlink(pmap -> heapFname); } #if NIX if (fsync(fd)) error(SYSTEM, "failed fsync in flushPhotonHeap()"); #endif pmap -> heapBufLen = 0; } #ifdef DEBUG_PMAP static int checkPhotonHeap (FILE *file) /* Check heap for nonsensical or duplicate photons */ { Photon p, lastp; int i, dup; rewind(file); memset(&lastp, 0, sizeof(lastp)); while (fread(&p, sizeof(p), 1, file)) { dup = 1; for (i = 0; i <= 2; i++) { if (p.pos [i] < thescene.cuorg [i] || p.pos [i] > thescene.cuorg [i] + thescene.cusize ) { sprintf(errmsg, "corrupt photon in heap at [%f, %f, %f]\n", p.pos [0], p.pos [1], p.pos [2] ); error(WARNING, errmsg); } dup &= p.pos [i] == lastp.pos [i]; } if (dup) { sprintf(errmsg, "consecutive duplicate photon in heap " "at [%f, %f, %f]\n", p.pos [0], p.pos [1], p.pos [2] ); error(WARNING, errmsg); } } return 0; } #endif int newPhoton (PhotonMap* pmap, const RAY* ray) { unsigned i; Photon photon; COLOR photonFlux; /* Account for distribution ratio */ if (!pmap || pmapRandom(pmap -> randState) > pmap -> distribRatio) return -1; /* Don't store on sources */ if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype)) return -1; /* Ignore photon if modifier in/outside exclude/include set */ if (ambincl != -1 && ray -> ro && ambincl != inset(ambset, ray->ro->omod) ) return -1; /* Store photon if within a region of interest (for ze Ecksperts!) */ if (!photonInROI(ray)) return -1; /* Adjust flux according to distribution ratio and ray weight */ copycolor(photonFlux, ray -> rcol); #if 0 /* Factored out ray -> rweight as deprecated (?) for pmap, and infact erroneously attenuates volume photon flux based on extinction, which is already factored in by photonParticipate() */ scalecolor( photonFlux, ray -> rweight / (pmap -> distribRatio ? pmap -> distribRatio : 1) ); #else scalecolor(photonFlux, 1.0 / (pmap -> distribRatio ? pmap -> distribRatio : 1) ); #endif setPhotonFlux(&photon, photonFlux); /* Set photon position and flags */ VCOPY(photon.pos, ray -> rop); photon.flags = 0; photon.caustic = PMAP_CAUSTICRAY(ray); /* Set photon's subprocess index (to consolidate contribution sources after photon distribution completes) */ photon.proc = PMAP_GETRAYPROC(ray); switch (pmap -> type) { case PMAP_TYPE_CONTRIB: if (pmap -> lastContribSrc.srcIdx < -1) { /* HACK: Contrib photon being precomputed; set contribution source from index passed in ray */ photon.aux.contribSrc = ray -> rsrc; } else { /* HACK: Contrib photon before precomputation (i.e. in forward pass); set contribution source from last index in contrib source array. Note the contribution source bin has already been set by newPhotonContribSrc(). */ photon.aux.contribSrc = pmap -> numContribSrc; /* Photon will be stored; set contribution source index, * thereby marking it as having spawned photon(s) */ pmap -> lastContribSrc.srcIdx = ray -> rsrc; } break; case PMAP_TYPE_TRANSIENT: #ifdef PMAP_PHOTONFLOW case PMAP_TYPE_TRANSLIGHTFLOW: #endif /* Set cumulative path length for transient photon, corresponding to its time of flight. The cumulative path length is obtained relative to the maximum path length photonMaxPathDist. The last intersection distance is subtracted as this isn't factored in until photonRay() generates a new photon ray. */ photon.aux.pathLen = photonMaxPathDist - ray -> rmax + ray -> rot; break; default: /* Set photon's path ID from ray serial num (supplemented by photon.proc below) */ photon.aux.pathID = ray -> rno; } /* Set normal */ for (i = 0; i <= 2; i++) photon.norm [i] = 127.0 * (isVolumePmap(pmap) #ifdef PMAP_PHOTONFLOW || isLightFlowPmap(pmap) #endif ? ray -> rdir [i] : ray -> ron [i] ); if (!pmap -> heapBuf) { /* Lazily allocate heap buffa */ #if NIX /* Randomise buffa size to temporally decorellate flushes in * multiprocessing mode */ srandom(randSeed + getpid()); pmap -> heapBufSize = PMAP_HEAPBUFSIZE * (0.5 + frandom()); #else /* Randomisation disabled for single processes on WIN; also useful * for reproducability during debugging */ pmap -> heapBufSize = PMAP_HEAPBUFSIZE; #endif if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)))) error(SYSTEM, "failed heap buffer allocation in newPhoton"); pmap -> heapBufLen = 0; } /* Photon initialised; now append to heap buffa */ memcpy(pmap -> heapBuf + pmap -> heapBufLen, &photon, sizeof(Photon)); if (++pmap -> heapBufLen >= pmap -> heapBufSize) /* Heap buffa full, flush to heap file */ flushPhotonHeap(pmap); pmap -> numPhotons++; /* Print photon attributes */ if (printPhoton) /* Non-const kludge */ printPhoton((RAY*)ray, &photon, pmap); return 0; } void buildPhotonMap (PhotonMap *pmap, double *photonFlux, PhotonContribSourceIdx *contribSrcOfs, unsigned nproc ) { unsigned long n, nCheck = 0; unsigned i; Photon *p; COLOR flux; char nuHeapFname [sizeof(PMAP_TMPFNAME)]; FILE *nuHeap; /* Need double here to reduce summation errors */ double avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0, avgPathLen = 0; FVECT d; if (!pmap) error(INTERNAL, "undefined photon map in buildPhotonMap()"); /* Get number of photons from heapfile size */ if (fseek(pmap -> heap, 0, SEEK_END) < 0) error(SYSTEM, "failed seek to end of photon heap in buildPhotonMap()"); pmap -> numPhotons = ftell(pmap -> heap) / sizeof(Photon); if (!pmap -> numPhotons) error(INTERNAL, "empty photon map in buildPhotonMap()"); if (!pmap -> heap) error(INTERNAL, "no heap in buildPhotonMap()"); #ifdef DEBUG_PMAP eputs("Checking photon heap consistency...\n"); checkPhotonHeap(pmap -> heap); sprintf(errmsg, "Heap contains %ld photons\n", pmap -> numPhotons); eputs(errmsg); #endif /* Allocate heap buffa */ if (!pmap -> heapBuf) { pmap -> heapBufSize = PMAP_HEAPBUFSIZE; pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)); if (!pmap -> heapBuf) error(SYSTEM, "failed to allocate postprocessed photon heap in buildPhotonMap()" ); } /* We REALLY don't need yet another @%&*! heap just to hold the scaled * photons, but can't think of a quicker fix... */ mktemp(strcpy(nuHeapFname, PMAP_TMPFNAME)); if (!(nuHeap = fopen(nuHeapFname, "w+b"))) error(SYSTEM, "failed to open postprocessed photon heap in buildPhotonMap()" ); rewind(pmap -> heap); #ifdef DEBUG_PMAP eputs("Postprocessing photons...\n"); #endif while (!feof(pmap -> heap)) { #ifdef DEBUG_PMAP printf("Reading %lu at %lu... ", pmap -> heapBufSize, ftell(pmap->heap) ); #endif pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufSize, pmap -> heap ); #ifdef DEBUG_PMAP printf("Got %lu\n", pmap -> heapBufLen); #endif if (ferror(pmap -> heap)) error(SYSTEM, "failed to read photon heap in buildPhotonMap()"); for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) { /* Update min and max pos and set photon flux */ for (i = 0; i <= 2; i++) { if (p -> pos [i] < pmap -> minPos [i]) pmap -> minPos [i] = p -> pos [i]; else if (p -> pos [i] > pmap -> maxPos [i]) pmap -> maxPos [i] = p -> pos [i]; /* Update centre of gravity with photon position */ CoG [i] += p -> pos [i]; } if (isContribPmap(pmap) && contribSrcOfs) /* Linearise contrib photon origin index from subprocess index * using the per-subprocess offsets in contribSrcOfs */ p -> aux.contribSrc += contribSrcOfs [p -> proc]; else if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) { /* Update min and max path length */ if (p -> aux.pathLen < pmap -> minPathLen) pmap -> minPathLen = p -> aux.pathLen; else if (p -> aux.pathLen > pmap -> maxPathLen) pmap -> maxPathLen = p -> aux.pathLen; /* Update avg photon path length */ avgPathLen += p -> aux.pathLen; } /* Scale photon's flux (hitherto normalised to 1 over RGB); in * case of a contrib photon map, this is done per light source, * and photonFlux is assumed to be an array */ getPhotonFlux(p, flux); if (photonFlux) { scalecolor(flux, photonFlux [isContribPmap(pmap) ? photonSrcIdx(pmap, p) : 0] ); setPhotonFlux(p, flux); } /* Update average photon flux; need a double here */ addcolor(avgFlux, flux); } /* Write modified photons to new heap */ fwrite(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufLen, nuHeap); if (ferror(nuHeap)) { error(SYSTEM, "failed postprocessing photon flux in buildPhotonMap()" ); /* Clean up temp files */ fclose(pmap -> heap); fclose(nuHeap); unlink(pmap -> heapFname); unlink(nuHeapFname); } nCheck += pmap -> heapBufLen; } + + /* !!! Flush pending writes to new heap !!! */ + fflush(nuHeap); #ifdef DEBUG_PMAP if (nCheck < pmap -> numPhotons) error(INTERNAL, "truncated photon heap in buildPhotonMap"); #endif /* Finalise average photon flux */ scalecolor(avgFlux, 1.0 / pmap -> numPhotons); copycolor(pmap -> photonFlux, avgFlux); /* Average photon positions to get centre of gravity */ for (i = 0; i < 3; i++) pmap -> CoG [i] = CoG [i] /= pmap -> numPhotons; /* Average photon path lengths */ pmap -> avgPathLen = avgPathLen /= pmap -> numPhotons; rewind(pmap -> heap); /* Compute average photon distance to centre of gravity */ while (!feof(pmap -> heap)) { pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufSize, pmap -> heap ); for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) { VSUB(d, p -> pos, CoG); CoGdist += DOT(d, d); if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) { /* Add temporal distance (in units of photon path length) for transient photons */ d [0] = p -> aux.pathLen - avgPathLen; CoGdist += d [0] * d [0]; } } } pmap -> CoGdist = CoGdist /= pmap -> numPhotons; /* Swap heaps, discarding unscaled photons */ fclose(pmap -> heap); unlink(pmap -> heapFname); pmap -> heap = nuHeap; strcpy(pmap -> heapFname, nuHeapFname); #ifdef PMAP_OOC OOC_BuildPhotonMap(pmap, nproc); #else if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) kdT_TransBuildPhotonMap(pmap); else kdT_BuildPhotonMap(pmap); #endif /* Trash heap and its buffa */ free(pmap -> heapBuf); fclose(pmap -> heap); unlink(pmap -> heapFname); pmap -> heap = NULL; pmap -> heapBuf = NULL; } /* PHOTON MAP SEARCH ROUTINES ------------------------------------------ */ /* Dynamic max photon search radius increase and reduction factors */ #define PMAP_MAXDIST_INC 4 #define PMAP_MAXDIST_DEC 0.9 /* Num successful lookups before reducing in max search radius */ #define PMAP_MAXDIST_CNT 1000 /* Threshold below which we assume increasing max radius won't help */ #define PMAP_SHORT_LOOKUP_THRESH 1 /* Coefficient for adaptive maximum search radius */ #define PMAP_MAXDIST_COEFF 100 void findPhotons (PhotonMap* pmap, const RAY* ray) { int redo = 0; const RREAL *norm = ray -> ron, *pos = ray -> rop; if (!pmap -> squeue.len) { /* Lazy init priority queue */ #ifdef PMAP_OOC OOC_InitFindPhotons(pmap); #else kdT_InitFindPhotons(pmap); #endif pmap -> minGathered = pmap -> maxGather; pmap -> maxGathered = pmap -> minGather; pmap -> totalGathered = 0; pmap -> numLookups = pmap -> numShortLookups = 0; pmap -> shortLookupPct = 0; pmap -> minError = FHUGE; pmap -> maxError = -FHUGE; pmap -> rmsError = 0; /* SQUARED max search radius limit is based on avg photon distance to * centre of gravity, unless fixed by user (maxDistFix > 0) */ pmap -> maxDist0 = pmap -> maxDist2Limit = maxDistFix > 0 ? maxDistFix * maxDistFix : PMAP_MAXDIST_COEFF * pmap -> squeue.len * pmap -> CoGdist / pmap -> numPhotons; } do { pmap -> squeue.tail = 0; pmap -> maxDist2 = pmap -> maxDist0; if (isVolumePmap(pmap) #ifdef PMAP_PHOTONFLOW || isLightFlowPmap(pmap) && sphericalIrrad #endif ) { /* Volume photons are not associated with a normal or intersection point; null the normal and set origin as lookup point. */ /* If a lightflow photon map is used and sphericalIrrad = 1, the mean spherical irradiance is evaluated, so normals are irrelevant and not filtered. */ norm = NULL; pos = ray -> rorg; } /* XXX/NOTE: status returned by XXX_FindPhotons() is currently ignored; if no photons are found, an empty queue is returned under the assumption all photons are too distant to contribute significant flux. */ #ifdef PMAP_OOC OOC_FindPhotons(pmap, pos, norm); #else if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) kdT_TransFindPhotons(pmap, pos, norm); else kdT_FindPhotons(pmap, pos, norm); #endif #ifdef PMAP_LOOKUP_INFO fprintf(stderr, "%d/%d %s photons found within radius %.3f " "at (%.2f,%.2f,%.2f) on %s\n", pmap -> squeue.tail, pmap -> squeue.len, pmapName [pmap -> type], sqrt(pmap -> maxDist2), ray -> rop [0], ray -> rop [1], ray -> rop [2], ray -> ro ? ray -> ro -> oname : "" ); #endif if (pmap -> squeue.tail < pmap -> squeue.len * pmap->gatherTolerance) { /* Short lookup; too few photons found */ if (pmap -> squeue.tail > PMAP_SHORT_LOOKUP_THRESH) { /* Ignore short lookups which return fewer than * PMAP_SHORT_LOOKUP_THRESH photons under the assumption there * really are no photons in the vicinity, and increasing the max * search radius therefore won't help */ #ifdef PMAP_LOOKUP_WARN sprintf(errmsg, "%d/%d %s photons found at (%.2f,%.2f,%.2f) on %s", pmap -> squeue.tail, pmap->squeue.len, pmapName [pmap->type], ray -> rop [0], ray -> rop [1], ray -> rop [2], ray -> ro ? ray -> ro -> oname : "" ); error(WARNING, errmsg); #endif /* Bail out after warning if maxDist is fixed */ if (maxDistFix > 0) return; if (pmap -> maxDist0 < pmap -> maxDist2Limit) { /* Increase max search radius if below limit & redo search */ pmap -> maxDist0 *= PMAP_MAXDIST_INC; #ifdef PMAP_LOOKUP_REDO redo = 1; #endif #ifdef PMAP_LOOKUP_WARN sprintf(errmsg, redo ? "restarting photon lookup with max radius %.1e" : "max photon lookup radius adjusted to %.1e", pmap -> maxDist0 ); error(WARNING, errmsg); #endif } #ifdef PMAP_LOOKUP_REDO else { sprintf(errmsg, "max photon lookup radius clamped to %.1e", pmap -> maxDist0 ); error(WARNING, errmsg); } #endif } /* Reset successful lookup counter */ pmap -> numLookups = 0; } else { /* Skip search radius reduction if maxDist is fixed */ if (maxDistFix > 0) return; /* Increment successful lookup counter and reduce max search radius if * wraparound */ pmap -> numLookups = (pmap -> numLookups + 1) % PMAP_MAXDIST_CNT; if (!pmap -> numLookups) pmap -> maxDist0 *= PMAP_MAXDIST_DEC; redo = 0; } } while (redo); } Photon *find1Photon (PhotonMap *pmap, const RAY* ray, Photon *photon) { /* Init (squared) search radius to avg photon dist to centre of gravity */ float maxDist2_0 = pmap -> CoGdist; int found = 0; #ifdef PMAP_LOOKUP_REDO #define REDO 1 #else #define REDO 0 #endif do { pmap -> maxDist2 = maxDist2_0; #ifdef PMAP_OOC found = OOC_Find1Photon(pmap, ray -> rop, ray -> ron, photon); #else found = kdT_Find1Photon(pmap, ray -> rop, ray -> ron, photon); #endif if (found < 0) { /* Expand search radius to retry */ maxDist2_0 *= 2; #ifdef PMAP_LOOKUP_WARN sprintf(errmsg, "failed 1-NN photon lookup" #ifdef PMAP_LOOKUP_REDO ", retrying with search radius %.2f", maxDist2_0 #endif ); error(WARNING, errmsg); #endif } } while (REDO && found < 0); /* Return photon buffer containing valid photon, else NULL */ return found < 0 ? NULL : photon; } void getPhoton (PhotonMap *pmap, PhotonIdx idx, Photon *photon) { #ifdef PMAP_OOC if (OOC_GetPhoton(pmap, idx, photon)) #else if (kdT_GetPhoton(pmap, idx, photon)) #endif error(INTERNAL, "failed photon lookup"); } Photon *getNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx) { #ifdef PMAP_OOC return OOC_GetNearestPhoton(squeue, idx); #else return kdT_GetNearestPhoton(squeue, idx); #endif } PhotonIdx firstPhoton (const PhotonMap *pmap) { #ifdef PMAP_OOC return OOC_FirstPhoton(pmap); #else return kdT_FirstPhoton(pmap); #endif } /* PHOTON MAP CLEANUP ROUTINES ------------------------------------------ */ void deletePhotons (PhotonMap* pmap) { #ifdef PMAP_OOC OOC_Delete(&pmap -> store); #else kdT_Delete(&pmap -> store); #endif free(pmap -> squeue.node); free(pmap -> biasCompHist); pmap -> numPhotons = pmap -> minGather = pmap -> maxGather = pmap -> squeue.len = pmap -> squeue.tail = 0; #ifdef PMAP_PATHFILT if (pmap -> pathLUT) { lu_done(pmap -> pathLUT); free(pmap -> pathLUT); } if (pmap -> pathLUTKeys) { unsigned k; for (k = 0; k < pmap->squeue.len + 1; free(pmap->pathLUTKeys [k++])); free(pmap -> pathLUTKeys); } #endif /* Contrib stuff */ if (isContribPmap(pmap) && pmap -> preCompContribTab) { + /* Parent contrib photon map; clean up child photon maps containing + precomputed contribs via lu_done() */ lu_done(pmap -> preCompContribTab); free(pmap -> preCompContribTab); pmap -> preCompContribTab = NULL; } } diff --git a/pmapdata.h b/pmapdata.h index 86d88f4..4fcff6d 100644 --- a/pmapdata.h +++ b/pmapdata.h @@ -1,436 +1,429 @@ /* 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] */ #ifdef PMAP_FLOAT_FLUX COLOR flux; #else COLR flux; #endif /* Auxiliary field; this is a multipurpose, type-specific field used by the following photon types (as identified by enum PhotonMapType in pmaptype.h): PMAP_TYPE_CONTRIB: Index into photon map's contrib origin array. PMAP_TYPE_TEMPLIGHTFLOW: Distance travelled by photon / time of flight All others: Photon path ID. */ PhotonAuxAttrib aux; } Photon; /* Define PMAP_FLOAT_FLUX to store photon flux as floats instead of * compact RGBE, which was found to improve accuracy in analytical * validation. */ #ifdef PMAP_FLOAT_FLUX #define setPhotonFlux(p,f) copycolor((p) -> flux, f) #define getPhotonFlux(p,f) copycolor(f, (p) -> flux) #else #define setPhotonFlux(p,f) setcolr((p)->flux, (f)[0], (f)[1], (f)[2]) #define getPhotonFlux(p,f) colr_color(f, (p) -> flux) #endif /* Define search queue and underlying data struct types */ #ifdef PMAP_OOC #include "pmapooc.h" #else #include "pmapkdt.h" #include "pmaptkdt.h" #endif /* Mean size of heapfile write buffer, in number of photons */ #define PMAP_HEAPBUFSIZE 1e6 /* Mean idle time between heap locking attempts, in usec */ #define PMAP_HEAPBUFSLEEP 2e6 /* Temporary filename for photon heaps */ #define PMAP_TMPFNAME TEMPLATE #define PMAP_TMPFNLEN (TEMPLEN + 1) - - /* Bias compensation history node */ - typedef struct { - COLOR irrad; - float weight; - } PhotonBiasCompNode; - - - /* Forward 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 *contribTab; /* LUT for binned contribs */ LUTAB *preCompContribTab; /* LUT for per-modifier precomp. contrib. 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 contribution photon bins prior to construction of store */ char contribHeapFname [sizeof(PMAP_TMPFNAME)]; /* ================================================================ * BIAS COMPENSATION STUFF * ================================================================ */ - PhotonBiasCompNode *biasCompHist; /* Bias compensation history */ + 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); /* Find single closest photon to ray -> rop with similar normal. Return NULL if none found, else the supplied Photon* buffer, indicating that it contains a valid photon. */ void getPhoton (PhotonMap *pmap, PhotonIdx idx, Photon *photon); /* Retrieve photon referenced by idx from pmap -> store */ Photon *getNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx); /* Retrieve photon from NN search queue after calling findPhotons() */ PhotonIdx firstPhoton (const PhotonMap *pmap); /* Index to first photon, to be passed to getPhoton(). Indices to * subsequent photons can be optained via increment operator (++) */ void deletePhotons (PhotonMap*); /* Free dem mammaries... */ #endif diff --git a/pmapio.c b/pmapio.c index 3a29abb..35e6bab 100644 --- a/pmapio.c +++ b/pmapio.c @@ -1,323 +1,383 @@ #ifndef lint static const char RCSid[] = "$Id: pmapio.c,v 2.13 2021/02/18 17:08:50 rschregle Exp $"; #endif /* ========================================================================= Photon map portable file I/O 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: pmapio.c,v 2.13 2021/02/18 17:08:50 rschregle Exp $ */ #include "pmapio.h" #include "pmapdiag.h" #include "pmapcontrib.h" #include "resolu.h" #include extern char *octname; void savePhotonMap (const PhotonMap *pmap, const char *fname, int argc, char **argv ) { - unsigned long i, j; - FILE *file; + unsigned j; + FILE *file; if (!pmap || !pmap -> numPhotons || !validPmapType(pmap -> type)) { error(INTERNAL, "attempt to save empty or invalid photon map"); return; } if (verbose) { sprintf(errmsg, "Saving %s (%ld photons)\n", fname, pmap -> numPhotons ); eputs(errmsg); fflush(stderr); } if (!(file = fopen(fname, "wb"))) { sprintf(errmsg, "can't open photon map file %s", fname); error(SYSTEM, errmsg); } /* Write header */ newheader("RADIANCE", file); /* Write command line */ printargs(argc, argv, file); /* Include statistics in info text */ fprintf(file, "NumPhotons\t= %ld\n" "AvgFlux\t\t= [%.2e, %.2e, %.2e]\n" "Bbox\t\t= [%.3f, %.3f, %.3f] [%.3f, %.3f, %.3f]\n" "CoG\t\t= [%.3f, %.3f, %.3f]\n" "MaxDist^2\t= %.3f\n" "Velocity\t= %g / sec\n" "Timespan\t= [%e, %e] sec\n" "AvgPathLen\t= %.3f\n", pmap -> numPhotons, pmap -> photonFlux [0], pmap -> photonFlux [1], pmap -> photonFlux [2], pmap -> minPos [0], pmap -> minPos [1], pmap -> minPos [2], pmap -> maxPos [0], pmap -> maxPos [1], pmap -> maxPos [2], pmap -> CoG [0], pmap -> CoG [1], pmap -> CoG [2], pmap -> CoGdist, pmap -> velocity, pmap -> minPathLen / pmap -> velocity, pmap -> maxPathLen / pmap -> velocity, pmap -> avgPathLen ); +#ifdef PMAP_OOC + if (isContribChildPmap(pmap) && pmap -> preCompContrib) { + /* Child contrib pmap containing precomputed contributions; + append binning params to header */ + const PreComputedContrib *preCompContrib = (PreComputedContrib*)( + pmap -> preCompContrib + ); + const unsigned long nBins = (unsigned long)preCompContrib -> scDim * + preCompContrib -> scDim; + + fprintf(file, + "Bin resolution\t= %d\n" + "Num bins\t= %d x %d\n" + "Num compressed\t= %d/%d\n" + "Compression\t= %.1f%%\n", + preCompContrib -> l, + preCompContrib -> scDim, preCompContrib -> scDim, + preCompContrib -> nCompressedBins, nBins, + 100 * (1 - (float)preCompContrib -> nCompressedBins / nBins) + ); + } +#endif + /* Write format, including human-readable file version */ fputformat((char *)pmapFormat [pmap -> type], file); fprintf(file, "VERSION=%s\n", PMAP_FILEVER); /* Empty line = end of header */ putc('\n', file); /* Write machine-readable file format version */ putstr(PMAP_FILEVER, file); /* Write number of photons as fixed size, which possibly results in * padding of MSB with 0 on some platforms. Unlike sizeof() however, * this ensures portability since this value may span 32 or 64 bits * depending on platform. */ putint(pmap -> numPhotons, PMAP_LONGSIZE, file); /* Write average photon flux */ for (j = 0; j < 3; j++) putflt(pmap -> photonFlux [j], file); /* Write max and min photon positions */ for (j = 0; j < 3; j++) { putflt(pmap -> minPos [j], file); putflt(pmap -> maxPos [j], file); } /* Write centre of gravity */ for (j = 0; j < 3; j++) putflt(pmap -> CoG [j], file); /* Write avg distance to centre of gravity */ putflt(pmap -> CoGdist, file); /* Save photon's velocity (speed of light in units of scene geometry), min/max/avg photon path length (= timespan and temporal CoG) */ putflt(pmap -> velocity, file); putflt(pmap -> minPathLen, file); putflt(pmap -> maxPathLen, file); putflt(pmap -> avgPathLen, file); + if (isContribChildPmap(pmap) && pmap -> preCompContrib) { + /* Child contrib photon map containing per-modifier precomputed + contributions; save bin resolution and num compressed bins */ + const PreComputedContrib *preCompContrib = (PreComputedContrib*)( + pmap -> preCompContrib + ); + + putint(preCompContrib -> l, sizeof(preCompContrib -> l), file); + putint(preCompContrib -> nCompressedBins, + sizeof(preCompContrib -> nCompressedBins), file + ); + } + /* Save photon storage */ - #ifdef PMAP_OOC - if (isContribPmap(pmap) && pmap -> preCompContribTab) - /* Save constituent per-modifier precomputed contribution pmaps */ - saveContribPhotonMap(pmap, argc, argv); - else if (OOC_SavePhotons(pmap, file)) { - #else - if (kdT_SavePhotons(pmap, file)) { - #endif - sprintf(errmsg, "error writing photon map file %s", fname); - error(SYSTEM, errmsg); - } +#ifdef PMAP_OOC + if (isContribPmap(pmap) && pmap -> preCompContribTab) + /* Parent contrib photon map; save child photon maps containing + per-modifier precomputed contributions */ + saveContribPhotonMap(pmap, argc, argv); + else if (OOC_SavePhotons(pmap, file)) { +#else + if (kdT_SavePhotons(pmap, file)) { +#endif + sprintf(errmsg, "error writing photon map file %s", fname); + error(SYSTEM, errmsg); + } fclose(file); } PhotonMapType loadPhotonMap (PhotonMap *pmap, const char *fname) { PhotonMapType ptype = PMAP_TYPE_NONE; char format [MAXFMTLEN]; unsigned long i, j; FILE *file; if (!pmap) return PMAP_TYPE_NONE; if ((file = fopen(fname, "rb")) == NULL) { sprintf(errmsg, "can't open photon map file %s", fname); error(SYSTEM, errmsg); } /* Get format string */ strcpy(format, PMAP_FORMAT_GLOB); if (checkheader(file, format, NULL) != 1) { sprintf(errmsg, "photon map file %s has an unknown format", fname); error(USER, errmsg); } /* Identify photon map type from format string */ for (ptype = 0; ptype < NUM_PMAP_TYPES && strcmp(pmapFormat [ptype], format); ptype++ ); if (!validPmapType(ptype)) { sprintf(errmsg, "file %s contains an unknown photon map type", fname); error(USER, errmsg); } initPhotonMap(pmap, ptype); /* Get file format version and check for compatibility */ if (strcmp(getstr(format, file), PMAP_FILEVER)) error(USER, "incompatible photon map file format"); /* Get number of photons as fixed size, which possibly results in * padding of MSB with 0 on some platforms. Unlike sizeof() however, * this ensures portability since this value may span 32 or 64 bits * depending on platform. */ pmap -> numPhotons = getint(PMAP_LONGSIZE, file); /* Get average photon flux */ for (j = 0; j < 3; j++) pmap -> photonFlux [j] = getflt(file); /* Get max and min photon positions */ for (j = 0; j < 3; j++) { pmap -> minPos [j] = getflt(file); pmap -> maxPos [j] = getflt(file); } /* Get centre of gravity */ for (j = 0; j < 3; j++) pmap -> CoG [j] = getflt(file); /* Get avg distance to centre of gravity */ pmap -> CoGdist = getflt(file); /* Get photon's velocity (speed of light in units of scene geometry), min/max/avg photon path length (= timespan) and temporal CoG */ pmap -> velocity = getflt(file); pmap -> minPathLen = getflt(file); pmap -> maxPathLen = getflt(file); pmap -> avgPathLen = getflt(file); + if (isContribChildPmap(pmap)) { + /* Per-modifier precomputed contribution (child) photon map; + allocate precomputed contribs, get bin resolution and num + compressed bins */ + PreComputedContrib *preCompContrib = malloc( + sizeof(PreComputedContrib) + ); + + if (!preCompContrib) { + sprintf(errmsg, + "cannot alloc precomputed contribs in photon map %s", fname + ); + error(SYSTEM, errmsg); + } + + pmap -> preCompContrib = (struct PreComputedContrib*)preCompContrib; + preCompContrib -> l = getint(sizeof(preCompContrib -> l), file); + preCompContrib -> scDim = PMAP_CONTRIB_SCDIM(preCompContrib -> l); + preCompContrib -> nCompressedBins = getint( + sizeof(preCompContrib -> nCompressedBins), file + ); + } + /* Load photon storage */ #ifdef PMAP_OOC if (OOC_LoadPhotons(pmap, file)) { #else if (kdT_LoadPhotons(pmap, file)) { #endif sprintf(errmsg, "error reading photon map file %s", fname); error(SYSTEM, errmsg); } fclose(file); return ptype; } void savePmaps (const PhotonMap **pmaps, int argc, char **argv) /* Wrapper to save all defined photon maps with specified command line */ { unsigned t; for (t = 0; t < NUM_PMAP_TYPES; t++) { if (pmaps [t]) savePhotonMap(pmaps [t], pmaps [t] -> fileName, argc, argv); } } void loadPmaps (PhotonMap **pmaps, const PhotonMapParams *parm) /* Wrapper to save all defined photon maps with specified command line */ { unsigned t; struct stat octstat, pmstat; PhotonMap *pm; PhotonMapType type; for (t = 0; t < NUM_PMAP_TYPES; t++) if (setPmapParam(&pm, parm + t)) { /* Check if photon map newer than octree */ if ( pm -> fileName && octname && !stat(pm -> fileName, &pmstat) && !stat(octname, &octstat) && octstat.st_mtime > pmstat.st_mtime ) { sprintf( errmsg, "photon map in file %s may be stale", pm -> fileName ); error(USER, errmsg); } /* Load photon map from file and get its type */ if ((type = loadPhotonMap(pm, pm -> fileName)) == PMAP_TYPE_NONE) error(USER, "failed loading photon map"); /* Assign to appropriate photon map type (deleting previously * loaded photon map of same type if necessary) */ if (pmaps [type]) { sprintf( errmsg, "multiple %s photon maps, dropping previous", pmapName [type] ); error(WARNING, errmsg); deletePhotons(pmaps [type]); free(pmaps [type]); } pmaps [type] = pm; /* Check for valid density estimate bandwidths */ if ( (pm -> minGather > 1 || pm -> maxGather > 1) && (type == PMAP_TYPE_PRECOMP) ) { /* Force bwidth to 1 for precomputed pmap */ error(WARNING, "ignoring bandwidth for precomp photon map"); pm -> minGather = pm -> maxGather = 1; } if ( (pm -> maxGather > pm -> minGather) && (type == PMAP_TYPE_VOLUME) ) { /* Biascomp for volume pmaps (see volumePhotonDensity() below) is considered redundant, and there's probably no point in recovering by using the lower bandwidth, since it's probably not what the user wants, so bail out. */ sprintf( errmsg, "bias compensation is not available with %s photon maps", pmapName [type] ); error(USER, errmsg); } if (pm -> maxGather > pm -> numPhotons) { /* Clamp lookup bandwidth to total number of photons (minus one, since density estimate gets extra photon to obtain averaged radius) */ sprintf( errmsg, "clamping density estimate bandwidth to %ld", pm -> numPhotons ); error(WARNING, errmsg); pm -> minGather = pm -> maxGather = pm -> numPhotons - 1; } } } diff --git a/pmapooc.c b/pmapooc.c index 864146e..bd20ee3 100644 --- a/pmapooc.c +++ b/pmapooc.c @@ -1,356 +1,362 @@ /* ====================================================================== Photon map interface to out-of-core octree Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ====================================================================== $Id: pmapooc.c,v 1.7 2021/03/22 23:00:00 rschregle Exp $ */ #include "pmapdata.h" /* Includes pmapooc.h */ #include "source.h" #include "otspecial.h" #include "oocsort.h" #include "oocbuild.h" #include "pmcontribsort.h" /* PHOTON MAP BUILD ROUTINES ------------------------------------------ */ -/* Returns photon position as sorting key for OOC_Octree & friends (notably - * for Morton code generation). +/* Returns photon position as sorting key for OOC_Octree & friends + * (notably for Morton code generation). * !!! XXX: Uses type conversion from float via TEMPORARY storage; * !!! THIS IS NOT THREAD SAFE! * !!! RETURNED KEY PERSISTS ONLY IF COPIED BEFORE NEXT CALL! */ RREAL *OOC_PhotonKey (const void *p) { static FVECT photonPos; /* Temp storage for type conversion */ VCOPY(photonPos, ((Photon*)p) -> pos); return photonPos; } #ifdef DEBUG_OOC static void OOC_CheckKeys (FILE *file, const OOC_Octree *oct) { Photon p, lastp; RREAL *key; MortonIdx zIdx, lastzIdx = 0; rewind(file); memset(&lastp, 0, sizeof(lastp)); while (fread(&p, sizeof(p), 1, file) > 0) { key = OOC_PhotonKey(&p); zIdx = OOC_KEY2MORTON(key, oct); if (zIdx < lastzIdx) error(INTERNAL, "photons not sorted"); if (zIdx == lastzIdx) { sprintf(errmsg, "identical key %021ld " "for [%.9f, %.9f, %.9f]\tand [%.9f, %.9f, %.9f]", zIdx, lastp.pos [0], lastp.pos [1], lastp.pos [2], p.pos [0], p.pos [1], p.pos [2] ); error(WARNING, errmsg); } lastzIdx = zIdx; memcpy(&lastp, &p, sizeof(p)); } } #endif void OOC_BuildPhotonMap (struct PhotonMap *pmap, unsigned numProc) { FILE *leafFile; - char leafFname [1024]; + char *leafFname; FVECT d, octOrg; int i; RREAL octSize = 0; /* Determine octree size and origin from pmap extent and init octree */ VCOPY(octOrg, pmap -> minPos); VSUB(d, pmap -> maxPos, pmap -> minPos); for (i = 0; i < 3; i++) if (octSize < d [i]) octSize = d [i]; if (octSize < FTINY) error(INTERNAL, "zero octree size in OOC_BuildPhotonMap"); /* Derive leaf filename from photon map and open file */ - strncpy(leafFname, pmap -> fileName, sizeof(leafFname)); - strncat(leafFname, PMAP_OOC_LEAFSUFFIX, - sizeof(leafFname) - strlen(leafFname) - 1 + leafFname = malloc( + strlen(pmap -> fileName) + strlen(PMAP_OOC_LEAFSUFFIX) + 1 ); + if (!leafFname) + error(SYSTEM, "failed leaf filename alloc in OOC_BuildPhotonMap"); + + strcpy(leafFname, pmap -> fileName); + strcat(leafFname, PMAP_OOC_LEAFSUFFIX); if (!(leafFile = fopen(leafFname, "w+b"))) error(SYSTEM, "failed to open leaf file in OOC_BuildPhotonMap"); #ifdef DEBUG_OOC eputs("Sorting photons by Morton code...\n"); #endif /* Sort photons in heapfile by Morton code and write to leaf file; - redirect to contrib sorting routine for contribution photons! */ - if (isContribPmap(pmap) + redirect to contrib sorting routine if contrib child pmap and + we have precomputed contributions, else we need a regular sort + before we can precompute */ + if (isContribChildPmap(pmap) && pmap -> preCompContrib ? contribPhotonSort(pmap, leafFile, octOrg, octSize, PMAP_OOC_NUMBLK, PMAP_OOC_BLKSIZE, numProc ) : OOC_Sort(pmap -> heap, leafFile, PMAP_OOC_NUMBLK, PMAP_OOC_BLKSIZE, numProc, sizeof(Photon), octOrg, octSize, OOC_PhotonKey ) ) error(INTERNAL, "failed out-of-core photon sort in OOC_BuildPhotonMap"); /* Init and build octree */ OOC_Init(&pmap -> store, sizeof(Photon), octOrg, octSize, OOC_PhotonKey, - leafFile + leafFile, leafFname ); #ifdef DEBUG_OOC eputs("Checking leaf file consistency...\n"); OOC_CheckKeys(leafFile, &pmap -> store); eputs("Building out-of-core octree...\n"); #endif if (!OOC_Build(&pmap -> store, PMAP_OOC_LEAFMAX, PMAP_OOC_MAXDEPTH)) error(INTERNAL, "failed out-of-core octree build in OOC_BuildPhotonMap"); #ifdef DEBUG_OOC eputs("Checking out-of-core octree consistency...\n"); if (OOC_Check(&pmap -> store, OOC_ROOT(&pmap -> store), octOrg, octSize, 0 )) error(INTERNAL, "inconsistent out-of-core octree; Time4Harakiri"); #endif } /* PHOTON MAP I/O ROUTINES ------------------------------------------ */ int OOC_SavePhotons (const struct PhotonMap *pmap, FILE *out) { return OOC_SaveOctree(&pmap -> store, out); } int OOC_LoadPhotons (struct PhotonMap *pmap, FILE *nodeFile) { FILE *leafFile; char leafFname [1024]; /* Derive leaf filename from photon map and open file */ strncpy(leafFname, pmap -> fileName, sizeof(leafFname)); strncat(leafFname, PMAP_OOC_LEAFSUFFIX, sizeof(leafFname) - strlen(leafFname) - 1 ); if (!(leafFile = fopen(leafFname, "r"))) error(SYSTEM, "failed to open leaf file in OOC_LoadPhotons"); if (OOC_LoadOctree(&pmap -> store, nodeFile, OOC_PhotonKey, leafFile)) return -1; #ifdef DEBUG_OOC /* Check octree for consistency */ if (OOC_Check(&pmap -> store, OOC_ROOT(&pmap -> store), pmap -> store.org, pmap -> store.size, 0 )) return -1; #endif return 0; } /* PHOTON MAP SEARCH ROUTINES ------------------------------------------ */ void OOC_InitFindPhotons (struct PhotonMap *pmap) { if (OOC_InitNearest(&pmap -> squeue, pmap -> maxGather + 1, sizeof(Photon) )) error(SYSTEM, "can't allocate photon search queue"); #ifdef PMAP_PATHFILT initPhotonPaths(pmap); #endif } static void OOC_InitPhotonCache (struct PhotonMap *pmap) /* Initialise OOC photon cache */ { static char warn = 1; if (!pmap -> store.cache && !pmap -> numDensity) { if (pmapCacheSize > 0) { const unsigned pageSize = pmapCachePageSize * pmap -> maxGather, numPages = pmapCacheSize / pageSize; /* Allocate & init I/O cache in octree */ pmap -> store.cache = malloc(sizeof(OOC_Cache)); if (!pmap -> store.cache || OOC_CacheInit(pmap -> store.cache, numPages, pageSize, sizeof(Photon) )) { error(SYSTEM, "failed OOC photon map cache init"); } } else if (warn) { error(WARNING, "OOC photon map cache DISABLED"); warn = 0; } } } int OOC_FindPhotons (struct PhotonMap *pmap, const FVECT pos, const FVECT norm ) { OOC_SearchFilterCallback filt; OOC_SearchFilterState filtState; #ifdef PMAP_PATHFILT OOC_SearchAttribCallback paths, *pathsPtr = NULL; #endif float res, n [3]; /* Lazily init OOC cache */ if (!pmap -> store.cache) OOC_InitPhotonCache(pmap); /* Set up filter callback */ if (norm) VCOPY(n, norm); filtState.pmap = pmap; filtState.pos = pos; filtState.norm = norm ? n : NULL; filt.state = &filtState; filt.filtFunc = filterPhoton; /* Get sought light source modifier from pmap -> lastContribSrc if precomputing contribution photon */ filtState.srcMod = isContribPmap(pmap) ? findmaterial(source [pmap -> lastContribSrc.srcIdx].so) : NULL; #ifdef PMAP_PATHFILT /* Set up path ID callback to filter for photons from unique paths */ paths.state = pmap; paths.findFunc = findPhotonPath; paths.addFunc = addPhotonPath; paths.delFunc = deletePhotonPath; paths.checkFunc = checkPhotonPaths; pathsPtr = &paths; resetPhotonPaths(pmap); res = OOC_FindNearest( &pmap -> store, OOC_ROOT(&pmap -> store), 0, pmap -> store.org, pmap -> store.size, pos, &filt, pathsPtr, &pmap -> squeue, pmap -> maxDist2 ); #else res = OOC_FindNearest( &pmap -> store, OOC_ROOT(&pmap -> store), 0, pmap -> store.org, pmap -> store.size, pos, &filt, NULL, &pmap -> squeue, pmap -> maxDist2 ); #endif /* PMAP_PATHFILT */ if (res < 0) error(INTERNAL, "failed k-NN photon lookup in OOC_FindPhotons"); /* Get (maximum distance)^2 from search queue head */ pmap -> maxDist2 = pmap -> squeue.node [0].dist2; /* Return success or failure (empty queue => none found) */ return pmap -> squeue.tail ? 0 : -1; } int OOC_Find1Photon (struct PhotonMap* pmap, const FVECT pos, const FVECT norm, Photon *photon ) { OOC_SearchFilterCallback filt; OOC_SearchFilterState filtState; float n [3], maxDist2; /* Lazily init OOC cache */ if (!pmap -> store.cache) OOC_InitPhotonCache(pmap); /* Set up filter callback */ if (norm) VCOPY(n, norm); filtState.pmap = pmap; filtState.norm = norm ? n : NULL; filtState.srcMod = NULL; filt.state = &filtState; filt.filtFunc = filterPhoton; maxDist2 = OOC_Find1Nearest(&pmap -> store, OOC_ROOT(&pmap -> store), 0, pmap -> store.org, pmap -> store.size,pos, &filt, photon, pmap -> maxDist2 ); if (maxDist2 < 0) error(INTERNAL, "failed 1-NN photon lookup in OOC_Find1Photon"); if (maxDist2 >= pmap -> maxDist2) /* No photon found => failed */ return -1; else { /* Set photon distance => success */ pmap -> maxDist2 = maxDist2; return 0; } } /* PHOTON ACCESS ROUTINES ------------------------------------------ */ int OOC_GetPhoton (struct PhotonMap *pmap, PhotonIdx idx, Photon *photon) { return OOC_GetData(&pmap -> store, idx, photon); } Photon *OOC_GetNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx) { return OOC_GetNearest(squeue, idx); } PhotonIdx OOC_FirstPhoton (const struct PhotonMap* pmap) { return 0; } diff --git a/pmaptype.c b/pmaptype.c index 3e18c39..8f2981d 100644 --- a/pmaptype.c +++ b/pmaptype.c @@ -1,55 +1,57 @@ #ifndef lint static const char RCSid[] = "$Id: pmaptype.c,v 2.5 2016/05/17 17:39:47 rschregle Exp $"; #endif /* ========================================================================= Photon map types and corresponding file format strings 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: pmaptype.c,v 2.5 2016/05/17 17:39:47 rschregle Exp $ */ #include "pmaptype.h" #ifdef PMAP_OOC #define PMAP_FMTSUFFIX "OOC_Photon_Map" #else #define PMAP_FMTSUFFIX "kdT_Photon_Map" #endif /* Format strings for photon map files corresponding to PhotonMapType */ const char *pmapFormat [NUM_PMAP_TYPES] = { "Radiance_Global_" PMAP_FMTSUFFIX, "Radiance_PreComp_" PMAP_FMTSUFFIX, "Radiance_Caustic_" PMAP_FMTSUFFIX, "Radiance_Volume_" PMAP_FMTSUFFIX, "Radiance_Direct_" PMAP_FMTSUFFIX, "Radiance_Contrib_" PMAP_FMTSUFFIX, - "Radiance_Transient_" PMAP_FMTSUFFIX + "Radiance_Transient_" PMAP_FMTSUFFIX, #ifdef PMAP_PHOTONFLOW - , "Radiance_LightFlow_" PMAP_FMTSUFFIX, - "Radiance_TransLightFlow_" PMAP_FMTSUFFIX + "Radiance_LightFlow_" PMAP_FMTSUFFIX, + "Radiance_TransLightFlow_" PMAP_FMTSUFFIX, #endif + "Radiance_Contrib_Child_" PMAP_FMTSUFFIX }; /* Photon map names per type */ const char *pmapName [NUM_PMAP_TYPES] = { "global", "precomp", "caustic", "volume", "direct", "contrib", - "transient" + "transient", #ifdef PMAP_PHOTONFLOW - , "lightflow", "transient lightflow" + "lightflow", "transient lightflow", #endif + "contrib child" }; diff --git a/pmaptype.h b/pmaptype.h index 8fdea06..d315e7e 100644 --- a/pmaptype.h +++ b/pmaptype.h @@ -1,56 +1,57 @@ /* RCSid $Id: pmaptype.h,v 2.5 2016/05/17 17:39:47 rschregle Exp $ */ /* ========================================================================= Photon map types and corresponding file format strings 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: pmaptype.h,v 2.5 2016/05/17 17:39:47 rschregle Exp $ */ /* TODO: Combine type-specific properties of pmap types as binary flags, e.g. PMAP_TYPE_TRANSLIGHTFLOW = PMAP_TYPE_VOLUME | PMAP_TYPE_TRANSIENT | PMAP_TYPE_LIGHTFLOW ??? */ #ifndef PMAPTYPE_H #define PMAPTYPE_H /* Photon map types */ typedef enum { - PMAP_TYPE_NONE = -1, PMAP_TYPE_GLOBAL, PMAP_TYPE_PRECOMP, - PMAP_TYPE_CAUSTIC, PMAP_TYPE_VOLUME, PMAP_TYPE_DIRECT, + PMAP_TYPE_NONE = -1, PMAP_TYPE_GLOBAL, PMAP_TYPE_PRECOMP, + PMAP_TYPE_CAUSTIC, PMAP_TYPE_VOLUME, PMAP_TYPE_DIRECT, PMAP_TYPE_CONTRIB, PMAP_TYPE_TRANSIENT, #ifdef PMAP_PHOTONFLOW PMAP_TYPE_LIGHTFLOW, PMAP_TYPE_TRANSLIGHTFLOW, #endif + PMAP_TYPE_CONTRIB_CHILD = 9, NUM_PMAP_TYPES } PhotonMapType; /* Check for valid photon map type */ #define validPmapType(t) ((t) >= 0 && (t) < NUM_PMAP_TYPES) /* Glob string for extracting photon map format from file header */ #define PMAP_FORMAT_GLOB "Radiance_*_Photon_Map" /* Format strings for photon map files corresponding to PhotonMapType */ extern const char *pmapFormat []; /* Photon map names per type */ extern const char *pmapName []; #endif diff --git a/pmcontrib3.c b/pmcontrib3.c index 268dc89..577321c 100644 --- a/pmcontrib3.c +++ b/pmcontrib3.c @@ -1,208 +1,195 @@ #ifndef lint static const char RCSid[] = "$Id$"; #endif /* ========================================================================= Photon map routines for precomputed light source contributions. These routines build and save precomputed contribution photon maps, 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$ */ #include "pmapcontrib.h" #include "pmapdiag.h" #include "pmapio.h" #include #if NIX #define __USE_XOPEN_EXTENDED /* Enables nftw() under Linux */ #include #endif #if NIX static int ftwFunc (const char *path, const struct stat *statData, int typeFlags, struct FTW *ftwData ) /* Callback for nftw() to delete a directory entry */ { return remove(path); /* = unlink()/rmdir() depending on type */ } #endif int buildPreCompContribPmap (const LUENT *preCompContribNode, void *p) /* Build per-modifier precomputed photon map. Returns 0 on success. */ { PhotonMap *preCompContribPmap = (PhotonMap*)( preCompContribNode -> data ); PreComputedContrib *preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib ); char dirName [1024]; struct stat dirStat; /* Flush photon/contrib heaps */ flushPhotonHeap(preCompContribPmap); fflush(preCompContribPmap -> contribHeap); if (verbose) { sprintf(errmsg, "Building precomputed contribution photon map " "for modifier %s\n", preCompContribNode -> key ); eputs(errmsg); #if NIX fflush(stderr); #endif } /* Set output directory for filename */ snprintf(dirName, sizeof(dirName), PMAP_CONTRIB_DIR, preCompContribPmap -> fileName ); /* Allocate & set output filenames to subdirectory and modifier */ preCompContribPmap -> fileName = malloc(strlen(dirName) + strlen(preCompContribNode -> key) + strlen(PMAP_CONTRIB_FILE) ); preCompContrib -> waveletFname = malloc(strlen(dirName) + strlen(preCompContribNode -> key) + strlen(PMAP_CONTRIB_WAVELETFILE) ); if (!preCompContribPmap -> fileName || !preCompContrib -> waveletFname) error(SYSTEM, "out of memory allocating filename in " "buildPreCompContribPmap()" ); snprintf(preCompContribPmap -> fileName, 1028, PMAP_CONTRIB_FILE, dirName, preCompContribNode -> key ); snprintf(preCompContrib -> waveletFname, 1029, PMAP_CONTRIB_WAVELETFILE, dirName, preCompContribNode -> key ); /* Cleanup previous directory contents if necessary */ if (!stat(dirName, &dirStat)) { /* File/dir exists */ if (S_ISREG(dirStat.st_mode)) /* Found regular file; delete and skip rest */ unlink(dirName); else if (S_ISDIR(dirStat.st_mode)) { /* Found subdirectory; delete its contents, skipping symlinks */ #if NIX if (nftw(dirName, ftwFunc, 1, FTW_DEPTH | FTW_MOUNT | FTW_PHYS) < 0 ) { sprintf(errmsg, "failed cleanup of output directory %s", dirName ); error(SYSTEM, errmsg); } #else /* Apparently there's no nftw() under Wind0z, so just skip * cleanup. Tuff luck, Wind0z Weenie! There's probably some * replacement for this, but we couldn't be buggered... */ #endif } else { /* Found neither regular file nor directory; whatever it is, * just stuff it and quit! */ sprintf(errmsg, "cannot remove existing output file %s", dirName ); error(SYSTEM, errmsg); } } /* Create new directory for per-modifier contribution photon maps */ if (mkdir(dirName, S_IFDIR | S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH) < 0 ) { sprintf(errmsg, "error creating output directory %s", dirName); error(SYSTEM, errmsg); } /* Open wavelet coefficient file; this is where the sorted contribs are written to by contribPhotonSort() (via buildPhotonMap()) */ if (!(preCompContrib -> waveletFile = fopen( preCompContrib -> waveletFname, "wb" )) ) { sprintf(errmsg, "can't open wavelet coefficient file %s", preCompContrib -> waveletFname ); error(SYSTEM, errmsg); } /* Rebuild underlying data structure, destroying heap, and sort contribs into wavelet coeff file */ buildPhotonMap(preCompContribPmap, NULL, NULL, 1); /* Free primary and transposed wavelet matrices; no longer needed */ free(preCompContrib -> waveletMatrix); freeWaveletMatrix(preCompContrib -> tWaveletMatrix, preCompContrib -> l); /* Free thresholded coefficients; no longer needed */ free(preCompContrib -> threshCoeffs); preCompContrib -> waveletMatrix = preCompContrib -> tWaveletMatrix = NULL; preCompContrib -> threshCoeffs = NULL; return 0; } static int savePreCompContrib (const LUENT *preCompContribNode, void *p) { PhotonMap *preCompContribPmap = (PhotonMap*)( preCompContribNode -> data ); - PreComputedContrib *preCompContrib = (PreComputedContrib*)( - preCompContribPmap -> preCompContrib - ); PhotonMapArgz *argz = (PhotonMapArgz*)p; /* Save child contrib photon map for current modifier; savePhotonMap() * will not recurse and call saveContribPhotonMap() again because * preCompContribPmap -> preCompContribTab == NULL */ savePhotonMap(preCompContribPmap, preCompContribPmap -> fileName, argz -> argc, argz -> argv ); - -#if 0 - /* Save wavelet coefficients */ - if (!(waveletFile = fopen(preCompContrib -> waveletFname, "wb"))) { - sprintf(errmsg, "can't open wavelet coefficient file %s", - preCompContrib -> waveletFname - ); - error(SYSTEM, errmsg); - } -#endif - + return 0; } void saveContribPhotonMap (const PhotonMap *pmap, int argc, char **argv) { PhotonMapArgz argz = {argc, argv}; lu_doall(pmap -> preCompContribTab, savePreCompContrib, &argz); } diff --git a/pmcontribsort.c b/pmcontribsort.c index eff2303..b6f0552 100644 --- a/pmcontribsort.c +++ b/pmcontribsort.c @@ -1,510 +1,536 @@ #ifndef lint static const char RCSid[] = "$Id$"; #endif /* ========================================================================= Support routines for N-way out-of-core merge sort of precomputed photon contributions. These routines are used adapted from OOC_Sort() to sort binned contributions in the same order as their associated precomputed photons. 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$ */ #if !defined(_WIN32) && !defined(_WIN64) || defined(PMAP_OOC) /* No Windoze 'cos no fork() */ #include "oocsort.h" #include "pmutil.h" #include "pmapcontrib.h" #include #include #include #include #include #define PHOTONSIZE sizeof(Photon) #define PMAPCONTRIB_KEY(p, kdat) Key2Morton3D( \ kdat -> key(p), kdat -> bbOrg, kdat -> mortonScale \ ) extern RREAL *OOC_PhotonKey (const void *p); static void contribPhotonQuickSort (Photon *photons, char *contribs, unsigned contribSize, const OOC_KeyData *keyData, unsigned long left, unsigned long right ) /* Partition photons and their binned contributions using a quicksort * algorithm. Returns index to pivot's position after partitioning. */ { - const unsigned long m = (right - left + 1) >> 1; /* Pivot in mid */ - unsigned long l = left, r = right; - const MortonIdx mKey = PMAPCONTRIB_KEY(photons + m, keyData); - char tmp [max(PHOTONSIZE, contribSize)]; + unsigned long l, r, m; + MortonIdx lKey, rKey, mKey, tKey; + + #define PMAPCONTRIB_KEYSWAP(i,j) (tKey = (i), (i) = (j), (j) = tKey); + void contribPhotonSwap (unsigned long i, unsigned long j) + { + /* TODO: swap photon indices instead of photons themselves? */ + char tmp [max(PHOTONSIZE, contribSize)]; + + memcpy(tmp, photons + i, PHOTONSIZE); + memcpy(photons + i, photons + j, PHOTONSIZE); + memcpy(photons + j, tmp, PHOTONSIZE); + memcpy(tmp, contribs + i, contribSize); + memcpy(contribs + i, contribs + j, contribSize); + memcpy(contribs + j, tmp, contribSize); + } + if (left < right) { + /* Select pivot from median-of-three and move to photons [right] + (a.k.a. Lomuto partitioning) */ + l = left; + r = right; + m = l + ((r - l) >> 1); /* Avoids overflow vs. (l+r) >> 1 */ + lKey = PMAPCONTRIB_KEY(photons + l, keyData); + rKey = PMAPCONTRIB_KEY(photons + r, keyData); + mKey = PMAPCONTRIB_KEY(photons + m, keyData); + if (mKey < lKey) { + contribPhotonSwap(l, m); + PMAPCONTRIB_KEYSWAP(lKey, mKey); + } + if (rKey < lKey) { + contribPhotonSwap(l, r); + PMAPCONTRIB_KEYSWAP(lKey, rKey); + } + if (mKey < rKey) { + contribPhotonSwap(m, r); + PMAPCONTRIB_KEYSWAP(mKey, rKey); + } + /* Pivot with key rKey is now in photons [right] */ + /* l & r converge, swapping photons (and their corresponding - * binned contributions) out of order with respect to pivot at - * photons [m]. The convergence point l = r is the pivot's final - * position */ - do { - while (l <= r && PMAPCONTRIB_KEY(photons + l, keyData) < mKey) + binned contributions) out of order with respect to pivot + (a.k.a. Hoare partitioning). The convergence point l = r + is the pivot's final position */ + while (l < r) { + while (l < r && PMAPCONTRIB_KEY(photons + l, keyData) < rKey) l++; - - while (r > l && PMAPCONTRIB_KEY(photons + r, keyData) >= mKey) + while (r > l && PMAPCONTRIB_KEY(photons + r, keyData) >= rKey) r--; - - /* Photons out of order; swap them and their corresponding - * binned contributions */ - /* TODO: swap photon indices instead of photons themselves? */ - memcpy(tmp, photons + l, PHOTONSIZE); - memcpy(photons + l, photons + r, PHOTONSIZE); - memcpy(photons + r, tmp, PHOTONSIZE); - memcpy(tmp, contribs + l, contribSize); - memcpy(contribs + l, contribs + r, contribSize); - memcpy(contribs + r, tmp, contribSize); - } while (l < r); + + /* Photons out of order; swap them and their contributions */ + if (l < r) + contribPhotonSwap(l, r); + }; - /* Now l = r; swap this photon with pivot, ditto for contribs */ - memcpy(tmp, photons + r, PHOTONSIZE); - memcpy(photons + r, photons + m, PHOTONSIZE); - memcpy(photons + m, tmp, PHOTONSIZE); - memcpy(tmp, contribs + r, contribSize); - memcpy(contribs + r, contribs + m, contribSize); - memcpy(contribs + m, tmp, contribSize); + /* Now l = r = pivot's final position; swap these photons */ + lKey = PMAPCONTRIB_KEY(photons + l, keyData); + contribPhotonSwap(l, right); - /* Recurse in left & right partitions */ - contribPhotonQuickSort(photons, contribs, contribSize, keyData, - left, m - 1 - ); - contribPhotonQuickSort(photons, contribs, contribSize, keyData, - m + 1, right - ); + /* Recurse in partitions flanking pivot */ + if (l > left) + contribPhotonQuickSort(photons, contribs, contribSize, + keyData, left, l - 1 + ); + if (l < right) + contribPhotonQuickSort(photons, contribs, contribSize, + keyData, l + 1, right + ); } } /* Active subprocess counter updated by parent process; must persist when * recursing into contribPhotonSortRecurse(), hence global */ static unsigned procCnt = 0; static int contribPhotonSortRecurse ( FILE *photonIn, FILE *contribIn, FILE* photonOut, FILE* contribOut, unsigned contribSize, unsigned numBlk, unsigned numProc, Photon *photonBuf, char *contribBuf, unsigned long bufLen, const OOC_KeyData *keyData, unsigned long blkLo, unsigned long blkHi ) /* Recursive part of contribPhotonSort(). Reads block of photons and * their binned contributions from input files photonIn resp. * contribOut, within the interval [blkLo, blkHi], and divides them into * numBlk blocks until the number of photons/contributions in a block * does not exceed bufLen, then quicksorts each block in-core into * temporary files. These files are then mergesorted via a priority * queue to the specified output files as the stack unwinds. * NOTE: Critical sections are prepended with '!!!' * * Parameters are as follows: * photonIn Opened primary input file containing unsorted photons; * their Morton codes determine the ordering * contribIn Opened secondary input file containing unsorted binned * contributions size contribSize * photonOut Opened primary output file for sorted photons * contribOut Opened secondary output file for sorted binned * contributions ordered by corresponding photons * contribSize Size of binned contributions per photon, in bytes * numBlk Number of blocks to divide into / merge from * numProc Number of parallel processes for in-core sort * photonBuf Preallocated in-core sort buffer for maxRec photons * contribBuf Preallocated in-core sort buffer for maxRec contribs * bufLen In-core sort buffer length, in number of photons/contribs * keyData Aggregate data for Morton code generation as sort keys * blkLo Start of current block in input files, in number of photons * blkHi End of current block in input files, in number of photons */ { const unsigned long blkLen = blkHi - blkLo + 1; if (!blkLen || blkHi < blkLo) return 0; if (blkLen > bufLen) { unsigned long blkLo2 = blkLo, blkHi2 = blkLo; const double blkLen2 = (double)blkLen / numBlk; FILE *pBlkFile [numBlk], *cBlkFile [numBlk]; Photon photon; char contrib [contribSize]; MortonIdx pri; int b, pid, stat; OOC_SortQueueNode pqueueNodes [numBlk]; OOC_SortQueue pqueue; #ifdef DEBUG_OOC_SORT unsigned long pqCnt = 0; #endif /* ====================================================== * Block too large for in-core sort -> divide into numBlk * subblocks and recurse * ====================================================== */ #ifdef DEBUG_OOC_SORT fprintf(stderr, "*DBG* contribPhotonSort: splitting block [%lu - %lu]\n", blkLo, blkHi ); #endif for (b = 0; b < numBlk; b++) { /* Open temporary file as output for subblock */ if (!(pBlkFile [b] = tmpfile()) || !(cBlkFile [b] = tmpfile())) { perror("contribPhotonSort: failed to open temp block file"); return -1; } /* Set new subblock interval [blkLo2, blkHi2] */ blkHi2 = blkLo + (b + 1) * blkLen2 - 1; if (blkHi2 - blkLo2 + 1 <= bufLen) { /* !!! Subblock [blkLo2, blkHi2] will be sorted in-core on * !!! recursion, so fork kid process. But don't fork more * !!! than numProc kids; wait if necessary */ while (procCnt >= numProc && wait(&stat) >= 0) { /* Kid process terminated; check status, update counter */ if (!WIFEXITED(stat) || WEXITSTATUS(stat)) return -1; procCnt--; } /* Now fork kid process */ if (!(pid = fork())) { /* Recurse on subblock with temp output files */ if (contribPhotonSortRecurse(photonIn, contribIn, pBlkFile [b], cBlkFile [b], contribSize, numBlk, numProc, photonBuf, contribBuf, bufLen, keyData, blkLo2, blkHi2 )) exit(-1); /* !!! Kid exits here. Apparently the parent's tmpfile * !!! isn't deleted (which is what we want), though * !!! some sources suggest using _Exit() instead; * !!! is this implementation specific? Anyway, we * !!! leave cleanup to the parent (see below) */ exit(0); } else if (pid < 0) { perror("contribPhotonSort: failed to fork subprocess"); return -1; } #ifdef DEBUG_OOC_FORK fprintf(stderr, "*DBG* contribPhotonSort: forking kid %d for block %d\n", procCnt, b ); #endif /* Parent continues here */ procCnt++; } else { /* Subblock will NOT be sorted in-core on recursion, so DON'T fork */ if (contribPhotonSortRecurse(photonIn, contribIn, pBlkFile [b], cBlkFile [b], contribSize, numBlk, numProc, photonBuf, contribBuf, bufLen, keyData, blkLo2, blkHi2 )) return -1; } /* Update offset for next subblock */ blkLo2 = blkHi2 + 1; } /* !!! Wait for any forked processes to terminate */ while (procCnt && wait(&stat) >= 0) { /* Check status, bail out on error or update process counter */ if (!WIFEXITED(stat) || WEXITSTATUS(stat)) return -1; procCnt--; } /* ============================================================= * Subblocks are now sorted; prepare priority queue * ============================================================= */ #ifdef DEBUG_OOC_SORT fprintf(stderr, "*DBG* contribPhotonSort: merging block [%lu - %lu]\n", blkLo, blkHi ); #endif /* Init priority queue */ pqueue.head = pqueueNodes; pqueue.tail = 0; pqueue.len = numBlk; /* Fill priority queue by peeking and enqueueing first photon from each subblock */ for (b = 0; b < numBlk; b++) { #ifdef DEBUG_OOC_SORT fseek(pBlkFile [b], 0, SEEK_END); fseek(cBlkFile [b], 0, SEEK_END); if (ftell(pBlkFile [b]) % PHOTONSIZE || ftell(cBlkFile [b]) % contribSize ) { fprintf(stderr, "contribPhotonSort: truncated record " "in temp block file %d\n", b ); return -1; } fprintf(stderr, "contribPhotonSort: temp block file %d " "contains %lu photons, %lu contribs\n", b, ftell(pBlkFile [b]) / PHOTONSIZE, ftell(cBlkFile [b]) / contribSize ); #endif rewind(pBlkFile [b]); rewind(cBlkFile [b]); if (OOC_SortPeek(pBlkFile [b], PHOTONSIZE, (char*)&photon) || OOC_SortPeek(cBlkFile [b], contribSize, contrib) ) { perror("contribPhotonSort: error reading block file"); return -1; } /* Enqueue photon's Morton code as priority */ pri = PMAPCONTRIB_KEY(&photon, keyData); if (OOC_SortPush(&pqueue, pri, b) < 0) { perror("contribPhotonSort: failed priority queue insertion"); return -1; } } /* ========================================================== * Subblocks now sorted and priority queue filled; * merge subblocks from temporary files * ========================================================== */ do { /* Get head of priority queue, read next photon/contrib in * corresponding subblock, and send to output */ if ((b = OOC_SortPop(&pqueue)) >= 0) { /* Priority queue not empty */ if (OOC_SortRead(pBlkFile [b], PHOTONSIZE, (char*)&photon) || OOC_SortRead(cBlkFile [b], contribSize, contrib) ) { /* Corresponding photon/contrib should be in subblock */ perror("contribPhotonSort: unexpected EOF in block file"); return -1; } if (OOC_SortWrite(photonOut, PHOTONSIZE, (char*)&photon) || OOC_SortWrite(contribOut, contribSize, contrib) ) { perror("contribPhotonSort: error writing output file"); return -1; } #ifdef DEBUG_OOC_SORT pqCnt++; #endif /* Peek next photon/contrib from same subblock and enqueue for next iteration */ if (!OOC_SortPeek(pBlkFile [b], PHOTONSIZE, (char*)&photon) && !OOC_SortPeek(cBlkFile [b], contribSize, contrib) ) { /* Subblock not exhausted */ pri = PMAPCONTRIB_KEY(&photon, keyData); if (OOC_SortPush(&pqueue, pri, b) < 0) { perror("contribPhotonSort: failed priority enqueue"); return -1; } } } } while (b >= 0); #ifdef DEBUG_OOC_SORT fprintf(stderr, "*DBG* contribPhotonSort: dequeued %lu rec\n", pqCnt ); fprintf(stderr, "*DBG* contribPhotonSort: merged file contains " "%lu photons, %lu contribs\n", ftell(photonOut) / PHOTONSIZE, ftell(contribOut) / contribSize ); #endif /* Priority queue now empty -> done; close temporary subblock files, * which deletes them (see note in kid kode above) */ for (b = 0; b < numBlk; b++) { fclose(pBlkFile [b]); fclose(cBlkFile [b]); } /* !!! Commit output file to disk before caller reads it; omitting * !!! this apparently leads to corrupted files (race condition?) * !!! when the caller tries to read them! */ fflush(photonOut); fflush(contribOut); fsync(fileno(photonOut)); fsync(fileno(contribOut)); } else { /* ========================================================== * Block is small enough for in-core sort * ========================================================== */ const unsigned long pBlkSize = blkLen * PHOTONSIZE, cBlkSize = blkLen * contribSize; int pIfd = fileno(photonIn), pOfd = fileno(photonOut), cIfd = fileno(contribIn), cOfd = fileno(contribOut); #ifdef DEBUG_OOC_SORT fprintf(stderr, "*DBG* contribPhotonSort: Proc %d (%d/%d) " "sorting block [%lu - %lu]\n", getpid(), procCnt, numProc - 1, blkLo, blkHi ); #endif /* Atomically seek and read block into in-core sort buffer */ /* !!! Unbuffered I/O via pread() avoids potential race conditions * !!! and buffer corruption which can occur with lseek()/fseek() * !!! followed by read()/fread(). */ if (pread(pIfd, photonBuf, pBlkSize, blkLo*PHOTONSIZE) != pBlkSize || pread(cIfd, contribBuf, cBlkSize, blkLo*contribSize) != cBlkSize ) { perror("contribPhotonSort: error reading block for sort"); return -1; } /* Quicksort block in-core and write to output file */ contribPhotonQuickSort(photonBuf, contribBuf, contribSize, keyData, blkLo, blkHi ); if (write(pOfd, photonBuf, pBlkSize) != pBlkSize || write(cOfd, contribBuf, cBlkSize) != cBlkSize ) { perror("contribPhotonSort: error writing sorted block"); return -1; } fsync(pOfd); fsync(cOfd); #ifdef DEBUG_OOC_SORT fprintf(stderr, "*DBG* contribPhotonSort: proc %d wrote " "%lu photons, %lu contribs\n", getpid(), lseek(pOfd, 0, SEEK_END) / PHOTONSIZE, lseek(cOfd, 0, SEEK_END) / contribSize ); #endif } return 0; } int contribPhotonSort ( const PhotonMap *pmap, FILE *leafFile, FVECT bbOrg, RREAL bbSize, unsigned numBlk, unsigned long blkSize, unsigned numProc ) /* Sort photons and their corresponding binned contributions by * subdividing input files into small blocks, sorting these in-core, and * merging them out-of-core via a priority queue. This code is adapted * from oocsort.{h,c} * * The sort is implemented as a recursive (numBlk)-way sort; the input * files are successively split into numBlk smaller blocks until these * are of size <= blkSize, i.e. small enough for in-core sorting, then * merging the sorted blocks as the stack unwinds. The in-core sort is * parallelised over numProc processes. * * Parameters are as follows: * pmap Photon Map with opened photon heap and precomputed * contributions * leafFile Opened output file for sorted photons (=OOC octree leaves) * numBlk Number of blocks to divide into / merge from * blkSize Max block size and size of in-core sort buffer, in bytes * numProc Number of parallel processes for in-core sort * bbOrg Origin of bounding box containing photons * bbSize Extent of bounding box containing photons */ { Photon *photonBuf = NULL; PreComputedContrib* preCompContrib = (PreComputedContrib*)( pmap ? pmap -> preCompContrib : NULL ); char *contribBuf = NULL; - unsigned long bufLen, contribSize; + unsigned long bufLen; int stat; OOC_KeyData keyData; if (numBlk < 1) numBlk = 1; if (!pmap || !pmap -> preCompContrib || !pmap -> numPhotons || !pmap -> heap || !pmap -> contribHeap || !preCompContrib -> nCompressedBins || !preCompContrib -> waveletFile || !numBlk || !blkSize ) error(INTERNAL, "contribPhotonSort: nothing to sort"); /* Figure out maximum buffer length (number of photons/contribs we can fit into the maximum block size) and allocate in-core sort buffers accordingly. Quit if block size too small */ - contribSize = preCompContrib -> nCompressedBins * sizeof(mRGBE); - - if (!(bufLen = blkSize / max(PHOTONSIZE, contribSize))) + bufLen = blkSize / max(PHOTONSIZE, preCompContrib -> contribSize); + if (!bufLen) error(INTERNAL, "contribPhotonSort: block size too small"); if (!(photonBuf = calloc(bufLen, PHOTONSIZE)) || - !(contribBuf = calloc(bufLen, contribSize)) + !(contribBuf = calloc(bufLen, preCompContrib -> contribSize)) ) error(SYSTEM, "contribPhotonSort: cannot allocate in-core sort buffer" ); /* Set key callback and data; reduces params on stack in recursion */ keyData.key = OOC_PhotonKey; keyData.mortonScale = MORTON3D_MAX / bbSize; VCOPY(keyData.bbOrg, bbOrg); stat = contribPhotonSortRecurse( pmap -> heap, pmap -> contribHeap, - leafFile, preCompContrib -> waveletFile, contribSize, - numBlk, numProc, photonBuf, contribBuf, bufLen, &keyData, + leafFile, preCompContrib -> waveletFile, + preCompContrib -> contribSize, numBlk, numProc, + photonBuf, contribBuf, bufLen, &keyData, 0, pmap -> numPhotons - 1 ); /* Cleanup */ free(photonBuf); free(contribBuf); return stat; } #endif /* NIX / PMAP_OOC */