libvmaf/ciede: avoid cross-platform mismatches

moves a few calculations to double precision
error should be below the logging precision now

libvmaf/src/feature/ciede.c

@@ -125,7 +125,7 @@ static float get_h_prime(const float x, const float y)
 {
     if ((x == 0.0) && (y == 0.0))
         return 0.0;
-    float hue_angle = atan2f(x, y);
+    float hue_angle = atan2(x, y);
     if (hue_angle < 0.0)
         hue_angle += 2. * M_PI;
     return hue_angle;
@@ -156,10 +156,10 @@ static float get_upcase_h_bar_prime(const float h_prime_1,
 static float get_upcase_t(const float upcase_h_bar_prime)
 {
     return 1.0 -
-           0.17 * cosf(upcase_h_bar_prime - M_PI / 6.0) +
-           0.24 * cosf(2.0 * upcase_h_bar_prime) +
-           0.32 * cosf(3.0 * upcase_h_bar_prime + M_PI / 30.0) -
-           0.20 * cosf(4.0 * upcase_h_bar_prime - 7.0 * M_PI / 20.0);
+           0.17 * cos(upcase_h_bar_prime - M_PI / 6.0) +
+           0.24 * cos(2.0 * upcase_h_bar_prime) +
+           0.32 * cos(3.0 * upcase_h_bar_prime + M_PI / 30.0) -
+           0.20 * cos(4.0 * upcase_h_bar_prime - 7.0 * M_PI / 20.0);
 }
 
 static float radians_to_degrees(const float radians)
@@ -219,7 +219,7 @@ static float ciede2000(LABColor color_1, LABColor color_2, KSubArgs ksub)
     const float h_prime_2 = get_h_prime(color_2.b, a_prime_2);
     const float delta_h_prime = get_delta_h_prime(c1, c2, h_prime_1, h_prime_2);
     const float delta_upcase_h_prime =
-            2.0 * sqrtf(c_prime_1 * c_prime_2) * sinf(delta_h_prime / 2.0);
+            2.0 * sqrt(c_prime_1 * c_prime_2) * sin(delta_h_prime / 2.0);
     const float upcase_h_bar_prime =
         get_upcase_h_bar_prime(h_prime_1, h_prime_2);
     const float upcase_t = get_upcase_t(upcase_h_bar_prime);
@@ -233,32 +233,32 @@ static float ciede2000(LABColor color_1, LABColor color_2, KSubArgs ksub)
                 pow(hue, 2) + r_sub_t * chroma * hue);
 }
 
-static float pow_2_4(float x)
+static double pow_2_4(double x)
 {
-    return powf(x, 2.4);
+    return pow(x, 2.4);
 }
 
-static float rgb_to_xyz_map(float c)
+static double rgb_to_xyz_map(double c)
 {
     if (c > 10. / 255.) {
-        const float A = 0.055;
-        const float D = 1.0 / 1.055;
+        const double A = 0.055;
+        const double D = 1.0 / 1.055;
         return pow_2_4((c + A) * D);
     } else {
-        const float D = 1.0 / 12.92;
+        const double D = 1.0 / 12.92;
         return (c * D);
     }
 }
 
-static float cbrt_approx(float c)
+static double cbrt_approx(double c)
 {
-    return powf(c, 1.0 / 3.0);
+    return pow(c, 1.0 / 3.0);
 }
 
-static float xyz_to_lab_map(float c)
+static float xyz_to_lab_map(double c)
 {
-    const float KAPPA = 24389.0 / 27.0;
-    const float EPSILON = 216.0 / 24389.0;
+    const double KAPPA = 24389.0 / 27.0;
+    const double EPSILON = 216.0 / 24389.0;
 
     if (c > EPSILON) {
         return cbrt_approx(c);
@@ -267,28 +267,28 @@ static float xyz_to_lab_map(float c)
     }
 }
 
