[babl] base: added IEEE 754-2008 half precision float
- From: Ãyvind KolÃs <ok src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [babl] base: added IEEE 754-2008 half precision float
- Date: Sat, 21 Apr 2012 23:34:47 +0000 (UTC)
commit 3b93d527598f1b4f9171e460c5be754b9ccf18de
Author: Ãyvind KolÃs <pippin gimp org>
Date: Sun Apr 22 01:31:05 2012 +0200
base: added IEEE 754-2008 half precision float
Using BSD licensed code that means applications using babl should display the
appropriate copyright header - or this code needs to be reimplemented so that
babl can retain the functionality without that additional requirement.
babl/babl-ids.h | 1 +
babl/base/Makefile.am | 1 +
babl/base/babl-base.c | 1 +
babl/base/babl-base.h | 1 +
babl/base/type-half.c | 351 +++++++++++++++++++++++++++++++++++++++++++++++++
5 files changed, 355 insertions(+), 0 deletions(-)
---
diff --git a/babl/babl-ids.h b/babl/babl-ids.h
index f0f6c61..47b5aa5 100644
--- a/babl/babl-ids.h
+++ b/babl/babl-ids.h
@@ -25,6 +25,7 @@ enum {
BABL_U8,
BABL_U16,
BABL_U32,
+ BABL_HALF,
BABL_FLOAT,
BABL_DOUBLE,
BABL_HALF_FLOAT,
diff --git a/babl/base/Makefile.am b/babl/base/Makefile.am
index 1c1da0f..74052c2 100644
--- a/babl/base/Makefile.am
+++ b/babl/base/Makefile.am
@@ -6,6 +6,7 @@ c_sources = \
babl-base.c \
formats.c \
type-float.c \
+ type-half.c \
type-u8.c \
type-u16.c \
type-u32.c \
diff --git a/babl/base/babl-base.c b/babl/base/babl-base.c
index 1271a43..0831eb1 100644
--- a/babl/base/babl-base.c
+++ b/babl/base/babl-base.c
@@ -50,6 +50,7 @@ babl_base_destroy (void)
static void
types (void)
{
+ babl_base_type_half ();
babl_base_type_float ();
babl_base_type_u8 ();
babl_base_type_u16 ();
diff --git a/babl/base/babl-base.h b/babl/base/babl-base.h
index bed3c64..f3bf6cd 100644
--- a/babl/base/babl-base.h
+++ b/babl/base/babl-base.h
@@ -23,6 +23,7 @@ void babl_base_init (void);
void babl_base_destroy (void);
void babl_formats_init (void);
+void babl_base_type_half (void);
void babl_base_type_float (void);
void babl_base_type_u8 (void);
void babl_base_type_u16 (void);
diff --git a/babl/base/type-half.c b/babl/base/type-half.c
new file mode 100644
index 0000000..987b88f
--- /dev/null
+++ b/babl/base/type-half.c
@@ -0,0 +1,351 @@
+/* babl - dynamically extendable universal pixel conversion library.
+ * Copyright (C) 2012, Ãyvind KolÃs.
+ *
+ * MATLAB (R) is a trademark of The Mathworks (R) Corporation
+ *
+ * Function: halfprecision
+ * Filename: halfprecision.c
+ * Programmer: James Tursa
+ * Version: 1.0
+ * Date: March 3, 2009
+ * Copyright: (c) 2009 by James Tursa, All Rights Reserved
+ *
+ * This code uses the BSD License:
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are
+ * met:
+ *
+ * * Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in
+ * the documentation and/or other materials provided with the distribution
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ *
+ * halfprecision converts the input argument to/from a half precision floating
+ * point bit pattern corresponding to IEEE 754r. The bit pattern is stored in a
+ * uint16 class variable. Please note that halfprecision is *not* a class. That
+ * is, you cannot do any arithmetic with the half precision bit patterns.
+ * halfprecision is simply a function that converts the IEEE 754r half precision
+ * bit pattern to/from other numeric MATLAB variables. You can, however, take
+ * the half precision bit patterns, convert them to single or double, do the
+ * operation, and then convert the result back manually.
+ *
+ * 1 bit sign bit
+ * 5 bits exponent, biased by 15
+ * 10 bits mantissa, hidden leading bit, normalized to 1.0
+ *
+ * Special floating point bit patterns recognized and supported:
+ *
+ * All exponent bits zero:
+ * - If all mantissa bits are zero, then number is zero (possibly signed)
+ * - Otherwise, number is a denormalized bit pattern
+ *
+ * All exponent bits set to 1:
+ * - If all mantissa bits are zero, then number is +Infinity or -Infinity
+ * - Otherwise, number is NaN (Not a Number)
+ */
+
+#include "config.h"
+#include <string.h>
+#include <assert.h>
+
+#include "babl.h"
+#include "babl-classes.h"
+#include "babl-ids.h"
+#include "babl-base.h"
+
+
+#define mwSize int
+#define INT16_TYPE short
+#define UINT16_TYPE unsigned short
+#define INT32_TYPE long
+#define UINT32_TYPE unsigned long
+
+static int next = 1; /* should be 0 for big endian */
+
+//-----------------------------------------------------------------------------
+
+static void doubles2halfp(void *target, void *source, mwSize numel)
+{
+ UINT16_TYPE *hp = (UINT16_TYPE *) target; // Type pun output as an unsigned 16-bit int
+ UINT32_TYPE *xp = (UINT32_TYPE *) source; // Type pun input as an unsigned 32-bit int
+ UINT16_TYPE hs, he, hm;
+ UINT32_TYPE x, xs, xe, xm;
+ int hes;
+
+ xp += next; // Little Endian adjustment if necessary
+
+ if( source == NULL || target == NULL ) { // Nothing to convert (e.g., imag part of pure real)
+ return;
+ }
+ while( numel-- ) {
+ x = *xp++; xp++; // The extra xp++ is to skip over the remaining 32 bits of the mantissa
+ if( (x & 0x7FFFFFFFu) == 0 ) { // Signed zero
+ *hp++ = (UINT16_TYPE) (x >> 16); // Return the signed zero
+ } else { // Not zero
+ xs = x & 0x80000000u; // Pick off sign bit
+ xe = x & 0x7FF00000u; // Pick off exponent bits
+ xm = x & 0x000FFFFFu; // Pick off mantissa bits
+ if( xe == 0 ) { // Denormal will underflow, return a signed zero
+ *hp++ = (UINT16_TYPE) (xs >> 16);
+ } else if( xe == 0x7FF00000u ) { // Inf or NaN (all the exponent bits are set)
+ if( xm == 0 ) { // If mantissa is zero ...
+ *hp++ = (UINT16_TYPE) ((xs >> 16) | 0x7C00u); // Signed Inf
+ } else {
+ *hp++ = (UINT16_TYPE) 0xFE00u; // NaN, only 1st mantissa bit set
+ }
+ } else { // Normalized number
+ hs = (UINT16_TYPE) (xs >> 16); // Sign bit
+ hes = ((int)(xe >> 20)) - 1023 + 15; // Exponent unbias the double, then bias the halfp
+ if( hes >= 0x1F ) { // Overflow
+ *hp++ = (UINT16_TYPE) ((xs >> 16) | 0x7C00u); // Signed Inf
+ } else if( hes <= 0 ) { // Underflow
+ if( (10 - hes) > 21 ) { // Mantissa shifted all the way off & no rounding possibility
+ hm = (UINT16_TYPE) 0u; // Set mantissa to zero
+ } else {
+ xm |= 0x00100000u; // Add the hidden leading bit
+ hm = (UINT16_TYPE) (xm >> (11 - hes)); // Mantissa
+ if( (xm >> (10 - hes)) & 0x00000001u ) // Check for rounding
+ hm += (UINT16_TYPE) 1u; // Round, might overflow into exp bit, but this is OK
+ }
+ *hp++ = (hs | hm); // Combine sign bit and mantissa bits, biased exponent is zero
+ } else {
+ he = (UINT16_TYPE) (hes << 10); // Exponent
+ hm = (UINT16_TYPE) (xm >> 10); // Mantissa
+ if( xm & 0x00000200u ) // Check for rounding
+ *hp++ = (hs | he | hm) + (UINT16_TYPE) 1u; // Round, might overflow to inf, this is OK
+ else
+ *hp++ = (hs | he | hm); // No rounding
+ }
+ }
+ }
+ }
+}
+
+#if 0
+//-----------------------------------------------------------------------------
+
+static void halfp2singles(void *target, void *source, mwSize numel)
+{
+ UINT16_TYPE *hp = (UINT16_TYPE *) source; // Type pun input as an unsigned 16-bit int
+ UINT32_TYPE *xp = (UINT32_TYPE *) target; // Type pun output as an unsigned 32-bit int
+ UINT16_TYPE h, hs, he, hm;
+ UINT32_TYPE xs, xe, xm;
+ INT32_TYPE xes;
+ int e;
+
+ if( source == NULL || target == NULL ) // Nothing to convert (e.g., imag part of pure real)
+ return;
+ while( numel-- ) {
+ h = *hp++;
+ if( (h & 0x7FFFu) == 0 ) { // Signed zero
+ *xp++ = ((UINT32_TYPE) h) << 16; // Return the signed zero
+ } else { // Not zero
+ hs = h & 0x8000u; // Pick off sign bit
+ he = h & 0x7C00u; // Pick off exponent bits
+ hm = h & 0x03FFu; // Pick off mantissa bits
+ if( he == 0 ) { // Denormal will convert to normalized
+ e = -1; // The following loop figures out how much extra to adjust the exponent
+ do {
+ e++;
+ hm <<= 1;
+ } while( (hm & 0x0400u) == 0 ); // Shift until leading bit overflows into exponent bit
+ xs = ((UINT32_TYPE) hs) << 16; // Sign bit
+ xes = ((INT32_TYPE) (he >> 10)) - 15 + 127 - e; // Exponent unbias the halfp, then bias the single
+ xe = (UINT32_TYPE) (xes << 23); // Exponent
+ xm = ((UINT32_TYPE) (hm & 0x03FFu)) << 13; // Mantissa
+ *xp++ = (xs | xe | xm); // Combine sign bit, exponent bits, and mantissa bits
+ } else if( he == 0x7C00u ) { // Inf or NaN (all the exponent bits are set)
+ if( hm == 0 ) { // If mantissa is zero ...
+ *xp++ = (((UINT32_TYPE) hs) << 16) | ((UINT32_TYPE) 0x7F800000u); // Signed Inf
+ } else {
+ *xp++ = (UINT32_TYPE) 0xFFC00000u; // NaN, only 1st mantissa bit set
+ }
+ } else { // Normalized number
+ xs = ((UINT32_TYPE) hs) << 16; // Sign bit
+ xes = ((INT32_TYPE) (he >> 10)) - 15 + 127; // Exponent unbias the halfp, then bias the single
+ xe = (UINT32_TYPE) (xes << 23); // Exponent
+ xm = ((UINT32_TYPE) hm) << 13; // Mantissa
+ *xp++ = (xs | xe | xm); // Combine sign bit, exponent bits, and mantissa bits
+ }
+ }
+ }
+}
+
+static void singles2halfp(void *target, void *source, mwSize numel)
+{
+ UINT16_TYPE *hp = (UINT16_TYPE *) target; // Type pun output as an unsigned 16-bit int
+ UINT32_TYPE *xp = (UINT32_TYPE *) source; // Type pun input as an unsigned 32-bit int
+ UINT16_TYPE hs, he, hm;
+ UINT32_TYPE x, xs, xe, xm;
+ int hes;
+
+ if( source == NULL || target == NULL ) { // Nothing to convert (e.g., imag part of pure real)
+ return;
+ }
+ while( numel-- ) {
+ x = *xp++;
+ if( (x & 0x7FFFFFFFu) == 0 ) { // Signed zero
+ *hp++ = (UINT16_TYPE) (x >> 16); // Return the signed zero
+ } else { // Not zero
+ xs = x & 0x80000000u; // Pick off sign bit
+ xe = x & 0x7F800000u; // Pick off exponent bits
+ xm = x & 0x007FFFFFu; // Pick off mantissa bits
+ if( xe == 0 ) { // Denormal will underflow, return a signed zero
+ *hp++ = (UINT16_TYPE) (xs >> 16);
+ } else if( xe == 0x7F800000u ) { // Inf or NaN (all the exponent bits are set)
+ if( xm == 0 ) { // If mantissa is zero ...
+ *hp++ = (UINT16_TYPE) ((xs >> 16) | 0x7C00u); // Signed Inf
+ } else {
+ *hp++ = (UINT16_TYPE) 0xFE00u; // NaN, only 1st mantissa bit set
+ }
+ } else { // Normalized number
+ hs = (UINT16_TYPE) (xs >> 16); // Sign bit
+ hes = ((int)(xe >> 23)) - 127 + 15; // Exponent unbias the single, then bias the halfp
+ if( hes >= 0x1F ) { // Overflow
+ *hp++ = (UINT16_TYPE) ((xs >> 16) | 0x7C00u); // Signed Inf
+ } else if( hes <= 0 ) { // Underflow
+ if( (14 - hes) > 24 ) { // Mantissa shifted all the way off & no rounding possibility
+ hm = (UINT16_TYPE) 0u; // Set mantissa to zero
+ } else {
+ xm |= 0x00800000u; // Add the hidden leading bit
+ hm = (UINT16_TYPE) (xm >> (14 - hes)); // Mantissa
+ if( (xm >> (13 - hes)) & 0x00000001u ) // Check for rounding
+ hm += (UINT16_TYPE) 1u; // Round, might overflow into exp bit, but this is OK
+ }
+ *hp++ = (hs | hm); // Combine sign bit and mantissa bits, biased exponent is zero
+ } else {
+ he = (UINT16_TYPE) (hes << 10); // Exponent
+ hm = (UINT16_TYPE) (xm >> 13); // Mantissa
+ if( xm & 0x00001000u ) // Check for rounding
+ *hp++ = (hs | he | hm) + (UINT16_TYPE) 1u; // Round, might overflow to inf, this is OK
+ else
+ *hp++ = (hs | he | hm); // No rounding
+ }
+ }
+ }
+ }
+}
+#endif
+
+
+//-----------------------------------------------------------------------------
+
+static void halfp2doubles(void *target, void *source, mwSize numel)
+{
+ UINT16_TYPE *hp = (UINT16_TYPE *) source; // Type pun input as an unsigned 16-bit int
+ UINT32_TYPE *xp = (UINT32_TYPE *) target; // Type pun output as an unsigned 32-bit int
+ UINT16_TYPE h, hs, he, hm;
+ UINT32_TYPE xs, xe, xm;
+ INT32_TYPE xes;
+ int e;
+
+ xp += next; // Little Endian adjustment if necessary
+
+ if( source == NULL || target == NULL ) // Nothing to convert (e.g., imag part of pure real)
+ return;
+ while( numel-- ) {
+ h = *hp++;
+ if( (h & 0x7FFFu) == 0 ) { // Signed zero
+ *xp++ = ((UINT32_TYPE) h) << 16; // Return the signed zero
+ } else { // Not zero
+ hs = h & 0x8000u; // Pick off sign bit
+ he = h & 0x7C00u; // Pick off exponent bits
+ hm = h & 0x03FFu; // Pick off mantissa bits
+ if( he == 0 ) { // Denormal will convert to normalized
+ e = -1; // The following loop figures out how much extra to adjust the exponent
+ do {
+ e++;
+ hm <<= 1;
+ } while( (hm & 0x0400u) == 0 ); // Shift until leading bit overflows into exponent bit
+ xs = ((UINT32_TYPE) hs) << 16; // Sign bit
+ xes = ((INT32_TYPE) (he >> 10)) - 15 + 1023 - e; // Exponent unbias the halfp, then bias the double
+ xe = (UINT32_TYPE) (xes << 20); // Exponent
+ xm = ((UINT32_TYPE) (hm & 0x03FFu)) << 10; // Mantissa
+ *xp++ = (xs | xe | xm); // Combine sign bit, exponent bits, and mantissa bits
+ } else if( he == 0x7C00u ) { // Inf or NaN (all the exponent bits are set)
+ if( hm == 0 ) { // If mantissa is zero ...
+ *xp++ = (((UINT32_TYPE) hs) << 16) | ((UINT32_TYPE) 0x7FF00000u); // Signed Inf
+ } else {
+ *xp++ = (UINT32_TYPE) 0xFFF80000u; // NaN, only the 1st mantissa bit set
+ }
+ } else {
+ xs = ((UINT32_TYPE) hs) << 16; // Sign bit
+ xes = ((INT32_TYPE) (he >> 10)) - 15 + 1023; // Exponent unbias the halfp, then bias the double
+ xe = (UINT32_TYPE) (xes << 20); // Exponent
+ xm = ((UINT32_TYPE) hm) << 10; // Mantissa
+ *xp++ = (xs | xe | xm); // Combine sign bit, exponent bits, and mantissa bits
+ }
+ }
+ xp++; // Skip over the remaining 32 bits of the mantissa
+ }
+}
+
+static long
+convert_double_half (char *src,
+ char *dst,
+ int src_pitch,
+ int dst_pitch,
+ long n)
+{
+ while (n--)
+ {
+ doubles2halfp (dst, src, 1);
+ dst += dst_pitch;
+ src += src_pitch;
+ }
+ return n;
+}
+
+static long
+convert_half_double (char *src,
+ char *dst,
+ int src_pitch,
+ int dst_pitch,
+ long n)
+{
+ while (n--)
+ {
+ halfp2doubles (dst, src, 1);
+ dst += dst_pitch;
+ src += src_pitch;
+ }
+ return n;
+}
+
+void
+babl_base_type_half (void)
+{
+ babl_type_new (
+ "half",
+ "id", BABL_HALF,
+ "bits", 16,
+ NULL);
+
+ babl_conversion_new (
+ babl_type_from_id (BABL_HALF),
+ babl_type_from_id (BABL_DOUBLE),
+ "plane", convert_half_double,
+ NULL
+ );
+
+ babl_conversion_new (
+ babl_type_from_id (BABL_DOUBLE),
+ babl_type_from_id (BABL_HALF),
+ "plane", convert_double_half,
+ NULL
+ );
+}
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
