Add log-based mean-square error (LogMSE) to exrmetrics

src/bin/exrmetrics/CMakeLists.txt

@@ -1,4 +1,4 @@
 # SPDX-License-Identifier: BSD-3-Clause
 # Copyright (c) Contributors to the OpenEXR Project.
 
-add_openexr_bin_program(exrmetrics SOURCES main.cpp exrmetrics.cpp)
+add_openexr_bin_program(exrmetrics SOURCES main.cpp exrmetrics.cpp distortionUtils.cpp)

src/bin/exrmetrics/distortionUtils.cpp

@@ -0,0 +1,57 @@
+// SPDX-License-Identifier: BSD-3-Clause
+// Copyright (c) Contributors to the OpenEXR Project.
+
+#include "distortionUtils.h"
+
+#include <cmath>
+#include <limits>
+
+using IMATH_NAMESPACE::half;
+
+template <>
+void
+accumLogMSE<half> (
+    const half* orig,
+    const half* reread,
+    uint64_t    pixelsInChannel,
+    double&     sumSq,
+    uint64_t&   count)
+{
+    static const double LN_HALF_DENORM_MIN = std::log(HALF_DENORM_MIN);
+    for (uint64_t px = 0; px < pixelsInChannel; ++px)
+    {
+        double a = static_cast<double> (orig[px]);
+        double b = static_cast<double> (reread[px]);
+        if (!std::isfinite (a)) { continue; }
+        if (!std::isfinite (b)) { count = 0; sumSq = 0.0; return; }
+
+        double diff = (a < 0 ? -1.0 : 1.0) * (std::log (std::abs (a) + HALF_DENORM_MIN) - LN_HALF_DENORM_MIN) -
+                        (b < 0 ? -1.0 : 1.0) * (std::log (std::abs (b) + HALF_DENORM_MIN) - LN_HALF_DENORM_MIN);
+        sumSq += diff * diff;
+        ++count;
+}
+}
+
+template <>
+void
+accumLogMSE<float> (
+    const float* orig,
+    const float* reread,
+    uint64_t     pixelsInChannel,
+    double&      sumSq,
+    uint64_t&    count)
+{
+    static const double ln_eps = std::log (std::numeric_limits<float>::denorm_min ());
+    static constexpr double eps = std::numeric_limits<float>::denorm_min ();
+    for (uint64_t px = 0; px < pixelsInChannel; ++px)
+    {
+        double a = static_cast<double> (orig[px]);
+        double b = static_cast<double> (reread[px]);
+        if (!std::isfinite (a)) { continue; }
+        if (!std::isfinite (b)) { count = 0; sumSq = 0.0; return; }
+        double diff = (a < 0 ? -1.0 : 1.0) * (std::log (std::abs (a) + eps) - ln_eps) -
+                        (b < 0 ? -1.0 : 1.0) * (std::log (std::abs (b) + eps) - ln_eps);
+        sumSq += diff * diff;
+        ++count;
+    }
+}

src/bin/exrmetrics/distortionUtils.h

@@ -0,0 +1,34 @@
+// SPDX-License-Identifier: BSD-3-Clause
+// Copyright (c) Contributors to the OpenEXR Project.
+
+#ifndef INCLUDED_DISTORTION_UTILS_H
+#define INCLUDED_DISTORTION_UTILS_H
+
+#include <half.h>
+#include <cstdint>
+
+template <typename T>
+void accumLogMSE (
+    const T*  orig,
+    const T*  reread,
+    uint64_t  pixelsInChannel,
+    double&   sumSq,
+    uint64_t& count);
+
+template <>
+void accumLogMSE<IMATH_NAMESPACE::half> (
+    const IMATH_NAMESPACE::half* orig,
+    const IMATH_NAMESPACE::half* reread,
+    uint64_t                     pixelsInChannel,
+    double&                      sumSq,
+    uint64_t&                    count);
+
+template <>
+void accumLogMSE<float> (
+    const float* orig,
+    const float* reread,
+    uint64_t     pixelsInChannel,
+    double&      sumSq,
+    uint64_t&    count);
+
+#endif

src/bin/exrmetrics/exrmetrics.cpp

@@ -5,6 +5,7 @@
 //
 
 #include "exrmetrics.h"
+#include "distortionUtils.h"
 
 #include "ImfChannelList.h"
 #include "ImfDeepFrameBuffer.h"
@@ -23,15 +24,20 @@
 #include "ImfTiledMisc.h"
 #include "ImfTiledOutputPart.h"
 
+#include <Imath/half.h>
+
 #include <chrono>
+#include <cmath>
 #include <ctime>
+#include <limits>
 #include <list>
 #include <stdexcept>
 #include <vector>
 #include <sys/stat.h>
 
 using namespace OPENEXR_IMF_NAMESPACE;
 using IMATH_NAMESPACE::Box2i;
+using IMATH_NAMESPACE::half;
 
 using std::cerr;
 using namespace std::chrono;
