Merge pull request #268 from zhujun98/implement_parallel_azimuthal_in…

…tegration

Implement parallel azimuthal integration

benchmarks/azimuthal_integration.ipynb

@@ -44,12 +44,22 @@
     "\n",
     "# img, mask = load_image(\"jf_ring.npy\")\n",
     "# cy, cx = -33, 1112\n",
+    "# pixel1, pixel2 = 75e-6, 75e-6  # pixel size (y, x)\n",
     "\n",
     "# img, mask = load_image(\"jf_ring_6modules.npy\")\n",
-    "# cy, cx = 537, 1132 \n",
+    "# cy, cx = 537, 1132\n",
+    "# pixel1, pixel2 = 75e-6, 75e-6  # pixel size (y, x)\n",
     "\n",
     "img, mask = load_image(\"lpd.npy\")\n",
-    "cy, cx = 606, 554"
+    "cy, cx = 606, 554\n",
+    "pixel1, pixel2 = 200e-6, 200e-6  # pixel size (y, x)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Integrate a single image"
    ]
   },
   {
@@ -60,7 +70,6 @@
    "source": [
     "dist = 1  # sample distance\n",
     "npt = 1024  # number of integration points\n",
-    "pixel1, pixel2 = 0.75e-6, 0.75e-6  # pixel size (y, x)\n",
     "poni1, poni2 = cy * pixel1, cx * pixel2  # integration center (y, x)"
    ]
   },
@@ -126,14 +135,54 @@
     "ax.legend(fontsize=16)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Integrate an array of images"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
-    "# %%timeit\n",
+    "import multiprocessing as mp\n",
+    "\n",
+    "print(mp.cpu_count())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "img_array = np.tile(img, (40, 1, 1))\n",
+    "print(img_array.shape)\n",
+    "\n",
+    "q_a, I_a = integrator.integrate1d(img_array, npt=npt)\n",
+    "np.testing.assert_array_equal(q_a, q)\n",
+    "np.testing.assert_array_equal(I_a[0], I)\n",
+    "np.testing.assert_array_equal(I_a[39], I)\n",
     "\n",
+    "%timeit integrator.integrate1d(img_array, npt=npt)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Concentric ring finding"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
     "min_count = 500\n",
     "prominence = 100\n",
     "distance = 10\n",
@@ -142,6 +191,15 @@
     "cx, cy = finder.search(img, cx, cy, min_count=min_count)"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%timeit finder.search(img, cx, cy, min_count=min_count)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -172,9 +230,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "foam-gcc6",
+   "display_name": "Python [conda env:foam-test]",
    "language": "python",
-   "name": "foam-gcc6"
+   "name": "conda-env-foam-test-py"
   },
   "language_info": {
    "codemirror_mode": {
@@ -186,7 +244,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.6"
+   "version": "3.7.7"
   }
  },
  "nbformat": 4,

docs/installation.rst

@@ -41,19 +41,29 @@ Install **EXtra-foam**
     $ git submodule update --init
 
     $ cd EXtra-foam
+    $ export CMAKE_PREFIX_PATH=${CONDA_PREFIX:-"$(dirname $(which conda))/../"}
+    $ pip install .
 
-    # optional
-    $ export FOAM_USE_TBB=0  # turn off intel TBB in extra-foam
-    $ export XTENSOR_USE_TBB=0  # turn off intel TBB in xtensor
-    $ export FOAM_USE_XSIMD=0  # turn off XSIMD in extra-foam
-    $ export XTENSOR_USE_XSIMD=0  # turn off XSIMD in xtensor
 
-    # Note: This step is also required if one wants to change the above
-    #       environmental parameters.
-    $ python setup.py clean  # alternatively "rm -r build"
+Install your own **EXtra-foam** kernel on the Maxwell cluster for offline analysis
+----------------------------------------------------------------------------------
 
