Add rygorous's float->srgb8 conversion routine to the stdlib.

Issue #230

docs/ispc.rst

@@ -147,6 +147,7 @@ Contents:
 
     * `Converting Between Array-of-Structures and Structure-of-Arrays Layout`_
     * `Conversions To and From Half-Precision Floats`_
+    * `Converting to sRGB8`_
 
   + `Systems Programming Support`_
 
@@ -3691,6 +3692,22 @@ precise.
     uniform int16 float_to_half_fast(uniform float f)
 
 
+Converting to sRGB8
+-------------------
+
+The sRGB color space is used in many applications in graphics and imaging;
+see the `Wikipedia page on sRGB`_ for more information.  The ``ispc``
+standard library provides two functions for converting floating-point color
+values to 8-bit values in the sRGB space.
+
+.. _Wikipedia page on sRGB: http://en.wikipedia.org/wiki/SRGB
+
+::
+
+    int float_to_srgb8(float v)
+    uniform int float_to_srgb8(uniform float v)
+
+
 Systems Programming Support
 ---------------------------
 

stdlib.ispc

@@ -3829,6 +3829,133 @@ static inline int16 float_to_half_fast(float f) {
     }
 }
 
+///////////////////////////////////////////////////////////////////////////
+// float -> srgb8
+
+// https://gist.github.com/2246678, from Fabian "rygorous" Giesen.
+//
+// The basic ideas are still the same, only this time, we squeeze
+// everything into the table, even the linear part of the range; since we
+// are approximating the function as piecewise linear anyway, this is
+// fairly easy.
+//
+// In the exact version of the conversion, any value that produces an
+// output float less than 0.5 will be rounded to an integer of
+// zero. Inverting the linear part of the transform, we get:
+//
+//   log2(0.5 / (255 * 12.92)) =~ -12.686
+//
+// which in turn means that any value smaller than about 2^(-12.687) will
+// return 0.  What this means is that we can adapt the clamping code to
+// just clamp to [2^(-13), 1-eps] and we're covered. This means our table
+// needs to cover a range of 13 different exponents from -13 to -1.
+//
+// The table lookup, storage and interpolation works exactly the same way
+// as in the code above.
+//
+// Max error for the whole function (integer-rounded result minus "exact"
+// value, as computed in floats using the official formula): 0.544403 at
+// 0x3e9f8000
+
+__declspec(safe)
+static inline int
+float_to_srgb8(float in)
+{
+    static const uniform unsigned int table[104] = {
+        0x0073000d, 0x007a000d, 0x0080000d, 0x0087000d,
+        0x008d000d, 0x0094000d, 0x009a000d, 0x00a1000d,
+        0x00a7001a, 0x00b4001a, 0x00c1001a, 0x00ce001a,
+        0x00da001a, 0x00e7001a, 0x00f4001a, 0x0101001a,
+        0x010e0033, 0x01280033, 0x01410033, 0x015b0033,
+        0x01750033, 0x018f0033, 0x01a80033, 0x01c20033,
+        0x01dc0067, 0x020f0067, 0x02430067, 0x02760067,
+        0x02aa0067, 0x02dd0067, 0x03110067, 0x03440067,
+        0x037800ce, 0x03df00ce, 0x044600ce, 0x04ad00ce,
+        0x051400ce, 0x057b00c5, 0x05dd00bc, 0x063b00b5,
+        0x06970158, 0x07420142, 0x07e30130, 0x087b0120,
+        0x090b0112, 0x09940106, 0x0a1700fc, 0x0a9500f2,
+        0x0b0f01cb, 0x0bf401ae, 0x0ccb0195, 0x0d950180,
+        0x0e56016e, 0x0f0d015e, 0x0fbc0150, 0x10630143,
+        0x11070264, 0x1238023e, 0x1357021d, 0x14660201,
+        0x156601e9, 0x165a01d3, 0x174401c0, 0x182401af,
+        0x18fe0331, 0x1a9602fe, 0x1c1502d2, 0x1d7e02ad,
+        0x1ed4028d, 0x201a0270, 0x21520256, 0x227d0240,
+        0x239f0443, 0x25c003fe, 0x27bf03c4, 0x29a10392,
+        0x2b6a0367, 0x2d1d0341, 0x2ebe031f, 0x304d0300,
+        0x31d105b0, 0x34a80555, 0x37520507, 0x39d504c5,
+        0x3c37048b, 0x3e7c0458, 0x40a8042a, 0x42bd0401,
+        0x44c20798, 0x488e071e, 0x4c1c06b6, 0x4f76065d,
+        0x52a50610, 0x55ac05cc, 0x5892058f, 0x5b590559,
+        0x5e0c0a23, 0x631c0980, 0x67db08f6, 0x6c55087f,
+        0x70940818, 0x74a007bd, 0x787d076c, 0x7c330723,
+    };
+
+    static const uniform unsigned int almost_one = 0x3f7fffff;
+    
+    // Clamp to [2^(-13), 1-eps]; these two values map to 0 and 1, respectively.
+    in = max(in, 0.0f);
+    in = min(in, floatbits(almost_one));
+
+    // Do the table lookup and unpack bias, scale
+    unsigned int tab = table[(intbits(in) - 0x39000000u) >> 20];
+    unsigned int bias = (tab >> 16) << 9;
+    unsigned int scale = tab & 0xffff;
+
+    // Grab next-highest mantissa bits and perform linear interpolation
+    unsigned int t = (intbits(in) >> 12) & 0xff;
+    return (bias + scale*t) >> 16;
+}
+
+
+__declspec(safe)
+static inline uniform int
+float_to_srgb8(uniform float in)
+{
+    static const uniform unsigned int table[104] = {
+        0x0073000d, 0x007a000d, 0x0080000d, 0x0087000d,
+        0x008d000d, 0x0094000d, 0x009a000d, 0x00a1000d,
+        0x00a7001a, 0x00b4001a, 0x00c1001a, 0x00ce001a,
+        0x00da001a, 0x00e7001a, 0x00f4001a, 0x0101001a,
+        0x010e0033, 0x01280033, 0x01410033, 0x015b0033,
+        0x01750033, 0x018f0033, 0x01a80033, 0x01c20033,
+        0x01dc0067, 0x020f0067, 0x02430067, 0x02760067,
+        0x02aa0067, 0x02dd0067, 0x03110067, 0x03440067,
+        0x037800ce, 0x03df00ce, 0x044600ce, 0x04ad00ce,
+        0x051400ce, 0x057b00c5, 0x05dd00bc, 0x063b00b5,
+        0x06970158, 0x07420142, 0x07e30130, 0x087b0120,
+        0x090b0112, 0x09940106, 0x0a1700fc, 0x0a9500f2,
+        0x0b0f01cb, 0x0bf401ae, 0x0ccb0195, 0x0d950180,
+        0x0e56016e, 0x0f0d015e, 0x0fbc0150, 0x10630143,
+        0x11070264, 0x1238023e, 0x1357021d, 0x14660201,
+        0x156601e9, 0x165a01d3, 0x174401c0, 0x182401af,
+        0x18fe0331, 0x1a9602fe, 0x1c1502d2, 0x1d7e02ad,
+        0x1ed4028d, 0x201a0270, 0x21520256, 0x227d0240,
+        0x239f0443, 0x25c003fe, 0x27bf03c4, 0x29a10392,
+        0x2b6a0367, 0x2d1d0341, 0x2ebe031f, 0x304d0300,
+        0x31d105b0, 0x34a80555, 0x37520507, 0x39d504c5,
+        0x3c37048b, 0x3e7c0458, 0x40a8042a, 0x42bd0401,
+        0x44c20798, 0x488e071e, 0x4c1c06b6, 0x4f76065d,
+        0x52a50610, 0x55ac05cc, 0x5892058f, 0x5b590559,
+        0x5e0c0a23, 0x631c0980, 0x67db08f6, 0x6c55087f,
+        0x70940818, 0x74a007bd, 0x787d076c, 0x7c330723,
+    };
+
+    static const uniform unsigned int almost_one = 0x3f7fffff;
+    
+    // Clamp to [2^(-13), 1-eps]; these two values map to 0 and 1, respectively.
+    in = max(in, 0.0f);
+    in = min(in, floatbits(almost_one));
+
+    // Do the table lookup and unpack bias, scale
+    uniform unsigned int tab = table[(intbits(in) - 0x39000000u) >> 20];
+    uniform unsigned int bias = (tab >> 16) << 9;
+    uniform unsigned int scale = tab & 0xffff;
+
+    // Grab next-highest mantissa bits and perform linear interpolation
+    uniform unsigned int t = (intbits(in) >> 12) & 0xff;
+    return (bias + scale*t) >> 16;
+}
+
 ///////////////////////////////////////////////////////////////////////////
 // RNG stuff