diff --git a/mrgbe.c b/mrgbe.c new file mode 100644 index 0000000..ccbd244 --- /dev/null +++ b/mrgbe.c @@ -0,0 +1,269 @@ +/* + ========================================================================= + 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 +#include +#include + + + +#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 \n", argv [0]); + return -1; + } + } +#endif diff --git a/mrgbe.h b/mrgbe.h new file mode 100644 index 0000000..3d8cbc3 --- /dev/null +++ b/mrgbe.h @@ -0,0 +1,76 @@ +/* + ========================================================================= + mRGBE (miniRGBE): reduced RGBE encoding of normalised RGB floats plus + associated integer payload data. By default, this encoding consists of: + + - three signed 5-bit mantissas for each RGB colour channel, + - a common 5-bit exponent, + - and an associated 12-bit integer datum. + + However, this allocation can be modified within a 32-bit envelope if more + precision is required by modifying mRGBE_MANTBITS, mRGBE_EXPBITS and + mRGBE_DATABITS. + + The encoded float RGB tuples are assumed to lie in the range [-1, 1], + with values in the interval [-mRGBE_MIN, mRGBE_MIN] clamped to 0, since + these cannot be resolved. The encoded values are therefore normalised + against a given normalisation factor (e.g. absolute maximum) in order to + maximise the useful range. + + 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$ +*/ + + +#ifndef _mRGBE_H + #define _mRGBE_H + + #include + + + /* Mantissa / exponent / payload data bits and their encoding ranges. + NOTE: binary shifts can overflow if the number of bits is increased; + in this case shifts must either be promoted to 64 bits (possibly at + the expense of portability), or replaced by much slower pow(2, ..) */ + #define mRGBE_MANTBITS 5 + #define mRGBE_EXPBITS 5 + #define mRGBE_DATABITS 12 + + /* Signed maximum for mantissa */ + #define mRGBE_MANTMAX (1 << (mRGBE_MANTBITS - 1)) + /* Absolute overflow limit for mantissa (sign flips if exceeded) */ + #define mRGBE_MANTOVF ((mRGBE_MANTMAX << 1) - 1) + #define mRGBE_MIN (1.0 / (1L << (1 << mRGBE_EXPBITS))) + #define mRGBE_MAX 1.0 + #define mRGBE_DATAMAX ((1U << mRGBE_DATABITS) - 1) + + + typedef union { + struct { + /* (mini)RGBE representation consisting of a 4-bit mantissa + sign + * per RGB, a common 5-bit exponent, and a 12-bit payload data */ + unsigned red: mRGBE_MANTBITS, + grn: mRGBE_MANTBITS, + blu: mRGBE_MANTBITS, + exp: mRGBE_EXPBITS, + dat: mRGBE_DATABITS; + }; + + uint32_t all; + } mRGBE; + + + 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). */ + + 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. */ +#endif