[devel] Fixes to rgb_to_gray and cHRM XYZ APIs

This commit is contained in:
John Bowler
2011-08-25 16:19:44 -05:00
committed by Glenn Randers-Pehrson
parent 0c03fc6f75
commit 736f40f459
13 changed files with 2055 additions and 514 deletions

339
png.c
View File

@@ -617,13 +617,13 @@ png_get_copyright(png_const_structp png_ptr)
#else
# ifdef __STDC__
return PNG_STRING_NEWLINE \
"libpng version 1.5.5beta06 - August 17, 2011" PNG_STRING_NEWLINE \
"libpng version 1.5.5beta06 - August 25, 2011" PNG_STRING_NEWLINE \
"Copyright (c) 1998-2011 Glenn Randers-Pehrson" PNG_STRING_NEWLINE \
"Copyright (c) 1996-1997 Andreas Dilger" PNG_STRING_NEWLINE \
"Copyright (c) 1995-1996 Guy Eric Schalnat, Group 42, Inc." \
PNG_STRING_NEWLINE;
# else
return "libpng version 1.5.5beta06 - August 17, 2011\
return "libpng version 1.5.5beta06 - August 25, 2011\
Copyright (c) 1998-2011 Glenn Randers-Pehrson\
Copyright (c) 1996-1997 Andreas Dilger\
Copyright (c) 1995-1996 Guy Eric Schalnat, Group 42, Inc.";
@@ -713,18 +713,9 @@ png_access_version_number(void)
#if defined(PNG_READ_SUPPORTED) || defined(PNG_WRITE_SUPPORTED)
# ifdef PNG_SIZE_T
/* Added at libpng version 1.2.6 */
PNG_EXTERN png_size_t PNGAPI png_convert_size PNGARG((size_t size));
png_size_t PNGAPI
png_convert_size(size_t size)
{
if (size > (png_size_t)-1)
PNG_ABORT(); /* We haven't got access to png_ptr, so no png_error() */
return ((png_size_t)size);
}
# endif /* PNG_SIZE_T */
/* png_convert_size: a PNGAPI but no longer in png.h, so deleted
* at libpng 1.5.5!
*/
/* Added at libpng version 1.2.34 and 1.4.0 (moved from pngset.c) */
# ifdef PNG_CHECK_cHRM_SUPPORTED
@@ -798,6 +789,326 @@ png_check_cHRM_fixed(png_structp png_ptr,
}
# endif /* PNG_CHECK_cHRM_SUPPORTED */
#ifdef PNG_cHRM_SUPPORTED
/* Added at libpng-1.5.5 to support read and write of true CIEXYZ values for
* cHRM, as opposed to using chromaticities. These internal APIs return
* non-zero on a parameter error. The X, Y and Z values are required to be
* positive and less than 1.0.
*/
int png_xy_from_XYZ(png_xy *xy, png_XYZ XYZ)
{
png_int_32 d, dwhite, whiteX, whiteY;
d = XYZ.redX + XYZ.redY + XYZ.redZ;
if (!png_muldiv(&xy->redx, XYZ.redX, PNG_FP_1, d)) return 1;
if (!png_muldiv(&xy->redy, XYZ.redY, PNG_FP_1, d)) return 1;
dwhite = d;
whiteX = XYZ.redX;
whiteY = XYZ.redY;
d = XYZ.greenX + XYZ.greenY + XYZ.greenZ;
if (!png_muldiv(&xy->greenx, XYZ.greenX, PNG_FP_1, d)) return 1;
if (!png_muldiv(&xy->greeny, XYZ.greenY, PNG_FP_1, d)) return 1;
dwhite += d;
whiteX += XYZ.greenX;
whiteY += XYZ.greenY;
d = XYZ.blueX + XYZ.blueY + XYZ.blueZ;
if (!png_muldiv(&xy->bluex, XYZ.blueX, PNG_FP_1, d)) return 1;
if (!png_muldiv(&xy->bluey, XYZ.blueY, PNG_FP_1, d)) return 1;
dwhite += d;
whiteX += XYZ.blueX;
whiteY += XYZ.blueY;
/* The reference white is simply the same of the end-point (X,Y,Z) vectors,
* thus:
*/
if (!png_muldiv(&xy->whitex, whiteX, PNG_FP_1, dwhite)) return 1;
if (!png_muldiv(&xy->whitey, whiteY, PNG_FP_1, dwhite)) return 1;
return 0;
}
int png_XYZ_from_xy(png_XYZ *XYZ, 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.)
*/
if (xy.redx < 0 || xy.redx > PNG_FP_1) return 1;
if (xy.redy < 0 || xy.redy > PNG_FP_1-xy.redx) return 1;
if (xy.greenx < 0 || xy.greenx > PNG_FP_1) return 1;
if (xy.greeny < 0 || xy.greeny > PNG_FP_1-xy.greenx) return 1;
if (xy.bluex < 0 || xy.bluex > PNG_FP_1) return 1;
if (xy.bluey < 0 || xy.bluey > PNG_FP_1-xy.bluex) return 1;
if (xy.whitex < 0 || xy.whitex > PNG_FP_1) return 1;
if (xy.whitey < 0 || xy.whitey > PNG_FP_1-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;
* (red,green,blue,white)x(x,y). This loses one degree of freedom and
* therefore an arbitrary ninth value has to be introduced to undo the
* original transformations.
*
* Think of the original end-points as points in (X,Y,Z) space. The
* chromaticity values (c) have the property:
*
* C
* c = ---------
* X + Y + Z
*
* For each c (x,y,z) from the corresponding original C (X,Y,Z). Thus the
* three chromaticity values (x,y,z) for each end-point obey the
* relationship:
*
* x + y + z = 1
*
* This describes the plane in (X,Y,Z) space that intersects each axis at the
* value 1.0; call this the chromaticity plane. Thus the chromaticity
* calculation has scaled each end-point so that it is on the x+y+z=1 plane
* and chromaticity is the intersection of the vector from the origin to the
* (X,Y,Z) value with the chromaticity plane.
*
* To fully invert the chromaticity calculation we would need the three
* end-point scale factors, (red-scale, green-scale, blue-scale), but these
* were not recorded. Instead we calculated the reference white (X,Y,Z) and
* recorded the chromaticity of this. The reference white (X,Y,Z) would have
* given all three of the scale factors since:
*
* color-C = color-c * color-scale
* white-C = red-C + green-C + blue-C
* = red-c*red-scale + green-c*green-scale + blue-c*blue-scale
*
* But cHRM records only white-x and white-y, so we have lost the white scale
* factor:
*
* white-C = white-c*white-scale
*
* To handle this the inverse transformation makes an arbitrary assumption
* about white-scale:
*
* Assume: white-Y = 1.0
* Hence: white-scale = 1/white-y
* Or: red-Y + green-Y + blue-Y = 1.0
*
* Notice the last statement of the assumption gives an equation in three of
* the nine values we want to calculate. 8 more equations come from the
* above routine as summarised at the top above (the chromaticity
* calculation):
*
* Given: color-x = color-X / (color-X + color-Y + color-Z)
* Hence: (color-x - 1)*color-X + color.x*color-Y + color.x*color-Z = 0
*
* This is 9 simultaneous equations in the 9 variables "color-C" and can be
* solved by Cramer's rule. Cramer's rule requires calculating 10 9x9 matrix
* determinants, however this is not as bad as it seems because only 28 of
* the total of 90 terms in the various matrices are non-zero. Nevertheless
* Cramer's rule is notoriously numerically unstable because the determinant
* calculation involves the difference of large, but similar, numbers. It is
* difficult to be sure that the calculation is stable for real world values
* and it is certain that it becomes unstable where the end points are close
* together.
*
* So this code uses the perhaps slighly less optimal but more understandable
* and totally obvious approach of calculating color-scale.
*
* This algorithm depends on the precision in white-scale and that is
* (1/white-y), so we can immediately see that as white-y approaches 0 the
* accuracy inherent in the cHRM chunk drops off substantially.
*
* libpng arithmetic: a simple invertion of the above equations
* ------------------------------------------------------------
*
* white_scale = 1/white-y
* white-X = white-x * white-scale
* white-Y = 1.0
* white-Z = (1 - white-x - white-y) * white_scale
*
* white-C = red-C + green-C + blue-C
* = red-c*red-scale + green-c*green-scale + blue-c*blue-scale
*
* This gives us three equations in (red-scale,green-scale,blue-scale) where
* all the coefficients are now known:
*
* red-x*red-scale + green-x*green-scale + blue-x*blue-scale
* = white-x/white-y
* red-y*red-scale + green-y*green-scale + blue-y*blue-scale = 1
* red-z*red-scale + green-z*green-scale + blue-z*blue-scale
* = (1 - white-x - white-y)/white-y
*
* In the last equation color-z is (1 - color-x - color-y) so we can add all
* three equations together to get an alternative third:
*
* red-scale + green-scale + blue-scale = 1/white-y = white-scale
*
* So now we have a Cramer's rule solution where the determinants are just
* 3x3 - far more tractible. Unfortunately 3x3 determinants still involve
* multiplication of three coefficients so we can't guarantee to avoid
* overflow in the libpng fixed point representation. Using Cramer's rule in
* floating point is probably a good choice here, but it's not an option for
* fixed point. Instead proceed to simplify the first two equations by
* eliminating what is likely to be the largest value, blue-scale:
*
* blue-scale = white-scale - red-scale - green-scale
*
* Hence:
*
* (red-x - blue-x)*red-scale + (green-x - blue-x)*green-scale =
* (white-x - blue-x)*white-scale
*
* (red-y - blue-y)*red-scale + (green-y - blue-y)*green-scale =
* 1 - blue-y*white-scale
*
* And now we can trivially solve for (red-scale,green-scale):
*
* green-scale =
* (white-x - blue-x)*white-scale - (red-x - blue-x)*red-scale
* -----------------------------------------------------------
* green-x - blue-x
*
* red-scale =
* 1 - blue-y*white-scale - (green-y - blue-y) * green-scale
* ---------------------------------------------------------
* red-y - blue-y
*
* Hence:
*
* red-scale =
* ( (green-x - blue-x) * (white-y - blue-y) -
* (green-y - blue-y) * (white-x - blue-x) ) / white-y
* -------------------------------------------------------------------------
* (green-x - blue-x)*(red-y - blue-y)-(green-y - blue-y)*(red-x - blue-x)
*
* green-scale =
* ( (red-y - blue-y) * (white-x - blue-x) -
* (red-x - blue-x) * (white-y - blue-y) ) / white-y
* -------------------------------------------------------------------------
* (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 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.
*
* Note that the values of the differences of the products of the
* chromaticities in the above equations tend to be small, for example for
* the sRGB chromaticities they are:
*
* red numerator: -0.04751
* green numerator: -0.08788
* denominator: -0.2241 (without white-y multiplication)
*
* The resultant Y coefficients from the chromaticities of some widely used
* color space definitions are (to 15 decimal places):
*
* sRGB
* 0.212639005871510 0.715168678767756 0.072192315360734
* Kodak ProPhoto
* 0.288071128229293 0.711843217810102 0.000085653960605
* Adobe RGB
* 0.297344975250536 0.627363566255466 0.075291458493998
* 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.
*/
if (!png_muldiv(&left, xy.greenx-xy.bluex, xy.redy - xy.bluey, 7)) return 2;
if (!png_muldiv(&right, xy.greeny-xy.bluey, xy.redx - xy.bluex, 7)) return 2;
denominator = left - right;
/* Now find the red numerator. */
if (!png_muldiv(&left, xy.greenx-xy.bluex, xy.whitey-xy.bluey, 7)) return 2;
if (!png_muldiv(&right, xy.greeny-xy.bluey, xy.whitex-xy.bluex, 7)) return 2;
/* 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) ||
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)) return 2;
if (!png_muldiv(&right, xy.redx-xy.bluex, xy.whitey-xy.bluey, 7)) return 2;
if (!png_muldiv(&green_inverse, xy.whitey, denominator, left-right) ||
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) return 1;
/* And fill in the png_XYZ: */
if (!png_muldiv(&XYZ->redX, xy.redx, PNG_FP_1, red_inverse)) return 1;
if (!png_muldiv(&XYZ->redY, xy.redy, PNG_FP_1, red_inverse)) return 1;
if (!png_muldiv(&XYZ->redZ, PNG_FP_1 - xy.redx - xy.redy, PNG_FP_1,
red_inverse))
return 1;
if (!png_muldiv(&XYZ->greenX, xy.greenx, PNG_FP_1, green_inverse)) return 1;
if (!png_muldiv(&XYZ->greenY, xy.greeny, PNG_FP_1, green_inverse)) return 1;
if (!png_muldiv(&XYZ->greenZ, PNG_FP_1 - xy.greenx - xy.greeny, PNG_FP_1,
green_inverse))
return 1;
if (!png_muldiv(&XYZ->blueX, xy.bluex, blue_scale, PNG_FP_1)) return 1;
if (!png_muldiv(&XYZ->blueY, xy.bluey, blue_scale, PNG_FP_1)) return 1;
if (!png_muldiv(&XYZ->blueZ, PNG_FP_1 - xy.bluex - xy.bluey, blue_scale,
PNG_FP_1))
return 1;
return 0; /*success*/
}
int png_XYZ_from_xy_checked(png_structp png_ptr, png_XYZ *XYZ, png_xy xy)
{
switch (png_XYZ_from_xy(XYZ, xy))
{
case 0: /* success */
return 1;
case 1:
/* The chunk may be technically valid, but we got png_fixed_point
* overflow while trying to get XYZ values out of it. This is
* entirely benign - the cHRM chunk is pretty extreme.
*/
png_chunk_benign_error(png_ptr,
"extreme cHRM chunk cannot be converted to tristimulus values");
break;
default:
/* libpng is broken; this should be a warning but if it happens we
* want error reports so for the moment it is an error.
*/
png_error(png_ptr, "internal error in png_XYZ_from_xy");
break;
}
/* ERROR RETURN */
return 0;
}
#endif
void /* PRIVATE */
png_check_IHDR(png_structp png_ptr,
png_uint_32 width, png_uint_32 height, int bit_depth,