Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F86179669
pmapcontrib.c
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Oct 4, 18:26
Size
61 KB
Mime Type
text/x-c
Expires
Sun, Oct 6, 18:26 (2 d)
Engine
blob
Format
Raw Data
Handle
21365536
Attached To
R10977 RADIANCE Photon Map
pmapcontrib.c
View Options
#ifndef lin
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"
#ifdef PMAP_CONTRIB
#include "pmapdiag.h"
#include "pmaprand.h"
#include "pmapmat.h"
#include "pmapsrc.h"
#include "pmutil.h"
#include "otspecial.h"
#include "otypes.h"
#include "lookup.h"
#if NIX
#include <sys/mman.h>
#include <sys/wait.h>
#endif
/* The following are convenient placeholders interfacing to mkpmap
and rcontrib. They can be externally set via initPmapContribTab()
and referenced within the contrib pmap modules. These variables
can then be safely ignored by rtrace/rpict/rvu, without annoying
linking errors. */
/* Global pointer to rcontrib's contribution binning LUT */
LUTAB
*
pmapContribTab
=
NULL
;
/* Contribution/coefficient mode flag */
int
*
pmapContribMode
;
#ifdef PMAP_CONTRIB_DBG
/* Sum/num of mRGBE encoding errors for wavelet coeffs */
static
WAVELET_COEFF
mrgbeDiffs
=
0
;
static
unsigned
long
mrgbeDiffCnt
=
0
;
#endif
static
int
ray2bin
(
const
RAY
*
ray
,
unsigned
scDim
)
/* Map ray dir (pointing away from origin) to its 1D bin in an
(scDim x scDim) Shirley-Chiu square.
Returns -1 if mapped bin is invalid (e.g. behind plane) */
{
static
int
scRHS
,
varInit
=
0
;
static
FVECT
scNorm
,
scUp
;
unsigned
scBin
[
2
];
FVECT
diskPlane
;
RREAL
dz
,
diskx
,
disky
,
rad
,
diskd2
,
scCoord
[
2
];
if
(
!
varInit
)
{
/* Lazy init Shirley-Chiu mapping orientation from func variables */
/* XXX: Note inversion of normal and up vectors to match rcontrib's
convention, since the input ray points away from its origin */
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 */
/* normalize(diskPlane); NEED TO NORMALISE? */
diskPlane
[
0
]
*=
scRHS
;
diskx
=
DOT
(
ray
->
rdir
,
diskPlane
);
disky
=
DOT
(
ray
->
rdir
,
scUp
);
/* Apply Shirley-Chiu mapping of (diskx, disky) to square */
disk2square
(
scCoord
,
diskx
,
disky
);
/* Map Shirley-Chiu square coords to 1D bin */
scBin
[
0
]
=
scCoord
[
0
]
*
scDim
;
scBin
[
1
]
=
scCoord
[
1
]
*
scDim
;
return
PMAP_CONTRIB_XY2LIN
(
scDim
,
scBin
[
0
],
scBin
[
1
]);
}
else
return
-
1
;
}
/* ------------------ CONTRIBSRC STUFF --------------------- */
#ifndef PMAP_CONTRIB_TEST
MODCONT
*
addContribModifier
(
LUTAB
*
contribTab
,
unsigned
*
numContribs
,
char
*
mod
,
char
*
binParm
,
int
binCnt
)
/* Add light source modifier mod to contribution lookup table contribTab,
and update numContribs. Return initialised contribution data for this
modifier. This code adapted from rcontrib.c:addmodifier(). */
{
LUENT
*
lutEntry
=
lu_find
(
contribTab
,
mod
);
MODCONT
*
contrib
;
EPNODE
*
eBinVal
;
unsigned
numCoeffs
;
/* Warn if potential coefficient index overflow in mRGBE encoding.
This requires getting the number of wavelet coefficients generated
by the transform a priori. */
if
(
!
(
numCoeffs
=
padWaveletXform2
(
NULL
,
NULL
,
sqrt
(
binCnt
),
NULL
)))
{
sprintf
(
errmsg
,
"can't determine number of wavelet coefficients "
"for modifier %s"
,
mod
);
error
(
INTERNAL
,
errmsg
);
}
if
(
numCoeffs
*
numCoeffs
>
PMAP_CONTRIB_MAXCOEFFS
)
{
sprintf
(
errmsg
,
"mRGBE data field may overflow for modifier %s; reduce -bn "
"and/or compression if contribution precomputation fails"
,
mod
);
error
(
WARNING
,
errmsg
);
}
if
(
lutEntry
->
data
)
{
/* Reject duplicate modifiers */
sprintf
(
errmsg
,
"duplicate light source modifier %s"
,
mod
);
error
(
USER
,
errmsg
);
}
if
(
*
numContribs
>=
MAXMODLIST
)
{
sprintf
(
errmsg
,
"too many modifiers (max. %d)"
,
MAXMODLIST
);
error
(
INTERNAL
,
errmsg
);
}
lutEntry
->
key
=
mod
;
if
(
binCnt
<=
0
)
{
sprintf
(
errmsg
,
"undefined/invalid bin count for modifier %s"
,
mod
);
error
(
USER
,
errmsg
);
}
/* Allocate and init contributions */
contrib
=
(
MODCONT
*
)
malloc
(
sizeof
(
MODCONT
)
+
sizeof
(
DCOLOR
)
*
(
binCnt
-
1
)
);
if
(
!
contrib
)
error
(
SYSTEM
,
"out of memory in addContribModifier()"
);
contrib
->
outspec
=
NULL
;
contrib
->
modname
=
mod
;
contrib
->
params
=
binParm
;
contrib
->
nbins
=
binCnt
;
contrib
->
binv
=
eBinVal
;
contrib
->
bin0
=
0
;
memset
(
contrib
->
cbin
,
0
,
sizeof
(
DCOLOR
)
*
binCnt
);
lutEntry
->
data
=
lutEntry
->
data
=
(
char
*
)
contrib
;
++
(
*
numContribs
);
return
(
contrib
);
}
void
addContribModfile
(
LUTAB
*
contribTab
,
unsigned
*
numContribs
,
char
*
modFile
,
char
*
binParm
,
int
binCnt
)
/* Add light source modifiers from file modFile to contribution lookup table
* contribTab, and update numContribs.
* NOTE: This code is adapted from rcontrib.c */
{
char
*
mods
[
MAXMODLIST
];
int
i
;
/* Find file and store strings */
i
=
wordfile
(
mods
,
MAXMODLIST
,
getpath
(
modFile
,
getrlibpath
(),
R_OK
));
if
(
i
<
0
)
{
sprintf
(
errmsg
,
"can't open modifier file %s"
,
modFile
);
error
(
SYSTEM
,
errmsg
);
}
if
(
*
numContribs
+
i
>=
MAXMODLIST
-
1
)
{
sprintf
(
errmsg
,
"too many modifiers (max. %d) in file %s"
,
MAXMODLIST
-
1
,
modFile
);
error
(
INTERNAL
,
errmsg
);
}
for
(
i
=
0
;
mods
[
i
];
i
++
)
/* Add each modifier */
addContribModifier
(
contribTab
,
numContribs
,
mods
[
i
],
binParm
,
binCnt
);
}
static
int
contribSourceBin
(
LUTAB
*
contribs
,
const
RAY
*
ray
)
/* Map contribution source ray to its bin for light source ray -> rsrc,
using Shirley-Chiu disk-to-square mapping.
Return invalid bin -1 if mapping failed. */
{
const
SRCREC
*
src
;
const
OBJREC
*
srcMod
;
const
MODCONT
*
srcCont
;
RAY
srcRay
;
int
bin
,
i
;
/* Check we have a valid ray and contribution LUT */
if
(
!
ray
||
!
contribs
)
return
-
1
;
src
=
&
source
[
ray
->
rsrc
];
srcMod
=
findmaterial
(
src
->
so
);
srcCont
=
(
MODCONT
*
)
lu_find
(
contribs
,
srcMod
->
oname
)
->
data
;
if
(
!
srcCont
)
/* Not interested in this source (modifier not in contrib LUT) */
return
-
1
;
/* Set up shadow ray pointing to source for disk2square mapping */
rayorigin
(
&
srcRay
,
SHADOW
,
NULL
,
NULL
);
srcRay
.
rsrc
=
ray
->
rsrc
;
VCOPY
(
srcRay
.
rorg
,
ray
->
rop
);
for
(
i
=
0
;
i
<
3
;
i
++
)
srcRay
.
rdir
[
i
]
=
-
ray
->
rdir
[
i
];
if
(
!
(
src
->
sflags
&
SDISTANT
?
sourcehit
(
&
srcRay
)
:
(
*
ofun
[
srcMod
->
otype
].
funp
)(
srcMod
,
&
srcRay
)
))
/* (Redundant?) sanity check for valid source ray? */
return
-
1
;
#if 0
worldfunc(RCCONTEXT, &srcRay);
#endif
set_eparams
((
char
*
)
srcCont
->
params
);
if
((
bin
=
ray2bin
(
&
srcRay
,
sqrt
(
srcCont
->
nbins
)))
<
0
)
error
(
WARNING
,
"Ignoring invalid bin in contribSourceBin()"
);
return
bin
;
}
PhotonContribSourceIdx
newPhotonContribSource
(
PhotonMap
*
pmap
,
const
RAY
*
contribSrcRay
,
FILE
*
contribSrcHeap
)
/* Add contribution source ray for emitted contribution photon and save
* light source index and binned direction. The current contribution
* source is stored in pmap -> lastContribSrc. If the previous contrib
* source spawned photons (i.e. has srcIdx >= 0), it's appended to
* contribSrcHeap. If contribSrcRay == NULL, the current contribution
* source is still flushed, but no new source is set. Returns updated
* contribution source counter pmap -> numContribSrc. */
{
if
(
!
pmap
||
!
contribSrcHeap
)
return
0
;
/* Check if last contribution source has spawned photons (srcIdx >= 0,
* see newPhoton()), in which case we save it to the heap file before
* clobbering it. (Note this is short-term storage, so we doan' need
* da portable I/O stuff here). */
if
(
pmap
->
lastContribSrc
.
srcIdx
>=
0
)
{
if
(
!
fwrite
(
&
pmap
->
lastContribSrc
,
sizeof
(
PhotonContribSource
),
1
,
contribSrcHeap
))
error
(
SYSTEM
,
"failed writing photon contrib source in "
"newPhotonContribSource()"
);
pmap
->
numContribSrc
++
;
if
(
!
pmap
->
numContribSrc
||
pmap
->
numContribSrc
>
PMAP_MAXCONTRIBSRC
)
error
(
INTERNAL
,
"contribution source overflow"
);
}
if
(
contribSrcRay
)
{
/* Mark this contribution source unused with a negative source
index until its path spawns a photon (see newPhoton()) */
pmap
->
lastContribSrc
.
srcIdx
=
-
1
;
/* Map ray to bin in anticipation that this contrib source will be
used, since the ray will be lost once a photon is spawned */
pmap
->
lastContribSrc
.
srcBin
=
contribSourceBin
(
pmapContribTab
,
contribSrcRay
);
if
(
pmap
->
lastContribSrc
.
srcBin
<
0
)
{
/* Warn if invalid bin, but trace photon nonetheless. It will
count as emitted to prevent bias, but will not be stored in
newPhoton(), as it contributes zero flux */
sprintf
(
errmsg
,
"invalid bin for light source %s, "
"contribution photons will be discarded"
,
source
[
contribSrcRay
->
rsrc
].
so
->
oname
);
error
(
WARNING
,
errmsg
);
}
}
return
pmap
->
numContribSrc
;
}
PhotonContribSourceIdx
buildContribSources
(
PhotonMap
*
pmap
,
FILE
**
contribSrcHeap
,
char
**
contribSrcHeapFname
,
PhotonContribSourceIdx
*
contribSrcOfs
,
unsigned
numHeaps
)
/* Consolidate per-subprocess contribution source heaps into array
* pmap -> contribSrc. Returns offset for contribution source index
* linearisation in pmap -> numContribSrc. The heap files in
* contribSrcHeap are closed on return. */
{
PhotonContribSourceIdx
heapLen
;
unsigned
heap
;
if
(
!
pmap
||
!
contribSrcHeap
||
!
contribSrcOfs
||
!
numHeaps
)
return
0
;
pmap
->
numContribSrc
=
0
;
for
(
heap
=
0
;
heap
<
numHeaps
;
heap
++
)
{
contribSrcOfs
[
heap
]
=
pmap
->
numContribSrc
;
if
(
fseek
(
contribSrcHeap
[
heap
],
0
,
SEEK_END
)
<
0
)
error
(
SYSTEM
,
"failed photon contrib source seek "
"in buildContribSources()"
);
pmap
->
numContribSrc
+=
heapLen
=
ftell
(
contribSrcHeap
[
heap
])
/
sizeof
(
PhotonContribSource
);
if
(
!
(
pmap
->
contribSrc
=
realloc
(
pmap
->
contribSrc
,
pmap
->
numContribSrc
*
sizeof
(
PhotonContribSource
)
)))
error
(
SYSTEM
,
"failed photon contrib source alloc "
"in buildContribSources()"
);
rewind
(
contribSrcHeap
[
heap
]);
if
(
fread
(
pmap
->
contribSrc
+
contribSrcOfs
[
heap
],
sizeof
(
PhotonContribSource
),
heapLen
,
contribSrcHeap
[
heap
]
)
!=
heapLen
)
error
(
SYSTEM
,
"failed reading photon contrib source "
"in buildContribSources()"
);
fclose
(
contribSrcHeap
[
heap
]);
unlink
(
contribSrcHeapFname
[
heap
]);
}
return
pmap
->
numContribSrc
;
}
/* ----------------- CONTRIB PMAP ENCODING STUFF -------------------- */
static
void
coeffSwap
(
PreComputedContribCoeff
*
c1
,
PreComputedContribCoeff
*
c2
)
{
PreComputedContribCoeff
tCoeff
;
tCoeff
.
coeff
=
c1
->
coeff
;
tCoeff
.
idx
=
c1
->
idx
;
c1
->
coeff
=
c2
->
coeff
;
c1
->
idx
=
c2
->
idx
;
c2
->
coeff
=
tCoeff
.
coeff
;
c2
->
idx
=
tCoeff
.
idx
;
}
static
int
coeffPartition
(
PreComputedContribCoeff
*
coeffs
,
unsigned
medianPos
,
unsigned
left
,
unsigned
right
)
/* REVERSE partition coefficients by magnitude, such that
coeffs [left..medianPos] >= coeffs [medianPos+1..right].
Returns median's position. */
{
unsigned
long
l
,
r
,
m
;
WAVELET_COEFF
lVal
,
rVal
,
mVal
,
tVal
;
#define COEFFVAL(c) (DOT((c) -> coeff, (c) -> coeff))
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 */
lVal
=
COEFFVAL
(
coeffs
+
l
);
rVal
=
COEFFVAL
(
coeffs
+
r
);
mVal
=
COEFFVAL
(
coeffs
+
m
);
if
(
mVal
>
lVal
)
{
coeffSwap
(
coeffs
+
m
,
coeffs
+
l
);
tVal
=
mVal
;
mVal
=
lVal
;
lVal
=
tVal
;
}
if
(
rVal
>
lVal
)
{
coeffSwap
(
coeffs
+
r
,
coeffs
+
l
);
tVal
=
rVal
;
rVal
=
lVal
;
lVal
=
tVal
;
}
if
(
mVal
>
rVal
)
{
coeffSwap
(
coeffs
+
m
,
coeffs
+
r
);
tVal
=
mVal
;
mVal
=
rVal
;
rVal
=
tVal
;
}
/* Pivot with key rVal is now in coeffs [right] */
/* l & r converge, swapping coefficients out of (reversed) order
with respect to pivot. The convergence point l = r is the
pivot's final position */
while
(
l
<
r
)
{
while
(
l
<
r
&&
COEFFVAL
(
coeffs
+
l
)
>
rVal
)
l
++
;
while
(
r
>
l
&&
COEFFVAL
(
coeffs
+
r
)
<=
rVal
)
r
--
;
/* Coeffs out of order, must swap */
if
(
l
<
r
)
coeffSwap
(
coeffs
+
l
,
coeffs
+
r
);
};
/* Now l == r is pivot's final position; swap these coeffs */
coeffSwap
(
coeffs
+
l
,
coeffs
+
right
);
/* Recurse in partition containing median */
if
(
l
>
medianPos
)
return
coeffPartition
(
coeffs
,
medianPos
,
left
,
l
-
1
);
else
if
(
l
<
medianPos
)
return
coeffPartition
(
coeffs
,
medianPos
,
l
+
1
,
right
);
else
/* l == medianPos, partitioning done */
return
l
;
}
else
return
left
;
}
static
int
coeffIdxCompare
(
const
void
*
c1
,
const
void
*
c2
)
/* Comparison function to sort thresholded coefficients by index */
/* TODO: scale coeffs by 2^-l to adapt to resolution/frequency l?
(this drops higher-frequency coeffs first). */
{
const
PreComputedContribCoeff
*
tcoeff1
=
c1
,
*
tcoeff2
=
c2
;
const
unsigned
v1
=
tcoeff1
->
idx
,
v2
=
tcoeff2
->
idx
;
if
(
v1
<
v2
)
return
-
1
;
else
if
(
v1
>
v2
)
return
1
;
else
return
0
;
}
static
int
thresholdContribs
(
PreComputedContrib
*
preCompContrib
)
/* Threshold wavelet detail coefficients in preCompContrib ->
waveletMatrix [approxDim..coeffDim-1] [approxDim..coeffDim-1] (where
approxDim = WAVELET_PADD2_APPROXDIM) by keeping the (preCompContrib ->
nCompressedCoeffs) largest of these and returning them in
preCompContrib -> compressedCoeffs along with their original matrix
indices.
NOTE: The wavelet approximation coefficients
preCompContrib -> waveletMatrix [0..approxDim-1] [0..approxDim-1]
are excluded from thresholding to minimise compression artefacts.
Returns 0 on success, else -1. */
{
unsigned
i
,
j
,
coeffDim
,
coeffIdx
,
nNzDetailCoeffs
,
nCompressedCoeffs
,
numThresh
;
WaveletMatrix2
waveletMatrix
;
PreComputedContribCoeff
*
threshCoeffs
,
*
threshCoeffPtr
;
WAVELET_COEFF
*
coeffPtr
;
if
(
!
preCompContrib
||
!
(
coeffDim
=
preCompContrib
->
coeffDim
)
||
!
(
threshCoeffs
=
preCompContrib
->
compressedCoeffs
)
||
!
(
waveletMatrix
=
preCompContrib
->
waveletMatrix
)
)
/* Should be initialised by caller! */
return
-
1
;
/* Set up coefficient buffer for compression (thresholding), skipping
* the approximation coefficients in the upper left of waveletMatrix,
* which are not thresholded. Also skip zero coefficients (resulting
* from padding), since these are already implicitly thresholded. The
* original 2D matrix indices are linearised to 1D and saved with each
* coefficient to restore the original sparse coefficient matrix. */
for
(
i
=
0
,
threshCoeffPtr
=
threshCoeffs
;
i
<
coeffDim
;
i
++
)
{
/* Calc row pointer (=mult) in outer loop, inc in inner */
coeffIdx
=
PMAP_CONTRIB_XY2LIN
(
coeffDim
,
i
,
0
);
for
(
j
=
0
;
j
<
coeffDim
;
j
++
,
coeffIdx
++
)
{
if
((
i
>=
WAVELET_PADD2_APPROXDIM
||
j
>=
WAVELET_PADD2_APPROXDIM
)
&&
!
coeffIsZero
(
waveletMatrix
[
i
]
[
j
])
)
{
/* Nonzero detail coefficient; set up pointer to coeff
(sorts faster than 3 doubles) and save original
(linearised) matrix index */
threshCoeffPtr
->
idx
=
coeffIdx
;
threshCoeffPtr
++
->
coeff
=
(
WAVELET_COEFF
*
)
&
(
waveletMatrix
[
i
]
[
j
]
);
}
}
}
/* Num of nonzero detail coeffs in buffer, number actually expected */
numThresh
=
threshCoeffPtr
-
threshCoeffs
;
nNzDetailCoeffs
=
WAVELET_PADD2_NUMDETAIL
(
preCompContrib
->
nNonZeroCoeffs
);
nCompressedCoeffs
=
preCompContrib
->
nCompressedCoeffs
;
/* If fewer nonzero detail coeff are in the threshold buffer than
* anticipated, the loop below fills the remainder of the threshold
* buffer with duplicates of a coefficient in the lower right of the
* matrix, which is padding and guaranteed to be zero. This condition
* can occur if the wavelet transform actually generated genuine zero
* detail coefficients. Infact it's quite common if the wavelet
* transformed contributions are quite sparse. */
i
=
j
=
coeffDim
-
1
;
coeffIdx
=
PMAP_CONTRIB_XY2LIN
(
coeffDim
,
i
,
j
);
for
(
coeffPtr
=
(
WAVELET_COEFF
*
)(
waveletMatrix
[
i
][
j
]);
numThresh
<
nNzDetailCoeffs
;
numThresh
++
)
{
threshCoeffPtr
->
idx
=
coeffIdx
;
threshCoeffPtr
++
->
coeff
=
coeffPtr
;
}
/* Partition coeffs in reverse order, such that all coeffs in
threshCoeffs [0..nCompressedCoeffs-1] are larger than those in
threshCoeffs [nCompressedCoeffs..nNzDetailCoeffs-1] */
coeffPartition
(
threshCoeffs
,
nCompressedCoeffs
-
1
,
0
,
nNzDetailCoeffs
-
1
);
#ifdef PMAP_CONTRIB_DBG
/* Check coefficient partitioning */
threshCoeffPtr
=
threshCoeffs
+
nCompressedCoeffs
-
1
;
for
(
i
=
0
;
i
<
nCompressedCoeffs
;
i
++
)
if
(
COEFFVAL
(
threshCoeffs
+
i
)
<
COEFFVAL
(
threshCoeffPtr
))
error
(
CONSISTENCY
,
"invalid wavelet coefficient partitioning "
"in thresholdContribs()"
);
for
(;
i
<
nNzDetailCoeffs
;
i
++
)
if
(
COEFFVAL
(
threshCoeffs
+
i
)
>
COEFFVAL
(
threshCoeffPtr
))
error
(
CONSISTENCY
,
"invalid wavelet coefficient partitioning "
"in thresholdContribs()"
);
#endif
/* Sort partition containing nCompressedCoeffs coefficients by index
* (ignoring the remaining coefficients since these are now dropped
* due to tresholding) */
qsort
(
threshCoeffs
,
nCompressedCoeffs
,
sizeof
(
PreComputedContribCoeff
),
coeffIdxCompare
);
return
0
;
}
static
int
encodeContribs
(
PreComputedContrib
*
preCompContrib
,
float
compressRatio
)
/* Apply wavelet transform to input matrix preCompContrib ->
waveletMatrix and compress according to compressRatio, storing
thresholded and mRGBE-encoded coefficients in preCompContrib ->
mrgbeCoeffs. Note that the approximation coefficients in the upper
left of the matrix are not encoded, and returned in
preCompContrib -> waveletMatrix
[0..WAVELET_PADD2_APPROXDIM-1] [0..WAVELET_PADD2_APPROXDIM-1].
Returns 0 on success. */
{
unsigned
i
,
j
,
k
,
scDim
,
lastCoeffIdx
;
WaveletMatrix2
waveletMatrix
,
tWaveletMatrix
;
PreComputedContribCoeff
*
threshCoeffs
;
mRGBERange
*
mrgbeRange
;
mRGBE
*
mrgbeCoeffs
;
WAVELET_COEFF
absCoeff
;
#ifdef PMAP_CONTRIB_DBG
WaveletCoeff3
decCoeff
;
unsigned
decIdx
;
#endif
if
(
!
preCompContrib
||
!
preCompContrib
->
mrgbeCoeffs
||
!
preCompContrib
->
compressedCoeffs
||
!
(
waveletMatrix
=
preCompContrib
->
waveletMatrix
)
||
!
(
tWaveletMatrix
=
preCompContrib
->
tWaveletMatrix
)
||
!
(
scDim
=
preCompContrib
->
scDim
)
)
/* Should be initialised by the caller! */
return
-
1
;
#ifdef PMAP_CONTRIB_DBG
for
(
i
=
0
;
i
<
scDim
;
i
++
)
{
for
(
j
=
0
;
j
<
scDim
;
j
++
)
{
for
(
k
=
0
;
k
<
3
;
k
++
)
{
#if 0
/* Set contributions to bins for debugging */
waveletMatrix [i] [j] [k] = PMAP_CONTRIB_VAL(
PMAP_CONTRIB_XY2LIN(scDim, i, j)
);
#elif 0
/* Replace contribs with "bump" function */
waveletMatrix
[
i
]
[
j
]
[
k
]
=
PMAP_CONTRIB_VAL
(
(
1.
-
fabs
(
i
-
scDim
/
2.
+
0.5
)
/
(
scDim
/
2.
-
0.5
))
*
(
1.
-
fabs
(
j
-
scDim
/
2.
+
0.5
)
/
(
scDim
/
2.
-
0.5
))
);
#elif 0
/* Inject reference contribs from rcontrib classic */
#include "rc-ref.c"
if
(
preCompContrib
->
nBins
!=
PMAP_CONTRIB_REFSIZE
)
{
sprintf
(
errmsg
,
"reference contribs require %d bins"
,
PMAP_CONTRIB_REFSIZE
);
error
(
USER
,
errmsg
);
}
waveletMatrix
[
i
]
[
j
]
[
k
]
=
PMAP_CONTRIB_VAL
(
refContrib
[
PMAP_CONTRIB_XY2LIN
(
scDim
,
i
,
j
)]
[
k
]
);
#endif
}
#if 0
/* Dump contribs prior to encoding for debugging */
printf("% 7.3g\t", colorAvg(waveletMatrix [i] [j]));
}
putchar('\n');
}
putchar('\n');
#else
}
}
#endif
#endif
/* Do 2D wavelet transform on preCompContrib -> waveletMatrix */
if
(
padWaveletXform2
(
waveletMatrix
,
tWaveletMatrix
,
scDim
,
NULL
)
!=
preCompContrib
->
coeffDim
)
error
(
INTERNAL
,
"failed wavelet transform in encodeContribs()"
);
/* Compress wavelet detail coeffs by thresholding */
if
(
thresholdContribs
(
preCompContrib
)
<
0
)
error
(
INTERNAL
,
"failed wavelet compression in encodeContribs()"
);
threshCoeffs
=
preCompContrib
->
compressedCoeffs
;
/* Init per-channel coefficient range for mRGBE encoding */
mrgbeRange
=
&
preCompContrib
->
mrgbeRange
;
setcolor
(
mrgbeRange
->
min
,
FHUGE
,
FHUGE
,
FHUGE
);
setcolor
(
mrgbeRange
->
max
,
0
,
0
,
0
);
/* Update per-channel coefficient range */
for
(
i
=
0
;
i
<
preCompContrib
->
nCompressedCoeffs
;
i
++
)
for
(
k
=
0
;
k
<
3
;
k
++
)
{
absCoeff
=
fabs
(
threshCoeffs
[
i
].
coeff
[
k
]);
if
(
absCoeff
<
mrgbeRange
->
min
[
k
])
mrgbeRange
->
min
[
k
]
=
absCoeff
;
if
(
absCoeff
>
mrgbeRange
->
max
[
k
])
mrgbeRange
->
max
[
k
]
=
absCoeff
;
}
if
(
preCompContrib
->
nCompressedCoeffs
==
1
)
/* Maximum compression with just 1 (!) compressed detail coeff (is
* this even useful?), so mRGBE range is undefined since min & max
* coincide; set minimum to 0, maximum to the single remaining
* coefficient */
setcolor
(
mrgbeRange
->
min
,
0
,
0
,
0
);
/* Init mRGBE coefficient normalisation from range */
if
(
!
mRGBEinit
(
mrgbeRange
,
mrgbeRange
->
min
,
mrgbeRange
->
max
))
return
-
1
;
mrgbeCoeffs
=
preCompContrib
->
mrgbeCoeffs
;
/* Encode wavelet detail coefficients to mRGBE */
for
(
i
=
lastCoeffIdx
=
0
;
i
<
preCompContrib
->
nCompressedCoeffs
;
lastCoeffIdx
=
threshCoeffs
[
i
++
].
idx
)
{
/* HACK: To reduce the likelihood of overflowing the mRGBE data
* field with the coefficient index, it is expressed incrementally
* w.r.t. the previously encoded coefficient's index, instead of
* absolutely. This implies threshCoeffs must be sorted by
* coefficient index to ensure increments are positive, and to
* minimise their magnitude.
* Note that an overflow cannot be predicted beforehand, e.g. by
* mkpmap when parsing the number of bins, as this depends on the
* sparseness of the wavelet coefficients (which in turn depends on
* the frequency distribution of the binned contributions), and the
* fraction of those that are dropped (i.e. the compression ratio).
* */
mrgbeCoeffs
[
i
]
=
mRGBEencode
(
threshCoeffs
[
i
].
coeff
,
mrgbeRange
,
threshCoeffs
[
i
].
idx
-
lastCoeffIdx
);
if
(
!
mrgbeCoeffs
[
i
].
all
)
error
(
INTERNAL
,
"failed mRGBE encoding in encodeContribs()"
);
#ifdef PMAP_CONTRIB_DBG
/* mRGBE encoding sanity check */
decIdx
=
mRGBEdecode
(
mrgbeCoeffs
[
i
],
mrgbeRange
,
decCoeff
);
mrgbeDiffs
+=
sqrt
(
dist2
(
decCoeff
,
threshCoeffs
[
i
].
coeff
));
mrgbeDiffCnt
++
;
if
(
/*decIdx > preCompContrib -> nCompressedCoeffs || */
decIdx
!=
threshCoeffs
[
i
].
idx
-
lastCoeffIdx
||
sqrt
(
dist2
(
decCoeff
,
threshCoeffs
[
i
].
coeff
))
>
0.1
*
colorAvg
(
mrgbeRange
->
max
)
)
{
sprintf
(
errmsg
,
"failed sanity check in encodeContribs()
\n
"
"Encoded: [%.3g %.3g %.3g %d]
\n
Decoded: [%.3g %.3g %.3g %d]"
,
threshCoeffs
[
i
].
coeff
[
0
],
threshCoeffs
[
i
].
coeff
[
1
],
threshCoeffs
[
i
].
coeff
[
2
],
threshCoeffs
[
i
].
idx
-
lastCoeffIdx
,
decCoeff
[
0
],
decCoeff
[
1
],
decCoeff
[
2
],
decIdx
);
error
(
CONSISTENCY
,
errmsg
);
}
for
(
k
=
0
;
k
<
3
;
k
++
)
{
absCoeff
=
fabs
(
threshCoeffs
[
i
].
coeff
[
k
])
>
1e-6
?
decCoeff
[
k
]
/
threshCoeffs
[
i
].
coeff
[
k
]
:
0
;
if
(
absCoeff
<
0
)
error
(
CONSISTENCY
,
"coefficient sign reversal in encodeContribs()"
);
}
#endif
}
return
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
*
modCont
=
(
MODCONT
*
)
lu_find
(
pmapContribTab
,
srcMod
->
oname
)
->
data
;
Photon
*
photon
;
COLOR
flux
;
DCOLOR
*
cbin
;
PhotonSearchQueueNode
*
sqn
;
float
r2
,
norm
;
unsigned
i
,
bin
,
numZeroBins
;
setcolor
(
irrad
,
0
,
0
,
0
);
if
(
!
modCont
)
/* Modifier not in LUT */
return
NULL
;
/* Zero bins for this source modifier */
memset
(
modCont
->
cbin
,
0
,
sizeof
(
DCOLOR
)
*
modCont
->
nbins
);
if
(
!
pmap
->
maxGather
)
/* Zero bandwidth */
return
NULL
;
/* Lookup photons; 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
:
"<null>"
,
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
));
#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
,
numZeroBins
=
modCont
->
nbins
;
i
<
pmap
->
squeue
.
tail
;
i
++
,
sqn
++
)
{
/* Get photon's contribution to density estimate */
photon
=
getNearestPhoton
(
&
pmap
->
squeue
,
sqn
->
idx
);
getPhotonFlux
(
photon
,
flux
);
scalecolor
(
flux
,
norm
);
#ifdef PMAP_EPANECHNIKOV
/* Apply Epanechnikov kernel to photon flux based on distance */
scalecolor
(
flux
,
1
-
sqn
->
dist2
/
r2
);
#endif
addcolor
(
irrad
,
flux
);
cbin
=
&
modCont
->
cbin
[
photonSrcBin
(
pmap
,
photon
)];
/* Bin was empty (zero) before accumulating flux;
decrement corresponding counter */
numZeroBins
-=
coeffIsZero
(
*
cbin
);
addcolor
(
*
cbin
,
flux
);
}
#ifdef PMAP_CONTRIB_WARNZEROBINS
/* Warn if too many empty bins */
if
(
numZeroBins
>
PMAP_CONTRIB_MAXZEROBINS
*
modCont
->
nbins
)
{
sprintf
(
errmsg
,
"empty bin ratio > %.2f at (%.3f, %.3f, %.3f) "
"for modifier %s; contributions may be biased"
,
PMAP_CONTRIB_MAXZEROBINS
,
ray
->
rop
[
0
],
ray
->
rop
[
1
],
ray
->
rop
[
2
],
modCont
->
modname
);
error
(
WARNING
,
errmsg
);
}
#endif
return
modCont
;
}
void
freePreCompContribNode
(
void
*
p
)
/* Clean up precomputed contribution LUT entry */
{
PhotonMap
*
preCompContribPmap
=
(
PhotonMap
*
)
p
;
PreComputedContrib
*
preCompContrib
;
if
(
preCompContribPmap
)
{
preCompContrib
=
(
PreComputedContrib
*
)(
preCompContribPmap
->
preCompContrib
);
if
(
preCompContrib
)
{
/* Free primary and transposed wavelet matrices */
freeWaveletMatrix2
(
preCompContrib
->
waveletMatrix
,
preCompContrib
->
coeffDim
);
freeWaveletMatrix2
(
preCompContrib
->
tWaveletMatrix
,
preCompContrib
->
coeffDim
);
/* Free thresholded coefficients */
free
(
preCompContrib
->
compressedCoeffs
);
/* Free mRGBE encoded coefficients */
free
(
preCompContrib
->
mrgbeCoeffs
);
free
(
preCompContrib
->
waveletFname
);
if
(
preCompContrib
->
cache
)
{
/* Free contribution cache */
OOC_DeleteCache
(
preCompContrib
->
cache
);
free
(
preCompContrib
->
cache
);
}
#if NIX
if
(
preCompContrib
->
numPhotonsShm
)
{
/* Release shared memory for photon counter */
munmap
(
preCompContrib
->
numPhotonsShm
,
sizeof
(
unsigned
long
));
close
(
preCompContrib
->
shmFile
);
unlink
(
preCompContrib
->
shmFname
);
}
#endif
#ifdef PMAP_CONTRIB_BINHISTO
free
(
preCompContrib
->
binHisto
);
#endif
}
/* Clean up precomputed contrib photon map */
deletePhotons
(
preCompContribPmap
);
free
(
preCompContribPmap
);
}
}
void
initPreComputedContrib
(
PreComputedContrib
*
preCompContrib
)
/* Initialise precomputed contribution container in photon map */
{
preCompContrib
->
waveletFname
=
NULL
;
preCompContrib
->
waveletFile
=
NULL
;
preCompContrib
->
waveletMatrix
=
preCompContrib
->
tWaveletMatrix
=
NULL
;
preCompContrib
->
compressedCoeffs
=
NULL
;
setcolor
(
preCompContrib
->
mrgbeRange
.
min
,
0
,
0
,
0
);
setcolor
(
preCompContrib
->
mrgbeRange
.
max
,
0
,
0
,
0
);
setcolor
(
preCompContrib
->
mrgbeRange
.
norm
,
0
,
0
,
0
);
preCompContrib
->
mrgbeCoeffs
=
NULL
;
preCompContrib
->
scDim
=
preCompContrib
->
nBins
=
preCompContrib
->
coeffDim
=
preCompContrib
->
nCoeffs
=
preCompContrib
->
nNonZeroCoeffs
=
preCompContrib
->
nCompressedCoeffs
=
preCompContrib
->
contribSize
=
0
;
preCompContrib
->
cache
=
NULL
;
preCompContrib
->
numPhotonsShm
=
NULL
;
preCompContrib
->
shmFile
=
-
1
;
#ifdef PMAP_CONTRIB_BINHISTO
preCompContrib
->
binHisto
=
NULL
;
#endif
}
static
int
allocPreCompContribPmap
(
const
LUENT
*
modContNode
,
void
*
p
)
/* Allocate & init child contrib photon map for modContNode's modifier
* and its contributions for parent photon map p. */
{
unsigned
scDim
,
nBins
,
coeffDim
,
nCoeffs
,
nCompressedCoeffs
,
nNzDetailCoeffs
;
PhotonMap
*
preCompContribPmap
,
*
parentPmap
=
(
PhotonMap
*
)
p
;
LUENT
*
preCompContribNode
;
PreComputedContrib
*
preCompContrib
;
const
MODCONT
*
modCont
=
(
MODCONT
*
)
modContNode
->
data
;
if
(
!
modCont
||
!
parentPmap
)
return
-
1
;
preCompContribNode
=
lu_find
(
parentPmap
->
preCompContribTab
,
modCont
->
modname
);
if
(
!
preCompContribNode
)
error
(
SYSTEM
,
"out of memory allocating LUT entry in "
"allocPreCompContribPmap()"
);
preCompContribNode
->
key
=
(
char
*
)
modCont
->
modname
;
preCompContribNode
->
data
=
(
char
*
)(
preCompContribPmap
=
malloc
(
sizeof
(
PhotonMap
))
);
if
(
!
preCompContribPmap
)
error
(
SYSTEM
,
"out of memory allocating photon map "
"in allocPreCompContribPmap()"
);
initPhotonMap
(
preCompContribPmap
,
PMAP_TYPE_CONTRIB_CHILD
);
/* Set output filename from parent photon map */
preCompContribPmap
->
fileName
=
parentPmap
->
fileName
;
/* Set contrib/coeff mode from parent photon map */
preCompContribPmap
->
contribMode
=
parentPmap
->
contribMode
;
preCompContrib
=
(
PreComputedContrib
*
)(
preCompContribPmap
->
preCompContrib
=
malloc
(
sizeof
(
PreComputedContrib
))
);
if
(
!
preCompContrib
)
error
(
SYSTEM
,
"out of memory allocating contributions "
"in allocPreCompContribPmap()"
);
initPreComputedContrib
(
preCompContrib
);
#if NIX
/* Allocate and map shared memory for photon counter */
strcpy
(
preCompContrib
->
shmFname
,
PMAP_TMPFNAME
);
preCompContrib
->
shmFile
=
mkstemp
(
preCompContrib
->
shmFname
);
if
(
preCompContrib
->
shmFile
<
0
||
ftruncate
(
preCompContrib
->
shmFile
,
sizeof
(
unsigned
long
))
<
0
)
error
(
SYSTEM
,
"failed opening shared memory file in allocPreCompContribPmap()"
);
preCompContrib
->
numPhotonsShm
=
mmap
(
NULL
,
sizeof
(
unsigned
long
),
PROT_READ
|
PROT_WRITE
,
MAP_SHARED
,
preCompContrib
->
shmFile
,
0
);
if
(
preCompContrib
->
numPhotonsShm
==
MAP_FAILED
)
error
(
SYSTEM
,
"failed mapping shared memory in allocPreCompContribPmap()"
);
#endif
/* Set Shirley-Chiu square & wavelet matrix dimensions
(number of bins resp. coefficients). */
preCompContrib
->
scDim
=
scDim
=
sqrt
(
preCompContrib
->
nBins
=
nBins
=
modCont
->
nbins
);
if
(
scDim
*
scDim
!=
nBins
)
/* Shouldn't happen; mkpmap took care of dis */
error
(
INTERNAL
,
"contribution bin count not integer square "
"in allocPreCompContribPmap()"
);
/* Allocate and initialise photon and contribution heaps */
initPhotonHeap
(
preCompContribPmap
);
if
(
nBins
>
1
)
{
/* BINNING ENABLED; set up wavelet xform & compression.
Get dimensions of wavelet coefficient matrix after
boundary padding, and number of nonzero coeffs. The
number of compressed (detail) coeffs is fixed and
based on the number of NONZERO coeffs (minus approx,
see NOTE below) since zero coeffs are considered
thresholded by default.
!!! NOTE: THE APPROX COEFFS IN THE UPPER LEFT OF
!!! THE MATRIX ARE _NEVER_ THRESHOLDED TO MINIMISE
!!! COMPRESSION ARTEFACTS! THEY ARE ESSENTIAL FOR
!!! RECONSTRUCTING THE ORIGINAL CONTRIBS. */
preCompContrib
->
coeffDim
=
coeffDim
=
padWaveletXform2
(
NULL
,
NULL
,
scDim
,
&
preCompContrib
->
nNonZeroCoeffs
);
preCompContrib
->
nCoeffs
=
nCoeffs
=
coeffDim
*
coeffDim
;
nNzDetailCoeffs
=
WAVELET_PADD2_NUMDETAIL
(
preCompContrib
->
nNonZeroCoeffs
);
nCompressedCoeffs
=
(
1
-
compressRatio
)
*
nNzDetailCoeffs
;
if
(
nCompressedCoeffs
<
1
)
{
error
(
WARNING
,
"maximum contribution compression, "
"clamping number of wavelet coefficients in "
"allocPreCompContribPmap()"
);
nCompressedCoeffs
=
1
;
}
if
(
nCompressedCoeffs
>=
preCompContrib
->
nCoeffs
)
{
error
(
WARNING
,
"minimum contribution compression, "
"clamping number of wavelet coefficients in "
"allocPreCompContribPmap()"
);
nCompressedCoeffs
=
nNzDetailCoeffs
;
}
preCompContrib
->
nCompressedCoeffs
=
nCompressedCoeffs
;
/* Lazily allocate primary & transposed wavelet coeff matrix */
preCompContrib
->
waveletMatrix
=
allocWaveletMatrix2
(
coeffDim
);
preCompContrib
->
tWaveletMatrix
=
allocWaveletMatrix2
(
coeffDim
);
if
(
!
preCompContrib
->
waveletMatrix
||
!
preCompContrib
->
tWaveletMatrix
)
error
(
SYSTEM
,
"out of memory allocating wavelet coefficient "
"matrix in allocPreCompContribPmap()"
);
/* Lazily allocate thresholded detail coefficients */
preCompContrib
->
compressedCoeffs
=
calloc
(
nNzDetailCoeffs
,
sizeof
(
PreComputedContribCoeff
)
);
/* Lazily alloc mRGBE-encoded compressed wavelet coeffs */
preCompContrib
->
mrgbeCoeffs
=
calloc
(
nCompressedCoeffs
,
sizeof
(
mRGBE
)
);
if
(
!
preCompContrib
->
compressedCoeffs
||
!
preCompContrib
->
mrgbeCoeffs
)
error
(
SYSTEM
,
"out of memory allocating compressed "
"contribution buffer in allocPreCompContribPmap()"
);
/* Set size of compressed contributions, in bytes */
preCompContrib
->
contribSize
=
PMAP_CONTRIB_ENCSIZE
(
nCompressedCoeffs
);
#ifdef PMAP_CONTRIB_BINHISTO
preCompContrib
->
binHisto
=
calloc
(
nBins
,
sizeof
(
unsigned
long
));
if
(
!
preCompContrib
->
binHisto
)
error
(
SYSTEM
,
"out of memory allocating contribution bin "
"histogram in allocPreCompContribPmap()"
);
else
memset
(
preCompContrib
->
binHisto
,
0
,
(
size_t
)
nBins
*
sizeof
(
unsigned
long
)
);
preCompContrib
->
minZero
=
nBins
;
preCompContrib
->
maxZero
=
preCompContrib
->
sumZero
=
0
;
#endif
}
}
static
int
sumNumPhotons
(
const
LUENT
*
preCompContribNode
,
void
*
p
)
/* Set number of photons for this modifier's child contrib photon map
* from the corresponding shared photon counter, and sum over all
* modifiers in (unsigned long*)p */
{
unsigned
long
*
numPhotons
=
(
unsigned
long
*
)
p
;
PhotonMap
*
preCompContribPmap
=
(
PhotonMap
*
)(
preCompContribNode
->
data
);
const
PreComputedContrib
*
preCompContrib
;
if
(
!
preCompContribPmap
)
return
0
;
else
preCompContrib
=
(
PreComputedContrib
*
)(
preCompContribPmap
->
preCompContrib
);
#if NIX
/* Lock shared photon counter file for reading */
shmLock
(
preCompContrib
->
shmFile
,
F_RDLCK
);
#endif
*
numPhotons
+=
(
preCompContribPmap
->
numPhotons
=
*
preCompContrib
->
numPhotonsShm
);
#if NIX
/* Release lock on shard photon counter file */
shmLock
(
preCompContrib
->
shmFile
,
F_UNLCK
);
#endif
return
0
;
}
void
preComputeContrib
(
PhotonMap
*
pmap
,
unsigned
numProc
)
/* Precompute contributions for a random subset of (finalGather *
pmap -> numPhotons) photons, and return the per-modifier precomputed
contribution photon maps in LUT preCompContribTab, discarding the
original photons. */
{
unsigned
long
p
,
numPreComp
;
unsigned
proc
;
int
stat
;
pid_t
pid
,
procPids
[
PMAP_MAXPROC
];
PhotonMap
nuPmap
;
LUTAB
lutInit
=
LU_SINIT
(
NULL
,
freePreCompContribNode
);
/* Init new parent photon map and set output filename */
initPhotonMap
(
&
nuPmap
,
pmap
->
type
);
nuPmap
.
fileName
=
pmap
->
fileName
;
/* Set contrib/coeff mode in new parent photon map */
nuPmap
.
contribMode
=
pmap
->
contribMode
;
/* Allocate and init LUT containing per-modifier child photon maps */
if
(
!
(
nuPmap
.
preCompContribTab
=
malloc
(
sizeof
(
LUTAB
))))
error
(
SYSTEM
,
"out of memory allocating LUT in preComputeContrib()"
);
memcpy
(
nuPmap
.
preCompContribTab
,
&
lutInit
,
sizeof
(
LUTAB
));
/* Allocate and init per-modifier child photon maps */
lu_doall
(
pmapContribTab
,
allocPreCompContribPmap
,
&
nuPmap
);
/* Record start time, baby */
repStartTime
=
repLastTime
=
time
(
NULL
);
repComplete
=
numPreComp
=
finalGather
*
pmap
->
numPhotons
;
repProgress
=
0
;
if
(
verbose
)
{
sprintf
(
errmsg
,
"
\n
Precomputing contributions for %ld photons
\n
"
,
numPreComp
);
eputs
(
errmsg
);
#if NIX
fflush
(
stderr
);
#endif
}
for
(
proc
=
0
;
proc
<
numProc
;
proc
++
)
{
#if NIX
if
(
!
(
pid
=
fork
()))
{
/* SUBPROCESS ENTERS HERE; opened and mmapped files inherited */
#else
if
(
1
)
{
/* No subprocess under Windoze */
#endif
int
photonOk
;
unsigned
i
,
j
,
k
,
coeffIdx
,
scDim
,
nBins
,
coeffDim
,
nCoeffs
,
nCompressedCoeffs
;
const
double
pInc
=
finalGather
>
FTINY
?
1
/
finalGather
:
0
,
numPhotonsPerProc
=
(
double
)
pmap
->
numPhotons
/
numProc
,
numPreCompPerProc
=
(
double
)
numPreComp
/
numProc
;
PhotonIdx
pIdx
;
Photon
photon
;
MODCONT
*
modCont
;
LUENT
*
preCompContribNode
;
PhotonMap
*
preCompContribPmap
;
PreComputedContrib
*
preCompContrib
;
WaveletMatrix2
waveletMatrix
;
WaveletCoeff3
*
coeffPtr
;
RAY
ray
;
unsigned
long
lastpIdx
=
0
;
photonRay
(
NULL
,
&
ray
,
PRIMARY
,
NULL
);
ray
.
ro
=
NULL
;
/* Decorrelate subprocess RNGs */
pmapSeed
(
randSeed
+
proc
,
pmap
->
randState
);
#ifdef PMAP_CONTRIB_DBG
/* Output child process PID after random delay to prevent
* scrambling console output with sibling subprocs */
usleep
(
1e6
*
pmapRandom
(
rouletteState
));
sprintf
(
errmsg
,
"*DBG* Proc %d: PID = %d "
"(waiting 10 sec to attach DEBUGGA...)
\n
"
,
proc
,
getpid
()
);
eputs
(
errmsg
);
/* Give debugga time to attach to subproc */
sleep
(
10
);
#endif
for
(
p
=
0
;
p
<
(
unsigned
long
)
numPreCompPerProc
;
p
++
)
{
do
{
/* Get random photon from stratified distribution in
* source heap to avoid duplicates and clustering. */
pIdx
=
firstPhoton
(
pmap
)
+
(
unsigned
long
)(
p
*
pInc
+
pmapRandom
(
pmap
->
randState
)
*
max
(
pInc
-
1
,
0
)
+
proc
*
numPhotonsPerProc
);
getPhoton
(
pmap
,
pIdx
,
&
photon
);
#ifdef PMAP_CONTRIB_DBG
if
(
pIdx
==
lastpIdx
)
{
sprintf
(
errmsg
,
"duplicate contrib photon "
"(proc %d, pIdx %d)"
,
proc
,
pIdx
);
error
(
WARNING
,
errmsg
);
}
lastpIdx
=
pIdx
;
#endif
/* Init dummy photon ray with intersection and normal at
* photon position. Set emitting light source index from
* origin. */
VCOPY
(
ray
.
rop
,
photon
.
pos
);
ray
.
rsrc
=
photonSrcIdx
(
pmap
,
&
photon
);
for
(
i
=
0
;
i
<
3
;
i
++
)
ray
.
ron
[
i
]
=
photon
.
norm
[
i
]
/
127.0
;
/* Get contributions at photon position; retry with
another photon if no contribs found */
modCont
=
getPhotonContrib
(
pmap
,
&
ray
,
ray
.
rcol
);
}
while
(
!
modCont
);
/* Get preallocated LUT node for this modifier */
preCompContribNode
=
lu_find
(
nuPmap
.
preCompContribTab
,
modCont
->
modname
);
if
(
!
preCompContribNode
||
!
preCompContribNode
->
key
||
!
preCompContribNode
->
data
)
error
(
CONSISTENCY
,
"missing LUT entry in preComputeContrib()"
);
preCompContribPmap
=
(
PhotonMap
*
)
preCompContribNode
->
data
;
preCompContrib
=
(
PreComputedContrib
*
)(
preCompContribPmap
->
preCompContrib
);
scDim
=
preCompContrib
->
scDim
;
nBins
=
preCompContrib
->
nBins
;
nCoeffs
=
preCompContrib
->
nCoeffs
;
nCompressedCoeffs
=
preCompContrib
->
nCompressedCoeffs
;
waveletMatrix
=
preCompContrib
->
waveletMatrix
;
if
(
nBins
>
1
)
{
/* Binning enabled */
char
contribBuf0
[
preCompContrib
->
contribSize
],
*
contribBuf
;
COLR
*
rgbe32
=
(
COLR
*
)
contribBuf0
;
/* copy binned contribs to wavelet
* input matrix before compressing and encoding */
for
(
i
=
0
;
i
<
scDim
;
i
++
)
{
/* Calculate row pointer (=multiplication) in outer
* loop, increment in inner loop */
coeffPtr
=
modCont
->
cbin
+
PMAP_CONTRIB_XY2LIN
(
scDim
,
i
,
0
);
#ifdef PMAP_CONTRIB_LOG
/* Logarithmic contribs; copy & apply per element */
for
(
j
=
0
;
j
<
scDim
;
j
++
,
coeffPtr
++
)
for
(
k
=
0
;
k
<
3
;
k
++
)
waveletMatrix
[
i
]
[
j
]
[
k
]
=
PMAP_CONTRIB_VAL
((
*
coeffPtr
)
[
k
]);
#else
/* Linear contribs; copy matrix rows */
memcpy
(
waveletMatrix
[
i
],
*
coeffPtr
,
scDim
*
WAVELET_COEFFSIZE
);
#endif
}
/* Compress and encode contribs */
if
(
encodeContribs
(
preCompContrib
,
compressRatio
))
error
(
INTERNAL
,
"failed contribution "
"compression/encoding in preComputeContrib()"
);
/* Encode wavelet approx coeffs in the upper left of
* the wavelet matrix to 32-bit RGBE, and buffer them.
* These coeffs are not thresholded to minimise
* compression artefacts. */
for
(
i
=
0
;
i
<
WAVELET_PADD2_APPROXDIM
;
i
++
)
{
for
(
j
=
0
;
j
<
WAVELET_PADD2_APPROXDIM
;
j
++
,
rgbe32
++
)
{
setcolr
(
*
rgbe32
,
fabs
(
waveletMatrix
[
i
]
[
j
]
[
0
]),
fabs
(
waveletMatrix
[
i
]
[
j
]
[
1
]),
fabs
(
waveletMatrix
[
i
]
[
j
]
[
2
])
);
/* HACK: depending on the boundary extension mode,
* some approx coeffs can be NEGATIVE (!), which
* 32-bit RGBE doesn't accommodate. To get around
* this, we sacrifice a bit in each mantissa for the
* sign. */
for
(
k
=
0
;
k
<
3
;
k
++
)
(
*
rgbe32
)
[
k
]
=
PMAP_CONTRIB_SET_RGBE32_SGN
(
(
*
rgbe32
)
[
k
],
waveletMatrix
[
i
]
[
j
]
[
k
]
);
}
}
/* Encode wavelet detail coeff range to 32-bit RGBE, and
* buffer it after approx coeffs.
* NOTE: mRGBE range maximum is buffered first to
* facilitate handling the special (but pointless) case of
* a single compressed detail coeff; if so, the mRGBE
* minimum is omitted, since implicitly zero. */
setcolr
(
*
(
rgbe32
++
),
preCompContrib
->
mrgbeRange
.
max
[
0
],
preCompContrib
->
mrgbeRange
.
max
[
1
],
preCompContrib
->
mrgbeRange
.
max
[
2
]
);
if
(
nCompressedCoeffs
>
1
)
setcolr
(
*
(
rgbe32
++
),
preCompContrib
->
mrgbeRange
.
min
[
0
],
preCompContrib
->
mrgbeRange
.
min
[
1
],
preCompContrib
->
mrgbeRange
.
min
[
2
]
);
contribBuf
=
(
char
*
)
rgbe32
;
/* Now buffa mRGBE encoded, compressed detail coeffs.
* NOTE: we can't copy to preCompContribPmap ->
* contribHeapBuf directly because it isn't allocated
* until the 1st call to newPhoton() below. IMHO this is
* cleaner and more consistent with the existing code than
* allocating it here. */
memcpy
(
contribBuf
,
preCompContrib
->
mrgbeCoeffs
,
nCompressedCoeffs
*
sizeof
(
mRGBE
)
);
/* Append precomputed photon to heap from ray, check if
* photon was accepted (i.e. newPhoton() returned 0).
* Note the locally buffered encoded contribs are passed
* as an extra parameter. */
ray
.
rsrc
=
proc
;
photonOk
=
!
newPhoton
(
preCompContribPmap
,
&
ray
,
contribBuf0
);
#ifdef PMAP_CONTRIB_BINHISTO
{
unsigned
long
numZero
=
0
;
/* Update bin histogram for last photon lookup */
for
(
i
=
1
;
i
<
pmap
->
squeue
.
tail
;
i
++
)
{
Photon
*
photon
=
getNearestPhoton
(
&
pmap
->
squeue
,
pmap
->
squeue
.
node
[
i
].
idx
);
++
preCompContrib
->
binHisto
[
photonSrcBin
(
pmap
,
photon
)
];
}
for
(
i
=
0
;
i
<
preCompContrib
->
nBins
;
i
++
)
numZero
+=
(
intens
(
modCont
->
cbin
[
i
])
<
FTINY
);
preCompContrib
->
sumZero
+=
numZero
;
if
(
numZero
<
preCompContrib
->
minZero
)
preCompContrib
->
minZero
=
numZero
;
if
(
numZero
>
preCompContrib
->
maxZero
)
preCompContrib
->
maxZero
=
numZero
;
}
#endif
}
else
{
/* Binning disabled; append precomputed photon to heap
* from ray _without_ encoded contribs, check if photon
* was accepted (i.e. newPhoton() returned 0). */
photonOk
=
!
newPhoton
(
preCompContribPmap
,
&
ray
,
NULL
);
}
/* Increment parent photon map's photon counter if photon was
* accepted by newPhoton(). */
nuPmap
.
numPhotons
+=
photonOk
;
#if NIX
/* Asynchronously increment shared photon counter for this
* modifier if photon was accepted by newPhoton().
* NOTE: write lock on memmapped file! */
if
(
photonOk
)
{
shmLock
(
preCompContrib
->
shmFile
,
F_WRLCK
);
*
preCompContrib
->
numPhotonsShm
+=
photonOk
;
shmLock
(
preCompContrib
->
shmFile
,
F_UNLCK
);
}
#else
/* Synchronously report progress on W33nD0z */
if
(
!
proc
&&
photonRepTime
>
0
&&
time
(
NULL
)
>=
repLastTime
+
photonRepTime
)
{
repProgress
=
nuPmap
.
numPhotons
;
pmapPreCompReport
();
}
#endif
}
/* End precomp. photon loop */
/* Flush heap buffa one final time to prevent data corruption */
flushPhotonHeap
(
preCompContribPmap
);
#ifdef PMAP_CONTRIB_BINHISTO
/* Print binning histogram/population count for 1st subproc */
if
(
!
proc
)
{
printf
(
"Binning histogram for modifier %s
\n
:"
"Avg photons / bin:
\n
"
,
preCompContribNode
->
key
);
for
(
i
=
0
;
i
<
preCompContrib
->
nBins
;
i
++
)
printf
(
"%f
\n
"
,
(
float
)
preCompContrib
->
binHisto
[
i
]
/
preCompContrib
->
nBins
/
nuPmap
.
numPhotons
);
printf
(
"Empty bin ratio (avg, min, max):
\t
%f
\t
%f
\t
%f
\n
"
,
(
float
)
preCompContrib
->
sumZero
/
preCompContrib
->
nBins
/
nuPmap
.
numPhotons
,
(
float
)
preCompContrib
->
minZero
/
preCompContrib
->
nBins
,
(
float
)
preCompContrib
->
maxZero
/
preCompContrib
->
nBins
);
}
#endif
#if NIX
#ifdef PMAP_CONTRIB_DBG
if
(
mrgbeDiffCnt
)
{
sprintf
(
errmsg
,
"*DBG* preComputeContrib: Avg mRGBE encoding diff = %g
\n
"
,
mrgbeDiffs
/
mrgbeDiffCnt
);
eputs
(
errmsg
);
}
sprintf
(
errmsg
,
"*DBG* Subproc %d exiting normally
\n
"
,
proc
);
eputs
(
errmsg
);
#endif
/* Terminate subprocess normally */
exit
(
0
);
#endif
}
/* if (!fork()) */
else
{
/* PARENT PROCESS CONTINUES FORKING LOOP HERE */
if
(
pid
<
0
)
error
(
SYSTEM
,
"failed to fork subprocess in preComputeContrib()"
);
else
/* Save child process ID */
procPids
[
proc
]
=
pid
;
}
}
/* End forking loop */
#if NIX
/* PARENT PROCESS CONTINUES HERE AFTER FORKING ALL SUBPROCS */
#ifdef SIGCONT
/* Enable progress report signal handler */
signal
(
SIGCONT
,
pmapPreCompReport
);
#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 & quit */
for
(
proc
=
0
;
proc
<
numProc
;
proc
++
)
kill
(
procPids
[
proc
],
SIGKILL
);
cleanUpPmaps
(
photonMaps
);
error
(
USER
,
"failed precomputing contribution photons"
);
}
--
proc
;
}
/* Nod off for a bit and update progress */
sleep
(
1
);
/* Update num photons from shared photon counters for all mods */
nuPmap
.
numPhotons
=
0
;
lu_doall
(
nuPmap
.
preCompContribTab
,
sumNumPhotons
,
&
nuPmap
.
numPhotons
);
repProgress
=
nuPmap
.
numPhotons
;
if
(
photonRepTime
>
0
&&
time
(
NULL
)
>=
repLastTime
+
photonRepTime
)
pmapPreCompReport
();
#ifdef SIGCONT
else
signal
(
SIGCONT
,
pmapPreCompReport
);
#endif
}
#endif
/* NIX */
#ifdef SIGCONT
/* Reset signal handler */
signal
(
SIGCONT
,
SIG_DFL
);
#endif
/* Trash original pmap and its leaf file, and replace with new parent,
which now acts as container for the per-modifier child pmaps */
unlink
(
pmap
->
store
.
leafFname
);
deletePhotons
(
pmap
);
memcpy
(
pmap
,
&
nuPmap
,
sizeof
(
PhotonMap
));
/* Bail out if no photons could be precomputed */
if
(
!
pmap
->
numPhotons
)
error
(
USER
,
"no contribution photons precomputed; try increasing -am"
);
/* !!! NOTE: buildPreCompContribPmap() bails out if any child contrib
* !!! pmaps are empty -- should this trigger a warning instead? */
/* 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 <stdio.h>
#include <random.h>
#include "func.h"
int
main
(
int
argc
,
char
*
argv
[])
{
unsigned
i
,
scdim
,
nsamp
,
numTheta
=
0
,
numPhi
=
0
;
double
t
,
p
;
RAY
ray
;
int
bin
;
if
(
argc
<
3
)
{
fprintf
(
stderr
,
"%s <scdim> <nsamp> [<var>=<value>; ..]
\n
"
,
argv
[
0
]
);
fprintf
(
stderr
,
"Default variable defs: %s
\n
"
,
PMAP_CONTRIB_SCDEFAULTS
);
fputs
(
"
\n
Missing Shirley-Chiu dimension scDim>1, "
"number of samples nsamp
\n
"
,
stderr
);
return
-
1
;
}
if
(
!
(
scdim
=
atoi
(
argv
[
1
])))
{
fputs
(
"Invalid Shirley-Chiu dimension
\n
"
,
stderr
);
return
-
1
;
}
if
(
!
(
nsamp
=
atoi
(
argv
[
2
])))
{
fputs
(
"Invalid num samples
\n
"
,
stderr
);
return
-
1
;
}
else
{
numTheta
=
(
int
)(
sqrt
(
nsamp
)
/
2
);
numPhi
=
4
*
numTheta
;
}
/* (Doan' need to) Init cal func routines for binning */
#if 0
initfunc();
calcontext(RCCONTEXT);
#endif
/* Compile default orientation variables for contribution binning */
scompile
(
PMAP_CONTRIB_SCDEFAULTS
,
NULL
,
0
);
/* Compile custom orientation variables from command line */
for
(
i
=
3
;
i
<
argc
;
i
++
)
scompile
(
argv
[
i
],
NULL
,
0
);
puts
(
"x
\t
y
\t
z
\t\t
theta
\t
phi
\t\t
bin"
);
for
(
i
=
0
;
i
<
nsamp
;
i
++
)
{
#if 0
/* Random */
t = 0.5 * PI * drand48();
p = 2 * PI * drand48();
#else
/* Stratified */
t
=
0.5
*
PI
*
((
i
%
numTheta
)
+
drand48
())
/
numTheta
;
p
=
2.0
*
PI
*
((
i
/
numTheta
)
+
drand48
())
/
numPhi
;
#endif
ray
.
rdir
[
0
]
=
sin
(
t
)
*
cos
(
p
);
ray
.
rdir
[
1
]
=
sin
(
t
)
*
sin
(
p
);
ray
.
rdir
[
2
]
=
cos
(
t
);
bin
=
ray2bin
(
&
ray
,
scdim
);
printf
(
"%.3f
\t
%.3f
\t
%.3f
\t
-->
\t
%.3f
\t
%.3f
\t
-->
\t
%d
\n
"
,
ray
.
rdir
[
0
],
ray
.
rdir
[
1
],
ray
.
rdir
[
2
],
t
,
p
,
bin
);
}
return
0
;
}
#endif
#endif
/* PMAP_CONTRIB */
Event Timeline
Log In to Comment