-    $ export CMAKE_PREFIX_PATH=${CONDA_PREFIX:-"$(dirname $(which conda))/../"}
-    $ pip install .
+For now, there is no documentation for the Python bindings of the C++ implementations in
+**EXtra-foam**. If you are interested in using those super fast C++ implementation to
+accelerate your analysis, please feel free to dig into the code and ask questions.
+
+.. code-block:: bash
+
+    # ssh to the Maxwell cluster and then
+    $ module load anaconda3
+    $ conda create --name extra-foam-offline -y
+    $ conda activate extra-foam-offline
+
+    # follow the previous steps to install EXtra-foam
+
+    $ conda install ipykernel nb_conda_kernels -y
+
+    # Now you should be able to load the newly created kernel on max-jhub.
 
 
 Install C++ API of **EXtra-foam** only

extra_foam/algorithms/tests/test_azimuthal_integ.py

@@ -37,8 +37,20 @@ def testIntegrate1D(self, dtype):
         assert q0 == q1
         assert s0 == s1
 
-        q10, s10 = integrator.integrate1d(img, npt=10, min_count=img.size)
-        assert not np.any(s10)
+        # test threshold
+        q10_0, s10_0 = integrator.integrate1d(img, npt=10, min_count=img.size)
+        assert not np.any(s10_0)
+
+        # test integrate an array of images
+        q10, s10 = integrator.integrate1d(img, npt=10)
+        img_a = np.tile(img, (4, 1, 1))
+        img_a[3, ...] = img - 1000
+        q10_a, s10_a = integrator.integrate1d(img_a, npt=10)
+        np.testing.assert_array_equal(q10_a, q10)
+        for i in range(3):
+            np.testing.assert_array_equal(s10_a[i], s10)
+        _, s10_2 = integrator.integrate1d(img - 1000, npt=10)
+        np.testing.assert_array_equal(s10_a[3], s10_2)
 
         q100, s100 = integrator.integrate1d(img, npt=999)
 

src/extra_foam/f_azimuthal_integrator.cpp

@@ -41,6 +41,19 @@ void declareAzimuthalIntegrator(py::module& m)
   AZIMUTHAL_INTEGRATE1D(float)
   AZIMUTHAL_INTEGRATE1D(uint16_t)
   AZIMUTHAL_INTEGRATE1D(int16_t)
+
+#define AZIMUTHAL_INTEGRATE1D_PARA(DTYPE)                                                             \
+  cls.def("integrate1d", (std::pair<foam::ReducedVectorTypeFromArray<xt::pytensor<value_type, 3>>,    \
+                                    foam::ReducedImageType<xt::pytensor<value_type, 3>>>              \
+                          (Integrator::*)(const xt::pytensor<DTYPE, 3>&, size_t, size_t,              \
+                                          foam::AzimuthalIntegrationMethod))                          \
+     &Integrator::template integrate1d<const xt::pytensor<DTYPE, 3>&>,                                \
+     py::arg("src").noconvert(), py::arg("npt"), py::arg("min_count")=1,                              \
+     py::arg("method")=foam::AzimuthalIntegrationMethod::HISTOGRAM);
+
+  AZIMUTHAL_INTEGRATE1D_PARA(float)
+  AZIMUTHAL_INTEGRATE1D_PARA(uint16_t)
+  AZIMUTHAL_INTEGRATE1D_PARA(int16_t)
 }
 
 void declareConcentricRingsFinder(py::module& m)

src/extra_foam/include/f_azimuthal_integrator.hpp

