Homec4science
Diffusion RADIANCE Photon Map (pmapmerge) Browse

RadiancePmap/mrgbe.c
246a688ea4c0pmapmerge

Last Change
Raw File

mrgbe.c
View Options

/*
=========================================================================
mRGBE (miniRGBE): reduced RGBE encoding of normalised RGB floats plus
associated integer payload data. See mrgbe.h for encoding details.
Compile with -DmRGBE_TEST to enable unit test (requires libm)
Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
(c) Lucerne University of Applied Sciences and Arts,
supported by the Swiss National Science Foundation
(SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment")
=========================================================================
$Id$
*/
#include "mrgbe.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b))
#define ABS(a) ((a) >= 0 ? (a) : -(a))
mRGBERange *mRGBEinit (mRGBERange *range,
double rgbMin [3], double rgbMax [3]
)
/* Initialise range state for mRGBEencode()/mRGBEdecode() with the specified
* float input ranges [rgbMin, rgbMax] in _absolute_ values per colour
* channel. This sets the offset and normalisation for the encoding. */
{
unsigned i;
double orgMin;
for (i = 0; i < 3; i++) {
range -> min [i] = orgMin = fabs(rgbMin [i]);
range -> max [i] = fabs(rgbMax [i]);
if (rgbMin [i] < rgbMax [i]) {
/* Iteratively reduce minimum and set normalisation factor until
* normalising the original minimum exceeds mRGBE_MIN (and therefore
* won't be clamped to zero). This avoids sign reversals when
* encoding negative values close or equal to the original minimum. */
do {
range -> min [i] = MAX(0, range -> min [i] - mRGBE_MIN);
range -> norm [i] = mRGBE_RANGE /
(range -> max [i] - range -> min [i]);
} while (range -> min [i] > 0 &&
(orgMin - range -> min [i]) * range -> norm [i] <= mRGBE_MIN
);
if (range -> norm [i] < mRGBE_MIN) {
fprintf(stderr,
"ERROR - Zero mRGBE normalisation / range overflow\n"
);
return NULL;
}
}
else if (rgbMin [i] > rgbMax [i]) {
/* Min/max reversed */
fprintf(stderr, "ERROR - invalid mRGBE range\n");
return NULL;
}
else {
/* Min == max (singular range) */
range -> norm [i] = 1;
}
}
return range;
}
mRGBE mRGBEencode (double rgb [3], const mRGBERange *range, unsigned data)
/* Encode signed float RGB tuple to mRGBE along with payload data. The
* supplied range state must be previously initialised with mRGBEinit().
* A warning is issued if rgb lies outside the initialised range. */
{
#define SGN(a, b) ((a) >= 0 ? (b) : -(b))
mRGBE mrgbe = {0, 0, 0, 0};
double rgbNorm [3], mantNorm;
int i, exp;
/* Check payload data */
if (data > mRGBE_DATAMAX) {
fprintf(stderr, "ERROR - mRGBE data overflow\n");
return mrgbe;
}
if (!range) {
fprintf(stderr, "ERROR - mRGBE range not initialised\n");
return mrgbe;
}
/* Offset and normalise RGB according to initialised range; note rgbNorm
* is in ABSOLUTE VALUES */
for (i = 0; i < 3; i++) {
rgbNorm [i] = fabs(rgb [i]);
/* Check range */
if (rgbNorm [i] < range -> min [i] || rgbNorm [i] > range -> max [i]) {
fprintf(stderr, "ERROR - mRGBE input outside range\n");
return mrgbe;
}
rgbNorm [i] = (rgbNorm [i] - range -> min [i]) * range -> norm [i];
if (fabs(rgbNorm [i] - mRGBE_MAX) < mRGBE_MIN)
/* Normalisation to maximum. This must be avoided by offsetting
by minimum, otherwise frexp() returns a positive exponent.
Normalisation to >mRGBE_MAX is caught by overflow check below */
rgbNorm [i] = mRGBE_RANGE;
}
/* Mantissa normalisation factor (= RGB maximum) */
mantNorm = MAX(MAX(rgbNorm [0], rgbNorm [1]), rgbNorm [2]);
/* Special handling for zero */
if (mantNorm < mRGBE_MIN)
mantNorm = exp = 0;
else
/* Get normalised mantissa and exponent to base 2 */
mantNorm = frexp(mantNorm, &exp) * mRGBE_MANTMAX / mantNorm;
/* Set RGB integer mantissas, adding excess (zero offset) and rounding
up or down, depending on sign. Set sign according to original value.
Clamp to mRGBE_MANTOVF to prevent rounding beyond overflow into
negative! */
mrgbe.red = MIN((int)(
SGN(rgb [0], rgbNorm [0] * mantNorm + mRGBE_ROUND) + mRGBE_MANTMAX
),
mRGBE_MANTOVF
);
mrgbe.grn = MIN((int)(
SGN(rgb [1], rgbNorm [1] * mantNorm + mRGBE_ROUND) + mRGBE_MANTMAX
),
mRGBE_MANTOVF
);
mrgbe.blu = MIN((int)(
SGN(rgb [2], rgbNorm [2] * mantNorm + mRGBE_ROUND) + mRGBE_MANTMAX
),
mRGBE_MANTOVF
);
/* Drop the exponent sign, since implicitly negative */
mrgbe.exp = abs(exp);
/* Set payload data */
mrgbe.dat = data;
return mrgbe;
#undef SGN
}
unsigned mRGBEdecode (mRGBE mrgbe, const mRGBERange *range, double rgb [3])
/* Decode mRGBE, returning signed float RGB tuple in array rgb. The
* associated payload data is decoded as return value. The supplied range
* state must be previously initialised with mRGBEinit(). */
{
#define SGN(a, b) ((a) >= mRGBE_MANTMAX ? (b) : -(b))
#ifdef mRGBE_ZEROJITTER
#define JITTER(m, e) (mRGBE_JITTER * drand48())
#else
#define JITTER(m, e) ((m) == mRGBE_MANTMAX && !(e) \
? 0 : mRGBE_JITTER * drand48() \
)
#endif
double mantNorm;
if (!range) {
fprintf(stderr, "ERROR - mRGBE range not initialised\n");
return 0;
}
/* Mantissa normalisation factor; exponent is implicitly negative */
mantNorm = ldexp(1.0, -mrgbe.exp) / mRGBE_MANTMAX;
/* Decode and set RGB components, subtracting excess (zero offset) and
* adding jitter to break up quantisation artefacts. Clamp to -mRGBE_MAX
* in case jitter underflows. (Note jitter is added/subtracted for
* negative/positive values to compensate rounding down/up during
* encoding). Finally denormalise and offset to original range with
* saved state. */
rgb [0] = SGN(mrgbe.red, 1) * (
ABS(mrgbe.red - mRGBE_MANTMAX + JITTER(mrgbe.red, mrgbe.exp)) *
mantNorm / range -> norm [0] +
range -> min [0]
);
rgb [1] = SGN(mrgbe.grn, 1) * (
ABS(mrgbe.grn - mRGBE_MANTMAX + JITTER(mrgbe.grn, mrgbe.exp)) *
mantNorm / range -> norm [1] +
range -> min [1]
);
rgb [2] = SGN(mrgbe.blu, 1) * (
ABS(mrgbe.blu - mRGBE_MANTMAX + JITTER(mrgbe.blu, mrgbe.exp)) *
mantNorm / range -> norm [2] +
range -> min [2]
);
/* Return payload data */
return mrgbe.dat;
#undef SGN
#undef JITTER
}
#ifdef mRGBE_TEST
/* Build standalone to run unit test */
#define RGBMIN 0
#define RGBMAX 1
int main (int argc, char **argv)
{
double rgb [3] = {0, 0, 0}, rgbDec [3],
rgbMin [3] = {0, 0, 0}, rgbMax [3] = {1, 1, 1},
rgbAvg [3] = {0, 0, 0}, rgbDecAvg [3] = {0, 0, 0},
rgbAbs, relDiff, rmse, rmseSum;
unsigned numTests, i, j, k, data, dataDec;
const short rndState0 [3] = {11, 23, 31};
short encFlag = 0, rndState [3];
mRGBERange mrgbeRange;
mRGBE mrgbe;
if (argc > 1 && (numTests = atoi(argv [1])) > 0) {
/* Test range checks */
mRGBEinit(&mrgbeRange, rgbMin, rgbMax);
puts("Check zero");
rgb [0] = rgb [1] = rgb [2] = 0;
mrgbe = mRGBEencode(rgb, &mrgbeRange, 0);
mRGBEdecode(mrgbe, &mrgbeRange, rgbDec);
printf("% 7.3f\t% 7.3f\t% 7.3f\t-->\t%02d:%02d:%02d:%02d\t-->"
"\t% 7.3f\t% 7.3f\t%7.3f\n",
rgb [0], rgb [1], rgb [2],
mrgbe.red, mrgbe.grn, mrgbe.blu, mrgbe.exp,
rgbDec [0], rgbDec [1], rgbDec [2]
);
puts("Check data overflow");
rgb [0] = rgb [1] = rgb [2] = mRGBE_MAX / 2;
mRGBEencode(rgb, &mrgbeRange, mRGBE_DATAMAX << 1);
puts("\nCheck overflow (positive)");
rgb [0] = rgb [1] = rgb [2] = 2 * mRGBE_MAX;
mRGBEencode(rgb, &mrgbeRange, 0);
puts("\nCheck overflow (negative)");
rgb [0] = rgb [1] = rgb [2] = -2 * mRGBE_MAX;
mRGBEencode(rgb, &mrgbeRange, 0);
puts("\nCheck singular range");
rgb [0] = rgb [1] = rgb [2] = RGBMAX * erand48(rndState);
mRGBEinit(&mrgbeRange, rgb, rgb);
mrgbe = mRGBEencode(rgb, &mrgbeRange, 0);
mRGBEdecode(mrgbe, &mrgbeRange, rgbDec);
printf("% 7.3f\t% 7.3f\t% 7.3f\t-->\t%02d:%02d:%02d:%02d\t-->"
"\t% 7.3f\t% 7.3f\t%7.3f\n",
rgb [0], rgb [1], rgb [2],
mrgbe.red, mrgbe.grn, mrgbe.blu, mrgbe.exp,
rgbDec [0], rgbDec [1], rgbDec [2]
);
puts("\nCheck invalid range");
mRGBEinit(&mrgbeRange, rgbMax, rgbMin);
putchar('\n');
rgbMin [0] = rgbMin [1] = rgbMin [2] = RGBMAX;
rgbMax [0] = rgbMax [1] = rgbMax [2] = RGBMIN;
/* Iterate twice with same random RGB tuples; gather minimum and
* maximum in 1st pass, then initialise and run encode/decode tests
* in 2nd pass. */
do {
/* Initialise RNG state so this can be restored to obtain the
same RGB tuple sequence in both passes. This state needs to
be kept separate from the RNG used for jittering in the
encode/decode routines. */
for (i = 0; i < 3; i++)
rndState [i] = rndState0 [i];
/* Random encode vs. decode */
for (k = 0, rmseSum = 0; k < numTests; k++) {
data = (unsigned)(mRGBE_DATAMAX * erand48(rndState));
for (i = 0; i < 3; i++) {
rgb [i] = (erand48(rndState) > 0.5 ? 1 : -1) *
((RGBMAX - RGBMIN) * erand48(rndState) + RGBMIN);
#if 1
/* Correlate RGB */
if (i > 0)
rgb [i] = rgb [0] * (0.1 + 0.9 * erand48(rndState));
#endif
#if 0
/* Check zeros */
rgb [i] = 0;
#endif
}
if (encFlag) {
printf("% 4d:\t", k);
/* Encode */
mrgbe = mRGBEencode(rgb, &mrgbeRange, data);
printf("% 7.3f\t% 7.3f\t% 7.3f\t%4d\t-->\t"
"%02d:%02d:%02d:%02d\t-->\t",
rgb [0], rgb [1], rgb [2], data,
mrgbe.red, mrgbe.grn, mrgbe.blu, mrgbe.exp
);
/* Decode */
dataDec = mRGBEdecode(mrgbe, &mrgbeRange, rgbDec);
for (j = rmse = 0; j < 3; j++) {
relDiff = fabs(rgb [j]) > 0
? 1 - rgbDec [j] / rgb [j] : 0;
rmse += relDiff * relDiff;
/* Update original and decoded averages */
rgbAvg [j] += rgb [j];
rgbDecAvg [j] += rgbDec [j];
}
rmseSum += rmse = sqrt(rmse / 3);
printf("% 7.3f\t% 7.3f\t% 7.3f\t%4d\tRMSE = %.1f%% %s\n",
rgbDec [0], rgbDec [1], rgbDec [2],
dataDec, 100 * rmse, rmse > 0.1 ? "!" : ""
);
#if 1
/* Stop at sign reversal; this error isn't acceptable
* under any circumstances! */
for (j = 0; j < 3; j++)
if (rgb [j] / rgbDec [j] < 0) {
puts("**** SIGN REVERSAL!!! ****\n");
return -1;
}
#endif
}
else {
/* Update RGB minimum and maximum */
for (i = 0; i < 3; i++) {
rgbAbs = fabs(rgb [i]);
if (rgbAbs > rgbMax [i])
rgbMax [i] = rgbAbs;
if (rgbAbs < rgbMin [i])
rgbMin [i] = rgbAbs;
#if 0
/* NOTE: Clamping range minimum to 0 increases RMSE
if RGBMIN > 0 */
rgbMin [i] = 0;
#endif
}
}
}
/* Dump results */
if (!encFlag) {
/* Initialise mRGBE encoding state with min/max */
if (!mRGBEinit(&mrgbeRange, rgbMin, rgbMax))
return -1;
printf(
"RGB Range = [%.3f, %.3f, %.3f] - [%.3f, %.3f, %.3f]\n\n",
rgbMin [0], rgbMin [1], rgbMin [2],
rgbMax [0], rgbMax [1], rgbMax [2]
);
}
else {
putchar(' ');
for (i = 0; i < 118; i++)
putchar('-');
printf("\n Avg:\t% 7.3f\t% 7.3f\t% 7.3f\t\t\t\t\t\t"
"% 7.3f\t% 7.3f\t% 7.3f\n",
rgbAvg [0] / numTests, rgbAvg [1] / numTests,
rgbAvg [2] / numTests, rgbDecAvg [0] / numTests,
rgbDecAvg [1] / numTests, rgbDecAvg [2] / numTests
);
printf("\nAvg RMSE = %.2f%%\n", 100 * rmseSum / numTests);
}
} while (encFlag = !encFlag);
return 0;
}
else {
printf("Usage: %s <numTests>\n", argv [0]);
return -1;
}
}
#endif
c4science · Help