Implement nanmean and nansum in C++

benchmarks/benchmark_statistics.py

@@ -0,0 +1,41 @@
+import time
+import pytest
+
+import numpy as np
+
+from extra_foam.algorithms import (
+    nanmean, nansum
+)
+
+
+def benchmark_nan(f_cpp, f_py, shape):
+    data = np.random.randn(*shape)
+    data[::2, ::2] = np.nan
+
+    t0 = time.perf_counter()
+    for i in range(100):
+        ret_cpp = f_cpp(data)
+    dt_cpp = time.perf_counter() - t0
+
+    t0 = time.perf_counter()
+    for i in range(100):
+        ret_py = f_py(data)
+    dt_py = time.perf_counter() - t0
+
+    assert ret_cpp == pytest.approx(ret_py)
+
+    print(f"\n{f_cpp.__name__} with {float} - \n"
+          f"dt (cpp): {dt_cpp:.4f}, "
+          f"dt (numpy): {dt_py:.4f}")
+
+
+if __name__ == "__main__":
+    print("*" * 80)
+    print("Benchmark image processing")
+    print("*" * 80)
+
+    s = (1096, 1120)
+
+    for f_cpp, f_py in [(nanmean, np.nanmean),
+                        (nansum, np.nansum)]:
+        benchmark_nan(f_cpp, f_py, s)

extra_foam/algorithms/statistics.py → extra_foam/algorithms/statistics_py.py

@@ -10,6 +10,7 @@
 import numpy as np
 
 from .imageproc_py import mask_image_data
+from .statistics import nanmean, nansum
 
 
 def find_actual_range(arr, range):

extra_foam/algorithms/tests/test_statistics.py

@@ -10,7 +10,15 @@
 )
 
 
-class TestMiscellaneous:
+class TestStatistics:
+
+    def testNanmean(self):
+        roi = np.array([[np.nan, 1, 2], [3, 6, np.nan]], dtype=np.float32)
+        assert 3 == nanmean(roi)
+
+    def testNansum(self):
+        roi = np.array([[np.nan, 1, 2], [3, 6, np.nan]], dtype=np.float32)
+        assert 12 == nansum(roi)
 
     def testNanhistWithStats(self):
         # case 1

extra_foam/pipeline/processors/image_roi.py

@@ -20,7 +20,9 @@
 from ...utils import profiler
 from ...config import AnalysisType, Normalizer, RoiCombo, RoiFom
 
-from extra_foam.algorithms import intersection, mask_image_data
+from extra_foam.algorithms import (
+    intersection, mask_image_data
+)
 
 
 class _RoiProcessorBase(_BaseProcessor):

setup.py

@@ -247,6 +247,7 @@ def finalize_options(self):
     def run(self):
         self.spawn(['python', 'benchmarks/benchmark_imageproc.py'])
         self.spawn(['python', 'benchmarks/benchmark_geometry.py'])
+        self.spawn(['python', 'benchmarks/benchmark_statistics.py'])
 
 
 class BinaryDistribution(Distribution):

src/extra_foam/CMakeLists.txt

@@ -9,6 +9,7 @@ set(extra-foam_MODULE_FILES
         f_imageproc.cpp
         f_datamodel.cpp
         f_geometry.cpp
+        f_statistics.cpp
 )
 
 if(UNIX)

src/extra_foam/f_geometry.cpp