@@ -32,7 +32,7 @@ namespace foam
 namespace
 {
 
-template<typename T, typename E, EnableIf<std::decay_t<E>, IsImage> = false>
+template<typename T, typename E>
 auto computeGeometry(E&& src, T poni1, T poni2, T pixel1, T pixel2)
 {
   auto shape = src.shape();
@@ -51,20 +51,13 @@ auto computeGeometry(E&& src, T poni1, T poni2, T pixel1, T pixel2)
   return geometry;
 }
 
-template<typename T, typename E, EnableIf<std::decay_t<E>, IsImage> = false>
-auto histogramAI(E&& src, const xt::xtensor<T, 2>& geometry, T q_min, T q_max,
-                 size_t n_bins, size_t min_count=1)
+template<typename E1, typename E2, typename E3, typename T>
+void histogramAI(E1&& src, const E2& geometry, E3& hist, T q_min, T q_max, size_t n_bins, size_t min_count)
 {
-  auto shape = src.shape();
-
-  using vector_type = ReducedVectorType<E, T>;
-
   T norm = T(1) / (q_max - q_min);
-
-  vector_type edges = xt::linspace<T>(q_min, q_max, n_bins + 1);
-  vector_type hist = xt::zeros<T>({ n_bins });
   xt::xtensor<size_t, 1> counts = xt::zeros<size_t>({ n_bins });
 
+  auto shape = src.shape();
   for (size_t i = 0; i < shape[0]; ++i)
   {
     for (size_t j = 0; j < shape[1]; ++j)
@@ -103,7 +96,21 @@ auto histogramAI(E&& src, const xt::xtensor<T, 2>& geometry, T q_min, T q_max,
     else
       hist(i) /= counts(i);
   }
+}
+
+template<typename T, typename E, EnableIf<std::decay_t<E>, IsImage> = false>
+auto histogramAI(E&& src, const xt::xtensor<T, 2>& geometry, T q_min, T q_max,
+                 size_t n_bins, size_t min_count=1)
+{
+  auto shape = src.shape();
 
+  using vector_type = ReducedVectorType<E, T>;
+
+  vector_type hist = xt::zeros<T>({ n_bins });
+
+  histogramAI(std::forward<E>(src), geometry, hist, q_min, q_max, n_bins, min_count);
+
+  vector_type edges = xt::linspace<T>(q_min, q_max, n_bins + 1);
   auto&& centers = 0.5 * (xt::view(edges, xt::range(0, -1)) + xt::view(edges, xt::range(1, xt::placeholders::_)));
 
   return std::make_pair<vector_type, vector_type>(centers, std::move(hist));
@@ -119,6 +126,43 @@ auto histogramAI(E&& src, T poni1, T poni2, T pixel1, T pixel2, size_t npt, size
   return histogramAI<T>(std::forward<E>(src), geometry, bounds[0], bounds[1], npt, min_count);
 }
 
+template<typename T, typename E, EnableIf<std::decay_t<E>, IsImageArray> = false>
+auto histogramAI(E&& src, const xt::xtensor<T, 2>& geometry, T q_min, T q_max,
+                 size_t n_bins, size_t min_count=1)
+{
+  size_t np = src.shape()[0];
+
+  using vector_type = ReducedVectorTypeFromArray<E, T>;
+  using image_type = ReducedImageType<E, T>;
+
+  image_type hist = xt::zeros<T>({ np, n_bins });
+
+#if defined(FOAM_USE_TBB)
+  tbb::parallel_for(tbb::blocked_range<int>(0, np),
+    [&src, &geometry, &hist, q_min, q_max, n_bins, min_count]
+    (const tbb::blocked_range<int> &block)
+    {
+      for(int k=block.begin(); k != block.end(); ++k)
+      {
+#else
+      for (size_t k = 0; k < np; ++k)
+      {
+#endif
+        auto hist_view = xt::view(hist, k, xt::all());
+        histogramAI(xt::view(src, k, xt::all(), xt::all()), geometry, hist_view,
+                    q_min, q_max, n_bins, min_count);
+      }
+#if defined(FOAM_USE_TBB)
+    }
+  );
+#endif
+
+  vector_type edges = xt::linspace<T>(q_min, q_max, n_bins + 1);
+  auto&& centers = 0.5 * (xt::view(edges, xt::range(0, -1)) + xt::view(edges, xt::range(1, xt::placeholders::_)));
+
+  return std::make_pair<vector_type, image_type>(centers, std::move(hist));
+}
+
 } //namespace
 
 enum class AzimuthalIntegrationMethod
@@ -186,6 +230,10 @@ class AzimuthalIntegrator
   template<typename E, EnableIf<std::decay_t<E>, IsImage> = false>
   auto integrate1d(E&& src, size_t npt, size_t min_count=1,
                    AzimuthalIntegrationMethod method=AzimuthalIntegrationMethod::HISTOGRAM);
+
+  template<typename E, EnableIf<std::decay_t<E>, IsImageArray> = false>
+  auto integrate1d(E&& src, size_t npt, size_t min_count=1,
+                   AzimuthalIntegrationMethod method=AzimuthalIntegrationMethod::HISTOGRAM);
 };
 
 template<typename T>
@@ -249,6 +297,34 @@ auto AzimuthalIntegrator<T>::integrate1d(E&& src,
   }
 }
 
+template<typename T>
+template<typename E, EnableIf<std::decay_t<E>, IsImageArray>>
+auto AzimuthalIntegrator<T>::integrate1d(E&& src,
+                                         size_t npt,
+                                         size_t min_count,
+                                         AzimuthalIntegrationMethod method)
+{
+  if (npt == 0) npt = 1;
+
+  auto src_shape = src.shape();
+  std::array<size_t, 2> q_shape = q_.shape();
+  if (!initialized_ || src_shape[1] != q_shape[0] || src_shape[2] != q_shape[1])
+  {
+    initQ(xt::view(src, 0, xt::all(), xt::all()));
+    initialized_ = true;
+  }
+
+  switch(method)
+  {
+    case AzimuthalIntegrationMethod::HISTOGRAM:
+    {
+      return histogramAI(std::forward<E>(src), q_, q_min_, q_max_, npt, min_count);
+    }
+    default:
+      throw std::runtime_error("Unknown azimuthal integration method");
+  }
+}
+
 /**
  * class for finding the center of concentric rings in an image.
  */

src/extra_foam/include/f_traits.hpp

@@ -82,6 +82,11 @@ using ReducedVectorType = decltype(xt::eval(xt::sum<std::conditional_t<std::is_s
                                                                        typename std::decay_t<E>::value_type,
                                                                        T>>(std::declval<E>(), {0})));
 
+template<typename E, typename T = void, EnableIf<std::decay_t<E>, IsImageArray> = false>
+using ReducedVectorTypeFromArray = decltype(xt::eval(xt::sum<std::conditional_t<std::is_same<T, void>::value,
+                                                                                typename std::decay_t<E>::value_type,
+                                                                                T>>(std::declval<E>(), {0, 1})));
+
 template<typename E, typename T = void, EnableIf<std::decay_t<E>, IsImageArray> = false>
 using ReducedImageType = decltype(xt::eval(xt::sum<std::conditional_t<std::is_same<T, void>::value,
                                                                       typename std::decay_t<E>::value_type,

test/test_azimuthal_integrator.cpp

@@ -49,6 +49,10 @@ TEST(TestAzimuthalIntegrator, TestDataType)
 TEST(TestAzimuthalIntegrator, TestIntegrator1D)
 {
   xt::xtensor<float, 2> src = xt::arange(1024).reshape({16, 128});
+  xt::xtensor<float, 2> src2 = src - 100;
+  auto src_a = xt::xtensor<float, 3>::from_shape({4, 16, 128});
+  for (size_t i = 0; i < 3; ++i) xt::view(src_a, i, xt::all(), xt::all()) = src;
+  xt::view(src_a, 3, xt::all(), xt::all()) = src2;
 
   double distance = 0.2;
   double pixel1 = 1e-4;
@@ -69,12 +73,19 @@ TEST(TestAzimuthalIntegrator, TestIntegrator1D)
   EXPECT_EQ(ret10.first, ret10_cut.first);
   EXPECT_THAT(ret10_cut.second, Each(Eq(0.)));
 
+  // test integrate an array of images
+  auto ret10_a = itgt.integrate1d(src_a, 10);
+  EXPECT_EQ(ret10.first, ret10_a.first);
+  for (size_t i = 0; i < 3; ++i) EXPECT_EQ(ret10.second, xt::view(ret10_a.second, i, xt::all()));
+  auto ret10_2 = itgt.integrate1d(src2, 10);
+  EXPECT_EQ(ret10_2.second, xt::view(ret10_a.second, 3, xt::all()));
+
   // big npt
   itgt.integrate1d(src, 999);
 
   // data has a single value
-  xt::xtensor<double, 2> src2 = xt::ones<double>({16, 128});
-  itgt.integrate1d(src2, 10);
+  xt::xtensor<double, 2> src_single_value = xt::ones<double>({16, 128});
+  itgt.integrate1d(src_single_value, 10);
 
   // integral source value type
   xt::xtensor<uint16_t, 2> src_int = xt::arange(1024).reshape({16, 128});