merge metrics to statistics, refactor

elsa/core/CMakeLists.txt

@@ -10,8 +10,7 @@ set(MODULE_HEADERS
     Utilities/DataContainerFormatter.hpp
     Utilities/FormatConfig.h
     Utilities/Math.hpp
-    Utilities/Metrics.hpp
-    Utilities/Statistics.hpp
+    Utilities/Statistics.h
     Utilities/TypeCasts.hpp
     Descriptors/DataDescriptor.h
     Descriptors/DescriptorUtils.h
@@ -44,6 +43,7 @@ set(MODULE_SOURCES
     Utilities/Assertions.cpp
     Descriptors/DataDescriptor.cpp
     Descriptors/VolumeDescriptor.cpp
+    Utilities/Statistics.cpp
     Descriptors/PlanarDetectorDescriptor.cpp
     Descriptors/RandomBlocksDescriptor.cpp
     Descriptors/DescriptorUtils.cpp

elsa/core/Utilities/Metrics.hpp

-#pragma once
-
-#include "elsaDefines.h"
-#include "DataContainer.h"
-#include "Statistics.hpp"
-
-#include <cmath>
-#include <numeric>
-#include <algorithm>
-#include <utility>
-
-namespace elsa
-{
-    /**
-     * @brief Compute the Relative Error between two given signals.
-     *
-     * @param dc1 DataContainer signal
-     * @param dc2 DataContainer signal
-     */
-    template <typename data_t = real_t>
-    constexpr auto relativeError(DataContainer<data_t> dc1, DataContainer<data_t> dc2)
-        -> long double
-    {
-        if (dc1.getDataDescriptor() != dc2.getDataDescriptor()) {
-            throw LogicError(
-                std::string("Metrics::relativeError: shapes of both signals should match"));
-        }
-
-        DataContainer<data_t> diff = dc1 - dc2;
-        return diff.l2Norm() / dc2.l2Norm();
-    }
-
-    /**
-     * @brief Compute the Peak Signal-to-Noise Ratio.
-     *
-     * @param dc1 DataContainer signal
-     * @param dc2 DataContainer signal
-     * @param dataRange The data range of the signals (distance between minimum and maximum possible
-     * values).
-     */
-    template <typename data_t = real_t>
-    constexpr auto peakSignalToNoiseRatio(DataContainer<data_t> dc1, DataContainer<data_t> dc2,
-                                          data_t dataRange) -> long double
-    {
-        if (dc1.getDataDescriptor() != dc2.getDataDescriptor()) {
-            throw LogicError(std::string(
-                "Metrics::peakSignalToNoiseRatio: shapes of both signals should match"));
-        }
-
-        long double err = meanSquaredError<data_t>(dc1, dc2);
-        return 10 * std::log10l((std::pow(dataRange, 2) / err));
-    }
-
-    /**
-     * @brief Compute the Peak Signal-to-Noise Ratio.
-     *
-     * @param dc1 DataContainer signal
-     * @param dc2 DataContainer signal
-     */
-    template <typename data_t = real_t>
-    constexpr auto peakSignalToNoiseRatio(DataContainer<data_t> dc1, DataContainer<data_t> dc2)
-        -> long double
-    {
-        if (dc1.getDataDescriptor() != dc2.getDataDescriptor()) {
-            throw LogicError(std::string(
-                "Metrics::peakSignalToNoiseRatio: shapes of both signals should match"));
-        }
-
-        data_t dataMin = std::numeric_limits<data_t>::min();
-        data_t dataMax = std::numeric_limits<data_t>::max();
-
-        data_t trueMin = std::min(minElement(dc1), minElement(dc2));
-        data_t trueMax = std::max(maxElement(dc1), maxElement(dc2));
-
-        if (trueMin < dataMin || trueMax > dataMax) {
-            throw LogicError(
-                std::string("Metrics::peakSignalToNoiseRatio: extreme values of the signals "
-                            "exceed the range expected for its data type"));
-        }
-
-        data_t dataRange;
-        if (trueMin >= 0) {
-            dataRange = dataMax;
-        } else {
-            dataRange = dataMax - dataMin;
-        }
-
-        return peakSignalToNoiseRatio(dc1, dc2, dataRange);
-    }
-} // namespace elsa

elsa/core/Utilities/Statistics.hpp → elsa/core/Utilities/Statistics.cpp

-#pragma once
-
-#include "elsaDefines.h"
-#include "DataContainer.h"
-
-#include <cmath>
-#include <numeric>
-#include <algorithm>
-#include <utility>
+#include "Statistics.h"
 
 namespace elsa
 {
-    /**
-     * @brief Compute the Mean Squared Error of two given signals.
-     *
-     * @param dc1 DataContainer signal
-     * @param dc2 DataContainer signal
-     */
-    template <typename data_t = real_t>
-    constexpr auto meanSquaredError(DataContainer<data_t> dc1, DataContainer<data_t> dc2)
-        -> long double
+    template <typename data_t>
+    data_t Statistics<data_t>::meanSquaredError(DataContainer<data_t> dc1,
+                                                DataContainer<data_t> dc2)
     {
         if (dc1.getDataDescriptor() != dc2.getDataDescriptor()) {
             throw LogicError(
@@ -26,16 +12,12 @@ namespace elsa
         }
 
         DataContainer<data_t> diff = dc1 - dc2;
-        DataContainer<data_t> squared = DataContainer<data_t>{square(diff)};
-        return squared.sum() / dc1.getSize();
+        return diff.squaredL2Norm() / dc1.getSize();
     }
 
-    /// @brief Calculate mean and standard deviation of a container
-    /// @param v Any container, such as `std::vector`
-    /// @tparam Container Container type of argument (e.g. `std::vector`)
-    /// @return a pair of mean and standard deviation (of type `Conatiner::value_type`)
+    template <typename data_t>
     template <typename Container>
-    constexpr auto calculateMeanStddev(Container v)
+    constexpr auto Statistics<data_t>::calculateMeanStddev(Container v)
         -> std::pair<typename Container::value_type, typename Container::value_type>
     {
         // value type of the vector
@@ -57,22 +39,9 @@ namespace elsa
         return std::make_pair(mean, stdev);
     }
 
-    /**
-     * @brief Compute the 95% confidence interval for a given number of samples `n`
-     * and the mean and standard deviation of the measurements.
-     *
-     * Compute it as \f$mean - c(n) * stddev, mean + c(n) * stddev\f$, where
-     * \f$c(n)\f$ is the n-th entry in the two tails T distribution table. For \f$n > 30\f$,
-     * it is assumed that \f$n = 1.96\f$.
-     *
-     * @param n Number of of samples
-     * @param mean mean of samples
-     * @param stddev standard deviation of samples
-     * @return pair of lower and upper bound of 95% confidence interval
-     */
-    template <typename data_t = real_t>
-    constexpr auto confidenceInterval95(std::size_t n, data_t mean, data_t stddev)
-        -> std::pair<real_t, real_t>
+    template <typename data_t>
+    std::pair<real_t, real_t> Statistics<data_t>::confidenceInterval95(std::size_t n, data_t mean,
+                                                                       data_t stddev)
     {
         // Do we run often enough to assume a large sample size?
         if (n > 30) {
@@ -96,4 +65,64 @@ namespace elsa
             return std::make_pair(lower, upper);
         }
     }
+
+    template <typename data_t>
+    data_t Statistics<data_t>::relativeError(DataContainer<data_t> dc1, DataContainer<data_t> dc2)
+    {
+        if (dc1.getDataDescriptor() != dc2.getDataDescriptor()) {
+            throw InvalidArgumentError(
+                std::string("Metrics::relativeError: shapes of both signals should match"));
+        }
+
+        DataContainer<data_t> diff = dc1 - dc2;
+        return diff.l2Norm() / dc2.l2Norm();
+    }
+
+    template <typename data_t>
+    data_t Statistics<data_t>::peakSignalToNoiseRatio(DataContainer<data_t> dc1,
+                                                      DataContainer<data_t> dc2, data_t dataRange)
+    {
+        if (dc1.getDataDescriptor() != dc2.getDataDescriptor()) {
+            throw InvalidArgumentError(std::string(
+                "Metrics::peakSignalToNoiseRatio: shapes of both signals should match"));
+        }
+
+        return data_t(10) * std::log10((std::pow(dataRange, 2) / meanSquaredError(dc1, dc2)));
+    }
+
+    template <typename data_t>
+    data_t Statistics<data_t>::peakSignalToNoiseRatio(DataContainer<data_t> dc1,
+                                                      DataContainer<data_t> dc2)
+    {
+        if (dc1.getDataDescriptor() != dc2.getDataDescriptor()) {
+            throw InvalidArgumentError(std::string(
+                "Metrics::peakSignalToNoiseRatio: shapes of both signals should match"));
+        }
+
+        data_t dataMin = std::numeric_limits<data_t>::min();
+        data_t dataMax = std::numeric_limits<data_t>::max();
+
+        data_t trueMin = std::min(dc1.minElement(), dc2.minElement());
+        data_t trueMax = std::max(dc1.maxElement(), dc2.maxElement());
+
+        if (trueMin < dataMin || trueMax > dataMax) {
+            throw InvalidArgumentError(
+                std::string("Metrics::peakSignalToNoiseRatio: extreme values of the signals "
+                            "exceed the range expected for its data type"));
+        }
+
+        data_t dataRange;
+        if (trueMin >= 0) {
+            dataRange = dataMax;
+        } else {
+            dataRange = dataMax - dataMin;
+        }
+
+        return peakSignalToNoiseRatio(dc1, dc2, dataRange);
+    }
+
+    // ------------------------------------------
+    // explicit template instantiation
+    template class Statistics<float>;
+    template class Statistics<double>;
 } // namespace elsa

elsa/core/Utilities/Statistics.h

+#pragma once
+
+#include "elsaDefines.h"
+#include "DataContainer.h"
+
+#include <cmath>
+#include <numeric>
+#include <algorithm>
+#include <utility>
+
+namespace elsa
+{
+    /**
+     * @brief Class representing ...
+     *
+     * This class represents an ...
+     *
+     * @tparam data_t data type for the domain and range of the problem, defaulting to real_t
+     *
+     * References:
+     */
+    template <typename data_t = real_t>
+    class Statistics
+    {
+        /**
+         * @brief Compute the Mean Squared Error of two given signals.
+         *
+         * @param dc1 DataContainer signal
+         * @param dc2 DataContainer signal
+         */
+        data_t meanSquaredError(DataContainer<data_t> dc1, DataContainer<data_t> dc2);
+
+        /// @brief Calculate mean and standard deviation of a container
+        /// @param v Any container, such as `std::vector`
+        /// @tparam Container Container type of argument (e.g. `std::vector`)
+        /// @return a pair of mean and standard deviation (of type `Container::value_type`)
+        template <typename Container>
+        constexpr auto calculateMeanStddev(Container v)
+            -> std::pair<typename Container::value_type, typename Container::value_type>;
+
+        /**
+         * @brief Compute the 95% confidence interval for a given number of samples `n`
+         * and the mean and standard deviation of the measurements.
+         *
+         * Compute it as \f$mean - c(n) * stddev, mean + c(n) * stddev\f$, where
+         * \f$c(n)\f$ is the n-th entry in the two tails T distribution table. For \f$n > 30\f$,
+         * it is assumed that \f$n = 1.96\f$.
+         *
+         * @param n Number of of samples
+         * @param mean mean of samples
+         * @param stddev standard deviation of samples
+         * @return pair of lower and upper bound of 95% confidence interval
+         */
+        std::pair<real_t, real_t> confidenceInterval95(std::size_t n, data_t mean, data_t stddev);
+
+        /**
+         * @brief Compute the Relative Error between two given signals.
+         *
+         * @param dc1 DataContainer signal
+         * @param dc2 DataContainer signal
+         */
+        data_t relativeError(DataContainer<data_t> dc1, DataContainer<data_t> dc2);
+
+        /**
+         * @brief Compute the Peak Signal-to-Noise Ratio.
+         *
+         * @param dc1 DataContainer signal
+         * @param dc2 DataContainer signal
+         * @param dataRange The data range of the signals (distance between minimum and maximum
+         * possible values).
+         */
+        data_t peakSignalToNoiseRatio(DataContainer<data_t> dc1, DataContainer<data_t> dc2,
+                                              data_t dataRange);
+
+        /**
+         * @brief Compute the Peak Signal-to-Noise Ratio.
+         *
+         * @param dc1 DataContainer signal
+         * @param dc2 DataContainer signal
+         */
+        data_t peakSignalToNoiseRatio(DataContainer<data_t> dc1, DataContainer<data_t> dc2);
+    };
+} // namespace elsa

elsa/core/tests/CMakeLists.txt

@@ -30,4 +30,4 @@ ELSA_DOCTEST(DataContainerFormatter)
 ELSA_DOCTEST(CartesianIndices)
 ELSA_DOCTEST(Bessel)
 ELSA_DOCTEST(Math)
-ELSA_DOCTEST(Metrics)
+ELSA_DOCTEST(Statistics)

elsa/core/tests/test_Metrics.cpp → elsa/core/tests/test_Statistics.cpp

 /**
- * @file test_Metrics.cpp
+ * @file test_Statistics.cpp
  *
- * @brief Tests for Metrics header
+ * @brief Tests for Statistics header
  *
  * @author Andi Braimllari
  */
 
-#include "Metrics.hpp"
+#include "Statistics.h"
 #include "DataContainer.h"
 #include "VolumeDescriptor.h"
 
@@ -17,8 +17,8 @@ using namespace elsa;
 using namespace doctest;
 
 TEST_SUITE_BEGIN("core");
-
-TEST_CASE_TEMPLATE("Metrics: Testing the metrics", TestType, float, double)
+// TODO do the Statistics' tests here
+TEST_CASE_TEMPLATE("Statistics: Testing the metrics", TestType, float, double)
 {
     GIVEN("a DataContainer")
     {
@@ -33,6 +33,8 @@ TEST_CASE_TEMPLATE("Metrics: Testing the metrics", TestType, float, double)
 
         WHEN("running the Relative Error")
         {
+            Statistics statistics;
+
             long double relErr = relativeError<TestType>(dataCont1, dataCont2);
             THEN("it produces the correct result")
             {

elsa/test_allmodules/TestDriver.h

 #pragma once
 
-#include "Statistics.hpp"
+#include "Statistics.h"
 #include "SetupHelpers.h"
 #include "NoiseGenerators.h"
 #include "LoggingHelpers.h"

examples/example_argparse.cpp

 #include "elsa.h"
 #include "LutProjector.h"
-#include "Utilities/Statistics.hpp"
+#include "Utilities/Statistics.h"
 #include "IO.h"
 #include "spdlog/fmt/bundled/core.h"