@@ -68,6 +68,8 @@ PYBIND11_MODULE(geometry, m)
 {
   xt::import_numpy();
 
+  m.doc() = "Detector geometry.";
+
   declare_1MGeometry<foam::LPD_1MGeometry>(m, "LPD");
 
   declare_1MGeometry<foam::DSSC_1MGeometry>(m, "DSSC");

src/extra_foam/f_imageproc.cpp

@@ -24,7 +24,7 @@ PYBIND11_MODULE(imageproc, m)
 
   using namespace foam;
 
-  m.doc() = "Calculate the mean of images, ignoring NaNs.";
+  m.doc() = "A collection of image processing functions.";
 
   m.def("nanmeanImageArray", [] (const xt::pytensor<double, 3>& src)
     { return nanmeanImageArray(src); },

src/extra_foam/f_statistics.cpp

@@ -0,0 +1,29 @@
+/**
+ * Distributed under the terms of the BSD 3-Clause License.
+ *
+ * The full license is in the file LICENSE, distributed with this software.
+ *
+ * Author: Jun Zhu <jun.zhu@xfel.eu>
+ * Copyright (C) European X-Ray Free-Electron Laser Facility GmbH.
+ * All rights reserved.
+ */
+#include "pybind11/pybind11.h"
+
+#include "f_statistics.hpp"
+#include "f_pyconfig.hpp"
+
+namespace py = pybind11;
+
+
+PYBIND11_MODULE(statistics, m)
+{
+  xt::import_numpy();
+
+  m.doc() = "A collection of statistics functions.";
+
+  m.def("nansum", [] (const xt::pytensor<double, 2>& src) { return foam::nansum(src); });
+  m.def("nansum", [] (const xt::pytensor<float, 2>& src) { return foam::nansum(src); });
+
+  m.def("nanmean", [] (const xt::pytensor<double, 2>& src) { return foam::nanmean(src); });
+  m.def("nanmean", [] (const xt::pytensor<float, 2>& src) { return foam::nanmean(src); });
+}

src/extra_foam/include/f_statistics.hpp

@@ -0,0 +1,44 @@
+/**
+ * Distributed under the terms of the BSD 3-Clause License.
+ *
+ * The full license is in the file LICENSE, distributed with this software.
+ *
+ * Author: Jun Zhu <jun.zhu@xfel.eu>
+ * Copyright (C) European X-Ray Free-Electron Laser Facility GmbH.
+ * All rights reserved.
+ */
+#ifndef EXTRA_FOAM_F_STATISTICS_HPP
+#define EXTRA_FOAM_F_STATISTICS_HPP
+
+#include <type_traits>
+
+#include "xtensor/xreducer.hpp"
+
+#if defined(FOAM_WITH_TBB)
+#include "tbb/parallel_for.h"
+#include "tbb/blocked_range2d.h"
+#include "tbb/blocked_range3d.h"
+#endif
+
+#include "f_traits.hpp"
+
+
+namespace foam
+{
+
+template<typename E, EnableIf<std::decay_t<E>, IsImage> = false>
+inline auto nansum(E&& src)
+{
+  return xt::nansum(std::forward<E>(src), xt::evaluation_strategy::immediate)[0];
+}
+
+template<typename E, EnableIf<std::decay_t<E>, IsImage> = false>
+inline auto nanmean(E&& src)
+{
+  return xt::nanmean(std::forward<E>(src), xt::evaluation_strategy::immediate)[0];
+}
+
+} // foam
+
+
+#endif //EXTRA_FOAM_F_STATISTICS_HPP

test/CMakeLists.txt

@@ -44,7 +44,8 @@ find_package(Threads REQUIRED)
 set(FOAM_TESTS
         test_tbb.cpp
         test_imageproc.cpp
-        test_geometry.cpp)
+        test_geometry.cpp
+        test_statistics.cpp)
 
 foreach(filename IN LISTS FOAM_TESTS)
     string(REPLACE ".cpp" "" targetname ${filename})

test/test_statistics.cpp

@@ -0,0 +1,47 @@
+/**
+ * Distributed under the terms of the BSD 3-Clause License.
+ *
+ * The full license is in the file LICENSE, distributed with this software.
+ *
+ * Author: Jun Zhu <jun.zhu@xfel.eu>
+ * Copyright (C) European X-Ray Free-Electron Laser Facility GmbH.
+ * All rights reserved.
+ */
+#include "gtest/gtest.h"
+#include "gmock/gmock.h"
+
+#include "xtensor/xio.hpp"
+#include "xtensor/xtensor.hpp"
+
+#include "f_statistics.hpp"
+
+namespace foam
+{
+namespace test
+{
+
+using ::testing::ElementsAre;
+using ::testing::ElementsAreArray;
+using ::testing::NanSensitiveFloatEq;
+using ::testing::FloatEq;
+
+TEST(TestNanmean, TestGeneral)
+{
+  auto nan = std::numeric_limits<float>::quiet_NaN();
+  auto nan_mt = NanSensitiveFloatEq(nan);
+
+  xt::xtensor<float, 2> img {{1.f, -1.f, 1.f}, {4.f, 5.f, nan}};
+  EXPECT_EQ(2.f, foam::nanmean(img));
+}
+
+TEST(TestNansum, TestGeneral)
+{
+  auto nan = std::numeric_limits<float>::quiet_NaN();
+  auto nan_mt = NanSensitiveFloatEq(nan);
+
+  xt::xtensor<float, 2> img {{1.f, -1.f, 1.f}, {4.f, 5.f, nan}};
+  EXPECT_EQ(10.f, foam::nansum(img));
+}
+
+} //test
+} //foam