@@ -988,7 +994,8 @@ exrmetrics (
     bool                               write,
     bool                               reread,
     PixelMode                          pixelMode,
-    bool                               verbose)
+    bool                               verbose,
+    bool                               computeDistortion)
 {
 
     if (verbose)
@@ -1081,7 +1088,7 @@ exrmetrics (
     if (!isinf (level) && level >= -1 && !compressionSet)
     {
         throw runtime_error (
-            "-l option only works for DWAA/DWAB,ZIP/ZIPS or ZSTD compression");
+            "-l option only works for DWAA/DWAB, ZIP/ZIPS or ZSTD compression");
     }
 
     vector<partData> parts (part == -1 ? in.parts () : 1);
@@ -1171,11 +1178,111 @@ exrmetrics (
         else { metrics.outputFileSize = fileSize; }
     }
 
+    //
+    // compute distortion vs. re-read data
+    //
+    if (computeDistortion && write && reread)
+    {
+        for (size_t p = 0; p < parts.size (); ++p)
+        {
+            metrics.stats[p].distortionCount = 0;
+            metrics.stats[p].distortion = std::numeric_limits<double>::quiet_NaN ();
+            string partType = outHeaders[p].type ();
+            if (partType != SCANLINEIMAGE && partType != TILEDIMAGE) continue;
+
+            // skip parts with mixed channel types or with no channels
+            {
+                auto chBegin = outHeaders[p].channels ().begin ();
+                auto chEnd   = outHeaders[p].channels ().end ();
+                if (chBegin == chEnd) continue;
+                PixelType firstType  = chBegin.channel ().type;
+                bool      allSameType = true;
+                for (auto i = chBegin; i != chEnd; ++i)
+                {
+                    if (i.channel ().type != firstType)
+                    {
+                        allSameType = false;
+                        break;
+                    }
+                }
+                if (!allSameType) continue;
+            }
+
+            Box2i    dw     = outHeaders[p].dataWindow ();
+            uint64_t width  = dw.max.x + 1 - dw.min.x;
+            uint64_t height = dw.max.y + 1 - dw.min.y;
+
+            double   sumSq = 0.0;
+            uint64_t count = 0;
+            int      channelIndex = 0;
+
+            for (ChannelList::ConstIterator i =
+                     outHeaders[p].channels ().begin ();
+                 i != outHeaders[p].channels ().end ();
+                 ++i, ++channelIndex)
+            {
+                if (i.channel ().type != HALF && i.channel ().type != FLOAT)
+                    continue;
+
+                uint64_t pixelsInChannel =
+                    (width / i.channel ().xSampling) *
+                    (height / i.channel ().ySampling);
+
+                const char* origData   = nullptr;
+                const char* rereadData = nullptr;
+
+                if (partType == SCANLINEIMAGE)
+                {
+                    origData =
+                        parts[p].readBuf.scanlinePixelData[channelIndex].data ();
+                    rereadData =
+                        parts[p].rereadBuf.scanlinePixelData[channelIndex].data ();
+                }
+                else
+                {
+                    if (parts[p].readBuf.tilePixelData.empty () ||
+                        parts[p].rereadBuf.tilePixelData.empty ())
+                        continue;
+                    origData =
+                        parts[p].readBuf.tilePixelData[0][channelIndex].data ();
+                    rereadData =
+                        parts[p].rereadBuf.tilePixelData[0][channelIndex].data ();
+                }
+
+                if (i.channel ().type == HALF)
+                {
+                    metrics.stats[p].metricKind = LOG_MSE_HALF;
+                    accumLogMSE (
+                        reinterpret_cast<const half*> (origData),
+                        reinterpret_cast<const half*> (rereadData),
+                        pixelsInChannel,
+                        sumSq,
+                        count);
+                }
+                else if (i.channel ().type == FLOAT)
+                {
+                    metrics.stats[p].metricKind = LOG_MSE_FLOAT;
+                    accumLogMSE (
+                        reinterpret_cast<const float*> (origData),
+                        reinterpret_cast<const float*> (rereadData),
+                        pixelsInChannel,
+                        sumSq,
+                        count);
+                }
+            }
+
+            metrics.stats[p].distortionCount = count;
+            metrics.stats[p].distortion =
+                count > 0 ? sumSq / count
+                          : std::numeric_limits<double>::quiet_NaN ();
+        }
+    }
+
     //
     // sum across all parts
     //
 
-    metrics.totalStats = metrics.stats[0];
+    metrics.totalStats      = metrics.stats[0];
     for (size_t i = 1; i < metrics.stats.size (); ++i)
     {
         accumulate (metrics.totalStats.readPerf, metrics.stats[i].readPerf);

src/bin/exrmetrics/exrmetrics.h

@@ -15,6 +15,7 @@
 
 #include "ImfCompression.h"
 
+#include <limits>
 #include <stdint.h>
 
 #include <vector>
@@ -43,6 +44,13 @@ struct partSizeData
     std::string partType = "";
 };
 
+enum DistortionMetric
+{
+    DISTORTION_METRIC_NONE,
+    LOG_MSE_HALF,
+    LOG_MSE_FLOAT,
+};
+
 struct partStats
 {
     std::vector<double>
@@ -57,6 +65,10 @@ struct partStats
         rereadPerf; // for deep, times reading the sample count, otherwise times reading the entire data
     uint64_t sizeOnDisk = 0; // record compressed size of part on disk.
 
+    DistortionMetric  metricKind  = DISTORTION_METRIC_NONE; // kind of distortion metric used
+    double   distortion      = std::numeric_limits<double>::quiet_NaN (); // distortion computed using the distortion metric
+    uint64_t distortionCount = 0; // number of samples used to compute the distortion
+
     partSizeData sizeData;
 };
 
@@ -78,6 +90,7 @@ fileMetrics exrmetrics (
     bool                               write,
     bool                               reread,
     PixelMode                          pixelMode,
-    bool                               verbose);
+    bool                               verbose,
+    bool                               computeDistortion = false);
 
 #endif