fix: Adjust colorspace calculations for ACES AP0

The subtracts in PNG_XYZ_from_xy are producing integer overflow with
some valid but extreme xy values.  This re-introduces the previous
checks, but with less limited bounds; sufficient to accomodate the
ACEScg end points (ACES AP1), but not for the ACES AP0 end points.
Those were not working anyway because libpng reads the cHRM parameters
as unsigned values, so they must always be at least 0.

A better solution requires recognizing reasonable negative values (ones
which violate the current spec) and allowing them too, at least on read.

Reviewed-by: Cosmin Truta <[email protected]>
Signed-off-by: John Bowler <[email protected]>
Signed-off-by: Cosmin Truta <[email protected]>

png.c

@@ -1180,22 +1180,66 @@ png_colorspace_sync(png_const_structrp png_ptr, png_inforp info_ptr)
 #endif /* GAMMA */
 
 #ifdef PNG_COLORSPACE_SUPPORTED
-static int
-png_safe_add(png_int_32 *addend0_and_result, png_int_32 addend1,
-      png_int_32 addend2) {
-   /* Safely add three integers.  Returns 0 on success, 1 on overlow.
+static png_int_32
+png_fp_add(png_int_32 addend0, png_int_32 addend1, int *error)
+{
+   /* Safely add two fixed point values setting an error flag and returning 0.5
+    * on overflow.
     * IMPLEMENTATION NOTE: ANSI requires signed overflow not to occur, therefore
     * relying on addition of two positive values producing a negative one is not
     * safe.
     */
-   int addend0 = *addend0_and_result;
-   if (0x7fffffff - addend0 < addend1)
-      return 1;
-   addend0 += addend1;
-   if (0x7fffffff - addend1 < addend2)
-      return 1;
-   *addend0_and_result = addend0 + addend2;
-   return 0;
+   if (addend0 > 0)
+   {
+      if (0x7fffffff - addend0 >= addend1)
+         return addend0+addend1;
+   }
+   else if (addend0 < 0)
+   {
+      if (-0x7fffffff - addend0 <= addend1)
+         return addend0+addend1;
+   }
+   else
+      return addend1;
+
+   *error = 1;
+   return PNG_FP_1/2;
+}
+
+static png_int_32
+png_fp_sub(png_int_32 addend0, png_int_32 addend1, int *error)
+{
+   /* As above but calculate addend0-addend1. */
+   if (addend1 > 0)
+   {
+      if (-0x7fffffff + addend1 <= addend0)
+         return addend0-addend1;
+   }
+   else if (addend1 < 0)
+   {
+      if (0x7fffffff + addend1 >= addend0)
+         return addend0+addend1;
+   }
+   else
+      return addend0;
+
+   *error = 1;
+   return PNG_FP_1/2;
+}
+
+static int
+png_safe_add(png_int_32 *addend0_and_result, png_int_32 addend1,
+      png_int_32 addend2)
+{
+   /* Safely add three integers.  Returns 0 on success, 1 on overflow.  Does not
+    * set the result on overflow.
+    */
+   int error = 0;
+   int result = png_fp_add(*addend0_and_result,
+                           png_fp_add(addend1, addend2, &error),
+                           &error);
+   if (!error) *addend0_and_result = result;
+   return error;
 }
 
 /* Added at libpng-1.5.5 to support read and write of true CIEXYZ values for
@@ -1266,6 +1310,29 @@ png_XYZ_from_xy(png_XYZ *XYZ, const png_xy *xy)
    png_fixed_point red_inverse, green_inverse, blue_scale;
    png_fixed_point left, right, denominator;
 
+   /* Check xy and, implicitly, z.  Note that wide gamut color spaces typically
+    * have end points with 0 tristimulus values (these are impossible end
+    * points, but they are used to cover the possible colors).  We check
+    * xy->whitey against 5, not 0, to avoid a possible integer overflow.
+    *
+    * The limits here will *not* accept ACES AP0, where bluey is -7700
+    * (-0.0770) because the PNG spec itself requires the xy values to be
+    * unsigned.  whitey is also required to be 5 or more to avoid overflow.
+    *
+    * Instead the upper limits have been relaxed to accomodate ACES AP1 where
+    * redz ends up as -600 (-0.006).  ProPhotoRGB was already "in range."
+    * The new limit accomodates the AP0 and AP1 ranges for z but not AP0 redy.
+    */
+   const png_fixed_point fpLimit = PNG_FP_1+(PNG_FP_1/10);
+   if (xy->redx   < 0 || xy->redx > fpLimit) return 1;
+   if (xy->redy   < 0 || xy->redy > fpLimit-xy->redx) return 1;
+   if (xy->greenx < 0 || xy->greenx > fpLimit) return 1;
+   if (xy->greeny < 0 || xy->greeny > fpLimit-xy->greenx) return 1;
+   if (xy->bluex  < 0 || xy->bluex > fpLimit) return 1;
+   if (xy->bluey  < 0 || xy->bluey > fpLimit-xy->bluex) return 1;
+   if (xy->whitex < 0 || xy->whitex > fpLimit) return 1;
+   if (xy->whitey < 5 || xy->whitey > fpLimit-xy->whitex) return 1;
+
    /* The reverse calculation is more difficult because the original tristimulus
     * value had 9 independent values (red,green,blue)x(X,Y,Z) however only 8
     * derived values were recorded in the cHRM chunk;
@@ -1409,18 +1476,23 @@ png_XYZ_from_xy(png_XYZ *XYZ, const png_xy *xy)
     *  (green-x - blue-x)*(red-y - blue-y)-(green-y - blue-y)*(red-x - blue-x)
     *
     * Accuracy:
-    * The input values have 5 decimal digits of accuracy.  The values are all in
-    * the range 0 < value < 1, so simple products are in the same range but may
-    * need up to 10 decimal digits to preserve the original precision and avoid
-    * underflow.  Because we are using a 32-bit signed representation we cannot
-    * match this; the best is a little over 9 decimal digits, less than 10.
+    * The input values have 5 decimal digits of accuracy.
+    *
+    * In the previous implementation the values were all in the range 0 < value
+    * < 1, so simple products are in the same range but may need up to 10
+    * decimal digits to preserve the original precision and avoid underflow.
+    * Because we are using a 32-bit signed representation we cannot match this;
+    * the best is a little over 9 decimal digits, less than 10.
+    *
+    * This range has now been extended to allow values up to 1.1, or 110,000 in
+    * fixed point.
     *
     * The approach used here is to preserve the maximum precision within the
     * signed representation.  Because the red-scale calculation above uses the
-    * difference between two products of values that must be in the range -1..+1
-    * it is sufficient to divide the product by 7; ceil(100,000/32767*2).  The
-    * factor is irrelevant in the calculation because it is applied to both
-    * numerator and denominator.
+    * difference between two products of values that must be in the range
+    * -1.1..+1.1 it is sufficient to divide the product by 8;
+    * ceil(121,000/32767*2).  The factor is irrelevant in the calculation
+    * because it is applied to both numerator and denominator.
     *
     * Note that the values of the differences of the products of the
     * chromaticities in the above equations tend to be small, for example for
@@ -1442,49 +1514,61 @@ png_XYZ_from_xy(png_XYZ *XYZ, const png_xy *xy)
     *  Adobe Wide Gamut RGB
     *    0.258728243040113 0.724682314948566 0.016589442011321
     */
