2 /* png.c - location for general purpose libpng functions
4 * Last changed in libpng 1.5.11 [June 14, 2012]
5 * Copyright (c) 1998-2012 Glenn Randers-Pehrson
6 * (Version 0.96 Copyright (c) 1996, 1997 Andreas Dilger)
7 * (Version 0.88 Copyright (c) 1995, 1996 Guy Eric Schalnat, Group 42, Inc.)
9 * This code is released under the libpng license.
10 * For conditions of distribution and use, see the disclaimer
11 * and license in png.h
16 /* Generate a compiler error if there is an old png.h in the search path. */
17 typedef png_libpng_version_1_5_12 Your_png_h_is_not_version_1_5_12;
19 /* Tells libpng that we have already handled the first "num_bytes" bytes
20 * of the PNG file signature. If the PNG data is embedded into another
21 * stream we can set num_bytes = 8 so that libpng will not attempt to read
22 * or write any of the magic bytes before it starts on the IHDR.
25 #ifdef PNG_READ_SUPPORTED
27 png_set_sig_bytes(png_structp png_ptr, int num_bytes)
29 png_debug(1, "in png_set_sig_bytes");
35 png_error(png_ptr, "Too many bytes for PNG signature");
37 png_ptr->sig_bytes = (png_byte)(num_bytes < 0 ? 0 : num_bytes);
40 /* Checks whether the supplied bytes match the PNG signature. We allow
41 * checking less than the full 8-byte signature so that those apps that
42 * already read the first few bytes of a file to determine the file type
43 * can simply check the remaining bytes for extra assurance. Returns
44 * an integer less than, equal to, or greater than zero if sig is found,
45 * respectively, to be less than, to match, or be greater than the correct
46 * PNG signature (this is the same behavior as strcmp, memcmp, etc).
49 png_sig_cmp(png_const_bytep sig, png_size_t start, png_size_t num_to_check)
51 png_byte png_signature[8] = {137, 80, 78, 71, 13, 10, 26, 10};
56 else if (num_to_check < 1)
62 if (start + num_to_check > 8)
63 num_to_check = 8 - start;
65 return ((int)(png_memcmp(&sig[start], &png_signature[start], num_to_check)));
68 #endif /* PNG_READ_SUPPORTED */
70 #if defined(PNG_READ_SUPPORTED) || defined(PNG_WRITE_SUPPORTED)
71 /* Function to allocate memory for zlib */
72 PNG_FUNCTION(voidpf /* PRIVATE */,
73 png_zalloc,(voidpf png_ptr, uInt items, uInt size),PNG_ALLOCATED)
76 png_structp p=(png_structp)png_ptr;
77 png_uint_32 save_flags=p->flags;
78 png_alloc_size_t num_bytes;
83 if (items > PNG_UINT_32_MAX/size)
85 png_warning (p, "Potential overflow in png_zalloc()");
88 num_bytes = (png_alloc_size_t)items * size;
90 p->flags|=PNG_FLAG_MALLOC_NULL_MEM_OK;
91 ptr = (png_voidp)png_malloc((png_structp)png_ptr, num_bytes);
97 /* Function to free memory for zlib */
99 png_zfree(voidpf png_ptr, voidpf ptr)
101 png_free((png_structp)png_ptr, (png_voidp)ptr);
104 /* Reset the CRC variable to 32 bits of 1's. Care must be taken
105 * in case CRC is > 32 bits to leave the top bits 0.
108 png_reset_crc(png_structp png_ptr)
110 /* The cast is safe because the crc is a 32 bit value. */
111 png_ptr->crc = (png_uint_32)crc32(0, Z_NULL, 0);
114 /* Calculate the CRC over a section of data. We can only pass as
115 * much data to this routine as the largest single buffer size. We
116 * also check that this data will actually be used before going to the
117 * trouble of calculating it.
120 png_calculate_crc(png_structp png_ptr, png_const_bytep ptr, png_size_t length)
124 if (PNG_CHUNK_ANCILLIARY(png_ptr->chunk_name))
126 if ((png_ptr->flags & PNG_FLAG_CRC_ANCILLARY_MASK) ==
127 (PNG_FLAG_CRC_ANCILLARY_USE | PNG_FLAG_CRC_ANCILLARY_NOWARN))
133 if (png_ptr->flags & PNG_FLAG_CRC_CRITICAL_IGNORE)
137 /* 'uLong' is defined as unsigned long, this means that on some systems it is
138 * a 64 bit value. crc32, however, returns 32 bits so the following cast is
139 * safe. 'uInt' may be no more than 16 bits, so it is necessary to perform a
142 if (need_crc && length > 0)
144 uLong crc = png_ptr->crc; /* Should never issue a warning */
148 uInt safeLength = (uInt)length;
150 safeLength = (uInt)-1; /* evil, but safe */
152 crc = crc32(crc, ptr, safeLength);
154 /* The following should never issue compiler warnings, if they do the
155 * target system has characteristics that will probably violate other
156 * assumptions within the libpng code.
159 length -= safeLength;
163 /* And the following is always safe because the crc is only 32 bits. */
164 png_ptr->crc = (png_uint_32)crc;
168 /* Check a user supplied version number, called from both read and write
169 * functions that create a png_struct
172 png_user_version_check(png_structp png_ptr, png_const_charp user_png_ver)
180 if (user_png_ver[i] != png_libpng_ver[i])
181 png_ptr->flags |= PNG_FLAG_LIBRARY_MISMATCH;
182 } while (png_libpng_ver[i++]);
186 png_ptr->flags |= PNG_FLAG_LIBRARY_MISMATCH;
188 if (png_ptr->flags & PNG_FLAG_LIBRARY_MISMATCH)
190 /* Libpng 0.90 and later are binary incompatible with libpng 0.89, so
191 * we must recompile any applications that use any older library version.
192 * For versions after libpng 1.0, we will be compatible, so we need
193 * only check the first digit.
195 if (user_png_ver == NULL || user_png_ver[0] != png_libpng_ver[0] ||
196 (user_png_ver[0] == '1' && user_png_ver[2] != png_libpng_ver[2]) ||
197 (user_png_ver[0] == '0' && user_png_ver[2] < '9'))
199 #ifdef PNG_WARNINGS_SUPPORTED
203 pos = png_safecat(m, sizeof m, pos, "Application built with libpng-");
204 pos = png_safecat(m, sizeof m, pos, user_png_ver);
205 pos = png_safecat(m, sizeof m, pos, " but running with ");
206 pos = png_safecat(m, sizeof m, pos, png_libpng_ver);
208 png_warning(png_ptr, m);
211 #ifdef PNG_ERROR_NUMBERS_SUPPORTED
219 /* Success return. */
223 /* Allocate the memory for an info_struct for the application. We don't
224 * really need the png_ptr, but it could potentially be useful in the
225 * future. This should be used in favour of malloc(png_sizeof(png_info))
226 * and png_info_init() so that applications that want to use a shared
227 * libpng don't have to be recompiled if png_info changes size.
229 PNG_FUNCTION(png_infop,PNGAPI
230 png_create_info_struct,(png_structp png_ptr),PNG_ALLOCATED)
234 png_debug(1, "in png_create_info_struct");
239 #ifdef PNG_USER_MEM_SUPPORTED
240 info_ptr = (png_infop)png_create_struct_2(PNG_STRUCT_INFO,
241 png_ptr->malloc_fn, png_ptr->mem_ptr);
243 info_ptr = (png_infop)png_create_struct(PNG_STRUCT_INFO);
245 if (info_ptr != NULL)
246 png_info_init_3(&info_ptr, png_sizeof(png_info));
251 /* This function frees the memory associated with a single info struct.
252 * Normally, one would use either png_destroy_read_struct() or
253 * png_destroy_write_struct() to free an info struct, but this may be
254 * useful for some applications.
257 png_destroy_info_struct(png_structp png_ptr, png_infopp info_ptr_ptr)
259 png_infop info_ptr = NULL;
261 png_debug(1, "in png_destroy_info_struct");
266 if (info_ptr_ptr != NULL)
267 info_ptr = *info_ptr_ptr;
269 if (info_ptr != NULL)
271 png_info_destroy(png_ptr, info_ptr);
273 #ifdef PNG_USER_MEM_SUPPORTED
274 png_destroy_struct_2((png_voidp)info_ptr, png_ptr->free_fn,
277 png_destroy_struct((png_voidp)info_ptr);
279 *info_ptr_ptr = NULL;
283 /* Initialize the info structure. This is now an internal function (0.89)
284 * and applications using it are urged to use png_create_info_struct()
289 png_info_init_3(png_infopp ptr_ptr, png_size_t png_info_struct_size)
291 png_infop info_ptr = *ptr_ptr;
293 png_debug(1, "in png_info_init_3");
295 if (info_ptr == NULL)
298 if (png_sizeof(png_info) > png_info_struct_size)
300 png_destroy_struct(info_ptr);
301 info_ptr = (png_infop)png_create_struct(PNG_STRUCT_INFO);
305 /* Set everything to 0 */
306 png_memset(info_ptr, 0, png_sizeof(png_info));
310 png_data_freer(png_structp png_ptr, png_infop info_ptr,
311 int freer, png_uint_32 mask)
313 png_debug(1, "in png_data_freer");
315 if (png_ptr == NULL || info_ptr == NULL)
318 if (freer == PNG_DESTROY_WILL_FREE_DATA)
319 info_ptr->free_me |= mask;
321 else if (freer == PNG_USER_WILL_FREE_DATA)
322 info_ptr->free_me &= ~mask;
326 "Unknown freer parameter in png_data_freer");
330 png_free_data(png_structp png_ptr, png_infop info_ptr, png_uint_32 mask,
333 png_debug(1, "in png_free_data");
335 if (png_ptr == NULL || info_ptr == NULL)
338 #ifdef PNG_TEXT_SUPPORTED
339 /* Free text item num or (if num == -1) all text items */
340 if ((mask & PNG_FREE_TEXT) & info_ptr->free_me)
344 if (info_ptr->text && info_ptr->text[num].key)
346 png_free(png_ptr, info_ptr->text[num].key);
347 info_ptr->text[num].key = NULL;
354 for (i = 0; i < info_ptr->num_text; i++)
355 png_free_data(png_ptr, info_ptr, PNG_FREE_TEXT, i);
356 png_free(png_ptr, info_ptr->text);
357 info_ptr->text = NULL;
358 info_ptr->num_text=0;
363 #ifdef PNG_tRNS_SUPPORTED
364 /* Free any tRNS entry */
365 if ((mask & PNG_FREE_TRNS) & info_ptr->free_me)
367 png_free(png_ptr, info_ptr->trans_alpha);
368 info_ptr->trans_alpha = NULL;
369 info_ptr->valid &= ~PNG_INFO_tRNS;
373 #ifdef PNG_sCAL_SUPPORTED
374 /* Free any sCAL entry */
375 if ((mask & PNG_FREE_SCAL) & info_ptr->free_me)
377 png_free(png_ptr, info_ptr->scal_s_width);
378 png_free(png_ptr, info_ptr->scal_s_height);
379 info_ptr->scal_s_width = NULL;
380 info_ptr->scal_s_height = NULL;
381 info_ptr->valid &= ~PNG_INFO_sCAL;
385 #ifdef PNG_pCAL_SUPPORTED
386 /* Free any pCAL entry */
387 if ((mask & PNG_FREE_PCAL) & info_ptr->free_me)
389 png_free(png_ptr, info_ptr->pcal_purpose);
390 png_free(png_ptr, info_ptr->pcal_units);
391 info_ptr->pcal_purpose = NULL;
392 info_ptr->pcal_units = NULL;
393 if (info_ptr->pcal_params != NULL)
396 for (i = 0; i < (int)info_ptr->pcal_nparams; i++)
398 png_free(png_ptr, info_ptr->pcal_params[i]);
399 info_ptr->pcal_params[i] = NULL;
401 png_free(png_ptr, info_ptr->pcal_params);
402 info_ptr->pcal_params = NULL;
404 info_ptr->valid &= ~PNG_INFO_pCAL;
408 #ifdef PNG_iCCP_SUPPORTED
409 /* Free any iCCP entry */
410 if ((mask & PNG_FREE_ICCP) & info_ptr->free_me)
412 png_free(png_ptr, info_ptr->iccp_name);
413 png_free(png_ptr, info_ptr->iccp_profile);
414 info_ptr->iccp_name = NULL;
415 info_ptr->iccp_profile = NULL;
416 info_ptr->valid &= ~PNG_INFO_iCCP;
420 #ifdef PNG_sPLT_SUPPORTED
421 /* Free a given sPLT entry, or (if num == -1) all sPLT entries */
422 if ((mask & PNG_FREE_SPLT) & info_ptr->free_me)
426 if (info_ptr->splt_palettes)
428 png_free(png_ptr, info_ptr->splt_palettes[num].name);
429 png_free(png_ptr, info_ptr->splt_palettes[num].entries);
430 info_ptr->splt_palettes[num].name = NULL;
431 info_ptr->splt_palettes[num].entries = NULL;
437 if (info_ptr->splt_palettes_num)
440 for (i = 0; i < (int)info_ptr->splt_palettes_num; i++)
441 png_free_data(png_ptr, info_ptr, PNG_FREE_SPLT, i);
443 png_free(png_ptr, info_ptr->splt_palettes);
444 info_ptr->splt_palettes = NULL;
445 info_ptr->splt_palettes_num = 0;
447 info_ptr->valid &= ~PNG_INFO_sPLT;
452 #ifdef PNG_UNKNOWN_CHUNKS_SUPPORTED
453 if (png_ptr->unknown_chunk.data)
455 png_free(png_ptr, png_ptr->unknown_chunk.data);
456 png_ptr->unknown_chunk.data = NULL;
459 if ((mask & PNG_FREE_UNKN) & info_ptr->free_me)
463 if (info_ptr->unknown_chunks)
465 png_free(png_ptr, info_ptr->unknown_chunks[num].data);
466 info_ptr->unknown_chunks[num].data = NULL;
474 if (info_ptr->unknown_chunks_num)
476 for (i = 0; i < info_ptr->unknown_chunks_num; i++)
477 png_free_data(png_ptr, info_ptr, PNG_FREE_UNKN, i);
479 png_free(png_ptr, info_ptr->unknown_chunks);
480 info_ptr->unknown_chunks = NULL;
481 info_ptr->unknown_chunks_num = 0;
487 #ifdef PNG_hIST_SUPPORTED
488 /* Free any hIST entry */
489 if ((mask & PNG_FREE_HIST) & info_ptr->free_me)
491 png_free(png_ptr, info_ptr->hist);
492 info_ptr->hist = NULL;
493 info_ptr->valid &= ~PNG_INFO_hIST;
497 /* Free any PLTE entry that was internally allocated */
498 if ((mask & PNG_FREE_PLTE) & info_ptr->free_me)
500 png_zfree(png_ptr, info_ptr->palette);
501 info_ptr->palette = NULL;
502 info_ptr->valid &= ~PNG_INFO_PLTE;
503 info_ptr->num_palette = 0;
506 #ifdef PNG_INFO_IMAGE_SUPPORTED
507 /* Free any image bits attached to the info structure */
508 if ((mask & PNG_FREE_ROWS) & info_ptr->free_me)
510 if (info_ptr->row_pointers)
513 for (row = 0; row < (int)info_ptr->height; row++)
515 png_free(png_ptr, info_ptr->row_pointers[row]);
516 info_ptr->row_pointers[row] = NULL;
518 png_free(png_ptr, info_ptr->row_pointers);
519 info_ptr->row_pointers = NULL;
521 info_ptr->valid &= ~PNG_INFO_IDAT;
526 mask &= ~PNG_FREE_MUL;
528 info_ptr->free_me &= ~mask;
531 /* This is an internal routine to free any memory that the info struct is
532 * pointing to before re-using it or freeing the struct itself. Recall
533 * that png_free() checks for NULL pointers for us.
536 png_info_destroy(png_structp png_ptr, png_infop info_ptr)
538 png_debug(1, "in png_info_destroy");
540 png_free_data(png_ptr, info_ptr, PNG_FREE_ALL, -1);
542 #ifdef PNG_HANDLE_AS_UNKNOWN_SUPPORTED
543 if (png_ptr->num_chunk_list)
545 png_free(png_ptr, png_ptr->chunk_list);
546 png_ptr->chunk_list = NULL;
547 png_ptr->num_chunk_list = 0;
551 png_info_init_3(&info_ptr, png_sizeof(png_info));
553 #endif /* defined(PNG_READ_SUPPORTED) || defined(PNG_WRITE_SUPPORTED) */
555 /* This function returns a pointer to the io_ptr associated with the user
556 * functions. The application should free any memory associated with this
557 * pointer before png_write_destroy() or png_read_destroy() are called.
560 png_get_io_ptr(png_structp png_ptr)
565 return (png_ptr->io_ptr);
568 #if defined(PNG_READ_SUPPORTED) || defined(PNG_WRITE_SUPPORTED)
569 # ifdef PNG_STDIO_SUPPORTED
570 /* Initialize the default input/output functions for the PNG file. If you
571 * use your own read or write routines, you can call either png_set_read_fn()
572 * or png_set_write_fn() instead of png_init_io(). If you have defined
573 * PNG_NO_STDIO or otherwise disabled PNG_STDIO_SUPPORTED, you must use a
574 * function of your own because "FILE *" isn't necessarily available.
577 png_init_io(png_structp png_ptr, png_FILE_p fp)
579 png_debug(1, "in png_init_io");
584 png_ptr->io_ptr = (png_voidp)fp;
588 # ifdef PNG_TIME_RFC1123_SUPPORTED
589 /* Convert the supplied time into an RFC 1123 string suitable for use in
590 * a "Creation Time" or other text-based time string.
592 png_const_charp PNGAPI
593 png_convert_to_rfc1123(png_structp png_ptr, png_const_timep ptime)
595 static PNG_CONST char short_months[12][4] =
596 {"Jan", "Feb", "Mar", "Apr", "May", "Jun",
597 "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"};
602 if (ptime->year > 9999 /* RFC1123 limitation */ ||
603 ptime->month == 0 || ptime->month > 12 ||
604 ptime->day == 0 || ptime->day > 31 ||
605 ptime->hour > 23 || ptime->minute > 59 ||
608 png_warning(png_ptr, "Ignoring invalid time value");
614 char number_buf[5]; /* enough for a four-digit year */
616 # define APPEND_STRING(string)\
617 pos = png_safecat(png_ptr->time_buffer, sizeof png_ptr->time_buffer,\
619 # define APPEND_NUMBER(format, value)\
620 APPEND_STRING(PNG_FORMAT_NUMBER(number_buf, format, (value)))
622 if (pos < (sizeof png_ptr->time_buffer)-1)\
623 png_ptr->time_buffer[pos++] = (ch)
625 APPEND_NUMBER(PNG_NUMBER_FORMAT_u, (unsigned)ptime->day);
627 APPEND_STRING(short_months[(ptime->month - 1)]);
629 APPEND_NUMBER(PNG_NUMBER_FORMAT_u, ptime->year);
631 APPEND_NUMBER(PNG_NUMBER_FORMAT_02u, (unsigned)ptime->hour);
633 APPEND_NUMBER(PNG_NUMBER_FORMAT_02u, (unsigned)ptime->minute);
635 APPEND_NUMBER(PNG_NUMBER_FORMAT_02u, (unsigned)ptime->second);
636 APPEND_STRING(" +0000"); /* This reliably terminates the buffer */
639 # undef APPEND_NUMBER
640 # undef APPEND_STRING
643 return png_ptr->time_buffer;
645 # endif /* PNG_TIME_RFC1123_SUPPORTED */
647 #endif /* defined(PNG_READ_SUPPORTED) || defined(PNG_WRITE_SUPPORTED) */
649 png_const_charp PNGAPI
650 png_get_copyright(png_const_structp png_ptr)
652 PNG_UNUSED(png_ptr) /* Silence compiler warning about unused png_ptr */
653 #ifdef PNG_STRING_COPYRIGHT
654 return PNG_STRING_COPYRIGHT
657 return PNG_STRING_NEWLINE \
658 "libpng version 1.5.12 - July 11, 2012" PNG_STRING_NEWLINE \
659 "Copyright (c) 1998-2012 Glenn Randers-Pehrson" PNG_STRING_NEWLINE \
660 "Copyright (c) 1996-1997 Andreas Dilger" PNG_STRING_NEWLINE \
661 "Copyright (c) 1995-1996 Guy Eric Schalnat, Group 42, Inc." \
664 return "libpng version 1.5.12 - July 11, 2012\
665 Copyright (c) 1998-2012 Glenn Randers-Pehrson\
666 Copyright (c) 1996-1997 Andreas Dilger\
667 Copyright (c) 1995-1996 Guy Eric Schalnat, Group 42, Inc.";
672 /* The following return the library version as a short string in the
673 * format 1.0.0 through 99.99.99zz. To get the version of *.h files
674 * used with your application, print out PNG_LIBPNG_VER_STRING, which
675 * is defined in png.h.
676 * Note: now there is no difference between png_get_libpng_ver() and
677 * png_get_header_ver(). Due to the version_nn_nn_nn typedef guard,
678 * it is guaranteed that png.c uses the correct version of png.h.
680 png_const_charp PNGAPI
681 png_get_libpng_ver(png_const_structp png_ptr)
683 /* Version of *.c files used when building libpng */
684 return png_get_header_ver(png_ptr);
687 png_const_charp PNGAPI
688 png_get_header_ver(png_const_structp png_ptr)
690 /* Version of *.h files used when building libpng */
691 PNG_UNUSED(png_ptr) /* Silence compiler warning about unused png_ptr */
692 return PNG_LIBPNG_VER_STRING;
695 png_const_charp PNGAPI
696 png_get_header_version(png_const_structp png_ptr)
698 /* Returns longer string containing both version and date */
699 PNG_UNUSED(png_ptr) /* Silence compiler warning about unused png_ptr */
701 return PNG_HEADER_VERSION_STRING
702 # ifndef PNG_READ_SUPPORTED
707 return PNG_HEADER_VERSION_STRING;
711 #ifdef PNG_HANDLE_AS_UNKNOWN_SUPPORTED
713 png_handle_as_unknown(png_structp png_ptr, png_const_bytep chunk_name)
715 /* Check chunk_name and return "keep" value if it's on the list, else 0 */
716 png_const_bytep p, p_end;
718 if (png_ptr == NULL || chunk_name == NULL || png_ptr->num_chunk_list <= 0)
719 return PNG_HANDLE_CHUNK_AS_DEFAULT;
721 p_end = png_ptr->chunk_list;
722 p = p_end + png_ptr->num_chunk_list*5; /* beyond end */
724 /* The code is the fifth byte after each four byte string. Historically this
725 * code was always searched from the end of the list, so it should continue
726 * to do so in case there are duplicated entries.
728 do /* num_chunk_list > 0, so at least one */
731 if (!png_memcmp(chunk_name, p, 4))
736 return PNG_HANDLE_CHUNK_AS_DEFAULT;
740 png_chunk_unknown_handling(png_structp png_ptr, png_uint_32 chunk_name)
742 png_byte chunk_string[5];
744 PNG_CSTRING_FROM_CHUNK(chunk_string, chunk_name);
745 return png_handle_as_unknown(png_ptr, chunk_string);
749 #ifdef PNG_READ_SUPPORTED
750 /* This function, added to libpng-1.0.6g, is untested. */
752 png_reset_zstream(png_structp png_ptr)
755 return Z_STREAM_ERROR;
757 return (inflateReset(&png_ptr->zstream));
759 #endif /* PNG_READ_SUPPORTED */
761 /* This function was added to libpng-1.0.7 */
763 png_access_version_number(void)
765 /* Version of *.c files used when building libpng */
766 return((png_uint_32)PNG_LIBPNG_VER);
771 #if defined(PNG_READ_SUPPORTED) || defined(PNG_WRITE_SUPPORTED)
772 /* png_convert_size: a PNGAPI but no longer in png.h, so deleted
776 /* Added at libpng version 1.2.34 and 1.4.0 (moved from pngset.c) */
777 # ifdef PNG_CHECK_cHRM_SUPPORTED
780 png_check_cHRM_fixed(png_structp png_ptr,
781 png_fixed_point white_x, png_fixed_point white_y, png_fixed_point red_x,
782 png_fixed_point red_y, png_fixed_point green_x, png_fixed_point green_y,
783 png_fixed_point blue_x, png_fixed_point blue_y)
786 unsigned long xy_hi,xy_lo,yx_hi,yx_lo;
788 png_debug(1, "in function png_check_cHRM_fixed");
793 /* (x,y,z) values are first limited to 0..100000 (PNG_FP_1), the white
794 * y must also be greater than 0. To test for the upper limit calculate
795 * (PNG_FP_1-y) - x must be <= to this for z to be >= 0 (and the expression
796 * cannot overflow.) At this point we know x and y are >= 0 and (x+y) is
797 * <= PNG_FP_1. The previous test on PNG_MAX_UINT_31 is removed because it
798 * pointless (and it produces compiler warnings!)
800 if (white_x < 0 || white_y <= 0 ||
801 red_x < 0 || red_y < 0 ||
802 green_x < 0 || green_y < 0 ||
803 blue_x < 0 || blue_y < 0)
806 "Ignoring attempt to set negative chromaticity value");
809 /* And (x+y) must be <= PNG_FP_1 (so z is >= 0) */
810 if (white_x > PNG_FP_1 - white_y)
812 png_warning(png_ptr, "Invalid cHRM white point");
816 if (red_x > PNG_FP_1 - red_y)
818 png_warning(png_ptr, "Invalid cHRM red point");
822 if (green_x > PNG_FP_1 - green_y)
824 png_warning(png_ptr, "Invalid cHRM green point");
828 if (blue_x > PNG_FP_1 - blue_y)
830 png_warning(png_ptr, "Invalid cHRM blue point");
834 png_64bit_product(green_x - red_x, blue_y - red_y, &xy_hi, &xy_lo);
835 png_64bit_product(green_y - red_y, blue_x - red_x, &yx_hi, &yx_lo);
837 if (xy_hi == yx_hi && xy_lo == yx_lo)
840 "Ignoring attempt to set cHRM RGB triangle with zero area");
846 # endif /* PNG_CHECK_cHRM_SUPPORTED */
848 #ifdef PNG_cHRM_SUPPORTED
849 /* Added at libpng-1.5.5 to support read and write of true CIEXYZ values for
850 * cHRM, as opposed to using chromaticities. These internal APIs return
851 * non-zero on a parameter error. The X, Y and Z values are required to be
852 * positive and less than 1.0.
854 int png_xy_from_XYZ(png_xy *xy, png_XYZ XYZ)
856 png_int_32 d, dwhite, whiteX, whiteY;
858 d = XYZ.redX + XYZ.redY + XYZ.redZ;
859 if (!png_muldiv(&xy->redx, XYZ.redX, PNG_FP_1, d)) return 1;
860 if (!png_muldiv(&xy->redy, XYZ.redY, PNG_FP_1, d)) return 1;
865 d = XYZ.greenX + XYZ.greenY + XYZ.greenZ;
866 if (!png_muldiv(&xy->greenx, XYZ.greenX, PNG_FP_1, d)) return 1;
867 if (!png_muldiv(&xy->greeny, XYZ.greenY, PNG_FP_1, d)) return 1;
869 whiteX += XYZ.greenX;
870 whiteY += XYZ.greenY;
872 d = XYZ.blueX + XYZ.blueY + XYZ.blueZ;
873 if (!png_muldiv(&xy->bluex, XYZ.blueX, PNG_FP_1, d)) return 1;
874 if (!png_muldiv(&xy->bluey, XYZ.blueY, PNG_FP_1, d)) return 1;
879 /* The reference white is simply the same of the end-point (X,Y,Z) vectors,
882 if (!png_muldiv(&xy->whitex, whiteX, PNG_FP_1, dwhite)) return 1;
883 if (!png_muldiv(&xy->whitey, whiteY, PNG_FP_1, dwhite)) return 1;
888 int png_XYZ_from_xy(png_XYZ *XYZ, png_xy xy)
890 png_fixed_point red_inverse, green_inverse, blue_scale;
891 png_fixed_point left, right, denominator;
893 /* Check xy and, implicitly, z. Note that wide gamut color spaces typically
894 * have end points with 0 tristimulus values (these are impossible end
895 * points, but they are used to cover the possible colors.)
897 if (xy.redx < 0 || xy.redx > PNG_FP_1) return 1;
898 if (xy.redy < 0 || xy.redy > PNG_FP_1-xy.redx) return 1;
899 if (xy.greenx < 0 || xy.greenx > PNG_FP_1) return 1;
900 if (xy.greeny < 0 || xy.greeny > PNG_FP_1-xy.greenx) return 1;
901 if (xy.bluex < 0 || xy.bluex > PNG_FP_1) return 1;
902 if (xy.bluey < 0 || xy.bluey > PNG_FP_1-xy.bluex) return 1;
903 if (xy.whitex < 0 || xy.whitex > PNG_FP_1) return 1;
904 if (xy.whitey < 0 || xy.whitey > PNG_FP_1-xy.whitex) return 1;
906 /* The reverse calculation is more difficult because the original tristimulus
907 * value had 9 independent values (red,green,blue)x(X,Y,Z) however only 8
908 * derived values were recorded in the cHRM chunk;
909 * (red,green,blue,white)x(x,y). This loses one degree of freedom and
910 * therefore an arbitrary ninth value has to be introduced to undo the
911 * original transformations.
913 * Think of the original end-points as points in (X,Y,Z) space. The
914 * chromaticity values (c) have the property:
920 * For each c (x,y,z) from the corresponding original C (X,Y,Z). Thus the
921 * three chromaticity values (x,y,z) for each end-point obey the
926 * This describes the plane in (X,Y,Z) space that intersects each axis at the
927 * value 1.0; call this the chromaticity plane. Thus the chromaticity
928 * calculation has scaled each end-point so that it is on the x+y+z=1 plane
929 * and chromaticity is the intersection of the vector from the origin to the
930 * (X,Y,Z) value with the chromaticity plane.
932 * To fully invert the chromaticity calculation we would need the three
933 * end-point scale factors, (red-scale, green-scale, blue-scale), but these
934 * were not recorded. Instead we calculated the reference white (X,Y,Z) and
935 * recorded the chromaticity of this. The reference white (X,Y,Z) would have
936 * given all three of the scale factors since:
938 * color-C = color-c * color-scale
939 * white-C = red-C + green-C + blue-C
940 * = red-c*red-scale + green-c*green-scale + blue-c*blue-scale
942 * But cHRM records only white-x and white-y, so we have lost the white scale
945 * white-C = white-c*white-scale
947 * To handle this the inverse transformation makes an arbitrary assumption
950 * Assume: white-Y = 1.0
951 * Hence: white-scale = 1/white-y
952 * Or: red-Y + green-Y + blue-Y = 1.0
954 * Notice the last statement of the assumption gives an equation in three of
955 * the nine values we want to calculate. 8 more equations come from the
956 * above routine as summarised at the top above (the chromaticity
959 * Given: color-x = color-X / (color-X + color-Y + color-Z)
960 * Hence: (color-x - 1)*color-X + color.x*color-Y + color.x*color-Z = 0
962 * This is 9 simultaneous equations in the 9 variables "color-C" and can be
963 * solved by Cramer's rule. Cramer's rule requires calculating 10 9x9 matrix
964 * determinants, however this is not as bad as it seems because only 28 of
965 * the total of 90 terms in the various matrices are non-zero. Nevertheless
966 * Cramer's rule is notoriously numerically unstable because the determinant
967 * calculation involves the difference of large, but similar, numbers. It is
968 * difficult to be sure that the calculation is stable for real world values
969 * and it is certain that it becomes unstable where the end points are close
972 * So this code uses the perhaps slightly less optimal but more
973 * understandable and totally obvious approach of calculating color-scale.
975 * This algorithm depends on the precision in white-scale and that is
976 * (1/white-y), so we can immediately see that as white-y approaches 0 the
977 * accuracy inherent in the cHRM chunk drops off substantially.
979 * libpng arithmetic: a simple invertion of the above equations
980 * ------------------------------------------------------------
982 * white_scale = 1/white-y
983 * white-X = white-x * white-scale
985 * white-Z = (1 - white-x - white-y) * white_scale
987 * white-C = red-C + green-C + blue-C
988 * = red-c*red-scale + green-c*green-scale + blue-c*blue-scale
990 * This gives us three equations in (red-scale,green-scale,blue-scale) where
991 * all the coefficients are now known:
993 * red-x*red-scale + green-x*green-scale + blue-x*blue-scale
995 * red-y*red-scale + green-y*green-scale + blue-y*blue-scale = 1
996 * red-z*red-scale + green-z*green-scale + blue-z*blue-scale
997 * = (1 - white-x - white-y)/white-y
999 * In the last equation color-z is (1 - color-x - color-y) so we can add all
1000 * three equations together to get an alternative third:
1002 * red-scale + green-scale + blue-scale = 1/white-y = white-scale
1004 * So now we have a Cramer's rule solution where the determinants are just
1005 * 3x3 - far more tractible. Unfortunately 3x3 determinants still involve
1006 * multiplication of three coefficients so we can't guarantee to avoid
1007 * overflow in the libpng fixed point representation. Using Cramer's rule in
1008 * floating point is probably a good choice here, but it's not an option for
1009 * fixed point. Instead proceed to simplify the first two equations by
1010 * eliminating what is likely to be the largest value, blue-scale:
1012 * blue-scale = white-scale - red-scale - green-scale
1016 * (red-x - blue-x)*red-scale + (green-x - blue-x)*green-scale =
1017 * (white-x - blue-x)*white-scale
1019 * (red-y - blue-y)*red-scale + (green-y - blue-y)*green-scale =
1020 * 1 - blue-y*white-scale
1022 * And now we can trivially solve for (red-scale,green-scale):
1025 * (white-x - blue-x)*white-scale - (red-x - blue-x)*red-scale
1026 * -----------------------------------------------------------
1030 * 1 - blue-y*white-scale - (green-y - blue-y) * green-scale
1031 * ---------------------------------------------------------
1037 * ( (green-x - blue-x) * (white-y - blue-y) -
1038 * (green-y - blue-y) * (white-x - blue-x) ) / white-y
1039 * -------------------------------------------------------------------------
1040 * (green-x - blue-x)*(red-y - blue-y)-(green-y - blue-y)*(red-x - blue-x)
1043 * ( (red-y - blue-y) * (white-x - blue-x) -
1044 * (red-x - blue-x) * (white-y - blue-y) ) / white-y
1045 * -------------------------------------------------------------------------
1046 * (green-x - blue-x)*(red-y - blue-y)-(green-y - blue-y)*(red-x - blue-x)
1049 * The input values have 5 decimal digits of accuracy. The values are all in
1050 * the range 0 < value < 1, so simple products are in the same range but may
1051 * need up to 10 decimal digits to preserve the original precision and avoid
1052 * underflow. Because we are using a 32-bit signed representation we cannot
1053 * match this; the best is a little over 9 decimal digits, less than 10.
1055 * The approach used here is to preserve the maximum precision within the
1056 * signed representation. Because the red-scale calculation above uses the
1057 * difference between two products of values that must be in the range -1..+1
1058 * it is sufficient to divide the product by 7; ceil(100,000/32767*2). The
1059 * factor is irrelevant in the calculation because it is applied to both
1060 * numerator and denominator.
1062 * Note that the values of the differences of the products of the
1063 * chromaticities in the above equations tend to be small, for example for
1064 * the sRGB chromaticities they are:
1066 * red numerator: -0.04751
1067 * green numerator: -0.08788
1068 * denominator: -0.2241 (without white-y multiplication)
1070 * The resultant Y coefficients from the chromaticities of some widely used
1071 * color space definitions are (to 15 decimal places):
1074 * 0.212639005871510 0.715168678767756 0.072192315360734
1076 * 0.288071128229293 0.711843217810102 0.000085653960605
1078 * 0.297344975250536 0.627363566255466 0.075291458493998
1079 * Adobe Wide Gamut RGB
1080 * 0.258728243040113 0.724682314948566 0.016589442011321
1082 /* By the argument, above overflow should be impossible here. The return
1083 * value of 2 indicates an internal error to the caller.
1085 if (!png_muldiv(&left, xy.greenx-xy.bluex, xy.redy - xy.bluey, 7)) return 2;
1086 if (!png_muldiv(&right, xy.greeny-xy.bluey, xy.redx - xy.bluex, 7)) return 2;
1087 denominator = left - right;
1089 /* Now find the red numerator. */
1090 if (!png_muldiv(&left, xy.greenx-xy.bluex, xy.whitey-xy.bluey, 7)) return 2;
1091 if (!png_muldiv(&right, xy.greeny-xy.bluey, xy.whitex-xy.bluex, 7)) return 2;
1093 /* Overflow is possible here and it indicates an extreme set of PNG cHRM
1094 * chunk values. This calculation actually returns the reciprocal of the
1095 * scale value because this allows us to delay the multiplication of white-y
1096 * into the denominator, which tends to produce a small number.
1098 if (!png_muldiv(&red_inverse, xy.whitey, denominator, left-right) ||
1099 red_inverse <= xy.whitey /* r+g+b scales = white scale */)
1102 /* Similarly for green_inverse: */
1103 if (!png_muldiv(&left, xy.redy-xy.bluey, xy.whitex-xy.bluex, 7)) return 2;
1104 if (!png_muldiv(&right, xy.redx-xy.bluex, xy.whitey-xy.bluey, 7)) return 2;
1105 if (!png_muldiv(&green_inverse, xy.whitey, denominator, left-right) ||
1106 green_inverse <= xy.whitey)
1109 /* And the blue scale, the checks above guarantee this can't overflow but it
1110 * can still produce 0 for extreme cHRM values.
1112 blue_scale = png_reciprocal(xy.whitey) - png_reciprocal(red_inverse) -
1113 png_reciprocal(green_inverse);
1114 if (blue_scale <= 0) return 1;
1117 /* And fill in the png_XYZ: */
1118 if (!png_muldiv(&XYZ->redX, xy.redx, PNG_FP_1, red_inverse)) return 1;
1119 if (!png_muldiv(&XYZ->redY, xy.redy, PNG_FP_1, red_inverse)) return 1;
1120 if (!png_muldiv(&XYZ->redZ, PNG_FP_1 - xy.redx - xy.redy, PNG_FP_1,
1124 if (!png_muldiv(&XYZ->greenX, xy.greenx, PNG_FP_1, green_inverse)) return 1;
1125 if (!png_muldiv(&XYZ->greenY, xy.greeny, PNG_FP_1, green_inverse)) return 1;
1126 if (!png_muldiv(&XYZ->greenZ, PNG_FP_1 - xy.greenx - xy.greeny, PNG_FP_1,
1130 if (!png_muldiv(&XYZ->blueX, xy.bluex, blue_scale, PNG_FP_1)) return 1;
1131 if (!png_muldiv(&XYZ->blueY, xy.bluey, blue_scale, PNG_FP_1)) return 1;
1132 if (!png_muldiv(&XYZ->blueZ, PNG_FP_1 - xy.bluex - xy.bluey, blue_scale,
1136 return 0; /*success*/
1139 int png_XYZ_from_xy_checked(png_structp png_ptr, png_XYZ *XYZ, png_xy xy)
1141 switch (png_XYZ_from_xy(XYZ, xy))
1143 case 0: /* success */
1147 /* The chunk may be technically valid, but we got png_fixed_point
1148 * overflow while trying to get XYZ values out of it. This is
1149 * entirely benign - the cHRM chunk is pretty extreme.
1151 png_warning(png_ptr,
1152 "extreme cHRM chunk cannot be converted to tristimulus values");
1156 /* libpng is broken; this should be a warning but if it happens we
1157 * want error reports so for the moment it is an error.
1159 png_error(png_ptr, "internal error in png_XYZ_from_xy");
1169 png_check_IHDR(png_structp png_ptr,
1170 png_uint_32 width, png_uint_32 height, int bit_depth,
1171 int color_type, int interlace_type, int compression_type,
1176 /* Check for width and height valid values */
1179 png_warning(png_ptr, "Image width is zero in IHDR");
1185 png_warning(png_ptr, "Image height is zero in IHDR");
1189 # ifdef PNG_SET_USER_LIMITS_SUPPORTED
1190 if (width > png_ptr->user_width_max)
1193 if (width > PNG_USER_WIDTH_MAX)
1196 png_warning(png_ptr, "Image width exceeds user limit in IHDR");
1200 # ifdef PNG_SET_USER_LIMITS_SUPPORTED
1201 if (height > png_ptr->user_height_max)
1203 if (height > PNG_USER_HEIGHT_MAX)
1206 png_warning(png_ptr, "Image height exceeds user limit in IHDR");
1210 if (width > PNG_UINT_31_MAX)
1212 png_warning(png_ptr, "Invalid image width in IHDR");
1216 if (height > PNG_UINT_31_MAX)
1218 png_warning(png_ptr, "Invalid image height in IHDR");
1222 if (width > (PNG_UINT_32_MAX
1223 >> 3) /* 8-byte RGBA pixels */
1224 - 48 /* bigrowbuf hack */
1225 - 1 /* filter byte */
1226 - 7*8 /* rounding of width to multiple of 8 pixels */
1227 - 8) /* extra max_pixel_depth pad */
1228 png_warning(png_ptr, "Width is too large for libpng to process pixels");
1230 /* Check other values */
1231 if (bit_depth != 1 && bit_depth != 2 && bit_depth != 4 &&
1232 bit_depth != 8 && bit_depth != 16)
1234 png_warning(png_ptr, "Invalid bit depth in IHDR");
1238 if (color_type < 0 || color_type == 1 ||
1239 color_type == 5 || color_type > 6)
1241 png_warning(png_ptr, "Invalid color type in IHDR");
1245 if (((color_type == PNG_COLOR_TYPE_PALETTE) && bit_depth > 8) ||
1246 ((color_type == PNG_COLOR_TYPE_RGB ||
1247 color_type == PNG_COLOR_TYPE_GRAY_ALPHA ||
1248 color_type == PNG_COLOR_TYPE_RGB_ALPHA) && bit_depth < 8))
1250 png_warning(png_ptr, "Invalid color type/bit depth combination in IHDR");
1254 if (interlace_type >= PNG_INTERLACE_LAST)
1256 png_warning(png_ptr, "Unknown interlace method in IHDR");
1260 if (compression_type != PNG_COMPRESSION_TYPE_BASE)
1262 png_warning(png_ptr, "Unknown compression method in IHDR");
1266 # ifdef PNG_MNG_FEATURES_SUPPORTED
1267 /* Accept filter_method 64 (intrapixel differencing) only if
1268 * 1. Libpng was compiled with PNG_MNG_FEATURES_SUPPORTED and
1269 * 2. Libpng did not read a PNG signature (this filter_method is only
1270 * used in PNG datastreams that are embedded in MNG datastreams) and
1271 * 3. The application called png_permit_mng_features with a mask that
1272 * included PNG_FLAG_MNG_FILTER_64 and
1273 * 4. The filter_method is 64 and
1274 * 5. The color_type is RGB or RGBA
1276 if ((png_ptr->mode & PNG_HAVE_PNG_SIGNATURE) &&
1277 png_ptr->mng_features_permitted)
1278 png_warning(png_ptr, "MNG features are not allowed in a PNG datastream");
1280 if (filter_type != PNG_FILTER_TYPE_BASE)
1282 if (!((png_ptr->mng_features_permitted & PNG_FLAG_MNG_FILTER_64) &&
1283 (filter_type == PNG_INTRAPIXEL_DIFFERENCING) &&
1284 ((png_ptr->mode & PNG_HAVE_PNG_SIGNATURE) == 0) &&
1285 (color_type == PNG_COLOR_TYPE_RGB ||
1286 color_type == PNG_COLOR_TYPE_RGB_ALPHA)))
1288 png_warning(png_ptr, "Unknown filter method in IHDR");
1292 if (png_ptr->mode & PNG_HAVE_PNG_SIGNATURE)
1294 png_warning(png_ptr, "Invalid filter method in IHDR");
1300 if (filter_type != PNG_FILTER_TYPE_BASE)
1302 png_warning(png_ptr, "Unknown filter method in IHDR");
1308 png_error(png_ptr, "Invalid IHDR data");
1311 #if defined(PNG_sCAL_SUPPORTED) || defined(PNG_pCAL_SUPPORTED)
1312 /* ASCII to fp functions */
1313 /* Check an ASCII formated floating point value, see the more detailed
1314 * comments in pngpriv.h
1316 /* The following is used internally to preserve the sticky flags */
1317 #define png_fp_add(state, flags) ((state) |= (flags))
1318 #define png_fp_set(state, value) ((state) = (value) | ((state) & PNG_FP_STICKY))
1321 png_check_fp_number(png_const_charp string, png_size_t size, int *statep,
1322 png_size_tp whereami)
1324 int state = *statep;
1325 png_size_t i = *whereami;
1330 /* First find the type of the next character */
1333 case 43: type = PNG_FP_SAW_SIGN; break;
1334 case 45: type = PNG_FP_SAW_SIGN + PNG_FP_NEGATIVE; break;
1335 case 46: type = PNG_FP_SAW_DOT; break;
1336 case 48: type = PNG_FP_SAW_DIGIT; break;
1337 case 49: case 50: case 51: case 52:
1338 case 53: case 54: case 55: case 56:
1339 case 57: type = PNG_FP_SAW_DIGIT + PNG_FP_NONZERO; break;
1341 case 101: type = PNG_FP_SAW_E; break;
1342 default: goto PNG_FP_End;
1345 /* Now deal with this type according to the current
1346 * state, the type is arranged to not overlap the
1347 * bits of the PNG_FP_STATE.
1349 switch ((state & PNG_FP_STATE) + (type & PNG_FP_SAW_ANY))
1351 case PNG_FP_INTEGER + PNG_FP_SAW_SIGN:
1352 if (state & PNG_FP_SAW_ANY)
1353 goto PNG_FP_End; /* not a part of the number */
1355 png_fp_add(state, type);
1358 case PNG_FP_INTEGER + PNG_FP_SAW_DOT:
1359 /* Ok as trailer, ok as lead of fraction. */
1360 if (state & PNG_FP_SAW_DOT) /* two dots */
1363 else if (state & PNG_FP_SAW_DIGIT) /* trailing dot? */
1364 png_fp_add(state, type);
1367 png_fp_set(state, PNG_FP_FRACTION | type);
1371 case PNG_FP_INTEGER + PNG_FP_SAW_DIGIT:
1372 if (state & PNG_FP_SAW_DOT) /* delayed fraction */
1373 png_fp_set(state, PNG_FP_FRACTION | PNG_FP_SAW_DOT);
1375 png_fp_add(state, type | PNG_FP_WAS_VALID);
1379 case PNG_FP_INTEGER + PNG_FP_SAW_E:
1380 if ((state & PNG_FP_SAW_DIGIT) == 0)
1383 png_fp_set(state, PNG_FP_EXPONENT);
1387 /* case PNG_FP_FRACTION + PNG_FP_SAW_SIGN:
1388 goto PNG_FP_End; ** no sign in fraction */
1390 /* case PNG_FP_FRACTION + PNG_FP_SAW_DOT:
1391 goto PNG_FP_End; ** Because SAW_DOT is always set */
1393 case PNG_FP_FRACTION + PNG_FP_SAW_DIGIT:
1394 png_fp_add(state, type | PNG_FP_WAS_VALID);
1397 case PNG_FP_FRACTION + PNG_FP_SAW_E:
1398 /* This is correct because the trailing '.' on an
1399 * integer is handled above - so we can only get here
1400 * with the sequence ".E" (with no preceding digits).
1402 if ((state & PNG_FP_SAW_DIGIT) == 0)
1405 png_fp_set(state, PNG_FP_EXPONENT);
1409 case PNG_FP_EXPONENT + PNG_FP_SAW_SIGN:
1410 if (state & PNG_FP_SAW_ANY)
1411 goto PNG_FP_End; /* not a part of the number */
1413 png_fp_add(state, PNG_FP_SAW_SIGN);
1417 /* case PNG_FP_EXPONENT + PNG_FP_SAW_DOT:
1420 case PNG_FP_EXPONENT + PNG_FP_SAW_DIGIT:
1421 png_fp_add(state, PNG_FP_SAW_DIGIT | PNG_FP_WAS_VALID);
1425 /* case PNG_FP_EXPONEXT + PNG_FP_SAW_E:
1428 default: goto PNG_FP_End; /* I.e. break 2 */
1431 /* The character seems ok, continue. */
1436 /* Here at the end, update the state and return the correct
1442 return (state & PNG_FP_SAW_DIGIT) != 0;
1446 /* The same but for a complete string. */
1448 png_check_fp_string(png_const_charp string, png_size_t size)
1451 png_size_t char_index=0;
1453 if (png_check_fp_number(string, size, &state, &char_index) &&
1454 (char_index == size || string[char_index] == 0))
1455 return state /* must be non-zero - see above */;
1457 return 0; /* i.e. fail */
1459 #endif /* pCAL or sCAL */
1461 #ifdef PNG_READ_sCAL_SUPPORTED
1462 # ifdef PNG_FLOATING_POINT_SUPPORTED
1463 /* Utility used below - a simple accurate power of ten from an integral
1467 png_pow10(int power)
1472 /* Handle negative exponent with a reciprocal at the end because
1473 * 10 is exact whereas .1 is inexact in base 2
1477 if (power < DBL_MIN_10_EXP) return 0;
1478 recip = 1, power = -power;
1483 /* Decompose power bitwise. */
1487 if (power & 1) d *= mult;
1495 /* else power is 0 and d is 1 */
1500 /* Function to format a floating point value in ASCII with a given
1504 png_ascii_from_fp(png_structp png_ptr, png_charp ascii, png_size_t size,
1505 double fp, unsigned int precision)
1507 /* We use standard functions from math.h, but not printf because
1508 * that would require stdio. The caller must supply a buffer of
1509 * sufficient size or we will png_error. The tests on size and
1510 * the space in ascii[] consumed are indicated below.
1513 precision = DBL_DIG;
1515 /* Enforce the limit of the implementation precision too. */
1516 if (precision > DBL_DIG+1)
1517 precision = DBL_DIG+1;
1519 /* Basic sanity checks */
1520 if (size >= precision+5) /* See the requirements below. */
1525 *ascii++ = 45; /* '-' PLUS 1 TOTAL 1 */
1529 if (fp >= DBL_MIN && fp <= DBL_MAX)
1531 int exp_b10; /* A base 10 exponent */
1532 double base; /* 10^exp_b10 */
1534 /* First extract a base 10 exponent of the number,
1535 * the calculation below rounds down when converting
1536 * from base 2 to base 10 (multiply by log10(2) -
1537 * 0.3010, but 77/256 is 0.3008, so exp_b10 needs to
1538 * be increased. Note that the arithmetic shift
1539 * performs a floor() unlike C arithmetic - using a
1540 * C multiply would break the following for negative
1543 (void)frexp(fp, &exp_b10); /* exponent to base 2 */
1545 exp_b10 = (exp_b10 * 77) >> 8; /* <= exponent to base 10 */
1547 /* Avoid underflow here. */
1548 base = png_pow10(exp_b10); /* May underflow */
1550 while (base < DBL_MIN || base < fp)
1552 /* And this may overflow. */
1553 double test = png_pow10(exp_b10+1);
1555 if (test <= DBL_MAX)
1556 ++exp_b10, base = test;
1562 /* Normalize fp and correct exp_b10, after this fp is in the
1563 * range [.1,1) and exp_b10 is both the exponent and the digit
1564 * *before* which the decimal point should be inserted
1565 * (starting with 0 for the first digit). Note that this
1566 * works even if 10^exp_b10 is out of range because of the
1567 * test on DBL_MAX above.
1570 while (fp >= 1) fp /= 10, ++exp_b10;
1572 /* Because of the code above fp may, at this point, be
1573 * less than .1, this is ok because the code below can
1574 * handle the leading zeros this generates, so no attempt
1575 * is made to correct that here.
1579 int czero, clead, cdigits;
1582 /* Allow up to two leading zeros - this will not lengthen
1583 * the number compared to using E-n.
1585 if (exp_b10 < 0 && exp_b10 > -3) /* PLUS 3 TOTAL 4 */
1587 czero = -exp_b10; /* PLUS 2 digits: TOTAL 3 */
1588 exp_b10 = 0; /* Dot added below before first output. */
1591 czero = 0; /* No zeros to add */
1593 /* Generate the digit list, stripping trailing zeros and
1594 * inserting a '.' before a digit if the exponent is 0.
1596 clead = czero; /* Count of leading zeros */
1597 cdigits = 0; /* Count of digits in list. */
1605 /* Use modf here, not floor and subtract, so that
1606 * the separation is done in one step. At the end
1607 * of the loop don't break the number into parts so
1608 * that the final digit is rounded.
1610 if (cdigits+czero-clead+1 < (int)precision)
1619 /* Rounding up to 10, handle that here. */
1623 if (cdigits == 0) --clead;
1628 while (cdigits > 0 && d > 9.0)
1632 if (exp_b10 != (-1))
1637 ch = *--ascii, ++size;
1638 /* Advance exp_b10 to '1', so that the
1639 * decimal point happens after the
1646 d = ch - 47; /* I.e. 1+(ch-48) */
1649 /* Did we reach the beginning? If so adjust the
1650 * exponent but take into account the leading
1653 if (d > 9.0) /* cdigits == 0 */
1655 if (exp_b10 == (-1))
1657 /* Leading decimal point (plus zeros?), if
1658 * we lose the decimal point here it must
1659 * be reentered below.
1664 ++size, exp_b10 = 1;
1666 /* Else lost a leading zero, so 'exp_b10' is
1673 /* In all cases we output a '1' */
1678 fp = 0; /* Guarantees termination below. */
1684 if (cdigits == 0) ++clead;
1689 /* Included embedded zeros in the digit count. */
1690 cdigits += czero - clead;
1695 /* exp_b10 == (-1) means we just output the decimal
1696 * place - after the DP don't adjust 'exp_b10' any
1699 if (exp_b10 != (-1))
1701 if (exp_b10 == 0) *ascii++ = 46, --size;
1702 /* PLUS 1: TOTAL 4 */
1705 *ascii++ = 48, --czero;
1708 if (exp_b10 != (-1))
1710 if (exp_b10 == 0) *ascii++ = 46, --size; /* counted
1715 *ascii++ = (char)(48 + (int)d), ++cdigits;
1718 while (cdigits+czero-clead < (int)precision && fp > DBL_MIN);
1720 /* The total output count (max) is now 4+precision */
1722 /* Check for an exponent, if we don't need one we are
1723 * done and just need to terminate the string. At
1724 * this point exp_b10==(-1) is effectively if flag - it got
1725 * to '-1' because of the decrement after outputing
1726 * the decimal point above (the exponent required is
1729 if (exp_b10 >= (-1) && exp_b10 <= 2)
1731 /* The following only happens if we didn't output the
1732 * leading zeros above for negative exponent, so this
1733 * doest add to the digit requirement. Note that the
1734 * two zeros here can only be output if the two leading
1735 * zeros were *not* output, so this doesn't increase
1738 while (--exp_b10 >= 0) *ascii++ = 48;
1742 /* Total buffer requirement (including the '\0') is
1743 * 5+precision - see check at the start.
1748 /* Here if an exponent is required, adjust size for
1749 * the digits we output but did not count. The total
1750 * digit output here so far is at most 1+precision - no
1751 * decimal point and no leading or trailing zeros have
1756 *ascii++ = 69, --size; /* 'E': PLUS 1 TOTAL 2+precision */
1758 /* The following use of an unsigned temporary avoids ambiguities in
1759 * the signed arithmetic on exp_b10 and permits GCC at least to do
1760 * better optimization.
1763 unsigned int uexp_b10;
1767 *ascii++ = 45, --size; /* '-': PLUS 1 TOTAL 3+precision */
1768 uexp_b10 = -exp_b10;
1776 while (uexp_b10 > 0)
1778 exponent[cdigits++] = (char)(48 + uexp_b10 % 10);
1783 /* Need another size check here for the exponent digits, so
1784 * this need not be considered above.
1786 if ((int)size > cdigits)
1788 while (cdigits > 0) *ascii++ = exponent[--cdigits];
1796 else if (!(fp >= DBL_MIN))
1798 *ascii++ = 48; /* '0' */
1804 *ascii++ = 105; /* 'i' */
1805 *ascii++ = 110; /* 'n' */
1806 *ascii++ = 102; /* 'f' */
1812 /* Here on buffer too small. */
1813 png_error(png_ptr, "ASCII conversion buffer too small");
1816 # endif /* FLOATING_POINT */
1818 # ifdef PNG_FIXED_POINT_SUPPORTED
1819 /* Function to format a fixed point value in ASCII.
1822 png_ascii_from_fixed(png_structp png_ptr, png_charp ascii, png_size_t size,
1825 /* Require space for 10 decimal digits, a decimal point, a minus sign and a
1826 * trailing \0, 13 characters:
1832 /* Avoid overflow here on the minimum integer. */
1834 *ascii++ = 45, --size, num = -fp;
1838 if (num <= 0x80000000) /* else overflowed */
1840 unsigned int ndigits = 0, first = 16 /* flag value */;
1845 /* Split the low digit off num: */
1846 unsigned int tmp = num/10;
1848 digits[ndigits++] = (char)(48 + num);
1849 /* Record the first non-zero digit, note that this is a number
1850 * starting at 1, it's not actually the array index.
1852 if (first == 16 && num > 0)
1859 while (ndigits > 5) *ascii++ = digits[--ndigits];
1860 /* The remaining digits are fractional digits, ndigits is '5' or
1861 * smaller at this point. It is certainly not zero. Check for a
1862 * non-zero fractional digit:
1867 *ascii++ = 46; /* decimal point */
1868 /* ndigits may be <5 for small numbers, output leading zeros
1869 * then ndigits digits to first:
1872 while (ndigits < i) *ascii++ = 48, --i;
1873 while (ndigits >= first) *ascii++ = digits[--ndigits];
1874 /* Don't output the trailing zeros! */
1880 /* And null terminate the string: */
1886 /* Here on buffer too small. */
1887 png_error(png_ptr, "ASCII conversion buffer too small");
1889 # endif /* FIXED_POINT */
1890 #endif /* READ_SCAL */
1892 #if defined(PNG_FLOATING_POINT_SUPPORTED) && \
1893 !defined(PNG_FIXED_POINT_MACRO_SUPPORTED)
1895 png_fixed(png_structp png_ptr, double fp, png_const_charp text)
1897 double r = floor(100000 * fp + .5);
1899 if (r > 2147483647. || r < -2147483648.)
1900 png_fixed_error(png_ptr, text);
1902 return (png_fixed_point)r;
1906 #if defined(PNG_READ_GAMMA_SUPPORTED) || \
1907 defined(PNG_INCH_CONVERSIONS_SUPPORTED) || defined(PNG__READ_pHYs_SUPPORTED)
1908 /* muldiv functions */
1909 /* This API takes signed arguments and rounds the result to the nearest
1910 * integer (or, for a fixed point number - the standard argument - to
1911 * the nearest .00001). Overflow and divide by zero are signalled in
1912 * the result, a boolean - true on success, false on overflow.
1915 png_muldiv(png_fixed_point_p res, png_fixed_point a, png_int_32 times,
1918 /* Return a * times / divisor, rounded. */
1921 if (a == 0 || times == 0)
1928 #ifdef PNG_FLOATING_ARITHMETIC_SUPPORTED
1934 /* A png_fixed_point is a 32-bit integer. */
1935 if (r <= 2147483647. && r >= -2147483648.)
1937 *res = (png_fixed_point)r;
1942 png_uint_32 A, T, D;
1943 png_uint_32 s16, s32, s00;
1946 negative = 1, A = -a;
1951 negative = !negative, T = -times;
1956 negative = !negative, D = -divisor;
1960 /* Following can't overflow because the arguments only
1961 * have 31 bits each, however the result may be 32 bits.
1963 s16 = (A >> 16) * (T & 0xffff) +
1964 (A & 0xffff) * (T >> 16);
1965 /* Can't overflow because the a*times bit is only 30
1968 s32 = (A >> 16) * (T >> 16) + (s16 >> 16);
1969 s00 = (A & 0xffff) * (T & 0xffff);
1971 s16 = (s16 & 0xffff) << 16;
1977 if (s32 < D) /* else overflow */
1979 /* s32.s00 is now the 64-bit product, do a standard
1980 * division, we know that s32 < D, so the maximum
1981 * required shift is 31.
1984 png_fixed_point result = 0; /* NOTE: signed */
1986 while (--bitshift >= 0)
1988 png_uint_32 d32, d00;
1991 d32 = D >> (32-bitshift), d00 = D << bitshift;
1998 if (s00 < d00) --s32; /* carry */
1999 s32 -= d32, s00 -= d00, result += 1<<bitshift;
2003 if (s32 == d32 && s00 >= d00)
2004 s32 = 0, s00 -= d00, result += 1<<bitshift;
2007 /* Handle the rounding. */
2008 if (s00 >= (D >> 1))
2014 /* Check for overflow. */
2015 if ((negative && result <= 0) || (!negative && result >= 0))
2027 #endif /* READ_GAMMA || INCH_CONVERSIONS */
2029 #if defined(PNG_READ_GAMMA_SUPPORTED) || defined(PNG_INCH_CONVERSIONS_SUPPORTED)
2030 /* The following is for when the caller doesn't much care about the
2034 png_muldiv_warn(png_structp png_ptr, png_fixed_point a, png_int_32 times,
2037 png_fixed_point result;
2039 if (png_muldiv(&result, a, times, divisor))
2042 png_warning(png_ptr, "fixed point overflow ignored");
2047 #ifdef PNG_READ_GAMMA_SUPPORTED /* more fixed point functions for gamma */
2048 /* Calculate a reciprocal, return 0 on div-by-zero or overflow. */
2050 png_reciprocal(png_fixed_point a)
2052 #ifdef PNG_FLOATING_ARITHMETIC_SUPPORTED
2053 double r = floor(1E10/a+.5);
2055 if (r <= 2147483647. && r >= -2147483648.)
2056 return (png_fixed_point)r;
2058 png_fixed_point res;
2060 if (png_muldiv(&res, 100000, 100000, a))
2064 return 0; /* error/overflow */
2067 /* A local convenience routine. */
2068 static png_fixed_point
2069 png_product2(png_fixed_point a, png_fixed_point b)
2071 /* The required result is 1/a * 1/b; the following preserves accuracy. */
2072 #ifdef PNG_FLOATING_ARITHMETIC_SUPPORTED
2073 double r = a * 1E-5;
2077 if (r <= 2147483647. && r >= -2147483648.)
2078 return (png_fixed_point)r;
2080 png_fixed_point res;
2082 if (png_muldiv(&res, a, b, 100000))
2086 return 0; /* overflow */
2089 /* The inverse of the above. */
2091 png_reciprocal2(png_fixed_point a, png_fixed_point b)
2093 /* The required result is 1/a * 1/b; the following preserves accuracy. */
2094 #ifdef PNG_FLOATING_ARITHMETIC_SUPPORTED
2099 if (r <= 2147483647. && r >= -2147483648.)
2100 return (png_fixed_point)r;
2102 /* This may overflow because the range of png_fixed_point isn't symmetric,
2103 * but this API is only used for the product of file and screen gamma so it
2104 * doesn't matter that the smallest number it can produce is 1/21474, not
2107 png_fixed_point res = png_product2(a, b);
2110 return png_reciprocal(res);
2113 return 0; /* overflow */
2115 #endif /* READ_GAMMA */
2117 #ifdef PNG_CHECK_cHRM_SUPPORTED
2118 /* Added at libpng version 1.2.34 (Dec 8, 2008) and 1.4.0 (Jan 2,
2119 * 2010: moved from pngset.c) */
2121 * Multiply two 32-bit numbers, V1 and V2, using 32-bit
2122 * arithmetic, to produce a 64-bit result in the HI/LO words.
2130 * where A and B are the high and low 16-bit words of V1,
2131 * C and D are the 16-bit words of V2, AD is the product of
2132 * A and D, and X || Y is (X << 16) + Y.
2136 png_64bit_product (long v1, long v2, unsigned long *hi_product,
2137 unsigned long *lo_product)
2142 a = (v1 >> 16) & 0xffff;
2144 c = (v2 >> 16) & 0xffff;
2147 lo = b * d; /* BD */
2148 x = a * d + c * b; /* AD + CB */
2149 y = ((lo >> 16) & 0xffff) + x;
2151 lo = (lo & 0xffff) | ((y & 0xffff) << 16);
2152 hi = (y >> 16) & 0xffff;
2154 hi += a * c; /* AC */
2156 *hi_product = (unsigned long)hi;
2157 *lo_product = (unsigned long)lo;
2159 #endif /* CHECK_cHRM */
2161 #ifdef PNG_READ_GAMMA_SUPPORTED /* gamma table code */
2162 #ifndef PNG_FLOATING_ARITHMETIC_SUPPORTED
2163 /* Fixed point gamma.
2165 * To calculate gamma this code implements fast log() and exp() calls using only
2166 * fixed point arithmetic. This code has sufficient precision for either 8-bit
2167 * or 16-bit sample values.
2169 * The tables used here were calculated using simple 'bc' programs, but C double
2170 * precision floating point arithmetic would work fine. The programs are given
2171 * at the head of each table.
2174 * This is a table of -log(value/255)/log(2) for 'value' in the range 128 to
2175 * 255, so it's the base 2 logarithm of a normalized 8-bit floating point
2176 * mantissa. The numbers are 32-bit fractions.
2182 for (i=128;i<256;++i) { .5 - l(i/255)/l(2)*65536*65536; }
2184 4270715492U, 4222494797U, 4174646467U, 4127164793U, 4080044201U, 4033279239U,
2185 3986864580U, 3940795015U, 3895065449U, 3849670902U, 3804606499U, 3759867474U,
2186 3715449162U, 3671346997U, 3627556511U, 3584073329U, 3540893168U, 3498011834U,
2187 3455425220U, 3413129301U, 3371120137U, 3329393864U, 3287946700U, 3246774933U,
2188 3205874930U, 3165243125U, 3124876025U, 3084770202U, 3044922296U, 3005329011U,
2189 2965987113U, 2926893432U, 2888044853U, 2849438323U, 2811070844U, 2772939474U,
2190 2735041326U, 2697373562U, 2659933400U, 2622718104U, 2585724991U, 2548951424U,
2191 2512394810U, 2476052606U, 2439922311U, 2404001468U, 2368287663U, 2332778523U,
2192 2297471715U, 2262364947U, 2227455964U, 2192742551U, 2158222529U, 2123893754U,
2193 2089754119U, 2055801552U, 2022034013U, 1988449497U, 1955046031U, 1921821672U,
2194 1888774511U, 1855902668U, 1823204291U, 1790677560U, 1758320682U, 1726131893U,
2195 1694109454U, 1662251657U, 1630556815U, 1599023271U, 1567649391U, 1536433567U,
2196 1505374214U, 1474469770U, 1443718700U, 1413119487U, 1382670639U, 1352370686U,
2197 1322218179U, 1292211689U, 1262349810U, 1232631153U, 1203054352U, 1173618059U,
2198 1144320946U, 1115161701U, 1086139034U, 1057251672U, 1028498358U, 999877854U,
2199 971388940U, 943030410U, 914801076U, 886699767U, 858725327U, 830876614U,
2200 803152505U, 775551890U, 748073672U, 720716771U, 693480120U, 666362667U,
2201 639363374U, 612481215U, 585715177U, 559064263U, 532527486U, 506103872U,
2202 479792461U, 453592303U, 427502463U, 401522014U, 375650043U, 349885648U,
2203 324227938U, 298676034U, 273229066U, 247886176U, 222646516U, 197509248U,
2204 172473545U, 147538590U, 122703574U, 97967701U, 73330182U, 48790236U,
2209 /* The following are the values for 16-bit tables - these work fine for the
2210 * 8-bit conversions but produce very slightly larger errors in the 16-bit
2211 * log (about 1.2 as opposed to 0.7 absolute error in the final value). To
2212 * use these all the shifts below must be adjusted appropriately.
2214 65166, 64430, 63700, 62976, 62257, 61543, 60835, 60132, 59434, 58741, 58054,
2215 57371, 56693, 56020, 55352, 54689, 54030, 53375, 52726, 52080, 51439, 50803,
2216 50170, 49542, 48918, 48298, 47682, 47070, 46462, 45858, 45257, 44661, 44068,
2217 43479, 42894, 42312, 41733, 41159, 40587, 40020, 39455, 38894, 38336, 37782,
2218 37230, 36682, 36137, 35595, 35057, 34521, 33988, 33459, 32932, 32408, 31887,
2219 31369, 30854, 30341, 29832, 29325, 28820, 28319, 27820, 27324, 26830, 26339,
2220 25850, 25364, 24880, 24399, 23920, 23444, 22970, 22499, 22029, 21562, 21098,
2221 20636, 20175, 19718, 19262, 18808, 18357, 17908, 17461, 17016, 16573, 16132,
2222 15694, 15257, 14822, 14390, 13959, 13530, 13103, 12678, 12255, 11834, 11415,
2223 10997, 10582, 10168, 9756, 9346, 8937, 8531, 8126, 7723, 7321, 6921, 6523,
2224 6127, 5732, 5339, 4947, 4557, 4169, 3782, 3397, 3014, 2632, 2251, 1872, 1495,
2229 PNG_STATIC png_int_32
2230 png_log8bit(unsigned int x)
2232 unsigned int lg2 = 0;
2233 /* Each time 'x' is multiplied by 2, 1 must be subtracted off the final log,
2234 * because the log is actually negate that means adding 1. The final
2235 * returned value thus has the range 0 (for 255 input) to 7.994 (for 1
2236 * input), return 7.99998 for the overflow (log 0) case - so the result is
2237 * always at most 19 bits.
2239 if ((x &= 0xff) == 0)
2242 if ((x & 0xf0) == 0)
2245 if ((x & 0xc0) == 0)
2248 if ((x & 0x80) == 0)
2251 /* result is at most 19 bits, so this cast is safe: */
2252 return (png_int_32)((lg2 << 16) + ((png_8bit_l2[x-128]+32768)>>16));
2255 /* The above gives exact (to 16 binary places) log2 values for 8-bit images,
2256 * for 16-bit images we use the most significant 8 bits of the 16-bit value to
2257 * get an approximation then multiply the approximation by a correction factor
2258 * determined by the remaining up to 8 bits. This requires an additional step
2259 * in the 16-bit case.
2261 * We want log2(value/65535), we have log2(v'/255), where:
2263 * value = v' * 256 + v''
2266 * So f is value/v', which is equal to (256+v''/v') since v' is in the range 128
2267 * to 255 and v'' is in the range 0 to 255 f will be in the range 256 to less
2268 * than 258. The final factor also needs to correct for the fact that our 8-bit
2269 * value is scaled by 255, whereas the 16-bit values must be scaled by 65535.
2271 * This gives a final formula using a calculated value 'x' which is value/v' and
2272 * scaling by 65536 to match the above table:
2274 * log2(x/257) * 65536
2276 * Since these numbers are so close to '1' we can use simple linear
2277 * interpolation between the two end values 256/257 (result -368.61) and 258/257
2278 * (result 367.179). The values used below are scaled by a further 64 to give
2279 * 16-bit precision in the interpolation:
2281 * Start (256): -23591
2285 PNG_STATIC png_int_32
2286 png_log16bit(png_uint_32 x)
2288 unsigned int lg2 = 0;
2290 /* As above, but now the input has 16 bits. */
2291 if ((x &= 0xffff) == 0)
2294 if ((x & 0xff00) == 0)
2297 if ((x & 0xf000) == 0)
2300 if ((x & 0xc000) == 0)
2303 if ((x & 0x8000) == 0)
2306 /* Calculate the base logarithm from the top 8 bits as a 28-bit fractional
2310 lg2 += (png_8bit_l2[(x>>8)-128]+8) >> 4;
2312 /* Now we need to interpolate the factor, this requires a division by the top
2313 * 8 bits. Do this with maximum precision.
2315 x = ((x << 16) + (x >> 9)) / (x >> 8);
2317 /* Since we divided by the top 8 bits of 'x' there will be a '1' at 1<<24,
2318 * the value at 1<<16 (ignoring this) will be 0 or 1; this gives us exactly
2319 * 16 bits to interpolate to get the low bits of the result. Round the
2320 * answer. Note that the end point values are scaled by 64 to retain overall
2321 * precision and that 'lg2' is current scaled by an extra 12 bits, so adjust
2322 * the overall scaling by 6-12. Round at every step.
2326 if (x <= 65536U) /* <= '257' */
2327 lg2 += ((23591U * (65536U-x)) + (1U << (16+6-12-1))) >> (16+6-12);
2330 lg2 -= ((23499U * (x-65536U)) + (1U << (16+6-12-1))) >> (16+6-12);
2332 /* Safe, because the result can't have more than 20 bits: */
2333 return (png_int_32)((lg2 + 2048) >> 12);
2336 /* The 'exp()' case must invert the above, taking a 20-bit fixed point
2337 * logarithmic value and returning a 16 or 8-bit number as appropriate. In
2338 * each case only the low 16 bits are relevant - the fraction - since the
2339 * integer bits (the top 4) simply determine a shift.
2341 * The worst case is the 16-bit distinction between 65535 and 65534, this
2342 * requires perhaps spurious accuracy in the decoding of the logarithm to
2343 * distinguish log2(65535/65534.5) - 10^-5 or 17 bits. There is little chance
2344 * of getting this accuracy in practice.
2346 * To deal with this the following exp() function works out the exponent of the
2347 * frational part of the logarithm by using an accurate 32-bit value from the
2348 * top four fractional bits then multiplying in the remaining bits.
2354 for (i=0;i<16;++i) { .5 + e(-i/16*l(2))*2^32; }
2356 /* NOTE: the first entry is deliberately set to the maximum 32-bit value. */
2357 4294967295U, 4112874773U, 3938502376U, 3771522796U, 3611622603U, 3458501653U,
2358 3311872529U, 3171459999U, 3037000500U, 2908241642U, 2784941738U, 2666869345U,
2359 2553802834U, 2445529972U, 2341847524U, 2242560872U
2363 /* Adjustment table; provided to explain the numbers in the code below. */
2365 for (i=11;i>=0;--i){ print i, " ", (1 - e(-(2^i)/65536*l(2))) * 2^(32-i), "\n"}
2366 11 44937.64284865548751208448
2367 10 45180.98734845585101160448
2368 9 45303.31936980687359311872
2369 8 45364.65110595323018870784
2370 7 45395.35850361789624614912
2371 6 45410.72259715102037508096
2372 5 45418.40724413220722311168
2373 4 45422.25021786898173001728
2374 3 45424.17186732298419044352
2375 2 45425.13273269940811464704
2376 1 45425.61317555035558641664
2377 0 45425.85339951654943850496
2380 PNG_STATIC png_uint_32
2381 png_exp(png_fixed_point x)
2383 if (x > 0 && x <= 0xfffff) /* Else overflow or zero (underflow) */
2385 /* Obtain a 4-bit approximation */
2386 png_uint_32 e = png_32bit_exp[(x >> 12) & 0xf];
2388 /* Incorporate the low 12 bits - these decrease the returned value by
2389 * multiplying by a number less than 1 if the bit is set. The multiplier
2390 * is determined by the above table and the shift. Notice that the values
2391 * converge on 45426 and this is used to allow linear interpolation of the
2395 e -= (((e >> 16) * 44938U) + 16U) >> 5;
2398 e -= (((e >> 16) * 45181U) + 32U) >> 6;
2401 e -= (((e >> 16) * 45303U) + 64U) >> 7;
2404 e -= (((e >> 16) * 45365U) + 128U) >> 8;
2407 e -= (((e >> 16) * 45395U) + 256U) >> 9;
2410 e -= (((e >> 16) * 45410U) + 512U) >> 10;
2412 /* And handle the low 6 bits in a single block. */
2413 e -= (((e >> 16) * 355U * (x & 0x3fU)) + 256U) >> 9;
2415 /* Handle the upper bits of x. */
2420 /* Check for overflow */
2422 return png_32bit_exp[0];
2424 /* Else underflow */
2429 png_exp8bit(png_fixed_point lg2)
2431 /* Get a 32-bit value: */
2432 png_uint_32 x = png_exp(lg2);
2434 /* Convert the 32-bit value to 0..255 by multiplying by 256-1, note that the
2435 * second, rounding, step can't overflow because of the first, subtraction,
2439 return (png_byte)((x + 0x7fffffU) >> 24);
2442 PNG_STATIC png_uint_16
2443 png_exp16bit(png_fixed_point lg2)
2445 /* Get a 32-bit value: */
2446 png_uint_32 x = png_exp(lg2);
2448 /* Convert the 32-bit value to 0..65535 by multiplying by 65536-1: */
2450 return (png_uint_16)((x + 32767U) >> 16);
2452 #endif /* FLOATING_ARITHMETIC */
2455 png_gamma_8bit_correct(unsigned int value, png_fixed_point gamma_val)
2457 if (value > 0 && value < 255)
2459 # ifdef PNG_FLOATING_ARITHMETIC_SUPPORTED
2460 double r = floor(255*pow(value/255.,gamma_val*.00001)+.5);
2463 png_int_32 lg2 = png_log8bit(value);
2464 png_fixed_point res;
2466 if (png_muldiv(&res, gamma_val, lg2, PNG_FP_1))
2467 return png_exp8bit(res);
2474 return (png_byte)value;
2478 png_gamma_16bit_correct(unsigned int value, png_fixed_point gamma_val)
2480 if (value > 0 && value < 65535)
2482 # ifdef PNG_FLOATING_ARITHMETIC_SUPPORTED
2483 double r = floor(65535*pow(value/65535.,gamma_val*.00001)+.5);
2484 return (png_uint_16)r;
2486 png_int_32 lg2 = png_log16bit(value);
2487 png_fixed_point res;
2489 if (png_muldiv(&res, gamma_val, lg2, PNG_FP_1))
2490 return png_exp16bit(res);
2497 return (png_uint_16)value;
2500 /* This does the right thing based on the bit_depth field of the
2501 * png_struct, interpreting values as 8-bit or 16-bit. While the result
2502 * is nominally a 16-bit value if bit depth is 8 then the result is
2503 * 8-bit (as are the arguments.)
2505 png_uint_16 /* PRIVATE */
2506 png_gamma_correct(png_structp png_ptr, unsigned int value,
2507 png_fixed_point gamma_val)
2509 if (png_ptr->bit_depth == 8)
2510 return png_gamma_8bit_correct(value, gamma_val);
2513 return png_gamma_16bit_correct(value, gamma_val);
2516 /* This is the shared test on whether a gamma value is 'significant' - whether
2517 * it is worth doing gamma correction.
2520 png_gamma_significant(png_fixed_point gamma_val)
2522 return gamma_val < PNG_FP_1 - PNG_GAMMA_THRESHOLD_FIXED ||
2523 gamma_val > PNG_FP_1 + PNG_GAMMA_THRESHOLD_FIXED;
2526 /* Internal function to build a single 16-bit table - the table consists of
2527 * 'num' 256-entry subtables, where 'num' is determined by 'shift' - the amount
2528 * to shift the input values right (or 16-number_of_signifiant_bits).
2530 * The caller is responsible for ensuring that the table gets cleaned up on
2531 * png_error (i.e. if one of the mallocs below fails) - i.e. the *table argument
2532 * should be somewhere that will be cleaned.
2535 png_build_16bit_table(png_structp png_ptr, png_uint_16pp *ptable,
2536 PNG_CONST unsigned int shift, PNG_CONST png_fixed_point gamma_val)
2538 /* Various values derived from 'shift': */
2539 PNG_CONST unsigned int num = 1U << (8U - shift);
2540 PNG_CONST unsigned int max = (1U << (16U - shift))-1U;
2541 PNG_CONST unsigned int max_by_2 = 1U << (15U-shift);
2544 png_uint_16pp table = *ptable =
2545 (png_uint_16pp)png_calloc(png_ptr, num * png_sizeof(png_uint_16p));
2547 for (i = 0; i < num; i++)
2549 png_uint_16p sub_table = table[i] =
2550 (png_uint_16p)png_malloc(png_ptr, 256 * png_sizeof(png_uint_16));
2552 /* The 'threshold' test is repeated here because it can arise for one of
2553 * the 16-bit tables even if the others don't hit it.
2555 if (png_gamma_significant(gamma_val))
2557 /* The old code would overflow at the end and this would cause the
2558 * 'pow' function to return a result >1, resulting in an
2559 * arithmetic error. This code follows the spec exactly; ig is
2560 * the recovered input sample, it always has 8-16 bits.
2562 * We want input * 65535/max, rounded, the arithmetic fits in 32
2563 * bits (unsigned) so long as max <= 32767.
2566 for (j = 0; j < 256; j++)
2568 png_uint_32 ig = (j << (8-shift)) + i;
2569 # ifdef PNG_FLOATING_ARITHMETIC_SUPPORTED
2570 /* Inline the 'max' scaling operation: */
2571 double d = floor(65535*pow(ig/(double)max, gamma_val*.00001)+.5);
2572 sub_table[j] = (png_uint_16)d;
2575 ig = (ig * 65535U + max_by_2)/max;
2577 sub_table[j] = png_gamma_16bit_correct(ig, gamma_val);
2583 /* We must still build a table, but do it the fast way. */
2586 for (j = 0; j < 256; j++)
2588 png_uint_32 ig = (j << (8-shift)) + i;
2591 ig = (ig * 65535U + max_by_2)/max;
2593 sub_table[j] = (png_uint_16)ig;
2599 /* NOTE: this function expects the *inverse* of the overall gamma transformation
2603 png_build_16to8_table(png_structp png_ptr, png_uint_16pp *ptable,
2604 PNG_CONST unsigned int shift, PNG_CONST png_fixed_point gamma_val)
2606 PNG_CONST unsigned int num = 1U << (8U - shift);
2607 PNG_CONST unsigned int max = (1U << (16U - shift))-1U;
2611 png_uint_16pp table = *ptable =
2612 (png_uint_16pp)png_calloc(png_ptr, num * png_sizeof(png_uint_16p));
2614 /* 'num' is the number of tables and also the number of low bits of the
2615 * input 16-bit value used to select a table. Each table is itself indexed
2616 * by the high 8 bits of the value.
2618 for (i = 0; i < num; i++)
2619 table[i] = (png_uint_16p)png_malloc(png_ptr,
2620 256 * png_sizeof(png_uint_16));
2622 /* 'gamma_val' is set to the reciprocal of the value calculated above, so
2623 * pow(out,g) is an *input* value. 'last' is the last input value set.
2625 * In the loop 'i' is used to find output values. Since the output is
2626 * 8-bit there are only 256 possible values. The tables are set up to
2627 * select the closest possible output value for each input by finding
2628 * the input value at the boundary between each pair of output values
2629 * and filling the table up to that boundary with the lower output
2632 * The boundary values are 0.5,1.5..253.5,254.5. Since these are 9-bit
2633 * values the code below uses a 16-bit value in i; the values start at
2634 * 128.5 (for 0.5) and step by 257, for a total of 254 values (the last
2635 * entries are filled with 255). Start i at 128 and fill all 'last'
2636 * table entries <= 'max'
2639 for (i = 0; i < 255; ++i) /* 8-bit output value */
2641 /* Find the corresponding maximum input value */
2642 png_uint_16 out = (png_uint_16)(i * 257U); /* 16-bit output value */
2644 /* Find the boundary value in 16 bits: */
2645 png_uint_32 bound = png_gamma_16bit_correct(out+128U, gamma_val);
2647 /* Adjust (round) to (16-shift) bits: */
2648 bound = (bound * max + 32768U)/65535U + 1U;
2650 while (last < bound)
2652 table[last & (0xffU >> shift)][last >> (8U - shift)] = out;
2657 /* And fill in the final entries. */
2658 while (last < (num << 8))
2660 table[last & (0xff >> shift)][last >> (8U - shift)] = 65535U;
2665 /* Build a single 8-bit table: same as the 16-bit case but much simpler (and
2666 * typically much faster). Note that libpng currently does no sBIT processing
2667 * (apparently contrary to the spec) so a 256-entry table is always generated.
2670 png_build_8bit_table(png_structp png_ptr, png_bytepp ptable,
2671 PNG_CONST png_fixed_point gamma_val)
2674 png_bytep table = *ptable = (png_bytep)png_malloc(png_ptr, 256);
2676 if (png_gamma_significant(gamma_val)) for (i=0; i<256; i++)
2677 table[i] = png_gamma_8bit_correct(i, gamma_val);
2679 else for (i=0; i<256; ++i)
2680 table[i] = (png_byte)i;
2683 /* Used from png_read_destroy and below to release the memory used by the gamma
2687 png_destroy_gamma_table(png_structp png_ptr)
2689 png_free(png_ptr, png_ptr->gamma_table);
2690 png_ptr->gamma_table = NULL;
2692 if (png_ptr->gamma_16_table != NULL)
2695 int istop = (1 << (8 - png_ptr->gamma_shift));
2696 for (i = 0; i < istop; i++)
2698 png_free(png_ptr, png_ptr->gamma_16_table[i]);
2700 png_free(png_ptr, png_ptr->gamma_16_table);
2701 png_ptr->gamma_16_table = NULL;
2704 #if defined(PNG_READ_BACKGROUND_SUPPORTED) || \
2705 defined(PNG_READ_ALPHA_MODE_SUPPORTED) || \
2706 defined(PNG_READ_RGB_TO_GRAY_SUPPORTED)
2707 png_free(png_ptr, png_ptr->gamma_from_1);
2708 png_ptr->gamma_from_1 = NULL;
2709 png_free(png_ptr, png_ptr->gamma_to_1);
2710 png_ptr->gamma_to_1 = NULL;
2712 if (png_ptr->gamma_16_from_1 != NULL)
2715 int istop = (1 << (8 - png_ptr->gamma_shift));
2716 for (i = 0; i < istop; i++)
2718 png_free(png_ptr, png_ptr->gamma_16_from_1[i]);
2720 png_free(png_ptr, png_ptr->gamma_16_from_1);
2721 png_ptr->gamma_16_from_1 = NULL;
2723 if (png_ptr->gamma_16_to_1 != NULL)
2726 int istop = (1 << (8 - png_ptr->gamma_shift));
2727 for (i = 0; i < istop; i++)
2729 png_free(png_ptr, png_ptr->gamma_16_to_1[i]);
2731 png_free(png_ptr, png_ptr->gamma_16_to_1);
2732 png_ptr->gamma_16_to_1 = NULL;
2734 #endif /* READ_BACKGROUND || READ_ALPHA_MODE || RGB_TO_GRAY */
2737 /* We build the 8- or 16-bit gamma tables here. Note that for 16-bit
2738 * tables, we don't make a full table if we are reducing to 8-bit in
2739 * the future. Note also how the gamma_16 tables are segmented so that
2740 * we don't need to allocate > 64K chunks for a full 16-bit table.
2743 png_build_gamma_table(png_structp png_ptr, int bit_depth)
2745 png_debug(1, "in png_build_gamma_table");
2747 /* Remove any existing table; this copes with multiple calls to
2748 * png_read_update_info. The warning is because building the gamma tables
2749 * multiple times is a performance hit - it's harmless but the ability to call
2750 * png_read_update_info() multiple times is new in 1.5.6 so it seems sensible
2751 * to warn if the app introduces such a hit.
2753 if (png_ptr->gamma_table != NULL || png_ptr->gamma_16_table != NULL)
2755 png_warning(png_ptr, "gamma table being rebuilt");
2756 png_destroy_gamma_table(png_ptr);
2761 png_build_8bit_table(png_ptr, &png_ptr->gamma_table,
2762 png_ptr->screen_gamma > 0 ? png_reciprocal2(png_ptr->gamma,
2763 png_ptr->screen_gamma) : PNG_FP_1);
2765 #if defined(PNG_READ_BACKGROUND_SUPPORTED) || \
2766 defined(PNG_READ_ALPHA_MODE_SUPPORTED) || \
2767 defined(PNG_READ_RGB_TO_GRAY_SUPPORTED)
2768 if (png_ptr->transformations & (PNG_COMPOSE | PNG_RGB_TO_GRAY))
2770 png_build_8bit_table(png_ptr, &png_ptr->gamma_to_1,
2771 png_reciprocal(png_ptr->gamma));
2773 png_build_8bit_table(png_ptr, &png_ptr->gamma_from_1,
2774 png_ptr->screen_gamma > 0 ? png_reciprocal(png_ptr->screen_gamma) :
2775 png_ptr->gamma/* Probably doing rgb_to_gray */);
2777 #endif /* READ_BACKGROUND || READ_ALPHA_MODE || RGB_TO_GRAY */
2781 png_byte shift, sig_bit;
2783 if (png_ptr->color_type & PNG_COLOR_MASK_COLOR)
2785 sig_bit = png_ptr->sig_bit.red;
2787 if (png_ptr->sig_bit.green > sig_bit)
2788 sig_bit = png_ptr->sig_bit.green;
2790 if (png_ptr->sig_bit.blue > sig_bit)
2791 sig_bit = png_ptr->sig_bit.blue;
2794 sig_bit = png_ptr->sig_bit.gray;
2796 /* 16-bit gamma code uses this equation:
2798 * ov = table[(iv & 0xff) >> gamma_shift][iv >> 8]
2800 * Where 'iv' is the input color value and 'ov' is the output value -
2803 * Thus the gamma table consists of up to 256 256-entry tables. The table
2804 * is selected by the (8-gamma_shift) most significant of the low 8 bits of
2805 * the color value then indexed by the upper 8 bits:
2807 * table[low bits][high 8 bits]
2809 * So the table 'n' corresponds to all those 'iv' of:
2811 * <all high 8-bit values><n << gamma_shift>..<(n+1 << gamma_shift)-1>
2814 if (sig_bit > 0 && sig_bit < 16U)
2815 shift = (png_byte)(16U - sig_bit); /* shift == insignificant bits */
2818 shift = 0; /* keep all 16 bits */
2820 if (png_ptr->transformations & (PNG_16_TO_8 | PNG_SCALE_16_TO_8))
2822 /* PNG_MAX_GAMMA_8 is the number of bits to keep - effectively
2823 * the significant bits in the *input* when the output will
2824 * eventually be 8 bits. By default it is 11.
2826 if (shift < (16U - PNG_MAX_GAMMA_8))
2827 shift = (16U - PNG_MAX_GAMMA_8);
2831 shift = 8U; /* Guarantees at least one table! */
2833 png_ptr->gamma_shift = shift;
2835 #ifdef PNG_16BIT_SUPPORTED
2836 /* NOTE: prior to 1.5.4 this test used to include PNG_BACKGROUND (now
2837 * PNG_COMPOSE). This effectively smashed the background calculation for
2838 * 16-bit output because the 8-bit table assumes the result will be reduced
2841 if (png_ptr->transformations & (PNG_16_TO_8 | PNG_SCALE_16_TO_8))
2843 png_build_16to8_table(png_ptr, &png_ptr->gamma_16_table, shift,
2844 png_ptr->screen_gamma > 0 ? png_product2(png_ptr->gamma,
2845 png_ptr->screen_gamma) : PNG_FP_1);
2847 #ifdef PNG_16BIT_SUPPORTED
2849 png_build_16bit_table(png_ptr, &png_ptr->gamma_16_table, shift,
2850 png_ptr->screen_gamma > 0 ? png_reciprocal2(png_ptr->gamma,
2851 png_ptr->screen_gamma) : PNG_FP_1);
2854 #if defined(PNG_READ_BACKGROUND_SUPPORTED) || \
2855 defined(PNG_READ_ALPHA_MODE_SUPPORTED) || \
2856 defined(PNG_READ_RGB_TO_GRAY_SUPPORTED)
2857 if (png_ptr->transformations & (PNG_COMPOSE | PNG_RGB_TO_GRAY))
2859 png_build_16bit_table(png_ptr, &png_ptr->gamma_16_to_1, shift,
2860 png_reciprocal(png_ptr->gamma));
2862 /* Notice that the '16 from 1' table should be full precision, however
2863 * the lookup on this table still uses gamma_shift, so it can't be.
2866 png_build_16bit_table(png_ptr, &png_ptr->gamma_16_from_1, shift,
2867 png_ptr->screen_gamma > 0 ? png_reciprocal(png_ptr->screen_gamma) :
2868 png_ptr->gamma/* Probably doing rgb_to_gray */);
2870 #endif /* READ_BACKGROUND || READ_ALPHA_MODE || RGB_TO_GRAY */
2873 #endif /* READ_GAMMA */
2874 #endif /* defined(PNG_READ_SUPPORTED) || defined(PNG_WRITE_SUPPORTED) */