Add a function to convert gain map metadata double values to fractions. (#1546)
Continued fractions (https://en.wikipedia.org/wiki/Continued_fraction) are used to find the best fractional approximation of a double value.
diff --git a/CMakeLists.txt b/CMakeLists.txt
index c893668..378cbf6 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -295,6 +295,9 @@
src/utils.c
src/write.c
)
+if(AVIF_ENABLE_EXPERIMENTAL_GAIN_MAP)
+ list(APPEND AVIF_SRCS src/gainmap.c)
+endif()
# Only applicable to macOS. In GitHub CI's macos-latest os image, this prevents using the libpng
# and libjpeg headers from /Library/Frameworks/Mono.framework/Headers instead of
@@ -709,6 +712,9 @@
endif()
if(NOT SKIP_INSTALL_HEADERS AND NOT SKIP_INSTALL_ALL)
install(FILES include/avif/avif.h DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}/avif")
+ if(AVIF_ENABLE_EXPERIMENTAL_GAIN_MAP)
+ install(FILES include/avif/gainmap.h DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}/avif")
+ endif()
endif()
# ---------------------------------------------------------------------------------------
diff --git a/include/avif/avif.h b/include/avif/avif.h
index 41d908b..1365f2e 100644
--- a/include/avif/avif.h
+++ b/include/avif/avif.h
@@ -522,7 +522,8 @@
// avifGainMap
// Gain Maps are a HIGHLY EXPERIMENTAL FEATURE. The format might still change and
// images containing a gain map encoded with the current version of libavif might
-// not decode with a feature future version of libavif. Use are your own risk.
+// not decode with a future version of libavif. The API is not guaranteed
+// to be stable, and might even be removed in the future. Use are your own risk.
// This is based on ISO/IEC JTC 1/SC 29/WG 3 m64379
// This product includes Gain Map technology under license by Adobe.
diff --git a/include/avif/gainmap.h b/include/avif/gainmap.h
new file mode 100644
index 0000000..b3b0376
--- /dev/null
+++ b/include/avif/gainmap.h
@@ -0,0 +1,44 @@
+// Copyright 2023 Google LLC
+// SPDX-License-Identifier: BSD-2-Clause
+
+// Utilities to deal with gain maps.
+// This API is experimental and may change or be removed in the future.
+
+#ifndef AVIF_GAINMAP_H
+#define AVIF_GAINMAP_H
+
+#include "avif/avif.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#if defined(AVIF_ENABLE_EXPERIMENTAL_GAIN_MAP)
+
+// Same as avifGainMapMetadata, but with fields of type double instead of uint32_t fractions.
+// Use avifGainMapMetadataDoubleToFractions() to convert this to a avifGainMapMetadata.
+// See avifGainMapMetadata in avif.h for detailed descriptions of fields.
+typedef struct avifGainMapMetadataDouble
+{
+ double gainMapMin[3];
+ double gainMapMax[3];
+ double gainMapGamma[3];
+ double offsetSdr[3];
+ double offsetHdr[3];
+ double hdrCapacityMin;
+ double hdrCapacityMax;
+ avifBool baseRenditionIsHDR;
+} avifGainMapMetadataDouble;
+
+// Converts a avifGainMapMetadataDouble to avifGainMapMetadata by converting double values
+// to the closest uint32_t fractions.
+// Returns AVIF_FALSE if some field values are < 0 or > UINT32_MAX.
+avifBool avifGainMapMetadataDoubleToFractions(avifGainMapMetadata * dst, const avifGainMapMetadataDouble * src);
+
+#endif // AVIF_ENABLE_EXPERIMENTAL_GAIN_MAP
+
+#ifdef __cplusplus
+} // extern "C"
+#endif
+
+#endif // ifndef AVIF_GAINMAP_H
diff --git a/include/avif/internal.h b/include/avif/internal.h
index a70b9d6..4d36701 100644
--- a/include/avif/internal.h
+++ b/include/avif/internal.h
@@ -86,6 +86,10 @@
avifBool avifFractionAdd(avifFraction a, avifFraction b, avifFraction * result);
avifBool avifFractionSub(avifFraction a, avifFraction b, avifFraction * result);
+// Creates a uint32 fraction that is approximately equal to 'v'.
+// Returns AVIF_FALSE if 'v' is < 0 or > UINT32_MAX or NaN.
+avifBool avifToUnsignedFraction(double v, uint32_t * numerator, uint32_t * denominator);
+
void avifImageSetDefaults(avifImage * image);
// Copies all fields that do not need to be freed/allocated from srcImage to dstImage.
void avifImageCopyNoAlloc(avifImage * dstImage, const avifImage * srcImage);
diff --git a/src/gainmap.c b/src/gainmap.c
new file mode 100644
index 0000000..4d0200c
--- /dev/null
+++ b/src/gainmap.c
@@ -0,0 +1,25 @@
+// Copyright 2023 Google LLC
+// SPDX-License-Identifier: BSD-2-Clause
+
+#include "avif/gainmap.h"
+#include "avif/avif.h"
+#include "avif/internal.h"
+
+#if defined(AVIF_ENABLE_EXPERIMENTAL_GAIN_MAP)
+
+avifBool avifGainMapMetadataDoubleToFractions(avifGainMapMetadata * dst, const avifGainMapMetadataDouble * src)
+{
+ for (int i = 0; i < 3; ++i) {
+ AVIF_CHECK(avifToUnsignedFraction(src->gainMapMin[i], &dst->gainMapMinN[i], &dst->gainMapMinD[i]));
+ AVIF_CHECK(avifToUnsignedFraction(src->gainMapMax[i], &dst->gainMapMaxN[i], &dst->gainMapMaxD[i]));
+ AVIF_CHECK(avifToUnsignedFraction(src->gainMapGamma[i], &dst->gainMapGammaN[i], &dst->gainMapGammaD[i]));
+ AVIF_CHECK(avifToUnsignedFraction(src->offsetSdr[i], &dst->offsetSdrN[i], &dst->offsetSdrD[i]));
+ AVIF_CHECK(avifToUnsignedFraction(src->offsetHdr[i], &dst->offsetHdrN[i], &dst->offsetHdrD[i]));
+ }
+ AVIF_CHECK(avifToUnsignedFraction(src->hdrCapacityMin, &dst->hdrCapacityMinN, &dst->hdrCapacityMinD));
+ AVIF_CHECK(avifToUnsignedFraction(src->hdrCapacityMax, &dst->hdrCapacityMaxN, &dst->hdrCapacityMaxD));
+ dst->baseRenditionIsHDR = src->baseRenditionIsHDR;
+ return AVIF_TRUE;
+}
+
+#endif // AVIF_ENABLE_EXPERIMENTAL_GAIN_MAP
diff --git a/src/utils.c b/src/utils.c
index deb12bc..64e823b 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -4,6 +4,7 @@
#include "avif/internal.h"
#include <assert.h>
+#include <float.h>
#include <math.h>
#include <string.h>
@@ -234,3 +235,48 @@
avifFractionSimplify(result);
return AVIF_TRUE;
}
+
+avifBool avifToUnsignedFraction(double v, uint32_t * numerator, uint32_t * denominator)
+{
+ if (isnan(v) || v < 0 || v > UINT32_MAX) {
+ return AVIF_FALSE;
+ }
+
+ // Maximum denominator: makes sure that both the numerator and denominator are <= UINT32_MAX.
+ const uint64_t maxD = (v <= 1) ? UINT32_MAX : (uint64_t)floor(UINT32_MAX / v);
+
+ // Find the best approximation of v as a fraction using continued fractions, see
+ // https://en.wikipedia.org/wiki/Continued_fraction
+ *denominator = 1;
+ uint32_t previousD = 0;
+ double currentV = v - floor(v);
+ int iter = 0;
+ // Set a maximum number of iterations to be safe. Most numbers should
+ // converge in less than ~20 iterations.
+ // The golden ratio is the worst case and takes 39 iterations.
+ const int maxIter = 39;
+ while (iter < maxIter) {
+ const double numeratorDouble = (double)(*denominator) * v;
+ assert(numeratorDouble <= UINT32_MAX);
+ *numerator = (uint32_t)round(numeratorDouble);
+ if (fabs(numeratorDouble - (*numerator)) == 0.0) {
+ return AVIF_TRUE;
+ }
+ currentV = 1.0 / currentV;
+ const double newD = previousD + floor(currentV) * (*denominator);
+ if (newD > maxD) {
+ // This is the best we can do with a denominator <= max_d.
+ return AVIF_TRUE;
+ }
+ previousD = *denominator;
+ assert(newD <= UINT32_MAX);
+ *denominator = (uint32_t)newD;
+ currentV -= floor(currentV);
+ ++iter;
+ }
+ // Maximum number of iterations reached, return what we've found.
+ // For max_iter >= 39 we shouldn't get here. max_iter can be set
+ // to a lower value to speed up the algorithm if needed.
+ *numerator = (uint32_t)round((double)(*denominator) * v);
+ return AVIF_TRUE;
+}
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 1722cee..4af8636 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -209,6 +209,11 @@
target_include_directories(aviftilingtest PRIVATE ${GTEST_INCLUDE_DIRS})
add_test(NAME aviftilingtest COMMAND aviftilingtest)
+ add_executable(avifutilstest gtest/avifutilstest.cc)
+ target_link_libraries(avifutilstest avif_internal ${GTEST_BOTH_LIBRARIES})
+ target_include_directories(avifutilstest PRIVATE ${GTEST_INCLUDE_DIRS})
+ add_test(NAME avifutilstest COMMAND avifutilstest)
+
add_executable(avify4mtest gtest/avify4mtest.cc)
target_link_libraries(avify4mtest aviftest_helpers ${GTEST_BOTH_LIBRARIES})
target_include_directories(avify4mtest PRIVATE ${GTEST_INCLUDE_DIRS})
diff --git a/tests/gtest/avifgainmaptest.cc b/tests/gtest/avifgainmaptest.cc
index 77f31bb..62e5450 100644
--- a/tests/gtest/avifgainmaptest.cc
+++ b/tests/gtest/avifgainmaptest.cc
@@ -1,9 +1,11 @@
// Copyright 2023 Google LLC
// SPDX-License-Identifier: BSD-2-Clause
+#include <cmath>
#include <fstream>
#include "avif/avif.h"
+#include "avif/gainmap.h"
#include "avif/internal.h"
#include "aviftest_helpers.h"
#include "gtest/gtest.h"
@@ -485,5 +487,60 @@
CheckGainMapMetadataMatches(decoded->gainMap.metadata, avifGainMapMetadata());
}
+#define EXPECT_FRACTION_NEAR(numerator, denominator, expected) \
+ EXPECT_NEAR(std::abs((double)numerator / denominator), expected, \
+ expected * 0.001);
+
+TEST(GainMapTest, Convert) {
+ avifGainMapMetadataDouble metadata_double = {};
+ metadata_double.gainMapMin[0] = 1.0;
+ metadata_double.gainMapMin[1] = 1.1;
+ metadata_double.gainMapMin[2] = 1.2;
+ metadata_double.gainMapMax[0] = 10.0;
+ metadata_double.gainMapMax[1] = 10.1;
+ metadata_double.gainMapMax[2] = 10.2;
+ metadata_double.gainMapGamma[0] = 1.0;
+ metadata_double.gainMapGamma[1] = 1.0;
+ metadata_double.gainMapGamma[2] = 1.2;
+ metadata_double.offsetSdr[0] = 1.0 / 32.0;
+ metadata_double.offsetSdr[1] = 1.0 / 64.0;
+ metadata_double.offsetSdr[2] = 1.0 / 128.0;
+ metadata_double.offsetHdr[0] = 0.004564;
+ metadata_double.offsetHdr[1] = 0.0;
+ metadata_double.hdrCapacityMin = 1.0;
+ metadata_double.hdrCapacityMax = 10.0;
+ metadata_double.baseRenditionIsHDR = AVIF_TRUE;
+
+ avifGainMapMetadata metadata = {};
+ ASSERT_TRUE(
+ avifGainMapMetadataDoubleToFractions(&metadata, &metadata_double));
+
+ for (int i = 0; i < 3; ++i) {
+ EXPECT_FRACTION_NEAR(metadata.gainMapMinN[i], metadata.gainMapMinD[i],
+ metadata_double.gainMapMin[i]);
+ EXPECT_FRACTION_NEAR(metadata.gainMapMaxN[i], metadata.gainMapMaxD[i],
+ metadata_double.gainMapMax[i]);
+ EXPECT_FRACTION_NEAR(metadata.gainMapGammaN[i], metadata.gainMapGammaD[i],
+ metadata_double.gainMapGamma[i]);
+ EXPECT_FRACTION_NEAR(metadata.offsetSdrN[i], metadata.offsetSdrD[i],
+ metadata_double.offsetSdr[i]);
+ EXPECT_FRACTION_NEAR(metadata.offsetHdrN[i], metadata.offsetHdrD[i],
+ metadata_double.offsetHdr[i]);
+ }
+ EXPECT_FRACTION_NEAR(metadata.hdrCapacityMinN, metadata.hdrCapacityMinD,
+ metadata_double.hdrCapacityMin);
+ EXPECT_FRACTION_NEAR(metadata.hdrCapacityMaxN, metadata.hdrCapacityMaxD,
+ metadata_double.hdrCapacityMax);
+ EXPECT_EQ(metadata.baseRenditionIsHDR, metadata_double.baseRenditionIsHDR);
+}
+
+TEST(GainMapTest, Invalid) {
+ avifGainMapMetadataDouble metadata_double = {};
+ metadata_double.gainMapGamma[0] = -42; // A negative value is invalid!
+ avifGainMapMetadata metadata = {};
+ ASSERT_FALSE(
+ avifGainMapMetadataDoubleToFractions(&metadata, &metadata_double));
+}
+
} // namespace
} // namespace libavif
diff --git a/tests/gtest/avifutilstest.cc b/tests/gtest/avifutilstest.cc
new file mode 100644
index 0000000..fe342f4
--- /dev/null
+++ b/tests/gtest/avifutilstest.cc
@@ -0,0 +1,129 @@
+// Copyright 2023 Google LLC
+// SPDX-License-Identifier: BSD-2-Clause
+
+#include <cmath>
+
+#include "avif/avif.h"
+#include "avif/internal.h"
+#include "aviftest_helpers.h"
+#include "gtest/gtest.h"
+
+namespace libavif {
+namespace {
+
+// Converts a double value to a fraction, and checks that the difference
+// between numerator/denominator and v is below relative_tolerance.
+void TestRoundTrip(double v, double relative_tolerance) {
+ uint32_t numerator, denominator;
+ ASSERT_TRUE(avifToUnsignedFraction(v, &numerator, &denominator)) << v;
+ const double reconstructed = (double)numerator / denominator;
+ const double tolerance = v * relative_tolerance;
+ EXPECT_NEAR(reconstructed, v, tolerance)
+ << "numerator " << (double)numerator << " denominator "
+ << (double)denominator;
+}
+
+constexpr double kLotsOfDecimals = 0.14159265358979323846;
+
+TEST(ToFractionUTest, RoundTrip) {
+ // Whole numbers and simple fractions should match perfectly.
+ constexpr double kPerfectTolerance = 0.0;
+ TestRoundTrip(0.0, kPerfectTolerance);
+ TestRoundTrip(1.0, kPerfectTolerance);
+ TestRoundTrip(42.0, kPerfectTolerance);
+ TestRoundTrip(102356.0, kPerfectTolerance);
+ TestRoundTrip(102356456.0f, kPerfectTolerance);
+ TestRoundTrip(UINT32_MAX / 2.0, kPerfectTolerance);
+ TestRoundTrip((double)UINT32_MAX - 1.0, kPerfectTolerance);
+ TestRoundTrip((double)UINT32_MAX, kPerfectTolerance);
+ TestRoundTrip(0.123, kPerfectTolerance);
+ TestRoundTrip(1.0 / 3.0, kPerfectTolerance);
+ TestRoundTrip(1.0 / 4.0, kPerfectTolerance);
+ TestRoundTrip(3.0 / 23.0, kPerfectTolerance);
+ TestRoundTrip(1253456.456, kPerfectTolerance);
+ TestRoundTrip(8598533.9, kPerfectTolerance);
+
+ // // Numbers with a lot of decimals or very large/small can show a small
+ // error.
+ constexpr double kSmallTolerance = 1e-9;
+ TestRoundTrip(0.0123456, kSmallTolerance);
+ TestRoundTrip(3 + kLotsOfDecimals, kSmallTolerance);
+ TestRoundTrip(sqrt(2.0), kSmallTolerance);
+ TestRoundTrip(exp(1.0), kSmallTolerance);
+ TestRoundTrip(exp(10.0), kSmallTolerance);
+ TestRoundTrip(exp(15.0), kSmallTolerance);
+ // The golden ratio, the irrational number that is the "most difficult" to
+ // approximate rationally according to Wikipedia.
+ const double kGoldenRatio = (1.0 + std::sqrt(5.0)) / 2.0;
+ TestRoundTrip(kGoldenRatio, kSmallTolerance); // Golden ratio.
+ TestRoundTrip(((double)UINT32_MAX) - 0.5, kSmallTolerance);
+ // Note that values smaller than this might have a larger relative error
+ // (e.g. 1.0e-10).
+ TestRoundTrip(4.2e-10, kSmallTolerance);
+}
+
+// Tests the max difference between the fraction-ified value and the original
+// value, for a subset of values between 0.0 and UINT32_MAX.
+TEST(ToFractionUTest, MaxDifference) {
+ double max_error = 0;
+ double max_error_v = 0;
+ double max_relative_error = 0;
+ double max_relative_error_v = 0;
+ for (uint64_t i = 0; i < UINT32_MAX; i += 1000) {
+ const double v = i + kLotsOfDecimals;
+ uint32_t numerator, denominator;
+ ASSERT_TRUE(avifToUnsignedFraction(v, &numerator, &denominator)) << v;
+ const double reconstructed = (double)numerator / denominator;
+ const double error = abs(reconstructed - v);
+ const double relative_error = error / v;
+ if (error > max_error) {
+ max_error = error;
+ max_error_v = v;
+ }
+ if (relative_error > max_relative_error) {
+ max_relative_error = relative_error;
+ max_relative_error_v = v;
+ }
+ }
+ EXPECT_LE(max_error, 0.5f) << max_error_v;
+ EXPECT_LT(max_relative_error, 1e-9) << max_relative_error_v;
+}
+
+// Tests the max difference between the fraction-ified value and the original
+// value, for a subset of values between 0 and 1.0/UINT32_MAX.
+TEST(ToFractionUTest, MaxDifferenceSmall) {
+ double max_error = 0;
+ double max_error_v = 0;
+ double max_relative_error = 0;
+ double max_relative_error_v = 0;
+ for (uint64_t i = 1; i < UINT32_MAX; i += 1000) {
+ const double v = 1.0 / (i + kLotsOfDecimals);
+ uint32_t numerator, denominator;
+ ASSERT_TRUE(avifToUnsignedFraction(v, &numerator, &denominator)) << v;
+ const double reconstructed = (double)numerator / denominator;
+ const double error = abs(reconstructed - v);
+ const double relative_error = error / v;
+ if (error > max_error) {
+ max_error = error;
+ max_error_v = v;
+ }
+ if (relative_error > max_relative_error) {
+ max_relative_error = relative_error;
+ max_relative_error_v = v;
+ }
+ }
+ EXPECT_LE(max_error, 1e-10) << max_error_v;
+ EXPECT_LT(max_relative_error, 1e-5) << max_relative_error_v;
+}
+
+TEST(ToFractionUTest, BadValues) {
+ uint32_t numerator, denominator;
+ // Negative value.
+ EXPECT_FALSE(avifToUnsignedFraction(-0.1, &numerator, &denominator));
+ // Too large.
+ EXPECT_FALSE(avifToUnsignedFraction(((double)UINT32_MAX) + 1.0, &numerator,
+ &denominator));
+}
+
+} // namespace
+} // namespace libavif