-   /* By the argument, above overflow should be impossible here. The return
-    * value of 2 indicates an internal error to the caller.
+   int error = 0;
+
+   /* By the argument above overflow should be impossible here, however the
+    * code now simply returns a failure code.  The xy subtracts in the arguments
+    * to png_muldiv are *not* checked for overflow because the checks at the
+    * start guarantee they are in the range 0..110000 and png_fixed_point is a
+    * 32-bit signed number.
     */
-   if (png_muldiv(&left, xy->greenx-xy->bluex, xy->redy - xy->bluey, 7) == 0)
+   if (png_muldiv(&left, xy->greenx-xy->bluex, xy->redy - xy->bluey, 8) == 0)
       return 1;
-   if (png_muldiv(&right, xy->greeny-xy->bluey, xy->redx - xy->bluex, 7) == 0)
+   if (png_muldiv(&right, xy->greeny-xy->bluey, xy->redx - xy->bluex, 8) == 0)
       return 1;
-   denominator = left - right;
+   denominator = png_fp_sub(left, right, &error);
+   if (error) return 1;
 
    /* Now find the red numerator. */
-   if (png_muldiv(&left, xy->greenx-xy->bluex, xy->whitey-xy->bluey, 7) == 0)
+   if (png_muldiv(&left, xy->greenx-xy->bluex, xy->whitey-xy->bluey, 8) == 0)
       return 1;
