Homec4science
Diffusion RADIANCE Photon Map (fixSpecular) Browse

RadiancePmap/pmcontribsort.c
b1706b2d7c9cfixSpecular

Last Change
Raw File

pmcontribsort.c
View Options

#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)
/* No Windoze 'cos no fork() */
#include "pmapcontrib.h"
#ifdef PMAP_CONTRIB
#include "oocsort.h"
#include "pmutil.h"
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <fcntl.h>
#include <sys/wait.h>
#define PHOTONSIZE sizeof(Photon)
#define PMAPCONTRIB_KEY(p, kdat) Key2Morton3D( \
kdat -> key(p), kdat -> bbOrg, kdat -> mortonScale \
)
#define PMAPCONTRIB_KEYSWAP(i,j) (tKey = (i), (i) = (j), (j) = tKey)
extern RREAL *OOC_PhotonKey (const void *p);
/* Active subprocess counter updated by parent process in
* contribPhotonSort(); must persist when recursing into
* contribPhotonSortRecurse(), hence global */
static unsigned procCnt = 0;
static void contribPhotonSwap (Photon *photons, char *contribs,
unsigned long i, unsigned long j, unsigned contribSize
)
{
/* TODO: swap photon indices instead of photons themselves? */
char tmp [max(PHOTONSIZE, contribSize)],
*ci = contribs + i * contribSize,
*cj = contribs + j * contribSize;
memcpy(tmp, photons + i, PHOTONSIZE);
memcpy(photons + i, photons + j, PHOTONSIZE);
memcpy(photons + j, tmp, PHOTONSIZE);
memcpy(tmp, ci, contribSize);
memcpy(ci, cj, contribSize);
memcpy(cj, tmp, contribSize);
}
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, */
{
unsigned long l, r, m;
MortonIdx lKey, rKey, mKey, tKey;
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(photons, contribs, l, m, contribSize);
PMAPCONTRIB_KEYSWAP(lKey, mKey);
}
if (rKey < lKey) {
contribPhotonSwap(photons, contribs, l, r, contribSize);
PMAPCONTRIB_KEYSWAP(lKey, rKey);
}
if (mKey < rKey) {
contribPhotonSwap(photons, contribs, m, r, contribSize);
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
(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) >= rKey)
r--;
/* Photons out of order; swap them and their contributions */
if (l < r)
contribPhotonSwap(photons, contribs, l, r, contribSize);
};
/* Now l = r = pivot's final position; swap these photons */
contribPhotonSwap(photons, contribs, l, right, contribSize);
/* 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
);
}
}
#ifdef PMAP_CONTRIB_DBG
static unsigned nCoeffs = 0, nCompressedCoeffs = 0;
#endif
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 num photons/contribs
* keyData Aggregate data for Morton code generation as sort keys
* blkLo Start of current block in input files, in num photons
* blkHi End of current block in input files, in num 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 PMAP_CONTRIB_DBG
unsigned long pqCnt = 0;
#endif
/* ======================================================
* Block too large for in-core sort -> divide into numBlk
* subblocks and recurse
* ====================================================== */
#ifdef PMAP_CONTRIB_DBG
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: cannot 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 proc terminated; check status, update counter */
if (!WIFEXITED(stat) || WEXITSTATUS(stat))
return -1;
procCnt--;
}
/* Now fork kid process */
if (!(pid = fork())) {
#if 0
fputs("waiting 10 sec to attach debugga", stderr);
sleep(10);
#endif
/* 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 PMAP_CONTRIB_DBG
fprintf(stderr, "*DBG* contribPhotonSort: forking kid %d "
"(PID %d) for block %d\n",
procCnt, pid, 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 count */
if (!WIFEXITED(stat) || WEXITSTATUS(stat))
return -1;
procCnt--;
}
/* ==========================================================
* Subblocks are now sorted; prepare priority queue
* ========================================================== */
#ifdef PMAP_CONTRIB_DBG
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 PMAP_CONTRIB_DBG
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 insert");
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 PMAP_CONTRIB_DBG
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 PMAP_CONTRIB_DBG
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 PMAP_CONTRIB_DBG
unsigned long contribIdx;
unsigned i;
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;
}
#ifdef PMAP_CONTRIB_DBG
/* Check compressed contribs for bullshit values BEFORE sort */
for (contribIdx = 0; contribIdx < blkLen; contribIdx++)
for (i = 0; i < nCompressedCoeffs; i++)
if (((mRGBE*)(contribBuf + contribIdx * contribSize +
(WAVELET_PADD4_NUMAPPROX + 1 +
(nCompressedCoeffs > 1)
) * sizeof(COLR)
)
) [i].dat >= nCoeffs
)
error(CONSISTENCY, "bullshit contribs BEFORE sort");
#endif
/* Quicksort block in-core and write to output file */
contribPhotonQuickSort(photonBuf, contribBuf, contribSize,
keyData, 0, blkLen - 1
);
#ifdef PMAP_CONTRIB_DBG
/* Check contribs for bullshit values AFTER sort */
for (contribIdx = 0; contribIdx < blkLen; contribIdx++)
for (i = 0; i < nCompressedCoeffs; i++)
if (((mRGBE*)(contribBuf + contribIdx * contribSize +
(WAVELET_PADD4_NUMAPPROX + 1 +
(nCompressedCoeffs > 1)
) * sizeof(COLR)
)
) [i].dat >= nCoeffs
)
error(CONSISTENCY, "bullshit contribs AFTER sort");
#endif
if (write(pOfd, photonBuf, pBlkSize) != pBlkSize ||
write(cOfd, contribBuf, cBlkSize) != cBlkSize
) {
perror("contribPhotonSort: error writing sorted block");
return -1;
}
fsync(pOfd);
fsync(cOfd);
#ifdef PMAP_CONTRIB_DBG
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;
int stat;
OOC_KeyData keyData;
if (numBlk < 1)
numBlk = 1;
if (!pmap || !preCompContrib || !pmap -> numPhotons ||
!numBlk || !blkSize
)
error(INTERNAL, "contribPhotonSort: nothing to sort");
if (preCompContrib -> nBins <= 1)
/* Binning disabled, fall back to photon sorting only */
return OOC_Sort(pmap -> heap, leafFile,
PMAP_OOC_NUMBLK, PMAP_OOC_BLKSIZE, numProc, sizeof(Photon),
bbOrg, bbSize, OOC_PhotonKey
);
#ifdef PMAP_CONTRIB_DBG
nCoeffs = preCompContrib -> nCoeffs;
nCompressedCoeffs = preCompContrib -> nCompressedCoeffs;
/* Check photon & contrib heap lengths before actual sort */
if (fseek(pmap -> heap, 0, SEEK_END) < 0 ||
ftell(pmap -> heap) / PHOTONSIZE != pmap -> numPhotons
)
error(INTERNAL,
"contribPhotonSort: photon heap length mismatch"
);
if (fseek(pmap -> contribHeap, 0, SEEK_END) < 0 ||
ftell(pmap -> contribHeap) / preCompContrib -> contribSize !=
pmap -> numPhotons
)
error(INTERNAL,
"contribPhotonSort: contribution heap length mismatch"
);
#endif
/* 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 */
bufLen = blkSize / max(PHOTONSIZE, preCompContrib -> contribSize);
if (!bufLen)
error(INTERNAL, "contribPhotonSort: block size too small");
if (!(photonBuf = calloc(bufLen, PHOTONSIZE)) ||
!(contribBuf = calloc(bufLen, preCompContrib -> contribSize))
)
error(SYSTEM,
"contribPhotonSort: cannot allocate in-core sort buffer"
);
/* Set key callback and data; reduces params on stack */
keyData.key = OOC_PhotonKey;
keyData.mortonScale = MORTON3D_MAX / bbSize;
VCOPY(keyData.bbOrg, bbOrg);
stat = contribPhotonSortRecurse(
pmap -> heap, pmap -> contribHeap,
leafFile, preCompContrib -> waveletFile,
preCompContrib -> contribSize, numBlk, numProc,
photonBuf, contribBuf, bufLen, &keyData,
0, pmap -> numPhotons - 1
);
/* Cleanup */
free(photonBuf);
free(contribBuf);
return stat;
}
#endif /* PMAP_CONTRIB */
#endif /* NIX */
c4science · Help