Homec4science
Diffusion RADIANCE Photon Map (pmapmerge) Browse

RadiancePmap/rgbecontrib.c
bca8d07c4d81pmapmerge

Last Change
Raw File

rgbecontrib.c
View Options

/*
=========================================================================
Routines to encode/decode RGBE representation of normalised binned
contributions. See rgbecontrib.h for encoding details.
Compile with -D_RGBE_CONTRIB_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 "rgbecontrib.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))
RGBEContrib RGBEContribEnc (double rgb [3], double norm [3], unsigned bin)
/* Normalise signed float RGB contributions against abs(norm), and encode
* to RGBE with associated bin. Issues warning if norm < abs(rgb). */
{
#define sgn(a) ((a) >= 0 ? 1.0 : -1.0)
double rgbNorm [3], mantNorm;
int i, exp;
RGBEContrib rgbeCont = {0, 0, 0, 0};
/* Check and set bin */
if (bin > RGBE_CONTRIB_BINMAX) {
fprintf(stderr, "WARNING - RGBEContrib bin overflow\n");
return rgbeCont;
}
else
rgbeCont.bin = bin;
/* Normalise RGB. */
for (i = 0; i < 3; i++)
if (norm [i] == 0) {
fprintf(
stderr, "WARNING - Zero RGBEContrib normalisation\n"
);
return rgbeCont;
}
else {
rgbNorm [i] = rgb [i] / norm [i];
if (fabs(fabs(rgbNorm [i]) - 1) < RGBE_CONTRIB_MIN)
/* Normalisation to 1. This must be avoided by offsetting by
minimum value, otherwise frexp() returns a positive exponent.
Normalisation to >1 is caught by overflow check below */
rgbNorm [i] = sgn(rgbNorm [i]) * (1 - RGBE_CONTRIB_MIN);
}
/* Mantissa normalisation factor (= RGB maximum) */
mantNorm = max(
max(fabs(rgbNorm [0]), fabs(rgbNorm [1])), fabs(rgbNorm [2])
);
/* Check RGB range */
if (mantNorm < RGBE_CONTRIB_MIN) {
/* Below encoding range, clamp to zero */
fprintf(stderr, "WARNING - RGBEContrib underflow, clamping to 0\n");
return rgbeCont;
}
if (mantNorm > RGBE_CONTRIB_MAX) {
/* Exceeds encoding range, clamp to maximum */
fprintf(
stderr, "WARNING - RGBEContrib overflow, clamping to +/-%.1f\n",
RGBE_CONTRIB_MAX
);
rgbeCont.red = (1 + sgn(rgbNorm [0])) * RGBE_CONTRIB_MANTMAX + 0.5;
rgbeCont.grn = (1 + sgn(rgbNorm [1])) * RGBE_CONTRIB_MANTMAX + 0.5;
rgbeCont.blu = (1 + sgn(rgbNorm [2])) * RGBE_CONTRIB_MANTMAX + 0.5;
return rgbeCont;
}
/* Get normalised mantissa and exponent to base 2 */
mantNorm = frexp(mantNorm, &exp) * RGBE_CONTRIB_MANTMAX / mantNorm;
/* Set RGB integer mantissas, adding excess (zero offset) and rounding
up or down, depending on sign. Clamp to RGBE_CONTRIB_MANTOVF to
prevent rounding beyond overflow into negative! */
rgbeCont.red = min(
(int)(rgbNorm [0] * mantNorm + RGBE_CONTRIB_MANTMAX + 0.5),
RGBE_CONTRIB_MANTOVF
);
rgbeCont.grn = min(
(int)(rgbNorm [1] * mantNorm + RGBE_CONTRIB_MANTMAX + 0.5),
RGBE_CONTRIB_MANTOVF
);
rgbeCont.blu = min(
(int)(rgbNorm [2] * mantNorm + RGBE_CONTRIB_MANTMAX + 0.5),
RGBE_CONTRIB_MANTOVF
);
/* Drop the exponent sign, since implicitly negative */
rgbeCont.exp = abs(exp);
return rgbeCont;
#undef sgn
}
unsigned RGBEContribDec (
RGBEContrib rgbeCont, double norm [3], double rgb [3]
)
/* Decode RGBE contribution, returning signed floating point contributions
* in array rgb after denormalisation against norm, and the associated
* bin as return value. */
{
#define sgn(a) ((a) >= RGBE_CONTRIB_MANTMAX ? 0.1 : -0.1)
/* Mantissa normalisation factor; exponent implicitly negative */
const double mantNorm = ldexp(1.0, -rgbeCont.exp) / RGBE_CONTRIB_MANTMAX;
/* Decode and set RGB components, subtracting excess (zero offset) and
* adding some jitter to break up quantisation artefacts.
* Clamp to -RGBE_CONTRIB_MAX (-1) in case jitter underflows.
* Denormalise to original range with user-specified factor. */
rgb [0] = norm [0] * max(
mantNorm * (
rgbeCont.red - RGBE_CONTRIB_MANTMAX + drand48() * sgn(rgbeCont.red)
),
-RGBE_CONTRIB_MAX
);
rgb [1] = norm [1] * max(
mantNorm * (
rgbeCont.grn - RGBE_CONTRIB_MANTMAX + drand48() * sgn(rgbeCont.grn)
),
-RGBE_CONTRIB_MAX
);
rgb [2] = norm [2] * max(
mantNorm * (
rgbeCont.blu - RGBE_CONTRIB_MANTMAX + drand48() * sgn(rgbeCont.blu)
),
-RGBE_CONTRIB_MAX
);
/* Return bin */
return rgbeCont.bin;
#undef sgn
}
#ifdef _RGBE_CONTRIB_TEST
/* Build standalone to run unit test */
#define RGBMAX 1
int main (int argc, char **argv)
{
double rgb [3] = {0, 0, 0}, rgbDec [3], rgbNorm [3] = {0, 0, 0},
relDiff, rmse, rmseSum;
const short rndState0 [3] = {11, 23, 31};
short encFlag = 0, rndState [3];
unsigned numTests, i, j, k, bin, binDec;
RGBEContrib rgbeCont;
/* Test range checks */
puts("Check bin overflow");
RGBEContribEnc(rgb, rgbNorm, RGBE_CONTRIB_BINMAX << 1);
puts("\nCheck zero normalisation");
RGBEContribEnc(rgb, rgbNorm, 0);
puts("\nCheck underflow");
rgbNorm [0] = rgbNorm [1] = rgbNorm [2] = 1;
rgb [0] = RGBE_CONTRIB_MIN / 2;
RGBEContribEnc(rgb, rgbNorm, 0);
puts("\nCheck overflow (positive)");
rgb [0] = 2;
RGBEContribEnc(rgb, rgbNorm, 0);
puts("\nCheck overflow (negative)");
rgb [0] = -2;
RGBEContribEnc(rgb, rgbNorm, 0);
putchar('\n');
if (argc > 1 && (numTests = atoi(argv [1])) > 0) {
rgbNorm [0] = rgbNorm [1] = rgbNorm [2] = 0;
/* Iterate twice with same random RGB tuples; gather normalisation
factors in 1st pass, then do actual encode/decode tests using
found normalisation 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++) {
bin = (unsigned)(RGBE_CONTRIB_BINMAX * erand48(rndState));
#if 1
for (i = 0; i < 3; i++) {
if (i > 0)
/* Correlate RGB */
rgb [i] = rgb [0] * (0.1 + 0.9 * erand48(rndState));
else
rgb [i] = RGBMAX * (2 * erand48(rndState) - 1);
}
#else
#if 1
rgb [0] = rgb [1] = rgb [2] = RGBMAX * (
2 * erand48(rndState) - 1
);
#else
rgb [0] = 0.9772;
rgb [1] = 0.4369;
rgb [2] = 0.1869;
#endif
#endif
if (encFlag) {
/* Encode */
/*
rgb [0] = rgbNorm [0];
rgb [1] = 0.0077;
rgb [2] = 0.0058;
*/
rgbeCont = RGBEContribEnc(rgb, rgbNorm, bin);
printf(
"% 7.4f\t% 7.4f\t% 7.4f\t%4d\t-->\t"
"%02d:%02d:%02d:%02d\t-->\t",
rgb [0], rgb [1], rgb [2], bin,
rgbeCont.red, rgbeCont.grn, rgbeCont.blu, rgbeCont.exp
);
/* Decode */
binDec = RGBEContribDec(rgbeCont, rgbNorm, rgbDec);
for (j = rmse = 0; j < 3; j++) {
relDiff = 1 - rgbDec [j] / rgb [j];
rmse += relDiff * relDiff;
}
rmseSum += rmse = sqrt(rmse / 3);
printf(
"% 7.4f\t% 7.4f\t% 7.4f\t%4d\tRMSE = %.1f%%\n",
rgbDec [0], rgbDec [1], rgbDec [2], binDec, 100 * rmse
);
}
else {
/* Update RGB normalisation factor */
for (i = 0; i < 3; i++)
if (fabs(rgb [i]) > rgbNorm [i])
rgbNorm [i] = fabs(rgb [i]);
}
}
/* Dump results */
if (!encFlag)
printf(
"RGB Normalisation = (%.3f, %.3f, %.3f)\n\n",
rgbNorm [0], rgbNorm [1], rgbNorm [2]
);
else
printf("\nAvg RMSE = %.1f%%\n", 100 * rmseSum / numTests);
} while (encFlag = !encFlag);
}
return 0;
}
#endif
c4science · Help