diff --git a/mrgbe.c b/mrgbe.c index 4542388..7eba4e1 100644 --- a/mrgbe.c +++ b/mrgbe.c @@ -1,423 +1,397 @@ -/* +/* ========================================================================= 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 + 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)) #define ABS(a) ((a) >= 0 ? (a) : -(a)) -mRGBERange *mRGBEinit (mRGBERange *range, +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++) { - if (rgbMin [i] >= rgbMax [i]) { - fprintf(stderr, "ERROR - invalid mRGBE range\n"); - return NULL; - } - range -> min [i] = orgMin = fabs(rgbMin [i]); range -> max [i] = fabs(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" + + 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)) + #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. + 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. + 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)) * + 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}, + 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}; + unsigned numTests, i, j, k, data, dataDec; + const short rndState0 [3] = {11, 23, 31}; short encFlag = 0, rndState [3]; mRGBERange mrgbeRange; mRGBE mrgbe; - #if 0 - rgb [0] = -1.481; - rgb [1] = -0.437; - rgb [2] = -1.392; - - rgbMin [0] = 1.003; - rgbMin [1] = 0.436; - rgbMin [2] = 0.299; - - rgbMax [0] = 14.863; - rgbMax [1] = 13.505; - rgbMax [2] = 13.359; - - mRGBEinit(&mrgbeRange, rgbMin, rgbMax); - mrgbe = mRGBEencode(rgb, &mrgbeRange, 0); - 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], 0, - 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; - } - 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 - ); - - return 0; - #endif - 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], + 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 underflow"); - rgb [0] = rgb [1] = rgb [2] = mRGBE_MIN / 2; - mRGBEencode(rgb, &mrgbeRange, 0); - + 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 empty range"); - rgbMax [0] = rgbMax [1] = rgbMax [2] = 0; - mRGBEinit(&mrgbeRange, rgbMin, rgbMax); + 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 + /* 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)); + data = (unsigned)(mRGBE_DATAMAX * erand48(rndState)); for (i = 0; i < 3; i++) { - - rgb [i] = (erand48(rndState) > 0.5 ? 1 : -1) * + + rgb [i] = (erand48(rndState) > 0.5 ? 1 : -1) * ((RGBMAX - RGBMIN) * erand48(rndState) + RGBMIN); #if 1 /* Correlate RGB */ - if (i > 0) + 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); + 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, + "%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 + 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], + 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 + /* 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, + 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 \n", argv [0]); return -1; } } #endif