[libpng16] New 'tools' directory containing tools used to generate libpng code.

contrib/tools/README.txt

@@ -0,0 +1,26 @@
+This directory (contrib/tools) contains tools used by the authors of libpng.
+
+Code and data placed in this directory is not required to build libpng,
+however the code in this directory has been used to generate data or code in
+the body of the libpng source.  The source code idenftifies where this has
+been done.  Code in this directory may not compile on all operating systems
+that libpng supports.
+
+NO COPYRIGHT RIGHTS ARE CLAIMED TO ANY OF THE FILES IN THIS DIRECTORY.
+
+To the extent possible under law, the authors have waived all copyright and
+related or neighboring rights to this work.  This work is published from:
+United States.
+
+The files may be used freely in any way.
+
+The source code and comments in this directory are the original work of the
+people named below.  No other person or organization has made contributions to
+the work in this directory.
+
+ORIGINAL AUTHORS
+    The following people have contributed to the code in this directory.  None
+    of the people below claim any rights with regard to the contents of this
+    directory.
+
+    John Bowler <jbowler@acm.org>

contrib/tools/cvtcolor.c

@@ -0,0 +1,188 @@
+/*-
+ * convert.c
+ *
+ * Last changed in libpng 1.6.0 [(PENDING RELEASE)]
+ *
+ * COPYRIGHT: Written by John Cunningham Bowler, 2011.
+ * To the extent possible under law, the author has waived all copyright and
+ * related or neighboring rights to this work.  This work is published from:
+ * United States.
+ *
+ * Convert 8-bit sRGB or 16-bit linear values to another format.
+ */
+#define _ISOC99_SOURCE 1
+
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <stdio.h>
+
+#include <fenv.h>
+
+#include "sRGB.h"
+
+static void
+usage(const char *prog)
+{
+   fprintf(stderr,
+      "%s: usage: %s [-linear|-sRGB] [-gray|-color] component{1,4}\n",
+      prog, prog);
+   exit(1);
+}
+
+unsigned long
+component(const char *prog, const char *arg, int issRGB)
+{
+   char *ep;
+   unsigned long c = strtoul(arg, &ep, 0);
+
+   if (ep <= arg || *ep || c > 65535 || (issRGB && c > 255))
+   {
+      fprintf(stderr, "%s: %s: invalid component value (%lu)\n", prog, arg, c);
+      usage(prog);
+   }
+
+   return c;
+}
+
+int
+main(int argc, const char **argv)
+{
+   const char *prog = *argv++;
+   int to_linear = 0, to_gray = 0, to_color = 0;
+   int channels = 0;
+   double c[4];
+
+   /* FE_TONEAREST is the IEEE754 round to nearest, preferring even, mode; i.e.
+    * everything rounds to the nearest value except that '.5' rounds to the
+    * nearest even value.
+    */
+   fesetround(FE_TONEAREST);
+
+   c[3] = c[2] = c[1] = c[0] = 0;
+
+   while (--argc > 0 && **argv == '-')
+   {
+      const char *arg = 1+*argv++;
+
+      if (strcmp(arg, "sRGB") == 0)
+         to_linear = 0;
+
+      else if (strcmp(arg, "linear") == 0)
+         to_linear = 1;
+
+      else if (strcmp(arg, "gray") == 0)
+         to_gray = 1, to_color = 0;
+
+      else if (strcmp(arg, "color") == 0)
+         to_gray = 0, to_color = 1;
+
+      else
+         usage(prog);
+   }
+
+   switch (argc)
+   {
+      default:
+         usage(prog);
+         break;
+
+      case 4:
+         c[3] = component(prog, argv[3], to_linear);
+         ++channels;
+      case 3:
+         c[2] = component(prog, argv[2], to_linear);
+         ++channels;
+      case 2:
+         c[1] = component(prog, argv[1], to_linear);
+         ++channels;
+      case 1:
+         c[0] = component(prog, argv[0], to_linear);
+         ++channels;
+         break;
+      }
+
+   if (to_linear)
+   {
+      int i;
+      int components = channels;
+
+      if ((components & 1) == 0)
+         --components;
+
+      for (i=0; i<components; ++i) c[i] = linear_from_sRGB(c[i] / 255);
+      if (components < channels)
+         c[components] = c[components] / 255;
+   }
+
+   else
+   {
+      int i;
+      for (i=0; i<4; ++i) c[i] /= 65535;
+
+      if ((channels & 1) == 0)
+      {
+         double alpha = c[channels-1];
+
+         if (alpha > 0)
+            for (i=0; i<channels-1; ++i) c[i] /= alpha;
+         else
+            for (i=0; i<channels-1; ++i) c[i] = 1;
+      }
+   }
+
+   if (to_gray)
+   {
+      if (channels < 3)
+      {
+         fprintf(stderr, "%s: too few channels (%d) for -gray\n",
+            prog, channels);
+         usage(prog);
+      }
+
+      c[0] = YfromRGB(c[0], c[1], c[2]);
+      channels -= 2;
+   }
+
+   if (to_color)
+   {
+      if (channels > 2)
+      {
+         fprintf(stderr, "%s: too many channels (%d) for -color\n",
+            prog, channels);
+         usage(prog);
+      }
+
+      c[3] = c[1]; /* alpha, if present */
+      c[2] = c[1] = c[0];
+   }
+
+   if (to_linear)
+   {
+      int i;
+      if ((channels & 1) == 0)
+      {
+         double alpha = c[channels-1];
+         for (i=0; i<channels-1; ++i) c[i] *= alpha;
+      }
+
+      for (i=0; i<channels; ++i) c[i] = nearbyint(c[i] * 65535);
+   }
+
+   else /* to sRGB */
+   {
+      int i = (channels+1)&~1;
+      while (--i >= 0)
+         c[i] = sRGB_from_linear(c[i]);
+
+      for (i=0; i<channels; ++i) c[i] = nearbyint(c[i] * 255);
+   }
+
+   {
+      int i;
+      for (i=0; i<channels; ++i) printf(" %g", c[i]);
+   }
+   printf("\n");
+
+   return 0;
+}