-static LABColor get_lab_color(float y, float u, float v, unsigned bpc)
+static LABColor get_lab_color(double y, double u, double v, unsigned bpc)
 {
-    const float scale = 1 << (bpc - 8);
+    const double scale = 1 << (bpc - 8);
 
     y = (y - 16.  * scale) * (1. / (219. * scale));
     u = (u - 128. * scale) * (1. / (224. * scale));
     v = (v - 128. * scale) * (1. / (224. * scale));
 
     // Assumes BT.709
-    float r = y + 1.28033 * v;
-    float g = y - 0.21482 * u - 0.38059 * v;
-    float b = y + 2.12798 * u;
+    double r = y + 1.28033 * v;
+    double g = y - 0.21482 * u - 0.38059 * v;
+    double b = y + 2.12798 * u;
 
     r = rgb_to_xyz_map(r);
     g = rgb_to_xyz_map(g);
     b = rgb_to_xyz_map(b);
 
-    float x = r * 0.4124564390896921 + g * 0.357576077643909 +
+    double x = r * 0.4124564390896921 + g * 0.357576077643909 +
               b * 0.18043748326639894;
           y = r * 0.21267285140562248 + g * 0.715152155287818 +
               b * 0.07217499330655958;
-    float z = r * 0.019333895582329317 + g * 0.119192025881303 +
+    double z = r * 0.019333895582329317 + g * 0.119192025881303 +
               b * 0.9503040785363677;
 
     x = xyz_to_lab_map(x * (1.0 / 0.95047));

libvmaf/test/meson.build

@@ -87,6 +87,12 @@ test_feature = executable('test_feature',
     include_directories : [libvmaf_inc, test_inc, include_directories('../src/')],
 )
 
+test_ciede = executable('test_ciede',
+    ['test.c', 'test_ciede.c'],
+    include_directories : [libvmaf_inc, test_inc, include_directories('../src/')],
+    link_with : get_option('default_library') == 'both' ? libvmaf.get_static_lib() : libvmaf,
+)
+
 test('test_picture', test_picture)
 test('test_feature_collector', test_feature_collector)
 test('test_thread_pool', test_thread_pool)
@@ -97,3 +103,4 @@ test('test_dict', test_dict)
 test('test_cpu', test_cpu)
 test('test_ref', test_ref)
 test('test_feature', test_feature)
+test('test_ciede', test_ciede)

libvmaf/test/test_ciede.c

@@ -0,0 +1,85 @@
+/**
+ *
+ *  Copyright 2016-2020 Netflix, Inc.
+ *
+ *     Licensed under the BSD+Patent License (the "License");
+ *     you may not use this file except in compliance with the License.
+ *     You may obtain a copy of the License at
+ *
+ *         https://opensource.org/licenses/BSDplusPatent
+ *
+ *     Unless required by applicable law or agreed to in writing, software
+ *     distributed under the License is distributed on an "AS IS" BASIS,
+ *     WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *     See the License for the specific language governing permissions and
+ *     limitations under the License.
+ *
+ */
+
+#include "test.h"
+#include "feature/ciede.c"
+
+static int close_enough(float a, float b)
+{
+    const float epsilon = 1e-9f;
+    return fabs(a - b) < epsilon;
+}
+
+static const KSubArgs default_ksub = { .l = 0.65, .c = 1.0, .h = 4.0 };
+
+static char *test_ciede()
+{
+    const LABColor color_1 = { .l = 0.052488625, .a = -0.587470829, .b = -8.98771572 };
+    const LABColor color_2 = { .l = 0.465437293, .a = 0.386364758, .b = -12.7648535 };
+
+    const float de00 = ciede2000(color_1, color_2, default_ksub);
+    mu_assert("de00 for this input should be 2.54780269",
+              close_enough(de00, 2.54780269));
+
+    return NULL;
+}
+
+static char *test_ciede2()
+{
+    const LABColor color_1 = { .l = 87.156334, .a = -12.049645, .b = -1.205325 };
+    const LABColor color_2 = { .l = 83.455727, .a = -9.040445, .b = -8.894289 };
+
+    const float de00 = ciede2000(color_1, color_2, default_ksub);
+    mu_assert("de00 for this input should be 4.22714281",
+              close_enough(de00, 4.22714281));
+
+    return NULL;
+}
+
+static char *test_ciede3()
+{
+    const LABColor color_1 = { .l = 79.718491, .a = 9.109915, .b = 13.727915 };
+    const LABColor color_2 = { .l = 78.717224, .a = 7.526546, .b = 5.597448 };
+
+    const float de00 = ciede2000(color_1, color_2, default_ksub);
+    mu_assert("de00 for this input should be 4.26012468",
+              close_enough(de00, 4.26012468));
+
+    return NULL;
+}
+
+static char *test_ciede4()
+{
+    const LABColor color_1 = { .l = 99.205299, .a = -3.339410, .b = 1.205873 };
+    const LABColor color_2 = { .l = 97.991730, .a = -2.497345, .b = 2.473533 };
+
+    const float de00 = ciede2000(color_1, color_2, default_ksub);
+    mu_assert("de00 for this input should be 1.26915979",
+              close_enough(de00, 1.26915979));
+
+    return NULL;
+}
+
+char *run_tests()
+{
+    mu_run_test(test_ciede);
+    mu_run_test(test_ciede2);
+    mu_run_test(test_ciede3);
+    mu_run_test(test_ciede4);
+    return NULL;
+}