Fix handling of floating point imprecision in unnormalizing samples

git-svn-id: http://svn.code.sf.net/p/netpbm/code/trunk@3514 9d0c8265-081b-0410-96cb-a4ca84ce46f8

12 files changed, 48 insertions, 20 deletions

diff --git a/doc/HISTORY b/doc/HISTORY
index 40fa133b..94eafeee 100644
--- a/doc/HISTORY
+++ b/doc/HISTORY
@@ -19,8 +19,8 @@ not yet  BJH  Release 10.86.00
               pamstretch-gen: Use -dropedge on the 'pamstretch' piece for
               better looking output.
 
-              libnetpbm: Fix ROUNDU so it properly rounds up fractions of
-              exactly half.  Affects many programs.
+              various: Fix unnormalizing code so a value exactly between two
+              sample values rounds consistently up.  Affects many programs.
 
               pamstretch: Reject very large scale factors instead of producing
               incorrect output.
diff --git a/generator/pgmkernel.c b/generator/pgmkernel.c
index ec634c16..37072c38 100644
--- a/generator/pgmkernel.c
+++ b/generator/pgmkernel.c
@@ -223,11 +223,13 @@ main(int argc, const char * argv[]) {
 
         unsigned int col;
         for (col = 0; col < (cmdline.cols +1) / 2; ++col) {
+            double const epsilon = 1e-15;
             double const dx2 = SQR(col - xcenter);
 
             double const normalized = t(dx2, dy2, cmdline.weight) / 2 / tMax;
 
-            gray const grayval = ROUNDU(cmdline.maxval * (0.5 + normalized));
+            gray const grayval =
+                ROUNDU(cmdline.maxval * (0.5 + normalized + epsilon));
 
             halfKernel[arow][col                   ] = grayval;
             halfKernel[arow][cmdline.cols - col - 1] = grayval;
diff --git a/generator/pgmmake.c b/generator/pgmmake.c
index 3843e316..ae706639 100644
--- a/generator/pgmmake.c
+++ b/generator/pgmmake.c
@@ -117,7 +117,7 @@ main(int argc, const char ** const argv) {
 
     parseCommandLine(argc, argv, &cmdline);
 
-    grayLevel = ROUNDU(cmdline.grayLevel * cmdline.maxval);
+    grayLevel = pgm_unnormalize(cmdline.grayLevel, cmdline.maxval);
 
     pgm_writepgminit(stdout, cmdline.cols, cmdline.rows, cmdline.maxval, 0);
 
diff --git a/lib/colorname.c b/lib/colorname.c
index 9400adf7..fe580cb9 100644
--- a/lib/colorname.c
+++ b/lib/colorname.c
@@ -247,9 +247,9 @@ pm_parse_dictionary_name(char    const colorname[],
 
     pm_parse_dictionary_namen(colorname, color);
 
-    r = ROUNDU(color[PAM_RED_PLANE] * maxval);
-    g = ROUNDU(color[PAM_GRN_PLANE] * maxval);
-    b = ROUNDU(color[PAM_BLU_PLANE] * maxval);
+    r = ppm_unnormalize(color[PAM_RED_PLANE], maxval);
+    g = ppm_unnormalize(color[PAM_GRN_PLANE], maxval);
+    b = ppm_unnormalize(color[PAM_BLU_PLANE], maxval);
 
     if (!closeOk) {
         if (maxval != PAM_COLORFILE_MAXVAL) {
diff --git a/lib/libpamcolor.c b/lib/libpamcolor.c
index f9add1a1..e1a24c66 100644
--- a/lib/libpamcolor.c
+++ b/lib/libpamcolor.c
@@ -366,9 +366,7 @@ pnm_parsecolor2(const char * const colorname,
 
     color = pnm_parsecolorn(colorname);
 
-    retval[PAM_RED_PLANE] = ROUNDU(color[PAM_RED_PLANE] * maxval);
-    retval[PAM_GRN_PLANE] = ROUNDU(color[PAM_GRN_PLANE] * maxval);
-    retval[PAM_BLU_PLANE] = ROUNDU(color[PAM_BLU_PLANE] * maxval);
+    pnm_unnormalizetuple(&pam, color, retval);
 
     if (!closeOk) {
         warnIfNotExact(colorname, retval, color, maxval, PAM_RED_PLANE);
diff --git a/lib/libpamn.c b/lib/libpamn.c
index 8ae57037..ae28283a 100644
--- a/lib/libpamn.c
+++ b/lib/libpamn.c
@@ -414,6 +414,19 @@ pnm_writepamn(struct pam * const pamP,
 
 
 
+samplen
+pnm_normalized_sample(struct pam * const pamP,
+                      sample       const sample) {
+    return (samplen)sample/pamP->maxval;
+}
+
+sample
+pnm_unnormalized_sample(struct pam * const pamP,
+                        samplen      const sampleVal) {
+    double const epsilon = 1e-6;
+    return (sample)((sampleVal + epsilon) * pamP->maxval + 0.5);
+}
+
 void
 pnm_normalizetuple(struct pam * const pamP,
                    tuple        const tuple,
@@ -422,7 +435,7 @@ pnm_normalizetuple(struct pam * const pamP,
     unsigned int plane;
 
     for (plane = 0; plane < pamP->depth; ++plane)
-        tuplen[plane] = (samplen)tuple[plane] / pamP->maxval;
+        tuplen[plane] = pnm_normalized_sample(pamP, tuple[plane]);
 }
 
 
@@ -435,7 +448,7 @@ pnm_unnormalizetuple(struct pam * const pamP,
     unsigned int plane;
 
     for (plane = 0; plane < pamP->depth; ++plane)
-        tuple[plane] = tuplen[plane] * pamP->maxval + 0.5;
+        tuple[plane] = pnm_unnormalized_sample(pamP, tuplen[plane]);
 }
 
 
diff --git a/lib/libppmcolor.c b/lib/libppmcolor.c
index 7079482c..c0a88dc2 100644
--- a/lib/libppmcolor.c
+++ b/lib/libppmcolor.c
@@ -553,9 +553,9 @@ ppm_color_from_hsv(struct hsv const hsv,
         }
     }
     PPM_ASSIGN(retval,
-               ROUNDU(R * maxval),
-               ROUNDU(G * maxval),
-               ROUNDU(B * maxval));
+               ppm_unnormalize(R, maxval),
+               ppm_unnormalize(G, maxval),
+               ppm_unnormalize(B, maxval));
 
     return retval;
 }
diff --git a/lib/pam.h b/lib/pam.h
index 74b20f46..0055858b 100644
--- a/lib/pam.h
+++ b/lib/pam.h
@@ -438,6 +438,13 @@ void
 pnm_writepamn(struct pam * const pamP,
               tuplen **    const tuplenarray);
 
+samplen
+pnm_normalized_sample(struct pam * const pamP,
+                      sample       const sample);
+
+sample
+pnm_unnormalized_sample(struct pam * const pamP,
+                        samplen      const sampleVal);
 
 void
 pnm_normalizetuple(struct pam * const pamP,
diff --git a/lib/pgm.h b/lib/pgm.h
index 81e86cfb..d4655239 100644
--- a/lib/pgm.h
+++ b/lib/pgm.h
@@ -42,6 +42,8 @@ typedef unsigned int gray;
 #define PGM_OVERALLMAXVAL 65535
 #define PGM_MAXMAXVAL 255
 
+#define pgm_unnormalize(value, maxval) \
+  ((gray)((value + 1e-6) * (maxval) + 0.5))
 
 /* Magic constants. */
 
diff --git a/lib/pnm.h b/lib/pnm.h
index 3b490552..f6f8137e 100644
--- a/lib/pnm.h
+++ b/lib/pnm.h
@@ -19,6 +19,7 @@ typedef pixel xel;
 typedef pixval xelval;
 #define PNM_OVERALLMAXVAL PPM_OVERALLMAXVAL
 #define PNM_MAXMAXVAL PPM_MAXMAXVAL
+#define pnm_unnormalize ppm_unnormalize
 #define PNM_GET1(x) PPM_GETB(x)
 #define PNM_GETR(x) PPM_GETR(x)
 #define PNM_GETG(x) PPM_GETG(x)
diff --git a/lib/ppm.h b/lib/ppm.h
index 7c483b59..9fc90bb3 100644
--- a/lib/ppm.h
+++ b/lib/ppm.h
@@ -23,6 +23,9 @@ typedef gray pixval;
 
 #define PPM_OVERALLMAXVAL PGM_OVERALLMAXVAL
 #define PPM_MAXMAXVAL PGM_MAXMAXVAL
+
+#define ppm_unnormalize pgm_unnormalize
+
 typedef struct {
     pixval r, g, b;
 } pixel;
diff --git a/lib/util/pm_c_util.h b/lib/util/pm_c_util.h
index 718be369..a093adb6 100644
--- a/lib/util/pm_c_util.h
+++ b/lib/util/pm_c_util.h
@@ -19,11 +19,13 @@
 #undef ROUND
 #define ROUND(X) (((X) >= 0) ? (int)((X)+0.5) : (int)((X)-0.5))
 #undef ROUNDU
-#define ROUNDU(X) ((unsigned int)((X)+0.5+1e-15))
-    /* We add the 1e-15 epsilon because for repeatability we want exactly half
-       to round up.  But the imprecision of the floating point arithmetic can
-       cause a number that is supposed to be exactly half to be slightly less
-       than half in reality.  Without this adjustment, it would round down.
+#define ROUNDU(X) ((unsigned int)((X)+0.5))
+    /* Note that imprecision in floating point arithmetic can make an exact
+       half fractional part round down instead of up.  What should be
+       1000.5 might actually be 1000.49999999999.
+
+       Use 'pnm_unnormalized_sample' instead of ROUNDU(samplen*maxval) to get
+       consistent rounding up when you are unnormalizing sample values.
     */
 
 /* ROUNDUP rounds up to a specified multiple.  E.g. ROUNDUP(22, 8) == 24 */