contrib/tools/intgamma.sh

@@ -0,0 +1,110 @@
+#!/bin/sh
+#
+# intgamma.sh
+#
+# Last changed in libpng 1.6.0 [(PENDING RELEASE)]
+#
+# COPYRIGHT: Written by John Cunningham Bowler, 2011.
+# To the extent possible under law, the author has waived all copyright and
+# related or neighboring rights to this work.  This work is published from:
+# United States.
+#
+# Shell script to generate png.c 8 and 16 bit log tables (see the code in png.c
+# for details).
+#
+# This script uses the "bc" arbitrary precision calculator to calculate 32 bit
+# fixed point values of logarithms appropriate to finding the log of an 8-bit
+# (0..255) value and a similar table for the exponent calculation.
+#
+# "bc" must be on the path when the script is executed, the math library (-l)
+# must be available
+#
+# function to print out a list of numbers as integers; the function truncates
+# the integers which must be one-per-line
+function print(){
+   awk 'BEGIN{
+      str = ""
+   }
+   {
+      sub("\\.[0-9]*$", "")
+      if ($0 == "")
+         $0 = "0"
+
+      if (str == "")
+         t = "   " $0 "U"
+      else
+         t = str ", " $0 "U"
+
+      if (length(t) >= 80) {
+         print str ","
+         str = "   " $0 "U"
+      } else
+         str = t
+   }
+   END{
+      print str
+   }'
+}
+#
+# The logarithm table.
+cat <<END
+/* 8-bit log table: png_8bit_l2[128]
+ * This is a table of -log(value/255)/log(2) for 'value' in the range 128 to
+ * 255, so it's the base 2 logarithm of a normalized 8-bit floating point
+ * mantissa.  The numbers are 32-bit fractions.
+ */
+static const png_uint_32
+png_8bit_l2[128] =
+{
+END
+#
+bc -lqws <<END | print
+f=65536*65536/l(2)
+for (i=128;i<256;++i) { .5 - l(i/255)*f; }
+END
+echo '};'
+echo
+#
+# The exponent table.
+cat <<END
+/* The 'exp()' case must invert the above, taking a 20-bit fixed point
+ * logarithmic value and returning a 16 or 8-bit number as appropriate.  In
+ * each case only the low 16 bits are relevant - the fraction - since the
+ * integer bits (the top 4) simply determine a shift.
+ *
+ * The worst case is the 16-bit distinction between 65535 and 65534, this
+ * requires perhaps spurious accuracty in the decoding of the logarithm to
+ * distinguish log2(65535/65534.5) - 10^-5 or 17 bits.  There is little chance
+ * of getting this accuracy in practice.
+ *
+ * To deal with this the following exp() function works out the exponent of the
+ * frational part of the logarithm by using an accurate 32-bit value from the
+ * top four fractional bits then multiplying in the remaining bits.
+ */
+static const png_uint_32
+png_32bit_exp[16] =
+{
+END
+#
+bc -lqws <<END | print
+f=l(2)/16
+for (i=0;i<16;++i) {
+   x = .5 + e(-i*f)*2^32;
+   if (x >= 2^32) x = 2^32-1;
+   x;
+}
+END
+echo '};'
+echo
+#
+# And the table of adjustment values.
+cat <<END
+/* Adjustment table; provided to explain the numbers in the code below. */
+#if 0
+END
+bc -lqws <<END | awk '{ printf "%5d %s\n", 12-NR, $0 }'
+for (i=11;i>=0;--i){
+   (1 - e(-(2^i)/65536*l(2))) * 2^(32-i)
+}
+END
+echo '#endif'

contrib/tools/makesRGB.c

@@ -0,0 +1,430 @@
+/* makesRGB.c -- build sRGB-to-linear and linear-to-sRGB conversion tables
+ *
+ * Last changed in libpng 1.6.0 [(PENDING RELEASE)]
+ *
+ * COPYRIGHT: Written by John Cunningham Bowler, 2011.
+ * To the extent possible under law, the author has waived all copyright and
+ * related or neighboring rights to this work.  This work is published from:
+ * United States.
+ *
+ * Make a table to convert 8-bit sRGB encoding values into the closest 16-bit
+ * linear value.
+ *
+ * Make two tables to take a linear value scaled to 255*65535 and return an
+ * approximation to the 8-bit sRGB encoded value.  Calculate the error in these
+ * tables and display it.
+ */
+#define _C99_SOURCE 1
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+
+/* pngpriv.h includes the definition of 'PNG_sRGB_FROM_LINEAR' which is required
+ * to verify the actual code.
+ */
+#include "../../pngpriv.h"
+
+#include "sRGB.h"
+
+/* The tables are declared 'const' in pngpriv.h, so this redefines the tables to
+ * be used.
+ */
+#define png_sRGB_table sRGB_table
+#define png_sRGB_base sRGB_base
+#define png_sRGB_delta sRGB_delta
+
+static png_uint_16 png_sRGB_table[256];
+static png_uint_16 png_sRGB_base[512];
+static png_byte png_sRGB_delta[512];
+
+static const unsigned int max_input = 255*65535;
+
+double
+fsRGB(double l)
+{
+   return sRGB_from_linear(l/max_input);
+}
+
+double
+sRGB(unsigned int i)
+{
+   return fsRGB(i);
+}
+
+double
+finvsRGB(unsigned int i)
+{
+   return 65535 * linear_from_sRGB(i/255.);
+}
+
+png_uint_16
+invsRGB(unsigned int i)
+{
+   unsigned int x = nearbyint(finvsRGB(i));
+
+   if (x > 65535)
+   {
+      fprintf(stderr, "invsRGB(%u) overflows to %u\n", i, x);
+      exit(1);
+   }
+
+   return (png_uint_16)x;
+}
+
+int
+main(int argc, char **argv)
+{
+   unsigned int i, i16, ibase;
+   double min_error = 0;
+   double max_error = 0;
+   double min_error16 = 0;
+   double max_error16 = 0;
+   double adjust;
+   double adjust_lo = 0.4, adjust_hi = 0.6, adjust_mid = 0.5;
+   unsigned int ec_lo = 0, ec_hi = 0, ec_mid = 0;
+   unsigned int error_count = 0;
+   unsigned int error_count16 = 0;
+   int test_only = 0;
+
+   if (argc > 1)
+      test_only = strcmp("--test", argv[1]) == 0;
+
+   /* Initialize the encoding table first. */
+   for (i=0; i<256; ++i)
+   {
+      png_sRGB_table[i] = invsRGB(i);
+   }
+
+   /* Now work out the decoding tables (this is where the error comes in because
+    * there are 512 set points and 512 straight lines between them.)
+    */
+   for (;;)
+   {
+      if (ec_lo == 0)
+         adjust = adjust_lo;
+
+      else if (ec_hi == 0)
+         adjust = adjust_hi;
+
+      else if (ec_mid == 0)
+         adjust = adjust_mid;
+
+      else if (ec_mid < ec_hi)
+         adjust = (adjust_mid + adjust_hi)/2;
+
+      else if (ec_mid < ec_lo)
+         adjust = (adjust_mid + adjust_lo)/2;
+
+      else
+      {
+         fprintf(stderr, "not reached: %u .. %u .. %u\n", ec_lo, ec_mid, ec_hi);
+         exit(1);
+      }
+
+      /* Calculate the table using the current 'adjust' */
+      for (i=0; i<=511; ++i)
+      {
+         double lo = 255 * sRGB(i << 15);
+         double hi = 255 * sRGB((i+1) << 15);
+         unsigned int calc;
+
+         calc = nearbyint((lo+adjust) * 256);
+         if (calc > 65535)
+         {
+            fprintf(stderr, "table[%d][0]: overflow %08x (%d)\n", i, calc,
+               calc);
+            exit(1);
+         }
+         png_sRGB_base[i] = calc;
+
+         calc = nearbyint((hi-lo) * 32);
+         if (calc > 255)
+         {
+            fprintf(stderr, "table[%d][1]: overflow %08x (%d)\n", i, calc,
+               calc);
+            exit(1);
+         }
+         png_sRGB_delta[i] = calc;
+      }
+
+      /* Check the 16-bit linear values alone: */
+      error_count16 = 0;
+      for (i16=0; i16 <= 65535; ++i16)
+      {
+         unsigned int i = 255*i16;
+         unsigned int iexact = nearbyint(255*sRGB(i));
+         unsigned int icalc = PNG_sRGB_FROM_LINEAR(i);
+
+         if (icalc != iexact)
+            ++error_count16;
+      }
+
+      /* Now try changing the adjustment. */
+      if (ec_lo == 0)
+         ec_lo = error_count16;
+
+      else if (ec_hi == 0)
+         ec_hi = error_count16;
+
+      else if (ec_mid == 0)
+      {
+         ec_mid = error_count16;
+         printf("/* initial error counts: %u .. %u .. %u */\n", ec_lo, ec_mid,
+            ec_hi);
+      }
+
+      else if (error_count16 < ec_mid)
+      {
+         printf("/* adjust (mid ): %f: %u -> %u */\n", adjust, ec_mid,
+            error_count16);
+         ec_mid = error_count16;
+         adjust_mid = adjust;
+      }
+
+      else if (adjust < adjust_mid && error_count16 < ec_lo)
+      {
+         printf("/* adjust (low ): %f: %u -> %u */\n", adjust, ec_lo,
+            error_count16);
+         ec_lo = error_count16;
+         adjust_lo = adjust;
+      }
+
+      else if (adjust > adjust_mid && error_count16 < ec_hi)
+      {
+         printf("/* adjust (high): %f: %u -> %u */\n", adjust, ec_hi,
+            error_count16);
+         ec_hi = error_count16;
+         adjust_hi = adjust;
+      }
+
+      else
+      {
+         adjust = adjust_mid;
+         printf("/* adjust: %f: %u */\n", adjust, ec_mid);
+         break;
+      }
+   }
+
+   /* For each entry in the table try to adjust it to minimize the error count
+    * in that entry.  Each entry corresponds to 128 input values.
+    */
+   for (ibase=0; ibase<65536; ibase+=128)
+   {
+      png_uint_16 base = png_sRGB_base[ibase >> 7], trybase = base, ob=base;
+      png_byte delta = png_sRGB_delta[ibase >> 7], trydelta = delta, od=delta;
+      unsigned int ecbase = 0, eco;
+
+      for (;;)
+      {
+         png_sRGB_base[ibase >> 7] = trybase;
+         png_sRGB_delta[ibase >> 7] = trydelta;
+
+         /* Check the 16-bit linear values alone: */
+         error_count16 = 0;
+         for (i16=ibase; i16 < ibase+128; ++i16)
+         {
+            unsigned int i = 255*i16;
+            unsigned int iexact = nearbyint(255*sRGB(i));
+            unsigned int icalc = PNG_sRGB_FROM_LINEAR(i);
+
+            if (icalc != iexact)
+               ++error_count16;
+         }
+
+         if (error_count16 == 0)
+            break;
+
+         if (ecbase == 0)
+         {
+            eco = ecbase = error_count16;
+            ++trybase; /* First test */
+         }
+
+         else if (error_count16 < ecbase)
+         {
+            if (trybase > base)
+            {
+               base = trybase;
+               ++trybase;
+            }
+            else if (trybase < base)
+            {
+               base = trybase;
+               --trybase;
+            }
+            else if (trydelta > delta)
+            {
+               delta = trydelta;
+               ++trydelta;
+            }
+            else if (trydelta < delta)
+            {
+               delta = trydelta;
+               --trydelta;
+            }
+            else
+            {
+               fprintf(stderr, "makesRGB: impossible\n");
+               exit(1);
+            }
+            ecbase = error_count16;
+         }
+
+         else
+         {
+            if (trybase > base)
+               trybase = base-1;
+            else if (trybase < base)
+            {
+               trybase = base;
+               ++trydelta;
+            }
+            else if (trydelta > delta)
+               trydelta = delta-1;
+            else if (trydelta < delta)
+               break; /* end of tests */
+         }
+      }
+
+      png_sRGB_base[ibase >> 7] = base;
+      png_sRGB_delta[ibase >> 7] = delta;
+      if (base != ob || delta != od)
+      {
+         printf("/* table[%u]={%u,%u} -> {%u,%u} %u -> %u errors */\n",
+            ibase>>7, ob, od, base, delta, eco, ecbase);
+      }
+      else if (0)
+         printf("/* table[%u]={%u,%u} %u errors */\n", ibase>>7, ob, od,
+            ecbase);
+   }
+
+   /* Only do the full (slow) test at the end: */
+   min_error = -.4999;
+   max_error = .4999;
+   error_count = 0;
+
+   for (i=0; i <= max_input; ++i)
+   {
+      unsigned int iexact = nearbyint(255*sRGB(i));
+      unsigned int icalc = PNG_sRGB_FROM_LINEAR(i);
+
+      if (icalc != iexact)
+      {
+         double err = 255*sRGB(i) - icalc;
+
+         if (err > (max_error+.001) || err < (min_error-.001))
+         {
+            printf(
+               "/* 0x%08x: exact: %3d, got: %3d [tables: %08x, %08x] (%f) */\n",
+               i, iexact, icalc, png_sRGB_base[i>>15],
+               png_sRGB_delta[i>>15], err);
+         }
+
+         ++error_count;
+         if (err > max_error)
+            max_error = err;
+         else if (err < min_error)
+            min_error = err;
+      }
+   }
+
+   /* Re-check the 16-bit cases too, including the warning if there is an error
+    * bigger than 1.
+    */
+   error_count16 = 0;
+   max_error16 = 0;
+   min_error16 = 0;
+   for (i16=0; i16 <= 65535; ++i16)
+   {
+      unsigned int i = 255*i16;
+      unsigned int iexact = nearbyint(255*sRGB(i));
+      unsigned int icalc = PNG_sRGB_FROM_LINEAR(i);
+
+      if (icalc != iexact)
+      {
+         double err = 255*sRGB(i) - icalc;
+
+         ++error_count16;
+         if (err > max_error16)
+            max_error16 = err;
+         else if (err < min_error16)
+            min_error16 = err;
+
+         if (abs(icalc - iexact) > 1)
+            printf(
+               "/* 0x%04x: exact: %3d, got: %3d [tables: %08x, %08x] (%f) */\n",
+               i16, iexact, icalc, png_sRGB_base[i>>15],
+               png_sRGB_delta[i>>15], err);
+      }
+   }
+
+   /* Check the round trip for each 8-bit sRGB value. */
+   for (i16=0; i16 <= 255; ++i16)
+   {
+      unsigned int i = 255 * png_sRGB_table[i16];
+      unsigned int iexact = nearbyint(255*sRGB(i));
+      unsigned int icalc = PNG_sRGB_FROM_LINEAR(i);
+
+      if (i16 != iexact)
+      {
+         fprintf(stderr, "8-bit rounding error: %d -> %d\n", i16, iexact);
+         exit(1);
+      }
+
+      if (icalc != i16)
+      {
+         double finv = finvsRGB(i16);
+
+         printf("/* 8-bit roundtrip error: %d -> %f -> %d(%f) */\n",
+            i16, finv, icalc, fsRGB(255*finv));
+      }
+   }
+
+
+   printf("/* error: %g - %g, %u (%g%%) of readings inexact */\n",
+      min_error, max_error, error_count, (100.*error_count)/max_input);
+   printf("/* 16-bit error: %g - %g, %u (%g%%) of readings inexact */\n",
+      min_error16, max_error16, error_count16, (100.*error_count16)/65535);
+
+   if (!test_only)
+   {
+      printf("PNG_CONST png_uint_16 png_sRGB_table[256] =\n{\n   ");
+      for (i=0; i<255; )
+      {
+         do
+         {
+            printf("%d,", png_sRGB_table[i++]);
+         }
+         while ((i & 0x7) != 0 && i<255);
+         if (i<255) printf("\n   ");
+      }
+      printf("%d\n};\n\n", png_sRGB_table[i]);
+
+
+      printf("PNG_CONST png_uint_16 png_sRGB_base[512] =\n{\n   ");
+      for (i=0; i<511; )
+      {
+         do
+         {
+            printf("%d,", png_sRGB_base[i++]);
+         }
+         while ((i & 0x7) != 0 && i<511);
+         if (i<511) printf("\n   ");
+      }
+      printf("%d\n};\n\n", png_sRGB_base[i]);
+
+      printf("PNG_CONST png_byte png_sRGB_delta[512] =\n{\n   ");
+      for (i=0; i<511; )
+      {
+         do
+         {
+            printf("%d,", png_sRGB_delta[i++]);
+         }
+         while ((i & 0xf) != 0 && i<511);
+         if (i<511) printf("\n   ");
+      }
+      printf("%d\n};\n\n", png_sRGB_delta[i]);
+   }
+
+   return 0;
+}

contrib/tools/sRGB.h

@@ -0,0 +1,48 @@
+/*-
+ * sRGB.h
+ *
+ * Last changed in libpng 1.6.0 [(PENDING RELEASE)]
+ *
+ * COPYRIGHT: Written by John Cunningham Bowler, 2011.
+ * To the extent possible under law, the author has waived all copyright and
+ * related or neighboring rights to this work.  This work is published from:
+ * United States.
+ *
+ * Utility file; not actually a header, this contains definitions of sRGB
+ * calculation functions for inclusion in those test programs that need them.
+ *
+ * All routines take and return a floating point value in the range
+ * 0 to 1.0, doing a calculation according to the sRGB specification
+ * (in fact the source of the numbers is the wikipedia article at
+ * http://en.wikipedia.org/wiki/SRGB).
+ */
+static double
+sRGB_from_linear(double l)
+{
+   if (l <= 0.0031308)
+      l *= 12.92;
+
+   else
+      l = 1.055 * pow(l, 1/2.4) - 0.055;
+
+   return l;
+}
+
+static double
+linear_from_sRGB(double s)
+{
+   if (s <= 0.04045)
+      return s / 12.92;
+
+   else
+      return pow((s+0.055)/1.055, 2.4);
+}
+
+static double
+YfromRGB(double r, double g, double b)
+{
+   /* Use the sRGB (rounded) coefficients for Rlinear, Glinear, Blinear to get
+    * the CIE Y value (also linear).
+    */
+   return 0.2126 * r + 0.7152 * g + 0.0722 * b;
+}