Page MenuHomec4science
Files F62645369

mrgbe.c
No OneTemporary
Actions

File Metadata

Created
Tue, May 14, 13:57

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))
mRGBE mRGBEencode (double rgb [3], double norm [3], unsigned data)
/* Normalise signed float RGB tuple by abs(norm), and encode to mRGBE
* with associated data. Issues warning if abs(norm) < abs(rgb). */
{
#define sgn(a) ((a) >= 0 ? 1.0 : -1.0)
double rgbNorm [3], mantNorm;
int i, exp;
mRGBE mrgbe = {0, 0, 0, 0};
/* Check and set payload data */
if (data > mRGBE_DATAMAX) {
fprintf(stderr, "WARNING - mRGBE data overflow\n");
return mrgbe;
}
else
mrgbe.dat = data;
/* Normalise RGB. */
for (i = 0; i < 3; i++)
if (norm [i] == 0) {
fprintf(stderr, "WARNING - Zero mRGBE normalisation\n");
return mrgbe;
}
else {
rgbNorm [i] = rgb [i] / norm [i];
if (fabs(fabs(rgbNorm [i]) - 1) < mRGBE_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 - mRGBE_MIN);
}
/* Mantissa normalisation factor (= RGB maximum) */
mantNorm = max(
max(fabs(rgbNorm [0]), fabs(rgbNorm [1])),
fabs(rgbNorm [2])
);
/* Check RGB range */
if (mantNorm < mRGBE_MIN) {
/* Below encoding range, clamp to zero */
fprintf(stderr, "WARNING - mRGBE underflow, clamping to 0\n");
return mrgbe;
}
if (mantNorm > mRGBE_MAX) {
/* Exceeds encoding range, clamp to maximum */
fprintf(stderr,
"WARNING - mRGBE overflow, clamping to +/-%.1f\n", mRGBE_MAX
);
mrgbe.red = (1 + sgn(rgbNorm [0])) * mRGBE_MANTMAX + 0.5;
mrgbe.grn = (1 + sgn(rgbNorm [1])) * mRGBE_MANTMAX + 0.5;
mrgbe.blu = (1 + sgn(rgbNorm [2])) * mRGBE_MANTMAX + 0.5;
return mrgbe;
}
/* 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. Clamp to mRGBE_MANTOVF to prevent
rounding beyond overflow into negative! */
mrgbe.red = min(
(int)(rgbNorm [0] * mantNorm + mRGBE_MANTMAX + 0.5), mRGBE_MANTOVF
);
mrgbe.grn = min(
(int)(rgbNorm [1] * mantNorm + mRGBE_MANTMAX + 0.5), mRGBE_MANTOVF
);
mrgbe.blu = min(
(int)(rgbNorm [2] * mantNorm + mRGBE_MANTMAX + 0.5), mRGBE_MANTOVF
);
/* Drop the exponent sign, since implicitly negative */
mrgbe.exp = abs(exp);
return mrgbe;
#undef sgn
}
unsigned mRGBEdecode (mRGBE mrgbe, double norm [3], double rgb [3])
/* Decode mRGBE, returning signed float RGB tuple in array rgb after
* denormalisation against abs(norm). The associated payload data is
* decoded as return value. */
{
#define sgn(a) ((a) >= mRGBE_MANTMAX ? 0.1 : -0.1)
/* Mantissa normalisation factor; exponent implicitly negative */
const double mantNorm = ldexp(1.0, -mrgbe.exp) / mRGBE_MANTMAX;
/* Decode and set RGB components, subtracting excess (zero offset) and
* adding some jitter to break up quantisation artefacts. Clamp to
* -mRGBE_MAX (=-1) in case jitter underflows. Denormalise to original
* range with user-specified factor. */
rgb [0] = norm [0] * max(
mantNorm * (mrgbe.red - mRGBE_MANTMAX + drand48() * sgn(mrgbe.red)),
-mRGBE_MAX
);
rgb [1] = norm [1] * max(
mantNorm * (mrgbe.grn - mRGBE_MANTMAX + drand48() * sgn(mrgbe.grn)),
-mRGBE_MAX
);
rgb [2] = norm [2] * max(
mantNorm * (mrgbe.blu - mRGBE_MANTMAX + drand48() * sgn(mrgbe.blu)),
-mRGBE_MAX
);
/* Return payload data */
return mrgbe.dat;
#undef sgn
}
#ifdef mRGBE_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, data, dataDec;
mRGBE mrgbe;
if (argc > 1 && (numTests = atoi(argv [1])) > 0) {
/* Test range checks */
puts("Check data overflow");
mRGBEencode(rgb, rgbNorm, mRGBE_DATAMAX << 1);
puts("\nCheck zero normalisation");
mRGBEencode(rgb, rgbNorm, 0);
puts("\nCheck underflow");
rgbNorm [0] = rgbNorm [1] = rgbNorm [2] = 1;
rgb [0] = mRGBE_MIN / 2;
mRGBEencode(rgb, rgbNorm, 0);
puts("\nCheck overflow (positive)");
rgb [0] = 2;
mRGBEencode(rgb, rgbNorm, 0);
puts("\nCheck overflow (negative)");
rgb [0] = -2;
mRGBEencode(rgb, rgbNorm, 0);
putchar('\n');
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++) {
data = (unsigned)(mRGBE_DATAMAX * 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;
*/
mrgbe = mRGBEencode(rgb, rgbNorm, data);
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], data,
mrgbe.red, mrgbe.grn, mrgbe.blu, mrgbe.exp
);
/* Decode */
dataDec = mRGBEdecode(mrgbe, 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],
dataDec, 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;
}
else {
printf("Usage: %s <numTests>\n", argv [0]);
return -1;
}
}
#endif

Event Timeline

c4science · Help