-   if (png_muldiv(&right, xy->greeny-xy->bluey, xy->whitex-xy->bluex, 7) == 0)
+   if (png_muldiv(&right, xy->greeny-xy->bluey, xy->whitex-xy->bluex, 8) == 0)
       return 1;
 
    /* Overflow is possible here and it indicates an extreme set of PNG cHRM
     * chunk values.  This calculation actually returns the reciprocal of the
     * scale value because this allows us to delay the multiplication of white-y
     * into the denominator, which tends to produce a small number.
     */
-   if (png_muldiv(&red_inverse, xy->whitey, denominator, left-right) == 0 ||
+   if (png_muldiv(&red_inverse, xy->whitey, denominator,
+                  png_fp_sub(left, right, &error)) == 0 || error ||
        red_inverse <= xy->whitey /* r+g+b scales = white scale */)
       return 1;
 
    /* Similarly for green_inverse: */
-   if (png_muldiv(&left, xy->redy-xy->bluey, xy->whitex-xy->bluex, 7) == 0)
+   if (png_muldiv(&left, xy->redy-xy->bluey, xy->whitex-xy->bluex, 8) == 0)
       return 1;
-   if (png_muldiv(&right, xy->redx-xy->bluex, xy->whitey-xy->bluey, 7) == 0)
+   if (png_muldiv(&right, xy->redx-xy->bluex, xy->whitey-xy->bluey, 8) == 0)
       return 1;
-   if (png_muldiv(&green_inverse, xy->whitey, denominator, left-right) == 0 ||
+   if (png_muldiv(&green_inverse, xy->whitey, denominator,
+                  png_fp_sub(left, right, &error)) == 0 || error ||
        green_inverse <= xy->whitey)
       return 1;
 
    /* And the blue scale, the checks above guarantee this can't overflow but it
     * can still produce 0 for extreme cHRM values.
     */
-   blue_scale = png_reciprocal(xy->whitey) - png_reciprocal(red_inverse) -
-       png_reciprocal(green_inverse);
-   if (blue_scale <= 0)
+   blue_scale = png_fp_sub(png_fp_sub(png_reciprocal(xy->whitey),
+                                      png_reciprocal(red_inverse), &error),
+                           png_reciprocal(green_inverse), &error);
+   if (error || blue_scale <= 0)
       return 1;
 
 
-   /* And fill in the png_XYZ: */
+   /* And fill in the png_XYZ.  Again the subtracts are safe because of the
+    * checks on the xy values at the start (the subtracts just calculate the
+    * corresponding z values.)
+    */
    if (png_muldiv(&XYZ->red_X, xy->redx, PNG_FP_1, red_inverse) == 0)
       return 1;
    if (png_muldiv(&XYZ->red_Y, xy->redy, PNG_FP_1, red_inverse) == 0)